Lattice class: carry the free class and decide it on what the pass realises

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 on everything below: the cell is refined with a = b, reflections
are predicted from it, and they are integrated at those predicted positions.
On a small cell the tolerance is far wider than the spot positions resolve, so
a genuinely orthorhombic crystal whose a and b differ by 2 % is integrated as
tetragonal and every decision downstream is read off the wreckage.

Measure the equality instead of assuming it. The indexer already refines each
candidate a second time with nothing held, so the first pass now reports how
far that free refinement leaves the two axes apart. A relative split of eps
displaces a reflection at radius r by eps/2 * r pixels: where that displacement
at the far corner of the detector stays inside the integration disc, imposing
the equality moves nothing out of its own box and the higher symmetry is kept
with no extra pass. Where it does not, both hypotheses are run as probe-only
passes - the promoted class, and the class the same walk carries when it is
granted no length equality at all - and whichever realises the lower held-out
positional residual is kept, the promoted class on a tie.

Two small-molecule sets whose axes differ by 2 % now index, refine and merge in
their own orthorhombic lattice instead of a tetragonal mean: one goes from a
cell 1 % wrong and P 1 at CC1/2 0.19 to the deposited cell within 0.2 % and
P 2 2 2 at CC1/2 0.95, the other from a tetragonal mean to a cell matching its
reference to 0.4 %. Protein sets whose symmetry is real (P41212, P4222, I23,
P6422, F4132) are unchanged: their free refinements leave the axes 0.02-0.10 %
apart, a few tenths of a pixel, so the question is never asked. Where it was
asked on a weak sweep whose free refinement diverged, the arms decided for the
higher symmetry and the output was identical.

No new threshold: the comparison is the integration radius the run already
integrates at, and the arms are judged by HeldOutResidualFell, as the geometry
walk's rounds are. postrefine_probe_only_ returns for the arms' sake - a pass
run only to measure stops before the scaling engine is built, so neither arm
pays for a merge or a space-group search.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT

common/IndexingSettings.cpp

@@ -223,6 +223,15 @@ float IndexingSettings::GetRotationIndexingAngularStride_deg() const {
     return rotation_indexing_angular_stride_deg;
 }
 
+IndexingSettings & IndexingSettings::LatticeClassDemoted(bool input) {
+    lattice_class_demoted = input;
+    return *this;
+}
+
+bool IndexingSettings::GetLatticeClassDemoted() const {
+    return lattice_class_demoted;
+}
+
 bool IndexingSettings::GetIndexIceRings() const {
     return index_ice_rings;
 }

common/IndexingSettings.h

@@ -62,6 +62,14 @@ private:
     bool blocking_behavior = true;
     bool index_ice_rings = false;
 
+    // Refuse the Bravais class the metric was promoted to on a LENGTH equality, and index in the
+    // class left when that equality is struck out. Not a user setting: rugnux runs one pass with it
+    // and one without, and keeps the one whose re-integrated held-out positional residual is lower
+    // (see RunAllPasses). The promotion is decided against a 3 % relative tolerance on a cell the
+    // spot positions place to ~0.1 %, so on a small cell a genuinely orthorhombic crystal is
+    // integrated with a = b imposed and nothing downstream can recover from it.
+    bool lattice_class_demoted = false;
+
     bool enable_rotation_indexing = false;
     float rotation_indexing_min_angular_range_deg = 20.0;
     float rotation_indexing_angular_stride_deg = 0.5;
@@ -89,6 +97,7 @@ public:
     IndexingSettings& RotationIndexingAngularStride_deg(float input);
     IndexingSettings& BlockingBehavior(bool input);
     IndexingSettings& MaxExtraLattices(int64_t input);
+    IndexingSettings& LatticeClassDemoted(bool input);
 
     [[nodiscard]] int64_t GetViableCellMinSpots() const;
     [[nodiscard]] IndexingAlgorithmEnum GetAlgorithm() const;
@@ -110,4 +119,5 @@ public:
     [[nodiscard]] float GetRotationIndexingAngularStride_deg() const;
     [[nodiscard]] bool GetBlockingBehavior() const;
     [[nodiscard]] int64_t GetMaxExtraLattices() const;
+    [[nodiscard]] bool GetLatticeClassDemoted() const;
 };

