mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 84 Wiki Activity

Merge broker-rugnux-live into rc176: rugnux in jfjoch_broker as a live thread during a DECTRIS rotation collection (off by default), the rugnux library without file I/O (results in ProcessResult::outputs, the CLI writes them through JFJochRugnuxOutput), END carrying the merged table, map coefficients and placed model

Browse Source 
Branch reviewed at e16db07ff (review: owner's decisions accepted - _maps.mtz instead of .ccp4 in the
viewer, -o "" computes everything incl. MR/validation, Ctrl+C/error paths may leave fewer files, the
writer links gemmi/zlib) and merged at 55184eb09 with its review fixes (END size guard, frame drop as
in Lite, outputs released before re-runs, pool priority, docs, frontend gating).

Conflicts: rugnux/Rugnux.h (kept rc176's host_memory_budget_bytes, output_prefix removed as on the
branch), rugnux/RugnuxPipeline.h (MtzTable future + rc176's GPUMerge pointer), MolecularReplacement.h/.cpp
(rc176's GPU cost model, packed transform and split functions with the branch's RugnuxModel signature;
search_model with rc176's other_budgets, result.model_name). tests/RugnuxLargeTest.cpp: rc176's
host-budget test writes its files through WriteRugnuxOutputs now that the library writes none.
Generated code (broker/gen, frontend client, Python client docs, Redoc) regenerated from the merged
jfjoch_api.yaml with update_version.sh (version unchanged).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-17 21:01:38 +02:00
2 parents 3e36fe9468 55184eb096
commit 78c2a69267
130 files changed
+7612 -1010
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/structure_refinement/ModelValidation.cpp
+91 -121
View File
@@ -21,16 +21,12 @@
#include <vector>
#include <unordered_map>
#include <gemmi/mmread_gz.hpp> // read_structure_gz
#include <gemmi/assembly.hpp> // expand_ncs
#include <gemmi/gz.hpp> // MaybeGzipped
#include <gemmi/it92.hpp> // IT92 x-ray form factors
#include <gemmi/dencalc.hpp> // DensityCalculator
#include <gemmi/fourier.hpp> // get_size_for_hkl, get_f_phi_on_grid (the maps)
#include <gemmi/solmask.hpp> // SolventMasker
#include <gemmi/scaling.hpp> // Scaling (bulk solvent + anisotropic B)
#include <gemmi/ccp4.hpp> // Ccp4 map I/O
#include <gemmi/mtz.hpp> // Mtz (map-coefficient output)
#include <gemmi/stats.hpp> // calculate_data_statistics
#include "../../common/CorrelationCoefficient.h"
@@ -82,22 +78,9 @@ gemmi::Grid<float> map_from_coefficients(gemmi::AsuData<std::complex<float>> &co
return MapFromFPhi(gemmi::get_f_phi_on_grid<float>(coef, size, true));
}
// A file this run may already have written - the validation in the model's setting writes over the
// first one's maps - is removed before it is written again, not truncated: XFS (and ext4) flush a file
// truncated and rewritten when it is closed, and three 150 MB maps forced out to a disk that way cost
// seconds where a fresh file costs nothing.
void remove_before_rewriting(const std::string &path) {
std::remove(path.c_str());
}
// Write a map as CCP4; return its RMS (the sigma the map is read in).
double write_ccp4(const gemmi::Grid<float> &map, const std::string &path) {
gemmi::Ccp4<float> ccp4;
ccp4.grid = map;
ccp4.update_ccp4_header(2);
remove_before_rewriting(path);
ccp4.write_ccp4_map(path);
return ccp4.hstats.rms;
// The RMS of a map, the sigma it is read in: the one its CCP4 header states (Ccp4::update_ccp4_header).
double map_rms(const gemmi::Grid<float> &map) {
return gemmi::calculate_data_statistics(map.data).rms;
}
// Cubic, not the default linear, for reading a map at a point. The maps are sampled every d_min/3,
@@ -538,8 +521,7 @@ namespace {
// validation either way.
ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
const UnitCell &cell,
const std::string &model_path,
const std::string &output_prefix,
const RugnuxModel &model,
Logger &logger,
const gemmi::SpaceGroup *data_space_group,
bool probe_indexing_ambiguity,
@@ -549,32 +531,29 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law,
bool rigid_body_gpu,
const gemmi::Structure *model) {
const gemmi::Structure *placed) {
ModelValidationResult result;
result.model_path = model_path;
result.model_name = model.name;
// --- read the atomic model ---
// --- the atomic model ---
ModelState mdl;
gemmi::Structure &st = mdl.st; // the real model, the one the maps and the report describe
try {
// Detect, not the default: without it GEMMI picks the format from the extension and only
// falls back to the content when it does not recognise one. A model arrives named however
// whoever produced it named it, so the file itself is the better authority.
st = model != nullptr ? *model : gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
st = placed != nullptr ? *placed : model.structure;
expand_strict_ncs(st);
} catch (const std::exception &e) {
result.failure_reason = fmt::format("cannot read model {}: {}", model_path, e.what());
result.failure_reason = fmt::format("cannot read model {}: {}", model.name, e.what());
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
if (st.models.empty() || !st.cell.is_crystal()) {
result.failure_reason = fmt::format("model {} has no atoms or no unit cell", model_path);
result.failure_reason = fmt::format("model {} has no atoms or no unit cell", model.name);
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
const gemmi::SpaceGroup *sg = st.find_spacegroup();
if (!sg) {
result.failure_reason = fmt::format("model {} has no usable space group", model_path);
result.failure_reason = fmt::format("model {} has no usable space group", model.name);
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
@@ -660,7 +639,7 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
// and their maps and reflection files can be compared directly. Moving the model would leave
// each dataset in whichever indexing its own run happened to pick.
// A model handed over already placed in the data's cell (molecular replacement) is on the data's axes.
if (data_cell.is_crystal() && model == nullptr) {
if (data_cell.is_crystal() && placed == nullptr) {
const gemmi::GroupOps model_gops = sg->operations();
const bool same_point_group =
data_space_group != nullptr && same_rotations(data_space_group->operations(), model_gops);
@@ -1734,23 +1713,27 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
mapfofc.v.push_back({terms[i].hkl, delfwt[i] * ph});
}
// A validation started before the one it follows had decided writes nothing until that one has,
// and nothing at all where it decided otherwise.
if (schedule.write_gate.valid() && !schedule.write_gate.get()) {
result.failure_reason = "superseded before its files were written";
// A validation started before the one it follows had decided makes no maps until that one has,
// and none at all where it decided otherwise.
if (schedule.maps_gate.valid() && !schedule.maps_gate.get()) {
result.failure_reason = "superseded before its maps were made";
return result;
}
// Where the validation in the model's setting, started beside this one, is going to be kept, it
// writes the maps and the map coefficients on the axes the files are written on, and this one's would
// be overwritten unread: they are not made at all.
// Where the validation in the model's setting, started beside this one, is going to be kept, its
// maps and map coefficients are on the axes the files are written on, and this one's would be
// thrown away unread: they are not made at all.
const bool maps_superseded = schedule.maps_superseded && schedule.maps_superseded(result);
std::shared_ptr<ModelMaps> maps;
if (!maps_superseded) {
// --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) ---
maps = std::make_shared<ModelMaps>();
// --- make the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) ---
// The three maps - this one, the difference map and the anomalous map - share nothing but what
// they read, so the other two are made and written on threads of their own beside this one.
// they read, so the other two are made on threads of their own beside this one. Nothing here
// reads the difference map, so it is made only where it is kept.
std::future<void> fofc_map = std::async(std::launch::async, [&] {
write_ccp4(map_from_coefficients(mapfofc, sf_gpu), output_prefix + "_fofc.ccp4");
if (schedule.keep_map_grids)
maps->map_fofc = map_from_coefficients(mapfofc, sf_gpu);
});
// --- anomalous difference map, where the merge kept the Bijvoet split ---
@@ -1764,13 +1747,15 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
result.indexing_op, fmodel);
result.anomalous_pairs = static_cast<int>(mapanom.v.size());
if (!mapanom.v.empty()) {
const gemmi::Grid<float> grid = map_from_coefficients(mapanom, sf_gpu);
read_anomalous_map(grid, write_ccp4(grid, output_prefix + "_anom.ccp4"), st, result);
gemmi::Grid<float> grid = map_from_coefficients(mapanom, sf_gpu);
read_anomalous_map(grid, map_rms(grid), st, result);
if (schedule.keep_map_grids)
maps->map_anom = std::move(grid);
}
});
const gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc, sf_gpu);
const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4");
gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc, sf_gpu);
const double rms2 = map_rms(grid2fofc);
{
double s = 0; int n = 0;
for (gemmi::Model &m : st.models)
@@ -1779,6 +1764,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; }
result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
}
if (schedule.keep_map_grids)
maps->map_2fofc = std::move(grid2fofc);
fofc_map.get();
anomalous_map.get();
}
@@ -1803,60 +1790,53 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
: "the model does not fit these data");
if (!maps_superseded) {
// --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. ---
try {
gemmi::Mtz mtz(true);
// The group the REFLECTIONS end up in, which is the model's only where its hand was adopted -
// AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes
// identically, so the coefficients are the same numbers either way and only the label moves;
// but the two groups have different screw translations, so a reader that expands symmetry out
// of this file works in the wrong one if the label disagrees with the .mtz beside it.
mtz.spacegroup = result.adopted_model_enantiomorph
? sg
// No data group means the caller merged in P1, and P1 is what the
// reflection files beside this one carry - not the model's group.
: (data_space_group != nullptr ? data_space_group
: gemmi::find_spacegroup_by_number(1));
mtz.set_cell_for_all(ucell);
mtz.add_dataset("model_validation");
mtz.datasets.back().wavelength = wavelength_A;
mtz.add_column("FP", 'F', -1, -1, false);
mtz.add_column("FC", 'F', -1, -1, false);
mtz.add_column("PHIC", 'P', -1, -1, false);
mtz.add_column("FWT", 'F', -1, -1, false);
mtz.add_column("PHWT", 'P', -1, -1, false);
mtz.add_column("DELFWT", 'F', -1, -1, false);
mtz.add_column("PHDELWT", 'P', -1, -1, false);
// The figure of merit the coefficients carry, so the weighting can be read off - and undone -
// from the file rather than having to be taken on trust.
mtz.add_column("FOM", 'W', -1, -1, false);
mtz.add_column("FREE", 'I', -1, -1, false);
std::vector<float> data;
for (size_t i = 0; i < terms.size(); i++) {
// [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg
// returns (-180, 180] and a P column is not supposed to.
const double phi_wrapped = terms[i].phi * 180.0 / PI;
const auto phi_deg = static_cast<float>(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped);
data.insert(data.end(), {static_cast<float>(terms[i].hkl[0]), static_cast<float>(terms[i].hkl[1]),
static_cast<float>(terms[i].hkl[2]),
static_cast<float>(terms[i].fo), static_cast<float>(terms[i].fc),
phi_deg,
fwt[i], phi_deg,
delfwt[i], phi_deg,
static_cast<float>(sigma_a.weight[i].m),
terms[i].free ? 0.0f : 1.0f});
}
mtz.nreflections = static_cast<int>(terms.size());
mtz.data = std::move(data);
remove_before_rewriting(output_prefix + "_maps.mtz");
mtz.write_to_file(output_prefix + "_maps.mtz");
} catch (const std::exception &e) {
logger.Warning("Model validation: could not write map MTZ: {}", e.what());
// --- the map coefficients, so the maps can be re-opened / rebuilt in Coot etc. ---
MtzTable &mtz = maps->coefficients;
// The group the REFLECTIONS end up in, which is the model's only where its hand was adopted -
// AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes
// identically, so the coefficients are the same numbers either way and only the label moves;
// but the two groups have different screw translations, so a reader that expands symmetry out
// of this file works in the wrong one if the label disagrees with the .mtz beside it.
mtz.space_group_name = result.adopted_model_enantiomorph
? sg->xhm()
// No data group means the caller merged in P1, and P1 is what the
// reflection files beside this one carry - not the model's group.
: (data_space_group != nullptr ? data_space_group->xhm()
: gemmi::find_spacegroup_by_number(1)->xhm());
mtz.unit_cell = UnitCell{static_cast<float>(ucell.a), static_cast<float>(ucell.b),
static_cast<float>(ucell.c), static_cast<float>(ucell.alpha),
static_cast<float>(ucell.beta), static_cast<float>(ucell.gamma)};
mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f},
{"model_validation", "model_validation", "model_validation",
static_cast<float>(wavelength_A)}};
// The figure of merit the coefficients carry, so the weighting can be read off - and undone -
// from the file rather than having to be taken on trust.
const std::pair<const char *, char> columns[] = {{"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"FP", 'F'},
{"FC", 'F'}, {"PHIC", 'P'}, {"FWT", 'F'}, {"PHWT", 'P'},
{"DELFWT", 'F'}, {"PHDELWT", 'P'}, {"FOM", 'W'},
{"FREE", 'I'}};
for (const auto &[label, type] : columns)
mtz.columns.push_back({label, type, type == 'H' ? 0 : 1, {}});
for (size_t i = 0; i < terms.size(); i++) {
// [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg
// returns (-180, 180] and a P column is not supposed to.
const double phi_wrapped = terms[i].phi * 180.0 / PI;
const auto phi_deg = static_cast<float>(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped);
const float row[] = {static_cast<float>(terms[i].hkl[0]), static_cast<float>(terms[i].hkl[1]),
static_cast<float>(terms[i].hkl[2]),
static_cast<float>(terms[i].fo), static_cast<float>(terms[i].fc),
phi_deg,
fwt[i], phi_deg,
delfwt[i], phi_deg,
static_cast<float>(sigma_a.weight[i].m),
terms[i].free ? 0.0f : 1.0f};
for (size_t c = 0; c < mtz.columns.size(); c++)
mtz.columns[c].values.push_back(row[c]);
}
}
result.ok = true;
result.maps_prefix = output_prefix;
result.maps = std::move(maps);
// The placed coordinates, for the caller to write out beside the maps once the frame is settled.
result.placed_model = std::make_shared<gemmi::Structure>(st);
logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) "
@@ -1885,12 +1865,7 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
"data were indexed in the wrong hand, and reindexing would hide which",
result.anomalous_deepest_site, result.anomalous_deepest_sigma);
if (maps_superseded)
logger.Info("Model validation: no maps made here - the validation in the model's setting writes them");
else
logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz",
output_prefix, output_prefix,
result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4",
output_prefix);
logger.Info("Model validation: no maps made here - the validation in the model's setting makes them");
if (result.fit_tested && !result.model_fits)
logger.Info("Model validation: the model was rejected, so nothing downstream moved - the "
"reflection files are byte for byte the ones a run with no model would have "
@@ -1975,8 +1950,7 @@ std::vector<size_t> OtherCardsTotalMemory(int device) {
ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &merged,
const UnitCell &cell,
const std::string &model_path,
const std::string &output_prefix,
const RugnuxModel &model,
Logger &logger,
const gemmi::SpaceGroup *data_space_group,
bool probe_indexing_ambiguity,
@@ -1985,18 +1959,18 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
const std::vector<float> &report_shell_d_min,
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law,
const gemmi::Structure *model) {
const gemmi::Structure *placed) {
#ifdef JFJOCH_USE_CUDA
// A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through:
// the validation reports why it did not finish, and a validation that did not finish decides nothing,
// so the reflection files are those of a run without a model.
try {
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group,
return Validate(merged, cell, model, logger, data_space_group,
probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law,
true, model);
true, placed);
} catch (const RigidBodyGPUFailure &e) {
ModelValidationResult failed;
failed.model_path = model_path;
failed.model_name = model.name;
failed.failure_reason = fmt::format("the GPU failed during the validation ({})", e.what());
logger.Error("Model validation: {}", failed.failure_reason);
cuda_throw_if_context_lost(); // a lost device ends the run here, not at its next user
@@ -2004,8 +1978,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
return failed;
}
#else
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity,
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, model);
return Validate(merged, cell, model, logger, data_space_group, probe_indexing_ambiguity,
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, placed);
#endif
}
@@ -2161,7 +2135,7 @@ void KeepModelVerdict(ModelValidationResult &remade, const ModelValidationResult
remade.model_space_group_number = decided.model_space_group_number;
}
std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model_path,
std::vector<MergedReflection> ModelReferenceIntensities(const RugnuxModel &model,
const std::optional<UnitCell> &cell,
const gemmi::SpaceGroup *space_group,
double d_min,
@@ -2172,19 +2146,15 @@ std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model
return out;
}
gemmi::Structure st;
gemmi::Structure st = model.structure;
try {
// Detect, not the default: without it GEMMI picks the format from the extension and only
// falls back to the content when it does not recognise one. A model arrives named however
// whoever produced it named it, so the file itself is the better authority.
st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
expand_strict_ncs(st);
} catch (const std::exception &e) {
logger.Error("Model reference: cannot read model {}: {}", model_path, e.what());
logger.Error("Model reference: cannot read model {}: {}", model.name, e.what());
return out;
}
if (st.models.empty() || !st.cell.is_crystal()) {
logger.Error("Model reference: model {} has no atoms or no unit cell", model_path);
logger.Error("Model reference: model {} has no atoms or no unit cell", model.name);
return out;
}
@@ -2207,7 +2177,7 @@ std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model
st.spacegroup_hm = space_group->xhm();
const gemmi::SpaceGroup *sg = st.find_spacegroup();
if (!sg) {
logger.Error("Model reference: model {} has no usable space group", model_path);
logger.Error("Model reference: model {} has no usable space group", model.name);
return out;
}
st.setup_cell_images();
@@ -2257,6 +2227,6 @@ std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model
}
logger.Info("Model reference: {} intensities computed from {} to {:.2f} A, space group {}",
out.size(), model_path, d_min, sg->short_name());
out.size(), model.name, d_min, sg->short_name());
return out;
}