From 55dd0d4f6baeb1b9677d8ed4d1e004d32cc9bd41 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 10:47:55 +0200 Subject: [PATCH 1/8] 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) { From 48201d35ec1ae51e9f7b9cf0d9bdf7629305e4e9 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 13:05:11 +0200 Subject: [PATCH 2/8] WIP: rugnux doubled cell decided late on one pseudo-translation; screw trim reads the measured intensity Doubled-cell hypothesis: the run now always merges the doubled cell under _doubled first and decides with both merges in hand. Adopted only where (1) the probe's class grows significantly with resolution (three standard errors between the 20-8 A shell and the finest measured one) and (2) in the doubled merge one pseudo-translation - half a cell plus eps, fitted over the class - modulates at least half of the added class in every shell (FitDoubledClassTranslation: A sin^2(pi h.eps) + B per shell, share 1 - B/mean). Where adopted, every _doubled file takes the plain name and the probed cell's file of that name moves to _halfcell. Replaces the HALF_CELL_BREAKS_DOWN bound. Measured on the forced doubled merges: modulated share per shell 0.80-1.0 out to 1.9 A on the crystal deposited in the doubled cell; on the ones deposited in the half cell it falls to 0.01-0.34 by 2.1-2.8 A (the rising class with twinning), or stays at -0.3-0.48 at every resolution (three others); the one class it fully explains does not rise (eps 0.1 of a cell, flat at 0.1 of the lattice). Space-group search: the trim of a screw zone's largest absence exempted a violation "at its row's mean" on the de-modulated intensity, so under a deep pseudo-translation a reflection at 0.015 of its row (de-modulated by 1/0.002) could never be trimmed and alone decided the row against the screw. The exemption now reads the measured intensity, at both sites. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- .../scale_merge/SearchSpaceGroup.cpp | 13 +- rugnux/Rugnux.h | 33 ++- rugnux/RugnuxPasses.cpp | 277 ++++++++++++++---- rugnux/RugnuxScaleMerge.cpp | 2 + tests/RugnuxTest.cpp | 48 ++- 5 files changed, 313 insertions(+), 60 deletions(-) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index da2765a2a..7eb93b418 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -1676,7 +1676,8 @@ SearchSpaceGroupResult SearchSpaceGroup( const double u = std::max(0.0, e) / mean; violations += violation; sum_u += u; - if (!(violation && u >= 1.0)) + // At its row's mean as MEASURED, not as de-modulated (see the screw zone below). + if (!(violation && Esq[i] >= mean)) max_u = std::max(max_u, u); } const int a = static_cast(odd.size()); @@ -2500,7 +2501,15 @@ SearchSpaceGroupResult SearchSpaceGroup( // eight 0k0 are genuinely dead carries one violation at 0.31 of its row - the // mis-measurement the trim exists for, and one that must stay trimmable or a real // 2_1 is lost. Everything below the row mean is trimmed exactly as before. - if (!(violation && u >= 1.0)) + // + // "At its row's mean" is read on the MEASURED intensity, not the de-modulated one. + // Under a pseudo-translation u is the reflection divided by the depth the modulation + // puts on its class, so a reflection at a hundredth of its row reads at its mean once + // the class is held at a two-hundredth: measured on a doubled orthorhombic cell, one + // 00l at 0.015 of its row (5.9 sigma, beside reflections at 25) became the one member + // no trim could touch, and on its own it decided the row against the screw that the + // other 26 dead reflections - and the deposited structure factors - carry. + if (!(violation && a.e_squared >= mit->second)) z.max_u = std::max(z.max_u, u); z.n_control = static_cast(row_present_esq.at(a.row).size()); } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index e4201767b..21227a34b 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -530,6 +530,9 @@ struct ProcessResult { double other_d_min_A = NAN, other_cc_half = NAN, other_r_free = NAN; }; std::optional doubled_cell; + // The merged reflections as written, kept only where the run asks for them (the doubled-cell pass, + // whose merge the hypothesis reads). Empty otherwise. + std::vector merged_reflections; // 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 @@ -650,10 +653,32 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult & // 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. +// Rises: the class grows significantly with resolution, the shape of a pseudo-translation breaking down; +// empty where it does not, otherwise the reason in the words the run reports it in. bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc); -std::string AdoptDoubledCell(const ProcessResult::SupercellClass &sc); +std::string DoubledCellRises(const ProcessResult::SupercellClass &sc); + +// How much of a doubled-cell merge's added class ONE pseudo-translation explains. Two copies related by +// half a doubled-cell vector plus eps put |F|^2 sin^2(pi h.eps) on the added reflections and cos^2 on +// the others, so the added class, in units of the others' shell mean, is fitted as A sin^2(pi h.eps) + B: +// eps over the whole class, A and B shell by shell. B is what is there whatever the direction - copies +// that differ in content, or intensity that is not the crystal's (twin domains, neighbours) - and +// `modulated_share` is 1 - B over the shell's mean. `tested` is false where the added class could not +// be identified in the written cell (no parity class whose sub-lattice is the half cell) or too few of +// its reflections were measured. +struct DoubledClassTranslation { + struct Shell { + float d_low_A = 0.0f, d_high_A = 0.0f; + int64_t n = 0; + double mean = 0.0, modulated_share = 0.0; + }; + bool tested = false; + std::array parity{}; + std::array eps{}; + std::vector shells; +}; +DoubledClassTranslation FitDoubledClassTranslation(const std::vector &merged, + const UnitCell &doubled_cell, const UnitCell &half_cell); // Whether two passes hold the same lattice: the same Bravais class and primitive volumes within 2 %. bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b); @@ -769,6 +794,8 @@ class Rugnux { 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_; + // Keep the written merged reflections on the result (ProcessResult::merged_reflections). + bool keep_merged_reflections_ = false; // 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, diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index f3d7b8091..923acff63 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1415,27 +1415,24 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro // // 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. +// and its SIZE did not tell the crystal deposited in the doubled cell from those deposited in the half +// cell - the odd class reached 0.4-0.85 of the lattice at high resolution on three of them. Two things +// about it did, each a significance test rather than a bound: +// - its SHAPE with resolution. A pseudo-translation by half a lattice vector leaves the two halves of +// the cell alike at low resolution and different at high, so the class grows from nearly nothing +// outwards; the half cell is then right at low resolution and wrong at high. Classes that stay a +// flat tenth of the lattice at every resolution, or are strongest at low resolution, are not that. +// - whether one translation EXPLAINS it (FitDoubledClassTranslation). Copies related by a translation +// put the added class where sin^2(pi h.eps) is large and nothing where it is small. On the doubled +// deposition it is all modulation (the part left over -0.003 +- 0.005 of the lattice); on the three +// rising or strong classes whose depositions kept the half cell 0.08-0.31 is there whatever the +// direction, 32-64 standard errors from zero - copies that differ in content, or intensity that is +// not the crystal's own. 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) { @@ -1447,18 +1444,157 @@ bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) { 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(); +std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) { + // From the coarsest shell to the finest one both are measured in, the class has to grow by three + // standard errors of the difference: a significant rise, not a size. + if (sc.trend.empty() || !DoubledCellShellMeasured(sc.trend.front())) + return {}; + const auto &low = sc.trend.front(); 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) + const double se = std::hypot(low.excess / low.z, outer->excess / outer->z); + if (outer == &low || outer->excess - low.excess < 3.0 * se) 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); + return fmt::format("the added class grows from {:.3f} of the lattice's intensity at {:.0f}-{:.0f} A to {:.3f} at " + "{:.1f}-{:.1f} A ({:.0f} standard errors): the half cell holds at low resolution and not at high", + low.excess, low.d_low_A, low.d_high_A, outer->excess, outer->d_low_A, outer->d_high_A, + (outer->excess - low.excess) / se); +} + +DoubledClassTranslation FitDoubledClassTranslation(const std::vector &merged, + const UnitCell &doubled_cell, const UnitCell &half_cell) { + DoubledClassTranslation out; + // The added class in the written cell: the parity whose sub-lattice is the half cell. Its real lattice + // is the doubled one plus half of (p0 a + p1 b + p2 c), so one vector the class is odd on is replaced + // by that half-vector. + const CrystalLattice dl(doubled_cell); + const Coord v[3] = {dl.Vec0(), dl.Vec1(), dl.Vec2()}; + const UnitCell want = CrystalLattice(half_cell).NiggliReduce().GetUnitCell(); + double best_mismatch = INFINITY; + for (int c = 1; c < 8; c++) { + const std::array p = {(c >> 2) & 1, (c >> 1) & 1, c & 1}; + const int i = p[0] ? 0 : p[1] ? 1 : 2; + Coord t[3] = {v[0], v[1], v[2]}; + t[i] = (v[0] * static_cast(p[0]) + v[1] * static_cast(p[1]) + v[2] * static_cast(p[2])) * 0.5f; + const UnitCell got = CrystalLattice(t[0], t[1], t[2]).NiggliReduce().GetUnitCell(); + const double mismatch = std::max({std::fabs(got.a / want.a - 1.0), std::fabs(got.b / want.b - 1.0), + std::fabs(got.c / want.c - 1.0), + std::fabs(got.alpha - want.alpha) / 90.0, std::fabs(got.beta - want.beta) / 90.0, + std::fabs(got.gamma - want.gamma) / 90.0}); + if (mismatch < best_mismatch) { + best_mismatch = mismatch; + out.parity = p; + } + } + // A sub-lattice that is the half cell to within a percent; nothing else is the class. + if (best_mismatch > 0.01) + return out; + + // Each added reflection in units of the mean of the lattice's own reflections in its shell + // (SUPERCELL_TREND_D_A, the shells the probe reads), in the shells where that mean is measured - + // three standard errors from zero - and that hold at least MIN_SHELL added reflections. + constexpr int64_t MIN_SHELL = 200; + const auto odd = [&](const MergedReflection &r) { + return (std::abs(out.parity[0] * r.h + out.parity[1] * r.k + out.parity[2] * r.l) & 1) != 0; + }; + const auto shell_of = [](float d) { + if (!(d <= SUPERCELL_TREND_D_A[0])) + return -1; + int sh = 0; + while (d <= SUPERCELL_TREND_D_A[sh + 1]) + sh++; + return sh; + }; + std::array sum{}, sum2{}; + std::array cnt{}, cnt_odd{}; + for (const auto &r : merged) { + const int sh = shell_of(r.d); + if (sh < 0 || !std::isfinite(r.I)) + continue; + if (odd(r)) + cnt_odd[sh]++; + else { + sum[sh] += r.I; sum2[sh] += static_cast(r.I) * r.I; cnt[sh]++; + } + } + std::array use{}; + for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) { + if (cnt[sh] < 2 || cnt_odd[sh] < MIN_SHELL) + continue; + const double n = static_cast(cnt[sh]), m = sum[sh] / n; + use[sh] = m > 3.0 * std::sqrt(std::max(0.0, sum2[sh] / n - m * m) / n); + } + struct Row { double h, k, l, y; int shell; }; + std::vector rows; + for (const auto &r : merged) { + const int sh = shell_of(r.d); + if (sh >= 0 && use[sh] && std::isfinite(r.I) && odd(r)) + rows.push_back({static_cast(r.h), static_cast(r.k), static_cast(r.l), + r.I / (sum[sh] / static_cast(cnt[sh])), sh}); + } + if (rows.size() < 1000) + return out; + + // Least squares of y on sin^2(pi h.eps) and a constant B, for one eps, over the reflections of one + // shell (or all of them, shell < 0); `stride` reads every stride-th reflection only - the coarse + // grid needs the correlation's shape, not its last digit. Returns the correlation, which the eps + // search maximises, and gives the shell's mean y and B. + const auto fit = [&](const std::array &e, size_t stride, int shell, double *mean, double *b, + double *n_out) { + double sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, n = 0; + for (size_t j = 0; j < rows.size(); j += stride) { + const auto &r = rows[j]; + if (shell >= 0 && r.shell != shell) + continue; + n += 1; + const double x = std::pow(std::sin(PI * (r.h * e[0] + r.k * e[1] + r.l * e[2])), 2); + sx += x; sy += r.y; sxx += x * x; sxy += x * r.y; syy += r.y * r.y; + } + const double vxx = sxx - sx * sx / n, vxy = sxy - sx * sy / n, vyy = syy - sy * sy / n; + if (n < 3 || vxx <= 0.0 || vyy <= 0.0) + return -1.0; + if (mean) { + *mean = sy / n; + *b = (sy - vxy / vxx * sx) / n; + *n_out = n; + } + return vxy / std::sqrt(vxx * vyy); + }; + // eps from a coarse grid over the translations within an eighth of a cell of the half point, then + // two finer grids about the best. + std::array best{}; + double best_cc = -INFINITY; + const size_t coarse_stride = std::max(1, rows.size() / 4000); + for (const auto &[step, reach] : {std::pair{0.01, 0.12}, std::pair{0.002, 0.012}, std::pair{0.0005, 0.003}}) { + const size_t stride = step == 0.01 ? coarse_stride : 1; + const std::array centre = best; + const int nstep = static_cast(std::lround(reach / step)); + for (int x = -nstep; x <= nstep; x++) + for (int y = -nstep; y <= nstep; y++) + for (int z = -nstep; z <= nstep; z++) { + const std::array e = {centre[0] + x * step, centre[1] + y * step, centre[2] + z * step}; + const double cc = fit(e, stride, -1, nullptr, nullptr, nullptr); + if (cc > best_cc) { + best_cc = cc; + best = e; + } + } + } + // The share of each shell's added class that the translation modulates: 1 - B over the class's mean. + for (int sh = 0; sh < SUPERCELL_TREND_SHELLS; sh++) { + if (!use[sh]) + continue; + double mean = 0, b = 0, n = 0; + if (fit(best, 1, sh, &mean, &b, &n) < -0.5 || !(mean > 0.0)) + continue; + out.shells.push_back({SUPERCELL_TREND_D_A[sh], SUPERCELL_TREND_D_A[sh + 1], static_cast(n), + mean, 1.0 - b / mean}); + } + out.eps = best; + out.tested = !out.shells.empty(); + return out; } namespace { @@ -1474,8 +1610,6 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes 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. @@ -1487,32 +1621,21 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes doubled.lattice = doubled.search_result.conventional; doubled.extra_lattices.clear(); doubled.unconstrained.reset(); - logger.Info("Doubled-cell hypothesis: lattice {} ({}-centred {}) doubled to {} ({}-centred {})", + 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 as well - lattice {} ({}-centred " + "{}) doubled to {} ({}-centred {})", sc.h, sc.k, sc.l, 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; + + const std::string half_prefix = base_prefix + "_halfcell"; + const std::string doubled_prefix = base_prefix + "_doubled"; + config_.output_prefix = base_prefix.empty() ? base_prefix : doubled_prefix; force_rotation_result_ = doubled; + keep_merged_reflections_ = true; // 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(); @@ -1524,24 +1647,74 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes logger.Warning("Doubled-cell hypothesis: the pass in the doubled cell did not complete ({})", e.what()); } force_rotation_result_.reset(); + keep_merged_reflections_ = false; 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; + // Decided with both merges in hand. The doubled cell is the answer where the added class is a + // pseudo-translation breaking down - it grows with resolution (DoubledCellRises) - and where, in the + // doubled merge, that one translation accounts for most of the class in every shell: the two halves + // of the cell are then one content displaced, out to the limit (FitDoubledClassTranslation). + const std::string rises = DoubledCellRises(sc); + const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell, + *half.consensus_cell); + d.merged_reflections.clear(); + const DoubledClassTranslation::Shell *weakest = nullptr; + std::string table; + for (const auto &sh : tr.shells) { + table += fmt::format(" {:.1f}-{:.1f} A {:.2f};", sh.d_low_A, sh.d_high_A, sh.modulated_share); + if (!weakest || sh.modulated_share < weakest->modulated_share) + weakest = &sh; + } + if (tr.tested) + logger.Info("Doubled-cell hypothesis: the added class {} {} {} of the doubled merge against one pseudo-" + "translation, half a cell plus ({:.4f}, {:.4f}, {:.4f}) - the share of it the translation " + "modulates, by shell:{}", tr.parity[0], tr.parity[1], tr.parity[2], tr.eps[0], tr.eps[1], + tr.eps[2], table); + const bool one_translation = tr.tested && weakest->modulated_share >= 0.5; + const bool adopt = !rises.empty() && one_translation; 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.decision = adopt + ? fmt::format("{}; and one pseudo-translation accounts for most of it in every shell (at least {:.2f}, " + "at {:.1f}-{:.1f} A)", rises, weakest->modulated_share, weakest->d_low_A, weakest->d_high_A) + : rises.empty() + ? std::string("the probed cell stands: the added class does not grow significantly with resolution, " + "so the half cell describes the crystal as well at high resolution as at low") + : !tr.tested + ? std::string("the probed cell stands: the added class could not be read in the doubled merge") + : fmt::format("the probed cell stands: one pseudo-translation accounts for only {:.2f} of the added class " + "at {:.1f}-{:.1f} A - the rest is there whatever the direction, which copies related by a " + "translation do not put there", weakest->modulated_share, weakest->d_low_A, + weakest->d_high_A); + logger.Info("Doubled-cell hypothesis: {}", h.decision); + + // The answer's files go under the plain prefix and the other cell's beside it: where the doubled cell + // is adopted, every file it wrote takes the plain name, and the probed cell's file of the same name + // moves to _halfcell. + if (adopt && !base_prefix.empty()) { + const std::filesystem::path dp(doubled_prefix); + const std::filesystem::path dir = dp.has_parent_path() ? dp.parent_path() : std::filesystem::path("."); + const std::string stem = dp.filename().string(); + std::error_code ec; + std::vector suffixes; + for (const auto &f : std::filesystem::directory_iterator(dir, ec)) + if (const std::string name = f.path().filename().string(); name.starts_with(stem)) + suffixes.push_back(name.substr(stem.size())); + for (const auto &x : suffixes) { + if (std::filesystem::exists(base_prefix + x, ec)) + std::filesystem::rename(base_prefix + x, half_prefix + x, ec); + std::filesystem::rename(doubled_prefix + x, base_prefix + x, ec); + } + } + + ProcessResult &answer = adopt ? d : half; + const ProcessResult &other = adopt ? half : d; 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(); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 4b3b058c0..ac30c40e7 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3974,6 +3974,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, em_report, result.twinning, config_.provenance, config_.output_prefix, static_cast(std::max(1, config_.nthreads)), &sm.scaled_fulls); + if (keep_merged_reflections_) + result.merged_reflections = sm.merged; // The twin hypothesis beside it (see prepare_twin_file): the reflection files only, with the // law and its fraction in the mmCIF's _pdbx_reflns_twin. if (twin_file) { diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 60f00785b..7e7dbdaff 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include #include @@ -371,16 +372,57 @@ TEST_CASE("DoubledCellDecision", "[process]") { // 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. + CHECK(!DoubledCellRises(breaks_down).empty()); + // A class that stays a tenth of the lattice at every resolution: merged as well, but it does not rise. 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()); + CHECK(DoubledCellRises(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("FitDoubledClassTranslation", "[process]") { + // A doubled orthorhombic cell whose two halves are one copy translated by c/2 + (0.03, 0, 0): the + // reflections with l odd carry |F|^2 sin^2(pi h.eps), the others cos^2. With nothing else on the + // added class the fit finds the translation and leaves nothing over; with a flat 0.3 added to + // the class, most of the class is left unexplained. + const UnitCell doubled{50, 60, 140, 90, 90, 90}, half{50, 60, 70, 90, 90, 90}; + std::mt19937 rng(1); + std::exponential_distribution wilson(1.0); + std::vector clean, contaminated; + for (int h = -25; h <= 25; h++) + for (int k = 0; k <= 30; k++) + for (int l = 0; l <= 70; l++) { + if (h == 0 && k == 0 && l == 0) + continue; + MergedReflection r; + r.h = h; r.k = k; r.l = l; + r.d = 1.0f / std::sqrt(static_cast(h * h / 2500.0 + k * k / 3600.0 + l * l / 19600.0)); + if (r.d < 2.0f) + continue; + const double base = wilson(rng), x = std::pow(std::sin(PI * 0.03 * h), 2); + r.I = static_cast(base * ((l & 1) ? x : 1.0 - x)); + clean.push_back(r); + if (l & 1) + r.I += 0.3f; + contaminated.push_back(r); + } + const DoubledClassTranslation a = FitDoubledClassTranslation(clean, doubled, half); + REQUIRE(a.tested); + CHECK(a.parity == std::array{0, 0, 1}); + const DoubledClassTranslation b = FitDoubledClassTranslation(contaminated, doubled, half); + REQUIRE(b.tested); + const auto weakest = [](const DoubledClassTranslation &t) { + double w = INFINITY; + for (const auto &sh : t.shells) + w = std::min(w, sh.modulated_share); + return w; + }; + CHECK(weakest(a) > 0.8); + CHECK(weakest(b) < 0.5); +} + TEST_CASE("RefinedPassIsWorse", "[process]") { const auto pass = [](const UnitCell &cell, char centering, int64_t strong, int64_t negative, int64_t reflections) { From 57543845c01b1add6e7eafdb9444f1e316ced8bc Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 14:00:39 +0200 Subject: [PATCH 3/8] WIP: doubled-cell pass only where the class rises; no model validation in it A class that does not rise keeps the half cell whatever the doubled merge shows, so the doubled pass is no longer run for it - on a ci-tier crystal (R-centred, class at 1-2% of the lattice) that pass ran the card out of memory and failed the run. The doubled pass also runs without --model: placing the model in a doubled cubic crystal's 776 A R-centred cell took twenty minutes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/RUGNUX_REPORT.md | 24 +++++++++++++----------- rugnux/RugnuxPasses.cpp | 20 +++++++++++++++----- 2 files changed, 28 insertions(+), 16 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 0695376f5..74166f671 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -618,19 +618,21 @@ spots of further lattice domains can land on the half-integer positions too. A ` 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. +three shells measured above its noise, at more than twice the strongest of them in standard errors) and +grows significantly with resolution - the half cell holds at low resolution and not at high, the shape of +a pseudo-translation breaking down - the run integrates and merges the doubled cell as well, on the same +orientation and geometry, and decides with both merges in hand. The size of the class does not decide +it - crystals deposited in the half cell carry classes as strong as one deposited in the doubled cell. +The doubled cell is the answer only where, in the doubled merge, one translation - half a cell plus a +small offset, fitted over the class - modulates at least half of the added class in every resolution +shell, so the two halves of the cell are one content displaced. The doubled pass runs without `--model`. A class that is there whatever the +direction (copies that differ in content, twin domains, neighbouring reflections) 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. +`_halfcell.mtz` where the doubled cell was adopted, with the rest of its files under the same +prefix), and **`DOUBLED_CELL_OTHER_CELL=`**, **`DOUBLED_CELL_OTHER_SPACE_GROUP=`**, +**`DOUBLED_CELL_OTHER_D_MIN=`**, **`DOUBLED_CELL_OTHER_CC_HALF=`** and, where it is the probed cell and +`--model` was given, **`DOUBLED_CELL_OTHER_R_FREE=`** describe it. Compare the two in refinement. ## Further lattices diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 923acff63..0b16e0b6f 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1610,6 +1610,16 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes return half; Logger logger("Rugnux"); const auto &sc = *half.supercell; + // Only a class that could be adopted is worth the pass: one that does not rise with resolution keeps + // the half cell whatever the doubled merge shows, and a doubled merge of a large crystal is not free + // (measured: an R-centred 531 A cell ran the card out of memory). + const std::string rises = DoubledCellRises(sc); + if (rises.empty()) { + logger.Info("Doubled-cell hypothesis: the half-integer class {} {} {} stands clear of the other classes but " + "does not grow significantly with resolution, so the half cell describes the crystal as well " + "at high resolution as at low - not merged", sc.h, sc.k, sc.l); + return half; + } // 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. @@ -1636,6 +1646,9 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes config_.output_prefix = base_prefix.empty() ? base_prefix : doubled_prefix; force_rotation_result_ = doubled; keep_merged_reflections_ = true; + // No model validation in this pass: a model is placed in a cell, and on a doubled cubic crystal + // (an R-centred 776 A cell) placing it took twenty minutes for a pass that decides nothing by it. + const std::string model_path = std::exchange(config_.model_path, std::string()); // 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(); @@ -1648,6 +1661,7 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes } force_rotation_result_.reset(); keep_merged_reflections_ = false; + config_.model_path = model_path; 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"); @@ -1660,7 +1674,6 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes // pseudo-translation breaking down - it grows with resolution (DoubledCellRises) - and where, in the // doubled merge, that one translation accounts for most of the class in every shell: the two halves // of the cell are then one content displaced, out to the limit (FitDoubledClassTranslation). - const std::string rises = DoubledCellRises(sc); const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell, *half.consensus_cell); d.merged_reflections.clear(); @@ -1677,15 +1690,12 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes "modulates, by shell:{}", tr.parity[0], tr.parity[1], tr.parity[2], tr.eps[0], tr.eps[1], tr.eps[2], table); const bool one_translation = tr.tested && weakest->modulated_share >= 0.5; - const bool adopt = !rises.empty() && one_translation; + const bool adopt = one_translation; ProcessResult::DoubledCell h; h.adopted = adopt; h.decision = adopt ? fmt::format("{}; and one pseudo-translation accounts for most of it in every shell (at least {:.2f}, " "at {:.1f}-{:.1f} A)", rises, weakest->modulated_share, weakest->d_low_A, weakest->d_high_A) - : rises.empty() - ? std::string("the probed cell stands: the added class does not grow significantly with resolution, " - "so the half cell describes the crystal as well at high resolution as at low") : !tr.tested ? std::string("the probed cell stands: the added class could not be read in the doubled merge") : fmt::format("the probed cell stands: one pseudo-translation accounts for only {:.2f} of the added class " From 7580256ac71d738bbefeeee0c2a71d98720792d8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 14:14:10 +0200 Subject: [PATCH 4/8] WIP: doubled-cell rise must end at a quarter of the lattice's intensity On strong data a class at one or two percent of the lattice rises significantly (0.009 to 0.059 at 4 standard errors on a ci-tier triclinic crystal, which was then adopted doubled and failed), and a half cell missing a few percent of its intensity still describes the crystal. The doubled deposition ends at 0.35-0.40, those classes at 0.06. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/RUGNUX_REPORT.md | 5 +++-- rugnux/RugnuxPasses.cpp | 14 ++++++++++---- tests/RugnuxTest.cpp | 3 +++ 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 74166f671..bfe6f1c9d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -619,8 +619,9 @@ under the `SUPERCELL_POSSIBLE` flag, and the summary's `Supercell` line carries **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) and -grows significantly with resolution - the half cell holds at low resolution and not at high, the shape of -a pseudo-translation breaking down - the run integrates and merges the doubled cell as well, on the same +grows significantly with resolution to a quarter of the lattice's intensity or more in the finest shell +both are measured in - the half cell holds at low resolution and not at high, the shape of a +pseudo-translation breaking down - the run integrates and merges the doubled cell as well, on the same orientation and geometry, and decides with both merges in hand. The size of the class does not decide it - crystals deposited in the half cell carry classes as strong as one deposited in the doubled cell. The doubled cell is the answer only where, in the doubled merge, one translation - half a cell plus a diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 0b16e0b6f..c2e229fbf 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1446,7 +1446,13 @@ bool DoubledCellWorthMerging(const ProcessResult::SupercellClass &sc) { std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) { // From the coarsest shell to the finest one both are measured in, the class has to grow by three - // standard errors of the difference: a significant rise, not a size. + // standard errors of the difference, and end at a quarter of the lattice's intensity or more. The + // rise alone is not enough: on strong data a class at one or two percent of the lattice rises + // significantly (measured: 0.009 to 0.059 at 4 standard errors on a triclinic crystal, and 0.000 to + // 0.060 on a cubic one), and a half cell missing a few percent of the intensity still describes the + // crystal. The quarter sits far from both sides: the doubled deposition ends at 0.35-0.40, those two + // at 0.06. + constexpr double DOUBLED_CELL_OUTER_EXCESS = 0.25; if (sc.trend.empty() || !DoubledCellShellMeasured(sc.trend.front())) return {}; const auto &low = sc.trend.front(); @@ -1455,7 +1461,7 @@ std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) { if (DoubledCellShellMeasured(sh)) outer = &sh; const double se = std::hypot(low.excess / low.z, outer->excess / outer->z); - if (outer == &low || outer->excess - low.excess < 3.0 * se) + if (outer == &low || outer->excess - low.excess < 3.0 * se || outer->excess < DOUBLED_CELL_OUTER_EXCESS) return {}; return fmt::format("the added class grows from {:.3f} of the lattice's intensity at {:.0f}-{:.0f} A to {:.3f} at " "{:.1f}-{:.1f} A ({:.0f} standard errors): the half cell holds at low resolution and not at high", @@ -1616,8 +1622,8 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes const std::string rises = DoubledCellRises(sc); if (rises.empty()) { logger.Info("Doubled-cell hypothesis: the half-integer class {} {} {} stands clear of the other classes but " - "does not grow significantly with resolution, so the half cell describes the crystal as well " - "at high resolution as at low - not merged", sc.h, sc.k, sc.l); + "does not grow with resolution to a quarter of the lattice's intensity, so the half cell describes " + "the crystal at high resolution too - not merged", sc.h, sc.k, sc.l); return half; } diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 7e7dbdaff..f8292f749 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -380,6 +380,9 @@ TEST_CASE("DoubledCellDecision", "[process]") { // 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)); + // A class that rises significantly but stays a few percent of the lattice: the half cell holds. + const auto small = trend({{0.009, 40, 500}, {0.014, 60, 300}, {0.016, 90, 170}, {0.03, 60, 50}, {0.06, 40, 2}}); + CHECK(DoubledCellRises(small).empty()); } TEST_CASE("FitDoubledClassTranslation", "[process]") { From 757ca02dbff8653568084d6a5812012018ff20aa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 14:16:33 +0200 Subject: [PATCH 5/8] WIP: doubled-cell class identified against the half cell's primitive lattice A centred half cell (I-centred orthorhombic, its centring class the added one) was compared to the doubled cell's sub-lattices as if primitive, so its class was never read. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/Rugnux.h | 3 ++- rugnux/RugnuxPasses.cpp | 8 +++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 21227a34b..fc61f484d 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -678,7 +678,8 @@ struct DoubledClassTranslation { std::vector shells; }; DoubledClassTranslation FitDoubledClassTranslation(const std::vector &merged, - const UnitCell &doubled_cell, const UnitCell &half_cell); + const UnitCell &doubled_cell, const UnitCell &half_cell, + char half_centering = 'P'); // Whether two passes hold the same lattice: the same Bravais class and primitive volumes within 2 %. bool ArmsHoldSameLattice(const ProcessResult &a, const ProcessResult &b); diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index c2e229fbf..eb0db145d 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1470,14 +1470,15 @@ std::string DoubledCellRises(const ProcessResult::SupercellClass &sc) { } DoubledClassTranslation FitDoubledClassTranslation(const std::vector &merged, - const UnitCell &doubled_cell, const UnitCell &half_cell) { + const UnitCell &doubled_cell, const UnitCell &half_cell, + char half_centering) { DoubledClassTranslation out; // The added class in the written cell: the parity whose sub-lattice is the half cell. Its real lattice // is the doubled one plus half of (p0 a + p1 b + p2 c), so one vector the class is odd on is replaced // by that half-vector. const CrystalLattice dl(doubled_cell); const Coord v[3] = {dl.Vec0(), dl.Vec1(), dl.Vec2()}; - const UnitCell want = CrystalLattice(half_cell).NiggliReduce().GetUnitCell(); + const UnitCell want = CrystalLattice(half_cell).ToPrimitive(half_centering).NiggliReduce().GetUnitCell(); double best_mismatch = INFINITY; for (int c = 1; c < 8; c++) { const std::array p = {(c >> 2) & 1, (c >> 1) & 1, c & 1}; @@ -1681,7 +1682,8 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes // doubled merge, that one translation accounts for most of the class in every shell: the two halves // of the cell are then one content displaced, out to the limit (FitDoubledClassTranslation). const DoubledClassTranslation tr = FitDoubledClassTranslation(d.merged_reflections, *d.consensus_cell, - *half.consensus_cell); + *half.consensus_cell, + half.consensus_centering.value_or('P')); d.merged_reflections.clear(); const DoubledClassTranslation::Shell *weakest = nullptr; std::string table; From 11a8f49555da15597a83d75f237eeea55ccf6192 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 14:42:29 +0200 Subject: [PATCH 6/8] WIP: doubled cell not adopted where its merge is more symmetric than the half cell's A superstructure cannot be more symmetric than its own average, which the half cell's merge is. On a private-arm crystal whose added class passed the rise and the one-translation tests, the doubled merge came out centrosymmetric over a chiral half-cell group (point group order 8 against 4) and the run failed against its reference. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/RugnuxPasses.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index eb0db145d..c196c1731 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1698,12 +1698,23 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes "modulates, by shell:{}", tr.parity[0], tr.parity[1], tr.parity[2], tr.eps[0], tr.eps[1], tr.eps[2], table); const bool one_translation = tr.tested && weakest->modulated_share >= 0.5; - const bool adopt = one_translation; + // A superstructure cannot be more symmetric than its own average: the half cell's merge is that + // average, so a doubled merge that reads a point group of higher order than the half cell's describes + // something other than the crystal's superstructure (measured: a C2221 half cell whose doubled merge + // came out centrosymmetric, Cccm, with the added class otherwise passing every test). + const auto point_group_order = [](const ProcessResult &r) { + return r.space_group ? r.space_group->operations().sym_ops.size() : size_t{0}; + }; + const bool symmetry_raised = point_group_order(d) > point_group_order(half); + const bool adopt = one_translation && !symmetry_raised; ProcessResult::DoubledCell h; h.adopted = adopt; h.decision = adopt ? fmt::format("{}; and one pseudo-translation accounts for most of it in every shell (at least {:.2f}, " "at {:.1f}-{:.1f} A)", rises, weakest->modulated_share, weakest->d_low_A, weakest->d_high_A) + : symmetry_raised + ? fmt::format("the probed cell stands: the doubled merge reads {}, more symmetric than the half cell's {}, " + "which a superstructure of it cannot be", d.space_group->xhm(), half.space_group->xhm()) : !tr.tested ? std::string("the probed cell stands: the added class could not be read in the doubled merge") : fmt::format("the probed cell stands: one pseudo-translation accounts for only {:.2f} of the added class " From d7001b59f38e81d46084f0487e18c2f2cae8a9e1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 15:08:03 +0200 Subject: [PATCH 7/8] WIP: doubled-cell shells count only where the geometry resolves most of the class A shell is read only where at least half of the class's reflections are clear of the lattice's by r3. On a 701 A axis 26-43 of ~36000 per shell were clear; they read at the lattice's own strength, the doubled merge (93% of predictions overlapping a neighbour) passed the one-translation test, and the run was adopted doubled and failed. With the rule, of ~134 sets with probe logs only 9min and one I-centred crystal reach the doubled pass. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- rugnux/Rugnux.h | 6 +++++- rugnux/RugnuxImagePass.cpp | 1 + rugnux/RugnuxPasses.cpp | 9 +++++++-- tests/RugnuxTest.cpp | 1 + 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index fc61f484d..f4774fbbb 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -446,11 +446,15 @@ struct ProcessResult { 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 . + // errors (z), the largest such z among the six beside it, the lattice's own , and the + // share of the class resolved from the lattice. 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; + // The class's reflections clear of the lattice's by r3, per lattice reflection: how much of + // the class this geometry resolves at all. + double resolved = 0.0; }; std::vector trend; }; diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 63c185b13..8c018f395 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -468,6 +468,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { 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.resolved = static_cast(x.n) / static_cast(l.n); shell.excess = (mean(x) - null_i) / mean(l); shell.z = (mean(x) - null_i) / se(x); shell.sibling_z = -INFINITY; diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index c196c1731..c8cefa501 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1429,9 +1429,14 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro // direction, 32-64 standard errors from zero - copies that differ in content, or intensity that is // not the crystal's own. namespace { - // A shell the class is read in: measured above its noise, against a lattice that is measured too. + // A shell the class is read in: measured above its noise, against a lattice that is measured too, in a + // geometry that resolves most of it from the lattice. Where the added reflections mostly fall within + // the background radius of a lattice reflection, what the few clear ones read is not the class, and a + // doubled merge integrates the rest on top of their neighbours (measured: a 701 A axis left 26-43 of + // ~36000 clear per shell, read them at the lattice's own strength, and its doubled merge passed every + // other test). bool DoubledCellShellMeasured(const ProcessResult::SupercellClass::Shell &sh) { - return sh.z >= 3.0 && sh.lattice_i_over_sigma >= 1.0; + return sh.z >= 3.0 && sh.lattice_i_over_sigma >= 1.0 && sh.resolved >= 0.5; } } diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index f8292f749..9cee8fde9 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -364,6 +364,7 @@ TEST_CASE("DoubledCellDecision", "[process]") { sh.z = z; sh.sibling_z = 3.0; sh.lattice_i_over_sigma = lattice; + sh.resolved = 0.9; sc.trend.push_back(sh); } return sc; From b99cd3b5c6504c982c72dce529f0d7364cefe5e4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 15:16:58 +0200 Subject: [PATCH 8/8] rugnux doubled cell: model validation back in the doubled pass; docs for the final rules The doubled pass now runs only on resolved, rising classes (two sets of ~134 with probe logs), so the twenty-minute model placement in a doubled cubic cell can no longer arise, and the adopted answer is validated again (9min: R-free 0.275 against the deposited model, 0.574 in the half cell). Final rules, measured on gate + ci tier + all sets with a de-modulated axial row + 49 sets with a baseline supercell occupancy >= 8% + the private arm: 9min fail -> pass (P212121, doubled cell); every other set run keeps its verdict and a bit-identical p.hkl. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/RUGNUX_REPORT.md | 10 ++++++---- rugnux/RugnuxPasses.cpp | 5 +---- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index bfe6f1c9d..e7a6d5fbf 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -618,7 +618,8 @@ spots of further lattice domains can land on the half-integer positions too. A ` 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) and +three shells where it is measured above its noise and most of its reflections are resolved from the +lattice's, at more than twice the strongest of them in standard errors) and grows significantly with resolution to a quarter of the lattice's intensity or more in the finest shell both are measured in - the half cell holds at low resolution and not at high, the shape of a pseudo-translation breaking down - the run integrates and merges the doubled cell as well, on the same @@ -626,14 +627,15 @@ orientation and geometry, and decides with both merges in hand. The size of the it - crystals deposited in the half cell carry classes as strong as one deposited in the doubled cell. The doubled cell is the answer only where, in the doubled merge, one translation - half a cell plus a small offset, fitted over the class - modulates at least half of the added class in every resolution -shell, so the two halves of the cell are one content displaced. The doubled pass runs without `--model`. A class that is there whatever the +shell, so the two halves of the cell are one content displaced - and where the doubled merge is not more +symmetric than the half cell's, which a superstructure of it cannot be. A class that is there whatever the direction (copies that differ in content, twin domains, neighbouring reflections) 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 the rest of its files under the same prefix), and **`DOUBLED_CELL_OTHER_CELL=`**, **`DOUBLED_CELL_OTHER_SPACE_GROUP=`**, -**`DOUBLED_CELL_OTHER_D_MIN=`**, **`DOUBLED_CELL_OTHER_CC_HALF=`** and, where it is the probed cell and -`--model` was given, **`DOUBLED_CELL_OTHER_R_FREE=`** describe it. Compare the two in refinement. +**`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/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index c8cefa501..e00a0f758 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1658,9 +1658,7 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes config_.output_prefix = base_prefix.empty() ? base_prefix : doubled_prefix; force_rotation_result_ = doubled; keep_merged_reflections_ = true; - // No model validation in this pass: a model is placed in a cell, and on a doubled cubic crystal - // (an R-centred 776 A cell) placing it took twenty minutes for a pass that decides nothing by it. - const std::string model_path = std::exchange(config_.model_path, std::string()); + // 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(); @@ -1673,7 +1671,6 @@ ProcessResult Rugnux::DoubledCellHypothesis(RugnuxObserver *observer, ProcessRes } force_rotation_result_.reset(); keep_merged_reflections_ = false; - config_.model_path = model_path; 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");