Rotation merge: scale from the fulls alone only where that leaves an error model

The per-frame scale was fitted on the partials from 50 rocking events per
frame and taken from the fulls alone below. The counts of small-molecule
sweeps (3-43) and proteins (5-900) overlap, and a weak protein at 36
events per frame scaled from its fulls had no resolved error model at all
(ISa undetermined, d_min 5.10 A) where its partials gave ISa 32 and 4.88 A.

The first merge of a pass that is not a space-group search is now made
both ways, and the count stands as the prior unless its merge has no
resolved error model while the other merge has one. Search merges (P1 /
subgroups) take the prior.

The two ISa values are deliberately not compared beyond that. Tried as
the arbiter, "higher ISa wins" (and weighted R_meas as fallback) agreed
with the external yardsticks on 10 of 11 crystals but chose the partials
on a 6-events-per-frame small-molecule sweep: ISa 10.5 vs 8.7, weighted
R_meas 0.102 vs 0.115, SHELXL R1 0.105 vs 0.062 - the partiality error
the partial scale absorbs is shared by symmetry mates at the same
rocking geometry, so their agreement cannot see it. Statistics of the
difference between the two arms' frame scales did not separate that
sweep from the proteins that want the partials either.

Effect: only crystals whose prior arm fails change; every small-molecule
set and every protein control keeps its arm. Private weak protein: ISa
undetermined -> 32.3, d_min 5.10 -> 4.88 A, weighted R_meas .343 -> .337.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB

docs/CHANGELOG.md

@@ -10,6 +10,7 @@
 * Rugnux claims a screw axis from a short axial row of a few weak reflections, so such crystals (e.g. P2_1 with a ~30 A unique axis) are no longer written without the screw.
 * Rugnux keeps the screw axes it found when a higher point group is adopted after the twin-law check, instead of writing the group without screws (e.g. P 4 2 2 for P4_2 2_1 2).
 * Rugnux drops a rotation reflection whose spot held a saturated (overloaded) pixel, as XDS does, and reports the count as `OBSERVATIONS_REJECTED_OVERLOAD=`.
+* Rugnux fits a rotation sweep's per-frame scale on the partial reflections, even with few reflections per frame, when scaling from the fulls alone leaves no measurable error model.
 
 ### 1.0.0-rc.173
 

docs/CPU_DATA_ANALYSIS_INTEGRATION.md

@@ -382,7 +382,9 @@ The combine groups each reflection's partials into rocking events (contiguous ru
 - **Capture-aware uncertainty.** A full captured incompletely ($f<1$) is extrapolated and biased high. The unobserved fraction is charged as an extra systematic uncertainty, $\sigma^2 \leftarrow \sigma^2 + \big(c\,(1-f)\,I\big)^2$, so the merge down-weights these extrapolated fulls and the error model treats their scatter as expected. The merge rebuilds every full's variance at the reflection's mean intensity (§10.4), and the capture term is rebuilt there too, as $\big(c\,(1-f)\,\langle I\rangle\big)^2$. It is enabled by default for the rotation path.
 - **Overloaded events.** An event in which any partial had a saturated pixel in its signal disk — or a pixel unreadable on that frame alone, beyond the run's pixel mask, which is how a detector that writes its error value for a pixel it could not count reports an overload — is dropped whole, as XDS drops an overloaded reflection. The brightest part of such a rocking curve is exactly what is missing, so neither the sum of the remaining partials nor their extrapolation by the partiality model measures the reflection: on a strongly diffracting small-molecule crystal these were the strongest low-order reflections, and they read 2–3× low. The integration keeps an overloaded partial, unfitted and flagged, only so that the event can be recognised; nothing else reads it. The count is `OBSERVATIONS_REJECTED_OVERLOAD=` in the report.
 
