From fd80583a5c8cecf0cce87a12a564c7a8169b7fc3 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 9 Nov 2026 09:59:11 +0100 Subject: [PATCH] rotation: the pre-pass keeps the plain per-frame fit; the prior's widths only from frames that fitted Two defects in the restrained per-frame refinement, both in what the geometry pre-pass measures. The spindle form (XtalOptimizerOnSpindle) was applied to every integrated rotation frame, the pre-pass's included, although the pre-pass is meant to keep the free per-frame fit. The pre-pass takes the run's geometry decisions on per-frame fits - the beam-centre arm comparison, the metric- symmetry check's held-out residuals, the post-refinement - and the spindle form changes how much of a wrong metric a frame's orientation absorbs. The metric-symmetry check then lost its margin: a tetragonal-promoted orthorhombic crystal (a and b 0.8% apart) kept a = b and ended in P1121, and a monoclinic crystal was demoted to triclinic. With the spindle form switched off alone both decisions return to the unrestrained ones; with the prior switched off alone they do not. The spindle form now goes with the prior: only the frames refined under it use it, and the pre-pass is the unrestrained one up to the integration of the pass after it. A pre-pass frame whose fit failed (too few spots in tolerance) keeps the sweep's geometry and was recorded as a zero shift. On sparse runs that is half the frames or more, and it pulled the running median, and the prior with it, toward zero - on a run where one frame in sixty indexes, to the 1e-5 rad clamp. Only frames whose fit succeeded are recorded now. Checked against the unrestrained build and 1f3c11f61, each with the first or last ten frames dropped as well as on the whole sweep: the orthorhombic set is P212121 on all three frame sets (1f3c11f61: P1121 on the whole sweep); the two-crystal I222 CCD set stays at ISa 15-19 on all three (1f3c11f61: 8.1 / 18.1 / 8.1; unrestrained 13.4-14.4); the ice-ring monoclinic set is P2/P21 on all three (unrestrained: P21, P21, P1). The small-molecule gains stay (kdp ISa 26-33 against 20-24 unrestrained, laspartic and lcystine R_meas below the unrestrained build's on every frame set). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK (cherry picked from commit 20f51f1cf8f05dc95cd8f845803a1f24f04abc2e) --- docs/CPU_DATA_ANALYSIS_INDEXING.md | 2 +- image_analysis/IndexAndRefine.cpp | 17 ++++++++++++----- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index b24608a31..6564911aa 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -356,7 +356,7 @@ Refinement is performed in stages with decreasing acceptance tolerance for inclu The loose first stage necessarily admits some spots that are not reflections of this lattice — the fraction of *randomly* placed spots inside a fractional-Miller tolerance $t$ is $\tfrac{4}{3}\pi t^3$, i.e. 11 % at $t=0.3$ — and an unweighted fit lets them pull the orientation. Each residual is therefore weighted by how strong its spot is **for its resolution**: the frame's spots are cut into equal-count resolution shells and each intensity is divided by its shell median, mapped to $w^2=r/(1+r)$. The shell normalisation is what makes this safe — genuine high-resolution spots are legitimately weaker and carry the cell and distance information, so an un-normalised intensity weight would suppress exactly the spots the fit needs. The weight is a property of the spot and never of the current residual, so it does not depend on how far the geometry is from convergence. -**A rotation frame's refinement is restrained toward the sweep.** A frame is integrated on the orientation and beam centre refined against its own spots, starting from the sweep's geometry. Where the crystal really moves the refined orientation follows it, and integrating on it is what keeps such a sweep together. Where a frame has few or poorly placed spots the same fit can instead jump on its own, by up to a few degrees, which drops strong reflections out of prediction on the frame they peak on. So after the geometry pre-pass, which refines its frames freely, each integrated frame is fitted twice: freely, for the scale $\sigma^2$ of its own weighted misses, and then from the sweep again with a Gaussian prior on how far it moves, minimising $\sum_i w_i^2\|\mathbf{r}_i\|^2/\sigma^2 + \sum_k \theta_k^2/\tau_k^2$ over the orientation and beam-centre shifts $\theta_k$ (the orientation's as the three components of the rotation from the sweep's orientation, so that its width means the same in every direction however the crystal sits in the laboratory). Both fits of an integrated frame are made with the spindle given, as over a sweep, so the beam centre's component along the spindle - which on a single frame nearly trades against the orientation - carries the same weak restraint toward the start that it carries in the sweep refinement. The widths $\tau$ - one for the orientation, one for the beam centre - are measured once per run on the pre-pass frames, as the root mean square of the seven-frame running median of their free shifts (the track's smooth part, which a jump on one frame does not reach). There is no gate: the data of each frame decide how far it moves, a frame with many well placed spots following a moving crystal and a sparse one staying near the sweep. Each frame is still refined on its own spots, independently of the others. The frames the first pass scores candidate lattices on and the pre-pass keep the free refinement, so the lattice and the geometry the run refines do not depend on the prior, and stills are not affected. +**A rotation frame's refinement is restrained toward the sweep.** A frame is integrated on the orientation and beam centre refined against its own spots, starting from the sweep's geometry. Where the crystal really moves the refined orientation follows it, and integrating on it is what keeps such a sweep together. Where a frame has few or poorly placed spots the same fit can instead jump on its own, by up to a few degrees, which drops strong reflections out of prediction on the frame they peak on. So after the geometry pre-pass, which refines its frames freely, each integrated frame is fitted twice: freely, for the scale $\sigma^2$ of its own weighted misses, and then from the sweep again with a Gaussian prior on how far it moves, minimising $\sum_i w_i^2\|\mathbf{r}_i\|^2/\sigma^2 + \sum_k \theta_k^2/\tau_k^2$ over the orientation and beam-centre shifts $\theta_k$ (the orientation's as the three components of the rotation from the sweep's orientation, so that its width means the same in every direction however the crystal sits in the laboratory). Both fits of such a frame are made with the spindle given, as over a sweep, so the beam centre's component along the spindle - which on a single frame nearly trades against the orientation - carries the same weak restraint toward the start that it carries in the sweep refinement. The widths $\tau$ - one for the orientation, one for the beam centre - are measured once per run on the pre-pass frames, as the root mean square of the seven-frame running median of the free shifts of the frames whose fit succeeded (a frame too sparse to fit keeps the sweep's geometry, which is not a measurement of how far it moves) (the track's smooth part, which a jump on one frame does not reach). There is no gate: the data of each frame decide how far it moves, a frame with many well placed spots following a moving crystal and a sparse one staying near the sweep. Each frame is still refined on its own spots, independently of the others. The frames the first pass scores candidate lattices on and the pre-pass keep the free refinement of a single frame, without the spindle form as well as without the prior, so the lattice, the geometry the run refines and every decision the pre-pass takes do not depend on either, and stills are not affected. **The rotation chain commits its best round, not the round it stops on.** After the winning candidate is selected, the refinement is run again — solve, re-accumulate the reciprocal-space cloud under the refined geometry, solve again — up to twenty times, and the loop stops on a test of the detector-tilt step. The chain is a trajectory and its last point is not always its best one: every solve ends by fitting only the spots inside its **tightest** gate, so a cell with a direction the data barely constrain (which is what a free cell whose metric is near a Bravais class has) can slide along it, pulling a core of spots tighter while the periphery falls out of the fit altogether. Measured on such a chain, the spots inside the tight gate rise over seventeen rounds while the spots inside the widest gate peak at round three and fall away — and round three is the round that merges at ISa 11.0 against 6.3 and $R_\mathrm{meas}$ 0.148 against 0.213. Nothing the run consulted could see it: indexed fraction, validation frames, validation spots and the tight-gate count all prefer the overfitted end. diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 3aa7cea70..2c4cf73b0 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -490,6 +490,11 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In data.crystal_system = gemmi::CrystalSystem::Hexagonal; const bool rotation_integration = integrating && experiment.IsRotationIndexing(); + // Restrained: the integrated rotation frames of the passes after the geometry pre-pass, once the run + // has its prior (SetRefinementPrior). Every other frame - the pre-pass's included - keeps the plain + // per-frame fit, so the pre-pass measures the frames, and decides everything it decides (beam centre, + // metric symmetry, post-refinement, the prior's own widths), exactly as without the restraint. + const bool restrained_frame = rotation_integration && refinement_prior_.has_value(); const CrystalLattice start_lattice = data.latt; const float start_beam_x = data.geom.GetBeamX_pxl(), start_beam_y = data.geom.GetBeamY_pxl(); switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) { @@ -502,9 +507,9 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In break; case GeomRefinementAlgorithmEnum::BeamCenter: { const XtalOptimizerData start = data; - // An integrated rotation frame is refined about the spindle (XtalOptimizerOnSpindle), both fits. + // A restrained frame is refined about the spindle (XtalOptimizerOnSpindle), both fits. const auto fit = [&](XtalOptimizerData &d) { - return rotation_integration + return restrained_frame ? XtalOptimizerOnSpindle(d, msg.spots, msg.number, *experiment.GetGoniometer()) : XtalOptimizer(d, msg.spots); }; @@ -514,7 +519,7 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In // the frame's own spots at the residual scale the free fit just found. A frame with few or noisy // spots then stays near the sweep and a well determined one follows a crystal that moves // (empirical Bayes: the widths are the run's, the evidence each frame's own). - if (refined && rotation_integration && refinement_prior_) { + if (refined && restrained_frame) { const double sigma_res = std::sqrt(data.fit_residual_var); XtalOptimizerData restrained = start; restrained.prior_lattice = start.latt; @@ -570,8 +575,10 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In outcome.lattice_candidate = data.latt; // The cell is held on rotation, so the refined vectors are the sweep's turned by one rotation R, - // A_refined = R A_start for the matrix A of the three real-space vectors. - if (rotation_integration && collect_frame_refinement_) { + // A_refined = R A_start for the matrix A of the three real-space vectors. Only a frame whose fit + // succeeded measured anything: one that failed keeps the sweep's geometry, and counted as a zero + // shift it would pull the prior toward a clamp on a run where many frames are too sparse to fit. + if (rotation_integration && collect_frame_refinement_ && outcome.beam_center_updated) { const auto matrix = [](const CrystalLattice &l) { Eigen::Matrix3d m; const Coord v[3] = {l.Vec0(), l.Vec1(), l.Vec2()};