diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index df3a89029..097512a26 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -539,16 +539,28 @@ 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 +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 e3b9ceb6f..fa065cd62 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1251,8 +1251,9 @@ 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, 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 c273f43ab..9fdf78de0 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1266,24 +1266,38 @@ 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; -// The twinned refinement is taken only where the law beats its null on the free set by this many standard errors. +// 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 - +// 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(); 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) @@ -1292,20 +1306,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; } } @@ -3028,11 +3049,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 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 -// 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) { @@ -3091,17 +3113,20 @@ 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 = 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; @@ -3112,6 +3137,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 1a85a59ba..faf34f585 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. @@ -2153,20 +2169,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 0d7d38991..94e3e27cd 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 @@ -4383,7 +4384,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) @@ -4395,7 +4399,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)); @@ -5380,6 +5385,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; // a copy: a rival taken below replaces *placed_model_out // In the group RefinableIn chose: the model's own where it has the rotations of the group written - // the screws are then the model's, as its fit already took them - so a model whose screws the merge diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index 1dd5b3990..799301cca 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -1540,6 +1540,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); @@ -1636,7 +1648,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