Rotation merge: a stretch removed whole is taken out of the data and the sweep scaled again without it

The delta-CC1/2 stretch rule decides last, on the corrected fulls: by the time
it removes a flagged stretch - a dead stretch at a tenth of the crystal's
scale, a second crystal or a slipped orientation taking over - every per-frame
scale, correction surface and the error model had been fitted with that
stretch in, and only the error model was refitted after it. So the first merge
that may be written (in the experiment's own group, not a search or pre-pass
merge) and removes a stretch whole now clears that stretch's reflections from
the outcomes and makes itself again from the start without them; every later
merge, and an engine made again after a change of basis, keeps them out
(RotationScaleMerge::TakenOut). The run's own batches are judged again on what
is left. The ledger keeps the stretch at the per-image scale, CC and cost it
was removed on; the quality guard's CC1/2 before corrections stays that of the
merge that still held it. A copy of the engine (a merge made ahead on another
card) takes nothing out; where such a merge removed a stretch the engine has
not taken out, it is made again on the engine.

The stretch-whole rule fires on 10 of 118 sets of the battery subset (78
open/in-house, 25 private, 15 multisweep); the other 108 are identical to the
rule alone. Against the rule alone, on the 10 (--no-refine): R_meas lower on 7
(an open set 0.425 -> 0.280, ISa 17.2 -> 19.9, REFMAC R-free 0.244 -> 0.243; a
private monoclinic set CC1/2 0.868 -> 0.922; the in-house half-image lysozyme
0.620 -> 0.559), about equal on 3. Costs: a powder-ring myoglobin sweep is
cut at 1.90 instead of 1.49 A (R_meas 1.48 -> 0.85), one private triclinic set
gives up another 5 points of completeness, and the ISa of two private sets
with very high ISa falls (60.8 -> 23.3, 45.3 -> 30.9) at equal R_meas. One
extra scaling per triggered run; wall time unchanged within noise.

Processing the sweep again as segments (branch stretch-v2-seg) recovers more
on two sets whose problem is integration-side, at 2-4x the wall time.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu

docs/RUGNUX_REPORT.md

@@ -617,6 +617,13 @@ measured on those batches alone; with fewer than three of them nothing is remove
 goes first and the rest are measured again against what is left. Weak stretches that agree with the
 rest of the run — a crystal that leaves the beam and comes back — cost nothing on this measure and stay.
 
+A stretch removed whole is then out of the data, not only out of this one merge: the per-frame scales,
+the correction surfaces and the error model it was judged on were all fitted with it in. So the first
+merge that may be written and removes one takes its frames' reflections out — they are in no later merge
+and in no file written from the integrated reflections — and makes itself again from the start without
+them; the run's own batches are then judged again on what is left, which can remove more. The stretch
+stays in this table at the scale, CC and ΔCC1/2 it was removed on.
+
 **What ΔCC1/2 cannot do**, because the report must not imply otherwise:
 
 - it says nothing about the **cause**: a shutter fault and a crystal that slipped have the identical

image_analysis/scale_merge/Merge.h

