Add utu.dem.plowman, the regularized DEM inversion of Plowman & Caspi (2020)

[roytsmart]

Summary

Adds utu.dem.plowman(), the regularized differential emission measure (DEM) inversion of Plowman & Caspi (2020). It is a Numba port of simple_reg_dem from EMToolKit, run in parallel over pixels, with a named-arrays front end.

dem, chi2 = utu.dem.plowman(
    intensity=intensity,        # DN/s, with a channel axis; every other axis is a pixel
    uncertainty=uncertainty,    # same units, broadcastable against intensity
    response=response,          # na.FunctionArray: temperature -> (temperature, channel) response
    axis_channel="channel",
    axis_temperature="temperature",
)

dem is an na.FunctionArray over the response's temperatures, in units of intensity over the units of the response (cm^-5 for AIA). chi2 is the reduced chi squared of each pixel, or -1 where the first step failed.

If the intensity, uncertainty or response is an na.UncertainScalarArray, the DEM and chi2 are too: the nominal values are inverted, and each sample of the distribution is inverted as one more pixel.

Arguments that can't be inverted are rejected before anything runs. That covers temperatures without units, not increasing, or not matching the responses; an intensity without the channel axis, or with the temperature axis; non-positive smoothness or floor; and an uncertainty or floor in units other than the intensity's.

The module is instrument-agnostic: responses and uncertainties come from the caller.

Accuracy and speed

Against EMToolKit itself:

	Agreement with EMToolKit	Speed, one core
Real AIA, 40k pixels, six channels	~1e-11 relative	~50k px/s against ~1,000
The random DEMs of the paper	~1e-10 relative	~33k px/s against ~1,000

A full 4096 x 4096 six-channel AIA image takes about 12 s on 24 cores, against about 3.6 h for the reference. EMToolKit's process-pool version is no faster than its serial one, because it pickles the whole cube for every pixel.

It cannot agree to the last bit. The regularization matrix is singular on its own, so where the data say little the linear systems are near-singular, and any two Cholesky implementations drift apart over the iteration. A port of the reference with numpy.linalg in place of scipy.linalg disagrees with it by as much. The worst cases are pixels with no signal, at ~1e-6.

The kernel allocates nothing inside the iteration, and compiles with only the contract and reassoc fast-math flags. They are worth about a third in speed and move results by ~1e-12; nnan and ninf are left off because they would let the compiler delete the checks that decide a pixel has failed. reassoc also simplifies exp(log(x)) to x, which drops the NaN of the logarithm of a negative number, so the initial guess checks for a negative ratio itself.

Departures from the reference

• Defaults follow the paper's code listing, not its text. The text describes 15 initial steps, steps of 0.1 and 0.75, a smoothness of 4, a tolerance of 1e-4, and picking the lower chi squared of the two trial steps. The code it lists, which produced its results, uses 5, (0.1, 0.5), 8, 0.1, and interpolates between the two steps. The docstring says so.
• floor, the least intensity the initial guess assumes in a channel, is a keyword that can be given in any units. The reference floors at a bare 0.01 in whatever units its data are in. That is the default here too, so it works for DN/s, DN/(pix s), counts or plain numbers.
• A negative uncertainty fails the pixel. The reference takes it as positive. NaN and zero fail the pixel in both.
• No exposure times. The reference takes counts and exposure times; this takes rates, so the uncertainty model stays with the caller.

Tests

The tests use no AIA data. Each pixel is independent, so testing needs only responses and intensities shaped like AIA's, and a copy of the algorithm to compare against:

• Responses: six Gaussians in log T, at the peaks of the AIA channels, with a cool secondary peak on the hottest one, as 94 and 131 have.
• DEMs: the random DEMs of Section 3.1 of the paper (sums of five log-normals), with Poisson noise.
• Reference: simple_reg_dem transcribed line for line into the test file, with numpy.linalg in place of scipy.linalg, so neither EMToolKit nor scipy is a test dependency.

They cover:

• agreement with the reference to 1e-6, for the default parameters and two other sets;
• recovery of known DEMs without noise;
• arbitrary pixel axes, and an uncertainty broadcast against them;
• a pixel with a NaN, zero or negative uncertainty failing, without disturbing its neighbours;
• negative responses failing every pixel as in the reference, compiled or not (this test fails on the kernel before the reassoc fix);
• a numerically singular system failing at the same step as in the reference, and the factorization itself;
• pixels sharing a chunk, and so its work arrays, giving bitwise the same results as pixels with one each, after a failed pixel too;
• the floor: converted from other units, matching the reference with the same floor, and changing the result where channels fall below it;
• intensities per pixel giving the same result, plain numbers in and out, and uncertain intensities;
• every validation error.

