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

ModelValidation: score the null, and the real fit against it, to 3.5 A [not exact - changes sigmas]

Browse Source 
QUALITY-AFFECTING, separate from the exact commits before it. The nine null replicates were fitted
at the data's full resolution: a full-resolution Fcalc and mask, 1 + n(indexing candidates) scale
fits over every reflection, the rigid body (which already stops at 3.5 A), then another Fcalc and
fit. Those full-resolution scale fits were most of the null's CPU. Now the replicates, the real
model's R-work and the real indexing margin they are compared with are all computed to
max(d_min, 3.5 A) - the same fits on the same reflections on both sides - but never on fewer than
1000 free reflections, since the margin is read in R-free: on a small cell the limit moves out to
where 1000 of them are, or to the data's own. The reported R-factors, maps, placement and files are
still made at full resolution; only the null test moves.

Why the free-reflection floor: without it the [ModelValidation] reindexing test (a 34x34x38 A cell,
1553 reflections) lost its decision - at 3.5 A its margin was read on ~30 free reflections, the
null margin sd went to 0.053 and the real margin to +2.6 sigma, under the gate of 3.

Measured against the floor-fixed build before this commit (battery commands, --model; CPU seconds
in the null - wall clock on this shared box moved with other jobs, 2-4x less where it was quiet):
  set   CPU s        R-work sigma     indexing margin sigma      verdict / decision / files
  9hnc  2307 ->  583  +10.4 -> +15.3   +106 -> +21 (reindexed)   same / same / identical
  8sa8  1382 ->  313  +328  -> +203    -                          same / same / identical
  8xtg   968 ->  496  +215  -> +172    +77 -> +27 (kept)          same / same / identical
  9fhc  1042 ->  737  +201  -> +139    +76 -> +41 (reindexed)     same / same / identical
  9bn8   393 ->  153  +378  -> +235    +137 -> +56 (reindexed)    same / same / identical
  5epe   631 ->  600  +0.01 -> -0.03   +26 -> +25 (not decided)   same / same / identical
