diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 564d97ee7..4d51dcdf9 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -14,6 +14,8 @@ * Rugnux keeps a second-pass lattice that corrects an axis-multiple supercell of the first pass, instead of reverting to the supercell because its centring differs. * Rugnux settles a disagreement between the file's and the measured beam centre over an axis harmonic in either direction, so a file centre off along the spindle no longer leaves the run on a doubled axis. * Rugnux commits a post-refined geometry only when the held-out residual falls by more than its own noise, instead of by a fixed 2 %. +* Rugnux reports the strongest index-2 superstructure class of a rotation run (`SUPERCELL_CLASS` and related keys) with the doubled cell to give with `-C`, and warns (`SUPERCELL_POSSIBLE`) where that class looks like Bragg reflections, so both settings can be tried; it does not change the lattice. +* Rugnux keeps refining the detector geometry between passes where the fit's new geometry puts more validation spots on the lattice, not only where the held-out residual falls. * Rugnux also reports `R_MEAS_WEIGHTED`, the R_meas with each observation weighted as the merge weights it, so weak frames kept in the merge at low weight can be told from real disagreement; `R_MEAS` is unchanged. * Rugnux no longer stops per-frame scaling as unsettled because of a single frame that is about to be dropped as blank. * `rugnux --polarization` is documented as the polarization degree (XDS `FRACTION_OF_POLARIZATION` = (1 + p)/2); the default 0.99 suits undulators. diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index 4ac0b252d..aa5b55c76 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -331,6 +331,8 @@ The refinement above (§7.2) runs per image against that image's spots. For rota **Whether the data determine the distance at all is asked, not assumed.** At a detector far enough away that no reflection reaches more than a few degrees of $2\theta$, a longer distance and a larger cell move every spot the same way to first order — the difference is of order $\sin^2\theta$ of the spot's own position, about 0.4 px rms per per cent of distance over a 2M detector at 820 mm against 2–3 px at the distances a crystal is usually collected at. The joint fit then finds a distance/cell pair that fits its own spot positions a little better than the header, commits it, and the pass re-integrated there finds the next pair: a walk along the degenerate direction that the realised residual never ratifies (measured on such a sweep: a header at 820 mm walked to 846 mm with the cell 3.3 % too large, the realised held-out residual flat at every round). So the same fit is asked once more with the distance **held at the header**, every other block as free as before — the nested hypothesis "the header distance is right" — and the two are compared on the one residual family that can tell them apart: the **excitation** residual. It never involves the detector, so it is blind to the distance itself; what it sees is the cell scale, and a held fit at a wrong header distance is forced into a wrong cell scale by the spot positions, which the rocking angles then refuse (measured: a header 1.4 % long leaves the held fit's excitation residual seventeen times the free fit's). Where freeing the distance lowers the held-out excitation residual below the held fit's by more than that residual's own standard error, the free fit is committed exactly as before; where it does not, the held fit is — header distance, refined beam, cell, orientation and axis — and the report says so (`POSTREFINE_DISTANCE_HELD`). The positional residual is deliberately not consulted for this: it is the family whose in-fit gain along the degenerate direction re-integration erases, and pooled with the excitation family it either drowns a decisive excitation gain in its own noise (a 54 % excitation gain read as 9 % pooled against a 9 % noise) or lends the degenerate direction a gain that is not there. Nothing is tuned here: the only input is the standard error of the residual itself, the same noise the geometry walk's rounds have to beat. The wavelength is never refined on a single crystal for the same reason in its exact form: it scales the spot positions and the rocking angles identically to the cell, so no sweep can tell the two apart at any $2\theta$. 3. **Pass 2** re-indexes de novo and re-integrates at the committed geometry. Only the **detector distance and beam centre** carry over: the refined cell, orientation and axis are what make the distance identifiable, but pass 2 re-indexes from scratch, so they are not propagated. Where that re-index indexes too few frames the run falls back to pass 1's lattice and integrates it at the refined geometry — and the cell is then **scaled to the distance it will be used at**, since a real-space cell is measured against the distance its spots were seen at, and carrying it across a distance change otherwise scales the whole cell by the ratio of the two. The orientation is untouched. + Pass 2 measures the post-refinement again, and where it still moves the geometry the run **walks**: it re-indexes and re-integrates at what the fit asks for, and repeats. A round is kept only for what it *realises*, not for what the fit predicts, and it can realise a gain in two ways, either of which has to beat its own noise: a lower held-out residual (the standard errors of the two means combined), or a larger share of the validation spots on the lattice beyond chance (the binomial noise of the two shares, z = 3.29). The residual alone misses exactly the errors that cost resolution: it is dominated by the low-resolution reflections, where a distance and the cell scale that compensates it move every spot alike, and its centroids are taken inside a disc centred on the prediction, so they follow the prediction part of the way; the high-resolution validation spots are the first to leave the lattice (measured: 0.6 % of distance read 0.74 of the residual's noise and 70.4 % against 58.1 % of the validation spots, and cost 0.07 Å of resolution). A move of a trust-region step or more starts the walk outright. A smaller one is first tried as two indexing probes on the validation frames — at the fit's geometry and at the one in hand, each stopping once the lattice is scored — and pays for a re-integrated round only where the fit's geometry scores higher. The run keeps the best round it reached. + The space group is determined **after** pass 2, on the geometry the run refined, and pass 1 does not search at all: a decision taken on the worse of the two passes and then carried forward is a constraint on the better one, and would have to be reconciled with what pass 2 later found. The guard that chooses which pass is written compares each pass's **first** merge — $P1$ on both sides, full resolution range, before the correction surfaces — which both passes produce anyway, so it never compares statistics computed in two different space groups. What it compares there is the **signal each pass measured**: the count of unique reflections merged at $I/\sigma \ge 2$. Self-consistency cannot do this job — against an external arbiter the signal count named the more accurate geometry on 23 of 27 arm-dataset pairs where $R_\mathrm{meas}$, $CC_{1/2}$ and ISa managed 13, a coin flip — and that merge's own $CC_{1/2}$ least of all, being pooled over the whole range, uncut and uncorrected, so the shells with no signal in them dominate it and they are exactly the shells a geometry move disturbs (it reads 0.13 on a crystal whose data merge at 0.995). The refined pass is sent back only where it merges more unique reflections than its cell can hold, where it measured decisively less signal (10 %, and only where it holds no more reflections either — a wider integration disk pulls weak reflections in and dilutes the strong fraction without measuring less), or where it lost the axial rows the systematic absences are read off. One index-time veto remains and is keyed to pass 1's **lattice** rather than its group: a centred pass-1 lattice against a primitive pass-2 one. Only pass 2 is written, as the canonical `_*` output. Pass 1's merge exists to give the guard something to judge pass 2 against, so it stops short of the parts of the merge that only fill in a file — the correction surfaces, the twinning and radiation-damage analyses, the R-free flags and the amplitudes — and writes no merged files of its own. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 4d09807a1..2c47cdf8c 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -148,7 +148,7 @@ and a run that produced a textbook data set read identically for their first thr `SYMMETRY_AMBIGUITY`, `CENTERING_UNTESTED`, `UNUSABLE_MERGE`, `LOW_COMPLETENESS`, `SWEEP_GAPS`, `GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`, `LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `CANCELLED`, `RESOLUTION_FIT`, - `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`. `NONE` when nothing fired. A code appears if and only if its + `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`, `SUPERCELL_POSSIBLE`. `NONE` when nothing fired. A code appears if and only if its warning fired, so the flags and the `WARNING:` lines are two renderings of one list — the closed type for machinery, the open sentence for a person. - **`WARNING_COUNT=`** and the **`WARNING:`** lines follow, in the same section. They are what they @@ -415,6 +415,34 @@ Rings are detected far more often than they are excluded: exclusion happens only found no usable lattice. The measurement is described in [CPU/GPU data analysis ▸ Resolution and ice-ring handling](CPU_DATA_ANALYSIS_IMAGE.md#33-resolution-and-ice-ring-handling). +## Index-2 superstructure + +On rotation data, section 4 reports whether the lattice the run adopted has intensity at half-integer +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. + +**`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. + +**`SUPERCELL_POSSIBLE=`** is `TRUE`, with a warning under the `SUPERCELL_POSSIBLE` flag, 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, 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. The +summary at the top of the report carries the same advice on its `Supercell` line, beside +pseudo-symmetry and twinning. + ## Translational pseudo-symmetry Two copies of the contents of the asymmetric unit related by a pure translation that is not a lattice diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index f4fca96c5..8dea3f1c9 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -640,6 +640,7 @@ void IndexAndRefine::ProbeSupercell(std::vector image_numbers) { std::sort(image_numbers.begin(), image_numbers.end()); supercell_probe_frames_ = std::move(image_numbers); supercell_probe_ = {}; + supercell_probe_primitive_.reset(); } SupercellProbe IndexAndRefine::GetSupercellProbe() { @@ -647,6 +648,11 @@ SupercellProbe IndexAndRefine::GetSupercellProbe() { return supercell_probe_; } +std::optional IndexAndRefine::GetSupercellProbePrimitive() { + const std::unique_lock ul(supercell_probe_mutex_); + return supercell_probe_primitive_; +} + void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction &prediction, const BraggIntegrateFn &integrate, const IndexingOutcome &outcome, const BraggPredictionSettings &settings) { @@ -684,6 +690,7 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction for (int i = 0; i < 8; i++) for (int j = 0; j < 2; j++) supercell_probe_[i][j] += frame[i][j]; + supercell_probe_primitive_ = prim; } std::optional diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 7d5537fe0..f3b30055c 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -123,6 +123,7 @@ class IndexAndRefine { // Supercell probe: the image numbers it reads, and what it has read (see ProbeSupercell). std::vector supercell_probe_frames_; SupercellProbe supercell_probe_{}; + std::optional supercell_probe_primitive_; std::mutex supercell_probe_mutex_; void ProbeSupercellFrame(const DataMessage &msg, BraggPrediction &prediction, const BraggIntegrateFn &integrate, const IndexingOutcome &outcome, @@ -189,6 +190,9 @@ 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(); + // The primitive lattice the parity classes are indexed in (the last 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(); // Index a single frame (no integration) with the current forced rotation lattice; used to score // first-pass sampling schemes on the real per-image path. Returns whether the frame indexed. bool IndexFrameOnly(DataMessage &msg, const SpotFindingSettings &settings); diff --git a/rugnux/ReportDocument.h b/rugnux/ReportDocument.h index be8690150..3b72acc0f 100644 --- a/rugnux/ReportDocument.h +++ b/rugnux/ReportDocument.h @@ -43,6 +43,7 @@ namespace PathologyCode { // in how exact it is. ANISOTROPY and RESOLUTION_FIT are already shared this way. constexpr const char *LATTICE_TRANSLATION = "LATTICE_TRANSLATION"; // a translation the cell does not declare constexpr const char *HARMONIC_CONTAMINATION = "HARMONIC_CONTAMINATION"; // the beam carries a higher harmonic + constexpr const char *SUPERCELL_POSSIBLE = "SUPERCELL_POSSIBLE"; // Bragg-like intensity at half-integer positions constexpr const char *GEOMETRY_NOT_CONVERGED = "GEOMETRY_NOT_CONVERGED"; // the geometry walk ran out of rounds constexpr const char *SCALING_NOT_CONVERGED = "SCALING_NOT_CONVERGED"; // the per-frame scales hit their cap still moving } diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 5f485d3de..f56b5e7cb 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -50,6 +50,15 @@ namespace { const char *BANNER = " ******************************************************************************"; + // Advisory, not a decision: the half-integer class is measured ( at least 0.5) and part of + // it rocks like Bragg reflections (at least 2% of the lattice's, three standard errors clear). On + // the battery's 105 rotation sets this names 8 - both crystals whose accepted cell is the doubled + // one, a pseudo-translation whose doubled description is an accepted alternative, and five correct + // sub-cells with real but weak half-integer intensity. + bool SupercellPossible(const ProcessResult::SupercellClass &sc) { + return sc.i_over_sigma >= 0.5 && sc.rock_pct >= 2.0 && sc.rock_pct > 3.0 * sc.rock_se_pct; + } + std::string CellString(const UnitCell &c) { return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c, c.alpha, c.beta, c.gamma); @@ -1075,6 +1084,43 @@ 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. + if (result.supercell) { + const auto &sc = *result.supercell; + Add(s, KeyText("SUPERCELL_CLASS", fmt::format("{} {} {}", sc.h, sc.k, sc.l))); + Add(s, KeyReal("SUPERCELL_OCCUPANCY_PCT", sc.occupancy_pct, "{:.1f}")); + Add(s, KeyReal("SUPERCELL_ROCK_PCT", sc.rock_pct, "{:.1f}")); + Add(s, KeyReal("SUPERCELL_ROCK_SE_PCT", sc.rock_se_pct, "{:.1f}")); + Add(s, KeyReal("SUPERCELL_I_OVER_SIGMA", sc.i_over_sigma, "{:.2f}")); + Add(s, KeyText("SUPERCELL_DOUBLED_CELL", CellString(sc.doubled_cell))); + const bool possible = SupercellPossible(sc); + Add(s, KeyBool("SUPERCELL_POSSIBLE", possible)); + if (possible) + Warn(doc, PathologyCode::SUPERCELL_POSSIBLE, + fmt::format("Half-integer reflections of parity {} {} {} carry Bragg-like intensity ({:.0f}% " + "of the lattice's, rocking part {:.1f}%): the true cell may be twice as large. " + "Check carefully, preferably by processing both settings - this cell and " + "-C \"{:.2f},{:.2f},{:.2f},{:.2f},{:.2f},{:.2f}\" - and comparing them in refinement", + sc.h, sc.k, sc.l, sc.occupancy_pct, sc.rock_pct, + sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c, + sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma)); + 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.")); + Add(s, Blank()); + } + // ---- powder contamination // A crystalline phase other than the crystal, diffracting as rings among its reflections. // Measured in the pre-scan on every run, from the spots found there for the width and the @@ -2009,6 +2055,16 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, : std::string("no indication")); if (result.twinning.l_test_pairs > 0) row("Twinning", TwinningVerdictLine(result.twinning)); + if (result.supercell) + row("Supercell", SupercellPossible(*result.supercell) + ? fmt::format("POSSIBLE - class {} {} {} rocks like Bragg reflections ({:.1f}%); " + "process also with -C \"{:.2f},{:.2f},{:.2f},{:.2f},{:.2f},{:.2f}\" " + "and compare", result.supercell->h, result.supercell->k, + result.supercell->l, result.supercell->rock_pct, + result.supercell->doubled_cell.a, result.supercell->doubled_cell.b, + result.supercell->doubled_cell.c, result.supercell->doubled_cell.alpha, + result.supercell->doubled_cell.beta, result.supercell->doubled_cell.gamma) + : std::string("no indication")); const double db = result.merge_statistics.radiation_damage_delta_b; if (std::isfinite(db)) row("Radiation damage", fmt::format("relative B {:+.2f} A^2 over the sweep", db)); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index fd6207325..4e4251246 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2505,9 +2505,10 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // geometry it did not come from, a wrong header does, and a pass that lost the crystal is thrown // out by the quality guard below, which returns the run to the header geometry. // - // What starts the walk is a move of more than one step: that is the fit saying the geometry is - // somewhere else entirely rather than a fraction of a percent away, and it leaves every run - // whose fit settles beside its header at the two passes it always had. Once the run IS walking + // What starts the walk outright is a move of more than one step: that is the fit saying the + // geometry is somewhere else entirely rather than a fraction of a percent away. A smaller move + // has to be preferred by the validation frames first (below), so a run whose fit settles beside + // the geometry in hand costs two indexing probes and no pass. Once the run IS walking // it keeps walking because the steps do not get smaller in proportion to what is left: measured // on a header 3.9 % long, the first round took a third of the error and each of the next took a // third of the rest. @@ -2525,6 +2526,53 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { // nothing ratifies it. A walk still paying when the rounds run out keeps its best round - the // last - and says it did not converge; the quality guard below judges it against the header // pass as it judges every refined pass. + // + // The residual is not the only thing a round realises. It is the mean over every reflection, + // measured by centroids taken inside a disc centred on the prediction - so they follow the + // prediction part of the way - and dominated by the low-resolution reflections, where a distance + // and the cell scale that compensates it move every spot alike. What does see such an error is + // the re-indexing: the share of the validation spots that land on the lattice, which the + // high-resolution spots lose first. Measured on a crystal whose canonical pass ran 0.6 % long, + // the round at the fit's distance put 70.4 % of the validation spots on the lattice against + // 58.1 %, while reading the residual only 0.74 of its noise lower - and the run the residual + // stopped kept the long distance and lost 0.05 A of resolution. So a round is kept when EITHER + // realises a gain beyond its noise: the residual (HeldOutResidualFell), or the validation + // frames (ValidationEvidencePrefers). + // + // And a move of less than one step no longer settles the geometry by being small: it is asked + // of the validation frames too, but cheaply first - an indexing probe at the fit's geometry + // against one at the geometry in hand, each scoring the lattice and stopping there - and only a + // move the probes prefer pays for a re-integrated round. A move of a step or more is run as a + // round directly, as it always was, and so is every round of the walk it starts. + const auto on_lattice_pct = [](const ValidationSpotEvidence &e) { + return 100.0 * static_cast(e.on_lattice - e.by_chance) + / static_cast(std::max(1, e.spots)); + }; + // An indexing probe at a geometry: the validation evidence of the lattice indexed there, and + // nothing else - the pass stops once the lattice is scored, and the experiment is put back. + const auto index_at_geometry = [&](const std::array &g) { + const DiffractionExperiment before_probe = experiment_; + const bool probe_was = postrefine_probe_; + const bool searched_before_probe = beam_center_searched_; + beam_center_searched_ = true; // a probe scores the geometry it is given, nothing else + set_geometry(g); + indexing_probe_only_ = true; + postrefine_probe_ = false; + ProcessResult r; + try { + r = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Info("Two-pass: the indexing probe at distance {:.3f} mm did not complete ({})", + g[2], e.what()); + } + indexing_probe_only_ = false; + postrefine_probe_ = probe_was; + beam_center_searched_ = searched_before_probe; + experiment_ = before_probe; + ++arm_passes; + return r.validation_evidence; + }; constexpr int MAX_GEOMETRY_ROUNDS = 8; int geometry_rounds = 0; // passes run at a geometry the walk moved to int walk_passes = 0; // every pass the walk ran, the return to its best round included @@ -2537,11 +2585,14 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { std::string walk_rounds; // distance / realised residual per round, for the pass decision if (pass2.post_refine) { best_fit = *pass2.post_refine; - walk_rounds = fmt::format("{:.3f} mm / {:.3e}", best_geometry[2], best_fit.held_out_before); + walk_rounds = fmt::format("{:.3f} mm / {:.3e} / {:.1f}%", best_geometry[2], best_fit.held_out_before, + on_lattice_pct(pass2.validation_evidence)); } + ValidationSpotEvidence best_evidence = pass2.validation_evidence; + std::optional probe_at_best; // the indexing probe at best_geometry, once run + bool large_walk = false; // a round of this walk was started by a move of a step or more while (!cancelled_ && prepass_detector_geometry_ - && pass2.post_refine && pass2.post_refine->detector_refined - && (pass2.post_refine->large_move || geometry_rounds > 0)) { + && pass2.post_refine && pass2.post_refine->detector_refined) { if (geometry_rounds == MAX_GEOMETRY_ROUNDS) { geometry_not_converged = true; walk_stop = fmt::format("the fit still moved the geometry after {} rounds, every one of " @@ -2552,6 +2603,23 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { } const PostRefineResult fit = *pass2.post_refine; const std::array g = *prepass_detector_geometry_; + large_walk = large_walk || fit.large_move; + if (!large_walk) { + if (!probe_at_best) + probe_at_best = index_at_geometry(best_geometry); + const ValidationSpotEvidence probe_at_fit = index_at_geometry(g); + const bool prefers = ValidationEvidencePrefers(*probe_at_best, probe_at_fit); + logger.Info("Two-pass: the post-refinement moves the geometry by less than a step (distance " + "{:.3f} -> {:.3f} mm, beam {:.2f},{:.2f} -> {:.2f},{:.2f} px); indexing probes put " + "{:.1f}% of the validation spots on the lattice beyond chance there against {:.1f}% " + "at the geometry in hand - {}", fit.distance_before_mm, fit.distance_after_mm, + fit.beam_x_before_px, fit.beam_y_before_px, fit.beam_x_after_px, fit.beam_y_after_px, + on_lattice_pct(probe_at_fit), on_lattice_pct(*probe_at_best), + prefers ? "worth a round" : "the geometry in hand stands"); + if (!prefers) + break; + probe_at_best = probe_at_fit; + } logger.Info("Two-pass: the post-refinement at the adopted geometry still moves it " "(distance {:.3f} -> {:.3f} mm) - re-integrating and re-indexing at what it " "asks for (round {} of at most {})", fit.distance_before_mm, @@ -2561,20 +2629,26 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { ++geometry_rounds; ++walk_passes; const double realised = pass2.post_refine ? pass2.post_refine->held_out_before : NAN; - walk_rounds += fmt::format(", {:.3f} mm / {:.3e}", g[2], realised); - if (!pass2.post_refine || !HeldOutResidualFell(best_fit, *pass2.post_refine)) { + walk_rounds += fmt::format(", {:.3f} mm / {:.3e} / {:.1f}%", g[2], realised, + on_lattice_pct(pass2.validation_evidence)); + const bool residual_fell = pass2.post_refine && HeldOutResidualFell(best_fit, *pass2.post_refine); + const bool evidence_prefers = ValidationEvidencePrefers(best_evidence, pass2.validation_evidence); + if (!pass2.post_refine || !(residual_fell || evidence_prefers)) { walk_stop = fmt::format( - "round {} did not lower the realised held-out residual by more than its noise " - "({:.3e} against {:.3e}, standard errors {:.1e} and {:.1e})", geometry_rounds, - realised, best_fit.held_out_before, + "round {} lowered neither the realised held-out residual by more than its noise " + "({:.3e} against {:.3e}, standard errors {:.1e} and {:.1e}) nor the share of " + "validation spots off the lattice ({:.1f}% on it beyond chance against {:.1f}%)", + geometry_rounds, realised, best_fit.held_out_before, pass2.post_refine ? pass2.post_refine->held_out_before_se : NAN, - best_fit.held_out_before_se); + best_fit.held_out_before_se, on_lattice_pct(pass2.validation_evidence), + on_lattice_pct(best_evidence)); logger.Info("Two-pass: {}", walk_stop); break; } best_round = geometry_rounds; best_geometry = g; best_fit = *pass2.post_refine; + best_evidence = pass2.validation_evidence; const int pushback = ReindexPushesCellBack(fit, *pass2.post_refine); if (pushback != 0 && pushback == last_pushback) { walk_stop = fmt::format( @@ -2615,7 +2689,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { if (geometry_rounds > 0) pass2.pass_decision = fmt::format( "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm " - "against the header's {:.3f}: {}; distance / realised held-out residual by round: {}", + "against the header's {:.3f}: {}; distance / realised held-out residual / validation spots on the lattice beyond chance by round: {}", geometry_rounds, geometry_rounds == 1 ? "" : "s", best_round, best_geometry[2], pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0, walk_stop.empty() ? "the fit at the kept geometry commits no further move" : walk_stop, @@ -5312,6 +5386,35 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } logger.Info("Supercell probe (the eight parity classes of the doubled primitive cell; , " "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. + 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(); + int best = 0; + for (int c = 1; c < 8; c++) + if (best == 0 || FitSupercellProbe(band(c)).mean_i > FitSupercellProbe(band(best)).mean_i) + 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]); + 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(); + result.supercell = sc; + } } } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 29e1c315a..77dd4c715 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -362,6 +362,16 @@ struct ProcessResult { // Higher-order contamination of the beam, read off the spots the lattice did not take. Measured // on the images themselves, so it is there whether or not anything merged. HarmonicContaminationResult harmonic; + // 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. + 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{}; + }; + std::optional supercell; // Powder contamination measured in the pre-scan: a crystalline phase other than the crystal, // diffracting as rings among its reflections. Measured and reported on every run. PowderRings powder;