The comparisons against EMToolKit and on real AIA data above were made outside the suite.

CI

Coverage cannot see inside a compiled function, so the coverage job now sets NUMBA_DISABLE_JIT=1 and runs the kernel as plain Python. The test matrix runs it compiled. Coverage of utu.dem is 100%.

Adds numba as a dependency, and intersphinx for scipy and numba.

🤖 Generated with Claude Code

https://claude.ai/code/session_01XR6qYUm91Mfxnsdgo3tcyg

[read-the-docs-community]

Documentation build overview

📚 utu | 🛠️ Build #34889616 | 📁 Comparing 6773ba7 against latest (3121255)

  🔍 Preview build  

10 files changed · + 3 added · ± 7 modified

+ Added

• _autosummary/utu.dem.html
• _autosummary/utu.dem.plowman.html
• _modules/utu/dem/_plowman.html

± Modified

• genindex.html
• py-modindex.html
• _autosummary/utu.html
• _autosummary/utu.rotation.html
• _autosummary/utu.spectrum.lines.html
• _modules/index.html
• tutorials/spectrum.html

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (3121255) to head (6773ba7).

Additional details and impacted files

@@            Coverage Diff             @@
##              main        #9    +/-   ##
==========================================
  Coverage   100.00%   100.00%            
==========================================
  Files           13        17     +4     
  Lines          501       966   +465     
==========================================
+ Hits           501       966   +465     

☔ 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]

Code review (Claude Code, /code-review xhigh) of 2e17e20

These findings come from an automated review of this branch at 2e17e20, diffed against main. The reviewer reproduced most of them with scratch scripts; the ones it verified directly are marked (verified). Nobody has triaged them yet. They are ordered by severity. The 16 new tests pass, both compiled and with NUMBA_DISABLE_JIT=1. Splitting pixels into chunks gives bitwise-identical results, so no state leaks between pixels today.

1. reassoc fast-math folds exp(log(x)) to x, so a failed pixel returns chi2 = NaN, not -1

utu/dem/_plowman_kernel.py#L24, together with L167-L170

With fastmath={"contract", "reassoc"}, LLVM rewrites np.exp(np.log(x)) as x. When the initial-guess ratio is negative, s is NaN but exp_s is a finite negative number. The Cholesky guard then passes, and the pixel iterates on garbage.

• Reproduction: descending temperatures, or any negative rvec.
  • The reference and the NUMBA_DISABLE_JIT coverage job give chi2 = [-1, -1, -1, -1].
  • The compiled kernel gives [nan, nan, nan, nan].
  • A one-line reassoc njit returns (nan, -2.0) for log and exp of -2 (verified).
• Consequences:
  • Callers that flag failures with chi2 == -1 miss these pixels.
  • The coverage job cannot see the bug, because it runs uncompiled.
  • The _flags docstring says only nnan/ninf could remove the failure checks, which is wrong.

2. The default floor=0.01 * u.DN / u.s breaks most unit systems

utu/dem/_plowman.py#L88

Any intensity whose unit doesn't convert to DN/s raises UnitConversionError when no floor is passed. That includes the SSW/aiapy-standard DN/(pix s), counts, and plain numbers.

• Reproduction: intensity in DN/(pix s) with responses in DN cm^5/(pix s) fails at L287 with "'DN / s' and 'DN / (pix s)' are not convertible" (verified).
• Suggested fix: default floor=None, meaning 0.01 in the intensity's own units, which is what the reference does.

3. The temperature grid is never validated, and unitless temperatures are silently treated as kelvin

utu/dem/_plowman.py#L286

The function doesn't check that the grid is increasing, positive, and at least 2 points long. Each of these fails silently:

• Descending response.inputs: every pixel returns a NaN DEM and a NaN chi2 with no exception (verified).
• Unitless log10 T inputs, which are common in DEM codes: read as K, giving chi2 of about 20-118 and a meaningless DEM with no warning (verified).
• smoothness=0, or a single temperature: gives a NaN regmat and fails every pixel.

4. The "intensity has no channel axis" check can be bypassed through uncertainty

utu/dem/_plowman.py#L265

The check runs on the shape after broadcasting with uncertainty. An intensity without the channel axis therefore passes whenever the uncertainty has one, and its single value is used for every channel.

• Reproduction: intensity axes ('pixel',) with uncertainty axes ('pixel', 'channel') raises no error. Every channel gets the same intensity, and the DEMs come back with chi2 of about 0.001-0.03 (verified).
• The error message also says it checks intensity alone.

5. Uncertain intensities crash

utu/dem/_plowman.py#L29

The signature accepts any na.AbstractScalar, but _ndarray calls .ndarray_aligned, which UncertainScalarArray lacks.

