tests: the FFT indexers on a frame of pure noise

The frame that works the indexer hardest is junk, not a crystal: a real lattice
prunes the FFT shortlist to four or five distinct directions, while noise leaves
dozens, which is what drives the length sort and the degenerate-plane fallback -
where a production broker segfaulted twice. There was no test in that regime.

Two cases, both over every FFT back-end the build has (GPU FFT under CUDA, CPU
FFTW always), so a non-CUDA job covers them as well:

- a cloud of 1500 reciprocal-space vectors of random direction and length, with
  no periodicity in it. Both back-ends return 70 candidate lattices from such a
  frame, which is the designed behaviour - Run() offers candidates and the caller
  scores them - so what the test pins is that none of them looks like a crystal:
  each takes about 1% of the cloud against the ~100% a real lattice takes, and
  the assertion allows 5%.
- the same noise flattened onto a tilted plane. A coplanar shortlist cannot close
  a cell, so this is the branch that goes looking for the missing row in a 3 deg
  cap (SearchCap), verified with a temporary probe to enter it on both back-ends
  and not to enter it on the isotropic frame. The SAME indexer object then has to
  index a clean lattice correctly, which is the regression test for the cap search
  putting the direction grid back.

The cloud is drawn from a fixed seed, and from the engine scaled by hand rather
than through uniform_real_distribution, whose output is not specified to be the
same in every standard library: a crash that needs one particular junk frame is
no use as a regression test if the frame is redrawn on each machine.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

tests/IndexingUnitTest.cpp

@@ -2,6 +2,8 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include <catch2/catch_all.hpp>
+#include <optional>
+#include <random>
 #include "../writer/HDF5Objects.h"
 #include "../image_analysis/indexing/IndexerFactory.h"
 #include "../image_analysis/indexing/PostIndexingRefinement.h"
@@ -595,4 +597,140 @@ TEST_CASE("FFBIDXIndexer_MultiLattice_TwoCrystals_BraggPrediction","[Indexing]")
 } */
 
 
