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

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 <algorithm>
 #include <cmath>
@@ -150,7 +151,8 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &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<MergedReflection> &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<size_t>(nthreads, 1);
+        const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
+        std::vector<std::string> 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<std::streamsize>(s.size()));
     }
 
     out << "#\n";
@@ -435,7 +462,8 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
 
 void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
                               const DiffractionExperiment &experiment,
-                              const std::string &filename) {
+                              const std::string &filename,
+                              size_t nthreads) {
     bool has_anom = true;
     const std::vector<MergedOutRow> rows = BuildMergedRows(reflections, experiment, has_anom);
 
@@ -459,20 +487,53 @@ void WriteShelxHklReflections(const std::vector<MergedReflection> &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<size_t>(std::clamp(n, 0, static_cast<int>(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<size_t>(nthreads, 1);
+        const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
+        std::vector<std::string> 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<std::streamsize>(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<std::streamsize>(tail.size()));
     out.close();
 }
 
@@ -482,11 +543,13 @@ void WriteReflections(const std::vector<MergedReflection> &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);
 }