-The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4). Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$.
+The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4).
+
+Before the combine a per-frame scale can also be fitted on the partials themselves. That needs a frame to hold many rocking events caught at different points of their curves: within one rocking curve a change of scale and an error of the partiality model are the same thing, and on a finely sliced sparse sweep the fit takes one for the other. The rocking events per frame are the prior (the partials are scaled from 50 events per frame), but the counts of small-molecule sweeps (3–43) and proteins (5–900) overlap. So the first merge of a run that is not a space-group search is made **both ways**, with the partials scaled and with the scale taken from the fulls alone, and the prior stands unless its merge has no resolved error model ($b$ not resolved from zero: the strong equivalents do not agree to within a measurable systematic error) while the other merge has one. The two ISa values are deliberately not compared beyond that: the partiality-model error a partial scale takes up is shared by symmetry mates measured at the same rocking geometry, so their agreement cannot see it — a small-molecule sweep with six events per frame read ISa 10.5 with its partials scaled against 8.7 from the fulls alone, and refined to $R_1$ 0.105 against 0.062. The log states the choice and both arms' ISa. Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$.
 
 How smooth that scale is over the rotation is left to the data rather than to a fixed window. Each round fits every frame on its own fulls against the current reference, giving a scale $G_f$ and its information $D_f=\sum w^2c^2$; the scale is then the curve $x=\log G$ minimising $\sum_f J_f\,(x_f-y_f)^2+\lambda\sum_f(\Delta^2 x)_f^2$, with $y_f$ the frame's own fit in log scale and $J_f$ its information carried there — a penalised (Whittaker–Eilers) smoother, solved as a five-band linear system. $\lambda$ is chosen by cross-validation: blocks of frames one rocking curve wide are left out in turn and predicted from the curve through the rest (neighbours closer than a rocking curve share their measurement, since a full sums those frames). A frame of hundreds of fulls is then followed frame by frame, a frame of two or three is carried by its neighbours, a stretch with none is bridged by a straight line, and a scale that falls a hundredfold over a few degrees — an absorbing crystal turning edge-on — is followed where a window would average across it. Once the curve settles, one free fit is shrunk toward it frame by frame by how much of each frame's deviation its neighbour shares (the lag-1 covariance), which hands back a real per-frame systematic and discards fit noise.
 

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -42,9 +42,9 @@ namespace {
     // takes one for the other - G swung 0.23..1.2 with a 180 deg period on a small-molecule sweep whose
     // own frame scale varies 0.79..0.99 (XDS), and the merge carried an hkl-dependent bias that the
     // equivalents cannot see but a refinement against the known structure does (R1 0.096 -> 0.062
-    // without it). Below this the per-frame scale comes from the fulls alone, as XDS takes it. Measured
-    // on the battery: the small-molecule sweeps hold 2.6 to 23 events per frame, the protein sets 84 to
-    // 900 - and on the proteins the fulls-only scale leaves model R-free unchanged (+-0.002).
+    // without it). Below this the per-frame scale comes from the fulls alone, as XDS takes it. This is
+    // the PRIOR: Run merges both ways and overrides it where its merge has no resolved error model and
+    // the other has (PartialScalingWins).
     constexpr double MIN_EVENTS_PER_FRAME_SCALE_PARTIALS = 50.0;
     constexpr int64_t MIN_REFLECTIONS_FOR_IMAGE_CC = 20;   // below this a frame's CC means nothing
     // The per-frame scaling loop has settled when the scales moved less than this between two
@@ -308,6 +308,7 @@ void RotationScaleMerge::MeasureIncidentFlux(const std::vector<double> &mean_bkg
 }
 
 void RotationScaleMerge::Ingest() {
+    scale_partials_choice.reset();
     n_frames = static_cast<int>(partials_out.size());
     partials_released = false;
     resident_ingest = false;
@@ -5779,8 +5780,60 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     return result;
 }
 
+namespace {
+    // Whether the per-frame scale is fitted on the partials, given the merge made both ways: with the
+    // partials scaled (with) and from the fulls alone (without). The prior - the count of rocking events
+    // per frame - stands unless the merge it picks has no resolved error model while the other has: an
+    // ISa whose b is not resolved from zero means the strong equivalents do not even agree to within a
+    // systematic error that can be measured, which is a scale that failed, not a judgement call.
+    // Comparing two resolved ISa values is NOT a test of the scale: on a small-molecule sweep of six
+    // rocking events per frame the partials' scale read ISa 10.5 against 8.7 from the fulls alone, and
+    // SHELXL R1 0.105 against 0.062 - the partiality-model error the partial scale takes up is shared by
+    // symmetry mates measured at the same rocking geometry, so their agreement cannot see it.
+    bool PartialScalingWins(const RotationScaleMerge::Result &with, const RotationScaleMerge::Result &without,
+                            bool prior) {
+        const bool w_ok = with.isa_resolved && std::isfinite(with.isa) && with.isa > 0.0;
+        const bool o_ok = without.isa_resolved && std::isfinite(without.isa) && without.isa > 0.0;
+        const bool prior_ok = prior ? w_ok : o_ok;
+        const bool other_ok = prior ? o_ok : w_ok;
+        return prior_ok || !other_ok ? prior : !prior;
+    }
+}   // namespace
+
 RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_stats,
                                                   bool measure_cc_before_corrections) {
+    // Whether a frame's scale can be fitted on its partials depends on whether the partials tell a change
+    // of scale from an error of the partiality model. The count of rocking events per frame is the prior
+    // (MIN_EVENTS_PER_FRAME_SCALE_PARTIALS), but the counts overlap between small molecules (2.6-43) and
+    // proteins (4.8-900): a weak protein at 36 events per frame scaled from its fulls alone had no
+    // resolved error model at all (ISa 21 with the partials scaled). So the first merge on this ingest
+    // that is not a space-group search is made both ways, and the prior stands unless its merge failed
+    // where the other did not (PartialScalingWins); later merges keep that choice. A search merge is in
+    // P1 or a subgroup, where a reflection has two or three observations, so it takes the prior. The arm
+    // returned is the one run last, so that every per-pass state the merge leaves behind is that arm's.
+    if (rocking_event_frames_at_start < 0)
+        rocking_event_frames_at_start = RockingEventFrames(&rocking_events_at_start);
+    const double events_per_frame = n_frames > 0 ? static_cast<double>(rocking_events_at_start) / n_frames : 0.0;
+    const bool prior = events_per_frame >= MIN_EVENTS_PER_FRAME_SCALE_PARTIALS;
+    if (for_search && !scale_partials_choice)
+        return RunArm(prior, for_search, full_stats, measure_cc_before_corrections);
+    if (!scale_partials_choice) {
+        Result without = RunArm(false, for_search, full_stats, measure_cc_before_corrections);
+        Result with = RunArm(true, for_search, full_stats, measure_cc_before_corrections);
+        scale_partials_choice = PartialScalingWins(with, without, prior);
+        logger.Info("Per-frame scale: {} ({:.1f} rocking events per frame) - ISa {:.1f}{} with the partials "
+                    "scaled, {:.1f}{} from the fulls alone",
+                    *scale_partials_choice ? "fitted on the partials" : "from the fulls alone",
+                    events_per_frame, with.isa, with.isa_resolved ? "" : " (unresolved)",
+                    without.isa, without.isa_resolved ? "" : " (unresolved)");
+        if (*scale_partials_choice)
+            return with;
+    }
+    return RunArm(*scale_partials_choice, for_search, full_stats, measure_cc_before_corrections);
+}
+
+RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool for_search, bool full_stats,
+                                                     bool measure_cc_before_corrections) {
     HKLKeyGenerator keygen(merge_friedel, x.GetSpaceGroupOrP1());
 
     // Start from the corr Ingest built, so this pass runs the scaling_iter iterations it was asked for
@@ -5830,16 +5883,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
                     "curves span {} frames)", smooth_window, smooth_window * osc_deg, smooth_g_deg,
                     rocking_event_frames_at_start);
     }