Sets where the null does not run (6oel, 5reo) are untouched by construction. The sigmas drop
because the null's sd at 3.5 A is about twice as wide; the smallest margin here is still 8x the gate,
but a borderline indexing decision would be decided less often, not more.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_f
2026-09-27 23:12:27 +02:00
co-authored by Claude Opus 5.5
parent 7eaf0f21b0
commit 3783d388ed
1 changed files with 83 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files
rugnux/ModelValidation.cpp
+83 -30
View File
@@ -9,6 +9,7 @@
#include <complex>
#include <array>
#include <atomic>
#include <functional>
#include <future>
#include <numeric>
#include <optional>
@@ -99,6 +100,11 @@ constexpr unsigned NULL_SEED = 20260902;
// measured in.
constexpr double MODEL_FIT_SIGMA = 3.0;
// The resolution the null and the real fit's side of it are scored to (see the null below), and the
// fewest free reflections it may leave the indexing margin, which is read in R-free.
constexpr double NULL_D_MIN = 3.5;
constexpr size_t NULL_MIN_FREE = 1000;
// Mean and sample standard deviation of a small sample.
std::pair<double, double> mean_sd(const std::vector<double> &v) {
if (v.size() < 2)
@@ -486,9 +492,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
// --- Fcalc (atomic) via electron density on a grid + FFT, plus a flat bulk-solvent mask -> Fmask.
// A lambda because the rigid-body step below moves the model and then needs both again, and it
// takes the state to work on so that a null replicate can run it on its own copy. ---
auto compute_model_factors = [&](ModelState &ms) {
auto compute_model_factors = [&](ModelState &ms, double to_d) {
gemmi::DensityCalculator<Table, float> dc;
dc.d_min = d_min;
dc.d_min = to_d;
dc.rate = 1.5;
dc.set_grid_cell_and_spacegroup(ms.st);
dc.set_refmac_compatible_blur(ms.st.models[0]);
@@ -505,7 +511,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
masker.put_mask_on_grid(mask_grid, ms.st.models[0]);
ms.fmask = MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0);
};
compute_model_factors(mdl);
compute_model_factors(mdl, d_min);
gemmi::GroupOps gops = sg->operations();
gemmi::ReciprocalAsu asu(sg);
@@ -611,7 +617,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
double rb_angle_deg = 0, rb_shift_A = 0, r_free_before_rb = 0;
};
auto place_and_fit = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in,
std::optional<Fit> fitted) -> Placement {
std::optional<Fit> fitted, double to_d) -> Placement {
Placement out;
out.fit = fitted ? std::move(*fitted) : fit_model(ms, obs_in);
@@ -625,7 +631,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
// exactly where it was read.
const std::vector<gemmi::Position> before = ModelPositions(ms.st.models[0]);
const RigidBodyRefineResult rb =
RefineRigidBody(ms.st.models[0], ucell, *sg, out.fit.fobs_work, d_min, logger, nthreads);
RefineRigidBody(ms.st.models[0], ucell, *sg, out.fit.fobs_work, to_d, logger, nthreads);
if (!rb.converged) {
// RefineRigidBody leaves the model wherever the solver left it, usable answer or not, so
// the restore cannot be conditional on the same flag the re-fit is.
@@ -633,7 +639,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
return out;
}
{
compute_model_factors(ms);
compute_model_factors(ms, to_d);
Fit moved = fit_model(ms, obs_in);
const bool commit = moved.r_free < out.fit.r_free;
logger.Info("Model validation: rigid body held-out R-free {:.4f} -> {:.4f} => {}",
@@ -657,13 +663,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
double margin = 0; // by how much in R-free it leads the runner-up
std::optional<Fit> identity_fit; // fit_model(ms, obs), for place_and_fit to reuse
};
auto probe_indexing = [&](ModelState &ms) {
auto probe_indexing = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in) {
IndexingProbe out;
out.identity_fit = fit_model(ms, obs);
out.identity_fit = fit_model(ms, obs_in);
std::vector<double> r_free{out.identity_fit->r_free}; // identity first, then the twin laws
double best_r_free = r_free.front();
for (const auto &op : reindex_ops) {
const double cand = fit_model(ms, ReindexReflections(obs, op)).r_free;
const double cand = fit_model(ms, ReindexReflections(obs_in, op)).r_free;
r_free.push_back(cand);
if (cand < best_r_free) {
best_r_free = cand;
@@ -683,7 +689,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
std::vector<MergedReflection> reindexed;
const std::vector<MergedReflection> *obs_model = &obs;
if (!reindex_ops.empty()) {
indexing = probe_indexing(mdl);
indexing = probe_indexing(mdl, obs);
if (!(indexing.op == gemmi::Op::identity())) {
reindexed = ReindexReflections(obs, indexing.op);
obs_model = &reindexed;
@@ -695,7 +701,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
if (obs_model == &obs)
identity_fit = std::move(indexing.identity_fit);
indexing.identity_fit.reset();
Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit));
Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit), d_min);
// --- is there anything for the model to decide? ---
// There are three: the space-group label, where the model asserts the other enantiomorph; the
@@ -746,6 +752,52 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
for (int i = 0; i < NULL_REPLICATES; i++)
null_rotation.push_back(random_rotation(rng));
}
// The null and the side of the real fit it is compared with are both scored to NULL_D_MIN,
// not to the data's resolution. The placement already stops there (RigidBodyRefine's ladder
// ends at 3.5 A), and whether a model is in the right orientation, or the data in its
// indexing, is decided at low resolution: the finer shells only add reflections to the scale
// fits, which were most of the null's cost. Both sides go through the same fits on the same
// reflections, so the comparison stays like for like. On a small cell 3.5 A leaves too few
// free reflections for a real indexing margin to be told from a random placement's, so the
// limit then moves out to where NULL_MIN_FREE of them are - to the data's own on the
// smallest, which are the cheap ones.
std::vector<float> free_d;
for (const MergedReflection &r : obs)
if (r.rfree_flag && !std::isnan(r.F))
free_d.push_back(r.d);
std::sort(free_d.begin(), free_d.end(), std::greater<>());
double null_d_min = std::max(d_min, NULL_D_MIN);
if (free_d.size() < NULL_MIN_FREE)
null_d_min = d_min;
else
null_d_min = std::min<double>(null_d_min, free_d[NULL_MIN_FREE - 1]);
std::vector<MergedReflection> obs_null;
for (const MergedReflection &r : obs)
if (r.d >= null_d_min)
obs_null.push_back(r);
// The real model's side of it: the indexing margin from the model as read, as a replicate's
// is, and the R-work of the model where its placement left it, against the reflections that
// placement was fitted to.
auto coarse_r_work = [&](const std::vector<MergedReflection> &placed_against) {
ModelState coarse;
coarse.st = mdl.st;
compute_model_factors(coarse, null_d_min);
std::vector<MergedReflection> in;
for (const MergedReflection &r : placed_against)
if (r.d >= null_d_min)
in.push_back(r);
return fit_model(coarse, in).r_work;
};
double real_margin = 0;
if (!reindex_ops.empty()) {
ModelState coarse;
coarse.st = mdl.st;
SetModelPositions(coarse.st.models[0], as_read);
compute_model_factors(coarse, null_d_min);
real_margin = probe_indexing(coarse, obs_null).margin;
}
double real_r_work = coarse_r_work(*obs_model);
// Each replicate works on its own copy of the model and of its structure factors, and the
// replicates share nothing else, so they run concurrently and the real model - which the
// maps and the atom-density readout below are still taken from - is never moved.
@@ -755,14 +807,14 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
rep.st = mdl.st;
SetModelPositions(rep.st.models[0], as_read);
reorient_about_centroid(rep.st.models[0], null_rotation[i]);
compute_model_factors(rep);
compute_model_factors(rep, null_d_min);
std::optional<Fit> identity_fit;
if (!reindex_ops.empty()) {
IndexingProbe probe = probe_indexing(rep);
IndexingProbe probe = probe_indexing(rep, obs_null);
null_margin[i] = probe.margin;
identity_fit = std::move(probe.identity_fit); // replicates are placed against obs
}
null_r_work[i] = place_and_fit(rep, obs, std::move(identity_fit)).fit.r_work;
null_r_work[i] = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min).fit.r_work;
};
// On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool
// worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of
@@ -785,7 +837,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
result.null_r_work_mean = null_mean;
result.null_r_work_sd = null_sd;
auto record_fit = [&](const Placement &p) {
auto record_fit = [&](double r_work) {
// The spread is floored, not trusted down to zero: sd == 0 and sd == 1e-7 are one ulp apart
// and would give sigmas of nothing and of millions. A null with no spread - a model too
// small or too symmetric for a rotation about its own centroid to move |Fcalc| - fits no
@@ -793,10 +845,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
// not a gate: a large, well-measured crystal has a null as narrow as 0.0006-0.0010
// (measured), and a real fit 0.3 below it is hundreds of sigma, not none.
constexpr double NULL_SD_FLOOR = 1e-3;
result.r_work_sigma = (null_mean - p.fit.r_work) / std::max(null_sd, NULL_SD_FLOOR);
result.r_work_sigma = (null_mean - r_work) / std::max(null_sd, NULL_SD_FLOOR);
result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA;
};
record_fit(real);
record_fit(real_r_work);
// --- does the model get to decide anything? ---
// The R-factors, the maps and the rigid-body placement are statements about the MODEL and are
@@ -806,12 +858,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
// wrong model must not be able to make.
if (!reindex_ops.empty()) {
result.indexing_probed = true;
result.indexing_margin = indexing.margin;
result.indexing_margin = real_margin;
const auto [margin_mean, margin_sd] = mean_sd(null_margin);
result.indexing_margin_null_mean = margin_mean;
result.indexing_margin_null_sd = margin_sd;
constexpr double MARGIN_SD_FLOOR = 1e-4;
result.indexing_margin_sigma = (indexing.margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR);
result.indexing_margin_sigma = (real_margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR);
// The margin is judged against the margin a random placement produces, not against a value.
// A random model also picks a winner, and on these data it picks one by a comparable lead
// (measured), so the raw margin says nothing on its own.
@@ -820,9 +872,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
if (decided)
result.indexing_op = indexing.op;
logger.Info("Model validation: probed {} indexing solution(s) against the model; the winner "
"leads the runner-up by {:.4f} in R-free, against {:.4f} +- {:.4f} for a random "
"placement of the same model ({:+.2f} sigma) => {}",
reindex_ops.size() + 1, indexing.margin, margin_mean, margin_sd,
"leads the runner-up by {:.4f} in R-free to {:.2f} A, against {:.4f} +- {:.4f} for a "
"random placement of the same model ({:+.2f} sigma) => {}",
reindex_ops.size() + 1, real_margin, null_d_min, margin_mean, margin_sd,
result.indexing_margin_sigma,
!decided ? "not decided; the data keep the indexing they were merged in"
: (result.indexing_op == gemmi::Op::identity()
@@ -831,14 +883,15 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
// The R-factors and the maps have to describe the reflections the files carry, and those
// are now the ones the merge produced, so the fit is remade on them.
SetModelPositions(st.models[0], as_read);
compute_model_factors(mdl);
real = place_and_fit(mdl, obs, std::nullopt);
record_fit(real);
compute_model_factors(mdl, d_min);
real = place_and_fit(mdl, obs, std::nullopt, d_min);
real_r_work = coarse_r_work(obs);
record_fit(real_r_work);
}
}
logger.Info("Model validation: R-work {:.4f} against a null of {:.4f} +- {:.4f} = {:+.2f} sigma "
"=> the model {} these data",
real.fit.r_work, null_mean, null_sd, result.r_work_sigma,
logger.Info("Model validation: R-work {:.4f} to {:.2f} A against a null of {:.4f} +- {:.4f} = "
"{:+.2f} sigma => the model {} these data",
real_r_work, null_d_min, null_mean, null_sd, result.r_work_sigma,
result.model_fits ? "FITS" : "DOES NOT FIT");
} else {
logger.Info("Model validation: the model asserts no other enantiomorph for these data and "
@@ -861,8 +914,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
result.indexing_op = gemmi::Op::identity();
result.r_work_sigma = 0.0;
SetModelPositions(st.models[0], as_read);
compute_model_factors(mdl);
real = place_and_fit(mdl, obs, std::nullopt);
compute_model_factors(mdl, d_min);
real = place_and_fit(mdl, obs, std::nullopt, d_min);
}
Fit &best = real.fit;