diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 84bbb1d45..80210a944 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -99,6 +100,11 @@ constexpr unsigned NULL_SEED = 20260902; // measured in. constexpr double MODEL_FIT_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. +constexpr double NULL_D_MIN = 3.5; +constexpr size_t NULL_MIN_FREE = 1000; + // Mean and sample standard deviation of a small sample. std::pair mean_sd(const std::vector &v) { if (v.size() < 2) @@ -486,9 +492,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // --- Fcalc (atomic) via electron density on a grid + FFT, plus a flat bulk-solvent mask -> Fmask. // A lambda because the rigid-body step below moves the model and then needs both again, and it // takes the state to work on so that a null replicate can run it on its own copy. --- - auto compute_model_factors = [&](ModelState &ms) { + auto compute_model_factors = [&](ModelState &ms, double to_d) { gemmi::DensityCalculator dc; - dc.d_min = d_min; + dc.d_min = to_d; dc.rate = 1.5; dc.set_grid_cell_and_spacegroup(ms.st); dc.set_refmac_compatible_blur(ms.st.models[0]); @@ -505,7 +511,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & masker.put_mask_on_grid(mask_grid, ms.st.models[0]); ms.fmask = MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0); }; - compute_model_factors(mdl); + compute_model_factors(mdl, d_min); gemmi::GroupOps gops = sg->operations(); gemmi::ReciprocalAsu asu(sg); @@ -611,7 +617,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & double rb_angle_deg = 0, rb_shift_A = 0, r_free_before_rb = 0; }; auto place_and_fit = [&](ModelState &ms, const std::vector &obs_in, - std::optional fitted) -> Placement { + std::optional fitted, double to_d) -> Placement { Placement out; out.fit = fitted ? std::move(*fitted) : fit_model(ms, obs_in); @@ -625,7 +631,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // exactly where it was read. const std::vector before = ModelPositions(ms.st.models[0]); const RigidBodyRefineResult rb = - RefineRigidBody(ms.st.models[0], ucell, *sg, out.fit.fobs_work, d_min, logger, nthreads); + RefineRigidBody(ms.st.models[0], ucell, *sg, out.fit.fobs_work, to_d, logger, nthreads); if (!rb.converged) { // RefineRigidBody leaves the model wherever the solver left it, usable answer or not, so // the restore cannot be conditional on the same flag the re-fit is. @@ -633,7 +639,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return out; } { - compute_model_factors(ms); + compute_model_factors(ms, to_d); Fit moved = fit_model(ms, obs_in); const bool commit = moved.r_free < out.fit.r_free; logger.Info("Model validation: rigid body held-out R-free {:.4f} -> {:.4f} => {}", @@ -657,13 +663,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & double margin = 0; // by how much in R-free it leads the runner-up std::optional identity_fit; // fit_model(ms, obs), for place_and_fit to reuse }; - auto probe_indexing = [&](ModelState &ms) { + auto probe_indexing = [&](ModelState &ms, const std::vector &obs_in) { IndexingProbe out; - out.identity_fit = fit_model(ms, obs); + out.identity_fit = fit_model(ms, obs_in); std::vector r_free{out.identity_fit->r_free}; // identity first, then the twin laws double best_r_free = r_free.front(); for (const auto &op : reindex_ops) { - const double cand = fit_model(ms, ReindexReflections(obs, op)).r_free; + const double cand = fit_model(ms, ReindexReflections(obs_in, op)).r_free; r_free.push_back(cand); if (cand < best_r_free) { best_r_free = cand; @@ -683,7 +689,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & std::vector reindexed; const std::vector *obs_model = &obs; if (!reindex_ops.empty()) { - indexing = probe_indexing(mdl); + indexing = probe_indexing(mdl, obs); if (!(indexing.op == gemmi::Op::identity())) { reindexed = ReindexReflections(obs, indexing.op); obs_model = &reindexed; @@ -695,7 +701,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & if (obs_model == &obs) identity_fit = std::move(indexing.identity_fit); indexing.identity_fit.reset(); - Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit)); + Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit), d_min); // --- is there anything for the model to decide? --- // There are three: the space-group label, where the model asserts the other enantiomorph; the @@ -746,6 +752,52 @@ ModelValidationResult ValidateAgainstModel(const std::vector & for (int i = 0; i < NULL_REPLICATES; i++) null_rotation.push_back(random_rotation(rng)); } + // The null and the side of the real fit it is compared with are both scored to NULL_D_MIN, + // not to the data's resolution. The placement already stops there (RigidBodyRefine's ladder + // ends at 3.5 A), and whether a model is in the right orientation, or the data in its + // indexing, is decided at low resolution: the finer shells only add reflections to the scale + // fits, which were most of the null's cost. Both sides go through the same fits on the same + // reflections, so the comparison stays like for like. On a small cell 3.5 A leaves too few + // free reflections for a real indexing margin to be told from a random placement's, so the + // limit then moves out to where NULL_MIN_FREE of them are - to the data's own on the + // smallest, which are the cheap ones. + std::vector free_d; + for (const MergedReflection &r : obs) + if (r.rfree_flag && !std::isnan(r.F)) + free_d.push_back(r.d); + std::sort(free_d.begin(), free_d.end(), std::greater<>()); + double null_d_min = std::max(d_min, NULL_D_MIN); + if (free_d.size() < NULL_MIN_FREE) + null_d_min = d_min; + else + null_d_min = std::min(null_d_min, free_d[NULL_MIN_FREE - 1]); + std::vector obs_null; + for (const MergedReflection &r : obs) + if (r.d >= null_d_min) + obs_null.push_back(r); + // The real model's side of it: the indexing margin from the model as read, as a replicate's + // is, and the R-work of the model where its placement left it, against the reflections that + // placement was fitted to. + auto coarse_r_work = [&](const std::vector &placed_against) { + ModelState coarse; + coarse.st = mdl.st; + compute_model_factors(coarse, null_d_min); + std::vector in; + for (const MergedReflection &r : placed_against) + if (r.d >= null_d_min) + in.push_back(r); + return fit_model(coarse, in).r_work; + }; + double real_margin = 0; + if (!reindex_ops.empty()) { + ModelState coarse; + coarse.st = mdl.st; + SetModelPositions(coarse.st.models[0], as_read); + compute_model_factors(coarse, null_d_min); + real_margin = probe_indexing(coarse, obs_null).margin; + } + double real_r_work = coarse_r_work(*obs_model); + // Each replicate works on its own copy of the model and of its structure factors, and the // replicates share nothing else, so they run concurrently and the real model - which the // maps and the atom-density readout below are still taken from - is never moved. @@ -755,14 +807,14 @@ ModelValidationResult ValidateAgainstModel(const std::vector & rep.st = mdl.st; SetModelPositions(rep.st.models[0], as_read); reorient_about_centroid(rep.st.models[0], null_rotation[i]); - compute_model_factors(rep); + compute_model_factors(rep, null_d_min); std::optional identity_fit; if (!reindex_ops.empty()) { - IndexingProbe probe = probe_indexing(rep); + IndexingProbe probe = probe_indexing(rep, obs_null); null_margin[i] = probe.margin; identity_fit = std::move(probe.identity_fit); // replicates are placed against obs } - null_r_work[i] = place_and_fit(rep, obs, std::move(identity_fit)).fit.r_work; + null_r_work[i] = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min).fit.r_work; }; // On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool // worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of @@ -785,7 +837,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.null_r_work_mean = null_mean; result.null_r_work_sd = null_sd; - auto record_fit = [&](const Placement &p) { + auto record_fit = [&](double r_work) { // The spread is floored, not trusted down to zero: sd == 0 and sd == 1e-7 are one ulp apart // and would give sigmas of nothing and of millions. A null with no spread - a model too // small or too symmetric for a rotation about its own centroid to move |Fcalc| - fits no @@ -793,10 +845,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // not a gate: a large, well-measured crystal has a null as narrow as 0.0006-0.0010 // (measured), and a real fit 0.3 below it is hundreds of sigma, not none. constexpr double NULL_SD_FLOOR = 1e-3; - result.r_work_sigma = (null_mean - p.fit.r_work) / std::max(null_sd, NULL_SD_FLOOR); + result.r_work_sigma = (null_mean - r_work) / std::max(null_sd, NULL_SD_FLOOR); result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA; }; - record_fit(real); + record_fit(real_r_work); // --- does the model get to decide anything? --- // The R-factors, the maps and the rigid-body placement are statements about the MODEL and are @@ -806,12 +858,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // wrong model must not be able to make. if (!reindex_ops.empty()) { result.indexing_probed = true; - result.indexing_margin = indexing.margin; + result.indexing_margin = real_margin; const auto [margin_mean, margin_sd] = mean_sd(null_margin); result.indexing_margin_null_mean = margin_mean; result.indexing_margin_null_sd = margin_sd; constexpr double MARGIN_SD_FLOOR = 1e-4; - result.indexing_margin_sigma = (indexing.margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR); + result.indexing_margin_sigma = (real_margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR); // 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. @@ -820,9 +872,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector & if (decided) result.indexing_op = indexing.op; logger.Info("Model validation: probed {} indexing solution(s) against the model; the winner " - "leads the runner-up by {:.4f} in R-free, against {:.4f} +- {:.4f} for a random " - "placement of the same model ({:+.2f} sigma) => {}", - reindex_ops.size() + 1, indexing.margin, margin_mean, margin_sd, + "leads the runner-up by {:.4f} in R-free to {:.2f} A, against {:.4f} +- {:.4f} for a " + "random placement of the same model ({:+.2f} sigma) => {}", + reindex_ops.size() + 1, real_margin, null_d_min, margin_mean, margin_sd, result.indexing_margin_sigma, !decided ? "not decided; the data keep the indexing they were merged in" : (result.indexing_op == gemmi::Op::identity() @@ -831,14 +883,15 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // The R-factors and the maps have to describe the reflections the files carry, and those // are now the ones the merge produced, so the fit is remade on them. SetModelPositions(st.models[0], as_read); - compute_model_factors(mdl); - real = place_and_fit(mdl, obs, std::nullopt); - record_fit(real); + compute_model_factors(mdl, d_min); + real = place_and_fit(mdl, obs, std::nullopt, d_min); + real_r_work = coarse_r_work(obs); + record_fit(real_r_work); } } - logger.Info("Model validation: R-work {:.4f} against a null of {:.4f} +- {:.4f} = {:+.2f} sigma " - "=> the model {} these data", - real.fit.r_work, null_mean, null_sd, result.r_work_sigma, + logger.Info("Model validation: R-work {:.4f} to {:.2f} A against a null of {:.4f} +- {:.4f} = " + "{:+.2f} sigma => the model {} these data", + real_r_work, null_d_min, null_mean, null_sd, result.r_work_sigma, result.model_fits ? "FITS" : "DOES NOT FIT"); } else { logger.Info("Model validation: the model asserts no other enantiomorph for these data and " @@ -861,8 +914,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.indexing_op = gemmi::Op::identity(); result.r_work_sigma = 0.0; SetModelPositions(st.models[0], as_read); - compute_model_factors(mdl); - real = place_and_fit(mdl, obs, std::nullopt); + compute_model_factors(mdl, d_min); + real = place_and_fit(mdl, obs, std::nullopt, d_min); } Fit &best = real.fit;