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

model validation: the model says, shell by shell, where the data still carry signal

Browse Source 
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
leonarski_f
2026-09-07 13:12:32 +02:00
co-authored by Claude Opus 5
parent cb08f63a52
commit 004fa5a781
10 changed files with 322 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/ModelValidationTest.cpp
+87
View File
@@ -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<float> 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<double> sum_i2(shells.size(), 0.0);
std::vector<int> count(shells.size(), 0);
for (const auto &r : obs) {
sum_i2[shell_of(r.d)] += static_cast<double>(r.I) * r.I;
++count[shell_of(r.d)];
}
std::mt19937 rng(20260907);
std::normal_distribution<double> 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<float>(obs[i].I + sd * normal(rng));
obs[i].sigma = static_cast<float>(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);
}