diff --git a/image_analysis/scale_merge/ErrorModel.cpp b/image_analysis/scale_merge/ErrorModel.cpp index 77999d826..a129cc3cf 100644 --- a/image_analysis/scale_merge/ErrorModel.cpp +++ b/image_analysis/scale_merge/ErrorModel.cpp @@ -121,6 +121,7 @@ bool ErrorModelUpdate(const std::vector &bs2, const std::vector if (!(Ass > 0.0)) return false; a_new = std::clamp(Bs / Ass, 0.25, 100.0); b2_new = 0.0; + fit.b_resolved = false; } else { if (!(det > 1e-10 * Ass * AII)) return false; const double a_fit = (Bs * AII - BI * AsI) / det; @@ -130,13 +131,18 @@ bool ErrorModelUpdate(const std::vector &bs2, const std::vector // Leverage is not enough: b^2 can have it and still come out at zero within its own error, // and 1/b then reads the noise in b as an I/sigma. The standard error of b^2 comes from the // bins' own scatter about the fitted line; this only decides what is printed. + // Nor is a standard error enough when a sits on a bound: the b^2 reported beside a clamped a + // is not the solution of any fit. A bin whose median deviation lies far below the others' - + // observations that repeat their mates to the bit - takes the whole fit through its 1/bd2^2 + // weight and pulls both terms towards zero; a stops at its bound, b^2 does not, and 1/b of + // the remainder reads as an I/sigma in the thousands or millions. double chi = 0.0; for (int bin = 0; bin < n_bins; ++bin) { if (!(bd2[bin] > 0.0)) continue; const double r = (bd2[bin] - a_fit * bs2[bin] - b2_fit * bI2[bin]) / bd2[bin]; chi += r * r; } - fit.b_resolved = b2_fit >= 2.0 * std::sqrt(chi / (n_bins - 2) * Ass / det); + fit.b_resolved = a_new == a_fit && b2_fit >= 2.0 * std::sqrt(chi / (n_bins - 2) * Ass / det); } const bool settled = std::fabs(a_new - a) <= 1e-4 * a && std::fabs(b2_new - b2) <= 1e-4 * b2; a = a_new; diff --git a/image_analysis/scale_merge/ErrorModel.h b/image_analysis/scale_merge/ErrorModel.h index 34ceb4a4b..d6ee315ec 100644 --- a/image_analysis/scale_merge/ErrorModel.h +++ b/image_analysis/scale_merge/ErrorModel.h @@ -35,7 +35,8 @@ struct ErrorModelFit { bool active = false; // the pool was large enough to fit double a = 1.0, b2 = 0.0; bool b_measured = true; // false: no bin reached where b could be seen; b is held at 0 - bool b_resolved = true; // b^2 stood clear of two of its standard errors + bool b_resolved = true; // b^2 stood clear of two of its standard errors, with a inside its bounds + // and b measured; the only state in which 1/b is reported as ISa }; // Scratch: the partitioned copy of the pool (reused between calls to save the allocation), with the rank diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index b9454e546..f422a2cb6 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -245,6 +245,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_a = 1.0; error_model_b = 0.0; error_model_chi2 = 0.0; + error_model_b_resolved = true; error_model_mean_I.clear(); reject_median_I.clear(); reject_up = 0.0; @@ -414,8 +415,21 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc const double det = Ass * AII - AsI * AsI; if (!(det > 1e-10 * Ass * AII)) return; - const double a = std::clamp((Bs * AII - BI * AsI) / det, 0.25, 100.0); - const double b2 = std::max((Ass * BI - AsI * Bs) / det, 0.0); + const double a_fit = (Bs * AII - BI * AsI) / det; + const double b2_fit = (Ass * BI - AsI * Bs) / det; + const double a = std::clamp(a_fit, 0.25, 100.0); + const double b2 = std::max(b2_fit, 0.0); + // Whether 1/b can be reported as ISa, decided as the rotation fit decides it (ErrorModelUpdate): + // b^2 two of its standard errors clear of zero, the error from the bins' scatter about the line, + // and a inside its bounds - beside a clamped a the b^2 above is not the solution of any fit, and + // it is a near-zero-scatter bin that puts a on its bound (see the weight floor above). + double chi = 0.0; + for (int bin = 0; bin < n_bins; ++bin) { + if (!(bd2[bin] > 0.0)) continue; + const double r = (bd2[bin] - a_fit * bs2[bin] - b2_fit * bI2[bin]) / bd2[bin]; + chi += r * r; + } + error_model_b_resolved = a == a_fit && b2_fit >= 2.0 * std::sqrt(chi / (n_bins - 2) * Ass / det); error_model_a = a; error_model_b = std::sqrt(b2); diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index cefedba66..bea32cae1 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -316,6 +316,7 @@ class MergeOnTheFly { double error_model_a = 1.0; double error_model_b = 0.0; double error_model_chi2 = 0.0; // achieved median reduced chi^2 (~1.0 = honestly calibrated sigmas) + bool error_model_b_resolved = true; // as ErrorModelFit::b_resolved; 1/b is reported as ISa only then // The (b*I)^2 term uses the reflection's *mean* intensity (constant over its // observations), so it inflates sigma without biasing the inverse-variance weights - // using the per-observation I_i instead would over-weight down-fluctuated points. @@ -366,6 +367,7 @@ public: [[nodiscard]] double ErrorModelA() const { return error_model_a; } [[nodiscard]] double ErrorModelB() const { return error_model_b; } [[nodiscard]] double ErrorModelChi2() const { return error_model_chi2; } + [[nodiscard]] bool ErrorModelBResolved() const { return error_model_b_resolved; } // Outlier rejection (driven by ScalingSettings::GetOutlierRejectNsigma) reports its count. [[nodiscard]] size_t RejectedCount() const { return reject_count; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 92e197e0e..b1722648c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -6634,7 +6634,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool pool.push_back(s); error_model_active = false; // a pool too small to fit leaves no model, not the run's own fit_error_model(pool); - if (error_model_active) + if (error_model_active && error_model_b_resolved) result.reference_isa = ToXdsErrorModel(error_model_a, error_model_b).isa; } diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index e3021e5dc..8d07a502e 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1507,6 +1507,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); result.error_model_isa = em.isa; + result.error_model_isa_resolved = merge_engine.ErrorModelBResolved(); result.error_model_a = em.a; result.error_model_b = em.b; } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 6b1114868..938e87c92 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2631,6 +2631,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); error_model_isa = em.isa; + error_model_isa_resolved = merge_engine.ErrorModelBResolved(); error_model_a = em.a; error_model_b = em.b; } diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index b7260c01f..1e9050001 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -534,10 +534,37 @@ TEST_CASE("ErrorModel_BNotMeasuredOnWeakData") { const auto fit = FitErrorModel(pool, scratch, 4); REQUIRE(fit.active); CHECK(!fit.b_measured); + CHECK(!fit.b_resolved); // b held at zero is not a b to report 1/b of CHECK(fit.b2 == 0.0); CHECK(fit.a == Catch::Approx(0.9).epsilon(0.03)); } +// Nearly half the observations repeat their mates to the bit (an input read twice): the bins' median deviations +// sit near zero, their 1/bd2^2 weights take the fit, and both terms are pulled towards zero. a stops on +// its lower bound while b^2 goes on to a remainder whose 1/b is an I/sigma in the hundreds or thousands. Beside a +// clamped a the b^2 is not the solution of any fit, so it is left unresolved and no ISa is printed. +TEST_CASE("ErrorModel_HalfTheMatesIdenticalLeavesISaUndetermined") { + std::mt19937 rng(1); + std::uniform_real_distribution u(0.0, 1.0); + std::normal_distribution z(0.0, 1.0); + std::vector pool; + for (int i = 0; i < 200000; ++i) { + const double s2 = std::pow(10.0, 2.0 * u(rng)); + const double snr = 300.0 * std::pow(10.0, -4.0 * u(rng)); + const double I2 = snr * snr * s2; + const double e = z(rng); + double dev2 = (1.3 * s2 + 0.03 * 0.03 * I2) * e * e; + if (u(rng) < 0.48) dev2 = 0.0; + pool.push_back({s2, I2, dev2, 2.0f}); + } + std::vector scratch; + const auto fit = FitErrorModel(pool, scratch, 4); + REQUIRE(fit.active); + CHECK(fit.a == 0.25); + CHECK(1.0 / std::sqrt(fit.b2) > 100.0); + CHECK(!fit.b_resolved); +} + #ifdef JFJOCH_USE_CUDA #include "../common/CUDAWrapper.h" #include "../image_analysis/scale_merge/ErrorModelGPU.h"