# Viewer reductions, float mean accuracy and gzip stack loads - 2026-10-11

## Question

Work done two weeks earlier on a laptop (a full-scan mean that delegates to the
native reduction, 6.2 s to 0.13 s on a 512 x 512 scan) had to be carried onto
the restructured code. Three things were open: is every reduction of a viewer
on the acquisition's own session, does a float32 acquisition give the float64
mean pattern, and which stage of opening an acquisition is slower than it
should be?

## Setup

Nine acquisitions: an Arina 512 x 512 scan of a 192 x 192 detector (uint16), an
EMPAD RAW scan and a background-subtracted EMPAD HDF5 (128 x 128 scans,
float32), three TCMEP stacks from Zenodo 15084123 (a simulation, twisted STO
and PCCO, one gzip chunk per pattern), the SrTiO3 dislocation scan and the
uint8 simulation from Zenodo 7464234, and the WSe2 simulation from Zenodo
15987625. Each was loaded, reduced, opened in `Show4DSTEM`, saved as `.qem` and
reopened on an Apple M5 Max (MPS) and on an RTX PRO 6000 Blackwell (CUDA; the
GPU and its disk were shared with other jobs, so its times are not isolated).
Stage times are from `tests/benchmarks/viewer_stages.py`, median of 3 on fresh
loads, before (378e777) and after on the same machine.

## Float mean patterns were a float32 running total

`FloatANSResident.mean_dp_device` added each decoded chunk to a float32 total.
Against NumPy float64 over the same loaded values, in units in the last place
of float32:

| Acquisition | float32 running total | 32-position reads summed in float32 (the laptop loop) | Now |
|---|---|---|---|
| EMPAD RAW 128 x 128 | 5 (MPS), 4 (CUDA) | 16 | 0 |
| EMPAD HDF5 128 x 128 | 4 | 18 | 0 |
| Twisted STO 100 x 100 | 7 | 17 | 0 |
| PCCO 100 x 100 | 4 | 20 | 0 |
| SrTiO3 dislocations 256 x 256 | 8 | 38 | 0 |
| WSe2 128 x 128 | 4 | 17 | 0 |

CUDA now sums each chunk in float64. Metal has no float64: a kernel
(`resident/mps/kernels/pattern_total.msl`) keeps each detector pixel's float32
sum with two orders of rounding residual (Kahan-Babuska-Klein) and the host
adds the three in float64. Scan-ROI sums and means take the same total. The
mean pattern equals the float64 reference bit for bit on all six, and CUDA and
MPS return the same bits. The virtual images were already exact: their kernel
compensates on both backends.

The first version waited for the GPU three times per chunk (decode, Torch,
kernel) and cost 6% on the mean and up to 25% on a scan ROI. Decoding the chunk
and summing it in one command, with the dark pattern subtracted and the
selected scans flagged in the kernel, removed that: SrTiO3 mean 153 to 147 ms,
scan ROI 37.5 to 22.4 ms.

## Bounded views read every scan position for some acquisitions

