Error model: ISa is reported only where b is the solution of a fit - not beside a clamped a, not when b was held at zero, and on stills under the same rule

The reported ISa is 1/b of the merge error model var = a*s2 + b^2*<I>^2. Three outputs
printed a number that was not the result of any fit:

* A bin whose median normalised deviation lies far below the others' - observations that
  repeat their symmetry mates to the bit, as an input read twice gives - takes the whole fit
  through its 1/bd2^2 relative weight and pulls both terms towards zero. a stops on its 0.25
  bound, b^2 does not, and the remainder's 1/b reads as an I/sigma in the hundreds to
  millions; the standard-error test let it through on about half of the synthetic pools.
  Beside a clamped a the b^2 is not the solution of any fit, so it is now left unresolved.
* When b had no leverage and was held at zero, b_resolved kept its default and the report
  printed ISA 0.00 where the log said ISa was not reported. Unresolved now.
* The stills fit (MergeOnTheFly::RefineErrorModel) had no resolution test at all, so any
  tiny positive b^2 went out as ISa; on pools with no systematic error one seed in twelve
  printed ISa ~1000. It now carries the same test and rugnux reports it the same way.
* The --report-resolution refit's ISa (REFRES_ISA) was written without the gate.

No a or b changes: the merge weights, outlier band and every value except the resolved flag
are untouched, so a merge whose a sits inside its bounds reports exactly what it did. The
flag also feeds the per-frame-scale arm choice (PartialScalingWins), which is the one
decision that can move - only on a merge whose a is on a bound.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu

image_analysis/scale_merge/ErrorModel.cpp

@@ -121,6 +121,7 @@ bool ErrorModelUpdate(const std::vector<double> &bs2, const std::vector<double>
         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<double> &bs2, const std::vector<double>
         // 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;

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

image_analysis/scale_merge/Merge.cpp

@@ -245,6 +245,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &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<IntegrationOutcome> &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);

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; }

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;
     }
 

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;
         }

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;
             }

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<double> u(0.0, 1.0);
+    std::normal_distribution<double> z(0.0, 1.0);
+    std::vector<ErrorModelSample> 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<ErrorModelBinned> 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"