diff --git a/common/IndexingSettings.cpp b/common/IndexingSettings.cpp index 2de46bcb..1cbbdd22 100644 --- a/common/IndexingSettings.cpp +++ b/common/IndexingSettings.cpp @@ -141,6 +141,17 @@ IndexingSettings &IndexingSettings::IndexingThreads(int64_t input) { return *this; } +int64_t IndexingSettings::GetRefineThreads() const { + return refine_threads; +} + +IndexingSettings &IndexingSettings::RefineThreads(int64_t input) { + check_min("Candidate-cell refinement thread count", input, 1); + check_max("Candidate-cell refinement thread count", input, 64); + refine_threads = input; + return *this; +} + IndexingSettings &IndexingSettings::UnitCellDistTolerance(float input) { check_min("Relative unit cell distance tolerance vs. reference", input, 0.0001); check_max("Relative unit cell distance tolerance vs. reference", input, 0.2001); diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 258e15ba..3fef55d9 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -23,6 +23,10 @@ class IndexingSettings { float unit_cell_dist_tolerance_vs_reference = 0.05; // relative static constexpr float unit_cell_angle_tolerance_deg = 5.0; // degree int64_t indexing_threads = 4; + // Threads splitting the candidate-cell refinement WITHIN one indexer call. 1 (the default) is the + // right answer whenever indexers already run one per image across all workers; it is raised only + // where few indexer threads exist and cores would otherwise sit idle. + int64_t refine_threads = 1; int64_t viable_cell_min_spots = 9; int64_t max_extra_lattices = 2; @@ -48,6 +52,7 @@ public: IndexingSettings& FFT_HighResolution_A(float input); IndexingSettings& Tolerance(float input); IndexingSettings& IndexingThreads(int64_t input); + IndexingSettings& RefineThreads(int64_t input); IndexingSettings& UnitCellDistTolerance(float input); IndexingSettings& GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum input); IndexingSettings& IndexIceRings(bool input); @@ -68,6 +73,7 @@ public: [[nodiscard]] float GetFFT_MinAngle_deg() const; [[nodiscard]] float GetFFT_MaxAngle_deg() const; [[nodiscard]] int64_t GetIndexingThreads() const; + [[nodiscard]] int64_t GetRefineThreads() const; [[nodiscard]] float GetUnitCellDistTolerance() const; [[nodiscard]] float GetUnitCellAngleTolerance_deg() const; [[nodiscard]] bool GetIndexIceRings() const; diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 610fdb3a..334445c5 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include #include @@ -39,12 +40,15 @@ std::string CurrentDateISO() { } /// Format a double with given decimal places; returns "?" for non-finite. +/// snprintf rather than an ostringstream: the reflection loop below calls this twelve times per +/// reflection, and building a stream (and its locale) per call dominated the time spent writing a +/// merged file. Same digits - both go through the C locale's %.*f. std::string Fmt(double val, int decimals = 4) { if (!std::isfinite(val)) return "?"; - std::ostringstream ss; - ss << std::fixed << std::setprecision(decimals) << val; - return ss.str(); + char buf[512]; + const int n = std::snprintf(buf, sizeof(buf), "%.*f", decimals, val); + return std::string(buf, n); } /// Quote a CIF string value; returns "?" for empty. diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 87b33ff8..8596507a 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -12,6 +12,7 @@ FFTIndexer::FFTIndexer(const IndexingSettings &settings) min_angle_deg(settings.GetFFT_MinAngle_deg()), max_angle_deg(settings.GetFFT_MaxAngle_deg()), nDirections(settings.GetFFT_NumVectors()), + refine_threads(static_cast(settings.GetRefineThreads())), result_fft(nDirections) { // Reciprocal-magnitude histogram in one_over_d = 1/d units (the internal convention - @@ -332,7 +333,8 @@ std::vector FFTIndexer::ReduceAndRefine(const std::vector .max_length_A = max_length_A, .min_angle_deg = min_angle_deg, .max_angle_deg = max_angle_deg, - .indexing_tolerance = indexing_tolerance + .indexing_tolerance = indexing_tolerance, + .refine_threads = refine_threads }; return Refine(coord, nspots, oCell, scores, parameters); diff --git a/image_analysis/indexing/FFTIndexer.h b/image_analysis/indexing/FFTIndexer.h index d699d4d2..b5061e3d 100644 --- a/image_analysis/indexing/FFTIndexer.h +++ b/image_analysis/indexing/FFTIndexer.h @@ -21,6 +21,7 @@ protected: const float max_angle_deg; const int nDirections; + const unsigned refine_threads; float histogram_spacing; int64_t histogram_size; std::vector direction_vectors; diff --git a/image_analysis/indexing/PostIndexingRefinement.cpp b/image_analysis/indexing/PostIndexingRefinement.cpp index e65f6486..7875accf 100644 --- a/image_analysis/indexing/PostIndexingRefinement.cpp +++ b/image_analysis/indexing/PostIndexingRefinement.cpp @@ -5,6 +5,7 @@ #include "PostIndexingRefinement.h" #include +#include namespace { struct config_ifssr final { @@ -145,7 +146,25 @@ std::vector Refine(const std::vector &in_spots, .min_spots = static_cast(p.viable_cell_min_spots) }; - RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr); + // Candidate cells refine independently - a block touches only its own scores(j) and cells rows, and + // holds its own scratch - so splitting them across threads gives the same numbers as one thread. + // Only worth it where few indexer threads run (the rotation first pass uses two, one per scheme, + // and leaves the rest of the machine idle); refine_threads stays 1 everywhere else. + const unsigned ncells = static_cast(scores.rows()); + const unsigned nblocks = std::max(1u, std::min(p.refine_threads, ncells)); + if (nblocks == 1) { + RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr); + } else { + std::vector workers; + workers.reserve(nblocks - 1); + for (unsigned b = 1; b < nblocks; b++) + workers.emplace_back([&, b] { + RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr, b, nblocks); + }); + RefineCandidateCells(spots.topRows(nspots), oCell, scores, cifssr, 0, nblocks); + for (auto &w : workers) + w.join(); + } std::vector candidates; diff --git a/image_analysis/indexing/PostIndexingRefinement.h b/image_analysis/indexing/PostIndexingRefinement.h index 205783f5..eb1d25d9 100644 --- a/image_analysis/indexing/PostIndexingRefinement.h +++ b/image_analysis/indexing/PostIndexingRefinement.h @@ -30,6 +30,7 @@ struct RefineParameters { float min_angle_deg; float max_angle_deg; float indexing_tolerance; + unsigned refine_threads = 1; }; std::vector Refine(const std::vector &in_spots, diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 431a3634..bb3f9c1c 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -375,16 +375,17 @@ SearchSpaceGroupResult SearchSpaceGroup( holohedry = HolohedryRotationSet(opt.lattice_system.value()); const auto point_groups = EnumeratePointGroups(holohedry); - // Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry - // equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces - // non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is - // more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well - // - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over - // the present (pass_cc) reflections. - auto chi2_under = [&](const std::vector& rotations) -> double { - struct Acc { double sw = 0.0, swI = 0.0; int n = 0; }; - std::unordered_map grp; - std::vector rep(n); + // Mapping every observation onto its symmetry representative under a candidate's rotations - one + // apply_to_hkl + Canonicalize per observation per operator - is the expensive half of BOTH + // 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 + }; + auto build_orbits = [&](const std::vector& rotations) -> Orbits { + Orbits orb; + orb.rep.resize(n); for (size_t i = 0; i < n; ++i) { if (!pass_cc[i] || !(Sigma[i] > 0.0)) continue; @@ -395,17 +396,29 @@ SearchSpaceGroupResult SearchSpaceGroup( if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l)) best = k2; } - rep[i] = best; - auto& g = grp[best]; + orb.rep[i] = best; + auto& g = orb.grp[best]; const double w = 1.0 / (Sigma[i] * Sigma[i]); g.sw += w; g.swI += w * I[i]; g.n += 1; } + return orb; + }; + + // Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry + // equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces + // non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is + // more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well + // - 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)) continue; - const auto& g = grp[rep[i]]; + const auto& g = grp.at(rep[i]); if (g.n < 2) continue; const double mean = g.swI / g.sw, dev = I[i] - mean; @@ -426,30 +439,14 @@ SearchSpaceGroupResult SearchSpaceGroup( // the systematic part of the scatter, which the fixed-sigma chi^2 ratio cannot: a genuine but // imperfectly-scaled high-symmetry merge and a twin can share a chi^2 ratio (~2) yet differ // sharply here (a genuine cubic step b x1.04 vs a merohedral twin b x1.6). - auto merge_systematic_b = [&](const std::vector& rotations) -> double { - struct Acc { double sw = 0.0, swI = 0.0; int n = 0; }; - std::unordered_map grp; - std::vector rep(n); - for (size_t i = 0; i < n; ++i) { - if (!pass_cc[i] || !(Sigma[i] > 0.0)) - continue; - HKLKey best = key[i]; - for (const auto& op : rotations) { - const auto m = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}}); - const HKLKey k2 = Canonicalize(m[0], m[1], m[2], opt.merge_friedel); - if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l)) - best = k2; - } - rep[i] = best; - auto& g = grp[best]; - const double w = 1.0 / (Sigma[i] * Sigma[i]); - g.sw += w; g.swI += w * I[i]; g.n += 1; - } + 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)) continue; - const auto& g = grp[rep[i]]; + const auto& g = grp.at(rep[i]); if (g.n < 2) continue; obs.push_back({I[i], Sigma[i], I[i] - g.swI / g.sw}); @@ -483,9 +480,13 @@ SearchSpaceGroupResult SearchSpaceGroup( const auto [present, min_class_cc] = point_group_present(pg.rotations); if (!present) continue; - const double ch = pg.rotations.empty() ? std::numeric_limits::quiet_NaN() - : chi2_under(pg.rotations); - const double be = pg.rotations.empty() ? 0.0 : merge_systematic_b(pg.rotations); + 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); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index fe837de5..8715ac06 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -523,8 +523,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // parallel (the main per-image loop then uses the forced lattice, no FFT). Size the pool to // those two so a rotation run doesn't pay to initialise (cuFFT plan + GPU alloc) indexers it // never uses. Stills fire the FFT once per image across all workers, so keep the default there. - if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) + if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) { indexing_settings.IndexingThreads(2); + // Two indexer threads leave the rest of the machine idle for the whole first pass, which is + // a third of the run. Split each scheme's candidate-cell refinement over -N/2 threads so the + // two together use -N. The split is by candidate cell, so the result does not depend on it. + indexing_settings.RefineThreads(std::max(1, config_.nthreads / 2)); + } indexer_pool = std::make_unique(indexing_settings, IndexerConstruction::OnFirstUse); indexer = std::make_unique(experiment_, indexer_pool.get()); if (!config_.reference_data.empty()) @@ -764,6 +769,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto coarse_settings = experiment_.GetIndexingSettings(); coarse_settings.FFT_HighResolution_A(3.5f); // low-order reflections only -> robust long axis coarse_settings.IndexingThreads(2); + coarse_settings.RefineThreads(std::max(1, config_.nthreads / 2)); IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- // resolution per-frame validation would be low (the coarse cell is metrically right but