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

rugnux --model: a wrong model confirms no shell, fits no hand, and fits the GPU reserve

Browse Source 
Three defects the battery's wrong-model null found (5 sets x 5 deposited models of unrelated proteins):

CC(model, data) significance. The per-shell sigma was Fisher's transform on independent reflections,
but the data and a scaled model share the falloff across a shell, the anisotropic B fitted to these
very data, the bulk solvent and the symmetry enhancement epsilon - whatever the model is. 21 of 24
wrong models read CC_MODEL_CONFIRMED_TO_D_MIN to the data's d_min (10-15 sigma in most shells).
The correlation is now formed on I/epsilon and |F_model|^2/epsilon, and its null is empirical: the
data at h against the model at h + d for the 26 nearest lattice offsets, which keeps everything smooth
in reciprocal space and loses what the model knows of h. SIGMA = (z - mean) / max(sd, 1/sqrt(N-3)) in
Fisher z over those 26, a new CC_NULL column shows the null, and a shell confirms only beside a coarser
shell that also clears 3 sigma (one shell alone did so once in 216 wrong-model shells). Prototyped on
the p_maps.mtz of a full open-arm run and the null runs: wrong models confirm 0/24 (was 21/24), the
deposited models confirm 197/200 (was 198/200 in the prototype's reading; the one lost has CC_MODEL
0.05-0.13 in every shell). The report no longer says a model cannot agree by accident.

MODEL_FIT gate. Accepted at 3.0 sample sds of a nine-draw null, as if the sd were known; an unrelated
P1 helical model on a P3(1) helical crystal was ACCEPTED at +3.05. Against 48 draws it reads +2.09 and
the best random placement fits as well as it does, and drawing the null about the cell axes instead
changes nothing (+2.15): the cause is the nine-draw sd, not the orientation. The gate is now 3 sigma
allowing for that: t(8) * sqrt(1 + 1/9) = 4.51. The deposited models of the open arm read +7.0 and up.
The indexing lead stays at 3 sigma (INDEXING_MARGIN_SIGMA): it is read only behind the fit gate, and
a deposited model settles one at +4.15.

Rigid-body GPU reserve. The finest zone's observations were counted by the data's d, but the zones
count them by the d of the index in the model's ASU, and a model's group need not keep the data's
metric (an I2(1)3 model on I222 data: 2320 in the zone against 1712 reserved), so the validation died
with "zone over the engine's reserve". Counted the zones' way, under every indexing probed.

