From 51f51ec5bf0ec6f506baeaf89a6d6e97d7b6cb6c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 17 Sep 2026 13:03:06 +0200 Subject: [PATCH] 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) --- tests/IndexingUnitTest.cpp | 140 ++++++++++++++++++++++++++++++++++++- 1 file changed, 139 insertions(+), 1 deletion(-) diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index 0d299f331..28de9a065 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -2,6 +2,8 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include +#include #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 \ No newline at end of file +#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 NoiseCloud(size_t count, uint32_t seed, + const std::optional &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(rng()) / static_cast(std::mt19937::max()) - 1.0f; + }; + + const Coord n = plane_normal ? plane_normal->Normalize() : Coord(0, 0, 0); + + std::vector 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 FFTAlgorithms() { + std::vector 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 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 &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(indexed) / static_cast(spots.size()); + } + + // The lattice of the FFTIndexer test above, as spots. + std::vector LatticeSpots(const CrystalLattice &cl) { + std::vector 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(algorithm)); + auto indexer = MakeFFTIndexer(algorithm); + REQUIRE(indexer); + + const auto result = indexer->Run(vec); + logger.Info("algorithm {}: {} candidate lattices from noise", static_cast(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(algorithm)); + auto indexer = MakeFFTIndexer(algorithm); + REQUIRE(indexer); + + const auto noise_result = indexer->Run(noise); + logger.Info("algorithm {}: {} candidate lattices from coplanar noise", + static_cast(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 lengths = {static_cast(found.a), static_cast(found.b), + static_cast(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)); + } +}