-    // Without partial scaling (MIN_EVENTS_PER_FRAME_SCALE_PARTIALS) every frame keeps G = 1 here and
-    // its scale comes from the fulls alone (scale-fulls below).
-    if (rocking_event_frames_at_start < 0)
-        rocking_event_frames_at_start = RockingEventFrames(&rocking_events_at_start);
-    const double events_per_frame = n_frames > 0 ? static_cast<double>(rocking_events_at_start) / n_frames : 0.0;
-    const bool scale_partials = events_per_frame >= MIN_EVENTS_PER_FRAME_SCALE_PARTIALS;
-    if (!scale_partials)
-        logger.Info("Per-frame scale from the fulls alone: {:.1f} rocking events per frame, under the {:.0f} "
-                    "a partial's scale needs to stay apart from the partiality model", events_per_frame,
-                    MIN_EVENTS_PER_FRAME_SCALE_PARTIALS);
+    // Without partial scaling (see Run) every frame keeps G = 1 here and its scale comes from the fulls
+    // alone (scale-fulls below).
     ScalingLoopOutcome partial_loop;
 #ifdef JFJOCH_USE_CUDA
     if (gpu_active_ && scale_partials) {

image_analysis/scale_merge/RotationScaleMerge.h

@@ -327,6 +327,9 @@ private:
     // the same for every Run on the same ingest; -1 until the first Run takes it.
     int rocking_event_frames_at_start = -1;
     int64_t rocking_events_at_start = 0;   // the events that walk counted
+    // Whether the per-frame scale is fitted on the partials, decided by the first Run on this ingest
+    // (PartialScalingWins) and kept for the rest; empty until then.
+    std::optional<bool> scale_partials_choice;
 
     // Raw-hkl ordering, built ONCE by Ingest and reused: `perm` lists partial indices sorted by
     // (raw h,k,l, image_number); each distinct raw hkl is a contiguous run [rawrun_start, +count) of it.
@@ -725,4 +728,8 @@ private:
     // / merge-accumulate / R_meas reductions run there (only per-group + samples come back).
     // full_stats: see Run().
     Result MergeAndStats(int n_groups, bool for_search, bool fulls_resident, bool full_stats);
+
+    // One whole pass of Run() with the per-frame scale fitted on the partials (scale_partials) or taken
+    // from the fulls alone.
+    Result RunArm(bool scale_partials, bool for_search, bool full_stats, bool measure_cc_before_corrections);
 };