From 3ba8ba1420fc4da231f4ad3e8440885bb5f82b2b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 11 Nov 2026 13:07:26 +0100 Subject: [PATCH] Model step: one set of parts for rotation and live stills Taken from rc186 (270ae432d, net of merge 2250e4019): rugnux/ModelStep, the rotation pipeline's model validation and refinement calls moved onto its parts (RugnuxScaleMerge.cpp), StillsLiveMerge ::SpaceGroup, SecondsSince in common/time_utc.h, and the two ModelStep cases of ModelValidationTest. 270ae432d states rotation output with a model is unchanged by the move (identical maps, models and blobs end to end). Left out: nothing of these files. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- common/time_utc.h | 4 + rugnux/CMakeLists.txt | 2 + rugnux/ModelStep.cpp | 256 ++++++++++++++++++++++++++++++++++ rugnux/ModelStep.h | 101 ++++++++++++++ rugnux/RugnuxScaleMerge.cpp | 52 ++----- rugnux/StillsLiveMerge.cpp | 9 +- rugnux/StillsLiveMerge.h | 2 + tests/ModelValidationTest.cpp | 145 +++++++++++++++++++ 8 files changed, 526 insertions(+), 45 deletions(-) create mode 100644 rugnux/ModelStep.cpp create mode 100644 rugnux/ModelStep.h diff --git a/common/time_utc.h b/common/time_utc.h index bec55fd75..2dfaedb20 100644 --- a/common/time_utc.h +++ b/common/time_utc.h @@ -12,6 +12,10 @@ #include #include +inline double SecondsSince(std::chrono::steady_clock::time_point t0) { + return std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); +} + inline std::string time_UTC(const std::chrono::time_point &input) { auto time_ms = std::chrono::duration_cast(input.time_since_epoch()).count(); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index bd1462509..55379b405 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -25,6 +25,8 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxLive.h StillsLiveMerge.cpp StillsLiveMerge.h + ModelStep.cpp + ModelStep.h RugnuxRefine.cpp RugnuxRefine.h Rugnux.h diff --git a/rugnux/ModelStep.cpp b/rugnux/ModelStep.cpp new file mode 100644 index 000000000..4d3d338bb --- /dev/null +++ b/rugnux/ModelStep.cpp @@ -0,0 +1,256 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ModelStep.h" + +#include +#include + +#include "spdlog/fmt/fmt.h" + +#include "../common/time_utc.h" +#include "../image_analysis/scale_merge/CrystalSetting.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" + +namespace { + bool Cut(const ModelStepSettings &settings) { + return settings.cut != nullptr && settings.cut->load(); + } + + void Step(const ModelStepSettings &settings, const std::string &title) { + if (settings.step) + settings.step(title); + } +} + +bool MergeIntoBasis(std::vector &merged, UnitCell &cell, gemmi::SpaceGroup &space_group, + const gemmi::Op &cob, bool merge_friedel, double rfree_fraction, size_t nthreads) { + const gemmi::SpaceGroup *moved = SpaceGroupInBasis(space_group, cob); + if (moved == nullptr) + return false; + merged = ReindexMergedIntoAsu(merged, HklOperator(cob), *moved, merge_friedel); + cell = CellInBasis(cell, cob); + AssignRfreeFlags(merged, *moved, rfree_fraction, 500, cell, nthreads); + space_group = *moved; + return true; +} + +ModelValidationResult ValidatePlacement(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + const gemmi::Structure &placed, double wavelength_A, + const std::vector &report_shell_d_min, size_t nthreads, + size_t gpu_resident_bytes, const ModelValidationResult *rigid_body, + Logger &logger) { + ModelValidationSchedule placed_elsewhere; + placed_elsewhere.fit_test = false; + placed_elsewhere.gpu_resident_bytes = gpu_resident_bytes; + ModelValidationResult v = ValidateAgainstModel(merged, cell, model, logger, &space_group, + /*probe_indexing_ambiguity=*/false, nthreads, wavelength_A, + report_shell_d_min, placed_elsewhere, nullptr, &placed); + if (v.ok && rigid_body != nullptr) { + v.fit_tested = rigid_body->fit_tested; + v.model_fits = rigid_body->model_fits; + v.null_r_work_mean = rigid_body->null_r_work_mean; + v.null_r_work_sd = rigid_body->null_r_work_sd; + v.r_work_sigma = rigid_body->r_work_sigma; + v.null_replicates = rigid_body->null_replicates; + } + return v; +} + +RefinementResult RefinePlacedModel(const gemmi::Structure &placed, const std::vector &merged, + const UnitCell &cell, const gemmi::SpaceGroup &space_group, double wavelength_A, + const RefinementSettings &settings, Logger &logger) { + try { + return RefineModel(placed, merged, gemmi::UnitCell(cell), space_group, wavelength_A, logger, settings); + } catch (const std::exception &e) { + RefinementResult failed; + failed.failure_reason = std::string("refinement failed: ") + e.what(); + return failed; + } +} + +ModelStepGroup ModelStepSpaceGroup(const std::optional &collection, const UnitCell &consensus, + const RugnuxModel &model) { + if (collection) + return {collection, ""}; + const gemmi::Structure &st = model.structure; + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + if (sg == nullptr || !st.cell.is_crystal()) + return {std::nullopt, "the collection states no space group, and the model gives no crystal's"}; + if (!CellsCorrespond(st.cell, gemmi::UnitCell(consensus))) + return {std::nullopt, fmt::format("the collection states no space group, and the model's cell {:.2f} {:.2f} " + "{:.2f} {:.1f} {:.1f} {:.1f} ({}) is not the run's {:.2f} {:.2f} {:.2f} " + "{:.1f} {:.1f} {:.1f}", st.cell.a, st.cell.b, st.cell.c, st.cell.alpha, + st.cell.beta, st.cell.gamma, sg->xhm(), consensus.a, consensus.b, + consensus.c, consensus.alpha, consensus.beta, consensus.gamma)}; + return {*sg, ""}; +} + +ModelStepResult RunModelStep(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + double wavelength_A, bool merge_friedel, + const std::vector &report_shell_d_min, + const ModelStepSettings &settings, const ModelStepResult *earlier, + Logger &logger) { + const auto t0 = std::chrono::steady_clock::now(); + ModelStepResult out; + out.merged = merged; + out.cell = cell; + out.space_group = space_group; + + // Warm: in the earlier step's frame, from where it placed the model + if (earlier && earlier->placed_model && earlier->validation) { + std::vector moved = out.merged; + UnitCell moved_cell = cell; + gemmi::SpaceGroup moved_group = space_group; + const bool setting_moved = !(earlier->to_model_setting == gemmi::Op::identity()); + if (setting_moved) + MergeIntoBasis(moved, moved_cell, moved_group, earlier->to_model_setting, merge_friedel, + settings.rfree_fraction, settings.nthreads); + const ReindexChoice c = ChooseReindex(moved, moved_cell, moved_group, + [&](const std::vector &m) { + const double cc = ReferenceIntensityCC(m, earlier->merged, moved_group); + return std::isfinite(cc) ? cc : -1.0; + }); + if (!c.is_identity) { + moved = ReindexMergedIntoAsu(moved, c.op, moved_group, merge_friedel); + moved_cell = CellInBasis(moved_cell, c.op); + } + Step(settings, "Model validation"); + ModelValidationResult v = ValidatePlacement(moved, moved_cell, earlier->space_group, model, + *earlier->placed_model->model, wavelength_A, report_shell_d_min, + settings.nthreads, settings.gpu_resident_bytes, nullptr, logger); + if (v.ok && v.placed_model) { + logger.Info("Model: placed as on the earlier merge{} - R-work {:.3f}, R-free {:.3f}", + c.is_identity ? "" : ", the merge reindexed into its frame", v.r_work, v.r_free); + KeepModelVerdict(v, *earlier->validation); + out.merged = std::move(moved); + out.cell = moved_cell; + out.space_group = earlier->space_group; + out.to_model_setting = earlier->to_model_setting; + out.placement = earlier->placement; + out.mr = earlier->mr; + out.warm = true; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, out.cell, out.space_group}; + out.validation = std::move(v); + } else + logger.Warning("Model: the earlier placement does not validate on this merge ({}) - placed again", + v.failure_reason); + } + + if (!out.warm) { + const bool predicted = PrepareSearchModel(model).predicted; + // An experimental structure as given, by its rigid body; a prediction has no placement of its own + if (!predicted) { + Step(settings, "Model validation"); + ModelValidationSchedule schedule; + schedule.gpu_resident_bytes = settings.gpu_resident_bytes; + ModelValidationResult v = ValidateAgainstModel(out.merged, cell, model, logger, &space_group, + /*probe_indexing_ambiguity=*/true, settings.nthreads, + wavelength_A, report_shell_d_min, schedule); + if (v.ok && v.model_fits) { + out.space_group = *AdoptModelFrame(v, out.merged, space_group, merge_friedel, logger); + // The cell moves with an indexing that is only a near-symmetry, and a model on other axes takes the merge + // to them; the maps and the placed model are then made again, there. As a rotation run does. + bool moved = false; + if (!(v.indexing_op == gemmi::Op::identity())) { + const UnitCell before = out.cell; + out.cell = CellInBasis(out.cell, HklOperator(v.indexing_op)); + moved = !(out.cell == before); + } + gemmi::Op to_model = v.change_of_basis_op; + to_model.tran = {0, 0, 0}; + if (!(to_model == gemmi::Op::identity()) + && MergeIntoBasis(out.merged, out.cell, out.space_group, to_model.inverse(), merge_friedel, + settings.rfree_fraction, settings.nthreads)) { + out.to_model_setting = to_model.inverse(); + moved = true; + } + if (moved) { + ModelValidationSchedule kept_verdict; + kept_verdict.fit_test = false; + kept_verdict.gpu_resident_bytes = settings.gpu_resident_bytes; + ModelValidationResult remade = ValidateAgainstModel(out.merged, out.cell, model, logger, &out.space_group, + /*probe_indexing_ambiguity=*/false, settings.nthreads, + wavelength_A, report_shell_d_min, kept_verdict); + if (remade.ok) { + KeepModelVerdict(remade, v); + v = std::move(remade); + } + } + out.placement = "RIGID_BODY"; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, out.cell, out.space_group}; + } + v.maps.reset(); + out.validation = std::move(v); + } + // Molecular replacement, where the rigid body did not place it + if (out.placement == "NOT_PLACED" && Cut(settings)) + out.left_out.emplace_back("molecular replacement"); + else if (out.placement == "NOT_PLACED") { + Step(settings, "Molecular replacement"); + MolecularReplacementResult mr = RunMolecularReplacement( + model, {{"MERGE", &out.merged, cell, space_group}}, logger, settings.nthreads, + settings.gpu_resident_bytes); + if (mr.chosen >= 0) { + const MRCandidate &c = mr.candidates[mr.chosen]; + if (mr.decision == "ENANTIOMORPH") + CheckEnantiomorphHand(mr, out.merged, space_group, logger, settings.nthreads, + settings.gpu_resident_bytes); + // The merge labelled in the group the model chose, where it chose one of those the merge cannot + // tell apart (the enantiomorph past its anomalous veto, or another screw) + if (mr.decision == "ENANTIOMORPH" || mr.decision == "SPACE_GROUP") { + out.space_group = c.space_group; + const gemmi::GroupOps ops = c.space_group.operations(); + std::erase_if(out.merged, [&](const MergedReflection &r) { + return ops.is_systematically_absent({{r.h, r.k, r.l}}); + }); + } + ModelValidationResult v = ValidatePlacement(out.merged, cell, out.space_group, model, *c.placed, + wavelength_A, report_shell_d_min, settings.nthreads, + settings.gpu_resident_bytes, + out.validation ? &*out.validation : nullptr, logger); + if (v.ok && v.placed_model) { + out.placement = "MR"; + out.maps = std::move(v.maps); + out.placed_model = PlacedModelOutput{v.placed_model, cell, out.space_group}; + out.validation = std::move(v); + } + } + out.mr = std::move(mr); + } + } + out.placement_s = SecondsSince(t0); + + if (!settings.refine || !out.placed_model) + return out; + if (Cut(settings)) { + out.left_out.emplace_back("refinement"); + return out; + } + Step(settings, "Refinement"); + const auto t1 = std::chrono::steady_clock::now(); + RefinementSettings rs; + rs.nthreads = settings.nthreads; + rs.gpu_resident_bytes = settings.gpu_resident_bytes; + rs.cancel = settings.cut; + // The free set is drawn by the merge: a model placed as given has most likely been refined against all of it + rs.as_given_competes = out.placement != "RIGID_BODY"; + out.refinement = RefinePlacedModel(*out.placed_model->model, out.merged, out.cell, out.space_group, wavelength_A, rs, + logger); + out.refinement_s = SecondsSince(t1); + return out; +} + +ProcessResult ModelStepReport(const ModelStepResult &result) { + ProcessResult r; + r.model_validation = result.validation; + r.model_prior = result.mr; + r.model_placement = result.placement; + r.refinement = result.refinement; + return r; +} diff --git a/rugnux/ModelStep.h b/rugnux/ModelStep.h new file mode 100644 index 000000000..1c0319eef --- /dev/null +++ b/rugnux/ModelStep.h @@ -0,0 +1,101 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include + +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/MolecularReplacement.h" +#include "../image_analysis/structure_refinement/Refinement.h" + +// The model step: a model placed against a merge - by its rigid body where it fits, by molecular replacement where +// it does not - then its maps and its refinement. No file is read or written here. +// +// A rotation run takes these parts one at a time, in Rugnux::ScaleMergeAndSymmetry, because each decision there moves what the +// run writes (the integrated observations, the twin hypothesis' merge, a rival refined beside it) and a space group, +// free set or merge it may search beside is its own. RunModelStep takes them in one go on a merge that is all there is +// of the run: the live stills of the broker, which make the step at each provisional merge and after the last frame. +// What to reuse from an earlier step and when to stop is the caller's. + +// The merge moved to the setting `cob` takes its group to: the indices, the cell, the group and, as the free set is a +// function of the written index (RfreeFlags.h), the free set. False, with nothing moved, where it reaches no setting. +bool MergeIntoBasis(std::vector &merged, UnitCell &cell, gemmi::SpaceGroup &space_group, + const gemmi::Op &cob, bool merge_friedel, double rfree_fraction, size_t nthreads); + +// The placement `placed` validated against the merge it was chosen on, for its maps, R-factors and CC(model, data), with +// nothing left for the validation to decide. `rigid_body`, where there is one, keeps its verdict on the model as given. +ModelValidationResult ValidatePlacement(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + const gemmi::Structure &placed, double wavelength_A, + const std::vector &report_shell_d_min, size_t nthreads, + size_t gpu_resident_bytes, const ModelValidationResult *rigid_body, + Logger &logger); + +// The refinement of a placed model. It decides nothing for the run, so one that fails is reported, in the result. +RefinementResult RefinePlacedModel(const gemmi::Structure &placed, const std::vector &merged, + const UnitCell &cell, const gemmi::SpaceGroup &space_group, double wavelength_A, + const RefinementSettings &settings, Logger &logger); + +// The group a merge is merged in where a model is given and no determination of the group was made (serial frames): the collection's, where it states +// one; else the model's own, where the model's cell describes the run's consensus cell on the same axes +// (CellsCorrespond); else none, and `reason` says why the model is left out. +struct ModelStepGroup { + std::optional space_group; + std::string reason; +}; +ModelStepGroup ModelStepSpaceGroup(const std::optional &collection, const UnitCell &consensus, + const RugnuxModel &model); + +struct ModelStepSettings { + size_t nthreads = 1; + size_t gpu_resident_bytes = 0; // held on the card meanwhile (the analysis's indexers) + double rfree_fraction = 0.05; // of a free set drawn again in the model's setting + bool refine = true; + // Set from any thread: the step is cut - a refinement that runs is cut short, keeping the trials that ran + // (RefinementSettings::cancel), and a molecular replacement or a refinement that has not started is left out. + const std::atomic *cut = nullptr; + // Called as each part starts: "Model validation", "Molecular replacement", "Refinement" + std::function step; +}; + +struct ModelStepResult { + // The merge in the model's frame - the indexing the model chose (AdoptModelFrame) on the axes the model is + // written on, labelled in the group the model was placed in (its enantiomorph, where it was taken) - and its cell. + std::vector merged; + UnitCell cell{}; + gemmi::SpaceGroup space_group; + gemmi::Op to_model_setting = gemmi::Op::identity(); // the change of basis that took the merge to those axes + std::string placement = "NOT_PLACED"; // RIGID_BODY | MR | NOT_PLACED + bool warm = false; // placed from an earlier step's placement + std::optional validation; // its maps moved to `maps` + std::optional mr; + // Where the model was placed: the maps of the placement and the model as placed, in the merge's frame + std::shared_ptr maps; + std::optional placed_model; + std::optional refinement; // its model, maps and blobs where it ran + std::vector left_out; // what a cut left out + double placement_s = 0, refinement_s = 0; +}; + +// The model placed against `merged` (in `space_group`, the cell `cell`): by its rigid body where it fits these +// data as given, its indexing ambiguity probed and the merge put in the model's indexing and on its axes; by molecular +// replacement where it does not (or is a prediction); then refined, where placed. `earlier`, the result of a +// step on an earlier merge of the same collection, is the warm start: the merge is put in its frame (the +// indexing a later merge happened to take is decided against it, ChooseReindex), the model goes in where it was +// placed there, and the verdict, the indexing and the molecular replacement's choice are its - no null, probe +// or search is made again. Where the warm placement does not validate, the step is made cold. +ModelStepResult RunModelStep(const std::vector &merged, const UnitCell &cell, + const gemmi::SpaceGroup &space_group, const RugnuxModel &model, + double wavelength_A, bool merge_friedel, + const std::vector &report_shell_d_min, + const ModelStepSettings &settings, const ModelStepResult *earlier, + Logger &logger); + +// The step's outcome in the terms of a ProcessResult (model_validation, model_prior, model_placement and +// refinement), for what reports one (LiveModel). +ProcessResult ModelStepReport(const ModelStepResult &result); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index dbf89d206..b1646a36a 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3,6 +3,7 @@ #include #include "Rugnux.h" +#include "ModelStep.h" #include "../image_analysis/structure_refinement/ModelValidation.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" @@ -4379,13 +4380,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { gemmi::Op to_model = ahead.change_of_basis_op; to_model.tran = {0, 0, 0}; const gemmi::Op cob = to_model.inverse(); - const gemmi::SpaceGroup *in = &sg; - if (const gemmi::SpaceGroup *moved = SpaceGroupInBasis(sg, cob)) { - merged = ReindexMergedIntoAsu(merged, HklOperator(cob), *moved, friedel); - cell = CellInBasis(cell, cob); - AssignRfreeFlags(merged, *moved, rfree_fraction, 500, cell, config_.nthreads); - in = moved; - } + MergeIntoBasis(merged, cell, sg, cob, friedel, rfree_fraction, config_.nthreads); ModelValidationSchedule schedule; schedule.maps_gate = gate; schedule.fit_test = false; // the verdict is the first validation's (KeepModelVerdict) @@ -4393,7 +4388,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { schedule.gpu_beside_bytes = beside; ahead.merged = std::move(merged); return ValidateAgainstModel(ahead.merged, cell, *config_.model, ahead.log, - in, /*probe_indexing_ambiguity=*/false, + &sg, /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), wavelength, report_shell_d_min, schedule); }); @@ -5482,27 +5477,15 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const ModelPriorFile *own = model_prior_file ? &*model_prior_file : nullptr; const gemmi::SpaceGroup data_sg = own ? own->space_group : experiment_.GetSpaceGroupOrP1(); const UnitCell cell = own ? own->cell : *result.consensus_cell; - ModelValidationSchedule placed_elsewhere; - placed_elsewhere.fit_test = false; - placed_elsewhere.gpu_resident_bytes = merge_on_card; phase("Validating against model (molecular replacement)"); - ModelValidationResult v = ValidateAgainstModel( - own ? own->merged : sm.merged, cell, *config_.model, logger, &data_sg, - /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), - experiment_.GetWavelength_A(), report_shell_d_min, placed_elsewhere, nullptr, c.placed.get()); + ModelValidationResult v = ValidatePlacement( + own ? own->merged : sm.merged, cell, data_sg, *config_.model, *c.placed, + experiment_.GetWavelength_A(), report_shell_d_min, static_cast(config_.nthreads), + merge_on_card, result.model_validation ? &*result.model_validation : nullptr, logger); // Its maps go beside the files of the merge it was validated against. take_maps(own ? result.outputs->model_prior_maps : result.outputs->maps, v); v.maps_suffix = result.model_prior_file_suffix; if (v.ok) { - if (result.model_validation) { - const ModelValidationResult &rb = *result.model_validation; - v.fit_tested = rb.fit_tested; - v.model_fits = rb.model_fits; - v.null_r_work_mean = rb.null_r_work_mean; - v.null_r_work_sd = rb.null_r_work_sd; - v.r_work_sigma = rb.r_work_sigma; - v.null_replicates = rb.null_replicates; - } // Refined where it was chosen: in the group of the model's choice, against that choice's merge. if (own) { model_prior_file->placed = v.placed_model; @@ -5672,15 +5655,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { refine_merged.push_back(r); } // processing-inert: a refinement that fails is reported, and the run goes on to its other outputs - RefinementResult refined; - try { - refined = RefineModel(*p.model, refine_merged.empty() ? *placed_merged : refine_merged, - gemmi::UnitCell(p.cell), refine_sg, experiment_.GetWavelength_A(), logger, - settings); - } catch (const std::exception &e) { - refined = RefinementResult(); - refined.failure_reason = std::string("refinement failed: ") + e.what(); - } + RefinementResult refined = RefinePlacedModel(*p.model, refine_merged.empty() ? *placed_merged : refine_merged, + p.cell, refine_sg, experiment_.GetWavelength_A(), settings, logger); // The rival hypothesis (see RivalPlacement), refined the same way, and the two decided on the refined // R-free: each refinement has made its own twin decision against its own null, so what is compared is // the best each hypothesis can do. The validation's choice stands unless the rival beats it by @@ -5689,14 +5665,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { phase("Refining the model (rival symmetry hypothesis)"); const RivalPlacement &rv = *rival; const std::vector &rv_merged = rv.prior ? rv.prior->merged : *rv.merged; - RefinementResult other; - try { - other = RefineModel(*rv.placed.model, rv_merged, gemmi::UnitCell(rv.placed.cell), - rv.placed.space_group, experiment_.GetWavelength_A(), logger, settings); - } catch (const std::exception &e) { - other = RefinementResult(); - other.failure_reason = std::string("refinement failed: ") + e.what(); - } + RefinementResult other = RefinePlacedModel(*rv.placed.model, rv_merged, rv.placed.cell, rv.placed.space_group, + experiment_.GetWavelength_A(), settings, logger); ProcessResult::RivalRefinement rr; rr.space_group = rv.placed.space_group.xhm(); rr.merge = rv.prior ? "TWIN_SUBGROUP_MERGE" : "MERGE"; diff --git a/rugnux/StillsLiveMerge.cpp b/rugnux/StillsLiveMerge.cpp index b4a6c68cb..5f9de2a3b 100644 --- a/rugnux/StillsLiveMerge.cpp +++ b/rugnux/StillsLiveMerge.cpp @@ -7,6 +7,7 @@ #include #include "../common/ParallelFor.h" +#include "../common/time_utc.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/ReferenceFreeReindex.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" @@ -14,10 +15,6 @@ #include "../image_analysis/scale_merge/SublatticeRelabel.h" namespace { - double SecondsSince(std::chrono::steady_clock::time_point t0) { - return std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - } - size_t ImagesWithReflections(const std::vector &outcomes) { size_t n = 0; for (const auto &o : outcomes) @@ -95,6 +92,10 @@ void StillsLiveMerge::Cell(const UnitCell &cell) { experiment.SetUnitCell(cell); } +void StillsLiveMerge::SpaceGroup(const gemmi::SpaceGroup &space_group) { + experiment.SetSpaceGroup(space_group); +} + void StillsLiveMerge::Add(int64_t image, IntegrationOutcome &&outcome) { ++added; if (outcome.reflections.empty()) diff --git a/rugnux/StillsLiveMerge.h b/rugnux/StillsLiveMerge.h index 43f2f3a46..ee8f841b9 100644 --- a/rugnux/StillsLiveMerge.h +++ b/rugnux/StillsLiveMerge.h @@ -50,6 +50,8 @@ public: // The cell merged in, where the experiment had none; before the first checkpoint. void Cell(const UnitCell &cell); + // The group merged in, where the experiment had none; before the first checkpoint. + void SpaceGroup(const gemmi::SpaceGroup &space_group); // Frame `image` as analysed: without reflections where it was not indexed. void Add(int64_t image, IntegrationOutcome &&outcome); // Frames added so far diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index e8986a690..60e676d4c 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -34,6 +34,7 @@ #include "../writer/WriteModel.h" #include "../writer/RugnuxModelFile.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../rugnux/ModelStep.h" namespace { // A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with @@ -2068,3 +2069,147 @@ TEST_CASE("ModelStructureFactorsGPU_MapMatchesCPU", "[ModelValidation][gpu]") { #endif + +// A model on other axes than the merge's takes the merge to them: its indices, cell, group and free set, as a rotation run +// takes its files. A later merge of the collection follows it there. +TEST_CASE("ModelStep_MergeTakesTheModelsAxes", "[ModelValidation][RugnuxLive]") { + Logger logger("ModelStep_MergeTakesTheModelsAxes"); + const auto data_model = WriteTemp("step_axes_data_test.pdb", ClusterPdb().c_str()); + auto obs = ModelReferenceIntensities(ReadModel(data_model), {}, {}, 2.5, logger); + REQUIRE(obs.size() > 500); + for (size_t i = 0; i < obs.size(); i++) { + obs[i].F = std::sqrt(std::max(0.0f, obs[i].I)); + obs[i].sigma = 0.02f * obs[i].I + 1.0f; + obs[i].sigmaF = 0.02f * obs[i].F + 0.1f; + obs[i].rfree_flag = (i % 20) == 0; + } + const UnitCell data_cell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup data_sg = *gemmi::find_spacegroup_by_name("P 21 21 21"); + + // a' = b, b' = c, c' = a + std::string permuted = "CRYST1 34.000 38.000 30.000 90.00 90.00 90.00 P 21 21 21 4\n"; + std::istringstream in(ClusterPdb()); + std::string line; + char buf[96]; + while (std::getline(in, line)) { + if (line.compare(0, 4, "ATOM") != 0) + continue; + std::snprintf(buf, sizeof buf, "%s%8.3f%8.3f%8.3f%s", line.substr(0, 30).c_str(), std::stod(line.substr(38, 8)), + std::stod(line.substr(46, 8)), std::stod(line.substr(30, 8)), line.substr(54).c_str()); + permuted += buf; + permuted += "\n"; + } + permuted += "END\n"; + const auto path = WriteTemp("step_axes_model_test.pdb", permuted.c_str()); + const RugnuxModel model = ReadModel(path); + + ModelStepSettings settings; + settings.refine = false; + const ModelStepResult first = RunModelStep(obs, data_cell, data_sg, model, 1.0, true, {}, settings, nullptr, logger); + CHECK(first.placement == "RIGID_BODY"); + CHECK_FALSE(first.to_model_setting == gemmi::Op::identity()); + CHECK(first.cell.a == Catch::Approx(34.0)); + CHECK(first.cell.b == Catch::Approx(38.0)); + CHECK(first.cell.c == Catch::Approx(30.0)); + REQUIRE(first.placed_model); + CHECK(first.placed_model->cell.a == Catch::Approx(34.0)); + REQUIRE(first.validation); + CHECK(first.validation->r_work < 0.15); + CHECK(first.validation->model_fits); + + const ModelStepResult second = RunModelStep(obs, data_cell, data_sg, model, 1.0, true, {}, settings, &first, logger); + CHECK(second.warm); + CHECK(second.to_model_setting == first.to_model_setting); + CHECK(second.cell.a == Catch::Approx(34.0)); + REQUIRE(second.validation); + CHECK(second.validation->r_work < 0.15); + + std::filesystem::remove(data_model); + std::filesystem::remove(path); +} + +// The model step of a live serial merge (ModelStep.h). A merge in the other indexing of a P4 crystal is put in the +// model's by the probe; a later merge of the same collection, again in the other indexing, is put in the earlier +// step's frame and placed from where the model was, with no null or probe made again; a cut leaves out what has not +// started. The group a merge without one is merged in is the model's where its cell is the run's. +TEST_CASE("ModelStep_PlacesWarmStartsAndCuts", "[ModelValidation][RugnuxLive]") { + Logger logger("ModelStep_PlacesWarmStartsAndCuts"); + const auto path = WriteTemp("stills_model_test.pdb", + ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n").c_str()); + const RugnuxModel model = ReadModel(path); + auto own = ModelReferenceIntensities(model, {}, {}, 2.5, logger); + REQUIRE(own.size() > 500); + for (size_t i = 0; i < own.size(); i++) { + own[i].F = std::sqrt(std::max(0.0f, own[i].I)); + own[i].sigma = 0.02f * own[i].I + 1.0f; + own[i].sigmaF = 0.02f * own[i].F + 0.1f; + own[i].rfree_flag = (i % 20) == 0; + } + const UnitCell cell{.a = 34, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup p4 = *gemmi::find_spacegroup_by_name("P 4"); + const auto laws = ReindexAmbiguityOperators(cell, p4); + REQUIRE(laws.size() == 1); + const auto other = ReindexMergedIntoAsu(own, laws.front(), p4, true); + + // The group: the collection's, else the model's on the run's cell, else none + CHECK(ModelStepSpaceGroup(*gemmi::find_spacegroup_by_name("P 4 2 2"), cell, model).space_group->xhm() == "P 4 2 2"); + CHECK(ModelStepSpaceGroup(std::nullopt, cell, model).space_group->xhm() == "P 4"); + const UnitCell longer{.a = 34, .b = 34, .c = 45, .alpha = 90, .beta = 90, .gamma = 90}; + const auto none = ModelStepSpaceGroup(std::nullopt, longer, model); + CHECK_FALSE(none.space_group); + CHECK_FALSE(none.reason.empty()); + + ModelStepSettings settings; + settings.refine = false; + const ModelStepResult first = RunModelStep(other, cell, p4, model, 1.0, true, {}, settings, nullptr, + logger); + REQUIRE(first.validation); + CHECK(first.placement == "RIGID_BODY"); + CHECK_FALSE(first.warm); + CHECK(first.validation->indexing_op == laws.front()); + REQUIRE(first.maps); + REQUIRE(first.placed_model); + // In the model's indexing: the model's own intensities + CHECK(ReferenceIntensityCC(first.merged, own, p4) > 0.99); + CHECK(first.validation->r_work < 0.15); + + // A later merge, in the other indexing again, and refined + settings.refine = true; + const ModelStepResult second = RunModelStep(other, cell, p4, model, 1.0, true, {}, settings, &first, + logger); + REQUIRE(second.validation); + CHECK(second.warm); + CHECK(second.placement == "RIGID_BODY"); + CHECK(second.validation->model_fits); + CHECK(ReferenceIntensityCC(second.merged, own, p4) > 0.99); + REQUIRE(second.refinement); + CHECK(second.refinement->ok); + CHECK(second.refinement->r_free < 0.2); + + // Cut before anything: the refinement is left out; and for data the model does not describe, the molecular + // replacement too + std::atomic cut{true}; + settings.cut = &cut; + const ModelStepResult cut_placed = RunModelStep(other, cell, p4, model, 1.0, true, {}, settings, + &first, logger); + CHECK_FALSE(cut_placed.refinement); + CHECK(cut_placed.left_out == std::vector{"refinement"}); + const auto unrelated_path = WriteTemp("stills_model_unrelated_test.pdb", + ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n", 7).c_str()); + auto unrelated = ModelReferenceIntensities(ReadModel(unrelated_path), {}, {}, 2.5, logger); + for (size_t i = 0; i < unrelated.size(); i++) { + unrelated[i].F = std::sqrt(std::max(0.0f, unrelated[i].I)); + unrelated[i].sigma = 0.02f * unrelated[i].I + 1.0f; + unrelated[i].sigmaF = 0.02f * unrelated[i].F + 0.1f; + unrelated[i].rfree_flag = (i % 20) == 0; + } + const ModelStepResult cut_unplaced = RunModelStep(unrelated, cell, p4, model, 1.0, true, {}, + settings, nullptr, logger); + CHECK(cut_unplaced.placement == "NOT_PLACED"); + CHECK_FALSE(cut_unplaced.mr); + CHECK(cut_unplaced.left_out == std::vector{"molecular replacement"}); + CHECK_FALSE(cut_unplaced.placed_model); + + std::filesystem::remove(path); + std::filesystem::remove(unrelated_path); +}