-#endif
+#endif
+
+namespace {
+    // A frame of pure noise: reciprocal-space vectors of random direction and length between 50 A
+    // and 2.5 A, with no periodicity of any kind in them. The seed is fixed, so the frame is the
+    // same on every run and on every machine - a crash that needs a particular junk frame is no use
+    // as a regression test if the frame is drawn afresh each time. With a plane normal given, the
+    // cloud is flattened onto the plane through the origin, which is what a degenerate net looks
+    // like to the indexer.
+    std::vector<Coord> NoiseCloud(size_t count, uint32_t seed,
+                                  const std::optional<Coord> &plane_normal = {}) {
+        // The engine alone, scaled by hand: std::uniform_real_distribution is not specified to give
+        // the same numbers in every standard library, and a fixed frame is the whole point here.
+        std::mt19937 rng(seed);
+        const auto draw = [&rng]() {
+            return 2.0f * static_cast<float>(rng()) / static_cast<float>(std::mt19937::max()) - 1.0f;
+        };
+
+        const Coord n = plane_normal ? plane_normal->Normalize() : Coord(0, 0, 0);
+
+        std::vector<Coord> vec;
+        while (vec.size() < count) {
+            Coord q(draw(), draw(), draw());
+            if (plane_normal)
+                q = q - n * (q * n);
+            const float len = q.Length();
+            if ((len < 0.05f) || (len > 1.0f))   // a ball, not the cube the three draws fill
+                continue;
+            vec.push_back(q * 0.4f);             // |q| in [0.02, 0.4] 1/A
+        }
+        return vec;
+    }
+
+    std::vector<IndexingAlgorithmEnum> FFTAlgorithms() {
+        std::vector<IndexingAlgorithmEnum> ret;
+#ifdef JFJOCH_USE_CUDA
+        ret.push_back(IndexingAlgorithmEnum::FFT);
+#endif
+#ifdef JFJOCH_USE_FFTW
+        ret.push_back(IndexingAlgorithmEnum::FFTW);
+#endif
+        return ret;
+    }
+
+    std::unique_ptr<Indexer> MakeFFTIndexer(IndexingAlgorithmEnum algorithm) {
+        DiffractionExperiment experiment;
+        IndexingSettings settings;
+        settings.Algorithm(algorithm).FFT_MaxUnitCell_A(250.0).FFT_HighResolution_A(2 * M_PI / 3.0);
+        experiment.ImportIndexingSettings(settings);
+        return CreateIndexer(experiment);
+    }
+
+    // Fraction of the cloud that sits on the lattice, the same test the indexer's own scoring uses.
+    float IndexedFraction(const std::vector<Coord> &spots, const CrystalLattice &lattice, float tolerance) {
+        const float tol_sq = tolerance * tolerance;
+        size_t indexed = 0;
+        for (const auto &q: spots) {
+            const float h = q * lattice.Vec0() - std::round(q * lattice.Vec0());
+            const float k = q * lattice.Vec1() - std::round(q * lattice.Vec1());
+            const float l = q * lattice.Vec2() - std::round(q * lattice.Vec2());
+            if (h * h + k * k + l * l < tol_sq)
+                indexed++;
+        }
+        return static_cast<float>(indexed) / static_cast<float>(spots.size());
+    }
+
+    // The lattice of the FFTIndexer test above, as spots.
+    std::vector<Coord> LatticeSpots(const CrystalLattice &cl) {
+        std::vector<Coord> vec;
+        for (int h = -2; h < 10; h++)
+            for (int k = -5; k < 10; k++)
+                for (int l = -3; l < 10; l++)
+                    vec.push_back(h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar());
+        return vec;
+    }
+}
+
+TEST_CASE("FFTIndexer_NoiseFrame", "[Indexing]") {
+    // The frame that works the indexer hardest is not a clean crystal but junk. A real lattice
+    // prunes the FFT shortlist to four or five distinct directions; noise leaves dozens, which is
+    // what drives the length sort and the degenerate-plane fallback - where a production broker
+    // segfaulted. The run has to come back, and what it brings back has to be recognisably nothing.
+    const auto vec = NoiseCloud(1500, 20260917);
+
+    Logger logger("FFTIndexer_NoiseFrame");
+
+    for (const auto algorithm: FFTAlgorithms()) {
+        INFO("algorithm " << static_cast<int>(algorithm));
+        auto indexer = MakeFFTIndexer(algorithm);
+        REQUIRE(indexer);
+
+        const auto result = indexer->Run(vec);
+        logger.Info("algorithm {}: {} candidate lattices from noise", static_cast<int>(algorithm),
+                    result.lattice.size());
+
+        // Candidates off a junk frame are spurious by construction, and the indexer is allowed to
+        // offer them - the caller scores them. What must not happen is one that looks like a real
+        // crystal: a true lattice takes most of its frame, and these take about 1%.
+        for (const auto &lattice: result.lattice)
+            CHECK(IndexedFraction(vec, lattice, 0.1f) < 0.05f);
+    }
+}
+
+TEST_CASE("FFTIndexer_CoplanarNoiseFrameLeavesTheIndexerUsable", "[Indexing]") {
+    // Noise confined to one plane: a shortlist that is coplanar cannot close a cell whatever is done
+    // with it, so the indexer goes looking for the missing row inside a 3 deg cap (SearchCap), which
+    // swaps the whole search grid out and back. The indexer is not thrown away afterwards - the
+    // broker's pool hands the same object the next image - so the frame after a junk one must index
+    // exactly as it would have on its own.
+    const auto noise = NoiseCloud(1500, 20260918, Coord(0.3f, -0.5f, 0.81f));
+
+    const UnitCell uc(39, 45, 78, 90, 90, 90);
+    const auto lattice_spots = LatticeSpots(CrystalLattice(uc));
+    Logger logger("FFTIndexer_CoplanarNoiseFrame");
+
+    for (const auto algorithm: FFTAlgorithms()) {
+        INFO("algorithm " << static_cast<int>(algorithm));
+        auto indexer = MakeFFTIndexer(algorithm);
+        REQUIRE(indexer);
+
+        const auto noise_result = indexer->Run(noise);
+        logger.Info("algorithm {}: {} candidate lattices from coplanar noise",
+                    static_cast<int>(algorithm), noise_result.lattice.size());
+        for (const auto &lattice: noise_result.lattice)
+            CHECK(IndexedFraction(noise, lattice, 0.1f) < 0.05f);
+
+        const auto result = indexer->Run(lattice_spots);
+        REQUIRE(result.lattice.size() == 1);
+        const auto found = result.lattice[0].GetUnitCell();
+        std::array<float, 3> lengths = {static_cast<float>(found.a), static_cast<float>(found.b),
+                                        static_cast<float>(found.c)};
+        std::sort(lengths.begin(), lengths.end());
+        CHECK(lengths[0] == Catch::Approx(uc.a).epsilon(0.01));
+        CHECK(lengths[1] == Catch::Approx(uc.b).epsilon(0.01));
+        CHECK(lengths[2] == Catch::Approx(uc.c).epsilon(0.01));
+    }
+}