From ef0be2e52ab909463891887eeba468a8716f62fc Mon Sep 17 00:00:00 2001 From: jungfrau Date: Sun, 23 Aug 2026 14:03:59 -0400 Subject: [PATCH] Bound the offline lattice refinement by iterations, and stop rotating the same axis three times Two things in the indexing path, one of them a reproducibility hole. XtalOptimizerData bounds a solve by iterations when it is told to and by WALL-CLOCK SECONDS when it is not, and its own header says why that matters: the same image refines to a different answer on a busier machine. The per-image refinement sets the iteration bound for exactly that reason. The rotation indexer never did, so its candidate-cell refinement ran under a one-second wall clock - three stages a candidate, up to eight candidates a scheme, twice a run. A run that has just been made reproducible from its prediction order to its accumulators was still free to pick a different lattice because the machine was loaded. It now takes the iteration bound offline and keeps the wall-clock one for a live acquisition, whose budget is real, which is the same split the per-image path already makes. The residual itself rotated the same axis three times over. It applies one orientation to three reciprocal-lattice vectors, and ceres::AngleAxisRotatePoint recomputes the angle, its sine, its cosine and the normalised axis on each call - and it does not inline at this optimisation level, so the compiler cannot notice. On a seventeen-parameter Jet each of those is a full dual-number evaluation. Computing the rotation once and applying it three times removes two hypots, two sines, two cosines, two divisions and six multiplies from every evaluation, which is about half the libm calls in it; hoisting a constant member's sine and cosine out of the same function takes two more. It runs everywhere the residual does - the indexer, the per-image refinement and the geometry refiner. Also lifts five SetParameterBlockConstant calls out of a per-observation loop in the detector solve, where they were executing once per observation to say the same thing. The rotation hoist was checked against the function it replaces on 200000 random dual numbers, including the small-angle branch, comparing the value and all seventeen derivative lanes: no difference in any component. Merged output is byte-identical on a 16 Mpx set and an ordinary one. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011n8riB6X59oRjkrSHzNPAU --- image_analysis/IndexAndRefine.cpp | 2 +- image_analysis/IndexAndRefine.h | 8 +- .../geom_refinement/XtalOptimizer.cpp | 7 +- image_analysis/geom_refinement/XtalResidual.h | 76 +++++++++++++++++-- .../rotation_indexer/RotationIndexer.cpp | 12 ++- .../rotation_indexer/RotationIndexer.h | 7 +- 6 files changed, 94 insertions(+), 18 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index c6455113..d5a1034b 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -71,7 +71,7 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool indexer_(indexer), rotation_indexer_counter(x) { if (indexer && x.IsRotationIndexing()) - rotation_indexer = std::make_unique(x, *indexer); + rotation_indexer = std::make_unique(x, *indexer, real_time); // Only retain the whole-run per-image reflections when a later scaling/merge pass will read them. if (retain_outcomes_) integration_outcome.resize(x.GetImageNum()); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index d79c2b0a..d10e66c7 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -92,10 +92,10 @@ class IndexAndRefine { std::optional RotationAngle(int64_t image) const; // mid-exposure angle for the indexer public: - // real_time: bound the per-image geometry refinement by WALL CLOCK, as online acquisition must - - // it has a real per-image budget. Offline (rugnux, the viewer) passes false and the refinement is - // bounded by iteration count instead, so the same file reprocesses to the same answer regardless of - // what else the machine was doing. + // real_time: bound the geometry refinements - the per-image one here and the candidate-cell ones + // in the rotation indexer - by WALL CLOCK, as online acquisition must, it having a real budget. + // Offline (rugnux, the viewer) passes false and they are bounded by iteration count instead, so the + // same file reprocesses to the same answer regardless of what else the machine was doing. IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer, bool retain_outcomes = true, bool real_time = false); diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 348cc9b0..1bd2a7a7 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -44,9 +44,10 @@ struct XtalResidualRotationOnlyPrecomp { T astar_rot[3], bstar_rot[3], cstar_rot[3]; - ceres::AngleAxisRotatePoint(rot_aa, astar_unrot, astar_rot); - ceres::AngleAxisRotatePoint(rot_aa, bstar_unrot, bstar_rot); - ceres::AngleAxisRotatePoint(rot_aa, cstar_unrot, cstar_rot); + const AngleAxisRotator rot(rot_aa); + rot.Rotate(astar_unrot, astar_rot); + rot.Rotate(bstar_unrot, bstar_rot); + rot.Rotate(cstar_unrot, cstar_rot); const Eigen::Matrix s_pred(T(h) * astar_rot[0] + T(k) * bstar_rot[0] + T(l) * cstar_rot[0], T(h) * astar_rot[1] + T(k) * bstar_rot[1] + T(l) * cstar_rot[1], diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index cf6a5a4d..93f8dcc3 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -13,6 +13,62 @@ #include "../../common/JFJochException.h" +// Rodrigues rotation with everything that depends only on the angle-axis worked out once. This is +// ceres::AngleAxisRotatePoint term for term - so the rotated point is bit-identical to it - but that +// function derives theta and its sine and cosine from the angle-axis on every call, which is pure +// repetition when one rotation is applied to several points, as it is to the three lattice columns +// below. On a Jet those trigonometric evaluations are the expensive part of the rotation. +template +struct AngleAxisRotator { + explicit AngleAxisRotator(const T *const angle_axis) + : aa{angle_axis[0], angle_axis[1], angle_axis[2]} { + using std::cos; + using std::fpclassify; + using std::hypot; + using std::sin; + + const T theta = hypot(aa[0], aa[1], aa[2]); + at_zero = (fpclassify(theta) == FP_ZERO); + if (at_zero) + return; + + costheta = cos(theta); + sintheta = sin(theta); + const T theta_inverse = T(1.0) / theta; + w[0] = aa[0] * theta_inverse; + w[1] = aa[1] * theta_inverse; + w[2] = aa[2] * theta_inverse; + } + + void Rotate(const T pt[3], T result[3]) const { + if (at_zero) { + // At zero the rotation is I + hat(angle_axis); the Taylor form keeps the derivatives + // meaningful where the angle-axis vanishes. + const T w_cross_pt[3] = {aa[1] * pt[2] - aa[2] * pt[1], + aa[2] * pt[0] - aa[0] * pt[2], + aa[0] * pt[1] - aa[1] * pt[0]}; + result[0] = pt[0] + w_cross_pt[0]; + result[1] = pt[1] + w_cross_pt[1]; + result[2] = pt[2] + w_cross_pt[2]; + return; + } + + const T w_cross_pt[3] = {w[1] * pt[2] - w[2] * pt[1], + w[2] * pt[0] - w[0] * pt[2], + w[0] * pt[1] - w[1] * pt[0]}; + const T tmp = (w[0] * pt[0] + w[1] * pt[1] + w[2] * pt[2]) * (T(1.0) - costheta); + + result[0] = pt[0] * costheta + w_cross_pt[0] * sintheta + w[0] * tmp; + result[1] = pt[1] * costheta + w_cross_pt[1] * sintheta + w[1] * tmp; + result[2] = pt[2] * costheta + w_cross_pt[2] * sintheta + w[2] * tmp; + } + + const T aa[3]; + T w[3]; + T costheta, sintheta; + bool at_zero; +}; + // Detector -> reciprocal geometry residual, shared by the per-image XtalOptimizer (one lattice, one // frame) and the offline GeometryRefiner (shared beam/distance/cell blocks, one orientation block per // frame). Parameter blocks: beam(2), distance_mm(1), detector_rot(2 = rot1,rot2), rotation_axis(3), @@ -32,7 +88,8 @@ struct XtalResidual { : obs_x(x), obs_y(y), inv_lambda(1.0/lambda), pixel_size(pixel_size), - rot3(rot3), + cos_rot3(std::cos(rot3)), + sin_rot3(std::sin(rot3)), exp_h(exp_h), exp_k(exp_k), exp_l(exp_l), @@ -67,9 +124,11 @@ struct XtalResidual { const T c2 = ceres::cos(rot2); const T s2 = ceres::sin(rot2); - // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0 - const T c3 = T(cos(rot3)); - const T s3 = T(sin(rot3)); + // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0. Its sine and cosine + // are taken in the constructor - they do not depend on any parameter, so recomputing them per + // evaluation is two libm calls per residual for nothing. + const T c3 = T(cos_rot3); + const T s3 = T(sin_rot3); // Detector coordinates in mm const T det_x = (T(obs_x) - beam[0]) * T(pixel_size); @@ -170,9 +229,10 @@ struct XtalResidual { T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2}; T col0_rot[3], col1_rot[3], col2_rot[3]; - ceres::AngleAxisRotatePoint(p0, col0_unrot, col0_rot); - ceres::AngleAxisRotatePoint(p0, col1_unrot, col1_rot); - ceres::AngleAxisRotatePoint(p0, col2_unrot, col2_rot); + const AngleAxisRotator rot_p0(p0); + rot_p0.Rotate(col0_unrot, col0_rot); + rot_p0.Rotate(col1_unrot, col1_rot); + rot_p0.Rotate(col2_unrot, col2_rot); const Eigen::Matrix A(col0_rot[0], col0_rot[1], col0_rot[2]); const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); @@ -205,7 +265,7 @@ struct XtalResidual { const double obs_x, obs_y; const double inv_lambda; const double pixel_size; - const double rot3; + const double cos_rot3, sin_rot3; const double exp_h; const double exp_k; const double exp_l; diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 1276ba2b..16965c50 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -39,6 +39,12 @@ namespace { // How far from a whole number the volume ratio may sit and still count as an axis multiple. constexpr double ROT_SUPERCELL_INTEGER_TOL = 0.15; + // Iteration bound for the candidate-cell refinements when the run is not real-time. Without one, + // XtalOptimizerData falls back to its wall-clock bound and the cell a candidate refines to - and + // therefore which candidate wins - depends on how busy the machine was. Same value as Ceres' own + // default iteration limit, so it only bites where the clock was biting before. + constexpr int ROT_REFINE_ITERATIONS = 50; + // Order of the lattice point group, so "lower symmetry" is a well-defined comparison // (gemmi's enum orders Trigonal after Tetragonal, which have 6 and 8 rotations). int LatticePointGroupOrder(gemmi::CrystalSystem s) { @@ -106,9 +112,11 @@ namespace { } } -RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) +RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer, + bool real_time) : experiment(x), index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), + real_time(real_time), v_(experiment.GetImageNum()), angle_deg_(experiment.GetImageNum()), axis_(x.GetGoniometer()), @@ -208,6 +216,7 @@ void RotationIndexer::RunIndexing() { .refine_detector_angles = true, .refine_rotation_axis = true, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), + .max_iterations = real_time ? 0 : ROT_REFINE_ITERATIONS, .axis = orig_axis }; if (d.crystal_system == gemmi::CrystalSystem::Trigonal) @@ -421,6 +430,7 @@ void RotationIndexer::RunIndexing() { .refine_unit_cell = false, .refine_rotation_axis = false, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), + .max_iterations = real_time ? 0 : ROT_REFINE_ITERATIONS, .axis = axis_ }; diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index 1fe5de58..d9b8dd51 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -34,6 +34,7 @@ class RotationIndexer { constexpr static size_t max_spots_per_image = 200; const bool index_ice_rings; + const bool real_time; // see the constructor std::vector> v_; // Per-image rotation angle (mid-exposure, deg) supplied by the caller; falls back to the @@ -53,7 +54,11 @@ class RotationIndexer { std::optional indexed_lattice; public: - RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); + // real_time: bound the candidate-cell refinements by WALL CLOCK, as online acquisition must - it + // has a real budget. Offline (rugnux, the viewer) passes false and they are bounded by iteration + // count instead, so the same file reprocesses to the same cell regardless of what else the machine + // was doing at the time. + RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer, bool real_time = false); // angle_deg is the image's mid-exposure rotation angle; if omitted, the goniometer angle at // `image` is used (only valid when `image` is the goniometer's own image index). void ProcessImage(int64_t image, const std::vector& spots,