From b2c0a5ea48640dd8fd0153a9e80abb3ff964157c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 2 Sep 2026 15:18:55 +0200 Subject: [PATCH] model validation: review fixes - the run's output, the origin gauge, and two claims From two independent code reviews of the model-validation work. Nothing here changes a verdict: the acceptance set still reads +17.69 / +1.70 / -0.37 sigma and the rejected runs still write files byte-identical to a run with no model. THE RUN'S OUTPUT. Model validation runs BEFORE the reflection files are written, and two paths through it could throw: the null's replicates (rotated models fed to a scaling path that fails outright on data it cannot pair up - the adversarial input for it), and the mmCIF coordinate writer. Either would have taken the .mtz, .cif, .hkl and _unmerged.mtz with it, after the merge had already been paid for. A null that cannot be built is a question that could not be put, which is the NOT_TESTED state this design already has; a coordinate file that cannot be written is a lost convenience. Both now degrade instead of aborting. THE ORIGIN GAUGE. Translating the whole cell content along a free-origin direction - all three in P1, the unique axis in a polar group - leaves every |F| exactly unchanged. The code said the LM damping and the R-free gate made that harmless between them. Neither does: the gate is a function of |F| and is blind to exactly this, and the gauge column of the Jacobian is not zero but noise divided by the difference step. It is now projected out after every zone, against the group's own common fixed subspace. P2_1 alone is a large share of deposited structures, and the reported shift was partly fiction in every one of them. TWO CLAIMS THAT WERE FALSE. The report told the user R-work carries the decision "because nothing was refined against it", six lines from where six placement parameters are refined against it; the real argument is that the null is placed the same way, so the optimism is common-mode and cancels. And the constant's own comment quoted a +4 vs +1 sigma gap where the measurement is +17.7 vs +1.7. Also: the map file's phase columns are back to [0, 360), the convention they carried before sigma_A weighting; with no data space group the map file follows the reflections into P1 rather than taking the model's group; the sigma is floored against a near-zero null spread rather than only an exactly-zero one; the model is restored whether or not the solver reported a usable answer; and the zone list is the ladder walked rather than the ladder planned. Co-Authored-By: Claude Opus 5 (1M context) --- docs/RUGNUX_REPORT.md | 6 ++-- rugnux/ModelValidation.cpp | 64 ++++++++++++++++++++++++++++++++------ rugnux/ResultReport.cpp | 5 +-- rugnux/RigidBodyRefine.cpp | 63 ++++++++++++++++++++++++++++--------- rugnux/WriteModel.cpp | 19 +++++++---- 5 files changed, 122 insertions(+), 35 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 4b2c89dac..f24ccc0b3 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -201,8 +201,10 @@ There is no threshold on R behind it. What a model that explains nothing reaches set depends on its atom count and B-factors as much as on the data, so the same model is refitted — and re-placed as a rigid body, exactly as the real one is — from `MODEL_FIT_NULL_REPLICATES` random orientations about its own centroid, and `MODEL_FIT_SIGMA` is how far the real fit sits above that -distribution. The statistic is R-work, not R-free: nothing is refined against the working set here, so -R-work carries no optimism, and it is decided on an order of magnitude more reflections. +distribution. The statistic is R-work, not R-free - not because nothing is refined against the working +set (the placement's six parameters are), but because every null replicate is placed the same way, so +what they buy is bought on both sides and cancels; and it is decided on an order of magnitude more +reflections than R-free. That null is **only built where the model claims one of the two things it could change** — the enantiomorph, or an indexing other than the one the data were merged in. A model already in the data's diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index d1bc187a4..8548eaf3a 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -88,8 +88,9 @@ constexpr unsigned NULL_SEED = 20260902; // How far above its own null a fit has to sit before the model is allowed to decide anything. The cut // is in sigma of that null and not in R: measured, the R a model that explains nothing reaches moves // with the model's atom count and B-factors as much as with the data, so no value of R separates the -// two on its own. A model that belongs to the crystal measures around +4 sigma here and an unrelated -// one under +1. +// 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; // Mean and sample standard deviation of a small sample. @@ -366,7 +367,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::vector before = ModelPositions(ms.st.models[0]); const RigidBodyRefineResult rb = RefineRigidBody(ms.st.models[0], ucell, *sg, fobs_work, d_min, logger); - if (rb.converged) { + 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. + SetModelPositions(ms.st.models[0], before); + return out; + } + { compute_model_factors(ms); Fit moved = fit_model(ms, obs_in); const bool commit = moved.r_free < out.fit.r_free; @@ -451,6 +458,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // computed and reported exactly as in every other case. const bool decision_pending = result.model_enantiomorph_candidate || !(indexing.op == gemmi::Op::identity()); + // A null that cannot be built is a question that could not be put, and this design already has a + // state for that. Without this the run dies here - model validation runs BEFORE the reflection + // files are written, so one failed replicate would take the .mtz, .cif, .hkl and _unmerged.mtz + // with it, after the merge has already been paid for. The replicates are rotated models fed to a + // scaling path that throws on data it cannot pair up, so this is the adversarial input for it. + try { if (decision_pending) { // --- the null: this model, this data, in random orientations --- // An R-factor on its own says nothing about whether a model belongs to a crystal. What a model @@ -461,9 +474,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // reoriented about its own centroid and run through the identical fit and rigid-body placement, // and the real fit is asked how far above the resulting distribution it sits. // - // R-work, not R-free: nothing is refined against the working set here - the scale has four - // parameters and the placement six - so R-work carries no optimism, and it is decided on an order - // of magnitude more reflections than R-free. + // R-work, not R-free. Not because nothing is refined against it - the placement's six parameters + // are, and the scale's four - but because every null replicate is placed and scaled the same + // way, so whatever that optimism is worth is bought on both sides and cancels in + // (mean - real)/sd. R-work is then decided on an order of magnitude more reflections. // // The replicates are fitted to the data as merged even where a reindexing won above: a reindexing // is a relabelling of the same intensities, and a model in a random orientation has no more to do @@ -501,7 +515,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.null_r_work_sd = null_sd; auto record_fit = [&](const Placement &p) { - result.r_work_sigma = null_sd > 0 ? (null_mean - p.fit.r_work) / null_sd : 0.0; + // Guarded at a floor, not at zero: sd == 0 and sd == 1e-7 are one ulp apart and land on + // opposite verdicts, and a null with no spread - a model too small or too symmetric for a + // rotation about its own centroid to move |Fcalc| - is the case that produces it. Measured + // nulls sit near 0.60 with sd ~0.009, so this is nowhere near anything genuine. + constexpr double NULL_SD_FLOOR = 1e-3; + result.r_work_sigma = + null_sd >= NULL_SD_FLOOR ? (null_mean - p.fit.r_work) / null_sd : 0.0; result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA; }; record_fit(real); @@ -518,8 +538,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector & 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 = - margin_sd > 0 ? (indexing.margin - margin_mean) / margin_sd : 0.0; + margin_sd >= MARGIN_SD_FLOOR ? (indexing.margin - margin_mean) / margin_sd : 0.0; // 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. @@ -553,6 +574,23 @@ ModelValidationResult ValidateAgainstModel(const std::vector & "reflections and no null was run; R-work {:.4f}, R-free {:.4f}", real.fit.r_work, real.fit.r_free); } + } catch (const std::exception &e) { + // The null could not be built, so no decision is licensed and none is taken - which is the + // NOT_TESTED state, not a failure of the run. If the indexing probe had already won on a + // relabelling, the fit that survives describes those reflections and the files will not carry + // them, so it is remade on the data as merged. + logger.Warning("Model validation: the null could not be built ({}), so the model decides " + "nothing; the reflections keep the group and the indexing they were merged in", + e.what()); + result.fit_tested = false; + result.model_fits = false; + result.indexing_probed = false; + 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); + } Fit &best = real.fit; result.rigid_body_applied = real.rb_applied; @@ -739,7 +777,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // of this file works in the wrong one if the label disagrees with the .mtz beside it. mtz.spacegroup = result.adopted_model_enantiomorph ? sg - : (data_space_group != nullptr ? data_space_group : sg); + // No data group means the caller merged in P1, and P1 is what the + // reflection files beside this one carry - not the model's group. + : (data_space_group != nullptr ? data_space_group + : gemmi::find_spacegroup_by_number(1)); mtz.set_cell_for_all(ucell); mtz.add_dataset("model_validation"); mtz.datasets.back().wavelength = wavelength_A; @@ -756,7 +797,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & mtz.add_column("FREE", 'I', -1, -1, false); std::vector data; for (size_t i = 0; i < terms.size(); i++) { - const auto phi_deg = static_cast(terms[i].phi * 180.0 / PI); + // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg + // returns (-180, 180] and a P column is not supposed to. + const double phi_wrapped = terms[i].phi * 180.0 / PI; + const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), static_cast(terms[i].hkl[2]), static_cast(terms[i].fo), static_cast(terms[i].fc), diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 6d7784ae7..16d4dc502 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -663,8 +663,9 @@ std::string RenderResultReport(const std::string &output_prefix, << " depends on its atom count and B-factors as much as on the data - so the same model was\n" << " refitted, and re-placed, from " << mv.null_replicates << " random orientations about its own centroid, and\n" << " MODEL_FIT_SIGMA is how far the real fit sits above that null. R-work carries the\n" - << " decision because nothing was refined against it and it has far more reflections than\n" - << " R-free. MODEL_DECISIONS_TAKEN names what the model was allowed to change about the\n" + << " decision because the null was placed the same way the real fit was, so what those six\n" + << " placement parameters buy is bought on both sides and cancels - and it has far more\n" + << " reflections than R-free. MODEL_DECISIONS_TAKEN names what the model was allowed to change about the\n" << " written reflections; R-factors, maps and the rigid-body placement are reported either\n" << " way, because they describe the model, not the data.\n"; else diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index 29f4fcc99..a2a725c55 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -43,8 +43,8 @@ constexpr double JACOBIAN_STEP_FRACTION = 0.01; // Parameters are carried as six lengths in angstroms - the first three are the angle-axis rotation // vector multiplied by the model's rms radius, so a unit of each of the six moves a typical atom by -// the same amount. That makes the Jacobian step and the trust region isotropic in something -// physical, rather than mixing radians with angstroms. +// the same amount. That makes the Jacobian step isotropic in something physical, rather than mixing +// radians with angstroms. struct Placement { gemmi::Position centre; // the model centroid: rotating about it decorrelates R from t double rms_radius = 1.0; // rms distance of the atoms from the centroid @@ -185,6 +185,28 @@ private: Evaluator &ev_; }; +// The directions in which this space group's origin is free. Translating the whole cell content along +// one of them multiplies every F by a phase and leaves every |F| EXACTLY unchanged, so the target +// cannot determine that component: all three directions in P1, the unique axis in a polar group. The +// R-free gate cannot stand in for this - it is a function of |F| too, so along such a direction it +// sees only grid noise and commits or not by coin flip, while the other five parameters carry the +// noise in with them. The free directions are the common fixed subspace of the group's rotation +// parts, and the projector onto it is simply their average. +gemmi::Mat33 GaugeProjector(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) { + const gemmi::GroupOps gops = sg.operations(); + double m[3][3] = {}; + const double n = static_cast(gops.sym_ops.size()) * gemmi::Op::DEN; + for (const gemmi::Op &op : gops.sym_ops) + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + m[i][j] += static_cast(op.rot[i][j]) / n; + const gemmi::Mat33 mean(m[0][0], m[0][1], m[0][2], + m[1][0], m[1][1], m[1][2], + m[2][0], m[2][1], m[2][2]); + // Fractional projector taken into orthogonal space, where the parameters live. + return cell.orth.mat.multiply(mean).multiply(cell.frac.mat); +} + } // namespace std::vector ModelPositions(const gemmi::Model &model) { @@ -209,11 +231,12 @@ void SetModelPositions(gemmi::Model &model, const std::vector & // bound ligand its own six parameters would refine against evidence this data does not separately // carry, and the ligand is what the difference map is meant to show rather than model away. // -// Some of the translation can be a gauge rather than a quantity - the origin is free in all three +// Some of the translation would be a gauge rather than a quantity - the origin is free in all three // directions in P1 and along the unique axis in a polar group, and |F| does not change when the whole -// content moves along it. Nothing is done about that beyond the Levenberg-Marquardt damping, which -// leaves an undetermined direction where it started, and the caller's R-free gate, which throws away -// a step that wandered along one. +// content moves along it - so that component is projected out after every zone. Neither of the two +// things that might look like they cover it actually does: the R-free gate is a function of |F| and +// therefore blind to exactly this, and the LM damping follows the gauge column of the Jacobian, which +// is not zero but noise divided by the difference step. RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, @@ -238,16 +261,18 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, if (!(placement.rms_radius > 0)) return result; + std::vector ladder; for (double zone : LADDER) if (zone >= d_min) - result.zones.push_back(zone); - if (result.zones.empty()) - result.zones.push_back(d_min); + ladder.push_back(zone); + if (ladder.empty()) + ladder.push_back(d_min); + const gemmi::Mat33 gauge = GaugeProjector(sg, cell); Evaluator ev(model, cell, sg, base, placement); double q[6] = {0, 0, 0, 0, 0, 0}; bool any_zone_solved = false; - for (double zone : result.zones) { + for (double zone : ladder) { gemmi::AsuData> zone_obs; zone_obs.unit_cell_ = fobs.unit_cell_; zone_obs.spacegroup_ = fobs.spacegroup_; @@ -256,6 +281,7 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, zone_obs.v.push_back(hv); if (zone_obs.v.size() < 50) continue; + result.zones.push_back(zone); // the ladder WALKED, which a thin zone drops out of ev.SetZone(zone_obs, zone); ceres::Problem problem; @@ -271,9 +297,12 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, const auto zone_t0 = std::chrono::steady_clock::now(); ceres::Solve(options, &problem, &summary); any_zone_solved = any_zone_solved || summary.IsSolutionUsable(); + const gemmi::Vec3 along = gauge.multiply(gemmi::Vec3(q[3], q[4], q[5])); + q[3] -= along.x; q[4] -= along.y; q[5] -= along.z; logger.Debug("Rigid body zone {:.1f} A: {} reflections, {} iterations, {} evaluations, {:.2f} s, " "rotation {:.3f} deg, translation {:.3f} A", zone, zone_obs.v.size(), - summary.iterations.size() - 1, ev.evaluations - evaluations_before, + summary.iterations.empty() ? 0 : summary.iterations.size() - 1, + ev.evaluations - evaluations_before, std::chrono::duration(std::chrono::steady_clock::now() - zone_t0).count(), std::sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / placement.rms_radius * 180.0 / PI, std::sqrt(q[3]*q[3] + q[4]*q[4] + q[5]*q[5])); @@ -287,9 +316,13 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, result.shift_A = std::sqrt(q[3] * q[3] + q[4] * q[4] + q[5] * q[5]); result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - logger.Info("Model validation: rigid body over {} resolution zone(s) down to {:.1f} A, " - "{} evaluations in {:.2f} s: rotation {:.3f} deg, translation {:.3f} A", - result.zones.size(), result.zones.back(), result.evaluations, result.seconds, - result.angle_deg, result.shift_A); + if (!result.zones.empty()) + logger.Info("Model validation: rigid body over {} resolution zone(s) down to {:.1f} A, " + "{} evaluations in {:.2f} s: rotation {:.3f} deg, translation {:.3f} A", + result.zones.size(), result.zones.back(), result.evaluations, result.seconds, + result.angle_deg, result.shift_A); + else + logger.Info("Model validation: rigid body had no resolution zone with enough reflections to " + "run in; the model is left where it arrived"); return result; } diff --git a/rugnux/WriteModel.cpp b/rugnux/WriteModel.cpp index b9e36afbf..e3719bb8b 100644 --- a/rugnux/WriteModel.cpp +++ b/rugnux/WriteModel.cpp @@ -29,12 +29,19 @@ void WritePlacedModel(const gemmi::Structure &placed, st.setup_cell_images(); // Fills in entity types and label_asym_id for a model read from a PDB, which carries neither. // Both are no-ops where the input already had them, i.e. for an mmCIF input. - gemmi::setup_entities(st); - - std::ofstream os(path); - gemmi::cif::write_cif_to_stream(os, gemmi::make_mmcif_document(st)); - if (!os) { - logger.Error("Model validation: cannot write the placed model to {}", path); + // setup_entities and make_mmcif_document can both throw, and this runs BEFORE the reflection files + // are written: a convenience deliverable must not be able to take the run's actual output with it. + // The map writer does the same. + try { + gemmi::setup_entities(st); + std::ofstream os(path); + gemmi::cif::write_cif_to_stream(os, gemmi::make_mmcif_document(st)); + if (!os) { + logger.Error("Model validation: cannot write the placed model to {}", path); + return; + } + } catch (const std::exception &e) { + logger.Warning("Model validation: could not write the placed model to {} ({})", path, e.what()); return; } logger.Info("Model validation: the model as placed against these data written to {} "