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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

image_analysis/scale_merge/TranslationalNCS.cpp

@@ -17,6 +17,7 @@
 #include <fftw3.h>
 
 #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<MergedReflection> &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<MergedReflectio
             // neighbouring one.
             double step = h_top > 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<double, 3>{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<double, 3> 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<double, NEIGHBOURS> 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;
             }

image_analysis/scale_merge/TranslationalNCS.h

@@ -119,6 +119,7 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
                                               const gemmi::UnitCell &cell,
                                               const gemmi::SpaceGroup *space_group = nullptr,
                                               int null_draws = 30,
-                                              uint32_t seed = 0);
+                                              uint32_t seed = 0,
+                                              size_t nthreads = 1);
 
 std::string TranslationalNCSToText(const TranslationalNCSResult &result);

rugnux/Rugnux.cpp

@@ -6776,7 +6776,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             // still carries (two thirds of them, and every axis step lands on one).
             TranslationalNCSResult pre_tncs;
             if (result.consensus_cell) {
-                pre_tncs = AnalyzeTranslationalNCS(sm.merged, gemmi::UnitCell(*result.consensus_cell));
+                pre_tncs = AnalyzeTranslationalNCS(sm.merged, gemmi::UnitCell(*result.consensus_cell), nullptr, 30, 0,
+                                                   config_.nthreads);
                 result.pre_promotion_tncs_measured = pre_tncs.measurable;
             }
             const char lattice_centring = end_msg.rotation_lattice_type.has_value()
@@ -7888,7 +7889,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             // it. Reporting only - no reflection and no merge depends on the outcome.
             if (result.consensus_cell)
                 result.tncs = AnalyzeTranslationalNCS(sm.merged, gemmi::UnitCell(*result.consensus_cell),
-                                                      twin_sg);
+                                                      twin_sg, 30, 0, config_.nthreads);
             if (result.tncs.detected)
                 logger.Warning("Translational pseudo-symmetry: Patterson off-origin peak {:.1f}% of the "
                                "origin at ({:.3f}, {:.3f}, {:.3f}), {:.1f} A - declare it to the "