@@ -157,6 +157,9 @@ struct SweepQuality {
     int frames_rejected = 0;
     float rejected_deg = 0.0f;
     std::vector<uint8_t> frame_disposition;   // FrameDisposition per processed image
+    // The stretches the per-image channels flagged that the merge removed whole (inclusive, processed-image
+    // ordinals); their reflections are then out of the data (RotationScaleMerge::TakenOut).
+    std::vector<std::pair<int, int>> removed_stretches;
 
     // The per-batch delta-CC1/2 curve the disposition was decided on, measured on the run as it stood
     // before anything was removed (NaN where a batch carries too few reflections to be judged). Same

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -558,9 +558,12 @@ void RotationScaleMerge::IngestHost() {
 #endif
 
     data_first = 0;
-    while (data_first < n_frames && partials_out[data_first].reflections.empty()) ++data_first;
+    const auto no_data = [&](int f) {
+        return partials_out[f].reflections.empty() && !(f < static_cast<int>(taken_out.frames.size()) && taken_out.frames[f]);
+    };
+    while (data_first < n_frames && no_data(data_first)) ++data_first;
     data_last = n_frames - 1;
-    while (data_last >= 0 && partials_out[data_last].reflections.empty()) --data_last;
+    while (data_last >= 0 && no_data(data_last)) --data_last;
 
     frame_cell_ok.assign(n_frames, 1);
     g_partial.assign(n_frames, 1.0);
@@ -999,6 +1002,7 @@ std::unique_ptr<RotationScaleMerge> RotationScaleMerge::CloneIngest(int device,
     c->logger = &log;
     c->gpu_device = device;
     c->write_back_per_frame_scale = false;
+    c->may_take_out = false;   // the outcomes are the original's (TakenOut)
 #ifdef JFJOCH_USE_CUDA
     if (gpu_)
         c->gpu_ = std::make_unique<RotationScaleMergeGPU>(device);
@@ -3101,6 +3105,12 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector<uint8_t> &partial
         if (have_cc && cc_n[f] >= MIN_REFLECTIONS_FOR_IMAGE_CC && std::isfinite(cc[f]))
             quality[f] = cc[f] / cc_typ;
     }
+    // A stretch taken out (TakenOut) reads as it did in the merge that removed it.
+    for (int f = 0; f < static_cast<int>(taken_out.frames.size()); ++f)
+        if (taken_out.frames[f]) {
+            scale[f] = taken_out.scale[f];
+            quality[f] = taken_out.cc[f];
+        }
 
     sweep_quality.measured = true;
     sweep_quality.sweep_deg = static_cast<float>(n_frames * osc_deg);
@@ -3417,6 +3427,9 @@ int RotationScaleMerge::MedianEventSpan(const std::vector<int32_t> &hist, int64_
 //    cycle and re-forming the reference in between is the authors' loop.
 void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
     frame_rejected.assign(n_frames, 0);
+    for (int f = 0; f < static_cast<int>(taken_out.frames.size()); ++f)
+        frame_rejected[f] = taken_out.frames[f];
+    removed_stretch_cost.clear();
     sweep_quality.delta_cc_half_batch.clear();
     sweep_quality.delta_cc_half_batch_deg = 0.0f;
     const auto gon = x.GetGoniometer();
@@ -3758,6 +3771,13 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
     // measured against the run as it now stands: the totals change only where a stretch is removed.
     std::vector<std::pair<double, double>> range_cost;   // (delta_cc, se) of each range as convicted
     std::vector<std::pair<int, int>> rejected_ranges;
+    // A stretch taken out (TakenOut) is in the ledger as removed, at the cost it was removed on.
+    for (size_t i = 0; i < taken_out.stretches.size(); ++i) {
+        rejected_ranges.push_back(taken_out.stretches[i]);
+        range_cost.push_back(taken_out.cost[i]);
+        sweep_quality.removed_stretches.push_back(taken_out.stretches[i]);
+        removed_stretch_cost.push_back(taken_out.cost[i]);
+    }
     std::vector<uint8_t> spent(candidates.size(), 0);
     int n_rejected = 0;
     while (true) {
@@ -3881,6 +3901,8 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) {
         }
         rejected_ranges.emplace_back(f0, f1);
         range_cost.emplace_back(stretch_delta[worst].delta_cc, stretch_delta[worst].se);
+        sweep_quality.removed_stretches.emplace_back(f0, f1);
+        removed_stretch_cost.emplace_back(stretch_delta[worst].delta_cc, stretch_delta[worst].se);
         logger->Info("delta-CC1/2: frames {}-{}, flagged by the per-image channels, removed whole - keeping "
                      "them lowered the plain CC1/2 by {:.4f} +/- {:.4f}", f0, f1, -stretch_delta[worst].delta_cc,
                      stretch_delta[worst].se);
@@ -7223,13 +7245,43 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
                                                   const gemmi::SpaceGroup *group) {
     const ScopedGpuPin pin(gpu_device);   // see IngestHost
     run_group = group;
+    struct ClearGroup { const gemmi::SpaceGroup *&g; ~ClearGroup() { g = nullptr; } } clear_group{run_group};
+    Result r = RunChoice(for_search, full_stats, measure_cc_before_corrections);
+    // See TakenOut. Once, and only on a merge in the experiment's own group that may be written.
+    if (may_take_out && taken_out.frames.empty() && group == nullptr && full_stats && !for_search && segments.empty()
+        && !r.statistics.sweep_quality.removed_stretches.empty()) {
+        taken_out.stretches = r.statistics.sweep_quality.removed_stretches;
+        taken_out.cost = removed_stretch_cost;
+        taken_out.scale = sweep_scale;
+        taken_out.cc = sweep_cc;
+        taken_out.frames.assign(n_frames, 0);
+        int64_t n = 0;
+        for (const auto &[f0, f1] : taken_out.stretches)
+            for (int f = f0; f <= f1; ++f) {
+                taken_out.frames[f] = 1;
+                n += static_cast<int64_t>(partials_out[f].reflections.size());
+                partials_out[f].reflections.clear();
+            }
+        logger->Info("The merge removed {} stretch(es) whole: their {} reflections are taken out of the data and "
+                     "the sweep is scaled again from the start without them", taken_out.stretches.size(), n);
+        Ingest();
+        run_group = group;
+        // The quality guard compares this CC1/2 with a pass that kept the stretch: it stays the first merge's.
+        const double cc_before = r.cc_half_before_corrections;
+        r = RunChoice(for_search, full_stats, measure_cc_before_corrections);
+        r.cc_half_before_corrections = cc_before;
+    }
+    return r;
+}
+
+RotationScaleMerge::Result RotationScaleMerge::RunChoice(bool for_search, bool full_stats,
+                                                        bool measure_cc_before_corrections) {
     // Pre-combined fulls: there are no partials to fit a scale on, and every event is one frame.
     if (!segments.empty()) {
         scale_partials_choice = false;
         rocking_event_frames_at_start = 1;
         rocking_events_at_start = n_partials_obs;
     }
-    struct ClearGroup { const gemmi::SpaceGroup *&g; ~ClearGroup() { g = nullptr; } } clear_group{run_group};
     // 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

image_analysis/scale_merge/RotationScaleMerge.h

@@ -218,6 +218,24 @@ public:
     Result Run(bool for_search, bool full_stats, bool measure_cc_before_corrections,
                const gemmi::SpaceGroup *group = nullptr);
 
+    // The stretches a merge removed whole (SweepQuality::removed_stretches) are out of the data from then on.
+    // The merge that removed them had fitted every scale, surface and error model with them in - a dead
+    // stretch at a tenth of the crystal's scale, a second crystal taking over - and only then judged them,
+    // so the first merge that may be written and removes one (not a search merge, nor a pre-pass merge)
+    // clears their frames' reflections from the outcomes and makes itself again from the start without
+    // them, as does every merge after it. Their frames keep, in the sweep's ledger, the per-image scale and
+    // CC and the cost they were removed on. A copy of the engine (CloneIngest) keeps what was taken out
+    // before it was made, and takes nothing out itself: the outcomes are not its own.
+    struct TakenOut {
+        std::vector<uint8_t> frames;                       // 1 = a frame of a stretch taken out
+        std::vector<std::pair<int, int>> stretches;        // inclusive, processed-image ordinals
+        std::vector<std::pair<double, double>> cost;       // (delta_cc, se) each was removed on
+        std::vector<double> scale, cc;                     // the ledger's two channels as they read then
+    };
+    [[nodiscard]] const TakenOut &GetTakenOut() const { return taken_out; }
+    // For an engine made again on the same outcomes (rugnux): set before Ingest().
+    void SetTakenOut(TakenOut t) { taken_out = std::move(t); }
+
     // Whether Run() writes the per-frame G / CC / mosaicity back onto the outcomes (on by default). Off
     // for a merge that is not the run's answer - the quality guard's merge of the pre-pass frames, the
     // finalist ledger's P1 merge - so the per-image table and the unmerged MTZ describe the merge that
@@ -1117,4 +1135,9 @@ private:
     // 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);
+    // Run() itself, without taking anything out.
+    Result RunChoice(bool for_search, bool full_stats, bool measure_cc_before_corrections);
+    TakenOut taken_out;                                    // see TakenOut
+    bool may_take_out = true;                              // false on a copy (CloneIngest)
+    std::vector<std::pair<double, double>> removed_stretch_cost;   // (delta_cc, se) per removed_stretches entry
 };

rugnux/RugnuxScaleMerge.cpp

@@ -1206,6 +1206,15 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
     const auto &rot_ss = experiment_.GetScalingSettings();
     const bool is_rotation = experiment_.IsRotationIndexing();   // rotation indexing -> rotation scaling/merge
     std::optional<RotationScaleMerge> rsm;
+    // The engine made again on the same frames (a change of basis, a corrected lattice): what the last one
+    // took out of the data stays out (RotationScaleMerge::TakenOut).
+    const auto remake_engine = [&] {
+        RotationScaleMerge::TakenOut taken_out = rsm->GetTakenOut();
+        rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
+                    static_cast<int>(config_.scaling_iter), config_.nthreads, logger, config_.observation_dump);
+        rsm->SetTakenOut(std::move(taken_out));
+        rsm->Ingest();
+    };
     std::optional<ScaleMergeResult> guard_merge;   // what the quality guard reads, where not the search merge
     std::optional<ScaleMergeResult> subset_merge;  // the pre-pass subset's, for a beam-centre arm (below)
     std::optional<PostRefineObservations> prepass_postrefine_obs;
@@ -3193,10 +3202,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                             .niggli_class = cand->niggli_class, .crystal_system = cand->system };
                         sg_search = s2;
                         if (rsm) {
-                            rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
-                                        static_cast<int>(config_.scaling_iter),
-                                        config_.nthreads, logger, config_.observation_dump);
-                            rsm->Ingest();
+                            remake_engine();
                             spec_same_ingest = false;
                         }
                         const auto &uc = *result.consensus_cell;
