Add utu.dem.emission_measure, the emission measure each temperature of a piecewise-linear DEM stands for

[roytsmart]

This adds utu.dem.emission_measure(dem, axis_temperature), the emission measure each temperature of a piecewise-linear DEM stands for, so that intensities predicted from a DEM of plowman are the ones its fit found.

Why

plowman takes the DEM D and each response R to be linear in log T between the temperatures, and fits the intensities with Rᵀ M D, where M is the mass matrix of those triangle functions. Code that predicts an intensity from the DEM afterwards usually sums D × R × step instead, which leaves out the off-diagonal terms of M. Three places do this now:

• The DEM tutorial of hinode: the ratio of the XRT intensity a DEM predicts to the one observed.
• The ESIS jet paper: the same ratio, in figures/_dem_regions.py. On its four regions, the sum is 0.3% to 1.4% away from Rᵀ M D.
• esis: esis.data.synth.scene_dem (Add synthetic scenes from the DEM of AIA images, and one for ESIS-II esis-mission/esis#109) wrote its own private _emission_measure, so that the AIA intensities of its synthetic scenes are the ones the fit found. This PR moves that function here.

What it does

emission_measure returns M D: Δ(Dᵢ₋₁ + 4Dᵢ + Dᵢ₊₁)/6 inside an even grid, and Δ(2D₀ + D₁)/6 and Δ(Dₙ₋₂ + 2Dₙ₋₁)/6 at the two ends.

• Grids: it takes uneven grids, a grid of its own for each pixel, and any other axes of the DEM, such as the two axes of a map.
• Uncertain arrays: an uncertain DEM, or uncertain temperatures, is weighted by the nominal values and by each sample alike.
• Memory: the result is accumulated in place, so the peak is twice the size of the result (5 times, before the review).
• Errors: it raises a ValueError if:
  • the inputs lack axis_temperature;
  • there are fewer than two temperatures;
  • the temperatures aren't positive or don't increase;
  • the outputs have a different number of temperatures from the inputs.
• Use with synthesize: the result can be given to utu.spectrum.synthesize as its emission measure per bin.

The docstring covers two limits:

• Ranges: summed over a range of temperatures, the result isn't the integral of the DEM over that range, because each temperature's weight reaches a step either side of it. The docstring shows the code for that integral: FunctionArray.integrate over log T, from the i-th temperature to the j-th, whose trapezoid rule is exact for a piecewise-linear DEM (as in Integrate the DEM tutorial's temperature ranges by the trapezoid rule solar-dynamics-observatory#30).

• Interpolation: it integrates the functions as interpolated linearly, which is off by about Δ²/8 times their curvature, not the functions themselves. For smooth functions the grid resolves, the plain sum of the product times the step is much closer to the true integral:

  • In the function's own example (a DEM 2 steps wide), it gives 9.30 DN/s against 8.754 for the integral of the two Gaussians, 6% high. The plain sum matches the integral to rounding error, and the example shows both numbers.
  • On the Gaussians of plowman's docstring example, the plain sum matches their integral to 0.000%, while linear interpolation is 14% high in the channels at 14.9 DN/s and 49% high in the two outer channels, which are nearly dark.

  So the function is for predicting intensities consistently with plowman's fit. plowman's example keeps its sum, and the docstring says to simulate test intensities from the functions themselves.

The docstrings of plowman (its returned DEM) and synthesize (its emission_measure) now point to the function.

Tests

utu/dem/_tests/test_emission_measure.py, 16 tests:

• Constant DEM: the result sums to the value times the span of log T, the ends carry half a step, and a DEM without the axis of temperature broadcasts along it.
• Integral: weighted by a random response, the result equals the integral of the product of the two linear interpolants on a grid of 200,001 points, to 1e-8. This runs on an even and an uneven grid.
• Mass matrix: the result equals the mass matrix plowman builds in _matrices, times the DEM, for every pixel of a 2 × 3 map, to 1e-14.
• Per-pixel grids: two pixels with different temperatures are each weighted by the mass matrix of their own grid.
• Uncertain arrays: an uncertain DEM (the nominal and three hand-picked samples), and uncertain temperatures (each sample a shifted grid), are each weighted by the matching mass matrix.
• Consistency with plowman: the χ² recomputed from the function's intensities equals the χ² plowman reports, to 1e-6. The test simulates its intensities from the functions themselves, as the docstring advises.
• Errors: seven invalid DEMs raise.

I checked that the tests can fail. On the final code:

Mutant	Tests that fail
The lumped mass matrix (the trapezoid rule)	7 of 16
Each end's own value in place of the other end's	7 of 16
One neighbor dropped	9 of 16

On the first version, the rectangle rule failed 7 of 10. I also checked that np.allclose on uncertain arrays compares their samples, and that the DEM passed in is left unchanged. The whole suite passes (124 tests), and black and ruff are clean. All three examples in the docstring run.

Code review

An xhigh review found no errors in the math. Its findings, all reproduced and fixed in 3fe870a and 4793426:

• A zero temperature gave NaN and inf, and a negative one the wrong message.
• A length-1 temperature axis in the outputs broadcast silently.
• The docstring was wrong about the trapezoid rule ("leaves out" the off-diagonal terms) and vague about the ends.
• The docstring understated the error of linear interpolation.
• It didn't show how to integrate over a range.
• The peak memory was 5 times the result.
• Some cases had no tests.

The review also pointed out that the mass matrix was written twice, here and in plowman's _matrices. 4793426 defines it once, in a private _mass, which gives the diagonal and the off-diagonal with the reference's arithmetic in the same order.

• emission_measure weights the DEM with those two vectors.
• _matrices builds its matrix from them, and it is bit for bit the matrix it built before. I checked this on the 0.05-dex grid, on an AIA-like grid and on a random uneven one.
• So a change to the basis now reaches both functions.

Not changed:

• Units: temperatures must convert to kelvin without equivalencies, as plowman requires of its responses.

After this

• Release utu 0.5.0.
• esis drops its private copy and requires utu~=0.5.
• PRs for the hinode tutorial (range maps with integrate, the XRT ratio with this function) and for the jet paper's ratio.

🤖 Generated with Claude Code

https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa

[read-the-docs-community]

Documentation build overview

📚 utu | 🛠️ Build #35198312 | 📁 Comparing 4793426 against latest (e0f889e)

  🔍 Preview build  

9 files changed · + 2 added · ± 7 modified

+ Added

• _autosummary/utu.dem.emission_measure.html
• _modules/utu/dem/_emission_measure.html

± Modified

• genindex.html
• _autosummary/utu.dem.html
• _autosummary/utu.dem.plowman.html
• _autosummary/utu.spectrum.synthesize.html
• _modules/index.html
• _modules/utu/dem/_plowman.html
• _modules/utu/spectrum/_synthesis.html

[codecov]

Codecov Report

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

Additional details and impacted files

@@            Coverage Diff             @@
##              main       #15    +/-   ##
==========================================
  Coverage   100.00%   100.00%            
==========================================
  Files           19        21     +2     
  Lines         1285      1414   +129     
==========================================
+ Hits          1285      1414   +129     

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