From ee350f1f74e3d3c04d253e34e879fd7c40395ded Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 11 Oct 2026 21:30:34 +0200 Subject: [PATCH] rugnux --model: a wrong model confirms no shell, fits no hand, and fits the GPU reserve 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) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/CPU_DATA_ANALYSIS_DECISIONS.md | 10 +- docs/RUGNUX_REPORT.md | 28 +++-- .../structure_refinement/ModelValidation.cpp | 106 +++++++++++++++--- .../structure_refinement/ModelValidation.h | 31 +++-- rugnux/ResultReport.cpp | 44 +++++--- tests/ModelValidationTest.cpp | 77 ++++++++++++- tests/ResultReportTest.cpp | 19 +++- 7 files changed, 255 insertions(+), 60 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 095a7a3ae..cc7bc70d5 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -180,7 +180,7 @@ a cut calibrated on one pair misjudges the next. In one measured arm an unrelate The only null that fits both the model and the data is therefore **made out of them**: the same model is re-oriented at random about its own centroid and run through the identical path — the same scaling, the same rigid-body placement, the same R — nine times, and the real fit is asked how far above the -resulting distribution it sits (`MODEL_FIT_SIGMA`, accepted at 3σ). The replicates are placed as well +resulting distribution it sits (`MODEL_FIT_SIGMA`, accepted at 4.51: see below). The replicates are placed as well as fitted, or the comparison would be between a placed model and unplaced nulls and the margin would be inflated by the placement rather than by the model. Measured on one rotation data set: the crystal's own model +15.0σ, an unrelated protein +1.8σ, and the correct model rigidly rotated 90° −1.0σ. @@ -217,7 +217,13 @@ The count is nine, and it is the *spread* that sets it rather than the mean: the sample that happens to come out narrow is what turns a model that does not fit into one that appears to. Measured on the case nearest the gate — an unrelated protein at 1.83σ against a threshold of 3 — the chance of reading over the gate on a different seed is 31 % at $n=3$, 17 % at $n=5$ and 6 % at -$n=9$. Nine is affordable only because the replicates run at once: the null costs the slowest of them +$n=9$. The gate itself allows for it: the intended cut is 3σ of the null (one-sided, 0.135 %), and a +new draw from a null known only through $n$ draws lies $t_{n-1}\sqrt{1+1/n}$ sample sds above their +mean at that probability — 4.51 for $n=9$. Read as if the sd were known, an unrelated helical protein +on a helical crystal was accepted at +3.05 against nine draws; against 48 it read +2.09, and the best +of the 48 random placements fitted as well as it did. The deposited models of the open arm read +7.0 +and up. The indexing lead (§14.6) stays at 3σ of its own null: it is only read for a model that has +passed this gate, and the deposited models that settle an indexing lead by +4.2σ and up. Nine is affordable only because the replicates run at once: the null costs the slowest of them rather than their sum, and the slowest of nine is barely above the slowest of five. The null costs one full fit and one rigid-body placement per replicate. The replicates share nothing diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 96c20d134..326ac6e9e 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -873,7 +873,7 @@ as usual. | `MODEL_FIT` | `ACCEPTED` \| `REJECTED` \| `NOT_TESTED` — whether the model may decide anything, or had nothing to decide | | `MODEL_FIT_STATISTIC` | What the verdict was taken on; `R_WORK` | | `MODEL_FIT_VALUE`, `MODEL_FIT_NULL_MEAN`, `MODEL_FIT_NULL_SD`, `MODEL_FIT_NULL_REPLICATES` | The real fit, and the null of the same model in random orientations. Absent when `NOT_TESTED` | -| `MODEL_FIT_SIGMA` | The real fit above that null, in its standard deviations. Signed. Absent when `NOT_TESTED` | +| `MODEL_FIT_SIGMA` | The real fit above that null, in its sample standard deviations. Signed. `ACCEPTED` from +4.51: 3 σ of the null, allowing for its mean and spread being estimated from nine draws (Student's t). Absent when `NOT_TESTED` | | `MODEL_DECISIONS_TAKEN` | `NONE` \| `ENANTIOMORPH` \| `INDEXING` \| `ENANTIOMORPH+INDEXING` | | `MODEL_ENANTIOMORPH_ADOPTED`, `MODEL_INDEXING_OPERATOR` | The two decisions individually; `x,y,z` is no reindexing | | `MODEL_CHANGE_OF_BASIS`, `MODEL_SETTING_AS_READ` | Present only where the model was written in another description of the lattice (other axes, other centring, or another point group) and was put into the data's to be scored. Where such a model fits, the data are then written in the model's setting instead (`SETTING_SOURCE=MODEL`) and these keys are absent, since the model no longer moves. An alternative indexing of the same point group is settled here only where a reference MTZ (`-z`) has already fixed the data's indexing; otherwise the data are reindexed into the model's, and that is `MODEL_INDEXING_OPERATOR` | @@ -933,19 +933,29 @@ observed value is the merged intensity itself rather than the French–Wilson |* that amplitude is a posterior mean under a Wilson prior, which pulls a weak reflection towards its shell mean and would show up as correlation in exactly the outer shells this number is read in. -Nothing here was refined against these reflections — the model is placed and scaled with eleven -parameters — so there is no work/free distinction to draw: the correlation is unbiased on **all** the -reflections of a shell, not only the few hundred free ones, and `SIGMA` is correspondingly sharp. +The model is placed and scaled with eleven parameters, so there is no work/free distinction to draw: +the correlation is formed on **all** the reflections of a shell, not only the few hundred free ones. +Both sides are divided by the reflection's symmetry enhancement ε, which they share whatever the model. | key | meaning | |---|---| | `CC_MODEL_OVERALL`, `CC_MODEL_REFLECTIONS` | The correlation over every reflection in the table, and how many — the `N` column sums to it. Like any overall correlation it is shell-weighted and can take any value between the best shell and the worst; the table is what to read | -| `CC_MODEL_CONFIRMED_TO_D_MIN` | The finest shell whose correlation reaches 3 σ, or `NONE`. A **lower bound** on the useful resolution | -| the `D_MIN / CC_MODEL / N / SIGMA` table | Per shell: the correlation, the reflections it was formed on, and how far above zero it sits (Fisher's transform, `atanh(CC)·√(N−3)`) | +| `CC_MODEL_CONFIRMED_TO_D_MIN` | The finest shell that, together with the shell coarser than it, reaches 3 σ, or `NONE`. A **lower bound** on the useful resolution | +| the `D_MIN / CC_MODEL / N / CC_NULL / SIGMA` table | Per shell: the correlation, the reflections it was formed on, its null, and how far above the null it sits | -**Read it in one direction only.** A shell whose correlation is significantly above zero carries -signal — a model cannot agree by accident with measurements it was never fitted to — so -`CC_MODEL_CONFIRMED_TO_D_MIN` is evidence for keeping *more* data. A shell whose correlation is near +**The null.** A correlation above zero is not evidence on its own. The data and a scaled model share +everything that is smooth in reciprocal space — the falloff across a shell, the anisotropic B (fitted +to these very data) and the bulk solvent — whatever the model is, so the deposited model of an +unrelated protein correlates with the data at 10–15 σ by Fisher's transform in most shells. `CC_NULL` +is the correlation of the data at *h* with the model at *h* + *d*, averaged over the 26 nearest lattice +offsets *d*: it keeps what is smooth and loses what the model knows of *h* itself. `SIGMA` is +(atanh `CC_MODEL` − mean atanh) / sd over those 26, the sd never narrower than Fisher's 1/√(N−3). A +shell confirms only beside a coarser shell that does as well: one shell alone clears 3 σ by chance +among a dozen, while signal runs on over neighbouring shells. + +**Read it in one direction only.** A shell whose correlation is significantly above its null carries +signal the model accounts for reflection by reflection, so `CC_MODEL_CONFIRMED_TO_D_MIN` is evidence +for keeping *more* data. A shell whose correlation is near zero says nothing about the data: the model may be incomplete, in the wrong hand, or simply wrong for this crystal, and cutting on it would be cutting because the model is poor. Nothing in the pipeline acts on these numbers; they are reported and no more. This is the same asymmetry cryo-EM works under, diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index e8ba4d186..0b971faa2 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -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 &merged, #ifdef JFJOCH_USE_CUDA std::unique_ptr 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 &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(std::max(nthreads, 1), NULL_REPLICATES) + 1, @@ -1257,7 +1283,7 @@ ModelValidationResult Validate(const std::vector &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 &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 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 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 shell_cc(report_shell_d_min.size()); + std::vector> shell_offset(report_shell_d_min.size()); std::vector 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 &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 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(); diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index ebd6b26cf..c4bbc1a0c 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -80,23 +80,32 @@ struct ModelValidationResult { // posterior mean under a Wilson prior, which pulls a weak reflection towards its shell mean, and // that is a correlation with the prior in exactly the outer shells this statistic is read in. // - // Nothing was refined against these reflections - the model is placed and scaled with a handful of - // parameters - so there is no free/work distinction to make here: the correlation is unbiased on - // all of them, which is worth having, because the free set alone in an outer shell is a few - // hundred reflections and the correlation on it correspondingly noisy. + // The model is placed and scaled with a handful of parameters, so there is no free/work + // distinction to make here and the correlation is formed on all the reflections: the free set + // alone in an outer shell is a few hundred reflections and the correlation on it correspondingly + // noisy. // - // The test it supports is ONE-SIDED. A correlation significantly above zero in a shell proves the - // shell carries signal, because a model cannot invent agreement with measurements it never saw, so - // it is evidence that a resolution limit could be pushed OUTWARDS. A correlation near zero proves - // nothing at all - the model may be the thing at fault - and must never pull a limit in. + // The test it supports is ONE-SIDED. A correlation significantly above its null in a shell is + // evidence that the shell carries signal, so it argues that a resolution limit could be pushed + // OUTWARDS. A correlation near zero proves nothing at all - the model may be the thing at fault - + // and must never pull a limit in. "Significantly" is the hard part: the data and a scaled model + // agree in everything smooth in reciprocal space - the falloff across a shell, the anisotropic B + // that was fitted to these data, the bulk solvent - whatever the model is, so a model of an + // unrelated protein correlates with them far beyond what independent reflections allow. The + // null is made from the same data and model for that reason (see sigma below). struct ModelDataCCShell { float d_min = 0.0f; // the shell's high-resolution bound, as the merge table labels its rows double cc = NAN; int n = 0; // reflections the correlation was formed from - // How far the correlation sits above zero, as a Fisher-z score: atanh(cc)*sqrt(n-3). Intensities - // are not bivariate normal, so this is the same approximation XDS makes when it calls a shell's - // CC1/2 significant, and it is a guide to the strength of the claim rather than an exact p. + // How far the correlation sits above its null, in Fisher-z: the null is the correlation of the + // data at h with the model at h + d for the 26 nearest lattice offsets d, which keeps + // everything the two share smoothly in reciprocal space and loses what the model knows of h + // itself. sigma = (atanh(cc) - mean) / sd over those 26, the sd never narrower than Fisher's + // 1/sqrt(n-3). A guide to the strength of the claim rather than an exact p. double sigma = NAN; + // The correlation those offset pairs reach on average: what an unrelated model scaled to these + // data would show in this shell. + double cc_null = NAN; // The shell's mean merged intensity over the mean |F_model|^2 of the SCALED model - the model // having already been given an overall scale, an anisotropic B and a bulk solvent fitted over // the whole range. That family is smooth and nearly exponential in s^2, so it can follow a diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 353a46c35..f57c05fab 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2375,14 +2375,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (!mv.cc_model_shells.empty()) { Add(s, KeyReal("CC_MODEL_OVERALL", mv.cc_model_overall, "{:.4f}")); Add(s, KeyInt("CC_MODEL_REFLECTIONS", mv.cc_model_n)); - // The finest shell in which the model confirms signal, at the significance XDS calls a - // shell's CC1/2 established at. A LOWER BOUND on the useful resolution and nothing else: - // shells past it are not thereby empty, only unconfirmed by this model. + // The finest shell in which the model confirms signal: it and the shell coarser than it + // both 3 sigma above their null. One shell alone clears that by chance among a dozen + // (measured on wrong models: once in 216 shells), while signal runs on over neighbouring + // shells. A LOWER BOUND on the useful resolution and nothing else: shells past it are not + // thereby empty, only unconfirmed by this model. constexpr double CC_MODEL_SIGNIFICANT = 3.0; float confirmed = 0.0f; - for (const auto &sh : mv.cc_model_shells) - if (sh.sigma >= CC_MODEL_SIGNIFICANT) - confirmed = sh.d_min; + for (size_t i = 1; i < mv.cc_model_shells.size(); i++) + if (mv.cc_model_shells[i].sigma >= CC_MODEL_SIGNIFICANT + && mv.cc_model_shells[i - 1].sigma >= CC_MODEL_SIGNIFICANT) + confirmed = mv.cc_model_shells[i].d_min; if (confirmed > 0.0f) Add(s, KeyReal("CC_MODEL_CONFIRMED_TO_D_MIN", confirmed, "{:.2f}")); else @@ -2453,8 +2456,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Blank()); Add(s, Prose(" CC(model, data): the correlation of the merged intensities with the placed, scaled\n" " model's |F_model|^2, on the same shells the merge table above reports, so a row here\n" - " can be read straight across from that shell's CC1/2. SIGMA is how far the correlation\n" - " sits above zero (Fisher's transform on N reflections). I/I_MODEL is the shell's\n" + " can be read straight across from that shell's CC1/2, both on I/epsilon. CC_NULL is what\n" + " the same data reach against the model at the 26 neighbouring lattice points - what\n" + " any model scaled to these data shares with them, through the falloff, the anisotropic\n" + " B and the solvent - and SIGMA is how far CC_MODEL sits above it, in units of its\n" + " spread over the 26 (Fisher's transform, never narrower than on N independent\n" + " reflections). I/I_MODEL is the shell's\n" " mean merged intensity over the model's, AFTER the model has been given its own\n" " overall scale, anisotropic B and bulk solvent - so it reads 1.000 where the two\n" " agree on how the data fall off with resolution, and drifts away from 1.000 where\n" @@ -2462,14 +2469,15 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " set of symmetry equivalents; R_meas, CC1/2 and chi^2 all cancel one.\n")); ReportEntry e; e.kind = ReportEntry::Kind::Table; - e.table.columns = {"D_MIN", "CC_MODEL", "N", "SIGMA", "I_OVER_MODEL"}; - e.table.text_header = " D_MIN CC_MODEL N SIGMA I/I_MODEL\n" - " -------- -------- --------- -------- ---------"; + e.table.columns = {"D_MIN", "CC_MODEL", "N", "CC_NULL", "SIGMA", "I_OVER_MODEL"}; + e.table.text_header = " D_MIN CC_MODEL N CC_NULL SIGMA I/I_MODEL\n" + " -------- -------- --------- -------- -------- ---------"; for (const auto &sh : mv.cc_model_shells) { if (sh.n == 0) continue; e.table.text_rows.push_back(fmt::format( - " {:8.2f} {:8.4f} {:9d} {:>8s} {:>9s}", sh.d_min, sh.cc, sh.n, + " {:8.2f} {:8.4f} {:9d} {:>8s} {:>8s} {:>9s}", sh.d_min, sh.cc, sh.n, + std::isfinite(sh.cc_null) ? fmt::format("{:.4f}", sh.cc_null) : std::string("-"), std::isfinite(sh.sigma) ? fmt::format("{:+.1f}", sh.sigma) : std::string("-"), std::isfinite(sh.i_over_model) ? fmt::format("{:.3f}", sh.i_over_model) : std::string("-"))); @@ -2477,16 +2485,20 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, row.push_back(KeyReal("", sh.d_min, "{:.2f}").value); row.push_back(KeyReal("", sh.cc, "{:.4f}").value); row.push_back(KeyInt("", sh.n).value); + row.push_back(KeyReal("", sh.cc_null, "{:.4f}").value); row.push_back(KeyReal("", sh.sigma, "{:+.1f}").value); row.push_back(KeyReal("", sh.i_over_model, "{:.3f}").value); e.table.cells.push_back(std::move(row)); } Add(s, std::move(e)); Add(s, Prose("\n" - " The table reads in ONE direction only. A shell whose correlation is significantly\n" - " above zero carries signal, because a model cannot agree by accident with measurements\n" - " it was never fitted to, so CC_MODEL_CONFIRMED_TO_D_MIN is a lower bound on the useful\n" - " resolution and an argument for keeping MORE data. A shell whose correlation is near\n" + " The table reads in ONE direction only. A correlation well above zero is no evidence\n" + " on its own: a model of an unrelated protein, scaled to these data, correlates with\n" + " them too, through everything the two share smoothly across a shell - which is what\n" + " CC_NULL measures. A shell whose correlation sits significantly above its null, beside\n" + " a coarser shell that does as well, carries signal the model accounts for reflection by\n" + " reflection, so CC_MODEL_CONFIRMED_TO_D_MIN is a lower bound on the useful resolution\n" + " and an argument for keeping MORE data. A shell whose correlation is near\n" " zero says nothing about the data: the model may be incomplete, in the wrong hand or\n" " simply wrong for this crystal, and cutting data on it would be cutting because the\n" " model is poor. Nothing in the pipeline acts on these numbers.")); diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 47720f21a..47b7ddf73 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -76,9 +76,9 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00 // The same cell in a polar group: b is the unique axis, so the origin is free along y alone. const char *kPolarCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 1 2 1 2\n"; - std::string ClusterPdb(const char *cryst = kCryst) { + std::string ClusterPdb(const char *cryst = kCryst, unsigned seed = 20260902) { std::string pdb = cryst; - std::mt19937 rng(20260902); + std::mt19937 rng(seed); std::uniform_real_distribution x(2, 14), y(2, 16), z(2, 18); char line[96]; for (int i = 1; i <= 150; i++) { @@ -720,6 +720,79 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") // The model path's structure factors come from FFTW; they must be gemmi's own transform to float // precision, in the same layout, so prepare_asu_data() reads the same reflections from either. + +// A wrong model may come in a group that does not keep the data's metric - a cubic group on an +// orthorhombic cell - and then its ASU holds an observation at another resolution than the data's own: +// more of them land inside the rigid body's finest zone than the data put there. The GPU engines are +// reserved for the zone as it is counted, so the validation runs instead of failing on the device. +TEST_CASE("ModelValidation_WrongModelInAGroupOffTheDataMetric", "[ModelValidation]") { + Logger logger("ModelValidation_WrongModelInAGroupOffTheDataMetric"); + const auto data_path = WriteTemp("off_metric_data.pdb", ClusterPdb().c_str()); + const auto wrong_path = WriteTemp("off_metric_wrong.pdb", + ClusterPdb("CRYST1 40.000 40.000 40.000 90.00 90.00 90.00 I 21 3 24\n", 7).c_str()); + auto obs = ModelReferenceIntensities(data_path, {}, {}, 2.0, logger); + REQUIRE(obs.size() > 1000); + for (size_t i = 0; i < obs.size(); i++) { + obs[i].sigma = std::max(1.0f, 0.05f * obs[i].I); + obs[i].F = std::sqrt(std::max(0.0f, obs[i].I)); + obs[i].rfree_flag = (i % 20) == 0; + } + const std::string prefix = (std::filesystem::temp_directory_path() / "off_metric_test").string(); + const auto result = + ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, + wrong_path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21")); + CHECK(result.ok); + CHECK(!result.adopted_model_enantiomorph); + + for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"}) + std::filesystem::remove(prefix + suffix); + std::filesystem::remove(data_path); + std::filesystem::remove(wrong_path); +} + +// A model of something else entirely, scaled to the data, shares their falloff across every shell, and +// by Fisher's transform on independent reflections that alone reads as a significant correlation. Against +// the null of the data at h with the model at h + d - which shares the falloff too - it must not. +TEST_CASE("ModelValidation_CCModelUnrelatedModelConfirmsNothing", "[ModelValidation]") { + Logger logger("ModelValidation_CCModelUnrelatedModelConfirmsNothing"); + const auto data_path = WriteTemp("cc_null_data.pdb", ClusterPdb().c_str()); + const auto wrong_path = WriteTemp("cc_null_wrong.pdb", ClusterPdb(kCryst, 7).c_str()); + auto obs = ModelReferenceIntensities(data_path, {}, {}, 2.0, logger); + std::erase_if(obs, [](const MergedReflection &r) { return r.d > 8.0f; }); + REQUIRE(obs.size() > 1000); + for (size_t i = 0; i < obs.size(); i++) { + // An anisotropic falloff, B 30 A^2 more along a* than along b* and c*: the validation fits an + // anisotropic B to the data and puts it on the model, whatever the model is. + const double s_a = obs[i].h / 30.0; + obs[i].I = static_cast(obs[i].I * std::exp(-30.0 * s_a * s_a / 2.0)); + obs[i].sigma = std::max(1.0f, 0.05f * obs[i].I); + obs[i].F = std::sqrt(std::max(0.0f, obs[i].I)); + obs[i].rfree_flag = (i % 20) == 0; + } + const std::vector shells{4.0f, 3.0f, 2.5f, 2.2f, 2.0f}; + const std::string prefix = (std::filesystem::temp_directory_path() / "cc_null_test").string(); + const auto result = + ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}, + wrong_path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), + /*probe_indexing_ambiguity=*/false, 1, 1.0, shells); + REQUIRE(result.ok); + REQUIRE(result.cc_model_shells.size() == shells.size()); + double fisher_max = 0; + for (const auto &sh : result.cc_model_shells) { + const double fisher = std::atanh(sh.cc) * std::sqrt(sh.n - 3.0); + fisher_max = std::max(fisher_max, fisher); + logger.Info("unrelated model {:.2f} A: CC {:.3f} (null {:.3f}) on {}: {:+.1f} sigma, {:+.1f} by Fisher alone", + sh.d_min, sh.cc, sh.cc_null, sh.n, sh.sigma, fisher); + CHECK(sh.sigma < 3.0); + } + // The trap is there: on independent reflections some shell would have read as significant. + CHECK(fisher_max > 3.0); + + for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"}) + std::filesystem::remove(prefix + suffix); + std::filesystem::remove(data_path); + std::filesystem::remove(wrong_path); +} TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") { // An even grid, and an odd one on every axis (FFTW and pocketfft split odd lengths differently). const auto size = GENERATE(std::array{20, 24, 30}, std::array{15, 21, 27}); diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index fd742ba43..757095e39 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -461,10 +461,11 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { CHECK(untested_text.find("was tried and REJECTED") == std::string::npos); // CC(model, data) by shell, with the deepest shell short of significance: the confirmed limit is - // the last shell that reached it, and it is a lower bound - the shell past it is not thereby empty. + // the last shell that reached it beside a coarser one that did, and it is a lower bound - the + // shell past it is not thereby empty. ModelValidationResult with_cc = good; - with_cc.cc_model_shells = {{3.20f, 0.9412, 4210, 91.4}, {2.10f, 0.5533, 3980, 40.2}, - {1.80f, 0.0412, 2110, 1.9}}; + with_cc.cc_model_shells = {{3.20f, 0.9412, 4210, 41.4, 0.1211}, {2.10f, 0.5533, 3980, 20.2, 0.0832}, + {1.80f, 0.0412, 2110, 1.9, 0.0105}}; with_cc.cc_model_overall = 0.8123; with_cc.cc_model_n = 10300; result.model_validation = with_cc; @@ -472,8 +473,8 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { CHECK(cc_text.find("\nCC_MODEL_OVERALL= 0.8123\n") != std::string::npos); CHECK(cc_text.find("\nCC_MODEL_REFLECTIONS= 10300\n") != std::string::npos); CHECK(cc_text.find("\nCC_MODEL_CONFIRMED_TO_D_MIN= 2.10\n") != std::string::npos); - CHECK(cc_text.find("D_MIN CC_MODEL N SIGMA") != std::string::npos); - CHECK(cc_text.find(" 1.80 0.0412 2110 +1.9") != std::string::npos); + CHECK(cc_text.find("D_MIN CC_MODEL N CC_NULL SIGMA") != std::string::npos); + CHECK(cc_text.find(" 1.80 0.0412 2110 0.0105 +1.9") != std::string::npos); CHECK(cc_text.find("The table reads in ONE direction only") != std::string::npos); // The comparable R rides on the same shells. R_WORK and R_FREE are what they always were - the @@ -510,6 +511,14 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { CHECK(RenderResultReport("p", "in.h5", x, result) .find("\nCC_MODEL_CONFIRMED_TO_D_MIN= NONE\n") != std::string::npos); + // One shell over the bar on its own, its neighbours not: chance among a dozen shells, not signal. + ModelValidationResult lone = with_cc; + lone.cc_model_shells[0].sigma = 0.4; + lone.cc_model_shells[1].sigma = 5.7; + result.model_validation = lone; + CHECK(RenderResultReport("p", "in.h5", x, result) + .find("\nCC_MODEL_CONFIRMED_TO_D_MIN= NONE\n") != std::string::npos); + result.model_validation = good; ModelValidationResult failed;