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()};