@@ -3433,10 +3439,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                         result.consensus_cell = adopted->GetUnitCell();
                         end_msg.unit_cell = result.consensus_cell;
                         if (rsm) {
-                            rsm.emplace(experiment_, indexer->GetIntegrationOutcome(),
-                                        result.consensus_cell, static_cast<int>(config_.scaling_iter),
-                                        config_.nthreads, logger, config_.observation_dump);
-                            rsm->Ingest();
+                            remake_engine();
                             spec_same_ingest = false;
                         }
                     }
@@ -3462,10 +3465,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                 end_msg.rotation_lattice_type = commit_lattice_type;
                 record_reindex(*commit_reindex);
                 if (rsm) {
-                    rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
-                                static_cast<int>(config_.scaling_iter),
-                                config_.nthreads, logger, config_.observation_dump);
-                    rsm->Ingest();
+                    remake_engine();
                     spec_same_ingest = false;
                 }
             }
@@ -3495,6 +3495,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                 }
                 mark.Stop();
                 m.log.ReplayInto(logger);
+                // A copy takes nothing out of the data (RotationScaleMerge::TakenOut): where its merge removed a
+                // stretch whole that this engine still holds, the merge is made again here, which does.
+                if (ahead && !ahead->first.statistics.sweep_quality.removed_stretches.empty()
+                    && rsm->GetTakenOut().frames.empty()) {
+                    logger.Debug("Final merge ahead: removed a stretch whole, which the engine takes out - not taken");
+                    ahead.reset();
+                }
                 if (ahead) {
                     m.rsm->WriteBackPerFrameScale();
                     rsm->TakeScalePartialsChoice(*m.rsm);
@@ -4031,10 +4038,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                 model_twin_choice.reset();
             determined_to_written = determined_to_written.combine(hkl);
             result.setting_operator = IndexTriplet(determined_to_written);
-            rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
-                        static_cast<int>(config_.scaling_iter),
-                        config_.nthreads, logger, config_.observation_dump);
-            rsm->Ingest();
+            remake_engine();
             phase("Re-merging in the reference's frame");
             sm = scale_and_merge(moved->short_name(), false);
             const auto &uc = *result.consensus_cell;

