RESULT-CHANGING: -C: an integer supercell of the given cell is stepped down to it, not refused

A candidate is kept only when it matches the given cell, and on sparse twinned data the FFT never
returns the cell itself: it returns the twins' coincidence lattice (a supercell of either domain) or a
cell built from rows taken at a multiple of their length. With the true cell given, such a run found
no lattice at all. Every node of a supercell is a node of the given lattice, so the given basis is a
combination of the candidate's rows with coefficients in (1/n)Z, n the volume ratio (2 to 10); those
combinations are enumerated and, of the triples that match the given cell, span its volume and contain
the candidate's lattice, the one indexing the most spots replaces the candidate (on a twin, the stronger
domain). Measured: a reticular-twin small-molecule sweep given its true P2_1 cell went from "no
lattice" to P2_1 on that cell.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqrcrEZy6XBUCMZvyT

docs/RUGNUX_CLI.md

@@ -256,7 +256,7 @@ For very long axes the direction grid's angular resolution binds well below 500
 | `--force-still` | Treat a rotation (goniometer) dataset as independent stills instead of rotation |
 | `-X, --indexing-algorithm <txt>` | `Auto` (default) \| `FFT` \| `cuFFT` \| `FFTW` \| `FFBIDX` \| `None`, case-insensitive. `FFT` is the de-novo FFT search, on the GPU (cuFFT) where there is one and `--no-gpu` is not given, otherwise on the CPU (FFTW); `cuFFT` and `FFTW` name one implementation each; `cuFFT` needs a GPU and is refused without one. `FFBIDX` (known cell, `-C`) runs on the GPU where there is one and otherwise on the CPU. `Auto` takes `FFBIDX` on stills with a known cell, and `FFT` otherwise; rotation data always use the FFT search |
 | `-F` | Short for `-X FFT`, kept for old command lines |
