diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ce7262e5c..e5e0c1523 100644 --- a/docs/RUGNUX_REPORT.md +++ b/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 diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 84c7c82ea..407eaa345 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -157,6 +157,9 @@ struct SweepQuality { int frames_rejected = 0; float rejected_deg = 0.0f; std::vector 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> 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 diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 42133bbbe..3c2ef36d4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/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(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::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(device); @@ -3101,6 +3105,12 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector &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(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(n_frames * osc_deg); @@ -3417,6 +3427,9 @@ int RotationScaleMerge::MedianEventSpan(const std::vector &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(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> range_cost; // (delta_cc, se) of each range as convicted std::vector> 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 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(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 diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 48d65db27..0671bf5e4 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/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 frames; // 1 = a frame of a stretch taken out + std::vector> stretches; // inclusive, processed-image ordinals + std::vector> cost; // (delta_cc, se) each was removed on + std::vector 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> removed_stretch_cost; // (delta_cc, se) per removed_stretches entry }; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 76f190707..f708d2d71 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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 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(config_.scaling_iter), config_.nthreads, logger, config_.observation_dump); + rsm->SetTakenOut(std::move(taken_out)); + rsm->Ingest(); + }; std::optional guard_merge; // what the quality guard reads, where not the search merge std::optional subset_merge; // the pre-pass subset's, for a beam-centre arm (below) std::optional 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(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(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(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(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; diff --git a/tests/RotationScaleMergeSegmentsTest.cpp b/tests/RotationScaleMergeSegmentsTest.cpp index 01f601b17..c5f20c7f9 100644 --- a/tests/RotationScaleMergeSegmentsTest.cpp +++ b/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(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()); +}