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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-11-11 13:39:47 +01:00
1 parent 6df0bd3a6a
commit 3ba8ba1420
8 files changed
+526 -45

No files matched your search

+4
View File
@@ -12,6 +12,10 @@
#include <ctime>
#include <cstdio>
inline double SecondsSince(std::chrono::steady_clock::time_point t0) {
return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
}
inline std::string time_UTC(const std::chrono::time_point<std::chrono::system_clock> &input) {
auto time_ms = std::chrono::duration_cast<std::chrono::milliseconds>(input.time_since_epoch()).count();
+2
View File
@@ -25,6 +25,8 @@ ADD_LIBRARY(JFJochRugnux STATIC
RugnuxLive.h
StillsLiveMerge.cpp
StillsLiveMerge.h
ModelStep.cpp
ModelStep.h
RugnuxRefine.cpp
RugnuxRefine.h
Rugnux.h
+256
View File
@@ -0,0 +1,256 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ModelStep.h"
#include <chrono>
#include <cmath>
#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<MergedReflection> &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<MergedReflection> &merged, const UnitCell &cell,
const gemmi::SpaceGroup &space_group, const RugnuxModel &model,
const gemmi::Structure &placed, double wavelength_A,
const std::vector<float> &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<MergedReflection> &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<gemmi::SpaceGroup> &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<MergedReflection> &merged, const UnitCell &cell,
const gemmi::SpaceGroup &space_group, const RugnuxModel &model,
double wavelength_A, bool merge_friedel,
const std::vector<float> &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<MergedReflection> 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<MergedReflection> &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;
}
+101
View File
@@ -0,0 +1,101 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <atomic>
#include <functional>
#include <optional>
#include <string>
#include <vector>
#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<MergedReflection> &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<MergedReflection> &merged, const UnitCell &cell,
const gemmi::SpaceGroup &space_group, const RugnuxModel &model,
const gemmi::Structure &placed, double wavelength_A,
const std::vector<float> &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<MergedReflection> &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<gemmi::SpaceGroup> space_group;
std::string reason;
};
ModelStepGroup ModelStepSpaceGroup(const std::optional<gemmi::SpaceGroup> &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<bool> *cut = nullptr;
// Called as each part starts: "Model validation", "Molecular replacement", "Refinement"
std::function<void(const std::string &)> 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<MergedReflection> 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<ModelValidationResult> validation; // its maps moved to `maps`
std::optional<MolecularReplacementResult> mr;
// Where the model was placed: the maps of the placement and the model as placed, in the merge's frame
std::shared_ptr<ModelMaps> maps;
std::optional<PlacedModelOutput> placed_model;
std::optional<RefinementResult> refinement; // its model, maps and blobs where it ran
std::vector<std::string> 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<MergedReflection> &merged, const UnitCell &cell,
const gemmi::SpaceGroup &space_group, const RugnuxModel &model,
double wavelength_A, bool merge_friedel,
const std::vector<float> &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);
+11 -41
View File
@@ -3,6 +3,7 @@
#include <cstdlib>
#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<size_t>(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<size_t>(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<size_t>(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<MergedReflection> &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";
+5 -4
View File
@@ -7,6 +7,7 @@
#include <ctime>
#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<double>(std::chrono::steady_clock::now() - t0).count();
}
size_t ImagesWithReflections(const std::vector<IntegrationOutcome> &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())
+2
View File
@@ -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
+145
View File
@@ -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<bool> 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<std::string>{"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<std::string>{"molecular replacement"});
CHECK_FALSE(cut_unplaced.placed_model);
std::filesystem::remove(path);
std::filesystem::remove(unrelated_path);
}