From 7632bd6a6696a5dbeacd398eefb4b4d727a8830e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 10:38:11 +0100 Subject: [PATCH 1/2] 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 From 7afb92d0eb75d2645607335166ce9037f0705c77 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 14:48:00 +0100 Subject: [PATCH 2/2] RESULT-CHANGING: Refinement twin decision: a law also taken at 2-3 standard errors over its null where the merge's L-test reads twinned Rule B of the twin-gate review. The null test alone (rule A, z >= 3) is short of power on small or borderline data; where the data's own intensity statistics read twinned - an independent piece of evidence - a lead of 2 standard errors is enough. The evidence is the verdict the report already gives: TWINNING_VERDICT INDICATED, carried by the L-test (ESTIMATED_TWIN_FRACTION_FROM L_TEST), so neither a low second moment alone nor an L-test a pseudo-translation left unreadable counts. No new threshold. The caller decides what the data say (RugnuxScaleMerge, which holds the merge's twinning analysis) and passes it as RefinementSettings::twin_data_indicated; the library applies the rule (AdoptTwinLaw), since the twinned refinement it would adopt exists only inside DecideTwin. The verdict is computed by one function, TwinningVerdictKey, which the TWINNING_VERDICT key now uses as well. --mode refine has no merge statistics (its input may be amplitudes only), so there the null alone decides. REFINE_TWIN_RULE= NULL | NULL+LTEST says which rule decided, where a law was tested to UNTWINNED or TWINNED. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 12 ++++++- docs/RUGNUX_REPORT.md | 3 +- .../structure_refinement/Refinement.cpp | 30 ++++++++++++------ .../structure_refinement/Refinement.h | 10 ++++++ rugnux/ResultReport.cpp | 31 +++++++++++-------- rugnux/ResultReport.h | 4 +++ rugnux/RugnuxScaleMerge.cpp | 7 +++++ tests/RefinementTest.cpp | 12 +++++++ 8 files changed, 85 insertions(+), 24 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 7ff6b6508..097512a26 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -550,7 +550,17 @@ standard error (`REFINE_TWIN_NULL_Z`). The R-free of the twinned refinement agai 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 +taken for a twin. + +**The data's own evidence.** A lead of two to three standard errors is enough where the merge's intensity statistics +read twinned as well - the L-test itself below 0.44 (`TWINNING_VERDICT` `INDICATED` with +`ESTIMATED_TWIN_FRACTION_FROM` `L_TEST`; not where only the second moment reads low, nor where a pseudo-translation +leaves the L-test unreadable): two independent pieces of evidence, the data's and the model's, then agree. +`REFINE_TWIN_RULE` says which rule decided. The data-side test is not made a gate - it misses twins a model sees, a +pseudo-merohedral law or a near-perfect twin among them - only further evidence where the null alone is short. +`--mode refine` has no merge statistics, so there the null alone decides. + +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 every law that is a rotation of the cell's own basis vectors; a law of a centred lattice that is not (its matrix fractional in the cell's basis) can be split by it, and so can a free set `--mode refine` reads from a file. Of several laws, the one the placed model asks the largest fraction of is refined. Each of the diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index b199804dc..4cee3fb0a 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1225,7 +1225,8 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md | `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 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_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, or 2 where the data read twinned (`REFINE_TWIN_RULE`) | +| `REFINE_TWIN_RULE` | Where `UNTWINNED` or `TWINNED`, the rule that decided ([§15.12](CPU_DATA_ANALYSIS_REFINEMENT.md#1512-twinned-data)): `NULL` (the null alone, at 3 standard errors) or `NULL+LTEST` (taken at 2-3 standard errors because the merge's L-test reads twinned as well: `TWINNING_VERDICT` `INDICATED` with `ESTIMATED_TWIN_FRACTION_FROM` `L_TEST`). In `--mode refine` there are no merge statistics, and it is always `NULL` | | `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 fb319e3ae..d070849ab 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1266,8 +1266,11 @@ void MakeMaps(Context &ctx, const std::vector &merged, double // 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. +// The twinned refinement is taken only where the law beats its null on the free set by this many standard errors - +// or by TWIN_NULL_Z_DATA where the data's own intensity statistics read twinned, a second, independent piece of +// evidence (AdoptTwinLaw). constexpr double TWIN_NULL_Z = 3.0; +constexpr double TWIN_NULL_Z_DATA = 2.0; void TwinNull(const Context &ctx, RefinementResult &result) { RefinementReflections null = ctx.refl; const size_t n = null.Size(); @@ -3040,11 +3043,12 @@ 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 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. +// (into `result`) where the law beats its own null (TwinNull) on the free reflections by TWIN_NULL_Z standard errors, +// or by TWIN_NULL_Z_DATA where the caller says the data read twinned (AdoptTwinLaw) - 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, const RefinementSettings &settings, RefinementResult &result) { @@ -3103,16 +3107,20 @@ void DecideTwin(const gemmi::Structure &placed, const std::vector 0 && result.twin_null_z >= TWIN_NULL_Z; + const bool adopt = AdoptTwinLaw(twinned.twin_fraction, result.twin_null_z, settings.twin_data_indicated); result.twin_decision = twinned.twin_fraction == 0 ? "ZERO_FRACTION" : adopt ? "TWINNED" : "UNTWINNED"; + if (twinned.twin_fraction > 0) + result.twin_rule = adopt && result.twin_null_z < TWIN_NULL_Z ? "NULL+LTEST" : "NULL"; 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 -> {}", + "(mates unrelated, fraction {:.3f}) R-free {:.4f}, the law ahead of it by {:.1f} standard errors{} -> {}", result.twin_law, result.twin_fraction, result.twin_r_work, result.twin_r_free, result.untwinned_r_free, - result.twin_fraction_null, result.twin_r_free_null, result.twin_null_z, result.twin_decision); + result.twin_fraction_null, result.twin_r_free_null, result.twin_null_z, + settings.twin_data_indicated ? ", the data read twinned" : "", result.twin_decision); if (!adopt) return; // the twinned refinement is the result, with the evidence that chose it twinned.twin_decision = result.twin_decision; + twinned.twin_rule = result.twin_rule; twinned.twin_laws = result.twin_laws; twinned.twin_r_work = result.twin_r_work; twinned.twin_r_free = result.twin_r_free; @@ -3123,6 +3131,10 @@ void DecideTwin(const gemmi::Structure &placed, const std::vector 0 && (null_z >= TWIN_NULL_Z || (data_indicated && null_z >= TWIN_NULL_Z_DATA)); +} + RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, Logger &logger, const RefinementSettings &settings) { diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index cb554f74f..6205c02ca 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -75,6 +75,9 @@ struct RefinementSettings { // A refinement against the data as twinned by each twin law the lattice admits (RefinementTwinLaws), beside the // untwinned one, and the model of whichever the evidence favours (RefineModel). bool twin = true; + // The data's own intensity statistics read twinned - the caller's verdict, from the merge (the L-test): the law + // is then taken from 2 standard errors over its null instead of 3 (AdoptTwinLaw). + bool twin_data_indicated = false; bool gpu = true; size_t nthreads = 1; size_t gpu_resident_bytes = 0; // held on the card by what runs beside @@ -219,8 +222,15 @@ struct RefinementResult { double untwinned_r_work = NAN, untwinned_r_free = NAN; double twin_r_free_null = NAN, twin_fraction_null = NAN; double twin_null_z = NAN; // the law's lead over the null on the free set, in standard errors + // Where tested to UNTWINNED or TWINNED, the rule that decided: NULL (the null alone, 3 standard errors), or + // NULL+LTEST (taken at 2-3 standard errors because the data read twinned, settings.twin_data_indicated). + std::string twin_rule; }; +// Whether the twinned refinement is taken (RefineModel): a fraction above 0, and the law ahead of its null on the +// free set by 3 standard errors - or by 2 where the data's own intensity statistics read twinned as well. +bool AdoptTwinLaw(double twin_fraction, double null_z, bool data_indicated); + // `placed`: the model in the data's cell and space group. With settings.solvent (and restrained, not score_only), // the weight trajectory is followed by up to SOLVENT_BUILD_ROUNDS rounds of the solvent build, each the trajectory // again from the model with its solvent changed, kept where it ends at a lower R-free (Refinement.cpp). `merged`: the diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index bed0cc371..6cfae066d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -595,6 +595,8 @@ namespace { Add(s, KeyReal("REFINE_TWIN_R_FREE_NULL", rf.twin_r_free_null, "{:.4f}")); Add(s, KeyReal("REFINE_TWIN_FRACTION_NULL", rf.twin_fraction_null, "{:.3f}", true)); Add(s, KeyReal("REFINE_TWIN_NULL_Z", rf.twin_null_z, "{:.1f}")); + if (!rf.twin_rule.empty()) + Add(s, KeyEnum("REFINE_TWIN_RULE", rf.twin_rule, {"NULL", "NULL+LTEST"})); Add(s, KeyReal("REFINE_UNTWINNED_R_WORK", rf.untwinned_r_work, "{:.4f}")); Add(s, KeyReal("REFINE_UNTWINNED_R_FREE", rf.untwinned_r_free, "{:.4f}")); } @@ -672,6 +674,20 @@ namespace { } } +// One verdict, in the words every other place uses (TwinningVerdictLine): INDICATED is a twin in a class that +// admits one; SYMMETRY_SUSPECT is the same low reading in a holohedral class, where it says an adopted operator +// averages unequal intensities. NOT_READABLE where <|L|> reads outside the range any twin fraction gives: then +// neither twinning nor a false operator is read from it. NO_TWIN_LAW: twin-like, but the lattice admits no twin law +// for the group (see TwinningAnalysisResult::twin_like_without_twin_law) - a question about the symmetry. +std::string TwinningVerdictKey(const ProcessResult &r) { + const auto &tw = r.twinning; + return TwinStatisticsAbnormal(r) ? "NOT_READABLE" + : tw.adopted_operators_suspect ? "SYMMETRY_SUSPECT" + : tw.twinning_suspected ? "INDICATED" + : tw.twin_like_without_twin_law ? "NO_TWIN_LAW" + : "NO_INDICATION"; +} + // ---------------------------------------------------------------------------- the emitter // Populates the document. Nothing here writes to a stream: the order of the FILE is the renderer's // business, which is what lets the SUMMARY sit above the sections whose warnings it reports. @@ -2134,20 +2150,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, : "UNAFFECTED", {"UNAFFECTED", "REPAIRED", "UNREADABLE"})); Add(s, KeyReal("SECOND_MOMENT_I", tw.second_moment, "{:.4f}")); - // One verdict, in the words every other place uses (TwinningVerdictLine): INDICATED is a twin - // in a class that admits one; SYMMETRY_SUSPECT is the same low reading in a holohedral class, - // where it says an adopted operator averages unequal intensities. - // NOT_READABLE where <|L|> reads outside the range any twin fraction gives: then neither - // twinning nor a false operator is read from it. const bool abnormal = TwinStatisticsAbnormal(result); - // NO_TWIN_LAW: twin-like, but the lattice admits no twin law for the group (see - // TwinningAnalysisResult::twin_like_without_twin_law) - a question about the symmetry. - Add(s, KeyEnum("TWINNING_VERDICT", - abnormal ? "NOT_READABLE" - : tw.adopted_operators_suspect ? "SYMMETRY_SUSPECT" - : tw.twinning_suspected ? "INDICATED" - : tw.twin_like_without_twin_law ? "NO_TWIN_LAW" - : "NO_INDICATION", + // One verdict (TwinningVerdictKey), the one refinement's twin decision reads too. + Add(s, KeyEnum("TWINNING_VERDICT", TwinningVerdictKey(result), {"NO_INDICATION", "INDICATED", "SYMMETRY_SUSPECT", "NO_TWIN_LAW", "NOT_READABLE"})); // The fraction comes from the statistic that carries the verdict, and the key beside it says // which. None is quoted where it cannot be measured: merged under a false operator, or from a diff --git a/rugnux/ResultReport.h b/rugnux/ResultReport.h index 5078a22cf..01da92200 100644 --- a/rugnux/ResultReport.h +++ b/rugnux/ResultReport.h @@ -59,6 +59,10 @@ inline constexpr const char *RUGNUX_ACKNOWLEDGEMENT = // The overall CC_anom above which the summary says the data carry anomalous signal. inline constexpr double CC_ANOM_SIGNAL = 0.15; +// TWINNING_VERDICT: NO_INDICATION, INDICATED, SYMMETRY_SUSPECT, NO_TWIN_LAW or NOT_READABLE, from the merge's +// intensity statistics (result.twinning, and the space-group search's merge where one was promoted from). +std::string TwinningVerdictKey(const ProcessResult &result); + // The emitter. Populates the result object every rendering reads from; it writes no text itself. ReportDocument BuildReportDocument(const std::string &output_prefix, const std::string &input_file, diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 70863dde0..d76e6efe7 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -12,6 +12,7 @@ #include "../image_analysis/SensorAbsorption.h" #include "DiagnosticOutput.h" #include "RugnuxFirstPass.h" +#include "ResultReport.h" #include #include @@ -5169,6 +5170,12 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { settings.gpu_resident_bytes = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; // A cut (the broker's deadline or cancel, Ctrl+C) cuts a running refinement short, keeping what it has settings.cancel = cancelled_.RunFlag(); + // The data-side evidence for a twin law, from this merge: the L-test itself reads twinned - the verdict + // the report gives (TWINNING_VERDICT INDICATED), carried by the L-test (ESTIMATED_TWIN_FRACTION_FROM + // L_TEST), not by a second moment alone or an L-test a pseudo-translation left unreadable. With it the + // refinement takes a law at a smaller lead over its null (AdoptTwinLaw). + settings.twin_data_indicated = TwinningVerdictKey(result) == "INDICATED" + && result.twinning.twin_fraction_source == TwinFractionSource::LTest; const PlacedModelOutput &p = *placed_model_out; // In the group the model was validated in: its own, where that has the rotations of the group // written - the screws are then the model's, as its fit already took them, with the reflections diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index b938c8a95..d47d569c8 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -1529,6 +1529,18 @@ TEST_CASE("RefinementTwin_laws_and_mates", "[Refinement]") { CHECK(SetTwinLaw(split, t.cell, t.sg, 2.2, t.law) > 0); } +TEST_CASE("RefinementTwin_rule_with_data_evidence", "[Refinement]") { + // the null alone: 3 standard errors + CHECK(AdoptTwinLaw(0.1, 3.0, false)); + CHECK_FALSE(AdoptTwinLaw(0.1, 2.9, false)); + // the data read twinned as well: 2 standard errors + CHECK(AdoptTwinLaw(0.1, 2.0, true)); + CHECK(AdoptTwinLaw(0.1, 2.5, true)); + CHECK_FALSE(AdoptTwinLaw(0.1, 1.9, true)); + // no fraction, no twin, whatever the evidence + CHECK_FALSE(AdoptTwinLaw(0.0, 10.0, true)); +} + TEST_CASE("RefinementTwin_gradient", "[Refinement]") { TwinCase t; gemmi::Structure truth = random_crystal(t.cell, "P 4", 240, 17);