image_analysis/lattice_search/LatticeSearch.cpp

@@ -476,6 +476,19 @@ std::optional<LatticeSearchResult> LatticeSearchForClass(const CrystalLattice &L
     return SearchCharacters(L, dist_tolerance, angle_tolerance, &only);
 }
 
+
+bool ClassImposesLengthEquality(gemmi::CrystalSystem system) {
+    switch (system) {
+        case gemmi::CrystalSystem::Tetragonal:
+        case gemmi::CrystalSystem::Trigonal:
+        case gemmi::CrystalSystem::Hexagonal:
+        case gemmi::CrystalSystem::Cubic:
+            return true;
+        default:
+            return false;
+    }
+}
+
 double MetricDeviation(const LatticeSearchResult &sr) {
     const UnitCell uc = sr.conventional.GetUnitCell();
     double angle = 0.0, length = 0.0;

image_analysis/lattice_search/LatticeSearch.h

@@ -48,6 +48,14 @@ std::optional<LatticeSearchResult> LatticeSearchForClass(const CrystalLattice& L
                                                          double dist_tolerance = LATTICE_SEARCH_DIST_TOLERANCE,
                                                          double angle_tolerance_deg = LATTICE_SEARCH_ANGLE_TOLERANCE_DEG);
 
+// Whether the conventional metric of this class asserts an equality between cell LENGTHS (a = b, and
+// for cubic a = b = c). Those are the promotions the relative tolerance above admits far more widely
+// than the spot positions resolve; every other class is carried by its angles alone. The walk run
+// with dist_tolerance = 0 grants no such equality, and so returns the class the metric carries on
+// its angles alone - the OTHER hypothesis about any of these promotions, which rugnux runs as its
+// own pass and judges on what it realises (see RunAllPasses).
+bool ClassImposesLengthEquality(gemmi::CrystalSystem system);
+
 // How far the matched lattice's metric sits from the IDEAL metric of its class, as a fraction of the
 // tolerance that admitted it: 0 is exactly the class, 1 is the edge of what the walk accepts. The
 // walk computes this and throws it away - a cell 0.02 deg from being C-centred monoclinic and one

image_analysis/rotation_indexer/RotationIndexer.cpp

@@ -160,8 +160,18 @@ void RotationIndexer::RunIndexing() {
         // 60 validation frames and an F-cubic one relabelled trigonal indexed 0 of 60, where the same
         // frames index at 36/60 and 51/60 without a group. The group is applied where it belongs - to
         // the scaling and the merge.
+        const bool demote_class = experiment.GetIndexingSettings().GetLatticeClassDemoted();
         auto build_sr = [&](const CrystalLattice &cand) -> LatticeSearchResult {
             auto ls = LatticeSearch(cand);
+            // The other hypothesis about a promotion that rests on a = b: walk the characters again
+            // granting NO length equality at all, and take the class the metric carries on its angles
+            // alone. Not the same as striking one class out - the tetragonal-P lattice also has a
+            // C-centred orthorhombic character, which is a different (doubled) cell and rests on the
+            // very equality in question. The caller runs this as its own pass and keeps whichever pass
+            // realises the lower held-out positional residual. The hexagonal re-expression below is a
+            // promotion of the same kind, so it is not applied here.
+            if (demote_class && ClassImposesLengthEquality(ls.system))
+                return LatticeSearch(cand, /*dist_tolerance=*/0.0);
             if (!IsHexagonalSystem(ls.system) && IsMetricallyHexagonal(ls.primitive_reduced)) {
                 ls.conventional = HexagonalConventional(ls.primitive_reduced);
                 ls.system = gemmi::CrystalSystem::Hexagonal;

rugnux/Rugnux.cpp

@@ -2086,6 +2086,106 @@ 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:
+        // the cell is refined with a = b, reflections are predicted from it and integrated at the
+        // predicted positions. On a small cell the tolerance is far wider than the positions resolve -
+        // a 5.8 A axis pair 2.4 % apart is read as tetragonal, and the mean the constraint substitutes
+        // mispredicts the outermost reflections by a good ten pixels against an integration disc of
+        // four - so a genuinely orthorhombic crystal is integrated as tetragonal and every decision
+        // downstream, up to the space group, is read off the wreckage.
+        //
+        // So the equality is measured rather than assumed. The indexer already refines the same
+        // candidate a second time with nothing held, and pass 1 reports how far that free refinement
+        // leaves the two axes apart (lattice_length_split). A relative split of eps displaces a
+        // reflection at radius r by eps/2 * r pixels, so the question "do these data resolve it" has
+        // an answer in the units the integration works in: where that displacement at the far corner
+        // of the detector stays inside the integration disc, imposing the equality moves nothing out
+        // of its own box and the higher symmetry is kept with no pass run. Where it does not, both
+        // hypotheses are run - the promoted class, and the class the same walk carries when it is
+        // granted no length equality at all - and whichever REALISES the lower held-out residual on a
+        // re-integrated, re-indexed pass is kept, the
+        // promoted class on a tie, because it is the one that asserts more. Both arms are probe-only
+        // passes (postrefine_probe_only_): the residual is all the comparison reads and it is in hand
+        // before the scaling engine is built, so neither arm pays for a merge, a space-group search or
+        // a set of correction surfaces, and only the class that wins is ever run as a canonical pass.
+        std::string lattice_arm;      // the arm decision, for the pass decision text
+        if (!cancelled_ && pass1.lattice_length_split && pass1.rotation_lattice_type) {
+            const double split = *pass1.lattice_length_split;
+            const double dx = std::max(experiment_.GetBeamX_pxl(),
+                                       static_cast<float>(experiment_.GetXPixelsNum()) - experiment_.GetBeamX_pxl());
+            const double dy = std::max(experiment_.GetBeamY_pxl(),
+                                       static_cast<float>(experiment_.GetYPixelsNum()) - experiment_.GetBeamY_pxl());
+            const double r_corner_px = std::hypot(dx, dy);
+            const double miss_px = 0.5 * split * r_corner_px;
+            const double disc_px = experiment_.GetBraggIntegrationSettings().GetR1();
+            if (miss_px > disc_px) {
+                logger.Info("Metric symmetry: {} was promoted on a = b, but the free refinement leaves "
+                            "the two axes {:.2f} % apart - {:.1f} px of misprediction at the far corner "
+                            "against an integration disc of {:.1f} px. Both classes are re-integrated "
+                            "and judged on the held-out residual each realises",
+                            gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system),
+                            100.0 * split, miss_px, disc_px);
+                // The arms are compared on the post-refinement's held-out residual, so they have to
+                // measure one: a canonical pass only does where the pre-pass moved the geometry.
+                const bool probe_was = postrefine_probe_;
+                postrefine_probe_ = true;
+                postrefine_probe_only_ = true;
+                const auto promoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false);
+                ++arm_passes;
+                IndexingSettings demoted_settings = experiment_.GetIndexingSettings();
+                demoted_settings.LatticeClassDemoted(true);
+                experiment_.ImportIndexingSettings(demoted_settings);
+                const auto demoted_arm = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/false);
+                ++arm_passes;
+                postrefine_probe_only_ = false;
+                postrefine_probe_ = probe_was;
+                const double realised_promoted = promoted_arm.post_refine ? promoted_arm.post_refine->held_out_before : NAN;
+                const double realised_demoted = demoted_arm.post_refine ? demoted_arm.post_refine->held_out_before : NAN;
+                const auto class_text = [](const ProcessResult &r) {
+                    return r.rotation_lattice_type
+                           ? fmt::format("{}-centred {}", r.rotation_lattice_type->centering,
+                                         gemmi::crystal_system_str(r.rotation_lattice_type->crystal_system))
+                           : std::string("no lattice");
+                };
+                if (!cancelled_ && promoted_arm.post_refine && demoted_arm.post_refine
+                    && HeldOutResidualFell(*promoted_arm.post_refine, *demoted_arm.post_refine)) {
+                    logger.Warning("Metric symmetry: re-integrated as {} the held-out residual is {:.3e} "
+                                   "against {:.3e} as {} - the equality is not this crystal's, and the "
+                                   "run goes on without it", class_text(demoted_arm), realised_demoted,
+                                   realised_promoted, class_text(promoted_arm));
+                    lattice_arm = fmt::format(
+                            "the metric symmetry was promoted to {} on axes the free refinement leaves "
+                            "{:.2f} % apart; re-integrated as {} the held-out residual is {:.3e} against "
+                            "{:.3e}, so the equality was dropped", class_text(promoted_arm), 100.0 * split,
+                            class_text(demoted_arm), realised_demoted, realised_promoted);
+                } else if (!cancelled_) {
+                    IndexingSettings promoted_settings = experiment_.GetIndexingSettings();
+                    promoted_settings.LatticeClassDemoted(false);
+                    experiment_.ImportIndexingSettings(promoted_settings);
+                    logger.Info("Metric symmetry: re-integrated as {} the held-out residual is {:.3e}, "
+                                "which {} does not beat ({:.3e}) - the higher symmetry stands",
+                                class_text(promoted_arm), realised_promoted, class_text(demoted_arm),
+                                realised_demoted);
+                    lattice_arm = fmt::format(
+                            "{} was run as well and realised a held-out residual of {:.3e} against {:.3e} "
+                            "as {}, so the metric symmetry stands", class_text(demoted_arm),
+                            realised_demoted, realised_promoted, class_text(promoted_arm));
+                }
+            } else {
+                logger.Info("Metric symmetry: {} asserts a = b and the free refinement leaves the two "
+                            "axes {:.2f} % apart - {:.1f} px at the far corner, inside the integration "
+                            "disc of {:.1f} px, so these data do not resolve the difference",
+                            gemmi::crystal_system_str(pass1.rotation_lattice_type->crystal_system),
+                            100.0 * split, miss_px, disc_px);
+            }
+        }
+
         logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
                     base_prefix);
         auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
