Add a roadmap for running MART on a GPU

[roytsmart]

Adds a Roadmap section to the docs with one page, docs/roadmap/cuda-mart.rst: the plan for running MartInverter on a GPU, and from there the ESIS Level-4 product.

The page covers:

• Where things stand: the pieces that already exist across regridding, ctis, optika and esis, and their status.
• What we learned: benchmarks of the MART loop on three problems. A single CSR operator makes a GPU iteration 29 to 223 times faster than ctis on the host, and 3.9 to 8.5 times faster than resampling one orthogonal element at a time on the GPU.
• Design:
  • the forward operator, and the conservative transpose of a per-channel operator for the backprojection, with broadcasting declared by shapes;
  • streaming two-pass assembly for ESIS-scale weights;
  • at ESIS scale, the sensor's charge diffusion applied as a convolution beside the matrix, with spectral lines kept on their own axis.
• Milestones: eight pull requests with their dependencies, and the gate each must pass. M7 and M8 reproduce tempest with diffusion off, within a tolerance agreed before the run.
• Memory at ESIS scale: the production grid needs about 48.5 GB in float32, so it needs an 80 GB GPU.
• Regenerate or cache: today a disk cache of the weights beats rebuilding them (about 40 s a run against 3 to 6 minutes), because optika's host-side preparation dominates the rebuild.
• Risks and open questions.

Notes:

• The three figures are static PNG renders.
• The benchmarks were rerun on 2026-10-12: all three problems back to back, each time the best of 3 runs on the host and 5 on the GPU.
• The regridding.Regridder prototype the page refers to is Add Regridder, a CSR operator for resampling on the host or a GPU regridding#70, a draft. The #69 it builds on is released in regridding 3.5.1.
• The diffusion plan follows the comment below and the optika sensor API in Spread the noiseless signal of a sensor with the kernel of its charge diffusion optika#255, which is planned for optika 4.0.
• The benchmark scripts are local and not part of this PR.

🤖 Generated with Claude Code

https://claude.ai/code/session_01GqceVSp8s2msxj1Dr9MrGU

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (b69e037) to head (dbafad7).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff            @@
##              main       #37   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           24        24           
  Lines          918       918           
=========================================
  Hits           918       918           

Flag	Coverage Δ
unittests	100.00% <ø> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[roytsmart]

Charge diffusion in the optika sensor, checked against this roadmap

