From 032bf9fe2bbf0e03f55b11a82dd0dcef2c63466a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 24 Sep 2026 11:34:06 +0200 Subject: [PATCH] Rugnux: walk the goniometer rotation scale to its fixed point, decided on the whole sweep The pass-1 post-refinement fits the rotation scale k only on the frames the stored angles still track, and a rate error is exactly what stops them tracking the rest: on a sweep whose stage turned ~3 % slow the fit read 0.979 over 94 deg, failed its leave-a-fifth-out test and was thrown away, leaving half the frames unscaled. The pass no longer decides. Between the passes, at the pass-2 detector geometry, the lattice is indexed (index-only probe) under the stored angles and under the fitted k and scored on the validation frames of the whole sweep: share of the spots on the lattice beyond the wrong-spindle null. k is adopted only where it scores higher by more than the binomial noise of the two (ValidationEvidencePrefers - the test the beam-centre arms already used, now one function); the run then integrates and post-refines at k (post-refine-only probe), fits again on top of it and repeats until the next k no longer scores better (WalkRotationScale). The stored angles are the first hypothesis. Measured: 1.000 25.1 %, 0.97874 58.9 %, 0.97041 90.0 %, 0.97006 90.7 % (not significant) -> 0.97041 adopted. Probes restore the experiment, the pass-2 geometry, pass-1 mosaicity and the beam-centre-search flag; a probe opens no beam-centre search. Forced pass-1 results get their axis scaled; the header revert drops the scale. GONIOMETER_ROTATION_SCALE reports the adopted k (SUSPECT = adopted). The leave-a-fifth-out figure stays in the log only. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- docs/CPU_DATA_ANALYSIS_INDEXING.md | 2 +- image_analysis/geom_refinement/PostRefine.cpp | 53 +---- image_analysis/geom_refinement/PostRefine.h | 4 +- rugnux/Rugnux.cpp | 212 +++++++++++++----- rugnux/Rugnux.h | 48 +++- tests/CMakeLists.txt | 1 + tests/RotationScaleWalkTest.cpp | 94 ++++++++ 7 files changed, 316 insertions(+), 98 deletions(-) create mode 100644 tests/RotationScaleWalkTest.cpp diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index 8f7e5cc7f..d219691ea 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -335,7 +335,7 @@ The space group is determined **after** pass 2, on the geometry the run refined, 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. -**Goniometer rotation scale (report only).** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. The excitation residual already measures it without a new degree of freedom: it rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under its own leave-a-fifth-of-the-sweep-out check — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it **under-reads** the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error. +**Goniometer rotation scale.** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. The excitation residual already measures it without a new degree of freedom: it rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor $k$ by which the stage actually turned, and normalising the axis throws it away. Pass 1 fits $k$ as a single parameter on its rocking events, with the crystal and the axis direction held at their committed values and the angle measured from the centre of the sweep. That fit **under-reads** a real error: it only sees the frames the stored angles still track, and a rate error is exactly what stops them tracking the rest. So it is not acted on directly. Between the passes, at the detector geometry pass 2 runs at, the lattice is indexed under the stored angles and under the fitted $k$, and each is scored on the validation frames spread over the whole sweep, as the share of their spots it puts on the lattice beyond what it puts there at a wrong spindle angle. The fitted $k$ is adopted only where it scores higher by more than the binomial noise of the two scores (z = 3.29); the run then integrates and post-refines at it, fits $k$ again on top of it, and repeats until the next $k$ no longer scores better - the fixed point of the fit. Otherwise the stored angles stand. The adopted $k$ drives every later pass (prediction, integration, scaling and the reported oscillation) and is reported as `GONIOMETER_ROTATION_SCALE`, with `GONIOMETER_ROTATION_SCALE_SUSPECT= TRUE`. `--rotation-scale ` asserts a calibration and skips all of this. ### 7.6 Detector geometry from powder rings diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index a2fc800f0..4acd73472 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -1080,50 +1080,19 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations, const double k_fit = solve_scale(-1); result.rotation_scale = k_fit; - // ----- Whether to COMMIT it. A stage fault is rare - 36 of 37 rotation datasets sit at 1.0000 - // on a direct scan - and a 1 % angle correction applied to a healthy dataset would damage it - // silently, so every test below has to pass. - // Preconditions: below these the fit is reported but never acted on. Under ~30 deg of sweep k - // entangles with the axis direction and 10-20 deg truncations of a perfect dataset wander by - // +-0.6 %; a screening wedge must not trigger a correction. - constexpr int MIN_SCALE_EVENTS = 5000; - constexpr double MIN_SCALE_SWEEP_DEG = 30.0; - // T1 significance: 0.5 % is 18 sigma on the between-dataset scatter of healthy stages - // (robust sd 2.8e-4) and still 3.5x below the one measured fault. - constexpr double ROTATION_SCALE_TOL = 0.005; - // T2 relevance: the misorientation the error produces at each end of the sweep. A large k over - // a short sweep moves nothing and is not worth correcting. - constexpr double MIN_SCALE_END_ERROR_DEG = 0.5; - // T3 uniformity: a stage error is a ramp present in EVERY part of the sweep, so dropping any - // fifth of it must leave the same k. A second lattice that dominates ONE END of the sweep - - // exactly what happens where the primary stops indexing - fakes a k indistinguishable from a - // real fault on T1 and T2, and is the reason this test is not optional. It replaces the - // hkl-hash split used elsewhere here, which cannot see it: both halves of that split sit at - // the same angles, so anything structured in phi survives in both folds. - constexpr double MIN_SCALE_JACKKNIFE_FRAC = 0.5; - const double end_error_deg = std::fabs(k_fit - 1.0) * sweep_deg / 2.0; - const bool enough_data = static_cast(n_events) >= MIN_SCALE_EVENTS - && sweep_deg >= MIN_SCALE_SWEEP_DEG; - const bool big_enough = enough_data && std::fabs(k_fit - 1.0) >= ROTATION_SCALE_TOL - && end_error_deg >= MIN_SCALE_END_ERROR_DEG; + // Whether the same k comes back with any fifth of the sweep left out, as the smallest share + // of the fitted excess the folds keep: a stage error is a ramp present in EVERY part of the + // sweep. Reported, not acted on - the fit only sees the frames the angles it was measured at + // still track, and a rate error is exactly what stops them tracking the rest, so the part + // it sees can be too short to agree with itself. What acts on k is the caller, which walks + // it to its fixed point and decides it on the whole sweep (rugnux WalkRotationScale). double jackknife = 1.0; - if (big_enough) - for (int f = 0; f < 5; ++f) - jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0)); - result.rotation_scale_suspect = big_enough && jackknife >= MIN_SCALE_JACKKNIFE_FRAC; + for (int f = 0; f < 5; ++f) + jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0)); logger.Info("Post-refine rotation SCALE: k = {:.5f} over {:.0f} deg of sweep centred on {:.1f} " - "deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f} => {}", - k_fit, sweep_deg, phi_c * 180.0 / PI, n_events, end_error_deg, jackknife, - result.rotation_scale_suspect ? "COMMIT" - : !enough_data ? "report only (too little sweep or too few events)" - : "reject (kept the stored angles)"); - if (result.rotation_scale_suspect) - logger.Warning("Goniometer rotation scale looks off by {:+.2f} % (fitted {:.5f}): the stage " - "appears to have turned {} than the angles stored in the file, which are the " - "COMMANDED values. This is a hardware calibration fault, not a data problem - " - "left uncorrected it inflates mosaicity, biases the cell and loses " - "high-resolution reflections", - 100.0 * (k_fit - 1.0), k_fit, k_fit > 1.0 ? "further" : "less far"); + "deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f}", + k_fit, sweep_deg, phi_c * 180.0 / PI, n_events, + std::fabs(k_fit - 1.0) * sweep_deg / 2.0, jackknife); // Assemble the committed geometry. result.distance_after_mm = dist[0]; diff --git a/image_analysis/geom_refinement/PostRefine.h b/image_analysis/geom_refinement/PostRefine.h index b5dcf079b..6ed478587 100644 --- a/image_analysis/geom_refinement/PostRefine.h +++ b/image_analysis/geom_refinement/PostRefine.h @@ -73,8 +73,8 @@ struct PostRefineResult { // the header). Fitted after the joint fit as a single free parameter, with the crystal and the axis // direction held at their committed values. Always the fitted value; 1.0 = header and stage agree. double rotation_scale = 1.0; - // Whether the fit passed every test needed to ACT on it: enough sweep and events, a significant and - // physically relevant size, and the same k from every fifth of the sweep. Only then is it applied. + // Whether the run ACTED on it: set by the caller where the scale, walked to its fixed point, was + // adopted on the evidence of the whole sweep (rugnux WalkRotationScale) - never by the fit itself. bool rotation_scale_suspect = false; }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6c5fdcc3f..05d9296ec 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -75,6 +75,48 @@ #include #include +bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate) { + const auto excess = [](const ValidationSpotEvidence &e) { + return static_cast(e.on_lattice - e.by_chance) / static_cast(std::max(1, e.spots)); + }; + const auto rate_var = [](const ValidationSpotEvidence &e) { + const double n = static_cast(std::max(1, e.spots)); + const double p = static_cast(e.on_lattice) / n; + return p * (1.0 - p) / n; + }; + return excess(candidate) - excess(current) + > SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(rate_var(candidate) + rate_var(current)); +} + +RotationScaleWalk WalkRotationScale(double first_fit, + const std::function &index_at, + const std::function(float)> &refit_at, + int max_rounds) { + RotationScaleWalk walk; + auto next = static_cast(first_fit); + if (next == walk.scale) + return walk; + const auto score = [](const ValidationSpotEvidence &e) { + return 100.0 * static_cast(e.on_lattice - e.by_chance) + / static_cast(std::max(1, e.spots)); + }; + walk.evidence = index_at(walk.scale); + walk.trail = fmt::format("{:.5f}: {:.1f}%", walk.scale, score(walk.evidence)); + for (int round = 0; round < max_rounds && next != walk.scale; ++round) { + const ValidationSpotEvidence e = index_at(next); + walk.trail += fmt::format(", {:.5f}: {:.1f}%", next, score(e)); + if (!ValidationEvidencePrefers(walk.evidence, e)) + break; + walk.scale = next; + walk.evidence = e; + const auto fit = refit_at(walk.scale); + if (!fit) + break; + next = static_cast(walk.scale * *fit); + } + return walk; +} + double MetricViolation(const UnitCell &uc, const gemmi::SpaceGroup &sg) { const double a = uc.a, b = uc.b, c = uc.c; const double d2r = PI / 180.0; @@ -2024,7 +2066,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { const DiffractionExperiment file_experiment = experiment_; const auto file_mosaicity = prepass_mosaicity_; const auto file_geometry = prepass_detector_geometry_; - const auto file_scale = prepass_rotation_scale_; const auto file_result = prepass_result_; experiment_ = experiment_before_first_pass; @@ -2036,7 +2077,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { experiment_.BeamX_pxl(alt_center[0]).BeamY_pxl(alt_center[1]); prepass_mosaicity_.clear(); prepass_detector_geometry_.reset(); - prepass_rotation_scale_.reset(); prepass_result_.reset(); ProcessResult alt; @@ -2113,7 +2153,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { experiment_ = file_experiment; prepass_mosaicity_ = file_mosaicity; prepass_detector_geometry_ = file_geometry; - prepass_rotation_scale_ = file_scale; prepass_result_ = file_result; } // The arm's own pass asks the same question again at its own centre; the run has already @@ -2152,10 +2191,90 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]) .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); } - // ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles - // in the file are the commanded ones, so a stage that ran fast is a geometry error like any other. + // The probe passes below and the lattice arms after them exist ONLY to measure - see there - so + // the count of passes the run made has to include them even though they write nothing. + int arm_passes = 0; + + // The goniometer rotation scale: the angles in the file are the commanded ones, so a stage that + // turned at the wrong rate is a geometry error like any other - but one pass 1's fit reads only + // in part (see WalkRotationScale), so it is walked to the fit's fixed point here and adopted + // only where the validation frames of the whole sweep prefer it. The walk runs at the detector + // geometry the second pass will, so every hypothesis, the stored angles included, is scored + // there, and at that geometry only: a probe whose angles lose the lattice must score what it + // scores, not start the beam-centre search that a poor first pass otherwise opens. Its passes + // only measure - an indexing probe stops once the lattice is scored, a refit stops once the + // post-refinement has measured - and they leave nothing behind: the experiment, the geometry + // the second pass is to run at, pass 1's mosaicity (a width in degrees fitted against the + // stored angles) and whether the run has searched for the beam centre are put back when the + // walk ends. + std::string rotation_scale_walk; + if (!cancelled_ && gonio_snapshot && pass1.post_refine && !config_.rotation_scale) { + const DiffractionExperiment before_walk = experiment_; + const auto geometry_before_walk = prepass_detector_geometry_; + const auto mosaicity_before_walk = prepass_mosaicity_; + const bool searched_before_walk = beam_center_searched_; + beam_center_searched_ = true; + const auto probe = [&](float k, bool index_only) { + experiment_ = before_walk; + experiment_.Goniometer(ScaleRotation(*before_walk.GetGoniometer(), k)); + prepass_rotation_scale_ = k; + prepass_mosaicity_.clear(); + indexing_probe_only_ = index_only; + postrefine_probe_ = postrefine_probe_only_ = !index_only; + ProcessResult r; + try { + r = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + // Angles under which nothing indexes have scored nothing, which is what an empty + // result reads as. + logger.Info("Rotation scale {:.5f}: the probe pass did not complete ({})", k, e.what()); + } + indexing_probe_only_ = postrefine_probe_ = postrefine_probe_only_ = false; + ++arm_passes; + return r; + }; + const auto index_at = [&](float k) { return probe(k, true).validation_evidence; }; + const auto refit_at = [&](float k) -> std::optional { + const ProcessResult r = probe(k, false); + if (!r.post_refine) + return std::nullopt; + logger.Info("Rotation scale {:.5f}: the post-refinement there fits {:.5f} on top of it, " + "held-out residual {:.3e}", k, r.post_refine->rotation_scale, + r.post_refine->held_out_before); + return r.post_refine->rotation_scale; + }; + constexpr int MAX_ROTATION_SCALE_ROUNDS = 8; + const RotationScaleWalk walk = WalkRotationScale(pass1.post_refine->rotation_scale, index_at, + refit_at, MAX_ROTATION_SCALE_ROUNDS); + experiment_ = before_walk; + prepass_detector_geometry_ = geometry_before_walk; + prepass_mosaicity_ = mosaicity_before_walk; + beam_center_searched_ = searched_before_walk; + prepass_rotation_scale_.reset(); + if (!walk.trail.empty()) { + rotation_scale_walk = fmt::format( + "goniometer rotation scale walked from the pass-1 fit {:.5f}, validation spots on " + "the lattice beyond chance by scale: {} - {}", pass1.post_refine->rotation_scale, + walk.trail, walk.scale == 1.0f ? "the stored angles stand" + : fmt::format("{:.5f} adopted", walk.scale)); + logger.Info("Two-pass: {}", rotation_scale_walk); + } + if (walk.scale != 1.0f) { + prepass_rotation_scale_ = walk.scale; + pass1.post_refine->rotation_scale = walk.scale; + pass1.post_refine->rotation_scale_suspect = true; + logger.Warning("Goniometer rotation scale {:.5f} ({:+.2f} %): the stage turned {} than the " + "angles stored in the file, which are the COMMANDED values. The second pass " + "integrates at the corrected angles; the fault is in the hardware and should " + "be fixed there", walk.scale, 100.0 * (walk.scale - 1.0), + walk.scale > 1.0f ? "further" : "less far"); + } + } + + // ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry. if (prepass_rotation_scale_ && gonio_snapshot) { - experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_)); + experiment_.Goniometer(ScaleRotation(*experiment_.GetGoniometer(), *prepass_rotation_scale_)); // The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has // just stopped using, and the override can only ever raise the second pass's own estimate (it // takes the larger of the two). Carrying it over would hold the second pass at the rocking @@ -2197,10 +2316,6 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]); }; - // The lattice arms below run passes that exist ONLY to measure - see there - so the - // count of passes the run made has to include them even though they write nothing. - int arm_passes = 0; - // The metric symmetry, asked of the spots rather than of the tolerance that admitted it. The // Bravais class is chosen by a walk that reads two axes as EQUAL when they agree to a fixed // relative tolerance, and the class carrying that equality is then imposed everywhere below: @@ -2548,6 +2663,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { if (const auto current = experiment_.GetGoniometer()) restored.Axis(current->GetAxis()); experiment_.Goniometer(restored); + prepass_rotation_scale_.reset(); + pass1.post_refine->rotation_scale_suspect = false; } config_.output_prefix = base_prefix; // This pass is the answer whatever it measures - the guard has had its one chance - @@ -2573,6 +2690,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { : "the post-refinement committed no geometry change, so this pass reproduces the first"; if (!lattice_arm.empty()) pass2.pass_decision = lattice_arm + "; " + pass2.pass_decision; + if (!rotation_scale_walk.empty()) + pass2.pass_decision = rotation_scale_walk + "; " + pass2.pass_decision; pass2.geometry_not_converged = geometry_not_converged; return pass2; } @@ -3234,11 +3353,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The pooled evidence for one candidate lattice: the validation frames' non-ice spots, how many // of them lie on the candidate, and how many a wrong spindle angle still puts there. - struct PooledEvidence { - int64_t spots = 0; - int64_t on_lattice = 0; - int64_t by_chance = 0; - }; + using PooledEvidence = ValidationSpotEvidence; // Displacements are a fixed fraction of the sweep, so the null is a function of the data alone // and the same file gives the same verdict every time. Sevenths: no crystallographic rotation @@ -3274,19 +3389,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b return static_cast(e.on_lattice - e.by_chance) > SPOT_BUDGET_SIGNIFICANCE_Z * sigma; }; - // The same quantity for COMPARING two lattices: the share of the pooled spots each one puts on - // itself over and above what a wrong spindle angle puts there. Each is scored against its own - // null, so a denser lattice is not credited for the spots it catches by accident - which is - // what makes the comparison fair between a cell and its axis harmonic. - auto excess = [](const PooledEvidence &e) { - return static_cast(e.on_lattice - e.by_chance) - / static_cast(std::max(1, e.spots)); - }; - auto rate_var = [](const PooledEvidence &e) { - const double n = static_cast(std::max(1, e.spots)); - const double p = static_cast(e.on_lattice) / n; - return p * (1.0 - p) / n; - }; + // Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the + // pooled spots each one puts on itself over and above what a wrong spindle angle puts there. + // Each is scored against its own null, which is what makes the comparison fair between a cell + // and its axis harmonic. // How deep into an image's intensity-ordered spot list this lattice is still being seen - the // measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the @@ -4033,10 +4139,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b restore_beam_center(header_x, header_y); if (best.result.has_value()) at_header = pooled_evidence(*indexer, *best.result); - const double margin = excess(at_measured) - excess(at_header); - const double margin_sigma = std::sqrt(rate_var(at_measured) + rate_var(at_header)); if (alt.result.has_value() && beats_chance(at_measured) - && margin > SPOT_BUDGET_SIGNIFICANCE_Z * margin_sigma) { + && ValidationEvidencePrefers(at_header, at_measured)) { try_beam_center(measured_x, measured_y); logger.Warning("Beam centre from the background: neither centre indexes a " "validation frame on its own, but the pooled spots put " @@ -4179,11 +4283,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const PooledEvidence at_header = pooled_evidence(*indexer, *best.result); try_beam_center(measured_x, measured_y); const PooledEvidence at_measured = pooled_evidence(*indexer, *alt.result); - const double margin = excess(at_measured) - excess(at_header); - const double margin_sigma = std::sqrt(rate_var(at_measured) - + rate_var(at_header)); adopted_measured = beats_chance(at_measured) - && margin > SPOT_BUDGET_SIGNIFICANCE_Z * margin_sigma; + && ValidationEvidencePrefers(at_header, at_measured); if (adopted_measured) { logger.Warning("Beam centre check: the larger cell is the measured " "centre's, and it is the one the spots are on - " @@ -4639,6 +4740,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b gemmi::crystal_system_str(prepass_result_->search_result.system), pc.a, pc.b, pc.c, pc.alpha, pc.beta, pc.gamma); best.result = *prepass_result_; + // Pass 1's result carries pass 1's goniometer; forcing it whole would put the + // uncorrected angles back (the same repair as the supercell re-run in RunAllPasses). + if (prepass_rotation_scale_ && best.result->axis) + best.result->axis = ScaleRotation(*best.result->axis, *prepass_rotation_scale_); // Pass 1's lattice was FITTED at pass 1's detector distance, and this pass // integrates at the post-refined one. Re-scoring it here is not the same as // re-fitting it: a real-space cell is measured against the distance the spots were @@ -4685,6 +4790,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b / static_cast(std::max(1, evidence.spots)); const double pooled_chance = static_cast(evidence.by_chance) / static_cast(std::max(1, evidence.spots)); + result.validation_evidence = evidence; + // A pass run only to score a lattice has its score, whatever it is: a lattice that does not + // beat chance is a result for the comparison that asked, not a reason to stop the run. + if (indexing_probe_only_) { + logger.Info("Indexing probe: scheme '{}', {}/{} validation frames, {}/{} validation spots " + "= {:.1f}% against {:.1f}% at a wrong spindle angle", best.name, best.score, + static_cast(validation.size()), evidence.on_lattice, evidence.spots, + 100.0 * pooled, 100.0 * pooled_chance); + return result; + } if (!beats_chance(evidence)) { { // Name the cell and Bravais class that was rejected. The commonest cause is a metric @@ -4821,7 +4936,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b "pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead", best.result->search_result.centering, v2, prepass_result_->search_result.centering, v1); - indexer->ForceRotationIndexerResult(*prepass_result_); + RotationIndexerResult forced = *prepass_result_; + if (prepass_rotation_scale_ && forced.axis) + forced.axis = ScaleRotation(*forced.axis, *prepass_rotation_scale_); + indexer->ForceRotationIndexerResult(forced); } } } @@ -5529,19 +5647,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Two-pass: geometry post-refine committed no detector change{}", geometry_prepass ? " - the second pass reproduces the first" : ""); } - // DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the - // scale is its own cross-validated fit and the crystals that have a stage fault are - // exactly the ones whose cell and detector steps do NOT pass, because the angle error - // is what their residual is made of. A calibration fault is rare (36 of 37 rotation - // datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset - // would silently damage it, so the asymmetry is deliberate: committed only when the - // fit is both cross-validated and outside the tolerance. A manual --rotation-scale is - // already on the goniometer and is left alone. - // Only the pre-pass fits the rotation scale. It is applied from the goniometer - // the file came with, so a scale measured again on angles that have already been - // corrected once is not a correction that can be applied on top of that one. - if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value() && geometry_prepass) - prepass_rotation_scale_ = static_cast(pr.rotation_scale); + // The goniometer rotation scale in pr is only a fit: RunAllPasses walks it to its + // fixed point and decides it on the validation frames of the whole sweep + // (WalkRotationScale). It does not ride on pr.ok - the crystals that have a stage + // fault are exactly the ones whose cell and detector steps do NOT pass, because the + // angle error is what their residual is made of. } } }; @@ -5559,8 +5669,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (postrefine_probe_only_) { result.processing_time_s = std::chrono::duration( std::chrono::steady_clock::now() - start_time).count(); - logger.Info("Lattice-arm probe: {} images integrated and post-refined in {:.2f} s " - "(the arm reads the held-out residual only, so this pass does not merge)", + logger.Info("Probe pass: {} images integrated and post-refined in {:.2f} s " + "(the comparison reads the post-refinement only, so this pass does not merge)", result.images_processed, result.processing_time_s); return result; } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index c74278625..d6375cd58 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -186,6 +187,42 @@ struct ProcessConfig { bool finalist_ledger = false; // --finalist-ledger; report-only symmetry evidence table }; +// A rotation lattice's score on the validation frames spread over the whole sweep: their spots (off +// the ice rings unless those are indexed), how many lie on the lattice, and how many still do at a +// wrong spindle angle - the median over displaced angles, which is what chance and the per-frame +// orientation polish give for free. +struct ValidationSpotEvidence { + int64_t spots = 0; + int64_t on_lattice = 0; + int64_t by_chance = 0; +}; + +// Whether `candidate` puts a larger share of its validation spots on its lattice, over and above what +// a wrong spindle angle puts there, than `current` does - by more than the binomial noise of the two +// shares (SPOT_BUDGET_SIGNIFICANCE_Z). Each is measured against its own null, so a denser lattice is +// not credited for the spots it catches by accident. +bool ValidationEvidencePrefers(const ValidationSpotEvidence ¤t, const ValidationSpotEvidence &candidate); + +// The goniometer rotation scale - the factor by which the stage turned relative to the angles stored +// in the file - walked to the fit's fixed point and decided on the whole sweep. The fit only sees the +// frames the angles it was measured at still track, and a stage at the wrong rate is exactly what +// stops them tracking the rest, so one fit reads only part of the error. So: index the lattice under +// the fitted scale and under the angles in hand; where the fitted one scores better on the validation +// frames (ValidationEvidencePrefers), adopt it, fit again there and repeat. The stored angles (k = 1) +// are the first hypothesis, and stand unless the evidence moves the run off them. +// index_at(k): the validation evidence of the lattice indexed with the stored angles scaled by k +// refit_at(k): the scale the post-refinement fits on reflections integrated at k, relative to k; +// empty where it fitted none +struct RotationScaleWalk { + float scale = 1.0f; // adopted; 1 = the stored angles stand + ValidationSpotEvidence evidence; // at the adopted scale + std::string trail; // every scale tried, with its validation score +}; +RotationScaleWalk WalkRotationScale(double first_fit, + const std::function &index_at, + const std::function(float)> &refit_at, + int max_rounds); + struct ProcessResult { bool cancelled = false; uint64_t images_processed = 0; @@ -233,6 +270,8 @@ struct ProcessResult { // cell - so two cells can only be compared for volume once this has brought both to primitive. std::optional consensus_centering; bool rotation_lattice_found = false; + // The rotation lattice's score on the validation frames; all zero where no lattice was scored. + ValidationSpotEvidence validation_evidence; // The metric symmetry the rotation lattice was classified in - the class GeometryRefiner holds the // cell to and the space-group search enumerates under. Empty when no rotation lattice was found. std::optional rotation_lattice_type; @@ -484,8 +523,9 @@ class Rugnux { // describes neither geometry. The two travel together or not at all. std::optional> prepass_detector_geometry_; - // Two-pass geometry pre-pass: the goniometer rotation scale the post-refine fitted and flagged as a - // stage fault, applied by Run() to the second pass's goniometer. Empty when the fit found nothing. + // The goniometer rotation scale the run adopted (WalkRotationScale, in RunAllPasses), applied to the + // second pass's goniometer - and, while the walk runs, the scale its current probe pass is at. + // Empty where the stored angles stand. std::optional prepass_rotation_scale_; // Two-pass geometry pre-pass: pass-1's FULL indexing result (correct lattice + orientation + refined @@ -591,6 +631,10 @@ class Rugnux { // Everything below that point - the merges, the space-group search, the correction surfaces, the // reports - is work no comparison looks at, and on a long sweep it is four fifths of the pass. bool postrefine_probe_only_ = false; + // Whether a pass exists only to index: it stops as soon as its first-pass indexing has scored the + // lattice on the validation frames (ProcessResult::validation_evidence), before a single image is + // integrated. The rotation-scale walk (see RunAllPasses) compares angle models on that score. + bool indexing_probe_only_ = false; // The file's detector distance, from before the first pass of the rotation two-pass: a pass whose // distance is still this one asks the post-refinement to test "the header distance is right" // (PostRefineSettings::distance_at_header); a pass that has walked off it does not. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 64cd2092d..24f86cd3e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -59,6 +59,7 @@ ADD_EXECUTABLE(jfjoch_test ResultReportTest.cpp DiagnosticOutputTest.cpp PostRefineTest.cpp + RotationScaleWalkTest.cpp RugnuxLargeTest.cpp TestData.h MovingAverageTest.cpp diff --git a/tests/RotationScaleWalkTest.cpp b/tests/RotationScaleWalkTest.cpp new file mode 100644 index 000000000..562d46a24 --- /dev/null +++ b/tests/RotationScaleWalkTest.cpp @@ -0,0 +1,94 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "../rugnux/Rugnux.h" + +namespace { + // A synthetic sweep whose stage turned `true_scale` times the stored angles. Scored at a scale k, + // the validation spots stay on the lattice as long as the angles track the rotation, and the share + // that does falls off with the relative rate error; a wrong spindle angle keeps half a percent. + struct SyntheticSweep { + double true_scale; + int64_t spots = 35000; + int index_calls = 0; + int refit_calls = 0; + + ValidationSpotEvidence IndexAt(float k) { + ++index_calls; + const double error = std::fabs(true_scale / k - 1.0); + const double on = 0.9 * std::max(0.0, 1.0 - 30.0 * error); + return ValidationSpotEvidence{spots, std::llround(on * spots), std::llround(0.005 * spots)}; + } + + // The post-refinement at k, relative to k. It reads only 70 % of the error that is left: it + // sees only the frames the angles at k still track. + std::optional RefitAt(float k) { + ++refit_calls; + return 1.0 + 0.7 * (true_scale / k - 1.0); + } + + RotationScaleWalk Walk(double first_fit) { + return WalkRotationScale(first_fit, [this](float k) { return IndexAt(k); }, + [this](float k) { return RefitAt(k); }, 8); + } + }; +} + +TEST_CASE("ValidationEvidencePrefers", "[RotationScale]") { + const ValidationSpotEvidence base{10000, 3000, 50}; + // 1 % more of the spots beyond chance is under the noise of two 30 % shares over 10000 spots + // (sqrt(2 * 0.3 * 0.7 / 10000) = 0.65 %, times 3.29); 5 % is well over it. + CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3100, 50})); + CHECK(ValidationEvidencePrefers(base, {10000, 3500, 50})); + // A candidate is judged against its own null: more spots on the lattice bought by a null that + // rose just as much is no gain. + CHECK_FALSE(ValidationEvidencePrefers(base, {10000, 3500, 550})); + // Never against itself, and nothing that scored nothing wins. + CHECK_FALSE(ValidationEvidencePrefers(base, base)); + CHECK_FALSE(ValidationEvidencePrefers(base, {})); + CHECK(ValidationEvidencePrefers({}, base)); +} + +TEST_CASE("WalkRotationScale_ReachesTheFixedPoint", "[RotationScale]") { + // A stage 3 % slow. The first fit reads 70 % of that; each refit at the adopted scale reads 70 % + // of what is left, and the walk goes on as long as the validation frames prefer the new scale. + SyntheticSweep sweep{0.97}; + const auto walk = sweep.Walk(sweep.RefitAt(1.0f).value()); + CHECK(walk.scale == Catch::Approx(0.97).margin(0.001)); + CHECK(walk.scale != 1.0f); + CHECK(walk.evidence.on_lattice > sweep.IndexAt(1.0f).on_lattice); + CHECK(sweep.refit_calls > 2); + CHECK_FALSE(walk.trail.empty()); +} + +TEST_CASE("WalkRotationScale_StoredAnglesStand", "[RotationScale]") { + SECTION("A healthy stage: a fit off by noise scores no better than the stored angles") { + SyntheticSweep sweep{1.0}; + const auto walk = sweep.Walk(1.0002); + CHECK(walk.scale == 1.0f); + CHECK(sweep.index_calls == 2); // the stored angles and the fit, nothing more + CHECK(sweep.refit_calls == 0); + } + SECTION("A real but small error the spots cannot resolve beyond their noise") { + SyntheticSweep sweep{0.999}; + sweep.spots = 400; + const auto walk = sweep.Walk(0.9993); + CHECK(walk.scale == 1.0f); + } + SECTION("A fit that tracks something other than the rotation scores worse, and is refused") { + SyntheticSweep sweep{1.0}; + const auto walk = sweep.Walk(0.98); + CHECK(walk.scale == 1.0f); + CHECK(sweep.refit_calls == 0); + } + SECTION("A fit of exactly one asks for no probe at all") { + SyntheticSweep sweep{1.0}; + const auto walk = sweep.Walk(1.0); + CHECK(walk.scale == 1.0f); + CHECK(walk.trail.empty()); + CHECK(sweep.index_calls == 0); + } +}