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