From 7632bd6a6696a5dbeacd398eefb4b4d727a8830e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 10:38:11 +0100 Subject: [PATCH] RESULT-CHANGING: Refinement twin decision: the null's standard error per twin pair, its mates from the same set, no R-free clause TwinNull: the paired difference |Fo-|F_null|| - |Fo-|F_law|| is summed over each free twin pair before its mean and variance are taken - both members share their two F_calc, so they were counted as two independent terms and the standard error was too small by up to sqrt(2). The null's unrelated mates are drawn from the reflection's own set (free with free, work with work), as a law's twin pairs are; each set is ordered by 1/d^2 with a twin pair's members adjacent and swapped 2 apart, past the reflection's own pair, because on the sparse free set the old step of 8 took the mate to another resolution and made the null worse than it is (measured on the synthetic untwinned case: null R-free 0.048 -> 0.065, z -0.5 -> 2.4 with a step of 8 within the free set; 0.060, 1.3 with the step of 2). DecideTwin: the clause R_free(twin) < R_free(untwinned) - 0.001 is dropped. It compares two refinement paths and never bound on the rc181 battery (smallest margin among z >= 3: 0.0029); the null test at z >= 3 is the decision. The H-test line (TWIN_LAW_n, key unchanged) no longer says "twin fraction X or more": untwinned crystals read 0.5 - H up to about 0.25, so it is stated as the fraction only where twinning alone makes the pairs agree. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 18 +++-- docs/RUGNUX_REPORT.md | 4 +- .../structure_refinement/Refinement.cpp | 75 ++++++++++++------- rugnux/RugnuxScaleMerge.cpp | 8 +- tests/RefinementTest.cpp | 2 + 5 files changed, 66 insertions(+), 41 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index df3a89029..7ff6b6508 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -539,14 +539,16 @@ is linear in the two intensities, so an unrelated partner takes a coefficient of least-squares fit of the twinned amplitudes instead took a fraction of 0.16 for unrelated partners of one model at R-free 0.31, the average of two predictions being a better prediction than either. -**The decision.** The twinned refinement is kept where its fraction is above 0, its R-free is below the -untwinned refinement's by more than 0.001, and the law beats its own **null** on the free reflections by at least -three standard errors - the null being the same twinned model at the same fraction with every reflection's mate -replaced by an unrelated reflection of the same resolution: what the averaging of two intensities buys whatever -they are (it lifts a weak reflection's model amplitude off zero). The comparison is paired, reflection by reflection, $|F_o - |F_\mathrm{null}|| - |F_o - |F_\mathrm{twin}||$ -over the free set, its mean against its standard error (`REFINE_TWIN_NULL_Z`): two refinements of one model -differ in R-free by a few thousandths from their trajectories alone, more than a fixed tie, and a fraction of one -or two percent is not resolved by them. A pseudo-symmetry the model already carries (non-crystallographic +**The decision.** The twinned refinement is kept where its fraction is above 0 and the law beats its own **null** +on the free reflections by at least three standard errors - the null being the same twinned model at the same +fraction with every reflection's mate replaced by an unrelated reflection of the same resolution and of the same +set (a free reflection's by a free one, as a law's are): what the averaging of two intensities buys whatever they +are (it lifts a weak reflection's model amplitude off zero). The comparison is paired, reflection by reflection, +$|F_o - |F_\mathrm{null}|| - |F_o - |F_\mathrm{twin}||$ over the free set, summed over each free twin pair (the two +mates share their two $F_\mathrm{calc}$, so they are one term, not two independent ones), its mean against its +standard error (`REFINE_TWIN_NULL_Z`). The R-free of the twinned refinement against the untwinned one's does not +enter: two refinements of one model differ in R-free by a few thousandths from their trajectories alone, and any +twin target lowers R by its averaging, so that comparison favours the twinned one with or without a twin. A pseudo-symmetry the model already carries (non-crystallographic symmetry parallel to the law) predicts $I(T\mathbf{h}) \approx I(\mathbf{h})$, so the law changes nothing and is not taken for a twin. A law is tested only where the free set keeps its twin pairs together, and is otherwise reported (`REFINE_TWIN` `FREE_SET_SPLIT`) rather than tested. The free set rugnux draws (§10.7) keeps together the pairs of diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 9dba7759b..b199804dc 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1224,8 +1224,8 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md | `REFINE_TWIN_LAW`, `REFINE_TWIN_FRACTION` | The twin law tested (an hkl operator in the data's basis) and the fraction it refined to | | `REFINE_TWIN_SPLIT_PAIRS` | Where `FREE_SET_SPLIT`: the twin pairs with one member in the working and one in the free set | | `REFINE_TWIN_R_WORK`, `REFINE_TWIN_R_FREE`, `REFINE_UNTWINNED_R_WORK`, `REFINE_UNTWINNED_R_FREE` | The twinned and the untwinned refinement, each where it ended (its solvent built) | -| `REFINE_TWIN_R_FREE_NULL` | The twinned model's R-free at its fraction with every mate replaced by an unrelated reflection of the same resolution: what averaging alone buys. `REFINE_TWIN_FRACTION_NULL` (*dev*) is the fraction those unrelated mates ask for, about 0 | -| `REFINE_TWIN_NULL_Z` | The law's lead over that null on the free reflections, paired reflection by reflection, in standard errors; `TWINNED` needs 3 | +| `REFINE_TWIN_R_FREE_NULL` | The twinned model's R-free at its fraction with every mate replaced by an unrelated reflection of the same resolution and set: what averaging alone buys, and the like-for-like comparison for `REFINE_TWIN_R_FREE`. `REFINE_TWIN_FRACTION_NULL` (*dev*) is the fraction those unrelated mates ask for, about 0 | +| `REFINE_TWIN_NULL_Z` | The law's lead over that null on the free reflections, paired reflection by reflection (a free twin pair as one term), in standard errors; `TWINNED` needs 3 | | `REFINE_BLOBS` | The blobs of the refined mFo-DFc map (above 3 σ, peaking above 4 σ) | | `REFINE_BLOBS_LIGAND`, `REFINE_BLOBS_WATER`, `REFINE_BLOBS_ION`, `REFINE_BLOBS_MODEL_ERROR` | How many of them are read as each ([§15.9](CPU_DATA_ANALYSIS_REFINEMENT.md#159-unexplained-density-blobs)); the rest are unassigned | | `REFINE_SOLVENT_SITES` | The waters and single-atom ions of the refined model the solvent audit scored ([§15.10](CPU_DATA_ANALYSIS_REFINEMENT.md#1510-the-solvent-audit)) | diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 11ea25145..fb319e3ae 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1260,24 +1260,35 @@ void MakeMaps(Context &ctx, const std::vector &merged, double } // The null of a twinned refinement, for the model as the last Refit() left it: its R-free at the same fraction with -// every reflection's mate replaced by an unrelated observed one of the same resolution - the reflections in order of -// 1/d^2 paired TWIN_NULL_SHIFT apart, beyond any run of equal 1/d^2 a twin law relates. What the averaging of two -// intensities buys on its own, whatever they are. -constexpr size_t TWIN_NULL_SHIFT = 8; +// every reflection's mate replaced by an unrelated observed one of the same resolution and the same set (free with +// free, work with work, as a law's twin pairs are). The reflections of each set are put in order of 1/d^2 with the two +// members of a twin pair next to each other, and paired TWIN_NULL_SHIFT apart - past the reflection's own pair, and no +// further, as the free set is sparse and a wider step takes the mate to another resolution, which a null pays for. +// What the averaging of two intensities buys on its own, whatever they are. +constexpr size_t TWIN_NULL_SHIFT = 2; // The twinned refinement is taken only where the law beats its null on the free set by this many standard errors. constexpr double TWIN_NULL_Z = 3.0; void TwinNull(const Context &ctx, RefinementResult &result) { RefinementReflections null = ctx.refl; const size_t n = null.Size(); null.mate_hkl.clear(); - std::vector order(n); - std::iota(order.begin(), order.end(), 0); - std::stable_sort(order.begin(), order.end(), [&](int a, int b) { return null.inv_d2[a] < null.inv_d2[b]; }); null.mate.assign(n, 0); - for (size_t r = 0; r < n; r++) { - // in blocks of twice the shift, the first half with the second - const size_t p = r % (2 * TWIN_NULL_SHIFT) < TWIN_NULL_SHIFT ? r + TWIN_NULL_SHIFT : r - TWIN_NULL_SHIFT; - null.mate[order[r]] = order[p < n ? p : r]; + for (const bool free : {false, true}) { + std::vector order; + for (size_t m = 0; m < n; m++) + if ((null.free[m] != 0) == free) + order.push_back(static_cast(m)); + // a twin pair is ordered at the 1/d^2 of its first member + auto pair_of = [&](int m) { return ctx.refl.mate[m] < static_cast(n) ? std::min(m, ctx.refl.mate[m]) : m; }; + std::stable_sort(order.begin(), order.end(), [&](int a, int b) { + const int pa = pair_of(a), pb = pair_of(b); + return std::make_pair(null.inv_d2[pa], pa) < std::make_pair(null.inv_d2[pb], pb); + }); + for (size_t r = 0; r < order.size(); r++) { + // in blocks of twice the shift, the first half with the second + const size_t p = r % (2 * TWIN_NULL_SHIFT) < TWIN_NULL_SHIFT ? r + TWIN_NULL_SHIFT : r - TWIN_NULL_SHIFT; + null.mate[order[r]] = order[p < order.size() ? p : r]; + } } // at the law's fraction: averaging any two intensities lifts the weak reflections' model amplitudes off zero // by itself (the fraction the unrelated mates ask for, reported beside it, is about 0) @@ -1286,20 +1297,27 @@ void TwinNull(const Context &ctx, RefinementResult &result) { const std::vector> f_null = TwinnedModel(null, rows, ctx.fcalc); result.twin_r_free_null = ComputeRFactors(null, f_null).r_free; // How far the law beats its null on the free reflections, paired: the mean of |Fo - |F_null|| - |Fo - |F_twin|| - // over its standard error. + // over its standard error. The two members of a free twin pair share their two F_calc, so their terms are summed + // into one: the standard error counts twin pairs, not reflections. const std::vector> f_law = TwinnedModel(ctx.refl, ctx.rows, ctx.fcalc); - double sum = 0, sum2 = 0; - int n_free = 0; + std::vector d(n, 0.0); for (size_t m = 0; m < n; m++) - if (ctx.refl.free[m]) { - const double d = std::fabs(ctx.refl.fo[m] - std::abs(f_null[m])) - std::fabs(ctx.refl.fo[m] - std::abs(f_law[m])); - sum += d; - sum2 += d * d; - n_free++; - } - if (n_free > 1) { - const double mean = sum / n_free, var = (sum2 - n_free * mean * mean) / (n_free - 1); - result.twin_null_z = var > 0 ? mean / std::sqrt(var / n_free) : 0.0; + if (ctx.refl.free[m]) + d[m] = std::fabs(ctx.refl.fo[m] - std::abs(f_null[m])) - std::fabs(ctx.refl.fo[m] - std::abs(f_law[m])); + double sum = 0, sum2 = 0; + int n_pairs = 0; + for (size_t m = 0; m < n; m++) { + const size_t t = ctx.refl.mate[m]; + if (!ctx.refl.free[m] || (t < n && t < m)) + continue; // the pair was counted at its first member + const double pair = t < n && t != m ? d[m] + d[t] : d[m]; + sum += pair; + sum2 += pair * pair; + n_pairs++; + } + if (n_pairs > 1) { + const double mean = sum / n_pairs, var = (sum2 - n_pairs * mean * mean) / (n_pairs - 1); + result.twin_null_z = var > 0 ? mean / std::sqrt(var / n_pairs) : 0.0; } } @@ -3022,10 +3040,10 @@ double FinalRWork(const RefinementResult &r) { // The twin hypothesis beside the untwinned refinement `result`: the placed model refined again, solvent build and // displacement model included, against the data as twinned by a law of the lattice, and the twinned refinement taken -// (into `result`) where its R-free is below the untwinned one's and the law beats its own null (TwinNull) on the free -// reflections by TWIN_NULL_Z standard errors - so that neither the drop in R any twin target buys by averaging, nor a -// pseudo-symmetry the model already carries (it predicts I(Th) ~ I(h), and the law then changes nothing) is read as -// a twin. A law is tested only where the free set keeps its twin pairs together; of several, the one the placed +// (into `result`) where the law beats its own null (TwinNull) on the free reflections by TWIN_NULL_Z standard errors - +// not on its R-free against the untwinned one's, which compares two refinement paths and so carries their noise - +// so that neither the drop in R any twin target buys by averaging, nor a pseudo-symmetry the model already carries +// (it predicts I(Th) ~ I(h), and the law then changes nothing) is read as a twin. A law is tested only where the free set keeps its twin pairs together; of several, the one the placed // model asks the largest fraction of. void DecideTwin(const gemmi::Structure &placed, const std::vector &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, Logger &logger, @@ -3085,8 +3103,7 @@ void DecideTwin(const gemmi::Structure &placed, const std::vector 0 && result.twin_r_free < result.untwinned_r_free - R_FREE_TIE && - result.twin_null_z >= TWIN_NULL_Z; + const bool adopt = twinned.twin_fraction > 0 && result.twin_null_z >= TWIN_NULL_Z; result.twin_decision = twinned.twin_fraction == 0 ? "ZERO_FRACTION" : adopt ? "TWINNED" : "UNTWINNED"; logger.Info("Refinement, twin law {}: fraction {:.3f}, R-work {:.4f}, R-free {:.4f}; untwinned R-free {:.4f}; null " "(mates unrelated, fraction {:.3f}) R-free {:.4f}, the law ahead of it by {:.1f} standard errors -> {}", diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 685d48615..70863dde0 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -4262,7 +4262,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { laws += (laws.empty() ? "" : " ; ") + t; // What each law reads on the merge: the disagreement H of the reflections it // relates (Yeates (1997) Methods Enzymol. 276, 344-358: a twin of fraction a reads - // (1-2a)/2, so a >= 0.5 - H), their correlation and their intensity-weighted R. + // (1-2a)/2), their correlation and their intensity-weighted R. 0.5 - H is the + // fraction only where nothing but twinning makes the pairs agree: the median over + // the stronger pairs and a pseudo-symmetry along the law both lower H, and + // untwinned crystals read 0.5 - H up to about 0.25 (measured over the battery). // Report-only, so a user sees which law a twin refinement would need. gemmi::Mat33 m; for (int i = 0; i < 3; ++i) @@ -4274,7 +4277,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { law_correlation.emplace(expanded_merge(), law_opts); if (const auto s = (*law_correlation)(m); s.has_value()) result.twin_law_evidence.push_back(fmt::format( - "{}: {} pairs, H {:.3f} (twin fraction {:.2f} or more), R {:.3f}, CC {:.3f}", + "{}: {} pairs, H {:.3f} (twin fraction {:.2f} if only twinning makes the pairs " + "agree; untwinned crystals read up to about 0.25), R {:.3f}, CC {:.3f}", t, s->n_pairs, s->h_stat, std::max(0.0, 0.5 - s->h_stat), s->r_stat, s->cc)); else result.twin_law_evidence.push_back(fmt::format("{}: too few pairs to read", t)); diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index 215da9d20..b938c8a95 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -1625,7 +1625,9 @@ TEST_CASE("Refinement_twin_decided", "[Refinement]") { CHECK(r.twin_fraction == Catch::Approx(fraction).margin(0.03)); CHECK(r.r_free < r.untwinned_r_free - 0.02); CHECK(r.twin_r_free < r.twin_r_free_null - 0.02); + CHECK(r.twin_null_z >= 3.0); } else { + CHECK((r.twin_decision == "ZERO_FRACTION" || r.twin_null_z < 3.0)); CHECK(r.twin_decision != "TWINNED"); CHECK_FALSE(r.twinned); // the untwinned refinement, as it is without the twin hypothesis