Tests: ModelValidation_CCModelUnrelatedModelConfirmsNothing (Fisher alone reads +4.7 on an unrelated
model with an anisotropic falloff; the new sigma stays under 1), ModelValidation_WrongModelInAGroupOffTheDataMetric
(fails with the old count), ResultReport_ModelValidationSection (CC_NULL column, lone-shell rule).

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
This commit is contained in:
leonarski_f
2026-10-11 21:30:34 +02:00
co-authored by Claude Opus 5.5
parent bd3c4efab6
commit ee350f1f74
7 changed files with 255 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/structure_refinement/ModelValidation.cpp
+91 -15
View File
@@ -142,7 +142,21 @@ constexpr int NULL_MAX_REDRAWS = 1000;
// two on its own. Measured at nine replicates, the crystal's own model reads +17.7 sigma and an
// unrelated protein +1.7, so the cut sits in a gap an order of magnitude wider than the sd it is
// measured in.
constexpr double MODEL_FIT_SIGMA = 3.0;
//
// The cut is 3 sigma of the null - one-sided, 0.135% - but the null's mean and sd are not known, only
// estimated from NULL_REPLICATES draws, and the fit is one more draw. At the same probability that
// draw lies t(n-1) * sqrt(1 + 1/n) sample sds above the sample mean: 4.28 * 1.054 = 4.51 for n = 9.
// Read as if the sd were known, a narrow sample is what lets a model that does not fit through:
// measured, an unrelated helical protein on a helical crystal read +3.05 against nine draws and +2.09
// against 48, whose best random placement fitted as well as it did. Across the open arm the deposited
// models read +7.0 and up.
constexpr double MODEL_FIT_SIGMA = 4.51;
// The lead an indexing must take over its runner-up, in sd of the lead random placements take. It is
// read only for a model that has already passed MODEL_FIT_SIGMA, so it is never all that stands between
// a wrong model and a decision; across the open arm the deposited models that settle an indexing lead
// by +4.2 sigma and up.
constexpr double INDEXING_MARGIN_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.
@@ -778,11 +792,23 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
#ifdef JFJOCH_USE_CUDA
std::unique_ptr<RigidBodyGPUPool> rigid_body_engines;
if (rigid_body_gpu) {
// Counted as the zones count them: by the resolution of the index once it is in the model's
// ASU, under each indexing the data may be placed in. The model's group need not keep the
// data's metric - a wrong model's cubic group on an orthorhombic cell sends h to an index of
// another resolution - so the data's own d can put a reflection on the wrong side of the zone.
const double finest_zone = RigidBodyLadder(d_min).back();
size_t zone_observations = 0;
for (const MergedReflection &r : obs)
if (!std::isnan(r.F) && r.d >= finest_zone)
++zone_observations;
const gemmi::GroupOps zone_gops = sg->operations();
const gemmi::ReciprocalAsu zone_asu(sg);
auto in_finest_zone = [&](const std::vector<MergedReflection> &v) {
size_t n = 0;
for (const MergedReflection &r : v)
if (!std::isnan(r.F) && ucell.calculate_d(zone_asu.to_asu({r.h, r.k, r.l}, zone_gops).first) >= finest_zone)
++n;
return n;
};
size_t zone_observations = in_finest_zone(obs);
for (const gemmi::Op &law : reindex_ops)
zone_observations = std::max(zone_observations, in_finest_zone(ReindexReflections(obs, law)));
// One engine for the real model and one for each replicate of the null placed beside it.
rigid_body_engines = RigidBodyGPUPool::Create(st.models[0], ucell, *sg, d_min, zone_observations,
std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES) + 1,
@@ -1257,7 +1283,7 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
// 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.
const bool decided = result.model_fits && result.indexing_margin_sigma >= MODEL_FIT_SIGMA;
const bool decided = result.model_fits && result.indexing_margin_sigma >= INDEXING_MARGIN_SIGMA;
result.indexing_decided = decided;
if (decided)
result.indexing_op = indexing.op;
@@ -1399,7 +1425,39 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
return i;
return report_shell_d_min.size();
};
// The correlation is formed on I/epsilon and |F_model|^2/epsilon: a reflection on a symmetry
// element is enhanced by epsilon in the data and in the model alike, and that shared factor
// is a correlation that any model in this group has, right or wrong.
//
// How far it sits above zero is NOT read off Fisher's transform alone. That transform treats
// the reflections of a shell as independent pairs, and they are not: the data and the model
// share the falloff across the shell, the anisotropic B (fitted to these very data) and the
// bulk solvent, so an unrelated model correlates with them well above what N independent
// pairs allow (measured: a wrong model of another protein reaches 10-15 sigma by Fisher's
// transform in most shells). The null is therefore made empirically, on the same data and the
// same model: the data at h against the model at h + d, for each of the 26 nearest lattice
// offsets d. Everything smooth in reciprocal space - the radial falloff, the anisotropy, the
// solvent - is the same at h and at h + d and stays in those correlations; what the model
// knows of the structure factor at h itself does not. The shell's score is how far its
// correlation sits above the mean of those 26, in units of their spread, never narrower than
// Fisher's 1/sqrt(N - 3).
struct OffsetSums {
CorrelationCoefficient cc;
int n = 0;
};
constexpr int OFFSETS = 26;
std::vector<gemmi::Miller> offsets;
for (int a = -1; a <= 1; a++)
for (int b = -1; b <= 1; b++)
for (int c = -1; c <= 1; c++)
if (a != 0 || b != 0 || c != 0)
offsets.push_back({a, b, c});
std::unordered_map<long, double> model_i; // |F_model|^2 / epsilon, keyed in the model's ASU
model_i.reserve(fmodel.v.size());
for (const auto &hv : fmodel.v)
model_i[hkl_key(hv.hkl)] = std::norm(hv.value) / gops.epsilon_factor(hv.hkl);
std::vector<CorrelationCoefficient> shell_cc(report_shell_d_min.size());
std::vector<std::array<OffsetSums, OFFSETS>> shell_offset(report_shell_d_min.size());
std::vector<int> shell_n(report_shell_d_min.size(), 0);
// Sums for the radial mismatch: the shell's mean observed intensity against the mean
// intensity of the already-scaled model. See ModelDataCCShell::i_over_model.
@@ -1415,23 +1473,41 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
if (bin >= shell_cc.size()) // finer than the finest shell the merge reported
continue;
const double Ic = std::norm(hv.value); // |F_model|^2
shell_cc[bin].Add(it->second.I, Ic);
const double eps = gops.epsilon_factor(hv.hkl);
const double I = it->second.I / eps;
shell_cc[bin].Add(I, Ic / eps);
for (int o = 0; o < OFFSETS; o++) {
const gemmi::Miller q{hv.hkl[0] + offsets[o][0], hv.hkl[1] + offsets[o][1], hv.hkl[2] + offsets[o][2]};
const auto m = model_i.find(hkl_key(asu.to_asu(q, gops).first));
if (m == model_i.end())
continue;
shell_offset[bin][o].cc.Add(I, m->second);
++shell_offset[bin][o].n;
}
shell_sum_I[bin] += it->second.I;
shell_sum_Im[bin] += Ic;
++shell_n[bin];
overall_cc.Add(it->second.I, Ic);
overall_cc.Add(I, Ic / eps);
++overall_n;
}
// Fisher's transform against a null of zero correlation. Below four reflections there is no
// score to give, and a correlation of exactly +-1 has no finite one.
auto fisher_sigma = [](double cc, int n) {
return (n > 3 && std::fabs(cc) < 1.0) ? std::atanh(cc) * std::sqrt(n - 3.0) : NAN;
};
result.cc_model_shells.reserve(report_shell_d_min.size());
for (size_t i = 0; i < report_shell_d_min.size(); i++) {
const double cc = shell_cc[i].GetCC();
result.cc_model_shells.push_back({report_shell_d_min[i], cc, shell_n[i],
fisher_sigma(cc, shell_n[i]),
std::vector<double> z_null;
for (const OffsetSums &o : shell_offset[i]) {
const double c = o.cc.GetCC();
if (o.n > 3 && std::fabs(c) < 1.0)
z_null.push_back(std::atanh(c));
}
double sigma = NAN, cc_null = NAN;
// Below four reflections there is no score to give, and a correlation of exactly +-1 has
// no finite one.
if (shell_n[i] > 3 && std::fabs(cc) < 1.0 && z_null.size() >= 2) {
const auto [z_mean, z_sd] = mean_sd(z_null);
cc_null = std::tanh(z_mean);
sigma = (std::atanh(cc) - z_mean) / std::max(z_sd, 1.0 / std::sqrt(shell_n[i] - 3.0));
}
result.cc_model_shells.push_back({report_shell_d_min[i], cc, shell_n[i], sigma, cc_null,
shell_sum_Im[i] > 0.0 ? shell_sum_I[i] / shell_sum_Im[i] : NAN});
}
result.cc_model_overall = overall_cc.GetCC();