unmerged MTZ: a batch header for every image the sweep spans

The header set was built from the batches that produced an observation, so an
image that indexed nothing left a hole in the phi series. AIMLESS starts a new
run at such a discontinuity: measured on a 360 degree sweep with 67 unindexed
images it made 23 runs of one sweep, its scale model diverged, Rmeas overflowed
to -1266 and the result no longer correlated with rugnux's own merge (Pearson
0.0018). Renumbering the batches contiguously does not help - the split is on
phi, not on numbering - so the fix is a header per image over the span the
observations cover.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N

docs/CHANGELOG.md

@@ -7,6 +7,7 @@
 * `rugnux` measures the beam centre on every run, and indexes with it when the file's value indexes nothing.
 * `rugnux` places a detector swung out on a 2theta arm where the file says it stands.
 * `rugnux` writes the unmerged MTZ by default, with a P1 merge beside it, so a wrong space group can be re-merged without reprocessing.
+* `rugnux` writes a batch header in the unmerged MTZ for every image the observations span, not only for the images that produced one, so a scaling program reads one run per sweep.
 * `rugnux` handles symmetry better: the lattice, the point group, the setting and the systematic absences.
 * `rugnux` fits the direction of the goniometer axis and not its length, so the cell chosen by the first pass is the one its own refinement scored.
 * `rugnux` reports the detector geometry it measured - the direct beam, the tilt and the beam centre - and what a single sweep can and cannot determine.

image_analysis/WriteReflections.cpp

@@ -12,7 +12,6 @@
 #include <cmath>
 #include <cstdio>
 #include <map>
-#include <set>
 #include <tuple>
 #include <fstream>
 #include <iomanip>
@@ -748,7 +747,13 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     // way to recover it. The partiality is NOT divided out - that is a scale, FRACTIONCALC carries
     // it, and every program this file is for wants to handle it its own way.
     gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
-    std::set<int> batch_numbers;
+    // The batch headers are written for every image the observations span, not only for the images
+    // that produced one. AIMLESS starts a new run wherever the phi series jumps, so an image that
+    // indexed nothing leaves a hole that cuts the sweep in two: measured on a 360 deg sweep with 67
+    // unindexed images, AIMLESS made 23 runs of it and its scale model diverged - CC(1/2) 0.811,
+    // Rmeas overflowed to -1266, and the result no longer correlated with rugnux's own merge
+    // (Pearson 0.0018). Filling the gaps gives one run and the numbers the sweep deserves.
+    int first_batch = 0, last_batch = -1;
 
     // Lattice-centring absences are integrated on purpose - prediction runs in P so the space-group
     // search can confirm or disprove the centring - but they are not written here. POINTLESS reads
@@ -766,7 +771,9 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
         // A summed full's image_number is its rocking-curve centroid, so this is the batch the
         // reflection is centred on - which is what a batch means for a full everywhere else.
         const int batch = 1 + static_cast<int>(std::lround(r.image_number));
-        batch_numbers.insert(batch);
+        if (last_batch < first_batch) { first_batch = batch; last_batch = batch; }
+        first_batch = std::min(first_batch, batch);
+        last_batch = std::max(last_batch, batch);
         mtz.data.push_back(static_cast<float>(hkl[0]));
         mtz.data.push_back(static_cast<float>(hkl[1]));
         mtz.data.push_back(static_cast<float>(hkl[2]));
@@ -861,7 +868,7 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
     batch.floats[115] = 1.0f;
     batch.floats[116] = static_cast<float>(experiment.GetYPixelsNum());
 
-    for (const int number : batch_numbers) {
+    for (int number = first_batch; number <= last_batch; number++) {
         batch.number = number;
         batch.floats[36] = phi_start_deg(static_cast<float>(number - 1));   // phistt
         batch.floats[37] = batch.floats[36] + wedge_deg;                    // phiend