diff --git a/docs/index.rst b/docs/index.rst index f3f5a52..b0347dd 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -234,6 +234,16 @@ Slide decks for presenting this package. * `Simulating and Inverting a CTIS Observation `__: a step-by-step walkthrough of :doc:`tutorials/simple-mart`. +Roadmap +======= + +Plans for the development of this package. + +.. toctree:: + :maxdepth: 1 + + roadmap/cuda-mart + Citation ======== diff --git a/docs/roadmap/cuda-mart.rst b/docs/roadmap/cuda-mart.rst new file mode 100644 index 0000000..fcf736c --- /dev/null +++ b/docs/roadmap/cuda-mart.rst @@ -0,0 +1,609 @@ +CUDA MART Roadmap +================= + +*As of 2026-10-12.* + +Every numerical piece of a CUDA implementation of MART now exists, +and a prototype ``regridding.Regridder`` joins most of them: +one sparse CSR operator with both of its transposes, assembled on the GPU. +What remains is a GPU path in :class:`~ctis.inverters.MartInverter` +and, for the ESIS Level-4 product, +a streaming build for weights too large to hold at once, +with the charge diffusion of the sensor applied beside the matrix. + +Where things stand +------------------ + +.. list-table:: + :header-rows: 1 + :widths: 22 22 18 38 + + * - Piece + - Where + - Status + - What it gives a CUDA MART + * - Device weights, ``regridding.weights(device="cuda")`` + - regridding + - Released + - Conservative weights built on the GPU; + the output grid must be a uniform, axis-aligned lattice + * - ``regridding.convolve_weights()`` on the device + - regridding 3.5.0 + - Released + - A PSF folded into the weights without leaving the GPU + * - ``regridding.regrid_from_weights()`` on the device + - regridding + - Released + - Applies device weights, one kernel launch per orthogonal element; + rejects ``torch`` tensors + * - ``regridding.transpose_weights_conservative()`` on the device + - `regridding #69 `__ + - Released in 3.5.1 on 2026-10-08 + - Conservative transposes of device weights; + its per-grid cell volumes give the CSR operator its conservative scaling + * - ``regridding.Regridder`` + - `regridding #70 `__, on main + - Draft prototype, 79 tests, CI green + - Every channel and wavelength in one CSR matrix, + with broadcasting declared by shapes, exact and conservative transposes, + and assembly on the host or GPU that agree bit for bit + * - :class:`~ctis.inverters.MartInverter` + - ctis 0.5.0 + - Released + - The reference algorithm, on the host through :mod:`named_arrays`, + with PSF support in :class:`~ctis.instruments.IdealInstrument` + * - ``torch`` ``Regridder`` + - ctis ``feature/parametric-inverter``, ``feature/elastic-net-inverter`` + - Unmerged + - Weights as one CSR matrix with an explicit transpose; + its CUDA product is that of ``torch``, which does not repeat exactly on long rows; + the regridding prototype generalizes it + * - Device weights in optika, ``LinearSystem.weights(device="cuda")`` + - `optika #229 `__ + - Released in 3.1.0 + - The weights of an optical model built on the GPU; + its distortion, vignetting and effective area are still prepared on the host + * - Sensor kernels, with charge diffusion + - `optika #255 `__, for optika 4.0 + - Draft + - The blur of the noise-free sensor as a kernel per wavelength, + its exact and scaled transposes, and a variance kernel for the noise model + * - Level-4 GPU MART scripts + - esis ``wip/level4-tempest-tooling`` + - Unmerged scripts + - The production run: 30 frames in about 10 minutes on one H100, + from CPU-built weights + * - MART benchmarks + - Local scripts + - Throwaway + - The MART loop of :class:`~ctis.inverters.MartInverter` on raw ``torch`` + tensors, with per-element device weights and with the CSR operator; + both match the solution of ctis to 2e-14 of its peak + +What we learned +--------------- + +One CSR operator makes a MART iteration on the GPU 29 to 223 times faster than ctis, +and 3.9 to 8.5 times faster than the per-element GPU path, +by replacing 160 kernel launches with two. +On the host alone it is 2.1 to 8 times faster than ctis. + +.. figure:: cuda-mart/iteration-times.png + :width: 100% + :alt: Milliseconds per MART iteration for three problems and four implementations, on a log scale. + + Milliseconds per MART iteration on an RTX 4090 and 48 host threads, + for the tutorial, the PSF tutorial with 21 by 21 PSFs, + and a 256 by 256 scene over 50 fixed iterations. + All three ran back to back on 2026-10-12, + with other applications keeping the GPU about 10% busy; + each time is the best of 3 runs on the host and 5 on the GPU. + Every implementation matched the solution of ctis to 2e-14 of its peak. + +- **Launches dominate the GPU.** + ``regrid_from_weights()`` launches one kernel per channel and wavelength, + 80 per regrid here, costing about 11 ms an iteration. + The tempest scripts (one flat table per channel and ``index_add_``) + and the ``torch`` ``Regridder`` (one CSR matrix) both avoid this. +- **Setup now dominates short runs.** + On the PSF tutorial the CSR GPU path spends 2.7 s building its operators + and 0.3 s iterating, + so end to end it is 4.0 times faster than ctis, + against 5.4 for the per-element path. + Building the backprojection straight from the weights would remove + one of its three assemblies and the transpose. +- **Short rows needed fewer threads.** + The kernel gave every row a 32-thread warp, + but without a PSF the rows of a backprojection hold about 3 entries each. + It now gives each row as many threads as an average row has entries, + rounded down to a power of two, up to a warp: + 2 for a backprojection without a PSF, that of ESIS included, + and 32 for a forward operator or a backprojection with the 21 by 21 PSFs. + At 256 by 256, the backprojection product went from 3.5 ms to under 0.5, + against 0.6 to 0.7 for the sparse product of ``torch`` (cuSPARSE). + Timed back to back, a MART iteration went from 6.1 ms to 2.8, + against 2.7 for ``torch``; + the gap left is about 0.2 ms of host time per product, + which only small problems notice. +- **Assembly is cheap once parallel.** + At 145M weights, assembling takes about 0.5 s on the GPU and on 48 host threads alike, + and the transpose another 0.5 s. + The host and GPU build the same matrix bit for bit. +- **Most of the host time of ctis is not regridding, except with large PSFs.** + The same algorithm on raw arrays, with the same per-element regridding, + takes 18 ms instead of 50 on the tutorial, + and 247 instead of 628 on the 256 by 256 scene. + The rest is :mod:`named_arrays` and :mod:`astropy` overhead: + units, uncertain arrays, the residual correlation. + On the PSF tutorial, where regridding 148M weights dominates, + it takes 184 ms instead of 199. +- **Building the weights is the expensive step at production scale.** + Tempest needed about 1.3 TB of cached CPU weights and 1.2 TB jobs per frame. + Built on the GPU, the 148M weights of the PSF tutorial and their transpose + take 0.95 s instead of 2.7 s. +- **MART needs the conservative transpose, which is a scaled exact one.** + The adjoint of the ``torch`` ``Regridder`` is the exact transpose, + but the backprojection of MART must conserve flux. + That is the exact transpose with each entry scaled by an input-cell factor + over an output-cell volume, so both share one set of indices. +- **Atomic adds are not deterministic.** + The device scatter adds in no fixed order, + and the ratios of MART amplified that to 1e-10 where both sides are near zero. + CSR gives identical results every run. +- **The sparse product of** ``torch`` **does not repeat exactly either.** + On rows of 4,500 entries, all 100 repeated products differed, + by 1 to 4 ulp in about half the values, + even under ``torch.use_deterministic_algorithms(True)``. + Short rows repeat, but the result for a column changes with the batch width. + The kernel of the regridding ``Regridder`` repeats exactly in every case. + The docstring of the ``torch`` ``Regridder`` in ctis calls its product + bitwise deterministic on CUDA, which is wrong for long rows. +- **Precision and memory trade off.** + Tempest stored values as float32 and indices as int32 + to fit the production grid on one 80 GB GPU. + The benchmarks used float64 throughout to match ctis exactly. + Measured, the CSR operator stores 12 bytes per weight per direction in float64 + and 8 in float32. +- **The GPU must be quiet to benchmark.** + On the PSF tutorial the per-element path took 21 ms an iteration on a quiet GPU, + 27 ms with other applications keeping it about 10% busy, + and 57 ms with video apps open. + Single runs on a shared GPU varied by up to 30%, + so the figure reports the best of several. +- **Stopping rules differ.** + The PSF tutorial stops at the first :math:`\langle \chi^2 \rangle \leq 1`; + tempest stops when :math:`\langle \chi^2 \rangle` improves by less than 1e-2; + ctis stops when it improves by less than 1e-3. + +Design +------ + +The CUDA MART is :class:`~ctis.inverters.MartInverter` running on ``regridding.Regridder``: +two sparse CSR matrices, the forward model and a conservative transpose for the backprojection, +assembled on the GPU from the device weights of regridding. +Only the merit, one number per channel, returns to the host each iteration, +until the solution does. + +.. figure:: cuda-mart/design.png + :width: 100% + :alt: The CUDA MART in ctis: weights built once on the GPU, then a five-step loop. + + Build once, then a five-step loop on the GPU. + +The loop is the host algorithm unchanged; +only the two products and where the arrays live differ. + +- **Operators as CSR, in regridding.** + ``regridding.Regridder`` holds every channel and wavelength in one matrix, + assembled from the device weights of regridding on the GPU, + with the ``-1`` slots dropped and no round trip through the host. + Each row sums in a fixed order, + so the host and GPU build the same matrix bit for bit, + and an iteration is two launches instead of 160. + Each row is summed by as many GPU threads as an average row has entries, + rounded down to a power of two, up to a warp, + so the short rows of a backprojection no longer waste most of a warp. +- **Broadcasting declared by shapes.** + The operator takes the shapes of the values going in and coming out. + Where the weights vary along an axis, + a length of 1 in the input means broadcast inside the matrix, + and a length of 1 in the output means summed. + Axes the weights lack are batch columns. + The forward model is one product from (wavelength, x, y) to (channel, x′, y′), + with the scene shared by the channels. + Its transpose would sum the channels, + but MART backprojects each channel on its own, + so the backprojection is the transpose of a second operator on the same weights, + which takes a scene per channel, (channel, wavelength, x, y). +- **Conservative values on the indices of the transpose.** + The conservative transpose scales each entry of the exact one + by an input-cell factor over an output-cell volume, + taken from the element of the weights the entry came from. + For the backprojection of MART, whose columns keep the channels apart, + scaling the rows and columns would do the same; + the prototype scales each entry, + so that the transpose of an operator sharing its input across channels is right too. + Each entry is scaled once at build time. + MART needs no exact values, + so its backprojection stores only the conservative ones: + 4 bytes per entry in float32. + A gradient inverter needs the exact values instead. +- **Inside** :class:`~ctis.inverters.MartInverter`\ **, behind a** ``device`` **option.** + The update rule and stopping tests stay shared with the host path. + :mod:`named_arrays` stays at the boundaries: + images in, a :class:`~named_arrays.FunctionArray` solution out. +- **float64 by default, float32 as an option.** + float64 matches the host to rounding. + float32 cuts each stored weight from 12 to 8 bytes per direction, + and the arrays of MART by half, + which the ESIS production grid needs. + Each row still sums in float64. + The products of the prototype read and return float64 whatever the matrix holds, + so M3 must let them read and return float32, + or the arrays of MART at production grow from 8.1 GB to 16.2. + Column indices are int32 whenever the matrix has fewer than :math:`2^{31}` columns. +- **The noise model of the instrument.** + :class:`~ctis.instruments.IdealInstrument` gets shot plus read noise, as on the host. + :class:`~ctis.instruments.OptikaInstrument` gets the per-wavelength model of its sensor + in a later milestone, + read from the kernels the sensor publishes from optika 4.0 + instead of measured with a noise probe. + The variance of an image is the variance kernel convolved with the photon image: + on the GPU, a second product with the amplitude of each wavelength in its values, + then one convolution. + The variance kernel starts as the variance-to-mean ratio times the kernel, + the model of tempest, + and can become exact later behind the same interface. +- **Streaming two-pass assembly at ESIS scale.** + The production weights are 35 to 53 GB on the device, + more than any GPU holds beside the matrix. + A builder takes them a block at a time, twice: + the first pass counts the entries of each row, which fixes the size of the matrix, + and the second places them. + Exact zeros, which the field stop leaves behind, are dropped. + Each block carries an explicit row and column offset and a scale, + so the tied lines of ESIS fold into their windows. +- **Blurs on the pixel grid beside the matrix, at ESIS scale.** + From optika 4.0, the noise-free model of a sensor includes charge diffusion: + the expected image is the photon image convolved, per wavelength, with a small kernel. + Folding that kernel into the weights, as ``convolve_weights()`` does for the PSF tutorial, + would grow the 3.5 pixels each scene cell touches at production + to about :math:`(\sqrt{3.5} + k - 1)^2` for a k by k kernel. + With the slab model of the e2v CCD97 in optika #255, + the 15 µm pixels of ESIS need 9 by 9 to hold all but 1e-6 of the charge + (5 by 5 leaves out 5e-4), + and the default kernel of the sensor is 11 by 11: + 28 to 40 times the entries. + The forward matrix would grow from 17.7 GB to 490 to 710 GB in float32. + Instead, the forward product is followed by the blur, + and the backprojection starts with its transpose: + two more launches an iteration, and no stored weights. + The blur is a gather, each pixel summing its neighbors in a fixed order, + so it repeats exactly like the products. + The sensor supplies both of its transposes, as the ``Regridder`` does: + the exact one for gradient inverters and the scaled one for MART. + An optical PSF of ESIS, if one is added, goes the same way. +- **Spectral lines on their own axis, at ESIS scale.** + The diffusion kernel depends on wavelength, but the forward product sums over it. + For the e2v CCD97 of optika #255 with 15 µm pixels, + the gain falls by a quarter across 580 to 640 Å, + but the shape of the kernel changes little: + by up to 1.1e-4, 0.11% of a nearest-neighbor value, + which moves 2.3e-4 of the charge, + and within the ±210 km/s window of one line it moves under 3e-6. + So the scene of ESIS keeps each line on its own axis, + the forward matrix sums only within a line, + and one kernel per channel and line blurs each image before the lines are summed. + The shape rules of the ``Regridder`` need nothing new for this. + The gain of each wavelength folds into the forward values, + and its inverse into the conservative ones, + so the GPU matches the per-wavelength model of :mod:`optika` to a few parts in a million. + At production this costs about 0.5 GB of row pointers and 0.25 GB of intermediate images. +- **Not chosen: a single-launch scatter, or a second matrix of transposed weights.** + The scatter would remove the launch overhead but stay nondeterministic; + CSR removes both. + A second matrix would repeat the indices of the transpose, + when only its values differ. +- **Not chosen: the sparse product of** ``torch`` **(cuSPARSE).** + Its results do not repeat exactly on long rows, + even in the deterministic mode of ``torch``, + and change with the batch width. + It cannot sum float32 entries in float64, + and its 64-bit row pointers force 64-bit column indices, + 8.8 GB more per direction at ESIS scale. + Our kernel now takes less time than it for the two products of MART together, + from 17M weights up. + +Milestones +---------- + +Eight pull requests take the GPU MART from prototype to the ESIS Level-4 product. +M5 is the first one users can run. +The ``Regridder`` prototype already covers most of M3; +M1 shipped in regridding 3.5.1, M2 and M4 can start now, +M3 is a draft pull request no longer waiting on M1, and M6 needs M3. + +.. figure:: cuda-mart/milestones.png + :width: 100% + :alt: Eight milestones and the ones each needs. + + Milestones M1 to M8, one pull request each, with the gate under the name. + Arrows point from a milestone to the ones that need it. + +M2 is off the path of MART, since the loop uses CSR products +rather than ``regrid_from_weights()``; +it fixes the ``torch`` bug the benchmark found. +The GPU transpose kernel of #69 is off the path too, +because the operator scales its own transpose, +but the cell-volume code of #69 is what that scaling uses. +M3 should be reviewed by the author of the ``torch`` ``Regridder``, +which it generalizes. +M6 is new: the production weights do not fit on any GPU at once, +so the operator has to be built from blocks of weights produced on demand. +M7 needs optika 4.0, whose sensor publishes the kernels of its blur and noise. + +Validation and benchmarks +------------------------- + +Every milestone is proven against the host path it replaces, then timed on an idle GPU. +The host results are the reference throughout, +since :class:`~ctis.inverters.MartInverter` is already validated. + +.. list-table:: + :header-rows: 1 + :widths: 10 50 40 + + * - Milestone + - Correct when + - Fast when + * - M1 + - The device tests of #69 pass on the GPU runner of regridding + - Device transpose at or under 166 ms for 35M weights + * - M2 + - A ``torch`` tensor resamples to the same result as a :mod:`numpy` array, + with no copy + - No target + * - M3 + - Block-diagonal products equal ``regrid_from_weights()`` + and the conservative transpose equals ``transpose_weights_conservative()``, + bit for bit; GPU assembly equals host assembly. + All three pass in the prototype + - One launch per product; + the two products together no slower than cuSPARSE from 17M weights up; + assembly in about 0.5 s per 145M weights. + The prototype meets all three + * - M4 + - Device weights match the host weights to 1e-12 in float64, + compared as dense matrices once the ``-1`` slots are dropped, + as the device tests of regridding compare them + - The weights of the PSF tutorial and their transpose in under 1 s (0.95 s in the prototype) + * - M5 + - The same iteration count as the host :class:`~ctis.inverters.MartInverter`, + and :math:`\langle \chi^2 \rangle` and the solution within 1e-10 of the host's, + on all three benchmark problems + - At least 6 times faster end to end on the PSF tutorial + * - M6 + - On weights that fit, the streamed build equals ``Regridder.from_weights()`` + bit for bit once the exact zeros are dropped from both; + each block's count agrees between the two passes + - Production assembly, beyond building the weights, + in under a minute on an 80 GB GPU + * - M7 + - With diffusion off (``diffusion=None`` on the material of the sensor), + ESIS frame 15 matches ``reproduce_level4.py`` of tempest + within a tolerance agreed before the run; + with it on, a science check + - The production grid fits one GPU; + tempest ran 30 frames in about 10 minutes on an H100 + * - M8 + - With diffusion off, the 30-frame cube matches the distributed Level-4 FITS files + within the tolerance of M7; + with it on, a science check + - No slower than tempest + +The benchmark problems are the tutorial, the PSF tutorial and the 256 by 256 scene. +ctis needs a GPU CI job before M4, +like the ``tests-cuda`` workflow of regridding, +which runs on the runner the ``GPU_RUNNER`` variable names. + +Memory at ESIS scale +-------------------- + +The ESIS production grid, 0.75″ pitch and 24 velocity bins, +needs about 48.5 GB in float32: +it fits an 80 GB GPU, not a 24 GB one. +Only the production row is measured, by the count of tempest of 2.21 billion weights; +the others scale it by grid. + +.. list-table:: + :header-rows: 1 + + * - Grid + - Weights per direction + - Forward (GB) + - Backprojection (GB) + - MART arrays (GB) + - Total, float32 (GB) + - Total, float64 (GB) + * - 0.75″, 24 bins (production) + - 2.21 billion + - 17.7 + - 22.7 + - 8.1 + - 48.5 + - 74 + * - 1.0″, 24 bins + - 1.72 billion + - 13.8 + - 16.6 + - 4.6 + - 35.0 + - 53 + * - 0.75″, 14 bins + - 1.29 billion + - 10.4 + - 13.2 + - 4.7 + - 28.3 + - 43 + * - 1.5″, 24 bins + - 1.29 billion + - 10.4 + - 11.6 + - 2.0 + - 24.0 + - 36 + * - 1.0″, 14 bins + - 1.00 billion + - 8.1 + - 9.7 + - 2.7 + - 20.4 + - 31 + * - 1.5″, 14 bins + - 0.75 billion + - 6.1 + - 6.8 + - 1.2 + - 14.0 + - 21 + +The costs are measured: 8 bytes per weight per direction in float32, +12 in float64, plus 8 per row. +The backprojection stores only the conservative values MART uses; +keeping the exact ones beside them would add 8.8 GB at production in float32. +The arrays of MART are counted in float32, +which the products of the prototype do not yet allow. +The backprojection keeps each channel apart, +so at production it has 623M rows, whose pointers alone take 5 GB. +The scene is 1139 cells across at 0.75″, +from the raytrace field of the distortion fit padded by 1.25, as tempest builds it. +Fewer velocity bins shrink the problem most; +coarser cells help less, since each one touches more detector pixels. +Keeping the spectral lines on their own axis, for the blur of the sensor, +adds about 0.75 GB at production. + +Regenerate or cache +------------------- + +The streamed build needs every block of weights twice. +A disk cache supplies them fastest today: about 40 s a run, +against 3.3 to 6.2 minutes for rebuilding. +What makes rebuilding slow is the host-side preparation in :mod:`optika`, +22.2 s per spectral line, not the GPU build, which takes 3.2 s. +Moving that preparation to the GPU would bring a rebuild within about 20 s of the cache. + +.. list-table:: + :header-rows: 1 + :widths: 44 22 34 + + * - Way to supply every block twice + - Per run, production grid + - Needs + * - Disk cache of the compacted forward weights + - 38 s + - 26.5 GB of NVMe, built once in about 3 minutes + * - The existing cache of tempest + - At least 79 s, if it were local + - 1.3 TB on the cluster; reads about 94 GB per pass + * - Build once on the GPU, park the weights in host RAM + - 3.3 minutes + - 26.5 GB of pinned host RAM + * - Rebuild each pass, keeping the preparation between passes + - 3.6 minutes + - 21 GB of host RAM + * - Rebuild everything each pass + - 6.2 minutes + - Nothing extra + +The rates are measured for one spectral line (four channels, 24 velocity bins) +on the flight 1 model of ESIS, and scaled to its seven lines, +on a workstation with an RTX 4090, an NVMe SSD reading 3.0 GB/s, +and pinned copies between host and GPU at 13 GB/s. +The preparation in :mod:`optika` (distortion, vignetting and effective area on the host) +takes 22.2 s; +the GPU build takes 3.2 s for 557M slots, and assembly 0.58 s per channel. +Linearizing the optics takes another 10.4 s per line once, which esis already caches. +Every option pays 16 s of assembly. + +Risks and open questions +------------------------ + +GPU memory is the main risk; +the rest are questions to settle before the milestone that needs them. + +- **GPU memory (M6, M8).** + The ESIS production grid needs about 48.5 GB in float32, + so an 80 GB GPU like the H100 tempest used. + A 24 GB card fits 14 velocity bins at 1.0″ only tightly; + see `Memory at ESIS scale`_. +- **A** ``device`` **option or a separate inverter (M5).** + The design assumes the GPU path goes inside :class:`~ctis.inverters.MartInverter` + behind a ``device`` option, rather than in a separate inverter. +- **Regenerate or cache (M6).** + A disk cache is fastest today, but needs about 27 GB of fast disk + and has to be invalidated when the optics, grids or code change; + moving the preparation in :mod:`optika` to the GPU would make rebuilding nearly as fast. + See `Regenerate or cache`_. +- **Sorting for CSR (resolved).** + The prototype sorts one block of weights at a time + and places it at the cursors of its rows, + so assembly needs one block's sort beyond the matrix, not a second copy of it. +- **optika weights on the GPU (M7).** + The device clipping needs a uniform, axis-aligned output lattice, + which the pixel grid of a sensor is, + and from optika 3.1.0 ``LinearSystem.weights()`` takes a ``device`` + and passes it to ``named_arrays.regridding.weights()`` + (`optika #229 `__). + ctis allows optika 3.1, + so :class:`~ctis.instruments.OptikaInstrument` needs only to pass its device on. + The distortion, vignetting and effective area are still prepared on the host, + 22.2 s per spectral line at production; + see `Regenerate or cache`_. +- **named-arrays cannot hold device arrays.** + The GPU path converts at its boundaries, which costs little. + A ``torch`` backend for :mod:`named_arrays` would remove that, + but is a far larger project. +- **Which stopping rule (M5).** + The GPU path should keep the host's: + today ctis stops when :math:`\langle \chi^2 \rangle` improves by less than 1e-3. + Tempest used 1e-2, and the PSF tutorial stops at :math:`\langle \chi^2 \rangle \leq 1`. +- **Coordination (M3, M8).** + Both build on unmerged work by a teammate, the ``torch`` ``Regridder`` + and the tempest scripts. + Ownership needs agreeing before either starts. +- **Where the ties of ESIS live (M6, M8).** + The weights of a tied line fold into the columns of its window with a fixed ratio, + and the backprojection of tempest blends the transpose of each line + by its share of the flux of the window. + The leaning is an esis-side loop over a generic builder, + so regridding knows nothing of ESIS. +- **GPU-built weights differ slightly (M4, M5, M7, M8).** + The GPU clipping keeps a few tiny overlaps the CPU build drops, + under 1e-12 in the device tests of regridding. + So M4 and M5 compare to a tolerance, + and a GPU-built Level-4 run will not match tempest's bit for bit: + M7 and M8 need a tolerance agreed before the run, and a science check. +- **The second product of the noise model (M7).** + The variance of an image needs a second forward product, + with the variance-to-mean ratio of each wavelength in its values. + Stored as a second values array on the forward indices, as tempest stores its ratio, + that is 8.8 GB at production in float32, + beyond the 48.5 GB of `Memory at ESIS scale`_. + The ratio is constant for each channel and wavelength, + so the product could compute it from the row and column of each entry instead. +- **Charge diffusion changes Level-4 (M7, M8).** + From optika 4.0, the noise-free model of a sensor includes charge diffusion, + so MART will also undo the diffusion blur, + and Level-4 reconstructions will differ from those of tempest. + ctis allows only optika below 4, so nothing changes until ctis raises that cap. + M7 and M8 therefore reproduce tempest with diffusion off, + and check the science with it on. +- **Host time per GPU call (M5).** + Each ``Regridder`` call spends about 0.2 ms on the host, + in the launch of ``numba``, its array wrappers and the layout in ``torch``. + That is why ``torch`` still wins a whole iteration on small problems: + timed back to back, 0.9 ms against 1.1 on the tutorial, + and 2.7 against 2.8 at 256 by 256. + Launching without the wrappers of ``numba`` could cut it; + at ESIS scale it is negligible. diff --git a/docs/roadmap/cuda-mart/design.png b/docs/roadmap/cuda-mart/design.png new file mode 100644 index 0000000..1ee531e Binary files /dev/null and b/docs/roadmap/cuda-mart/design.png differ diff --git a/docs/roadmap/cuda-mart/iteration-times.png b/docs/roadmap/cuda-mart/iteration-times.png new file mode 100644 index 0000000..c52856b Binary files /dev/null and b/docs/roadmap/cuda-mart/iteration-times.png differ diff --git a/docs/roadmap/cuda-mart/milestones.png b/docs/roadmap/cuda-mart/milestones.png new file mode 100644 index 0000000..d29fb58 Binary files /dev/null and b/docs/roadmap/cuda-mart/milestones.png differ