`BoundedDetectorCompute` reduced on the acquisition's session only when it was
an encoded `Dataset4dstemGPU`. A GPU tensor held as a dense acquisition (the
widget's tutorial preview), a bare tensor and a paired acquisition were reduced
by 32-position reads. Every one has a backend of its own, so the view now takes
its owner's session for all of them and the read loops are gone. CUDA, GPU
shared:

| Acquisition viewed | Stage | Reads | Session | Values |
|---|---|---|---|---|
| Dense uint16, 192 x 192 x 96 x 96 | mean / BF / scan-ROI mean / viewer | 4.03 / 7.47 / 5.91 / 12.4 s | 0.071 / 0.044 / 0.019 / 0.26 s | identical |
| Dense float32, same shape | mean / BF / viewer | 4.35 / 5.51 / 10.7 s | 0.016 / 0.029 / 0.040 s | mean now 0 instead of 1 unit in the last place off float64; BF identical |
| Paired uint16, 512 x 512 x 192 x 192 | mean / BF / scan-ROI mean / viewer | 84.6 / 69.5 / 8.28 / 65.4 s | 0.014 / 0.036 / 0.166 / 0.14 s | identical |

On MPS the dense uint16 viewer went from 0.46 s to 0.11 s and the float32 one
from 0.39 s to 0.02 s. `TorchBackend` float sums moved to float64 accuracy with
the loops (float64 on CUDA, the Metal kernel on MPS), so no precision was lost.

### Signed data and scan regions of dense owners

Routing dense owners to their backends exposed what those backends did with
data that is not counts. A review of the branch found, and tests now cover:

- int16 data with negative totals (dark-subtracted frames): `ArrayBackend`,
  `TorchBackend.reduce_frames` and the CuPy `reduce_frames` and mean summed in
  uint64 or cast the total to it, so a negative total came back near 1.8e19,
  and the Torch maximum started from 0. Signed integers now accumulate in
  int64 and a maximum starts from the first chunk. Unsigned results are
  unchanged.
- A scan region's mean gathers its patterns: `ArrayBackend` and the CuPy
  backend took the whole selection at once (10.5 GB for a 400 x 400 region of
  128 x 128 float32 patterns). Both now gather in chunks.
- `ArrayBackend.reduce_frames(..., "mean")` of float32 data used a float32
  mean; it now divides the float64 total once.

## Metal scan-ROI reductions summed on the host

MPS `CountDetectorCompute.reduce_frames` decoded each run of selected scans,
copied it to the host and summed in NumPy: 203 ms for the 8245 patterns of a
51-scan-radius circle (98 ms NumPy sum, 79 ms in 103 decode commands, 28 ms
copy). A kernel (`streamed_counts_selected_reduce`) now adds the selected scans
per detector pixel on the device, skipping stream blocks without one: 30 ms for
that circle, 89 ms for 70681 positions, the same sums, maxima and means.

## gzip stacks inflated on one core

A stack of one gzip chunk per pattern is read by HDF5 one chunk at a time:
3.9 s of a 5.8 s load. The stored chunks are now read as they are and inflated
by zlib on up to 8 threads. Every loaded value is unchanged (the loader's
logical SHA-256 is the same).

| Stack | MPS before | MPS after | CUDA host before | CUDA host after |
|---|---|---|---|---|
| Twisted STO, 10000 x 256 x 256 | 3.69 s | 1.74 s | 4.83 s | 2.60 s |
| PCCO, 10000 x 128 x 128 | 2.21 s | 0.76 s | 2.54 s | 1.03 s |
| TCMEP simulation, 676 x 128 x 128 | 0.084 s | 0.049 s | 0.100 s | 0.063 s |

## Stage times on MPS, before and after

Median ms of 3; the float mean and scan-ROI values changed as described, every
other output is identical.

| Acquisition | load | mean | viewer | detector change | scan-ROI mean |
|---|---|---|---|---|---|
| Arina 512 x 512 x 192 x 192 uint16 | 1588 / 1626 | 130 / 131 | 141 / 141 | 6.2 / 6.2 | 449 / 46 |
| EMPAD RAW | 562 / 562 | 42 / 41 | 47 / 49 | 3.0 / 3.0 | 10.5 / 7.6 |
| EMPAD HDF5 | 534 / 544 | 42 / 43 | 51 / 51 | 10.9 / 10.8 | 10.4 / 7.9 |
| TCMEP simulation (676 patterns) | 86 / 49 | 4.2 / 4.5 | 11.7 / 12.6 | 2.2 / 2.2 | 1.1 / 1.7 |
| Twisted STO | 3672 / 1743 | 76 / 64 | 90 / 76 | 20.9 / 20.8 | 20.0 / 10.7 |
| PCCO | 2183 / 759 | 35 / 28 | 44 / 41 | 8.1 / 8.0 | 7.0 / 5.3 |
| SrTiO3 dislocations | 2188 / 2218 | 153 / 147 | 167 / 166 | 30.0 / 30.0 | 37.5 / 22.4 |
| SrTiO3 simulation uint8 | 247 / 246 | 14.3 / 14.3 | 19.2 / 19.0 | 2.3 / 2.2 | 29.3 / 5.3 |
| WSe2 simulation | 556 / 554 | 43 / 41 | 49 / 49 | 4.8 / 4.8 | 10.5 / 7.8 |

A cursor move takes 0.9 to 1.8 ms on every acquisition, before and after. On the
676-pattern simulation the scan-ROI mean is 0.6 ms slower: the fused path's fixed
cost (one window, one flag upload) exceeds what it saves on 13 patterns.

## Rejected

- A pairwise TwoSum tree in Torch operations for the Metal float total: about
  54 operations per 32 MiB window, several times slower than the decode.
- First-order compensation only (one residual): its error grows as
  n times eps squared, about 2e-10 of the total after 65536 patterns, which
  rounds a few pixels in every hundred to the neighbouring float32.
- Queuing the Metal sums of every window and waiting once: a decoded window is
  recycled by the next decode before a queued command has read it.
- Copying decoded patterns to the host and summing in float64 there.
- A tolerance on the float mean test: the comparison is exact.

## Left

- A `.qem` save writes its body twice (a temporary file, then the published
  file) and hashes it on one thread: for 4.3 GB on MPS 1.28 s final write,
  1.24 s SHA-256, 1.02 s temporary write and 0.36 s read back of a 3.9 s save.
- The exact uint64 products (`reduce_frames_exact`, `masked_sum_exact`) are defined
  for counts; signed data with a negative total has no uint64 value there.
- `TorchBackend.masked_sum` of float data sums in float32 (1 unit in the last
  place on the dense float32 stack above); the residents' virtual images are
  exact.
- The compensated Metal kernel costs more than a float32 sum on a dense float32
  stack: 12.8 ms against 2.8 ms for a float32 running sum of 1.36 GB.