From ea383f9e83c122d40a18ee613be7d7c29c5d2987 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 23 Oct 2026 16:35:25 +0200 Subject: [PATCH] FFBIDXIndexerCPU: the fast-feedback indexer on the CPU, used where there is no GPU A port of the fast-feedback indexer's CUDA kernels to plain C++/Eigen: the same half-sphere and rotation samples, the same scores and the same 32 output cells for Refine(). The loops are reorganised so the inner one runs over samples and vectorises (GCC/Clang, no intrinsics), and the cell search separates the rotation about the candidate vector out of the per-spot terms. The vector refinement reproduces the library's warp sums - partial sums its shuffle reduction leaves in the other lanes carry into the next round - and its float arithmetic (FMA contraction, the Fibonacci points in float), because the candidates searched depend on it: on the library's 27 test spot lists and 4 synthetic stills the CPU and GPU best cells agree to ~1e-6 relative, the top-32 sets in 30 of 31 cases completely. FFBIDX now resolves to the GPU library where a GPU is visible and to the CPU port otherwise (no GPU, --no-gpu, or a build without CUDA); Auto takes FFBIDX for stills whenever a cell is known, also without a GPU, instead of FFTW. With a GPU nothing changes (a stills set and a rotation set reproduce the baseline exactly). Without a GPU the stills indexer pool gets one indexer per worker instead of eight, since the CPU indexer holds no card memory. Single core, synthetic still, loaded 5950X: ~8 ms per call at 30 spots, ~15 ms at 80, ~35 ms at 250; 16 threads ~1300/700/330 calls/s. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit 9a8eb8669d1be6c3ec16db6786280affd6830ee5) --- common/DiffractionExperiment.cpp | 9 +- frontend/src/components/IndexingSettings.tsx | 2 +- image_analysis/indexing/CMakeLists.txt | 7 + image_analysis/indexing/FFBIDXIndexerCPU.cpp | 431 ++++++++++++++++++ image_analysis/indexing/FFBIDXIndexerCPU.h | 30 ++ image_analysis/indexing/IndexerFactory.cpp | 9 +- image_analysis/indexing/IndexerThreadPool.cpp | 7 +- rugnux/Rugnux.cpp | 8 +- rugnux/Rugnux.h | 3 + rugnux/RugnuxStills.cpp | 2 +- rugnux/rugnux_cli.cpp | 12 +- tests/CMakeLists.txt | 12 + tests/FFBIDXIndexerCPUTest.cpp | 264 +++++++++++ viewer/widgets/JFJochViewerSettingsDock.cpp | 8 +- 14 files changed, 781 insertions(+), 23 deletions(-) create mode 100644 image_analysis/indexing/FFBIDXIndexerCPU.cpp create mode 100644 image_analysis/indexing/FFBIDXIndexerCPU.h create mode 100644 tests/FFBIDXIndexerCPUTest.cpp diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index a21b1ca1e..488a8f1fb 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1714,11 +1714,10 @@ IndexingAlgorithmEnum DiffractionExperiment::GetIndexingAlgorithm() const { return IndexingAlgorithmEnum::None; return IndexingAlgorithmEnum::FFBIDX; case IndexingAlgorithmEnum::Auto: - if (get_gpu_count() == 0) - return IndexingAlgorithmEnum::FFTW; - if (!cell) - return IndexingAlgorithmEnum::FFT; - return IndexingAlgorithmEnum::FFBIDX; + // FFBIDX runs on the GPU where there is one and on the CPU otherwise + if (cell) + return IndexingAlgorithmEnum::FFBIDX; + return get_gpu_count() == 0 ? IndexingAlgorithmEnum::FFTW : IndexingAlgorithmEnum::FFT; case IndexingAlgorithmEnum::FFT: // No silent switch to FFTW, too slow for live stills: a missing GPU turns indexing off return get_gpu_count() == 0 ? IndexingAlgorithmEnum::None : IndexingAlgorithmEnum::FFT; diff --git a/frontend/src/components/IndexingSettings.tsx b/frontend/src/components/IndexingSettings.tsx index 93dbbb849..adae9df31 100644 --- a/frontend/src/components/IndexingSettings.tsx +++ b/frontend/src/components/IndexingSettings.tsx @@ -104,7 +104,7 @@ function IndexingSettings({s: serverS, status}: MyProps) { value={s.algorithm} onChange={handleAlgorithmChange} > - FFBIDX (known unit cell) - GPU + FFBIDX (known unit cell) - GPU, or CPU without one FFT (guess unit cell) - GPU with CuFFT library FFT (guess unit cell) - CPU with FFTW library Auto (FFBIDX if unit cell provided, FFT otherwise) diff --git a/image_analysis/indexing/CMakeLists.txt b/image_analysis/indexing/CMakeLists.txt index 5f7451c55..212810ce4 100644 --- a/image_analysis/indexing/CMakeLists.txt +++ b/image_analysis/indexing/CMakeLists.txt @@ -69,6 +69,13 @@ ELSE() TARGET_LINK_LIBRARIES(JFJochIndexing Eigen3::Eigen) ENDIF() +# The CPU fast-feedback indexer is always built. Without errno and FP traps to preserve, GCC and Clang +# vectorise its scoring loops (sqrt, rint, min/max); neither option changes a computed value. +TARGET_SOURCES(JFJochIndexing PRIVATE FFBIDXIndexerCPU.cpp FFBIDXIndexerCPU.h) +IF (NOT MSVC) + SET_SOURCE_FILES_PROPERTIES(FFBIDXIndexerCPU.cpp PROPERTIES COMPILE_OPTIONS "-fno-math-errno;-fno-trapping-math") +ENDIF() + # FFTW (fftw3f) is always available via FetchContent -> the CPU FFT indexer is always built. TARGET_SOURCES(JFJochIndexing PRIVATE FFTIndexerCPU.cpp FFTIndexerCPU.h) TARGET_LINK_LIBRARIES(JFJochIndexing fftw3f) diff --git a/image_analysis/indexing/FFBIDXIndexerCPU.cpp b/image_analysis/indexing/FFBIDXIndexerCPU.cpp new file mode 100644 index 000000000..244c149a8 --- /dev/null +++ b/image_analysis/indexing/FFBIDXIndexerCPU.cpp @@ -0,0 +1,431 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// Ported from the CUDA kernels of the fast-feedback indexer, +// https://github.com/paulscherrerinstitute/fast-feedback-indexer (indexer/src/indexer_gpu.cu) +// (c) 2022 Paul Scherrer Institute, BSD-3-Clause +// +// The samples and the scores are the library's; the loops are reorganised for a CPU. Each block of +// samples is scored against one spot at a time, so the inner loop runs over samples and vectorises, +// and in the cell search the rotation about the candidate vector is separated out of the spot terms. + +#include "FFBIDXIndexerCPU.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "PostIndexingRefinement.h" + +namespace { + // The GPU library's runtime defaults (config_runtime in ffbidx/indexer.h) and rugnux's persistent + // configuration of it (FFBIDXIndexer.h). Keep them in step: the CPU vs GPU test compares the two. + constexpr int NUM_HALFSPHERE_POINTS = 32 * 1024; + constexpr int NUM_ANGLE_POINTS = 1024; // the library's own choice for 32k half-sphere points + constexpr int NUM_CANDIDATE_VECTORS = 32; + constexpr int NUM_OUTPUT_CELLS = 32; + constexpr int MIN_SPOTS = 6; // vector refinement + constexpr float TRIML = 0.001f; + constexpr float TRIMH = 0.3f; + constexpr float DELTA = 0.1f; + constexpr float DIST1 = 0.2f; + constexpr float DIST3 = 0.15f; + constexpr float LENGTH_THRESHOLD = 1e-9f; + + // Half-sphere directions are scored in blocks of this many, so a block's accumulators stay in L1 + constexpr int BLOCK = 1024; + static_assert(NUM_HALFSPHERE_POINTS % BLOCK == 0); + + // Every factor of a score's product is in [TRIML + DELTA, TRIMH + DELTA]; after 8 of them the + // product is still a normal float + constexpr int RENORMALISE_EVERY = 8; + + struct Tables { + std::array, 3> dir; // Fibonacci points on the half sphere z >= 0 + std::vector cos_alpha, sin_alpha; // rotation angles about a cell vector + }; + + // (a, rest) = a + x * y + rest, a Kahan sum, with the product fused into the first addition as + // nvcc compiles the library's + inline void KahanAdd(float &a, float &rest, float x, float y = 1.0f) { + const float s = std::fma(x, y, rest); + const float t = a; + a = t + s; + rest = s - (a - t); + } + + Tables MakeTables() { + constexpr double pi2 = 2.0 * std::numbers::pi; + const double dl = 3.0 - std::sqrt(5.0); + + Tables t; + for (auto &v: t.dir) + v.resize(NUM_HALFSPHERE_POINTS); + std::array dl2n; + for (int n = 0; n < 32; n++) + dl2n[n] = static_cast(std::fmod(dl * std::ldexp(1.0, n), pi2)); + for (int i = 0; i < NUM_HALFSPHERE_POINTS; i++) { + const float z = 1.0f - static_cast(i) / static_cast(NUM_HALFSPHERE_POINTS - 1); + const float r = std::sqrt(std::fma(-z, z, 1.0f)); + float l = 0, rest = 0; + for (int n = 0; n < 32; n++) + if (i & (1 << n)) + KahanAdd(l, rest, dl2n[n]); + t.dir[0][i] = std::cos(l) * r; + t.dir[1][i] = std::sin(l) * r; + t.dir[2][i] = z; + } + for (int q = 0; q < NUM_ANGLE_POINTS; q++) { + const double alpha = q * pi2 / NUM_ANGLE_POINTS; + t.cos_alpha.push_back(static_cast(std::cos(alpha))); + t.sin_alpha.push_back(static_cast(std::sin(alpha))); + } + return t; + } + + const Tables &GetTables() { + static const Tables t = MakeTables(); + return t; + } + + inline float Dist2Int(float x) { + return std::abs(x - std::rint(x)); + } + + // One spot's share of a sample's score: a factor of the product, and whether it counts as good + inline void AddSpot(float d, float good_limit, float &product, int &good) { + good += static_cast(d < good_limit); + product *= std::clamp(d, TRIML, TRIMH) + DELTA; + } + + // Moves the power of two of the product into the exponent, leaving the product in [1, 2). The + // library takes a log2 every 8 spots instead; this is the same sum, and it vectorises. + inline void Renormalise(float &product, int &exponent) { + const auto bits = std::bit_cast(product); + exponent += static_cast(bits >> 23) - 127; + product = std::bit_cast((bits & 0x007fffffu) | 0x3f800000u); + } + + // exp2(mean log2(trim(d) + delta)) - delta - (good + 1): lower is better. A NaN, from a degenerate + // candidate, goes last instead of breaking the sort. + inline float Score(float log2_sum, int good, int nspots) { + const float s = std::exp2(log2_sum / static_cast(nspots)) - DELTA - static_cast(good + 1); + return std::isnan(s) ? std::numeric_limits::infinity() : s; + } + + // Indices of the k lowest scores, best first + std::vector Best(const std::vector &score, int k) { + std::vector idx(score.size()); + std::iota(idx.begin(), idx.end(), 0); + auto lower = [&score](int a, int b) { return score[a] < score[b] || (score[a] == score[b] && a < b); }; + std::nth_element(idx.begin(), idx.begin() + k, idx.end(), lower); + idx.resize(k); + std::sort(idx.begin(), idx.end(), lower); + return idx; + } + + // Stage 1: every half-sphere direction as the direction of a cell vector of this length, scored on + // how close the spots' projections on it come to integers + std::vector FindCandidates(const Tables &t, const std::vector &spots, float length) { + const int nspots = static_cast(spots.size()); + std::vector score(NUM_HALFSPHERE_POINTS); + std::array product; + std::array exponent, good; + + for (int b = 0; b < NUM_HALFSPHERE_POINTS; b += BLOCK) { + const float *vx = t.dir[0].data() + b; + const float *vy = t.dir[1].data() + b; + const float *vz = t.dir[2].data() + b; + product.fill(1.0f); + exponent.fill(0); + good.fill(0); + for (int i = 0; i < nspots; i++) { + const Eigen::Vector3f s = length * spots[i]; + for (int k = 0; k < BLOCK; k++) + AddSpot(Dist2Int(s.x() * vx[k] + s.y() * vy[k] + s.z() * vz[k]), DIST1, product[k], good[k]); + if ((i + 1) % RENORMALISE_EVERY == 0) + for (int k = 0; k < BLOCK; k++) + Renormalise(product[k], exponent[k]); + } + for (int k = 0; k < BLOCK; k++) + score[b + k] = Score(static_cast(exponent[k]) + std::log2(product[k]), good[k], nspots); + } + return Best(score, NUM_CANDIDATE_VECTORS); + } + + // The library's dot product: rounding matters here, as a spot on the edge of a tolerance decides + // which spots the candidate is refined against + inline float Dot(const Eigen::Vector3f &a, const Eigen::Vector3f &b) { + return std::fma(a[0], b[0], std::fma(a[1], b[1], a[2] * b[2])); + } + + // A sum as the library forms it on a warp of 32 lanes: spot i is added on lane i % 32, and the + // lanes are reduced by shuffles that leave partial sums behind in every lane but the first. Those + // stay in the lanes and enter the next round's sum. Reproduced as it is, because the refined + // candidates - and so the cells searched - depend on it. + struct WarpSum { + std::array value{}, rest{}; + + void Add(size_t spot, float x, float y) { + KahanAdd(value[spot % 32], rest[spot % 32], x, y); + } + + float Reduce() { + std::array r{}; + for (int d = 1; d < 32; d *= 2) { + const auto old_value = value; + const auto old_r = r; + for (int l = 0; l < 32; l++) { + const int from = (l % (2 * d) < d) ? l + d : l; + r[l] = old_r[from]; + KahanAdd(value[l], r[l], old_value[from]); + } + } + return value[0]; + } + }; + + // x = A^-1 v for symmetric A given by its upper triangle (00 01 02 11 12 22), by cofactors, as + // the library solves it (products fused as nvcc compiles them) + Eigen::Vector3f SymSolve(const std::array &A, const std::array &v) { + auto cofactor = [](float a, float b, float c, float d) { return std::fma(a, b, -(c * d)); }; + const float co00 = cofactor(A[3], A[5], A[4], A[4]); + const float co01 = cofactor(A[1], A[5], A[2], A[4]); + const float co02 = cofactor(A[1], A[4], A[2], A[3]); + const float co11 = cofactor(A[0], A[5], A[2], A[2]); + const float co12 = cofactor(A[0], A[4], A[1], A[2]); + const float co22 = cofactor(A[0], A[3], A[1], A[1]); + // a0 * b0 + a1 * b1 + a2 * b2 + auto sum3 = [](float a0, float b0, float a1, float b1, float a2, float b2) { + float t = a0 * b0; + float rest = 0; + KahanAdd(t, rest, a1, b1); + KahanAdd(t, rest, a2, b2); + return t; + }; + const float r_det = 1.0f / sum3(A[0], co00, -A[1], co01, A[2], co02); + return {sum3(co00, v[0], -co01, v[1], co02, v[2]) * r_det, + sum3(-co01, v[0], co11, v[1], -co12, v[2]) * r_det, + sum3(co02, v[0], -co12, v[1], co22, v[2]) * r_det}; + } + + // Stage 2: a candidate refined by least squares against the spots it nearly indexes, in four rounds + // of tightening tolerance; the Miller indices stay those of the unrefined vector + Eigen::Vector3f RefineCandidate(const std::vector &spots, const Eigen::Vector3f &z0) { + std::array STS; + std::array STm; + Eigen::Vector3f z = z0; + for (const float beta: {0.1f, 0.05f, 0.025f, 0.01f}) { + int good = 0; + for (size_t i = 0; i < spots.size(); i++) { + const Eigen::Vector3f &s = spots[i]; + const float m = std::rint(Dot(s, z0)); + if (std::abs(Dot(s, z) - m) >= beta) + continue; + good++; + int e = 0; + for (int j = 0; j < 3; j++) { + STm[j].Add(i, s[j], m); + for (int k = j; k < 3; k++) + STS[e++].Add(i, s[j], s[k]); + } + } + if (good < MIN_SPOTS) + break; + std::array A; + std::array v; + for (int e = 0; e < 3; e++) + v[e] = STm[e].Reduce(); + for (int e = 0; e < 6; e++) + A[e] = STS[e].Reduce(); + z = SymSolve(A, v); + } + return z; + } + + // The input cell turned so its vector cell_vec lies along c, in the frame (x, y, z = c/|c|) the + // rotation about c is sampled in: a(alpha) = a_z z + a_xy (cos alpha x + sin alpha y), and b the + // same at alpha + delta. + struct CellFrame { + int cell_vec; + Eigen::Vector3f x, y, z; + float c_length; + float a_z, a_xy, b_z, b_xy; + float cos_delta, sin_delta; + + Eigen::Vector3f A(float cos_a, float sin_a) const { + return a_z * z + a_xy * (cos_a * x + sin_a * y); + } + + Eigen::Vector3f B(float cos_a, float sin_a) const { + const float c = cos_a * cos_delta - sin_a * sin_delta; + const float s = sin_a * cos_delta + cos_a * sin_delta; + return b_z * z + b_xy * (c * x + s * y); + } + }; + + CellFrame MakeFrame(const Eigen::Matrix3f &cell, int cell_vec, const Eigen::Vector3f &c) { + CellFrame f; + f.cell_vec = cell_vec; + f.c_length = c.norm(); + f.z = c / f.c_length; + + // Mirroring on the bisector of the input vector and z takes one onto the other + const Eigen::Vector3f v = cell.row(cell_vec); + Eigen::Vector3f a = cell.row((cell_vec + 1) % 3); + Eigen::Vector3f b = cell.row((cell_vec + 2) % 3); + Eigen::Vector3f t = v + v.norm() * f.z; + if (t.norm() >= std::numeric_limits::epsilon()) + t.normalize(); + else + t = f.z.cross(a).normalized(); + a = 2.0f * a.dot(t) * t - a; + b = 2.0f * b.dot(t) * t - b; + + f.a_z = a.dot(f.z); + f.x = (a - f.a_z * f.z).normalized(); + f.a_xy = a.dot(f.x); + f.y = f.z.cross(f.x); + + f.b_z = b.dot(f.z); + const Eigen::Vector3f b_dir = (b - f.b_z * f.z).normalized(); + f.b_xy = b.dot(b_dir); + f.cos_delta = b_dir.dot(f.x); + f.sin_delta = b_dir.dot(f.y); + return f; + } + + // Stage 3: the cell of one frame at every rotation angle, scored on how close each spot comes to + // a lattice point. The spot's projections on z, x and y are taken once; c does not turn at all. + void ScoreAngles(const Tables &t, const CellFrame &f, const std::vector &spots, float *score) { + const int nspots = static_cast(spots.size()); + const float *cos_a = t.cos_alpha.data(); + const float *sin_a = t.sin_alpha.data(); + std::array product; + std::array exponent, good; + product.fill(1.0f); + exponent.fill(0); + good.fill(0); + + for (int i = 0; i < nspots; i++) { + const float zs = f.z.dot(spots[i]); + const float xs = f.x.dot(spots[i]); + const float ys = f.y.dot(spots[i]); + const float dc = Dist2Int(f.c_length * zs); + const float a0 = f.a_z * zs; + const float a1 = f.a_xy * xs; + const float a2 = f.a_xy * ys; + const float b0 = f.b_z * zs; + const float b1 = f.b_xy * (f.cos_delta * xs + f.sin_delta * ys); + const float b2 = f.b_xy * (f.cos_delta * ys - f.sin_delta * xs); + for (int q = 0; q < NUM_ANGLE_POINTS; q++) { + const float da = Dist2Int(a0 + a1 * cos_a[q] + a2 * sin_a[q]); + const float db = Dist2Int(b0 + b1 * cos_a[q] + b2 * sin_a[q]); + AddSpot(std::sqrt(dc * dc + da * da + db * db), DIST3, product[q], good[q]); + } + if ((i + 1) % RENORMALISE_EVERY == 0) + for (int q = 0; q < NUM_ANGLE_POINTS; q++) + Renormalise(product[q], exponent[q]); + } + // The library starts this sum at 1, not 0; kept, so the scores agree with the GPU + for (int q = 0; q < NUM_ANGLE_POINTS; q++) + score[q] = Score(1.0f + static_cast(exponent[q]) + std::log2(product[q]), good[q], nspots); + } +} + +// Following the fast-feedback indexer (H.-C. Stadler, PSI) and TORO, Gasparotto et al. (2024) +// J. Appl. Cryst. 57, 931-944 +FFBIDXCells FFBIDXFindCells(const Eigen::Matrix3f &cell, const std::vector &coord, size_t nspots) { + const Tables &t = GetTables(); + + std::vector spots(nspots); + for (size_t i = 0; i < nspots; i++) + spots[i] = Eigen::Vector3f(coord[i].x, coord[i].y, coord[i].z); + + // Cell vectors of the same length share one set of candidates + std::array group{}; + std::vector group_length; + for (int v = 0; v < 3; v++) { + const float length = cell.row(v).norm(); + auto it = std::find_if(group_length.begin(), group_length.end(), + [length](float l) { return std::abs(l - length) < LENGTH_THRESHOLD; }); + group[v] = static_cast(it - group_length.begin()); + if (it == group_length.end()) + group_length.push_back(length); + } + + std::vector> candidates(group_length.size()); + for (size_t g = 0; g < group_length.size(); g++) { + for (const int i: FindCandidates(t, spots, group_length[g])) { + const Eigen::Vector3f dir(t.dir[0][i], t.dir[1][i], t.dir[2][i]); + candidates[g].push_back(RefineCandidate(spots, group_length[g] * dir)); + } + } + + // Every cell vector along every candidate of its length, both ways round + std::vector frames; + for (int v = 0; v < 3; v++) { + for (const auto &c: candidates[group[v]]) { + frames.push_back(MakeFrame(cell, v, c)); + frames.push_back(MakeFrame(cell, v, -c)); + } + } + + std::vector score(frames.size() * NUM_ANGLE_POINTS); + for (size_t i = 0; i < frames.size(); i++) + ScoreAngles(t, frames[i], spots, score.data() + i * NUM_ANGLE_POINTS); + + FFBIDXCells ret; + ret.cells.resize(3 * NUM_OUTPUT_CELLS, 3); + ret.scores.resize(NUM_OUTPUT_CELLS); + const auto best = Best(score, NUM_OUTPUT_CELLS); + for (int k = 0; k < NUM_OUTPUT_CELLS; k++) { + const CellFrame &f = frames[best[k] / NUM_ANGLE_POINTS]; + const int q = best[k] % NUM_ANGLE_POINTS; + ret.scores(k) = score[best[k]]; + ret.cells.row(3 * k + f.cell_vec) = f.c_length * f.z; + ret.cells.row(3 * k + (f.cell_vec + 1) % 3) = f.A(t.cos_alpha[q], t.sin_alpha[q]); + ret.cells.row(3 * k + (f.cell_vec + 2) % 3) = f.B(t.cos_alpha[q], t.sin_alpha[q]); + } + return ret; +} + +void FFBIDXIndexerCPU::SetupUnitCell(const std::optional &cell) { + if (!cell.has_value()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "FFBIDX requires unit cell"); + reference_unit_cell = cell; + + CrystalLattice l(cell.value()); + input_cell << l.Vec0().x, l.Vec0().y, l.Vec0().z, + l.Vec1().x, l.Vec1().y, l.Vec1().z, + l.Vec2().x, l.Vec2().y, l.Vec2().z; +} + +std::vector FFBIDXIndexerCPU::RunInternal(const std::vector &coord, size_t nspots) { + if (nspots > coord.size()) + nspots = coord.size(); + + if (nspots < viable_cell_min_spots) + return {}; + + auto found = FFBIDXFindCells(input_cell, coord, nspots); + + RefineParameters parameters{ + .viable_cell_min_spots = viable_cell_min_spots, + .dist_tolerance_vs_reference = dist_tolerance_vs_reference, + .reference_unit_cell = reference_unit_cell, + .min_length_A = 1, // doesn't matter + .max_length_A = 1000, // doesn't matter + .min_angle_deg = 30, + .max_angle_deg = 150, + .indexing_tolerance = indexing_tolerance + }; + + return Refine(coord, nspots, found.cells, found.scores, parameters); +} diff --git a/image_analysis/indexing/FFBIDXIndexerCPU.h b/image_analysis/indexing/FFBIDXIndexerCPU.h new file mode 100644 index 000000000..2136bbf15 --- /dev/null +++ b/image_analysis/indexing/FFBIDXIndexerCPU.h @@ -0,0 +1,30 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/Coord.h" +#include "Indexer.h" + +// The best cells of the fast-feedback indexer for one input cell, before Refine(): three rows per +// cell (the real-space cell vectors), scores ascending - lower is better. The layout of the GPU +// library's oCellM()/oScoreV(), so both feed Refine() the same way. +struct FFBIDXCells { + Eigen::MatrixX3 cells; + Eigen::VectorX scores; +}; + +// The fast-feedback indexer on the CPU: same samples, scores and output as the GPU library. +// cell holds the real-space cell vectors as rows; the first nspots of spots are used. +FFBIDXCells FFBIDXFindCells(const Eigen::Matrix3f &cell, const std::vector &spots, size_t nspots); + +class FFBIDXIndexerCPU : public Indexer { + Eigen::Matrix3f input_cell; + + void SetupUnitCell(const std::optional &cell) override; +public: + std::vector RunInternal(const std::vector &coord, size_t nspots) override; +}; diff --git a/image_analysis/indexing/IndexerFactory.cpp b/image_analysis/indexing/IndexerFactory.cpp index cf8adecb2..52be026d3 100644 --- a/image_analysis/indexing/IndexerFactory.cpp +++ b/image_analysis/indexing/IndexerFactory.cpp @@ -2,6 +2,8 @@ // SPDX-License-Identifier: GPL-3.0-only #include "IndexerFactory.h" +#include "FFBIDXIndexerCPU.h" +#include "../../common/CUDAWrapper.h" #ifdef JFJOCH_USE_CUDA #include "FFBIDXIndexer.h" @@ -17,9 +19,12 @@ std::unique_ptr CreateIndexer(const DiffractionExperiment& experiment) switch (experiment.GetIndexingAlgorithm()) { case IndexingAlgorithmEnum::FFBIDX: #ifdef JFJOCH_USE_CUDA - indexer = std::make_unique(); - indexer->Setup(experiment); + if (get_gpu_count() > 0) + indexer = std::make_unique(); #endif + if (!indexer) + indexer = std::make_unique(); + indexer->Setup(experiment); break; case IndexingAlgorithmEnum::FFTW: #ifdef JFJOCH_USE_FFTW diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp index bd8e2696d..ea4150583 100644 --- a/image_analysis/indexing/IndexerThreadPool.cpp +++ b/image_analysis/indexing/IndexerThreadPool.cpp @@ -11,6 +11,7 @@ #include "../common/CUDAWrapper.h" #include "../common/Logger.h" +#include "FFBIDXIndexerCPU.h" #ifdef JFJOCH_USE_CUDA #include "FFBIDXIndexer.h" @@ -173,6 +174,8 @@ static std::unique_ptr MakeIndexer(IndexingAlgorithmEnum algorithm, con return std::make_unique(); } #endif + if (algorithm == IndexingAlgorithmEnum::FFBIDX) + return std::make_unique(); #ifdef JFJOCH_USE_FFTW if (algorithm == IndexingAlgorithmEnum::FFTW) return std::make_unique(settings); @@ -446,7 +449,7 @@ IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, co const bool servable = (algorithm == IndexingAlgorithmEnum::FFT && have_gpu && (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFT)) || - (algorithm == IndexingAlgorithmEnum::FFBIDX && have_gpu && + (algorithm == IndexingAlgorithmEnum::FFBIDX && (requested == IndexingAlgorithmEnum::Auto || requested == IndexingAlgorithmEnum::FFBIDX)) || (algorithm == IndexingAlgorithmEnum::FFTW && fftw_built && ((requested == IndexingAlgorithmEnum::Auto && !have_gpu) || requested == IndexingAlgorithmEnum::FFTW)); @@ -457,8 +460,6 @@ IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, co "always uses the GPU FFT indexer here; select FFT or Auto, or run FFTW on a CPU-only node."; else if (algorithm == IndexingAlgorithmEnum::FFT && !have_gpu) msg = "FFT is the GPU indexer but no GPU is available. Select FFTW or Auto for CPU indexing."; - else if (algorithm == IndexingAlgorithmEnum::FFBIDX && !have_gpu) - msg = "FFBIDX is a GPU indexer but no GPU is available. Select FFTW or Auto for CPU indexing."; else if (algorithm == IndexingAlgorithmEnum::FFTW) msg = "FFTW (CPU) indexing was requested but this build has no FFTW indexer."; else diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a1401382a..e7aa66f59 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -593,6 +593,12 @@ int Rugnux::FirstPassIndexerThreads() const { return gpu_work_beside_ && gpu_takes_turns() ? 1 : 2; } +int Rugnux::StillsIndexerThreads() const { + if (experiment_.GetIndexingAlgorithm() == IndexingAlgorithmEnum::FFBIDX && get_gpu_count() == 0) + return std::clamp(config_.nthreads, 1, 64); + return std::clamp(config_.nthreads, 1, 8); +} + int Rugnux::FirstPassSpotWorkers() const { const int gpus = get_gpu_count(); if (gpus == 0) @@ -1054,7 +1060,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Measured on a 16 Mpx serial set, de-novo FFT: 83 frames/s with 4, 125 with 8, 133 with 16 - // by eight the card is what limits it. Each GPU FFT indexer holds ~0.5 GB of the card // (sixteen of them took the run to 12.7 GB), so eight rather than one per worker. - indexing_settings.IndexingThreads(std::clamp(config_.nthreads, 1, 8)); + indexing_settings.IndexingThreads(StillsIndexerThreads()); } indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A(); indexer_pool = std::make_unique(indexing_settings, IndexerConstruction::OnFirstUse, diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 20d350226..9b539dfc2 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -1489,6 +1489,9 @@ class Rugnux { // An FFT indexer holds 0.3-0.6 GB of the card, and which one indexes a frame does not change what it // finds. [[nodiscard]] int FirstPassIndexerThreads() const; + // The indexers of a stills pool: eight where each holds card memory; one per worker for the CPU + // fast-feedback indexer, which holds next to nothing. + [[nodiscard]] int StillsIndexerThreads() const; // The first pass's spot-finding workers, each with GPU engines of its own where there is a card. [[nodiscard]] int FirstPassSpotWorkers() const; // What a probe pass that only indexes holds of a card at its most: its first pass's workers and FFT diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index 26d4f3d08..29872e86b 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -841,7 +841,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // As many indexer threads as the image pass has (RunPipeline): the hypotheses are indexed side by side, and // with the default four the workers waited on the pool. IndexingSettings pool_settings = experiment_.GetIndexingSettings(); - pool_settings.IndexingThreads(std::clamp(config_.nthreads, 1, 8)); + pool_settings.IndexingThreads(StillsIndexerThreads()); IndexerThreadPool pool(pool_settings, IndexerConstruction::OnFirstUse, gpu_work_beside_); // A hypothesis is indexed once: the searches below meet the same one more than once - the decided centre diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index c6c8468be..9a9f45e01 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -138,7 +138,7 @@ void print_usage() { std::cout << " -N, --threads Number of threads (default: every CPU the run may use - its affinity and cgroup quota - of which the " "per-image loop uses at most 16 per GPU; an explicit value is used as given)" << std::endl; std::cout << " --no-gpu Run on the CPU only, as on a machine without a GPU: CUDA is not " - "touched and every GPU step takes its CPU path (-X Auto and FFT index with FFTW)" << std::endl; + "touched and every GPU step takes its CPU path (FFT indexes with FFTW, FFBIDX with its CPU implementation)" << std::endl; std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; std::cout << " -e, --end-image End image number, exclusive: -s 0 -e 100 processes images 0-99 (default: all)" << std::endl; std::cout << " -t, --stride Image stride (default: 1)" << std::endl; @@ -197,7 +197,9 @@ void print_usage() { std::cout << " --prepass-fraction Integrate the first fraction f of the sweep, but at least 120 deg of it, in the two-pass rotation geometry pre-pass (0.05-1, default 0.5; 1 = the whole sweep); pass 2 reads every frame. The pre-scan and the first-pass indexing sample the same part of the sweep" << std::endl; std::cout << " -X, --indexing-algorithm Indexing algorithm (Auto|FFT|cuFFT|FFTW|FFBIDX|None, default Auto). 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 and FFBIDX (known cell) need a GPU" << std::endl; + "CPU (FFTW); cuFFT and FFTW name one implementation each, and cuFFT needs a GPU. FFBIDX is the known-cell " + "fast-feedback indexer for stills, likewise on the GPU where there is one and otherwise on the CPU; " + "Auto picks it for stills whenever a cell is known" << std::endl; std::cout << " -S, --space-group Space group number (96) or symbol (P43212) - for indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; std::cout << " --fft-min-unit-cell Shortest cell axis the FFT search accepts, in A (default: 10). A candidate with a shorter axis is discarded, but a de-novo rotation run also tries a second first-pass hypothesis with the floor at 5 A, so a small-molecule cell is indexed without this flag; -C lowers it on its own to cover the cell given" << std::endl; @@ -1789,10 +1791,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } } - if ((indexing_algorithm == IndexingAlgorithmEnum::FFT || indexing_algorithm == IndexingAlgorithmEnum::FFBIDX) - && get_gpu_count() == 0) { - logger.Error("-X {} runs on the GPU only, and {}; -X FFT or Auto index on the CPU", - indexing_algorithm == IndexingAlgorithmEnum::FFT ? "cuFFT" : "FFBIDX", + if (indexing_algorithm == IndexingAlgorithmEnum::FFT && get_gpu_count() == 0) { + logger.Error("-X cuFFT runs on the GPU only, and {}; -X FFT or Auto index on the CPU", no_gpu ? "--no-gpu is given" : "no GPU is visible"); exit(EXIT_FAILURE); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 70744345f..cecfe9c5a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,6 +35,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL LoadReferenceMtzFreeFlagsTest.cpp LoadReferenceSfCifTest.cpp RotationIndexerTest.cpp + FFBIDXIndexerCPUTest.cpp ) target_link_libraries(jfjoch_portable_test Catch2WithMain JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter JFJochCommon) @@ -44,6 +45,13 @@ IF (JFJOCH_PORTABLE_ONLY AND JFJOCH_CUDA_AVAILABLE AND TARGET CUDA::cufft_static SET_TARGET_PROPERTIES(jfjoch_portable_test PROPERTIES CUDA_RESOLVE_DEVICE_SYMBOLS ON) ENDIF() +# The CPU vs GPU test of the fast-feedback indexer reads the library's own test spot lists +IF (JFJOCH_CUDA_AVAILABLE) + FetchContent_GetProperties(fast-indexer SOURCE_DIR FFBIDX_SOURCE_DIR) + SET(FFBIDX_SIMPLE_DATA_DIR "${FFBIDX_SOURCE_DIR}/data/simple/files") + TARGET_COMPILE_DEFINITIONS(jfjoch_portable_test PRIVATE FFBIDX_SIMPLE_DATA_DIR="${FFBIDX_SIMPLE_DATA_DIR}") +ENDIF() + IF (JFJOCH_PORTABLE_ONLY) RETURN() ENDIF() @@ -67,6 +75,7 @@ ADD_EXECUTABLE(jfjoch_test JFJochReceiverIntegrationTest.cpp AcquisitionCountersTest.cpp IndexingUnitTest.cpp + FFBIDXIndexerCPUTest.cpp ThreadSafeFIFOSetTest.cpp ReceiveRingBufferTest.cpp FPGANetworkTest.cpp @@ -215,6 +224,9 @@ target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJ JFJochRugnuxOutput JFJochRugnux JFJochModelWriter JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview jfjoch_xds_plugin ZLIB::ZLIB httplib::httplib) target_include_directories(jfjoch_test PRIVATE .) +IF (JFJOCH_CUDA_AVAILABLE) + TARGET_COMPILE_DEFINITIONS(jfjoch_test PRIVATE FFBIDX_SIMPLE_DATA_DIR="${FFBIDX_SIMPLE_DATA_DIR}") +ENDIF() find_package(Threads REQUIRED) diff --git a/tests/FFBIDXIndexerCPUTest.cpp b/tests/FFBIDXIndexerCPUTest.cpp new file mode 100644 index 000000000..3c91c1f1e --- /dev/null +++ b/tests/FFBIDXIndexerCPUTest.cpp @@ -0,0 +1,264 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "../common/CUDAWrapper.h" +#include "../common/DiffractionExperiment.h" +#include "../image_analysis/indexing/FFBIDXIndexerCPU.h" + +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/indexing/FFBIDXIndexer.h" +#endif + +namespace { + // A still of a test cell in a fixed random orientation: the reciprocal lattice points to 2 A that lie + // within 0.001 1/A of the Ewald sphere of a 1 A beam along z, and then random spots in the same + // resolution range. The seed is fixed, and the engine scaled by hand, so the frame is the same + // with every standard library. + std::vector SyntheticStill(const UnitCell &uc, uint32_t seed, size_t noise_spots) { + std::mt19937 rng(seed); + auto draw = [&rng]() { + return 2.0f * static_cast(rng()) / static_cast(std::mt19937::max()) - 1.0f; + }; + const Eigen::Matrix3f rot = Eigen::Quaternionf(draw(), draw(), draw(), draw()).normalized().toRotationMatrix(); + + const CrystalLattice l(uc); + auto rotate = [&rot](const Coord &v) { + const Eigen::Vector3f x = rot * Eigen::Vector3f(v.x, v.y, v.z); + return Coord(x.x(), x.y(), x.z()); + }; + const Coord astar = rotate(l.Astar()), bstar = rotate(l.Bstar()), cstar = rotate(l.Cstar()); + + const Coord k0(0, 0, 1.0f); + std::vector spots; + for (int h = -40; h <= 40; h++) + for (int k = -40; k <= 40; k++) + for (int m = -40; m <= 40; m++) { + const Coord r = h * astar + k * bstar + m * cstar; + if (r.Length() > 0.5f || r.Length() < 0.02f) + continue; + if (std::abs((r + k0).Length() - 1.0f) < 0.001f) + spots.push_back(r); + } + while (noise_spots > 0) { + const Coord r(0.5f * draw(), 0.5f * draw(), 0.5f * draw()); + if (r.Length() > 0.5f || r.Length() < 0.02f) + continue; + spots.push_back(r); + noise_spots--; + } + return spots; + } + + // Fraction of spots within 0.1 of a lattice point (in Miller index units) + float IndexedFraction(const std::vector &spots, const CrystalLattice &l) { + size_t n = 0; + for (const auto &s: spots) { + const float h = s * l.Vec0(), k = s * l.Vec1(), m = s * l.Vec2(); + const float dh = h - std::round(h), dk = k - std::round(k), dm = m - std::round(m); + if (dh * dh + dk * dk + dm * dm < 0.01f) + n++; + } + return static_cast(n) / static_cast(spots.size()); + } + + bool SameLengths(const UnitCell &a, const UnitCell &b, float rel_tol) { + std::array x = {a.a, a.b, a.c}, y = {b.a, b.b, b.c}; + std::sort(x.begin(), x.end()); + std::sort(y.begin(), y.end()); + for (int i = 0; i < 3; i++) + if (std::abs(x[i] - y[i]) > rel_tol * y[i]) + return false; + return true; + } + + // The cell vectors as rows + Eigen::Matrix3f CellMatrix(const UnitCell &uc) { + const CrystalLattice l(uc); + Eigen::Matrix3f m; + m << l.Vec0().x, l.Vec0().y, l.Vec0().z, l.Vec1().x, l.Vec1().y, l.Vec1().z, + l.Vec2().x, l.Vec2().y, l.Vec2().z; + return m; + } + + const std::vector &TestCells() { + static const std::vector cells = { + {30, 40, 50, 90, 90, 90}, + {80, 80, 90, 90, 90, 120}, + {40, 45, 80, 90, 82.5, 90}, + {60, 60, 60, 90, 90, 90} + }; + return cells; + } +} + +TEST_CASE("Portable_FFBIDXIndexerCPU_SyntheticStill", "[portable][Indexing]") { + for (size_t i = 0; i < TestCells().size(); i++) { + const UnitCell &uc = TestCells()[i]; + const auto lattice_spots = SyntheticStill(uc, 1234 + i, 0); + auto spots = SyntheticStill(uc, 1234 + i, lattice_spots.size() / 5); + INFO("cell " << i << ", " << lattice_spots.size() << " lattice spots"); + REQUIRE(lattice_spots.size() > 30); + + DiffractionExperiment experiment; + experiment.SetUnitCell(uc); + FFBIDXIndexerCPU indexer; + indexer.Setup(experiment); + const auto result = indexer.Run(spots); + REQUIRE(!result.lattice.empty()); + CHECK(SameLengths(result.lattice[0].GetUnitCell(), uc, 0.01f)); + CHECK(IndexedFraction(lattice_spots, result.lattice[0]) > 0.95f); + } +} + +// Time per call of the CPU indexer, alone and with Refine(), on one thread and on 16 at once +TEST_CASE("FFBIDXIndexerCPU_Benchmark", "[.][ffbidx_bench]") { + const UnitCell uc(40, 45, 80, 90, 82.5, 90); + DiffractionExperiment experiment; + experiment.SetUnitCell(uc); + + for (const size_t n: {30, 80, 250}) { + auto spots = SyntheticStill(uc, 99, 0); + const auto lattice_spots = spots.size(); + spots = SyntheticStill(uc, 99, n > lattice_spots ? n - lattice_spots : 0); + spots.resize(n); + + const int calls = 20; + auto t0 = std::chrono::steady_clock::now(); + for (int i = 0; i < calls; i++) + FFBIDXFindCells(CellMatrix(uc), spots, n); + const double find_ms = std::chrono::duration(std::chrono::steady_clock::now() - t0).count() / calls; + + FFBIDXIndexerCPU indexer; + indexer.Setup(experiment); + t0 = std::chrono::steady_clock::now(); + for (int i = 0; i < calls; i++) + indexer.Run(spots); + const double run_ms = std::chrono::duration(std::chrono::steady_clock::now() - t0).count() / calls; + + const int threads = 16; + t0 = std::chrono::steady_clock::now(); + std::vector workers; + for (int t = 0; t < threads; t++) + workers.emplace_back([&] { + FFBIDXIndexerCPU w; + w.Setup(experiment); + for (int i = 0; i < calls; i++) + w.Run(spots); + }); + for (auto &w: workers) + w.join(); + const double wall_s = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + + WARN(n << " spots: FindCells " << find_ms << " ms, with Refine " << run_ms << " ms; " + << threads << " threads " << threads * calls / wall_s << " calls/s"); + } +} + +#ifdef JFJOCH_USE_CUDA + +namespace { + // The GPU library's raw output, before Refine(), with the configuration FFBIDXIndexer uses + FFBIDXCells RunGPU(const Eigen::Matrix3f &cell, const std::vector &spots) { + fast_feedback::config_runtime crt{}; + fast_feedback::config_persistent cpers{ + .max_output_cells = 32, + .max_input_cells = 1, + .max_spots = static_cast(spots.size()), + .num_candidate_vectors = 32, + .redundant_computations = true, + }; + fast_feedback::refine::indexer indexer{cpers, crt}; + indexer.iCellM() = cell; + for (size_t i = 0; i < spots.size(); i++) { + indexer.spotX(i) = spots[i].x; + indexer.spotY(i) = spots[i].y; + indexer.spotZ(i) = spots[i].z; + } + indexer.index(1, spots.size()); + return {indexer.oCellM(), indexer.oScoreV()}; + } + + // Cell (rows) and spots of one of the library's test files: the first line is the cell, then a + // spot per line; '#' starts a comment + std::pair> ReadSimpleData(const std::string &path) { + std::ifstream f(path); + REQUIRE(f.good()); + std::vector v; + std::string line; + while (std::getline(f, line)) { + if (line.empty() || line[0] == '#') + continue; + std::istringstream s(line); + float x; + while (s >> x) + v.push_back(x); + } + REQUIRE(v.size() >= 12); + Eigen::Matrix3f cell; + cell << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]; + std::vector spots; + for (size_t i = 9; i + 2 < v.size(); i += 3) + spots.emplace_back(v[i], v[i + 1], v[i + 2]); + return {cell, spots}; + } + + float RelativeDifference(const Eigen::Matrix3f &a, const Eigen::Matrix3f &b) { + return (a - b).norm() / b.norm(); + } + + // CPU and GPU agree when the best scores agree and the CPU's best cell is one of the GPU's best + // cells - the samples are the same, so a matching cell agrees to rounding, and two cells within + // rounding of each other in score may come out in either order + void CompareWithGPU(const std::string &name, const Eigen::Matrix3f &cell, const std::vector &spots) { + const auto gpu = RunGPU(cell, spots); + const auto cpu = FFBIDXFindCells(cell, spots, spots.size()); + + float best_match = INFINITY; + for (int k = 0; k < gpu.scores.size(); k++) + best_match = std::min(best_match, RelativeDifference(cpu.cells.block<3, 3>(0, 0), gpu.cells.block<3, 3>(3 * k, 0))); + + int top_matched = 0; + for (int j = 0; j < cpu.scores.size(); j++) + for (int k = 0; k < gpu.scores.size(); k++) + if (RelativeDifference(cpu.cells.block<3, 3>(3 * j, 0), gpu.cells.block<3, 3>(3 * k, 0)) < 1e-3f) { + top_matched++; + break; + } + + UNSCOPED_INFO(name << ": " << spots.size() << " spots, score CPU " << cpu.scores(0) << " GPU " << gpu.scores(0) + << ", best cell to GPU " << best_match << ", " << top_matched << "/32 cells matched"); + CHECK(std::abs(cpu.scores(0) - gpu.scores(0)) <= 1e-3f * std::abs(gpu.scores(0))); + CHECK(best_match < 1e-3f); + CHECK(top_matched >= 24); + } +} + +TEST_CASE("FFBIDXIndexerCPU_MatchesGPU", "[Indexing][GPU]") { + if (get_gpu_count() == 0) + SKIP("no GPU"); + + for (int i = 0; i < 9; i++) { + for (const std::string kind: {"local", "peakfinder8", "radial"}) { + const std::string name = "image" + std::to_string(i) + "_" + kind + ".txt"; + const auto [cell, spots] = ReadSimpleData(std::string(FFBIDX_SIMPLE_DATA_DIR) + "/" + name); + CompareWithGPU(name, cell, spots); + } + } + + for (size_t i = 0; i < TestCells().size(); i++) { + const auto spots = SyntheticStill(TestCells()[i], 1234 + i, 10); + CompareWithGPU("synthetic " + std::to_string(i), CellMatrix(TestCells()[i]), spots); + } +} + +#endif diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index b953a1d5d..0b66f75ca 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -371,7 +371,7 @@ QWidget *JFJochViewerSettingsDock::BuildMXPage() { algo_->setSizeAdjustPolicy(QComboBox::AdjustToMinimumContentsLengthWithIcon); algo_->setMinimumContentsLength(12); algo_->addItem("Auto", static_cast(IndexingAlgorithmEnum::Auto)); - algo_->addItem("FFBIDX (GPU, known cell)", static_cast(IndexingAlgorithmEnum::FFBIDX)); + algo_->addItem("FFBIDX (known cell)", static_cast(IndexingAlgorithmEnum::FFBIDX)); algo_->addItem("FFT (GPU, de-novo)", static_cast(IndexingAlgorithmEnum::FFT)); algo_->addItem("FFTW (CPU, de-novo)", static_cast(IndexingAlgorithmEnum::FFTW)); algo_->addItem("None", static_cast(IndexingAlgorithmEnum::None)); @@ -765,14 +765,14 @@ void JFJochViewerSettingsDock::UpdateAlgorithmDescription() { QString text; switch (a) { case IndexingAlgorithmEnum::FFBIDX: - text = "GPU, needs a known cell — best for sparse serial stills"; break; + text = "needs a known cell, GPU or CPU — best for sparse serial stills"; break; case IndexingAlgorithmEnum::FFT: text = "GPU, de-novo — best for strong rotation data"; break; case IndexingAlgorithmEnum::FFTW: text = "CPU, de-novo — no GPU needed"; break; case IndexingAlgorithmEnum::Auto: { - const IndexingAlgorithmEnum r = !gpu ? IndexingAlgorithmEnum::FFTW - : (cell_known ? IndexingAlgorithmEnum::FFBIDX : IndexingAlgorithmEnum::FFT); + const IndexingAlgorithmEnum r = cell_known ? IndexingAlgorithmEnum::FFBIDX + : (gpu ? IndexingAlgorithmEnum::FFT : IndexingAlgorithmEnum::FFTW); const QString rn = r == IndexingAlgorithmEnum::FFBIDX ? "FFBIDX" : r == IndexingAlgorithmEnum::FFT ? "FFT" : "FFTW"; text = QString("resolves to %1 here (%2, cell %3)")