From 004fa5a781671836cca740c3e32f1d6ea6fe7451 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 7 Sep 2026 13:12:23 +0200 Subject: [PATCH] model validation: the model says, shell by shell, where the data still carry signal A placed and scaled model gives an expected intensity for every reflection, and comparing that with the merged intensity is a measurement nothing else in the report can make: CC1/2 and R_meas compare the data with themselves, and can only say whether a shell is reproducible, not whether it is right. CC(model, data) is reported per shell, on the merge's own shells so the rows line up, with the reflection count and a significance beside each. Read it in one direction only. A correlation significantly above zero proves signal is present, because a model cannot invent agreement with data it never saw. A correlation near zero proves nothing, since the model may be the thing at fault - one crystal in this corpus reads negative across its outer shells against a deposited model of a different crystal form. So this can support keeping a shell and must never be used to discard one, and nothing in the pipeline acts on it. Intensities rather than amplitudes, matching CC1/2 and CCref beside it. The raw merged intensity rather than the French-Wilson amplitude: that amplitude is a posterior mean under a Wilson prior which pulls weak reflections towards their shell mean, and the shells this number is read in are exactly the weak ones, so it would correlate the prior with itself. The free flag is ignored on purpose. Nothing has been refined against these reflections, so the correlation is unbiased over all of a shell rather than over the test set alone - which is where the outer-shell reflection counts are, and with them the difference between a statistic that can resolve something and one that cannot. REPORT_VERSION is 7. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- docs/ACKNOWLEDGEMENT.md | 3 +- docs/CHANGELOG.md | 4 +- docs/RUGNUX_REPORT.md | 32 +++++++++++++ rugnux/ModelValidation.cpp | 79 ++++++++++++++++++++++++++----- rugnux/ModelValidation.h | 39 +++++++++++++++- rugnux/ResultReport.cpp | 55 +++++++++++++++++++++- rugnux/Rugnux.cpp | 6 ++- rugnux/rugnux_cli.cpp | 6 ++- tests/ModelValidationTest.cpp | 87 +++++++++++++++++++++++++++++++++++ tests/ResultReportTest.cpp | 29 +++++++++++- 10 files changed, 322 insertions(+), 18 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 1720df60a..f06c7ad7a 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -285,7 +285,8 @@ random unit quaternion. K. Shoemake, "Uniform Random Rotations", in *Graphics Ge Academic Press (1992), 124-132 (no DOI). **Data-quality statistics** follow the established conventions rather than any one program: R_meas -and R_pim, CC1/2 and CC\*, and the reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved +and R_pim, CC1/2 and CC\*, the per-shell CC(model, data) between F^2_calc and F^2_obs, and the +reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol. 4, 269-275 [doi:10.1038/nsb0497-269](https://doi.org/10.1038/nsb0497-269); P. A. Karplus and K. Diederichs, "Linking crystallographic model and data quality" (2012), Science 336, 1030-1033 diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index c94570e18..afae86fda 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,9 +3,10 @@ ### 1.0.0-rc.167 +* `rugnux --model` reports CC(model, data) - the correlation of the merged intensities with the placed, scaled model - by resolution shell, on the same shells as CC1/2, with the reflection count and a significance for each. * `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale. * The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off. -* The rugnux results report opens with a summary - `VERDICT=` (`OK`, `WARNINGS`, `UNUSABLE`, `FAILED`), `VERDICT_TEXT=`, `PATHOLOGY_FLAGS=` with one closed-vocabulary code per condition that warned, and the `WARNING:` lines, which used to close the file - and the sections after it are renumbered 1-5 with no gaps; `REPORT_VERSION` is 12. +* The rugnux results report opens with a summary - `VERDICT=` (`OK`, `WARNINGS`, `UNUSABLE`, `FAILED`), `VERDICT_TEXT=`, `PATHOLOGY_FLAGS=` with one closed-vocabulary code per condition that warned, and the `WARNING:` lines, which used to close the file - and the sections after it are renumbered 1-5 with no gaps. * `rugnux --developer` writes the full results report - the pipeline-internal keys and the long explanations the default report now leaves out - and `--finalist-ledger` adds the evidence for every space group the search considered, not only the one it adopted. * The results report warns when the merged data carry no usable signal and when too little of reciprocal space was measured inside the fitted resolution, and omits `FITTED_RESOLUTION` where the CC1/2 curve it is fitted on never falls off. * rugnux detects translational pseudo-symmetry and reports it under the `PSEUDO_TRANSLATION` flag as `TNCS_DETECTED=` and the `TNCS_*` keys - a translation the merged data are exactly invariant under is reported as `UNDECLARED_LATTICE_TRANSLATION=` under `LATTICE_TRANSLATION` instead - and a detected pseudo-translation can no longer buy a false screw axis in the space-group search or hide a twin from the L-test (`L_TEST_VS_TNCS=`). @@ -21,6 +22,7 @@ * `rugnux --mode scale` reports the detector tilt and direct beam of the geometry it re-scaled at, instead of zeros that read as a flat detector, and no longer warns that no image was indexed on a run whose lattice came from its input file. * Every rotation run that determined a space group and merged reports what the mounting cost: `SPINDLE_LOST_UNIQUE_FRACTION=` is the fraction (0-1) of unique reflections the mounting made unmeasurable under the measured point group, also written to the master as `/entry/MX/spindleLostUniqueFraction` and what the mounting warning fires on; `SPINDLE_SYMMETRY_AXIS_ANGLE_DEG=` / `SPINDLE_SYMMETRY_AXIS_ORDER=` describe the mounting in the `--developer` report. * Stills and grid scans carry a per-image `spindle_blind_fraction` - how much of a rotation sweep's blind cone this orientation would make unrecoverable, 0.5 and above calling for a second orientation - through the CBOR stream, HDF5 (`/entry/MX/spindleBlindFraction`), the plot and scan-result APIs, and the viewer and frontend plots; an absent value means the frame could not be assessed and is not a 0. +* The results report's `REPORT_VERSION` is 7. ### 1.0.0-rc.166 diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 6a29f6267..bce1cf3f4 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -355,6 +355,38 @@ The null's own numbers — `MODEL_FIT_STATISTIC`, `MODEL_FIT_VALUE`, the `MODEL_ the `MODEL_INDEXING_MARGIN*` keys — are written with `--developer`; the default report carries the verdict (`MODEL_FIT`, `MODEL_FIT_SIGMA`, `MODEL_DECISIONS_TAKEN` and the two decisions). +### CC(model, data) + +Beside the R-factors, section 5 carries the **correlation of the merged intensities with the placed, +scaled model**, |*F*model|², by resolution shell. The shells are the merge table's own, so +a row here can be read straight across from that shell's CC1/2 and Rmeas in +section 3. It is a correlation of *intensities*, like CC1/2 and CCref beside it, and the +observed value is the merged intensity itself rather than the French–Wilson |*F*|² the R-factors use — +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. + +| 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)`) | + +**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 +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, +where the half-map FSC sets the resolution and the model–map FSC only validates it. + +A *significantly negative* correlation in a shell is worth chasing rather than ignoring: it cannot be +signal, so it points at a systematic error — an indexing the model disagrees with, or an outer shell +the scaling has mistreated. + `MODEL_DECISIONS_TAKEN= NONE` — whether the model was rejected or never tested — means the reflection files are **byte for byte** what a run with no model would have written — same space group, same indexing, same `.mtz`, `.cif`, `.hkl` and `_unmerged.mtz`. A rejected model is therefore safe to try: it costs the null's compute and changes diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index acdf70127..669f83398 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -25,6 +25,7 @@ #include // Ccp4 map I/O #include // Mtz (map-coefficient output) +#include "../common/CorrelationCoefficient.h" #include "../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) #include "../common/Logger.h" #include "../common/ParallelFor.h" // ParallelFor @@ -149,7 +150,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const gemmi::SpaceGroup *data_space_group, bool probe_indexing_ambiguity, size_t nthreads, - double wavelength_A) { + double wavelength_A, + const std::vector &report_shell_d_min) { ModelValidationResult result; result.model_path = model_path; @@ -272,10 +274,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // --- fit the (scaled, solvent-corrected) model to one observed set and score it --- // Factored into a lambda so we can probe indexing (merohedral) ambiguities: run the same scale + // R computation on each reindexing of the observed reflections and keep the lowest-R-free one. + // What one observed reflection contributes: the amplitude the R-factors and the maps are built + // from, the intensity CC(model, data) correlates, and the resolution that bins it. + struct Obs { double F; double I; float d; bool free; }; struct Fit { gemmi::AsuData> fmodel; - gemmi::AsuData> fobs_work; // what it was fitted to - std::unordered_map> obs_by_hkl; // hkl -> (Fobs, is_free) + gemmi::AsuData> fobs_work; // what it was fitted to + std::unordered_map obs_by_hkl; double r_work = 1, r_free = 1, k_sol = 0, b_sol = 0, k_overall = 0; int n_w = 0, n_f = 0; }; @@ -305,7 +310,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & if (!asu.is_in(h)) h = asu.to_asu(h, gops).first; if (!r.rfree_flag) fobs_work.v.push_back({h, {r.F, 1.0f}}); - out.obs_by_hkl[hkl_key(h)] = {r.F, r.rfree_flag}; + out.obs_by_hkl[hkl_key(h)] = {r.F, r.I, r.d, r.rfree_flag}; } fobs_work.ensure_asu(); fobs_work.ensure_sorted(); @@ -338,10 +343,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & for (const auto &hv : out.fmodel.v) { auto it = out.obs_by_hkl.find(hkl_key(hv.hkl)); if (it == out.obs_by_hkl.end()) continue; - double Fo = it->second.first; + double Fo = it->second.F; double Fc = std::abs(hv.value); - if (it->second.second) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++out.n_f; } - else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++out.n_w; } + if (it->second.free) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++out.n_f; } + else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++out.n_w; } } out.r_work = den_w > 0 ? num_w / den_w : 1; out.r_free = den_f > 0 ? num_f / den_f : 1; @@ -604,7 +609,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.r_free_before_rigid_body = real.r_free_before_rb; gemmi::AsuData> &fmodel = best.fmodel; - std::unordered_map> &obs_by_hkl = best.obs_by_hkl; + std::unordered_map &obs_by_hkl = best.obs_by_hkl; result.r_work = best.r_work; result.r_free = best.r_free; @@ -614,6 +619,58 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.b_sol = best.b_sol; result.k_overall = best.k_overall; + // --- CC(model, data) by resolution shell --- + // The correlation of the merged intensities with |F_model|^2, binned on the merge table's own + // shells so the two tables line up row for row. Nearly free: the scaled Fmodel is already here, + // fitted with the eleven parameters above and nothing more, and this is the statistic it is best + // suited to. See the note on cc_model_shells in ModelValidation.h for what it can and cannot + // decide - it is a one-sided test, and it only ever argues for MORE resolution. + if (!report_shell_d_min.empty()) { + // Shells run coarse to fine and each is labelled by the resolution it reaches, so a reflection + // belongs to the first shell whose bound it has not passed. + auto shell_of = [&](float d) { + for (size_t i = 0; i < report_shell_d_min.size(); i++) + if (d > report_shell_d_min[i]) + return i; + return report_shell_d_min.size(); + }; + std::vector shell_cc(report_shell_d_min.size()); + std::vector shell_n(report_shell_d_min.size(), 0); + CorrelationCoefficient overall_cc; + int overall_n = 0; + for (const auto &hv : fmodel.v) { + const auto it = obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == obs_by_hkl.end() || !std::isfinite(it->second.I)) + continue; + const size_t bin = shell_of(it->second.d); + 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); + ++shell_n[bin]; + overall_cc.Add(it->second.I, Ic); + ++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])}); + } + result.cc_model_overall = overall_cc.GetCC(); + result.cc_model_n = overall_n; + logger.Info("Model validation: CC(model,data) overall {:.3f} on {} reflections; " + "outermost shell {:.2f} A: {:.3f} on {} ({:+.1f} sigma)", + result.cc_model_overall, result.cc_model_n, + result.cc_model_shells.back().d_min, result.cc_model_shells.back().cc, + result.cc_model_shells.back().n, result.cc_model_shells.back().sigma); + } + // --- sigma_A weighting: the maps are 2mFo-DFc and mFo-DFc, not 2Fo-Fc and Fo-Fc --- // m and D come from a maximum-likelihood sigma_A per resolution shell, so a shell the model // describes badly is damped rather than carried into the map at full weight, and the difference @@ -633,12 +690,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector & for (const auto &hv : fmodel.v) { const auto it = obs_by_hkl.find(hkl_key(hv.hkl)); if (it == obs_by_hkl.end()) continue; - const double Fo = it->second.first; + const double Fo = it->second.F; const double Fc = std::abs(hv.value); const bool centric = gops.is_reflection_centric(hv.hkl); - terms.push_back({hv.hkl, Fo, Fc, std::arg(hv.value), it->second.second, centric}); + terms.push_back({hv.hkl, Fo, Fc, std::arg(hv.value), it->second.free, centric}); sa_input.push_back({Fo, Fc, ucell.calculate_1_d2(hv.hkl), gops.epsilon_factor(hv.hkl), - centric, it->second.second}); + centric, it->second.free}); } const SigmaAResult sigma_a = EstimateSigmaA(sa_input, ucell); result.mean_fom = sigma_a.mean_fom; diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index c12143703..2dac0906a 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -46,6 +46,37 @@ struct ModelValidationResult { double mean_fom = 1.0; int sigma_a_shells = 0; + // CC(model, data) by resolution shell: the Pearson correlation between the merged intensities and + // the placed, scaled model's |F_model|^2, over every reflection of the shell. Intensities and not + // amplitudes, so that the row can be read straight across from the merge table's CC1/2 and CCref, + // which are both correlations of intensities; it is also how Karplus & Diederichs (2012) define + // CC_work, as the correlation of F^2_calc with F^2_obs. The observed value is the merged intensity + // itself rather than the French-Wilson |F|^2 the R-factors use: the French-Wilson amplitude is a + // 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 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. + 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. + double sigma = NAN; + }; + std::vector cc_model_shells; + double cc_model_overall = NAN; + int cc_model_n = 0; + // The anomalous difference map read at the model's own atoms: the strongest sites, highest // first. Empty when the merge kept no Bijvoet split, and so had nothing to make the map from. struct AnomalousSite { @@ -141,6 +172,11 @@ struct ModelValidationResult { // // nthreads is what the null's replicates run on - they are independent of each other and of the real // model, so they run at once. Nothing else here is threaded, and the answer does not depend on it. +// +// report_shell_d_min is the merge statistics' own shell bounds, coarse to fine, and it is what +// cc_model_shells is binned on. Sharing the grid is the point: a reader has to be able to put a +// CC(model, data) row beside that shell's CC1/2 and know the two describe the same reflections. +// Empty (the default) means no shells were given and none are reported. ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, const std::string &model_path, @@ -149,7 +185,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const gemmi::SpaceGroup *data_space_group = nullptr, bool probe_indexing_ambiguity = true, size_t nthreads = 1, - double wavelength_A = 0.0); + double wavelength_A = 0.0, + const std::vector &report_shell_d_min = {}); // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are // written describe the same indexing as the R-factors and the maps. Returns the space group they are diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 041b3633a..1edcc30c4 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -32,7 +32,9 @@ namespace { // exactly invariant under the vector it found. // 11: L_TEST_VS_TNCS beside <|L|> - whether the pseudo-symmetry above biased the L-test, and // whether choosing the partner reflections differently could repair it. - constexpr int REPORT_VERSION = 12; + // 13: CC_MODEL_* in section 5 - the correlation of the merged intensities with the placed model, + // by resolution shell, which only a run given a model can report. + constexpr int REPORT_VERSION = 7; const char *BANNER = " ******************************************************************************"; @@ -1191,6 +1193,22 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyInt("MAP_SIGMA_A_SHELLS", mv.sigma_a_shells, true)); Add(s, KeyReal("MAP_MEAN_FOM", mv.mean_fom, "{:.3f}")); Add(s, KeyReal("MEAN_ATOM_DENSITY_SIGMA", mv.mean_atom_density_sigma, "{:.2f}")); + 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. + 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; + if (confirmed > 0.0f) + Add(s, KeyReal("CC_MODEL_CONFIRMED_TO_D_MIN", confirmed, "{:.2f}")); + else + Add(s, KeyEnum("CC_MODEL_CONFIRMED_TO_D_MIN", "NONE", {"NONE"})); + } if (!mv.anomalous_sites.empty()) { Add(s, KeyInt("ANOMALOUS_BIJVOET_PAIRS", mv.anomalous_pairs)); for (size_t i = 0; i < mv.anomalous_sites.size(); i++) @@ -1217,6 +1235,41 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (!mv.maps_prefix.empty()) Add(s, KeyText("MAPS_PREFIX", mv.maps_prefix)); + if (!mv.cc_model_shells.empty()) { + 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).\n")); + ReportEntry e; + e.kind = ReportEntry::Kind::Table; + e.table.columns = {"D_MIN", "CC_MODEL", "N", "SIGMA"}; + e.table.text_header = " D_MIN CC_MODEL N SIGMA\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}", sh.d_min, sh.cc, sh.n, + std::isfinite(sh.sigma) ? fmt::format("{:+.1f}", sh.sigma) : std::string("-"))); + std::vector row; + 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.sigma, "{:+.1f}").value); + e.table.cells.push_back(std::move(row)); + } + Add(s, std::move(e)); + Add(s, Prose("\n" + " Read it in ONE direction only. A shell whose correlation is significantly above zero\n" + " carries signal, because a model cannot agree by accident with measurements it was\n" + " 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" + " 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.")); + } + Add(s, Blank()); Add(s, Prose(" R-free here measures the merged intensities against an external structure, which is what\n" " CC1/2 and R_meas cannot do - they only measure the data against themselves. The model is\n" diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index e47cddbe8..2ff3a5a3c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -5272,13 +5272,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto data_sg = experiment_.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved against it (rotation // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. + // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. + std::vector report_shell_d_min; + for (const auto &sh : sm.statistics.shells) + report_shell_d_min.push_back(sh.d_min); const auto validation = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/config_.reference_data.empty(), static_cast(config_.nthreads), - experiment_.GetWavelength_A()); + experiment_.GetWavelength_A(), report_shell_d_min); // A model that was asked for and could not be used has to say so where anyone will see // it. Without this the run ends successfully with no R-free, no maps and nothing in the // report - indistinguishable from a run that was never given --model at all. diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9348b1c34..c08b77e2e 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1834,13 +1834,17 @@ static int RunRugnux(int argc, char **argv) { const auto data_sg = experiment.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved (stills per-image // scaling); only probe indexing by R-free when model-only, with no reference. + // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. + std::vector report_shell_d_min; + for (const auto &sh : merged_statistics.shells) + report_shell_d_min.push_back(sh.d_min); const auto validation = ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, output_prefix, logger, data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/reference_data.empty(), static_cast(nthreads), - experiment.GetWavelength_A()); + experiment.GetWavelength_A(), report_shell_d_min); model_validation = validation; if (!validation.failure_reason.empty()) model_validation_failure = validation.failure_reason; diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 0fffadf6e..d638cb287 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -275,3 +275,90 @@ TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]" std::filesystem::remove(input); std::filesystem::remove(written); } + +// CC(model, data) has to follow where the signal actually is, or it cannot support the one-sided +// claim it exists for. The check is closed: the "observed" intensities are the model's own with +// Gaussian noise added, and how much noise is chosen per shell - almost none in the first, some in +// the second, enough to bury the signal in the third - so the answer is known before the run. +TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") { + Logger logger("ModelValidation_CCModelFollowsTheSignalByShell"); + + const auto path = WriteTemp("cc_model_test.pdb", ClusterPdb().c_str()); + auto obs = ModelReferenceIntensities(path, {}, {}, 2.5, logger); + REQUIRE(obs.size() > 1000); + + // The shells the correlation is reported on, coarse to fine, and the noise each one gets as a + // multiple of the r.m.s. intensity of that shell. Nothing coarser than the first shell is kept: + // the reference intensities carry the bulk solvent at fixed constants while the validation fits + // its own, and below about 6 A that difference is a large part of |F| and would decorrelate a + // shell this test needs to be clean. + const std::vector shells{5.0f, 3.2f, 2.5f}; + const double noise[3] = {0.02, 1.0, 30.0}; + + auto shell_of = [&](float d) { + for (size_t s = 0; s < shells.size(); s++) + if (d > shells[s]) return s; + return shells.size(); + }; + std::erase_if(obs, [&](const MergedReflection &r) { return r.d > 6.0f || shell_of(r.d) >= shells.size(); }); + REQUIRE(obs.size() > 500); + + std::vector sum_i2(shells.size(), 0.0); + std::vector count(shells.size(), 0); + for (const auto &r : obs) { + sum_i2[shell_of(r.d)] += static_cast(r.I) * r.I; + ++count[shell_of(r.d)]; + } + + std::mt19937 rng(20260907); + std::normal_distribution normal(0.0, 1.0); + for (size_t i = 0; i < obs.size(); i++) { + const size_t bin = shell_of(obs[i].d); + const double sd = noise[bin] * std::sqrt(sum_i2[bin] / count[bin]); + obs[i].I = static_cast(obs[i].I + sd * normal(rng)); + obs[i].sigma = static_cast(std::max(1.0, sd)); + 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() / "cc_model_test").string(); + const auto result = + ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, + .alpha = 90, .beta = 90, .gamma = 90}, + 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()); + int n_total = 0; + for (size_t s = 0; s < shells.size(); s++) { + const auto &sh = result.cc_model_shells[s]; + logger.Info("CC(model,data) {:.2f} A: {:.3f} on {} refl, {:+.1f} sigma", + sh.d_min, sh.cc, sh.n, sh.sigma); + CHECK(sh.d_min == shells[s]); + CHECK(sh.n > 20); + n_total += sh.n; + } + CHECK(n_total == result.cc_model_n); + // Essentially noiseless: the model is the data, so the correlation is high and hugely significant. + CHECK(result.cc_model_shells[0].cc > 0.9); + CHECK(result.cc_model_shells[0].sigma > 10.0); + // Noise at the shell's own r.m.s. still leaves plenty to see. + CHECK(result.cc_model_shells[1].cc > 0.25); + CHECK(result.cc_model_shells[1].sigma > 5.0); + // Buried: the shell must NOT come out significant, or the one-sided test would fire on noise. + CHECK(std::fabs(result.cc_model_shells[2].cc) < 0.15); + CHECK(std::fabs(result.cc_model_shells[2].sigma) < 4.0); + + // No shells asked for, none reported: a run that did not measure it writes no key. + const auto no_shells = + ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38, + .alpha = 90, .beta = 90, .gamma = 90}, + path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"), + false, 1, 1.0); + CHECK(no_shells.ok); + CHECK(no_shells.cc_model_shells.empty()); + + for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"}) + std::filesystem::remove(prefix + suffix); + std::filesystem::remove(path); +} diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 1a9bf2791..e93069571 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -70,7 +70,7 @@ TEST_CASE("ResultReport_Render", "[Diagnostics]") { // The stable keys a consumer greps for. The version is pinned on purpose: a key added to the // report is a contract change, and this line is where it has to be acknowledged. - CHECK(text.find("\nREPORT_VERSION= 12\n") != std::string::npos); + CHECK(text.find("\nREPORT_VERSION= 7\n") != std::string::npos); // SOHNCKE_SPACE_GROUP names the best group a chiral crystal could have, and it comes from the // space-group SEARCH. This fixture is given its group rather than searching for one, so there is // no Sohncke candidate to name and the key is absent - which is the honest behaviour and the @@ -313,6 +313,9 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { CHECK(dev_text.find("\nMODEL_FIT_NULL_REPLICATES= 5\n") != std::string::npos); } CHECK(text.find("\nMODEL_FIT_SIGMA= +4.12\n") != std::string::npos); + // A run that measured no CC(model, data) writes no key for it. + CHECK(text.find("CC_MODEL_OVERALL=") == std::string::npos); + CHECK(text.find("CC_MODEL_CONFIRMED_TO_D_MIN=") == std::string::npos); CHECK(text.find("\nMODEL_DECISIONS_TAKEN= ENANTIOMORPH\n") != std::string::npos); // The hand is ASSUMED from the model, never determined: merged intensities cannot see it. CHECK(text.find("\nSPACE_GROUP_ENANTIOMORPH= ASSUMED_FROM_MODEL\n") != std::string::npos); @@ -346,6 +349,30 @@ TEST_CASE("ResultReport_ModelValidationSection", "[Diagnostics]") { CHECK(untested_text.find("MODEL_FIT_NULL_MEAN=") == std::string::npos); 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. + 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_overall = 0.8123; + with_cc.cc_model_n = 10300; + result.model_validation = with_cc; + const auto cc_text = RenderResultReport("p", "in.h5", x, result); + 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("Read it in ONE direction only") != std::string::npos); + + // Nothing significant anywhere: the key still has to be written, saying so. + ModelValidationResult no_signal = with_cc; + for (auto &sh : no_signal.cc_model_shells) + sh.sigma = 1.0; + result.model_validation = no_signal; + 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;