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()) {