From c244d384b057d886e718ecfec091dee9d7bada23 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 26 Sep 2026 18:08:36 +0200 Subject: [PATCH] tNCS: evaluate the class-vector ladder's neighbours in parallel Where a pseudo-translation is detected, the class vector is re-refined up a ladder on the whole data by a greedy sweep over the 26 neighbours of the current vector, each a cosine correlation over up to every acentric reflection - seconds of one thread on the main path of the run (3.4 s and 1.5 s for the two calls on a P3_121 16M set). The neighbours of the current vector are now evaluated together; the first improvement in the sweep's own order is taken, and the neighbours after it are evaluated again about the new vector - exactly the moves the serial sweep makes. Each correlation is still one serial sum, so the result is bit-identical. AnalyzeTranslationalNCS takes the thread count. md5-identical output and report; 47.1 -> 44.4 s on that set. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- .../scale_merge/TranslationalNCS.cpp | 51 +++++++++++-------- image_analysis/scale_merge/TranslationalNCS.h | 3 +- rugnux/Rugnux.cpp | 5 +- 3 files changed, 36 insertions(+), 23 deletions(-) diff --git a/image_analysis/scale_merge/TranslationalNCS.cpp b/image_analysis/scale_merge/TranslationalNCS.cpp index 5c8ba4779..ed2176762 100644 --- a/image_analysis/scale_merge/TranslationalNCS.cpp +++ b/image_analysis/scale_merge/TranslationalNCS.cpp @@ -17,6 +17,7 @@ #include #include "../../common/FFTWPlannerLock.h" +#include "../../common/ParallelFor.h" namespace { // The Patterson is measured in a band, not over the whole dataset: below ~20 A the bulk solvent @@ -151,7 +152,7 @@ namespace { TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup *space_group, - int null_draws, uint32_t seed) { + int null_draws, uint32_t seed, size_t nthreads) { TranslationalNCSResult result; if (cell.a <= 0.0 || cell.b <= 0.0 || cell.c <= 0.0) { result.refusal = "no unit cell"; @@ -677,28 +678,38 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector 0 ? 0.25 / h_top : std::max({step0[0], step0[1], step0[2]}); double best = CosineCorrelation(cls_hkl, cls_e2, u, n); + // All 26 neighbours, not the 6 along the axes: the phase is h.u, so a vector whose class + // comes mostly from one combination of the axes - h/3 + 2k/3, say - sits in a valley that + // runs diagonally, and a search that only moves one component at a time stalls on its + // wall (measured: 0.005 short in two components at once, a quarter of a turn by |h| = 60). + constexpr int NEIGHBOURS = 26; + const auto neighbour = [&](int j) { + const int k = j < 13 ? j : j + 1; // the 27 offsets in dx, dy, dz order, centre skipped + const int dx = k / 9 - 1, dy = k / 3 % 3 - 1, dz = k % 3 - 1; + return std::array{u[0] + dx * step, u[1] + dy * step, u[2] + dz * step}; + }; for (int round = 0; round < 8; ++round) { for (bool improved = true; improved;) { improved = false; - // All 26 neighbours, not the 6 along the axes: the phase is h.u, so a vector - // whose class comes mostly from one combination of the axes - h/3 + 2k/3, say - - // sits in a valley that runs diagonally, and a search that only moves one - // component at a time stalls on its wall (measured: 0.005 short in two - // components at once, a quarter of a turn by |h| = 60). - for (int dx = -1; dx <= 1; ++dx) - for (int dy = -1; dy <= 1; ++dy) - for (int dz = -1; dz <= 1; ++dz) { - if (dx == 0 && dy == 0 && dz == 0) - continue; - const std::array v{u[0] + dx * step, u[1] + dy * step, - u[2] + dz * step}; - const double f = CosineCorrelation(cls_hkl, cls_e2, v, n); - if (f > best + 1e-9) { - best = f; - u = v; - improved = true; - } - } + // The sweep takes each neighbour of the current vector in turn and moves to it + // when it improves, so the neighbours after a move are those of the new vector. + // They are evaluated together, the first improvement in sweep order is taken, + // and the rest are evaluated again about it - the serial sweep's moves exactly. + std::array f{}; + for (int j = 0; j < NEIGHBOURS;) { + ParallelFor(NEIGHBOURS - j, nthreads, [&](int i) { + f[j + i] = CosineCorrelation(cls_hkl, cls_e2, neighbour(j + i), n); + }); + int k = j; + while (k < NEIGHBOURS && !(f[k] > best + 1e-9)) + ++k; + if (k == NEIGHBOURS) + break; + best = f[k]; + u = neighbour(k); + improved = true; + j = k + 1; + } } step *= 0.5; } diff --git a/image_analysis/scale_merge/TranslationalNCS.h b/image_analysis/scale_merge/TranslationalNCS.h index dea828cdd..4b8c3f6b6 100644 --- a/image_analysis/scale_merge/TranslationalNCS.h +++ b/image_analysis/scale_merge/TranslationalNCS.h @@ -119,6 +119,7 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector