From 55dd0d4f6baeb1b9677d8ed4d1e004d32cc9bd41 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 10:47:55 +0200 Subject: [PATCH] rugnux: doubled-cell hypothesis - the supercell probe reads its class to the limit, and a class that stands clear is merged in the doubled cell as well Supercell probe, log-only part: the class it reads over 20-3 A is followed out to the run's own limit (a second prediction of the same supercell on the same 60 frames, so the 20-3 A band and the SUPERCELL_* keys are unchanged). Half-integer reflections are counted only where no lattice reflection predicted on the same frame lies within r3, and each resolution shell is read against the median of the six other parity classes - which a single doubling leaves empty - as an excess in units of the lattice's intensity and as a z in the class's own standard errors, the strongest sibling's z beside it ("Supercell probe to the limit" in the log). Doubled-cell hypothesis: where that class stands clear of the six empty ones (measured above its noise and at more than twice the strongest of them, in at least three shells), the run integrates and merges the doubled cell too, forcing the last pass's indexing result with its lattice doubled along the class (orientation, geometry and axis kept). One of the two is written as the answer, the other's merged .mtz/.cif/.hkl beside it as _doubled or _halfcell; the report names both (DOUBLED_CELL_*). Which one is the answer is a policy behind one constant: ADDED_CLASS_SIGNIFICANT (the doubled cell wherever the class stands clear) or HALF_CELL_BREAKS_DOWN (default: only where the class reaches a quarter of the lattice's intensity in the finest shell both are measured in). The doubled lattice is classed with a near-zero tolerance: the default lattice search promoted a triclinic doubled cell with an 88.2 deg angle to orthorhombic C, and no frame indexed. Measured on forced doubled-cell merges (rc175 b530c2d base): the open-arm crystal deposited in the doubled cell merges at R-free 0.275 against its deposited model, against 0.574 in the half cell (odd/even 0.02 at >8 A rising to 0.65 at 1.95 A, odd half-set CC1/2 0.93 at 8-5 A). Five crystals deposited in the half cell read classes the doubled merge holds as well (odd CC1/2 0.73-0.98 at 8-4 A): two rise with resolution like the first (to 0.4-0.5 at 2.5 A), one is strong from low resolution (0.3-0.5), two stay flat at 0.1-0.2. Nothing measured on intensities - the trend, the neighbour correlation, the separation from the nearest lattice reflection after the r3 filter - separated the doubled deposition from the rising half-cell ones, so the policy is the owner's call. Also: the SUPERCELL_POSSIBLE prose no longer says correct cells "read the same", and the warning no longer advises -C with the doubled cell, which the indexer reduces back to the half cell. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX_REPORT.md | 42 ++++++--- image_analysis/IndexAndRefine.cpp | 73 +++++++++++++- image_analysis/IndexAndRefine.h | 18 ++++ rugnux/ResultReport.cpp | 51 ++++++---- rugnux/Rugnux.cpp | 10 +- rugnux/Rugnux.h | 38 +++++++- rugnux/RugnuxImagePass.cpp | 61 ++++++++++-- rugnux/RugnuxPasses.cpp | 152 +++++++++++++++++++++++++++++- tests/IndexingUnitTest.cpp | 17 ++++ tests/RugnuxTest.cpp | 32 +++++++ 10 files changed, 441 insertions(+), 53 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index fd76f8b88..0695376f5 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -595,30 +595,42 @@ On rotation data, section 4 reports whether the lattice the run adopted has inte positions it does not index. On 60 frames spread over the sweep, after each frame's own integration, the lattice is predicted doubled along all three primitive axes and integrated to 3 Å; the reflections split into eight parity classes of h, k and l, of which `0 0 0` is the lattice itself and each of the other -seven is one index-2 superstructure. Nothing is decided on it: whether a superstructure belongs in the -cell is as much the depositor's call as the data's, and several crystals whose accepted cell is the -sub-cell carry one. +seven is one index-2 superstructure. The class read there is also followed out to the run's own limit, +counting only half-integer reflections that no lattice reflection on the same frame lies within the +background radius of, each resolution shell against the median of the six other classes - which a single +doubling leaves empty, so they measure what sits at a half-integer node without being a reflection. The +log prints that table (`Supercell probe to the limit`). **`SUPERCELL_CLASS=`** is the parity class with the most intensity over 20–3 Å. **`SUPERCELL_OCCUPANCY_PCT=`** is its mean intensity against the lattice's own reflections on the same frames, and **`SUPERCELL_ROCK_PCT=`** (± **`SUPERCELL_ROCK_SE_PCT=`**) the part of it that follows the partiality the way a Bragg reflection does, from a fit I = a + b p over the class, on the same scale. **`SUPERCELL_I_OVER_SIGMA=`** is its mean I/σ. A class near zero on both is empty. One that is occupied -and rocks is a superstructure whose reflections this run did not integrate; **`SUPERCELL_DOUBLED_CELL=`** -is the Niggli-reduced cell the lattice would double to, to give with `-C` to process on it. One that is -occupied but hardly rocks is diffuse or disordered intensity rather than Bragg reflections. +and rocks holds reflections of a superstructure this cell does not index; **`SUPERCELL_DOUBLED_CELL=`** +is the Niggli-reduced cell the lattice doubles to. One that is occupied but hardly rocks is diffuse or +disordered intensity rather than Bragg reflections. **`SUPERCELL_POSSIBLE=`** is `TRUE` where the class is measured (`SUPERCELL_I_OVER_SIGMA=` at least 0.5) and part of it rocks like Bragg reflections (`SUPERCELL_ROCK_PCT=` at least 2 %, three standard -errors clear of zero). It is advice to check, not a finding: the same numbers come from a real doubled -cell and from a correct cell with weak ordered intensity between its reflections - or with further -lattice domains whose spots land on the half-integer positions - and which of the two a structure is -decided by refinement. Process both settings - the run's cell, and `SUPERCELL_DOUBLED_CELL=` given with -`-C` - and compare them there. A `TRUE` is also a warning under the `SUPERCELL_POSSIBLE` flag, worded as a -prompt to check: on a battery of rotation data sets most of the crystals it named refine normally in the -sub-cell (rocking parts of 2-16 %), but a real doubled cell can read inside that range too. The summary's -`Supercell` line carries the same advice, and says where further lattice domains were found that may put -spots on the half-integer positions. +errors clear of zero). It says the doubled cell is there in the data, not that it is the cell to refine +in: a crystal whose two halves are nearly alike (a pseudo-translation) is described by either cell, and +spots of further lattice domains can land on the half-integer positions too. A `TRUE` is also a warning +under the `SUPERCELL_POSSIBLE` flag, and the summary's `Supercell` line carries the same advice. + +**Both cells merged.** Where the class stands clear of the six empty ones out to the limit (in at least +three shells measured above its noise, at more than twice the strongest of them in standard errors), the +run integrates and merges the doubled cell as well, on the same orientation and geometry, and writes one +of the two as its answer and the other's merged files beside it. Which one is the answer is a policy, +because nothing measured in the intensities separates a crystal deposited in the doubled cell from one +deposited in the half cell: the doubled cell is adopted where the class reaches a quarter of the lattice's +intensity in the finest shell both are measured in - where the half cell no longer approximates the +crystal - and a class that stays a small fraction of the lattice at every resolution keeps the half cell. +**`DOUBLED_CELL_ADOPTED=`** says which was written as the answer and **`DOUBLED_CELL_DECISION=`** why; +**`DOUBLED_CELL_OTHER_FILE=`** names the other hypothesis's merged file (`_doubled.mtz`, or +`_halfcell.mtz` where the doubled cell was adopted, with its `.cif` and `.hkl`), and +**`DOUBLED_CELL_OTHER_CELL=`**, **`DOUBLED_CELL_OTHER_SPACE_GROUP=`**, **`DOUBLED_CELL_OTHER_D_MIN=`**, +**`DOUBLED_CELL_OTHER_CC_HALF=`** and, with `--model`, **`DOUBLED_CELL_OTHER_R_FREE=`** describe it. +Compare the two in refinement. ## Further lattices diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 1e717dde5..fec1c4bb7 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include "IndexAndRefine.h" #include "bragg_integration/CalcISigma.h" @@ -618,6 +619,15 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, ProbeSupercellFrame(msg, prediction, probe_integrate, outcome, settings_prediction); } +CrystalLattice DoubledLattice(const CrystalLattice &primitive, const std::array &parity) { + const Coord v[3] = {primitive.Vec0(), primitive.Vec1(), primitive.Vec2()}; + const int i = parity[0] ? 0 : parity[1] ? 1 : 2; + Coord t[3]; + for (int j = 0; j < 3; j++) + t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast(parity[j]); + return {t[0], t[1], t[2]}; +} + SupercellProbeFit FitSupercellProbe(const SupercellProbeClass &c) { SupercellProbeFit f; if (c.n < 3) @@ -692,8 +702,69 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction c.sum_pi += r.partiality * r.I; c.sum_ii += static_cast(r.I) * r.I; } + + // The same supercell again, to the run's own limit (SupercellTrend). Its own prediction, so the band + // above integrates exactly what it always has. + BraggPredictionSettings s_trend = settings; + s_trend.centering = 'P'; + ApplyPredictionRange(s_trend, experiment, super); + const auto ntrend = prediction.Calc(outcome.experiment, super, s_trend); + const std::vector &predicted = prediction.GetReflections(); + + // The lattice's own reflections on this frame, binned on a grid of r3 squares so that every one + // within r3 of a point is in the 3x3 squares around it. + const float r3 = experiment.GetBraggIntegrationSettings().GetR3(); + const auto square = [r3](float x, float y) { + return (static_cast(std::floor(x / r3)) << 32) + | static_cast(static_cast(std::floor(y / r3))); + }; + std::unordered_map>> lattice_nodes; + for (int i = 0; i < ntrend; i++) + if (!((predicted[i].h | predicted[i].k | predicted[i].l) & 1)) + lattice_nodes[square(predicted[i].predicted_x, predicted[i].predicted_y)].emplace_back( + predicted[i].predicted_x, predicted[i].predicted_y); + const auto clear_of_lattice = [&](float x, float y) { + for (int dx = -1; dx <= 1; dx++) + for (int dy = -1; dy <= 1; dy++) { + const auto it = lattice_nodes.find(square(x + dx * r3, y + dy * r3)); + if (it == lattice_nodes.end()) + continue; + for (const auto &[nx, ny] : it->second) + if ((nx - x) * (nx - x) + (ny - y) * (ny - y) < r3 * r3) + return false; + } + return true; + }; + + SupercellTrend trend{}; + for (const auto &r : integrate(predicted, ntrend, msg.number)) { + if (!(r.d <= SUPERCELL_TREND_D_A[0]) || !(r.sigma > 0.0f) || !std::isfinite(r.I)) + continue; + const int c = 4 * (r.h & 1) + 2 * (r.k & 1) + (r.l & 1); + if (c != 0 && !clear_of_lattice(r.predicted_x, r.predicted_y)) + continue; + int sh = 0; + while (r.d <= SUPERCELL_TREND_D_A[sh + 1]) + sh++; + auto &t = trend[c][sh]; + t.n++; + t.sum_i += r.I; + t.sum_i_over_sigma += r.I / r.sigma; + t.sum_ii += static_cast(r.I) * r.I; + } + const std::unique_lock ul(supercell_probe_mutex_); - supercell_probe_by_frame_.insert_or_assign(msg.number, ProbedFrame{frame, prim}); + supercell_probe_by_frame_.insert_or_assign(msg.number, ProbedFrame{frame, trend, prim}); +} + +SupercellTrend IndexAndRefine::GetSupercellTrend() { + const std::unique_lock ul(supercell_probe_mutex_); + SupercellTrend sum{}; + for (const auto &[image, f] : supercell_probe_by_frame_) + for (int i = 0; i < 8; i++) + for (int j = 0; j < SUPERCELL_TREND_SHELLS; j++) + sum[i][j] += f.trend[i][j]; + return sum; } std::optional diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index a5eb15a28..87ac0300a 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -50,6 +50,22 @@ struct SupercellProbeClass { // [class][shell], the shells 20-5 A and 5-3 A. using SupercellProbe = std::array, 8>; +// The supercell probe's reach to the run's resolution limit: the eight parity classes again, in fixed +// resolution shells from 20 A down to the prediction limit (SUPERCELL_TREND_D_A gives the shell edges), +// with a half-integer reflection counted only where no lattice reflection predicted on the same frame +// lies within the background radius r3 of it - so what is summed is its own intensity, not the tail of +// a neighbour's. A superstructure from a pseudo-translation is weakest at low resolution, where the two +// halves of the cell look alike, and grows towards the limit, which the 20-3 A band cannot show. +inline constexpr float SUPERCELL_TREND_D_A[] = {20.0f, 8.0f, 5.0f, 4.0f, 3.3f, 2.8f, 2.4f, 2.1f, 1.9f, 1.7f, + 1.5f, 1.3f, 1.1f, 0.0f}; +inline constexpr int SUPERCELL_TREND_SHELLS = std::size(SUPERCELL_TREND_D_A) - 1; +using SupercellTrend = std::array, 8>; + +// The index-2 supercell of a primitive lattice whose added reflections are the half-integer class of +// parity (h,k,l): the vectors t with (h,k,l).t even - twice one axis the class is odd on, and each +// other axis plus that one where the class is odd on it too. Not reduced. +CrystalLattice DoubledLattice(const CrystalLattice &primitive, const std::array &parity); + // One parity class of the supercell probe, fitted I = a + b p over partiality p: b (with its standard // error) is what rocks like a Bragg reflection, a what does not. struct SupercellProbeFit { @@ -126,6 +142,7 @@ class IndexAndRefine { std::vector supercell_probe_frames_; struct ProbedFrame { SupercellProbe probe; + SupercellTrend trend; CrystalLattice primitive; }; std::map supercell_probe_by_frame_; @@ -195,6 +212,7 @@ public: // The probe's integrations are the caller's to keep out of anything else it counts. void ProbeSupercell(std::vector image_numbers); [[nodiscard]] SupercellProbe GetSupercellProbe(); + [[nodiscard]] SupercellTrend GetSupercellTrend(); // The primitive lattice the parity classes are indexed in (the highest probed frame's; only its metric // is meant to be read), so a class can be turned into the cell it would double to. [[nodiscard]] std::optional GetSupercellProbePrimitive(); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index d8f50ed55..55ab95b5a 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1426,10 +1426,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } // ---- index-2 superstructure - // Report only, like the harmonic above: on 60 frames the probe integrates the lattice doubled - // along all three primitive axes and reads the seven half-integer parity classes. Whether a - // class that is there belongs in the cell is the depositor's call as much as the data's - - // several crystals whose accepted cell is the sub-cell carry one - so no verdict is printed. + // On 60 frames the probe integrates the lattice doubled along all three primitive axes and reads + // the seven half-integer parity classes. Whether a class that is there belongs in the cell is the + // depositor's call as much as the data's - several crystals whose accepted cell is the sub-cell + // carry one - so where the run merged both cells (DOUBLED_CELL_*) the report names both. if (result.supercell) { const auto &sc = *result.supercell; Add(s, KeyText("SUPERCELL_CLASS", fmt::format("{} {} {}", sc.h, sc.k, sc.l))); @@ -1445,25 +1445,40 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Warn(doc, PathologyCode::SUPERCELL_POSSIBLE, fmt::format("Check the cell: half-integer reflections of parity {} {} {} carry Bragg-like intensity " "({:.0f}% of the lattice's, rocking part {:.1f}%), so the true cell may be twice as " - "large - or this is weak ordered intensity of a correct cell{}. Process both settings - " - "this cell and -C \"{:.2f},{:.2f},{:.2f},{:.2f},{:.2f},{:.2f}\" - and compare them in refinement", + "large - or the crystal's two halves are nearly alike and either cell describes it{}. {}", sc.h, sc.k, sc.l, sc.occupancy_pct, sc.rock_pct, - domains ? ", or spots of the further lattice domains found landing on those positions" : "", - sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c, - sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma)); + domains ? ", or spots of the further lattice domains found land on those positions" : "", + result.doubled_cell + ? fmt::format("Both cells were merged; the other one is in {}.mtz - compare them in refinement", + result.doubled_cell->other_prefix) + : fmt::format("The doubled cell {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} was not merged", + sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c, + sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma))); + if (result.doubled_cell) { + const auto &h = *result.doubled_cell; + Add(s, KeyBool("DOUBLED_CELL_ADOPTED", h.adopted)); + Add(s, KeyText("DOUBLED_CELL_DECISION", h.decision)); + Add(s, KeyText("DOUBLED_CELL_OTHER_FILE", h.other_prefix + ".mtz")); + Add(s, KeyText("DOUBLED_CELL_OTHER_CELL", CellString(h.other_cell))); + Add(s, KeyText("DOUBLED_CELL_OTHER_SPACE_GROUP", h.other_space_group)); + Add(s, KeyReal("DOUBLED_CELL_OTHER_D_MIN", h.other_d_min_A, "{:.2f}")); + Add(s, KeyReal("DOUBLED_CELL_OTHER_CC_HALF", h.other_cc_half, "{:.4f}")); + if (std::isfinite(h.other_r_free)) + Add(s, KeyReal("DOUBLED_CELL_OTHER_R_FREE", h.other_r_free, "{:.4f}")); + } Add(s, Blank()); Add(s, Prose( " SUPERCELL_CLASS is the half-integer class (parity of h, k, l in the primitive cell of the\n" " lattice) with the most intensity over 20-3 A, read on 60 frames. OCCUPANCY is its mean\n" " intensity against the lattice's own reflections; ROCK is the part of it that follows the\n" " partiality the way a Bragg reflection does, on the same scale. A class near zero on both is\n" - " empty. One that is occupied and rocks is an index-2 superstructure the lattice does not\n" - " index: its reflections are not integrated. To process on the doubled lattice, give\n" - " SUPERCELL_DOUBLED_CELL with -C. An occupied class that does not rock is diffuse or\n" - " disordered intensity rather than Bragg reflections. SUPERCELL_POSSIBLE is TRUE, and a\n" - " warning to check, where the class is measured and rocks; correct cells with weak ordered\n" - " intensity, or with further lattice domains whose spots land on the half-integer\n" - " positions, read the same, so only refinement of both cells decides.")); + " empty. One that is occupied and rocks holds reflections of an index-2 superstructure the\n" + " lattice does not index, and they are not integrated in this cell. An occupied class that\n" + " does not rock is diffuse or disordered intensity rather than Bragg reflections.\n" + " SUPERCELL_POSSIBLE is TRUE, and a warning to check, where the class is measured and rocks.\n" + " It says the doubled cell is there in the data, not that it is the cell to refine in: a\n" + " crystal whose two halves are nearly alike (a pseudo-translation) is described in either,\n" + " and spots of further lattice domains can also land on the half-integer positions.")); Add(s, Blank()); } @@ -2783,8 +2798,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, checked_row("Supercell", {PathologyCode::SUPERCELL_POSSIBLE}, "no indication"); else checked_row("Supercell", {PathologyCode::SUPERCELL_POSSIBLE}, fmt::format("POSSIBLE, check - class {} {} {} rocks like Bragg reflections " - "(rocking part {:.1f}%); correct cells with weak ordered intensity " - "read the same{}; {}", sc.h, sc.k, sc.l, sc.rock_pct, + "(rocking part {:.1f}%); a nearly exact pseudo-translation is described " + "in either cell{}; {}", sc.h, sc.k, sc.l, sc.rock_pct, domains ? ", and the further lattice domains found may put " "spots there" : "", doubled)); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 659039ccc..20793707f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -861,12 +861,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // at the half cell by every arbiter. The probe reads the question in intensities instead, on the // same frames the first pass validates on (see IndexAndRefine::ProbeSupercell). // - // REPORT ONLY - it decides nothing. Measured on the battery, a weak real class is not yet told apart - // from what sits at the half-integer nodes of crystals whose cell is right: a class at a fifth of - // the lattice's intensity on a crystal whose axis is truly doubled, against up to a tenth on - // correct cells (diffuse scattering and the tails of long-axis neighbours, which do not rock like - // Bragg reflections - but neither, on that crystal, does its real class). The line it logs is the - // population the decision has to be calibrated on. + // It decides nothing in this pass. Read on its own it does not tell a weak real class from what sits + // at the half-integer nodes of crystals whose cell is right, so the run carries the doubled cell as a + // second hypothesis where the class stands clear out to the limit (DoubledCellHypothesis), and merges + // both. if (full && indexer && experiment_.IsRotationIndexing()) { std::vector frames; for (const int o : select_equally_spaced_image_ordinals(images_to_process, 60)) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0f2806303..e4201767b 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -439,11 +439,20 @@ struct ProcessResult { // The strongest index-2 superstructure class the supercell probe read (20-3 A, 60 frames): which // parity class of the primitive cell, its mean intensity and the part that rocks like a Bragg // reflection, both against the lattice's own, and the cell the lattice would double to were it - // real. Report only - nothing decides on it. + // real. What reads it is the doubled-cell hypothesis (Rugnux::DoubledCellHypothesis). struct SupercellClass { int h = 0, k = 0, l = 0; // parity of the primitive indices double occupancy_pct = 0.0, rock_pct = 0.0, rock_se_pct = 0.0, i_over_sigma = 0.0; UnitCell doubled_cell{}; + // The class out to the run's limit (SupercellTrend), by resolution shell, less the median of the + // six other classes: its mean intensity against the lattice's (excess) and in its own standard + // errors (z), the largest such z among the six beside it, and the lattice's own . + struct Shell { + float d_low_A = 0.0f, d_high_A = 0.0f; + int64_t n = 0; + double excess = 0.0, z = 0.0, sibling_z = 0.0, lattice_i_over_sigma = 0.0; + }; + std::vector trend; }; std::optional supercell; // Rotation: what the spots the crystal's lattice leaves over index to (split-crystal domains, twin @@ -508,6 +517,20 @@ struct ProcessResult { }; std::optional twin_subgroup_merge; + // The doubled-cell hypothesis (Rugnux::DoubledCellHypothesis): where the supercell probe's class + // stands clear of the crystal's own null out to the limit, the run integrates and merges in the + // doubled cell as well, and one of the two cells is the answer. What the other one is, and where its + // merged files are. Unset where the hypothesis was not run. + struct DoubledCell { + bool adopted = false; // the answer is the doubled cell, the other one the probed cell + std::string decision; // why, in the words the run reports it in + std::string other_prefix; // _halfcell or _doubled + UnitCell other_cell{}; + std::string other_space_group; + double other_d_min_A = NAN, other_cc_half = NAN, other_r_free = NAN; + }; + std::optional doubled_cell; + // What --model found, when one was given. Set even when the validation could not run - its // failure_reason then says why - so the report can report the failure rather than leave a run // that was asked to validate looking like one that was not. Unset when no model was asked for. @@ -625,6 +648,13 @@ constexpr float OVERLOAD_LOW_RES_A = 4.0f; // reports it in, or an empty string when the refined pass stands. std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &pass2); +// The doubled-cell hypothesis (Rugnux::DoubledCellHypothesis), read off the supercell probe's class out +// to the limit (SupercellClass::trend). Worth merging: the class stands clear of the crystal's own null. +// Adopt: the doubled cell is the answer, under DOUBLED_CELL_POLICY; empty when the probed cell stands, +// otherwise why the doubled one is adopted. +bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc); +std::string AdoptDoubledCell(const ProcessResult::SupercellClass &sc); + // Whether two passes hold the same lattice: the same Bravais class and primitive volumes within 2 %. bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b); @@ -737,6 +767,8 @@ class Rugnux { // Set by Run()'s supercell-collapse guard to force pass-1's full result into the second-pass indexer // (RunPipeline calls indexer->ForceRotationIndexerResult when this is set); reset right after that re-run. std::optional force_rotation_result_; + // The indexing result the last full pass integrated on, which the doubled-cell hypothesis doubles. + std::optional last_rotation_result_; // The space group the USER fixed (-S on the CLI, the settings panel in the viewer), captured before // anything can clear it - a pass that determines a group of its own writes it onto the experiment, @@ -1032,6 +1064,10 @@ class Rugnux { }; // The refined passes of the rotation two-pass, from the second pass on. Empty where they ran beside // the rotation-scale walk and it overturned them. + // After the refined passes: where the supercell probe's class stands clear, integrate and merge the + // doubled cell too, and return whichever of the two is the answer (ProcessResult::doubled_cell). + ProcessResult DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half, + const std::string &base_prefix); std::optional RefinedPasses(RugnuxObserver *observer, ProcessResult pass1, const RefinedPassStart &start); // The rotation-scale walk, run on the run while the refined passes run on a copy of it; `check` is diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index f0d159506..63c185b13 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -417,9 +417,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { "mean intensity, and the fit I = a + b p against the lattice's own):{}", table); // The class with the most intensity over 20-3 A goes into the report, with the cell it would - // double the lattice to: the vectors t of the primitive lattice with (h,k,l).t even - twice - // one axis the class is odd on, and each other axis plus that one where the class is odd - // on it too - Niggli-reduced. + // double the lattice to (DoubledLattice), Niggli-reduced. const auto band = [&](int c) { SupercellProbeClass x = probe[c][0]; x += probe[c][1]; return x; }; const SupercellProbeFit f0 = FitSupercellProbe(band(0)); const auto prim = indexer->GetSupercellProbePrimitive(); @@ -429,19 +427,61 @@ void Rugnux::ImagePass(PipelineLocals &p) { best = c; if (prim && f0.mean_i > 0.0 && f0.b != 0.0) { const SupercellProbeFit f = FitSupercellProbe(band(best)); - const int p[3] = {(best >> 2) & 1, (best >> 1) & 1, best & 1}; - const Coord v[3] = {prim->Vec0(), prim->Vec1(), prim->Vec2()}; - const int i = p[0] ? 0 : p[1] ? 1 : 2; - Coord t[3]; - for (int j = 0; j < 3; j++) - t[j] = j == i ? v[i] * 2.0f : v[j] + v[i] * static_cast(p[j]); + const std::array p = {(best >> 2) & 1, (best >> 1) & 1, best & 1}; ProcessResult::SupercellClass sc; sc.h = p[0]; sc.k = p[1]; sc.l = p[2]; sc.occupancy_pct = 100.0 * f.mean_i / f0.mean_i; sc.rock_pct = 100.0 * f.b / f0.b; sc.rock_se_pct = 100.0 * f.b_se / f0.b; sc.i_over_sigma = f.mean_i_over_sigma; - sc.doubled_cell = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell(); + sc.doubled_cell = DoubledLattice(*prim, p).NiggliReduce().GetUnitCell(); + + // The same class out to the limit, half-integer reflections clear of the lattice's only + // (SupercellTrend). The six other classes, which a single doubling leaves empty, are the + // crystal's own null: measured on the same frames through the same integration, they + // carry whatever sits at a half-integer node without being a reflection, so the class + // is read against their median - in units of the lattice's intensity (excess) and of + // its own standard error (z), the strongest sibling's z beside it. + const SupercellTrend trend = indexer->GetSupercellTrend(); + const auto mean = [](const SupercellProbeClass &c) { return c.sum_i / static_cast(c.n); }; + const auto se = [&](const SupercellProbeClass &c) { + const double m = mean(c); + return std::sqrt(std::max(0.0, c.sum_ii / static_cast(c.n) - m * m) + / static_cast(c.n)); + }; + std::string trend_table; + for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) { + const SupercellProbeClass &l = trend[0][sh], &x = trend[best][sh]; + std::vector siblings; + for (int c = 1; c < 8; c++) + if (c != best && trend[c][sh].n > 1) + siblings.push_back(c); + if (l.n == 0 || x.n < 2 || siblings.size() < 6 || l.sum_i <= 0.0) + continue; + std::vector sib_mean; + for (const int c : siblings) + sib_mean.push_back(mean(trend[c][sh])); + std::ranges::sort(sib_mean); + const double null_i = (sib_mean[2] + sib_mean[3]) / 2.0; + ProcessResult::SupercellClass::Shell shell; + shell.d_low_A = SUPERCELL_TREND_D_A[sh]; + shell.d_high_A = SUPERCELL_TREND_D_A[sh + 1]; + shell.n = x.n; + shell.lattice_i_over_sigma = l.sum_i_over_sigma / static_cast(l.n); + shell.excess = (mean(x) - null_i) / mean(l); + shell.z = (mean(x) - null_i) / se(x); + shell.sibling_z = -INFINITY; + for (const int c : siblings) + shell.sibling_z = std::max(shell.sibling_z, (mean(trend[c][sh]) - null_i) / se(trend[c][sh])); + sc.trend.push_back(shell); + trend_table += fmt::format("\n {:.1f}-{:.1f} A: lattice n {} {:.2f}; class n {} excess {:.3f} " + "z {:.1f}, other classes z up to {:.1f}", shell.d_low_A, shell.d_high_A, + l.n, shell.lattice_i_over_sigma, x.n, shell.excess, shell.z, + shell.sibling_z); + } + logger.Info("Supercell probe to the limit, class {}{}{} (its reflections clear of the lattice's by r3, " + "less the median of the six other classes: mean intensity against the lattice's own, " + "and in standard errors):{}", p[0], p[1], p[2], trend_table); result.supercell = sc; } } @@ -524,6 +564,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { if (full && !cancelled_) { if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) { + last_rotation_result_ = *rot; end_msg.rotation_lattice = rot->lattice; // The refined geometry, not the nominal StartMessage values. Two things about these five // numbers that are not obvious from their names. The beam centre is the PONI - the foot of diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 1ad46836d..f3d7b8091 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -592,7 +593,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { if (!check->overturns) { if (output_error) std::rethrow_exception(output_error); - return *result; + return output->DoubledCellHypothesis(observer, std::move(*result), base_prefix); } } @@ -613,7 +614,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { logger.Info("Two-pass: the rotation-scale walk changed the spot finding the refined passes " "start from, so they run again after it"); } - return *RefinedPasses(observer, std::move(pass1), start); + return DoubledCellHypothesis(observer, *RefinedPasses(observer, std::move(pass1), start), base_prefix); } return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); } @@ -1407,3 +1408,150 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro pass2.geometry_not_converged = geometry_not_converged; return pass2; } + +// The supercell probe reads one half-integer class out to the limit (SupercellClass::trend), each +// shell less the median of the six other classes - which a single doubling leaves empty, so they are the +// null measured on the same frames, through the same integration and the same neighbour filter. +// +// Measured on the forced doubled-cell merges of six crystals whose class the probe reads: a class read +// this way is a class the doubled merge holds (odd-class half-set CC1/2 0.7-0.98 at 8-4 A on all six), +// but nothing read from intensities told the crystal whose deposited cell is the doubled one from those +// deposited in the half cell - the odd/even ratio rose with resolution alike on both kinds, and the +// odd class correlated with its lattice neighbours alike. Which of the two cells to write is therefore +// a policy (DOUBLED_CELL_POLICY), and both are written. +namespace { + // A shell the class is read in: measured above its noise, against a lattice that is measured too. + bool DoubledCellShellMeasured(const ProcessResult::SupercellClass::Shell &sh) { + return sh.z >= 3.0 && sh.lattice_i_over_sigma >= 1.0; + } + + // ADDED_CLASS_SIGNIFICANT: the doubled cell wherever the class stands clear of the null - it + // describes every Bragg reflection the crystal has, where the half cell leaves a measured class + // out. + // HALF_CELL_BREAKS_DOWN: the doubled cell only where the class reaches DOUBLED_CELL_BREAKDOWN_EXCESS + // of the lattice's intensity in the finest shell both are measured in - where the half cell no + // longer approximates the crystal. A class that stays a small fraction of the lattice at every + // resolution (two halves that differ slightly throughout) keeps the half cell, as the depositions + // of such crystals do. + enum class DoubledCellPolicy { AddedClassSignificant, HalfCellBreaksDown }; + constexpr DoubledCellPolicy DOUBLED_CELL_POLICY = DoubledCellPolicy::HalfCellBreaksDown; + constexpr double DOUBLED_CELL_BREAKDOWN_EXCESS = 0.25; +} + +bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) { + // Clear of the null: in at least three shells the class is measured and stands at more than twice + // the strongest of the six empty classes, in standard errors. + int clear = 0; + for (const auto &sh : sc.trend) + clear += DoubledCellShellMeasured(sh) && sh.z > 2.0 * sh.sibling_z; + return clear >= 3; +} + +std::string AdoptDoubledCell(const ProcessResult::SupercellClass &sc) { + if (DOUBLED_CELL_POLICY == DoubledCellPolicy::AddedClassSignificant) + return DoubledCellWorthMerging(sc) ? "the added class stands clear of the other classes" : std::string(); + const ProcessResult::SupercellClass::Shell *outer = nullptr; + for (const auto &sh : sc.trend) + if (DoubledCellShellMeasured(sh)) + outer = &sh; + if (!outer || outer->excess < DOUBLED_CELL_BREAKDOWN_EXCESS) + return {}; + return fmt::format("the added class reaches {:.2f} of the lattice's intensity at {:.1f}-{:.1f} A, the finest " + "shell both are measured in, so the half cell no longer approximates the crystal", + outer->excess, outer->d_low_A, outer->d_high_A); +} + +namespace { + std::string CellText(const UnitCell &c) { + return fmt::format("{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}", c.a, c.b, c.c, c.alpha, c.beta, c.gamma); + } +} + +ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half, + const std::string &base_prefix) { + if (cancelled_ || user_cell_ || user_fixed_sg_ || !half.supercell || !last_rotation_result_ + || !half.has_merge_statistics || !half.consensus_cell || !DoubledCellWorthMerging(*half.supercell)) + return half; + Logger logger("Rugnux"); + const auto &sc = *half.supercell; + const std::string adopt_why = AdoptDoubledCell(sc); + const bool adopt = !adopt_why.empty(); + + // The lattice the probed cell was integrated on, doubled along the class, with everything else of + // that indexing - orientation, geometry, axis - kept, so the two merges differ by the cell alone. + RotationIndexerResult doubled = *last_rotation_result_; + const CrystalLattice prim = doubled.lattice.ToPrimitive(doubled.search_result.centering); + // Classed only on what the refined metric already holds exactly: a looser search promotes a doubled + // cell whose angle merely lies near 90 deg (measured: 88.2 deg, held to 90, and no frame indexed). + doubled.search_result = LatticeSearch(DoubledLattice(prim, {sc.h, sc.k, sc.l}), 0.001, 0.05); + doubled.lattice = doubled.search_result.conventional; + doubled.extra_lattices.clear(); + doubled.unconstrained.reset(); + logger.Info("Doubled-cell hypothesis: lattice {} ({}-centred {}) doubled to {} ({}-centred {})", + CellText(last_rotation_result_->lattice.GetUnitCell()), last_rotation_result_->search_result.centering, + gemmi::crystal_system_str(last_rotation_result_->search_result.system), + CellText(doubled.lattice.GetUnitCell()), doubled.search_result.centering, + gemmi::crystal_system_str(doubled.search_result.system)); + + // The answer's files go under the plain prefix and the other cell's merged files beside them, so + // where the doubled cell is adopted the probed cell's merged files move aside before it writes. + const std::string half_prefix = base_prefix + "_halfcell"; + const std::string doubled_prefix = base_prefix + "_doubled"; + const auto move_merged = [&](const std::string &from, const std::string &to) { + for (const char *ext : {".mtz", ".cif", ".hkl"}) + if (std::error_code ec; std::filesystem::exists(from + ext, ec)) + std::filesystem::rename(from + ext, to + ext, ec); + }; + logger.Info("Doubled-cell hypothesis: the half-integer class {} {} {} stands clear of the other classes " + "out to the limit, so the run integrates and merges the doubled cell {:.2f} {:.2f} {:.2f} " + "{:.2f} {:.2f} {:.2f} as well; {}", sc.h, sc.k, sc.l, sc.doubled_cell.a, sc.doubled_cell.b, + sc.doubled_cell.c, sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma, + adopt ? "it is the answer - " + adopt_why : std::string("the probed cell stays the answer")); + if (observer) + observer->OnStep("Doubled-cell hypothesis"); + if (adopt && !base_prefix.empty()) + move_merged(base_prefix, half_prefix); + config_.output_prefix = base_prefix.empty() ? base_prefix : adopt ? base_prefix : doubled_prefix; + force_rotation_result_ = doubled; + // The quality guard compares a refined pass against the header pass in the same cell; this pass is + // in another cell and is judged here instead. + quality_guard_pass1_.reset(); + ProcessResult d; + try { + d = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Warning("Doubled-cell hypothesis: the pass in the doubled cell did not complete ({})", e.what()); + } + force_rotation_result_.reset(); + config_.output_prefix = base_prefix; + if (!d.has_merge_statistics || !d.consensus_cell) { + logger.Warning("Doubled-cell hypothesis: the doubled cell merged nothing - the probed cell stays the answer"); + if (adopt && !base_prefix.empty()) + move_merged(half_prefix, base_prefix); + return half; + } + d.pass_number = d.pass_count = half.pass_count + 1; + half.pass_count = d.pass_count; + + ProcessResult &answer = adopt ? d : half; + const ProcessResult &other = adopt ? half : d; + ProcessResult::DoubledCell h; + h.adopted = adopt; + h.decision = adopt ? adopt_why : DOUBLED_CELL_POLICY == DoubledCellPolicy::HalfCellBreaksDown + ? fmt::format("the probed cell stands: the added class stays below {:.2f} of the lattice's intensity " + "where both are measured", DOUBLED_CELL_BREAKDOWN_EXCESS) + : std::string("the probed cell stands"); + h.other_prefix = adopt ? half_prefix : doubled_prefix; + h.other_cell = *other.consensus_cell; + h.other_space_group = other.space_group ? other.space_group->xhm() : std::string(); + h.other_d_min_A = other.merge_statistics.overall.d_min; + h.other_cc_half = other.merge_statistics.overall.cc_half; + if (other.model_validation) + h.other_r_free = other.model_validation->r_free; + answer.doubled_cell = h; + // The answer's report describes the probe of the probed cell, which is what the decision read. + answer.supercell = half.supercell; + answer.pass_decision = fmt::format("doubled-cell hypothesis merged ({}); {}", h.decision, half.pass_decision); + return answer; +} diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index 2bd7e1be3..535d0dcbd 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -817,6 +817,23 @@ TEST_CASE("FitSupercellProbe", "[SupercellProbe]") { CHECK(FitSupercellProbe(SupercellProbeClass{}).b == 0.0); } +TEST_CASE("DoubledLattice", "[SupercellProbe]") { + const CrystalLattice prim(Coord(40, 0, 0), Coord(0, 50, 0), Coord(0, 0, 60)); + // Class 0 0 1: the c axis doubled, the others kept. + const CrystalLattice c2 = DoubledLattice(prim, {0, 0, 1}); + CHECK(c2.Vec0().x == Catch::Approx(40.0)); + CHECK(c2.Vec1().y == Catch::Approx(50.0)); + CHECK(c2.Vec2().z == Catch::Approx(120.0)); + // Class 1 1 0: a centred doubling of the ab plane - twice a, and b plus a - so every vector of the + // supercell has h + k even, at twice the volume. + const CrystalLattice ab = DoubledLattice(prim, {1, 1, 0}); + CHECK(std::fabs(ab.CalcVolume()) == Catch::Approx(2.0 * 40 * 50 * 60)); + CHECK(ab.Vec0().x == Catch::Approx(80.0)); + CHECK(ab.Vec1().x == Catch::Approx(40.0)); + CHECK(ab.Vec1().y == Catch::Approx(50.0)); + CHECK(ab.Vec2().z == Catch::Approx(60.0)); +} + TEST_CASE("PostIndexingRefinement_CoplanarCandidateIsRejected", "[Indexing]") { // A candidate whose three rows lie in one plane passes every length and angle test - three rows // 120 deg apart in a plane meet the 30-150 deg window - and every spot along the plane normal diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 3a99f3c2d..60f00785b 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -349,6 +349,38 @@ TEST_CASE("SpotWidth_Gaussian", "[process]") { // The rotation two-pass quality guard. The refined pass is sent back to the header geometry only for // a decisive loss of signal - reflections at I/sigma >= 2 net of the noise tail at <= -2 - over no // more of reciprocal space. +TEST_CASE("DoubledCellDecision", "[process]") { + // Shells of the supercell probe's trend: d range, excess over the null, z, the strongest null z and + // the lattice's . + const auto trend = [](std::vector> excess_z_lattice) { + ProcessResult::SupercellClass sc; + float d = 20.0f; + for (const auto &[excess, z, lattice] : excess_z_lattice) { + ProcessResult::SupercellClass::Shell sh; + sh.d_low_A = d; + sh.d_high_A = d *= 0.7f; + sh.excess = excess; + sh.z = z; + sh.sibling_z = 3.0; + sh.lattice_i_over_sigma = lattice; + sc.trend.push_back(sh); + } + return sc; + }; + // A pseudo-translation that breaks down: weak at low resolution, a third of the lattice where the + // lattice is still measured. + const auto breaks_down = trend({{0.02, 17, 18}, {0.10, 38, 7}, {0.26, 36, 6}, {0.35, 45, 2.4}, {0.46, 33, 0.8}}); + CHECK(DoubledCellWorthMerging(breaks_down)); + CHECK(!AdoptDoubledCell(breaks_down).empty()); + // A class that stays a tenth of the lattice at every resolution: merged as well, the half cell kept. + const auto flat = trend({{0.09, 10, 70}, {0.11, 23, 28}, {0.10, 24, 33}, {0.11, 24, 18}, {0.10, 23, 8}}); + CHECK(DoubledCellWorthMerging(flat)); + CHECK(AdoptDoubledCell(flat).empty()); + // A class no stronger than the null: nothing. + const auto null_class = trend({{0.01, 4, 20}, {0.02, 5, 10}, {0.03, 2, 5}, {0.05, 1, 2}}); + CHECK(!DoubledCellWorthMerging(null_class)); +} + TEST_CASE("RefinedPassIsWorse", "[process]") { const auto pass = [](const UnitCell &cell, char centering, int64_t strong, int64_t negative, int64_t reflections) {