From e5442e7a0745aa49b5af12e48aba406d263935e1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 24 Aug 2026 21:19:59 +0200 Subject: [PATCH] Read the cell surface's twenty bytes, and give each candidate group a thread Three unrelated costs in the tail, each measured before and after. The cell surface fits and scores by reaching into the fulls for sixteen bytes of an eighty-byte record, twenty-two times over. That is four times the traffic of the data it uses, and it was the whole of the cost: the arithmetic never was. One fused compaction of the intensity, sigma, correction, cell term and group, built in the passes that were already being made, and every later pass walks the compact array instead. Fit accumulation falls to a third, scoring to a third. What is left is the term build and the scatter, not the fit. The estimator is untouched. The space-group search hashed a reflection key per observation per candidate. The orbit representative is now interned to a dense index when the orbits are built, so the two tests that follow index an array. Both candidate loops also run a thread per candidate - the point-group loop, and the space-group loop, which was the larger of the two by far: a full pass over the merge with three absence tests and a map insert per reflection, once for every candidate group. The operator cache is filled by the serial pass that precedes them, so each candidate still sees the same operators, does its own arithmetic unchanged, and appends in the same order. The search runs in a third of the time. Writing the reflections was seventeen stream insertions per row for a quarter of a million rows. The rows are formatted in parallel blocks and written in order, in the same widths and precisions as before: the mmCIF in a tenth of the time, the hkl in a quarter. Faster on twelve of twelve matched pairs, eight per cent on the sum of minima, and region timers account for the wall clock to within five per cent. Every reflection file over the whole rotation test set is byte-identical, in both passes and all three formats. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_016NNnL26LAvruQ9eLUUWvrJ --- image_analysis/WriteReflections.cpp | 135 +++++++++++++----- image_analysis/WriteReflections.h | 11 +- .../scale_merge/RotationScaleMerge.cpp | 93 ++++++++---- .../scale_merge/SearchSpaceGroup.cpp | 88 ++++++++---- rugnux/Rugnux.cpp | 3 +- rugnux/rugnux_cli.cpp | 2 +- 6 files changed, 235 insertions(+), 97 deletions(-) diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 73ddedce..20012b5e 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -5,6 +5,7 @@ #include "scale_merge/Merge.h" #include "scale_merge/HKLKey.h" #include "scale_merge/TwinningAnalysis.h" +#include "../common/ParallelFor.h" #include #include @@ -150,7 +151,8 @@ void WriteMmcifReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::string &filename) { + const std::string &filename, + size_t nthreads) { std::ofstream out(filename); if (!out) @@ -332,25 +334,50 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_refln.status_free\n"; out << "_refln.status\n"; - for (const auto& r : reflections) { - out << std::setw(5) << r.h << " " - << std::setw(5) << r.k << " " - << std::setw(5) << r.l << " " - << std::setw(14) << Fmt(r.I, 4) << " " - << std::setw(14) << Fmt(r.sigma, 4) << " " - << std::setw(14) << Fmt(r.I_plus, 4) << " " - << std::setw(14) << Fmt(r.sigma_plus, 4) << " " - << std::setw(14) << Fmt(r.I_minus, 4) << " " - << std::setw(14) << Fmt(r.sigma_minus, 4) << " " - << std::setw(14) << Fmt(r.F, 4) << " " - << std::setw(14) << Fmt(r.sigmaF, 4) << " " - << std::setw(14) << Fmt(r.F_plus, 4) << " " - << std::setw(14) << Fmt(r.sigmaF_plus, 4) << " " - << std::setw(14) << Fmt(r.F_minus, 4) << " " - << std::setw(14) << Fmt(r.sigmaF_minus, 4) << " " - << (r.rfree_flag ? 1 : 0) << " " - << "o" // 'o' = observed - << "\n"; + // One row per unique reflection, twelve formatted floats each - tens of megabytes on a crowded + // crystal, and the largest single-threaded stretch left in a run. Nothing about a row depends on + // any other, so each worker formats its own block into its own string and the blocks go to the + // file in order. The columns are written exactly as the stream wrote them: the same "%.4f" (or + // "?" where the value is not finite), right-aligned in the same width. + { + const size_t nrow = reflections.size(); + const size_t nw = std::max(nthreads, 1); + const int nch = static_cast(ThreadsForWork(nrow, nw, 4096)); + std::vector block(nch); + ParallelChunks(nch, nw, [&](int tlo, int thi) { + for (int t = tlo; t < thi; ++t) { + const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch; + std::string &s = block[t]; + s.reserve((hi - lo) * 208); + const auto column = [&s](const std::string &v, size_t width) { + if (v.size() < width) s.append(width - v.size(), ' '); + s.append(v); + s.push_back(' '); + }; + for (size_t i = lo; i < hi; ++i) { + const auto &r = reflections[i]; + column(std::to_string(r.h), 5); + column(std::to_string(r.k), 5); + column(std::to_string(r.l), 5); + column(Fmt(r.I, 4), 14); + column(Fmt(r.sigma, 4), 14); + column(Fmt(r.I_plus, 4), 14); + column(Fmt(r.sigma_plus, 4), 14); + column(Fmt(r.I_minus, 4), 14); + column(Fmt(r.sigma_minus, 4), 14); + column(Fmt(r.F, 4), 14); + column(Fmt(r.sigmaF, 4), 14); + column(Fmt(r.F_plus, 4), 14); + column(Fmt(r.sigmaF_plus, 4), 14); + column(Fmt(r.F_minus, 4), 14); + column(Fmt(r.sigmaF_minus, 4), 14); + s.push_back(r.rfree_flag ? '1' : '0'); + s.append(" o\n"); // 'o' = observed + } + } + }); + for (const std::string &s : block) + out.write(s.data(), static_cast(s.size())); } out << "#\n"; @@ -435,7 +462,8 @@ void WriteMtzReflections(const std::vector &reflections, void WriteShelxHklReflections(const std::vector &reflections, const DiffractionExperiment &experiment, - const std::string &filename) { + const std::string &filename, + size_t nthreads) { bool has_anom = true; const std::vector rows = BuildMergedRows(reflections, experiment, has_anom); @@ -459,20 +487,53 @@ void WriteShelxHklReflections(const std::vector &reflections, std::ofstream out(filename); if (!out) throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename); - out << std::fixed << std::setprecision(2); - const auto emit = [&out, scale](int h, int k, int l, float I, float sigma) { - out << std::setw(4) << h << std::setw(4) << k << std::setw(4) << l - << std::setw(8) << scale * I << std::setw(8) << scale * sigma << "\n"; + // Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker + // blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in + // the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write. + const auto column = [](std::string &s, const std::string &v, size_t width) { + if (v.size() < width) s.append(width - v.size(), ' '); + s.append(v); }; - for (const auto& r : rows) { - const bool plus = usable(r.Ip, r.sIp); - const bool minus = usable(r.Im, r.sIm); - if (plus) emit(r.h, r.k, r.l, r.Ip, r.sIp); - if (minus) emit(-r.h, -r.k, -r.l, r.Im, r.sIm); - if (!plus && !minus && usable(r.Imean, r.sImean)) - emit(r.h, r.k, r.l, r.Imean, r.sImean); + const auto num2 = [](double v) { + char b[64]; + const int n = std::snprintf(b, sizeof b, "%.2f", v); + return std::string(b, static_cast(std::clamp(n, 0, static_cast(sizeof b) - 1))); + }; + const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) { + column(s, std::to_string(h), 4); + column(s, std::to_string(k), 4); + column(s, std::to_string(l), 4); + column(s, num2(scale * I), 8); + column(s, num2(scale * sigma), 8); + s.push_back('\n'); + }; + { + const size_t nrow = rows.size(); + const size_t nw = std::max(nthreads, 1); + const int nch = static_cast(ThreadsForWork(nrow, nw, 4096)); + std::vector block(nch); + ParallelChunks(nch, nw, [&](int tlo, int thi) { + for (int t = tlo; t < thi; ++t) { + const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch; + std::string &s = block[t]; + s.reserve((hi - lo) * 2 * 29); + for (size_t i = lo; i < hi; ++i) { + const auto &r = rows[i]; + const bool plus = usable(r.Ip, r.sIp); + const bool minus = usable(r.Im, r.sIm); + if (plus) emit(s, r.h, r.k, r.l, r.Ip, r.sIp); + if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm); + if (!plus && !minus && usable(r.Imean, r.sImean)) + emit(s, r.h, r.k, r.l, r.Imean, r.sImean); + } + } + }); + for (const std::string &s : block) + out.write(s.data(), static_cast(s.size())); } - emit(0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker + std::string tail; + emit(tail, 0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker + out.write(tail.data(), static_cast(tail.size())); out.close(); } @@ -482,11 +543,13 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::string &filename) { + const std::string &filename, + size_t nthreads) { // Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as // the SHELXC / ANODE substructure-solution input). WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); - WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, filename + ".cif"); - WriteShelxHklReflections(reflections, experiment, filename + ".hkl"); + WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, + filename + ".cif", nthreads); + WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 3cf42c64..fb3f86c5 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -24,13 +24,15 @@ struct ErrorModelReport { std::string b; }; +// nthreads: workers for the per-reflection row formatting, which is the bulk of the file. void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::string &filename); + const std::string &filename, + size_t nthreads); void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, @@ -38,9 +40,11 @@ void WriteMtzReflections(const std::vector &reflections, const std::string &filename); // SHELX HKLF-4 text file (h k l I sigma(I), Bijvoet mates separate) for SHELXC / ANODE. +// nthreads: workers for the per-reflection row formatting, as for the mmCIF. void WriteShelxHklReflections(const std::vector &reflections, const DiffractionExperiment &experiment, - const std::string &filename); + const std::string &filename, + size_t nthreads); void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, @@ -48,4 +52,5 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::string &filename); \ No newline at end of file + const std::string &filename, + size_t nthreads); \ No newline at end of file diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e1198609..34f3e437 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1977,16 +1977,58 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int auto usable = [&](const Obs &o) { return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && o.partiality >= min_partiality; }; - // The observations any of the passes below can use, in fulls order, split by frame parity. Every - // pass then walks one list instead of re-testing (and re-reading) the whole of fulls - the - // cross-validated halves cost half a pass rather than a whole one with a skip. The lists keep the - // fulls order, so each sum below is formed from exactly the same terms in exactly the same order. + // The observations any of the passes below can use, in fulls order, with the four fields those + // passes read copied out beside them. Around twenty passes follow, and each one used to reach back + // into `fulls` for sixteen bytes of an eighty-byte Obs - so every pass paid four times the memory + // traffic of the data it touched, over hundreds of megabytes. Copying once turns them into + // sequential walks of a compact array. The copy keeps fulls order, so each sum below is formed + // from exactly the same terms in exactly the same order. + struct Term { float I, sigma, corr; int32_t cell, group; }; + std::vector term; + std::vector term_parity; // frame parity, read only by the group-ordered copy below + // Which terms each pass walks, as positions in `term`. The cross-validated halves then cost half a + // pass rather than a whole one with a skip; the full fit walks the array in order. std::vector idx_all, idx_even, idx_odd; - idx_all.reserve(fulls.size()); - for (size_t i = 0; i < fulls.size(); ++i) { - if (!usable(fulls[i]) || cell[i] < 0) continue; - idx_all.push_back(static_cast(i)); - (fulls[i].frame & 1 ? idx_odd : idx_even).push_back(static_cast(i)); + { + // Count what each chunk keeps, and how much of that sits on an even frame, so the second pass + // can write straight into its own span of every list. Chunk boundaries depend on nothing but + // the length and the worker count, so the lists come out in fulls order on every run. + const size_t nf = fulls.size(); + const int nch = static_cast(ThreadsForWork(nf, nthreads)); + auto cpart = [nf, nch](int t) { return nf * static_cast(t) / static_cast(nch); }; + std::vector off(nch + 1, 0), eoff(nch + 1, 0); + ParallelChunks(nch, nthreads, [&](int tlo, int thi) { + for (int t = tlo; t < thi; ++t) { + size_t keep = 0, even = 0; + for (size_t i = cpart(t); i < cpart(t + 1); ++i) { + if (!usable(fulls[i]) || cell[i] < 0) continue; + ++keep; + if (!(fulls[i].frame & 1)) ++even; + } + off[t + 1] = keep; eoff[t + 1] = even; + } + }); + for (int t = 0; t < nch; ++t) { off[t + 1] += off[t]; eoff[t + 1] += eoff[t]; } + term.resize(off[nch]); + term_parity.resize(off[nch]); + idx_all.resize(off[nch]); + idx_even.resize(eoff[nch]); + idx_odd.resize(off[nch] - eoff[nch]); + ParallelChunks(nch, nthreads, [&](int tlo, int thi) { + for (int t = tlo; t < thi; ++t) { + size_t w = off[t], e = eoff[t], d = off[t] - eoff[t]; + for (size_t i = cpart(t); i < cpart(t + 1); ++i) { + const Obs &o = fulls[i]; + if (!usable(o) || cell[i] < 0) continue; + const uint8_t parity = static_cast(o.frame & 1); + term[w] = {o.I, o.sigma, o.corr, cell[i], o.group}; + term_parity[w] = parity; + idx_all[w] = static_cast(w); + (parity ? idx_odd[d++] : idx_even[e++]) = static_cast(w); + ++w; + } + } + }); } auto subset = [&](int parity) -> const std::vector & { return parity < 0 ? idx_all : (parity ? idx_odd : idx_even); @@ -2000,22 +2042,21 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int // thread count, so a sum split this way gives the same answer on every run. auto part = [nt](size_t n, int t) { return n * static_cast(t) / static_cast(nt); }; - // The same observations once more, in ASU-GROUP order, with the fields the reference pass reads - // copied out beside them. That pass sums per group, so group order is the only layout that lets a - // thread own a whole group's sum - and it still adds that group's terms in fulls order, because the - // counting sort is stable (`fill` walks idx_all, which is ascending). Copying the fields out is what - // makes it worth doing: the alternative is chasing every term back through `fulls` at random over a - // few hundred megabytes of Obs, which is the locality this pass used to stay serial to keep. + // The same terms once more, in ASU-GROUP order. The reference pass sums per group, so group order + // is the only layout that lets a thread own a whole group's sum - and it still adds that group's + // terms in fulls order, because the counting sort is stable (it walks `term`, which is in fulls + // order). A separate copy rather than a permutation: chasing every term through an index at random + // over a few hundred megabytes is the locality this pass used to stay serial to keep. struct RefTerm { float I, sigma, corr; int32_t cell; uint8_t parity; }; std::vector g_start(n_groups + 1, 0); - for (const int32_t i : idx_all) ++g_start[fulls[i].group + 1]; + for (const Term &t : term) ++g_start[t.group + 1]; for (int g = 0; g < n_groups; ++g) g_start[g + 1] += g_start[g]; - std::vector gterm(idx_all.size()); + std::vector gterm(term.size()); { std::vector fill(g_start.begin(), g_start.end() - 1); - for (const int32_t i : idx_all) { - const Obs &o = fulls[i]; - gterm[fill[o.group]++] = {o.I, o.sigma, o.corr, cell[i], static_cast(o.frame & 1)}; + for (size_t k = 0; k < term.size(); ++k) { + const Term &t = term[k]; + gterm[fill[t.group]++] = {t.I, t.sigma, t.corr, t.cell, term_parity[k]}; } } @@ -2072,14 +2113,13 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int std::fill(xcross.begin(), xcross.end(), 0.0); std::fill(xref2.begin(), xref2.end(), 0.0); for (size_t k = part(sel.size(), t); k < part(sel.size(), t + 1); ++k) { - const int32_t i = sel[k]; - const Obs &o = fulls[i]; + const Term &o = term[sel[k]]; if (sw[o.group] <= 0.0) continue; - const double Iref = swI[o.group] / sw[o.group], a = A[cell[i]]; + const double Iref = swI[o.group] / sw[o.group], a = A[o.cell]; const double Is = static_cast(o.I) * o.corr * a, sc = static_cast(o.sigma) * o.corr * a; if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0) || !(sc > 0.0)) continue; const double w = 1.0 / (sc * sc); - xcross[cell[i]] += w * Is * Iref; xref2[cell[i]] += w * Iref * Iref; + xcross[o.cell] += w * Is * Iref; xref2[o.cell] += w * Iref * Iref; } } }); @@ -2111,10 +2151,9 @@ void RotationScaleMerge::ApplyCellSurface(const std::vector &cell, int for (int t = tlo; t < thi; ++t) { double num = 0.0, den = 0.0; for (size_t k = part(sel.size(), t); k < part(sel.size(), t + 1); ++k) { - const int32_t i = sel[k]; - const Obs &o = fulls[i]; + const Term &o = term[sel[k]]; if (sw[o.group] <= 0.0) continue; - const double a = A[cell[i]], Is = static_cast(o.I) * o.corr * a; + const double a = A[o.cell], Is = static_cast(o.I) * o.corr * a; const double Iref = swI[o.group] / sw[o.group]; if (!std::isfinite(Iref) || Iref <= 0.0) continue; num += std::abs(Is - Iref); den += Iref; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index e58abb91..21b3fe61 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -459,12 +459,16 @@ SearchSpaceGroupResult SearchSpaceGroup( // quantities below, and both need exactly the same mapping. Build it once per point group. struct Acc { double sw = 0.0, swI = 0.0; int n = 0; }; struct Orbits { - std::unordered_map grp; // representative -> inverse-variance accumulator - std::vector rep; // observation -> its representative + std::vector acc; // one inverse-variance accumulator per orbit, in first-seen order + std::vector orbit; // observation -> its orbit, -1 for one neither quantity below uses }; auto build_orbits = [&](const std::vector& rotations) -> Orbits { Orbits orb; - orb.rep.resize(n); + orb.orbit.assign(n, -1); + // The representative is interned to a dense index right here, so the two quantities below + // index an array rather than hashing a key per observation - each of them is a pass over the + // whole merge, and the lookup was the larger half of both. Same orbits, same order, same sums. + std::unordered_map orbit_id; for (size_t i = 0; i < n; ++i) { if (!pass_cc[i] || !(Sigma[i] > 0.0)) continue; @@ -475,8 +479,11 @@ SearchSpaceGroupResult SearchSpaceGroup( if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l)) best = k2; } - orb.rep[i] = best; - auto& g = orb.grp[best]; + const auto [it, fresh] = orbit_id.emplace(best, static_cast(orb.acc.size())); + if (fresh) + orb.acc.emplace_back(); + orb.orbit[i] = it->second; + Acc& g = orb.acc[it->second]; const double w = 1.0 / (Sigma[i] * Sigma[i]); g.sw += w; g.swI += w * I[i]; g.n += 1; } @@ -490,20 +497,18 @@ SearchSpaceGroupResult SearchSpaceGroup( // - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over // the present (pass_cc) reflections. auto chi2_under = [&](const Orbits& orb) -> double { - const auto& grp = orb.grp; - const auto& rep = orb.rep; double chi2 = 0.0; long dof = 0; for (size_t i = 0; i < n; ++i) { - if (!pass_cc[i] || !(Sigma[i] > 0.0)) + if (orb.orbit[i] < 0) continue; - const auto& g = grp.at(rep[i]); + const Acc& g = orb.acc[orb.orbit[i]]; if (g.n < 2) continue; const double mean = g.swI / g.sw, dev = I[i] - mean; chi2 += dev * dev / (Sigma[i] * Sigma[i]); } - for (const auto& [k, g] : grp) + for (const Acc& g : orb.acc) if (g.n >= 2) dof += g.n - 1; return dof > 0 ? chi2 / static_cast(dof) : std::numeric_limits::quiet_NaN(); @@ -523,13 +528,11 @@ SearchSpaceGroupResult SearchSpaceGroup( // convention; converting them to any other silently squares the ratios and makes the absolute // floor a-dependent, on a quantity that has no a. Leave it alone. auto merge_systematic_b = [&](const Orbits& orb) -> double { - const auto& grp = orb.grp; - const auto& rep = orb.rep; std::vector> obs; // I, sigma, deviation-from-orbit-mean for (size_t i = 0; i < n; ++i) { - if (!pass_cc[i] || !(Sigma[i] > 0.0)) + if (orb.orbit[i] < 0) continue; - const auto& g = grp.at(rep[i]); + const Acc& g = orb.acc[orb.orbit[i]]; if (g.n < 2) continue; obs.push_back({I[i], Sigma[i], I[i] - g.swI / g.sw}); @@ -559,20 +562,37 @@ SearchSpaceGroupResult SearchSpaceGroup( std::string refused_pg_hm, refused_why; std::vector pg_cands; double chi2_ref = std::numeric_limits::infinity(); + // Which point groups their own operators confirm. Serial: this fills the operator cache, and it + // is only a handful of cache lookups per group either way. + std::vector confirmed; + std::vector confirmed_cc; for (const auto& pg : point_groups) { const auto [present, min_class_cc] = point_group_present(pg.rotations); if (!present) continue; - double ch = std::numeric_limits::quiet_NaN(); - double be = 0.0; - if (!pg.rotations.empty()) { - const Orbits orb = build_orbits(pg.rotations); - ch = chi2_under(orb); - be = merge_systematic_b(orb); - } - pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, min_class_cc, ch, be}); - if (!pg.rotations.empty() && std::isfinite(ch)) - chi2_ref = std::min(chi2_ref, ch); + confirmed.push_back(&pg); + confirmed_cc.push_back(min_class_cc); + } + // Merging under a candidate costs three passes over the whole merge, and there are a handful of + // candidates, each reading nothing but the shared reflection arrays - so give each one a thread. + // Every candidate's sums are still formed on one thread, over the same reflections in the same + // order, so the numbers this hands back do not depend on the split. + std::vector cand_chi2(confirmed.size(), std::numeric_limits::quiet_NaN()); + std::vector cand_b(confirmed.size(), 0.0); + ParallelFor(static_cast(confirmed.size()), + std::min(opt.nthreads, confirmed.size()), [&](int i) { + if (confirmed[i]->rotations.empty()) + return; + const Orbits orb = build_orbits(confirmed[i]->rotations); + cand_chi2[i] = chi2_under(orb); + cand_b[i] = merge_systematic_b(orb); + }); + for (size_t i = 0; i < confirmed.size(); ++i) { + const PointGroupInfo& pg = *confirmed[i]; + pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, confirmed_cc[i], + cand_chi2[i], cand_b[i]}); + if (!pg.rotations.empty() && std::isfinite(cand_chi2[i])) + chi2_ref = std::min(chi2_ref, cand_chi2[i]); } // Choose the largest point group that is both operator-confirmed AND self-consistent (its merge @@ -779,10 +799,18 @@ SearchSpaceGroupResult SearchSpaceGroup( return result; } - for (const auto& sg : gemmi::spacegroup_tables::main) { - if (!sg.is_sohncke() || !sg.is_reference_setting() || RotationSetOf(sg) != best_pg->rotation_set) - continue; - + // The candidate space groups of the chosen point group. Scoring one is a pass over the whole + // merge with three absence tests per reflection, and there are up to a dozen of them; they read + // nothing but the shared reflection arrays, so give each one a thread. Every candidate's own pass + // is unchanged and they are appended in table order, so the ranking below sees what it saw before. + std::vector sg_cands; + for (const auto& sg : gemmi::spacegroup_tables::main) + if (sg.is_sohncke() && sg.is_reference_setting() && RotationSetOf(sg) == best_pg->rotation_set) + sg_cands.push_back(&sg); + std::vector sg_scored(sg_cands.size()); + ParallelFor(static_cast(sg_cands.size()), + std::min(opt.nthreads, sg_cands.size()), [&](int ci) { + const gemmi::SpaceGroup& sg = *sg_cands[ci]; const gemmi::GroupOps gops = sg.operations(); SpaceGroupCandidateScore s{.space_group = sg}; double absent_sum = 0, present_sum = 0; @@ -948,8 +976,10 @@ SearchSpaceGroupResult SearchSpaceGroup( const bool screw_ok = screw_absent == 0 || screw_violations <= opt.max_absent_violation_fraction * screw_absent; s.consistent = centering_ok && screw_ok; + sg_scored[ci] = std::move(s); + }); + for (auto& s : sg_scored) result.candidates.push_back(std::move(s)); - } // A candidate is eligible when its absences are confirmed and there are enough of them to // trust (the symmorphic group, with no absences, is always eligible as the fallback). Rank diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d0cff4c9..57bbd2ff 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2862,7 +2862,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(), result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()}; WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, - em_report, result.twinning, config_.output_prefix); + em_report, result.twinning, config_.output_prefix, + static_cast(std::max(1, config_.nthreads))); // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS // comparison. The offline self-scaling result is otherwise not exposed (process.h5's // per-image arrays are only filled on the online per-image path). Sourced from the diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 8fcad8ef..63c445a4 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1527,7 +1527,7 @@ static int RunRugnux(int argc, char **argv) { error_model_a > 0 ? fmt::format("{:.3f}", error_model_a) : std::string(), error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, - em_report, twinning, output_prefix); + em_report, twinning, output_prefix, static_cast(nthreads)); } if (!output_prefix.empty() && !model_pdb.empty()) {