tests/RotationScaleMergeSegmentsTest.cpp

@@ -268,3 +268,40 @@ TEST_CASE("RotationScaleMerge_ReportedCCHalfIsThePlainPearson", "[RotationScale]
     CHECK(r.statistics.overall.cc_half_internal == Catch::Approx(weighted.GetCC()).epsilon(1e-9));
     CHECK(std::fabs(r.statistics.overall.cc_half - r.statistics.overall.cc_half_internal) > 1e-6);
 }
+
+// A sweep whose second half is another crystal, weaker: the per-image channels flag that half (low scale,
+// low CC to the merge), and the merge removes it whole - it costs the plain CC1/2 far beyond its noise and
+// beyond the spread of the first half's batches - and says which frames it removed so.
+TEST_CASE("RotationScaleMerge_SecondCrystalStretchRemovedWhole", "[RotationScale]") {
+    auto data = MakeSweeps({{600, 1.0, 0.0}}, 21);
+    auto other = MakeSweeps({{400, 0.3, 0.0}}, 22);   // another seed: other intensities
+    for (auto &o: other.frames) {
+        for (auto &r: o.reflections)
+            r.image_number += 600.0f;
+        data.frames.push_back(std::move(o));
+    }
+    const auto r = Merge(data, SyntheticExperiment(true), {});
+    const auto &sq = r.statistics.sweep_quality;
+    const auto &removed = sq.removed_stretches;
+    REQUIRE(!removed.empty());
+    int in_first = 0, in_second = 0;
+    for (const auto &[f0, f1]: removed)
+        for (int f = f0; f <= f1; ++f)
+            (f < 600 ? in_first : in_second)++;
+    CHECK(in_second >= 380);
+    CHECK(in_first <= 10);
+    // ... and the merge was made again without them (RotationScaleMerge::TakenOut): their reflections are out
+    // of the data, their frames rejected in the ledger.
+    for (const auto &[f0, f1]: removed)
+        for (int f = f0; f <= f1; ++f) {
+            CHECK(data.frames[f].reflections.empty());
+            CHECK(sq.frame_disposition[f] == static_cast<uint8_t>(FrameDisposition::Rejected));
+        }
+}
+
+// The same sweep, one crystal throughout and the second half only weaker: nothing to remove.
+TEST_CASE("RotationScaleMerge_WeakStretchOfOneCrystalKept", "[RotationScale]") {
+    auto data = MakeSweeps({{600, 1.0, 0.0}, {400, 0.3, 0.0}}, 21);
+    const auto r = Merge(data, SyntheticExperiment(true), {});
+    CHECK(r.statistics.sweep_quality.removed_stretches.empty());
+}