Context. optika 3.0 (sun-data/optika#244 and a follow-up PR) will make the noise-free sensor model include charge diffusion. Today ImagingSensor.expose(noise=False) only applies a gain per pixel, while noise=True diffuses the charge, so the expected image is not the mean of the noisy ones.

After the change:

• The expected image is the photon image convolved, per wavelength, with a kernel. The kernel is the electrons per absorbed photon collected in each pixel around the one it was absorbed in, weighted by charge collection, and its sum is the old gain.
• The sensor stays separate from the optics. LinearSystem composes its weights with the sensor's operator, and the sensor gains transposes for the inverters.

The plan fits this roadmap, with the adjustments below.

What already lines up

• MART is unchanged, on the host and on the GPU. OptikaInstrument is unchanged too. The change is inside optika's image_from_weights and backproject_from_weights.

• Same backprojection convention. The sensor's backprojection is the kernel's transpose divided by the gain of each channel and wavelength. That's a scaled transpose like the conservative one here, so the division folds into the stored conservative values at build time.

• The kernel is the same object as IdealInstrument.psf:

  • it holds the fraction of the light landing in a pixel that reaches each pixel around it
  • it has odd size and is centered
  • it broadcasts over channel and wavelength

  optika's kernel can be passed to ctis as data. IdealInstrument(psf=...) against OptikaInstrument then makes a cross-check.

• The noise model matches M7. The per-wavelength variance ratio of tempest, constant for each channel and velocity slot, is optika's variance-to-mean ratio.

• The weights cache stays valid. With the sensor outside the weights, changing a diffusion model doesn't invalidate the weights cache. The kernels are cheap host computations, so optika needs no device path for them.

Adjustments

1. Don't fold diffusion into the CSR matrix at ESIS scale.
convolve_weights() is fine for the PSF tutorial. At production scale the backprojection rows average 3.5 entries. Folding in a k×k kernel grows a compact footprint to roughly (√3.5 + k − 1)² entries: about 4× for 3×3 and 10× for 5×5. The 15 µm pixels of ESIS need 5×5 to hold all but 1e-6 of the charge, and the exact slab model planned for optika, whose tails are heavier, may need more. That takes the forward matrix alone from 17.7 GB to roughly 70 to 180 GB in float32.

Applying the kernel as a small convolution beside the matrix costs two launches per iteration and no stored weights:

• forward: the matrix product, then the blur
• backprojection: the blur's transpose, then the transposed matrix

The same holds for an optical PSF of ESIS, if one is ever added.

2. Keep spectral lines on their own axis.
The forward matrix sums over wavelength inside the product, but the kernel depends on wavelength. Measured with optika#244 (ee4d7de) for e2v_ccd97 with 15 µm pixels and a 5×5 kernel:

• Across 580–640 Å: the kernel varies by up to 8e-4, which is 0.7% of a nearest-neighbor value.
• Within one line's velocity window: roughly 1e-5.

So the forward matrix should sum within each line and keep the lines on their own axis. One kernel per channel and line is then applied before the lines are summed. With the Regridder's shape rules that only needs the line to be its own axis of the scene grid. It costs about 0.5 GB of row pointers and 0.25 GB of intermediate images at production.

The per-wavelength gain folds into the forward values, and its inverse into the conservative values. The GPU path then matches optika's per-wavelength host model to about 1e-5.

3. The M7 and M8 gates.
Tempest was built without diffusion in the noise-free model. With optika 3.0, MART will also remove the diffusion blur, so Level-4 reconstructions will change. The gates need:

• a run with diffusion off (material.replace(diffusion=None)) to reproduce tempest
• a separate science check with diffusion on, like the one already planned for GPU-built weights

4. Exact and scaled transposes of the sensor.
ElasticNetInverter (#33) and the parametric inverter use the exact transpose. Their gradients need the exact transpose of the blur, while MART uses the scaled one. The optika sensor will provide both, mirroring the Regridder.

5. Determinism.
The blur should be computed as a gather, each output pixel summing its neighbors in a fixed order: depthwise conv2d in torch, never index_add_ or atomic adds.

6. The noise model as published data, instead of a probe.
The optika sensor will publish its kernel and a variance kernel. The variance of an image is the variance kernel convolved with the photon image:

• First, the variance kernel is just the variance-to-mean ratio times the kernel, today's model.
• Later, it becomes the exact variance kernel, with the same interface.

On the GPU that is a second matrix product, with the per-wavelength amplitude in the values, followed by one convolution, which is the variance-ratio pattern of this page. M7 can read these directly instead of measuring them with a noise probe.

Suggested edits to this page

1. Design: blurs on the sensor grid (diffusion, and an ESIS PSF if added) are applied as small convolutions beside the CSR matrix at ESIS scale, not folded into it, with the memory estimate above.
2. Design: ESIS scenes keep spectral lines on their own axis, so the forward matrix sums within lines.
3. M7: build the operator from the sensor's published kernel and variance kernel instead of a noise probe.
4. M7 and M8 gates: reproduce tempest with diffusion off, plus a science check with diffusion on.
5. Risks: from optika 3.0 the noise-free forward model includes diffusion, so Level-4 reconstructions change.

The planned optika sensor API, for reference

Not final until the follow-up PR:

AbstractImagingSensor.kernel(wavelength, direction, num=None)           # electrons per absorbed photon, per offset
AbstractImagingSensor.kernel_variance(wavelength, direction, num=None)  # variance contributed per photon
AbstractImagingSensor.expose(..., noise=False)                          # the kernel applied to the photons
AbstractImagingSensor.expose_transposed(image, ...)                     # exact transpose, for gradient inverters
AbstractImagingSensor.backproject(image, ...)                           # scaled transpose, for MART

num=None sizes the kernel to hold all but 1e-6 of the charge of the widest cloud.

🤖 Generated with Claude Code