Merge branch 'd-lattice-class' into pool-C

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