• Reproduction: passing na.UncertainScalarArray(...) as the intensity raises AttributeError: 'UncertainScalarArray' object has no attribute 'ndarray_aligned' (verified).
• Suggested fix: treat the distribution axis as one more pixel axis.

6. An intensity carrying axis_temperature runs the full inversion before failing

utu/dem/_plowman.py#L273

shape_pixel isn't checked for axis_temperature.

• Reproduction: intensity axes ('temperature', 'channel') invert every pixel, about 12 s for a full AIA frame. Only then does L326 raise "Each axis name must be unique" (verified).
• Suggested fix: reject this up front.

7. A negative uncertainty is silently used as its absolute value

utu/dem/_plowman.py#L162

The docstring says a non-positive uncertainty gives chi2 = -1 and a NaN DEM.

• What actually happens: uncertainty = -sigma gives exactly the same chi2 (0.9957, ...) as +sigma (verified). Only an uncertainty of exactly zero fails.
• Consequence: a sign error in the uncertainty produces a normal-looking fit, with no failure flag.

8. A unitless uncertainty is accepted when intensity has units

utu/dem/_plowman.py#L284

_ndarray assumes plain numbers are already in unit. The docstring says the uncertainty must be "in the same units".

• Reproduction: an intensity in DN/s with a bare-ndarray uncertainty runs without error (verified). If that ndarray is in DN/min, or is a fraction, the fit is weighted wrongly.
• Asymmetry: the reverse case, a unitless intensity with a Quantity uncertainty, raises.

9. The test helper _reference doesn't match the reference implementation it transliterates

utu/dem/_tests/test_plowman.py#L109

The helper calls np.linalg.cholesky, which raises LinAlgError on a finite matrix that isn't positive-definite. The real simple_reg_dem instead catches every cho_factor exception and breaks out of the loop.

• Consequence: a pixel whose linearized system becomes singular crashes the test instead of being compared. One way this happens is exp(s) underflowing so that amat is about regmat.
• Gap: the kernel's _cholesky returning False on a finite matrix is never checked against the reference behavior.

10. Reusing work arrays across pixels within a chunk is never tested

utu/dem/_plowman.py#L307

The tests use at most 40 pixels, while num_chunks = min(num_pixel, 64 * num_threads). Every chunk therefore holds one pixel, so the main path at AIA scale (one chunk, many pixels) never runs in either CI mode.

• Risk: a future edit that reads a work array before writing it, such as the upper triangle of a or solution, would leak state between neighboring pixels in production but pass every test.
• Suggested test: pass more pixels than chunks, or call the kernel with num_chunks=1.

11. test_plowman_recovers asserts looser bounds than its docstring promises

utu/dem/_tests/test_plowman.py#L234

	Docstring	Assertion
Peak position	within a quarter of a decade	< 0.3 dex
Emission measure	within 40 percent	rtol=0.5

A regression that moves peaks by 0.28 dex or emission measures by 48% still passes.

12. A mismatched temperature count fails deep inside _matrices

utu/dem/_plowman.py#L257

Nothing checks that response.outputs has as many temperatures as response.inputs.

• Reproduction: 30 input temperatures with 31 outputs raises "matmul: Input operand 1 has a mismatch in its core dimension 0 ... (size 30 is different from 31)" (verified).
• Suggested fix: raise a clear ValueError here, like the other shape checks do.

13. A per-channel uncertainty is expanded into a full copy

utu/dem/_plowman.py#L32

_ndarray makes every broadcast input contiguous, so a single uncertainty per channel becomes a full (pixel, channel) float64 array before the kernel runs.

• Cost: about 805 MB extra for a 4096^2 x 6 AIA cube.
• Suggested fix: pass the stride-0 broadcast view (numba accepts 'A'-layout arrays), or pass a per-channel errors array.

14. unit or u.dimensionless_unscaled is repeated three times

utu/dem/_plowman.py#L314

The pattern appears at L31, L287 and L314-L316, and re-implements na.unit_normalized. na.unit_normalized(intensity) and na.unit_normalized(response.outputs) would give the same result; the None check for unitless output still needs na.unit once.

• The output unit isn't normalized either: an intensity in DN/min gives a DEM in s / (min cm5).

15. _matrices uses raw numpy instead of named-arrays

utu/dem/_plowman.py#L60

The group's guideline is to use named_arrays rather than numpy directly, except inside the packages named-arrays itself depends on (ndfilters, colorsynth, regridding). utu is not one of those. The numba kernel itself has to take ndarrays, but _matrices is ordinary wrapper code built from np.diag, np.concatenate and np.matmul.

🤖 Generated with Claude Code