-| `-C, --unit-cell <cell>` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it |
+| `-C, --unit-cell <cell>` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it. A candidate that is an integer supercell of the cell given (a twin's coincidence lattice, a row taken at a multiple of its length) is stepped down to it |
 | `--fft-min-unit-cell <num>` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but a de-novo rotation run also runs a **second first-pass hypothesis with the floor lowered to 5 Å by default**, adopted when the standing cell turns out to be an integer supercell of what it finds, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given (and the second hypothesis then stays out of the way). On stills, where neither applies, a crystal below the floor cannot be indexed unless this is lowered |
 | `--min-indexed-spots <num>` | Spots a frame must have on the lattice before it counts as indexed (default: 6, minimum 4). It sets the reported indexing rate and the count the rotation first pass ranks candidate lattices by; whether the run has a lattice at all is decided on the pooled spots instead ([§4.1](CPU_DATA_ANALYSIS_INDEXING.md#41-indexed-spot-decision-inlier-test)), and integration is not gated by it |
 | `-S, --space-group <num\|symbol>` | Space group number (`96`) or Hermann-Mauguin symbol (`P43212`) — for indexing and scaling. On serial stills whose group admits two indexings, it also lets the crystals be put on one from the data ([the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)) |

image_analysis/indexing/PostIndexingRefinement.cpp

@@ -7,6 +7,7 @@
 #include <iostream>
 #include <limits>
 #include <future>
+#include <optional>
 
 namespace {
     struct config_ifssr final {
@@ -142,6 +143,106 @@ namespace {
             cells.block(3u * j, 0u, 3u, 3u) = cell.transpose();
         }
     }
+
+    // The reference test, on a candidate's three rows: sorted lengths within the relative tolerance and
+    // sorted angles within REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG of any of the reference cells. Each
+    // angle is folded to its acute complement (min(x,180-x)) so the obtuse/acute setting choice is
+    // irrelevant. Guards against a right-edges/wrong-angle cell (a pseudo-symmetric near-metric, e.g. a
+    // monoclinic beta refined to the wrong value) passing on lengths.
+    bool MatchesReference(const Eigen::Matrix3f &cell_rows, const std::vector<UnitCell> &reference_cells,
+                          float dist_tolerance) {
+        std::array<float, 3> obs = {cell_rows.row(0).norm(), cell_rows.row(1).norm(), cell_rows.row(2).norm()};
+        std::sort(obs.begin(), obs.end());
+
+        auto fold = [](float deg) { return std::min(deg, 180.0f - deg); };
+        auto row_angle = [&](int i, int j) {
+            return std::acos(std::clamp(cell_rows.row(i).normalized().dot(cell_rows.row(j).normalized()),
+                                        -1.0f, 1.0f)) * 180.0f / PI;
+        };
+        std::array<float, 3> obs_ang = {fold(row_angle(1, 2)), fold(row_angle(0, 2)), fold(row_angle(0, 1))};
+        std::sort(obs_ang.begin(), obs_ang.end());
+
+        auto matches = [&](const UnitCell &reference) {
+            std::array<float, 3> ref = {reference.a, reference.b, reference.c};
+            std::sort(ref.begin(), ref.end());
+            for (int k = 0; k < 3; ++k) {
+                const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON);
+                if (std::abs(obs[k] - ref[k]) / denom > dist_tolerance)
+                    return false;
+            }
+
+            std::array<float, 3> ref_ang = {fold(reference.alpha), fold(reference.beta), fold(reference.gamma)};
+            std::sort(ref_ang.begin(), ref_ang.end());
+            for (int k = 0; k < 3; ++k) {
+                if (std::abs(obs_ang[k] - ref_ang[k]) > REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG)
+                    return false;
+            }
+            return true;
+        };
+        return std::any_of(reference_cells.begin(), reference_cells.end(), matches);
+    }
+
+    // A candidate that is an integer SUPERCELL of the reference lattice: the FFT builds one from rows it
+    // reports at a multiple of their length, and on a twin by reticular merohedry the cell indexing both
+    // domains is their coincidence lattice, a supercell of either. Every node of the candidate is then a
+    // node of the reference lattice, so each reference basis vector is a combination of the candidate's
+    // rows with coefficients in (1/n)Z, n the volume ratio. Those combinations are enumerated, and of the
+    // triples that match the reference, span its volume and contain the candidate's lattice, the one
+    // that indexes the most spots is returned (on a twin, the stronger domain).
+    std::optional<Eigen::Matrix3f> StepDownToReference(const Eigen::Matrix3f &cell_rows,
+                                                       const std::vector<UnitCell> &reference_cells,
+                                                       float dist_tolerance,
+                                                       const Eigen::Ref<const Eigen::MatrixX3<float>> &spots,
+                                                       float indexing_tolerance) {
+        constexpr int MAX_INDEX = 10;
+        const float vol = std::abs(cell_rows.determinant());
+        std::optional<Eigen::Matrix3f> best;
+        int64_t best_count = -1;
+        for (const auto &ref: reference_cells) {
+            const float ref_vol = CrystalLattice(ref).CalcVolume();
+            if (!(ref_vol > 0.0f))
+                continue;
+            const float ratio = vol / ref_vol;
+            const int n = static_cast<int>(std::lround(ratio));
+            if (n < 2 || n > MAX_INDEX || std::abs(ratio - static_cast<float>(n)) > 0.15f * static_cast<float>(n))
+                continue;
+
+            const std::array<float, 3> lengths = {ref.a, ref.b, ref.c};
+            std::array<std::vector<Eigen::RowVector3f>, 3> near;
+            for (int i = -n; i <= n; i++)
+                for (int j = -n; j <= n; j++)
+                    for (int k = -n; k <= n; k++) {
+                        const Eigen::RowVector3f v = (static_cast<float>(i) * cell_rows.row(0)
+                                                      + static_cast<float>(j) * cell_rows.row(1)
+                                                      + static_cast<float>(k) * cell_rows.row(2))
+                                                     / static_cast<float>(n);
+                        for (int a = 0; a < 3; a++)
+                            if (std::abs(v.norm() - lengths[a]) <= dist_tolerance * lengths[a])
+                                near[a].push_back(v);
+                    }
+
+            for (const auto &v0: near[0])
+                for (const auto &v1: near[1])
+                    for (const auto &v2: near[2]) {
+                        Eigen::Matrix3f sub;
+                        sub << v0, v1, v2;
+                        if (std::abs(std::abs(sub.determinant()) - ref_vol) > 0.15f * ref_vol)
+                            continue;
+                        const Eigen::Matrix3f m = cell_rows * sub.inverse();
+                        if ((m - m.array().round().matrix()).cwiseAbs().maxCoeff() > 0.1f)
+                            continue;
+                        if (!MatchesReference(sub, {ref}, dist_tolerance))
+                            continue;
+                        int64_t count = 0;
+                        ComputeIndexedMask(spots, sub.transpose(), indexing_tolerance, count);
+                        if (count > best_count) {
+                            best_count = count;
+                            best = sub;
+                        }
+                    }
+        }
+        return best;
+    }
 }
 
 std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
@@ -216,41 +317,15 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
         Eigen::Vector3f row_norms = cell_rows.rowwise().norm();
 
         if (!reference_cells.empty()) {
-            std::array<float, 3> obs = {row_norms(0), row_norms(1), row_norms(2)};
-            std::sort(obs.begin(), obs.end());
-
-            // The angles are compared as well as the lengths. Each is folded to its acute complement
-            // (min(x,180-x)) so the obtuse/acute setting choice is irrelevant, then the sorted
-            // triples are compared. Guards against a right-edges/wrong-angle cell (a pseudo-symmetric
-            // near-metric, e.g. a monoclinic beta refined to the wrong value) passing on lengths.
-            auto fold = [](float deg) { return std::min(deg, 180.0f - deg); };
-            auto row_angle = [&](int i, int j) {
-                return std::acos(std::clamp(cell_rows.row(i).normalized().dot(cell_rows.row(j).normalized()),
-                                            -1.0f, 1.0f)) * 180.0f / PI;
-            };
-            std::array<float, 3> obs_ang = {fold(row_angle(1, 2)), fold(row_angle(0, 2)), fold(row_angle(0, 1))};
-            std::sort(obs_ang.begin(), obs_ang.end());
-
-            auto matches = [&](const UnitCell &reference) {
-                std::array<float, 3> ref = {reference.a, reference.b, reference.c};
-                std::sort(ref.begin(), ref.end());
-                for (int k = 0; k < 3; ++k) {
-                    const float denom = std::max(ref[k], REFINE_MIN_REFERENCE_LENGTH_EPSILON);
-                    if (std::abs(obs[k] - ref[k]) / denom > p.dist_tolerance_vs_reference)
-                        return false;
-                }
-
-                std::array<float, 3> ref_ang = {fold(reference.alpha), fold(reference.beta), fold(reference.gamma)};
-                std::sort(ref_ang.begin(), ref_ang.end());
-                for (int k = 0; k < 3; ++k) {
-                    if (std::abs(obs_ang[k] - ref_ang[k]) > REFINE_ANGLE_TOLERANCE_VS_REFERENCE_DEG)
-                        return false;
-                }
-                return true;
-            };
-
-            if (std::none_of(reference_cells.begin(), reference_cells.end(), matches))
-                continue;
+            if (!MatchesReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference)) {
+                const auto sub = StepDownToReference(cell_rows, reference_cells, p.dist_tolerance_vs_reference,
+                                                     spots.topRows(nspots), p.indexing_tolerance);
+                if (!sub)
+                    continue;
+                cell_rows = *sub;
+                cell_cols = cell_rows.transpose();
+                row_norms = cell_rows.rowwise().norm();
+            }
         } else {
             if (row_norms.minCoeff() < p.min_length_A || row_norms.maxCoeff() > p.max_length_A)
                 continue;

tests/IndexingUnitTest.cpp

@@ -280,6 +280,46 @@ TEST_CASE("FFTIndexer_AxisLongerThanTheBound","[Indexing]") {
     CHECK(out[2] == Catch::Approx(360.0f).epsilon(0.01));
 }
 
+TEST_CASE("PostIndexingRefinement_SupercellStepsDownToReference","[Indexing]") {
+    // A candidate that is a 2x supercell of the given cell (one axis doubled, as the coincidence lattice
+    // of a twin or a row reported at its first multiple gives) comes back as the given cell, not refused.
+    const UnitCell uc(6.0, 8.0, 11.0, 90, 100, 90);
+    const CrystalLattice cl(uc);
+
+    std::vector<Coord> spots;
+    for (int h = -5; h <= 5; h++)
+        for (int k = -6; k <= 6; k++)
+            for (int l = -8; l <= 8; l++)
+                if (h != 0 || k != 0 || l != 0)
+                    spots.push_back(h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar());
+
+    Eigen::MatrixX3<float> cells(3, 3);
+    const Coord a2 = cl.Vec0() * 2.0f, b = cl.Vec1(), c = cl.Vec2();
+    cells << a2.x, a2.y, a2.z, b.x, b.y, b.z, c.x, c.y, c.z;
+    Eigen::VectorX<float> scores(1);
+    scores(0) = 0.2f;
+
+    RefineParameters p{
+        .viable_cell_min_spots = 6,
+        .dist_tolerance_vs_reference = 0.05f,
+        .reference_unit_cell = uc,
+        .min_length_A = 2.5f,
+        .max_length_A = 500.0f,
+        .min_angle_deg = 30.0f,
+        .max_angle_deg = 150.0f,
+        .indexing_tolerance = 0.1f
+    };
+    const auto out = Refine(spots, spots.size(), cells, scores, p);
+    REQUIRE(out.size() == 1);
+    CHECK(std::abs(out[0].CalcVolume()) == Catch::Approx(cl.CalcVolume()).epsilon(0.01));
+    const auto uc_out = out[0].GetUnitCell();
+    std::array<float, 3> lengths = {uc_out.a, uc_out.b, uc_out.c};
+    std::sort(lengths.begin(), lengths.end());
+    CHECK(lengths[0] == Catch::Approx(6.0f).epsilon(0.01));
+    CHECK(lengths[1] == Catch::Approx(8.0f).epsilon(0.01));
+    CHECK(lengths[2] == Catch::Approx(11.0f).epsilon(0.01));
+}
+
 TEST_CASE("FFTIndexer_CentredReferenceCell","[Indexing]") {
     // A reference cell is quoted the way a deposit states it, which for a centred lattice is the
     // conventional cell - and the indexer's candidates are Niggli-reduced primitive cells, whose