@@ -2205,8 +2305,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
             pass1.post_refine->distance_held = false;
         }
         pass2.post_refine = pass1.post_refine;
-        pass2.pass_number = 2 + walk_passes;
-        pass2.pass_count = 2 + walk_passes;
+        pass2.pass_number = 2 + arm_passes + walk_passes;
+        pass2.pass_count = 2 + arm_passes + walk_passes;
         if (geometry_rounds > 0)
             pass2.pass_decision = fmt::format(
                     "post-refined geometry walked {} round{} and kept round {}, detector distance {:.3f} mm "
@@ -2341,6 +2441,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
             pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_)
                     ? "post-refined geometry adopted"
                     : "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;
         pass2.geometry_not_converged = geometry_not_converged;
         return pass2;
     }
@@ -4869,6 +4971,28 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             };
             result.rotation_lattice_found = true;
             result.rotation_lattice_type = end_msg.rotation_lattice_type;
+            // What the equality this class asserts is worth on a REFINED cell. The indexer refines the
+            // same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that
+            // cell over whenever the constrained one won; it was seeded from the primitive reduced cell
+            // and comes back in that basis, so the class's own reindex re-expresses it in the
+            // conventional one and a and b can be read off it. Measured here, decided in RunAllPasses,
+            // which is where there is a pass to compare it against.
+            if (rot->unconstrained && ClassImposesLengthEquality(rot->search_result.system)) {
+                const UnitCell free_uc =
+                        rot->unconstrained->lattice.Multiply(rot->search_result.reindex).GetUnitCell();
+                double split = std::fabs(free_uc.a - free_uc.b) / (0.5 * (free_uc.a + free_uc.b));
+                if (rot->search_result.system == gemmi::CrystalSystem::Cubic)
+                    split = std::max(split, std::fabs(free_uc.b - free_uc.c)
+                                            / (0.5 * (free_uc.b + free_uc.c)));
+                if (std::isfinite(split)) {
+                    result.lattice_length_split = split;
+                    logger.Info("Rotation lattice: {}-centred {} asserts a = b; the free refinement of "
+                                "the same candidate leaves them {:.2f} % apart ({:.4f} / {:.4f} A)",
+                                rot->search_result.centering,
+                                gemmi::crystal_system_str(rot->search_result.system),
+                                100.0 * split, free_uc.a, free_uc.b);
+                }
+            }
             // The consensus cell below is this lattice's conventional cell, so this is the centring
             // it is written in - the pair the volume comparisons in Run() need.
             result.consensus_centering = rot->search_result.centering;
@@ -5214,8 +5338,20 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
         // before the engine is built rather than beside it - the gathered observations and their
         // rocking events are gigabytes on a long axis, and beside the engine's own arrays they set
         // the pass's memory high-water mark.
-        if (postrefine_probe_ && !geometry_prepass && is_rotation)
+        if (postrefine_probe_ && !geometry_prepass && is_rotation) {
             post_refine_geometry();
+            // A pass the lattice arms run to compare two Bravais classes has now produced the only
+            // thing the comparison reads - the held-out residual of the fit just measured. Stop here
+            // rather than scale, search, fit surfaces and merge for a result nothing looks at.
+            if (postrefine_probe_only_) {
+                result.processing_time_s = std::chrono::duration<double>(
+                        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)",
+                            result.images_processed, result.processing_time_s);
+                return result;
+            }
+        }
         if (is_rotation) {
             if (experiment_.GetRefineRotationWedgeInScaling()
                 || rot_ss.GetRotationWedgeForScaling().has_value())

rugnux/Rugnux.h

@@ -232,6 +232,13 @@ struct ProcessResult {
     // 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<LatticeMessage> rotation_lattice_type;
+    // Where that class asserts a = b: how far apart the FREE (triclinic) refinement of the same
+    // candidate - which the indexer runs anyway, as its pseudo-symmetry guard - leaves those two
+    // axes, relative. The class was decided against a fixed relative tolerance on an UNREFINED FFT
+    // cell; this is the same equality measured on a refined one, and it is what RunAllPasses turns
+    // into pixels to ask whether this crystal's spot positions resolve the equality at all. Empty
+    // where the class asserts no length equality, or where the free refinement is not available.
+    std::optional<double> lattice_length_split;
     MeanProcessingTime mean_processing_time{};
     std::optional<std::string> written_master_path;
     std::string merge_statistics_text;   // populated when scaling/merging ran
@@ -555,6 +562,12 @@ class Rugnux {
     // whether it has finished walking is a question only reflections integrated after the move can
     // answer; RunAllPasses turns this on for the passes that answer it.
     bool postrefine_probe_ = false;
+    // Whether such a pass exists ONLY for that measurement. The lattice-class arms (see RunAllPasses)
+    // read the held-out residual the probe returns and nothing else, so a pass run for one of them
+    // stops as soon as the probe has measured - which is before the scaling engine is built.
+    // 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;
     // 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.

tests/LatticeSearchTest.cpp

@@ -639,3 +639,32 @@ TEST_CASE("SymmetrizeMetric puts a cell onto the metric its group fixes") {
     CHECK(turned_fixed.Vec0() * turned.Vec0()
           == Catch::Approx(turned_fixed.Vec0().Length() * turned.Vec0().Length()).epsilon(1e-4));
 }
+
+TEST_CASE("LatticeSearch - a pseudo-tetragonal cell is tetragonal at the walk's tolerance and "
+          "orthorhombic with no length equality granted") {
+    // a and b 2.4 % apart: inside LATTICE_SEARCH_DIST_TOLERANCE, far outside what spot positions
+    // resolve. The two hypotheses rugnux carries forward are exactly these two answers.
+    const double a = 5.795, b = 5.933, c = 12.362;
+    CrystalLattice L(Coord(a, 0, 0), Coord(0, b, 0), Coord(0, 0, c));
+
+    const auto promoted = LatticeSearch(L);
+    CHECK(promoted.system == gemmi::CrystalSystem::Tetragonal);
+    CHECK(promoted.centering == 'P');
+    CHECK(ClassImposesLengthEquality(promoted.system));
+
+    const auto free_class = LatticeSearch(L, /*dist_tolerance=*/0.0);
+    CHECK(free_class.system == gemmi::CrystalSystem::Orthorhombic);
+    CHECK(free_class.centering == 'P');
+    const auto uc = free_class.conventional.GetUnitCell();
+    check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-2, 1e-2);
+}
+
+TEST_CASE("ClassImposesLengthEquality names the classes whose metric asserts a = b") {
+    CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Tetragonal));
+    CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Trigonal));
+    CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Hexagonal));
+    CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Cubic));
+    CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Orthorhombic));
+    CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Monoclinic));
+    CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Triclinic));
+}