From 69dcfb247ee351dc255635a62f9724a12324958d Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 16 Oct 2026 10:34:43 -0600 Subject: [PATCH 1/3] Add `utu.dem.emission_measure`, the emission measure each temperature of a piecewise-linear DEM stands for plowman fits the intensities with R^T M D, where M is the mass matrix of the triangle functions the DEM and the responses are interpolated with. emission_measure returns M D, so a response times it, summed over temperature, is the intensity the fit found, and it can be given to synthesize as the emission measure of each bin. It moves esis's private _emission_measure here, takes uneven grids and uncertain DEMs, and checks its axis and temperatures. The docstrings of plowman and synthesize point to it. Its own says it is not the integral over a range of temperatures, which FunctionArray.integrate gives, and that it integrates the interpolated functions rather than the functions themselves. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa --- utu/dem/__init__.py | 2 + utu/dem/_emission_measure.py | 154 +++++++++++++++++ utu/dem/_plowman.py | 5 +- utu/dem/_tests/test_emission_measure.py | 213 ++++++++++++++++++++++++ utu/spectrum/_synthesis.py | 5 +- 5 files changed, 377 insertions(+), 2 deletions(-) create mode 100644 utu/dem/_emission_measure.py create mode 100644 utu/dem/_tests/test_emission_measure.py diff --git a/utu/dem/__init__.py b/utu/dem/__init__.py index de88ed2..77127d4 100644 --- a/utu/dem/__init__.py +++ b/utu/dem/__init__.py @@ -1,7 +1,9 @@ """Differential emission measures, how much plasma there is at each temperature.""" +from ._emission_measure import emission_measure from ._plowman import plowman __all__ = [ + "emission_measure", "plowman", ] diff --git a/utu/dem/_emission_measure.py b/utu/dem/_emission_measure.py new file mode 100644 index 0000000..b06c462 --- /dev/null +++ b/utu/dem/_emission_measure.py @@ -0,0 +1,154 @@ +"""The emission measure each temperature of a piecewise-linear DEM stands for.""" + +from typing import cast + +import astropy.units as u +import named_arrays as na +import numpy as np + +__all__ = [ + "emission_measure", +] + + +def emission_measure( + dem: na.FunctionArray[na.AbstractScalar, na.AbstractScalar], + axis_temperature: str, +) -> na.AbstractScalar: + r""" + The emission measure each temperature of a DEM stands for, if the DEM is + linear in :math:`\log_{10} T` between its temperatures, as the DEMs of + :func:`utu.dem.plowman` are. + + A DEM which is linear in :math:`\log_{10} T` between its temperatures is + a sum of triangle functions, one at each temperature. If a temperature + response or a contribution function is taken to be linear between the + same temperatures too, as :func:`utu.dem.plowman` takes its responses, + the integral of the product of the two is :math:`R^T M D`, where + :math:`D` and :math:`R` are their values at the temperatures and + :math:`M` is the mass matrix of the triangles. This function returns + :math:`M D`. At a temperature whose neighbors are :math:`\Delta_-` below + and :math:`\Delta_+` above it, in :math:`\log_{10} T`, that is + + .. math:: + + (MD)_i = \frac{\Delta_- (D_{i-1} + 2 D_i) + \Delta_+ (2 D_i + D_{i+1})}{6}, + + where the :math:`\Delta` beyond either end is zero. For an even step + :math:`\Delta`, it is :math:`\Delta (D_{i-1} + 4 D_i + D_{i+1}) / 6` + inside the grid, and half as much at either end. + + This is the integral :func:`utu.dem.plowman` fits the intensities with, + so the intensity of a channel predicted from these emission measures and + its response, ``(response.outputs * result).sum(axis_temperature)``, is + the intensity the inversion found. Neither the sum of the DEM times the + step nor the trapezoid rule applied to the product is: both leave out + the off-diagonal terms of :math:`M`, which matter wherever the DEM and + the response both change across a step. The result can also be given to + :func:`utu.spectrum.synthesize` as the emission measure of each bin. + + This is the integral of the functions as interpolated between the + temperatures, not of the functions themselves, and the two differ where + either curves on the scale of one step, as the far tails of narrow peaks + do. For smooth functions which the grid resolves, the plain sum of their + product times the step can be the closer of the two to their integral. + So to simulate the intensities of a known DEM, to test an inversion + against, integrate the functions themselves rather than this model of + them. + + Summed over a range of temperatures, the result is not the integral of + the DEM over that range, since each temperature stands for plasma up to + one step either side of it. For that, use + :meth:`~named_arrays.AbstractFunctionArray.integrate` on the DEM with + :math:`\log_{10} T` as its inputs, whose trapezoid rule is exact for a + DEM which is linear between its temperatures. + + Parameters + ---------- + dem + The DEM per unit :math:`\log_{10} T`, as :func:`utu.dem.plowman` + returns it: its inputs are temperatures along `axis_temperature`, + increasing and with units of temperature. Its outputs may be + uncertain, and have any other axes, such as the two axes of a map. + axis_temperature + The name of the axis along the temperatures of `dem`. + + Returns + ------- + : + The emission measure each temperature stands for, in the units of + the outputs of `dem`, along `axis_temperature` and every other axis + of `dem`. + + Examples + -------- + A Gaussian DEM holding :math:`10^{27}\,\mathrm{cm^{-5}}` of plasma in + all, and the intensity it gives in a channel whose response peaks at a + hotter temperature. + + .. jupyter-execute:: + + import astropy.units as u + import named_arrays as na + import numpy as np + import utu + + temperature = 10 ** na.linspace(5.5, 7, axis="temperature", num=31) * u.K + logt = np.log10(temperature / u.K) + + width = 0.1 + dem = na.FunctionArray( + inputs=temperature, + outputs=1e27 / u.cm**5 * np.exp(-np.square((logt - 6.2) / width) / 2) + / (np.sqrt(2 * np.pi) * width), + ) + + emission_measure = utu.dem.emission_measure(dem, axis_temperature="temperature") + + # the total emission measure, 1e27 cm^-5 + emission_measure.sum("temperature") + + .. jupyter-execute:: + + response = 1e-25 * np.exp(-np.square((logt - 6.5) / 0.15)) * u.DN * u.cm**5 / u.s + + # the intensity, as utu.dem.plowman models it + (response * emission_measure).sum("temperature") + """ + temperature = dem.inputs + shape = na.shape(temperature) + if axis_temperature not in shape: + raise ValueError( + f"the inputs of `dem` have no axis {axis_temperature!r}, " + f"got {tuple(shape)}" + ) + if shape[axis_temperature] < 2: + raise ValueError("`dem` needs at least two temperatures") + + logt = np.log10(temperature.to(u.K).value) + + lower = {axis_temperature: slice(None, -1)} + upper = {axis_temperature: slice(1, None)} + + # the width of each interval between neighboring temperatures + step = logt[upper] - logt[lower] + if not np.all(step > 0): + raise ValueError("the temperatures of `dem` must be increasing") + + d = na.broadcast_to(dem.outputs, na.shape_broadcasted(temperature, dem.outputs)) + + # the width of the interval below and above each temperature, zero + # beyond either end + zero = 0 * step[{axis_temperature: slice(0, 1)}] + below = na.concatenate([zero, step], axis=axis_temperature) + above = na.concatenate([step, zero], axis=axis_temperature) + + # the neighbor below and above each temperature, any value beyond either + # end, since its interval there is zero + first = d[{axis_temperature: slice(0, 1)}] + last = d[{axis_temperature: slice(-1, None)}] + d_below = na.concatenate([first, d[lower]], axis=axis_temperature) + d_above = na.concatenate([d[upper], last], axis=axis_temperature) + + result = (below * (d_below + 2 * d) + above * (2 * d + d_above)) / 6 + return cast(na.AbstractScalar, result) diff --git a/utu/dem/_plowman.py b/utu/dem/_plowman.py index cf92e77..7a6b58a 100644 --- a/utu/dem/_plowman.py +++ b/utu/dem/_plowman.py @@ -188,7 +188,10 @@ def plowman( :math:`\log_{10} T`, in units of ``intensity`` over the units of the outputs of ``response``: :math:`\mathrm{cm^{-5}}` for responses in :math:`\mathrm{DN\,cm^5\,s^{-1}}` and intensities in - :math:`\mathrm{DN\,s^{-1}}`. + :math:`\mathrm{DN\,s^{-1}}`. The fit takes it to be linear in + :math:`\log_{10} T` between the temperatures, and + :func:`utu.dem.emission_measure` weights it the same way, so that + the intensities it predicts from the DEM are the ones the fit found. chi2 The reduced :math:`\chi^2` of each pixel. It is :math:`-1` where the iteration failed before its first step, in which case the DEM is not a diff --git a/utu/dem/_tests/test_emission_measure.py b/utu/dem/_tests/test_emission_measure.py new file mode 100644 index 0000000..3808362 --- /dev/null +++ b/utu/dem/_tests/test_emission_measure.py @@ -0,0 +1,213 @@ +import astropy.units as u +import named_arrays as na +import numpy as np +import pytest + +import utu +from utu.dem._plowman import _matrices + +_logt = np.round(np.arange(5.5, 7.0 + 1e-9, 0.05), 2) + +_logt_uneven = np.array([5.5, 5.7, 5.8, 5.85, 5.9, 6.0, 6.2, 6.5, 7.0]) +"""An uneven grid of temperatures, closest together near 1 MK.""" + + +def _mass(logt: np.ndarray) -> np.ndarray: + """The mass matrix :func:`utu.dem.plowman` fits with, from its own code.""" + return _matrices(logt, np.eye(logt.size), smoothness=1)[0] + + +def _temperature(logt: np.ndarray) -> na.ScalarArray: + """The temperatures of a DEM, from their logarithms.""" + return na.ScalarArray(10**logt * u.K, axes="temperature") + + +def _dem( + logt: np.ndarray, + dem: np.ndarray, + axes: tuple[str, ...] = ("temperature",), +) -> na.FunctionArray: + """A DEM in units of :math:`\\mathrm{cm^{-5}}` from plain arrays.""" + return na.FunctionArray( + inputs=_temperature(logt), + outputs=na.ScalarArray(dem / u.cm**5, axes=axes), + ) + + +@pytest.mark.parametrize( + argnames="outputs", + argvalues=[ + na.ScalarArray(np.full(_logt.size, 3e27) / u.cm**5, axes="temperature"), + # without the axis of temperature, broadcast along it + na.ScalarArray(3e27 / u.cm**5), + ], +) +def test_emission_measure_constant(outputs: na.ScalarArray) -> None: + """ + A constant DEM holds its value times the span of :math:`\\log_{10} T`, + and its ends stand for half a step. + """ + dem = na.FunctionArray(inputs=_temperature(_logt), outputs=outputs) + + result = utu.dem.emission_measure(dem, axis_temperature="temperature") + + unit = 1 / u.cm**5 + assert result.shape == dict(temperature=_logt.size) + assert np.allclose(result.sum("temperature"), 3e27 * unit * 1.5) + assert np.allclose(result[{"temperature": 0}], 3e27 * unit * 0.025) + assert np.allclose(result[{"temperature": 5}], 3e27 * unit * 0.05) + assert np.allclose(result[{"temperature": -1}], 3e27 * unit * 0.025) + + +@pytest.mark.parametrize( + argnames="logt", + argvalues=[_logt, _logt_uneven], +) +def test_emission_measure_integral(logt: np.ndarray) -> None: + """ + Weighted by a response, it is the integral of the product of the + response and the DEM, each linear between the temperatures. + """ + rng = np.random.default_rng(0) + dem = rng.uniform(1e26, 1e28, logt.size) + response = rng.uniform(0, 1e-25, logt.size) + + result = utu.dem.emission_measure(_dem(logt, dem), axis_temperature="temperature") + intensity = (na.ScalarArray(response, axes="temperature") * result).sum( + "temperature" + ) + + # the trapezoid rule on a grid fine enough to be exact to about 1e-10 + fine = np.linspace(logt[0], logt[-1], 200_001) + product = np.interp(fine, logt, dem) * np.interp(fine, logt, response) + expected = np.trapezoid(product, fine) / u.cm**5 + + assert np.allclose(intensity.ndarray, expected, rtol=1e-8) + + +def test_emission_measure_matrix() -> None: + """ + It is the mass matrix of :func:`utu.dem.plowman` times the DEM, for + every pixel of a map. + """ + rng = np.random.default_rng(1) + dem = rng.uniform(1e26, 1e28, (2, 3, _logt_uneven.size)) + + result = utu.dem.emission_measure( + _dem(_logt_uneven, dem, axes=("y", "x", "temperature")), + axis_temperature="temperature", + ) + + expected = na.ScalarArray( + dem @ _mass(_logt_uneven) / u.cm**5, + axes=("y", "x", "temperature"), + ) + assert result.shape == dict(y=2, x=3, temperature=_logt_uneven.size) + assert np.allclose(result, expected, rtol=1e-14) + + +def test_emission_measure_uncertain() -> None: + """The nominal DEM and each sample of its distribution are weighted alike.""" + axis = na.UncertainScalarArray.axis_distribution + nominal = 1e27 * np.exp(-np.square((_logt - 6.2) / 0.2) / 2) + samples = np.stack( + [ + 0.5 * nominal, + 1e27 * np.exp(-np.square((_logt - 6.0) / 0.1) / 2), + np.full(_logt.size, 2e26), + ] + ) + dem = na.FunctionArray( + inputs=_temperature(_logt), + outputs=na.UncertainScalarArray( + nominal=na.ScalarArray(nominal / u.cm**5, axes="temperature"), + distribution=na.ScalarArray( + samples / u.cm**5, + axes=(axis, "temperature"), + ), + ), + ) + + result = utu.dem.emission_measure(dem, axis_temperature="temperature") + + mass = _mass(_logt) + expected = na.UncertainScalarArray( + nominal=na.ScalarArray(nominal @ mass / u.cm**5, axes="temperature"), + distribution=na.ScalarArray( + samples @ mass / u.cm**5, + axes=(axis, "temperature"), + ), + ) + assert isinstance(result, na.UncertainScalarArray) + assert np.allclose(result, expected, rtol=1e-14) + + +def test_emission_measure_plowman() -> None: + """ + The intensities it predicts from the DEM :func:`utu.dem.plowman` finds + give the reduced :math:`\\chi^2` the inversion reports. + """ + peaks = np.array([5.7, 5.95, 6.15, 6.3, 6.5, 6.9]) + response = na.FunctionArray( + inputs=_temperature(_logt), + outputs=na.ScalarArray( + ndarray=1e-26 + * np.exp(-np.square((_logt[:, None] - peaks) / 0.12)) + * u.DN + * u.cm**5 + / u.s, + axes=("temperature", "channel"), + ), + ) + + amplitude = na.ScalarArray(np.array([1e27, 3e27, 1e28]), axes="pixel") + logt = na.ScalarArray(_logt, axes="temperature") + dem_true = na.FunctionArray( + inputs=response.inputs, + outputs=amplitude * np.exp(-np.square((logt - 6.2) / 0.15) / 2) / u.cm**5, + ) + intensity = ( + response.outputs * utu.dem.emission_measure(dem_true, "temperature") + ).sum("temperature") + + # noise, so that the fit stops short of the data + rng = np.random.default_rng(0) + uncertainty = 0.05 * intensity + noise = na.ScalarArray(rng.normal(size=(6, 3)), axes=("channel", "pixel")) + intensity = intensity + noise * uncertainty + + dem, chi2 = utu.dem.plowman( + intensity=intensity, + uncertainty=uncertainty, + response=response, + axis_channel="channel", + axis_temperature="temperature", + ) + + model = (response.outputs * utu.dem.emission_measure(dem, "temperature")).sum( + "temperature" + ) + expected = np.square((intensity - model) / uncertainty).mean("channel") + + assert np.all(chi2 >= 0) + assert np.allclose(chi2, expected, rtol=1e-6) + + +@pytest.mark.parametrize( + argnames="dem,match", + argvalues=[ + ( + na.FunctionArray( + inputs=na.ScalarArray(10**_logt * u.K, axes="t"), + outputs=na.ScalarArray(np.ones(_logt.size) / u.cm**5, axes="t"), + ), + "no axis 'temperature'", + ), + (_dem(_logt[:1], np.ones(1)), "at least two temperatures"), + (_dem(_logt[::-1], np.ones(_logt.size)), "must be increasing"), + ], +) +def test_emission_measure_invalid(dem: na.FunctionArray, match: str) -> None: + """A DEM whose temperatures cannot be integrated over is an error.""" + with pytest.raises(ValueError, match=match): + utu.dem.emission_measure(dem, axis_temperature="temperature") diff --git a/utu/spectrum/_synthesis.py b/utu/spectrum/_synthesis.py index 5d733f5..79c6aec 100644 --- a/utu/spectrum/_synthesis.py +++ b/utu/spectrum/_synthesis.py @@ -94,7 +94,10 @@ def synthesize( This is the emission measure of each bin, not the differential emission measure: a DEM per unit :math:`\log_{10} T`, such as :func:`utu.dem.plowman` returns, has the same units, and has to be - multiplied by the width of each bin in :math:`\log_{10} T` first. + converted first: multiplied by the width of each bin in + :math:`\log_{10} T`, or, to integrate it the way + :func:`utu.dem.plowman` fits the intensities, weighted with + :func:`utu.dem.emission_measure`. velocity The vertices of the bins of velocity, along `axis_velocity`. Positive velocities are away from the observer, to longer wavelengths. From 3fe870a1e9205e7edcb63c69cb5a9e0605e3f6ef Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 16 Oct 2026 11:09:09 -0600 Subject: [PATCH 2/3] Check the temperatures of emission_measure, accumulate it in place, and correct its docstring From the xhigh code review of #15: - Raise for temperatures that are not positive, which gave NaN and inf, or the wrong message, and for outputs whose axis of temperature does not match the inputs, which a length-1 axis broadcast past. - Accumulate the result in place: 2 times the size of the result at most, rather than 5. - The docstring gives the weights of the two ends, says the sum times the step and the trapezoid rule move the off-diagonal terms of M onto its diagonal rather than leave them out, and shows how to integrate the DEM over a range of temperatures. - It also says how much linear interpolation costs: in its own example, emission_measure gives an intensity 6% above the integral of the two Gaussians, which the plain sum matches, and the example shows both. - Tests of per-pixel grids, uncertain temperatures and the new errors, and the chi2 test simulates its intensities from the functions themselves, as the docstring advises. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa --- utu/dem/_emission_measure.py | 104 +++++++++++++++--------- utu/dem/_tests/test_emission_measure.py | 84 ++++++++++++++++++- 2 files changed, 148 insertions(+), 40 deletions(-) diff --git a/utu/dem/_emission_measure.py b/utu/dem/_emission_measure.py index b06c462..4c83618 100644 --- a/utu/dem/_emission_measure.py +++ b/utu/dem/_emission_measure.py @@ -34,42 +34,60 @@ def emission_measure( (MD)_i = \frac{\Delta_- (D_{i-1} + 2 D_i) + \Delta_+ (2 D_i + D_{i+1})}{6}, - where the :math:`\Delta` beyond either end is zero. For an even step - :math:`\Delta`, it is :math:`\Delta (D_{i-1} + 4 D_i + D_{i+1}) / 6` - inside the grid, and half as much at either end. + where the :math:`\Delta` beyond either end is zero, so that the first + and last temperatures stand for :math:`\Delta_+ (2 D_0 + D_1) / 6` and + :math:`\Delta_- (D_{n-2} + 2 D_{n-1}) / 6`. For an even step + :math:`\Delta`, a temperature inside the grid stands for + :math:`\Delta (D_{i-1} + 4 D_i + D_{i+1}) / 6`. This is the integral :func:`utu.dem.plowman` fits the intensities with, so the intensity of a channel predicted from these emission measures and its response, ``(response.outputs * result).sum(axis_temperature)``, is the intensity the inversion found. Neither the sum of the DEM times the - step nor the trapezoid rule applied to the product is: both leave out - the off-diagonal terms of :math:`M`, which matter wherever the DEM and - the response both change across a step. The result can also be given to - :func:`utu.spectrum.synthesize` as the emission measure of each bin. - - This is the integral of the functions as interpolated between the - temperatures, not of the functions themselves, and the two differ where - either curves on the scale of one step, as the far tails of narrow peaks - do. For smooth functions which the grid resolves, the plain sum of their - product times the step can be the closer of the two to their integral. - So to simulate the intensities of a known DEM, to test an inversion - against, integrate the functions themselves rather than this model of - them. + step nor the trapezoid rule applied to the product is: both weight each + temperature by the sum of its row of :math:`M`, which moves the + off-diagonal terms of :math:`M` onto its diagonal, and that matters + wherever the DEM and the response both change across a step. The result + can also be given to :func:`utu.spectrum.synthesize` as the emission + measure of each bin. + + This is the integral of the functions as interpolated linearly between + the temperatures, not of the functions themselves, and linear + interpolation is off by about :math:`\Delta^2 / 8` times the curvature + of each. For smooth functions which the grid resolves, the plain sum of + their product times the step is much the closer of the two to their + integral: in the example below, this function gives an intensity 6% + higher than the integral of the two Gaussians, which the plain sum + matches to rounding error. So use this function to predict intensities + consistently with the fit of :func:`utu.dem.plowman`, and, to simulate + the intensities of a known DEM to test an inversion against, integrate + the functions themselves rather than this model of them. Summed over a range of temperatures, the result is not the integral of the DEM over that range, since each temperature stands for plasma up to - one step either side of it. For that, use - :meth:`~named_arrays.AbstractFunctionArray.integrate` on the DEM with - :math:`\log_{10} T` as its inputs, whose trapezoid rule is exact for a - DEM which is linear between its temperatures. + one step either side of it. To integrate the DEM from its ``i``-th to + its ``j``-th temperature, give it :math:`\log_{10} T` as its inputs and + use :meth:`~named_arrays.AbstractFunctionArray.integrate`, whose + trapezoid rule is exact for a DEM which is linear between its + temperatures: + + .. code-block:: python + + dem_logt = dem.replace(inputs=np.log10(dem.inputs / u.K)) + dem_logt[{axis_temperature: slice(i, j + 1)}].integrate(axis_temperature) + + A range whose ends fall between the temperatures needs the DEM + interpolated at its ends first. Parameters ---------- dem The DEM per unit :math:`\log_{10} T`, as :func:`utu.dem.plowman` returns it: its inputs are temperatures along `axis_temperature`, - increasing and with units of temperature. Its outputs may be - uncertain, and have any other axes, such as the two axes of a map. + positive, increasing and with units of temperature. Its outputs may + be uncertain, and have any other axes, such as the two axes of a map. + If they vary along `axis_temperature`, they must have one value for + each temperature. axis_temperature The name of the axis along the temperatures of `dem`. @@ -114,6 +132,13 @@ def emission_measure( # the intensity, as utu.dem.plowman models it (response * emission_measure).sum("temperature") + + The plain sum of the product times the step, 6% less, which for these + smooth functions is the integral of the functions themselves. + + .. jupyter-execute:: + + (response * dem.outputs).sum("temperature") * 0.05 """ temperature = dem.inputs shape = na.shape(temperature) @@ -122,10 +147,20 @@ def emission_measure( f"the inputs of `dem` have no axis {axis_temperature!r}, " f"got {tuple(shape)}" ) - if shape[axis_temperature] < 2: + num = shape[axis_temperature] + if num < 2: raise ValueError("`dem` needs at least two temperatures") + num_outputs = na.shape(dem.outputs).get(axis_temperature, num) + if num_outputs != num: + raise ValueError( + f"the outputs of `dem` have {num_outputs} temperatures, " + f"but its inputs have {num}" + ) - logt = np.log10(temperature.to(u.K).value) + temperature = temperature.to(u.K) + if not np.all(temperature > 0 * u.K): + raise ValueError("the temperatures of `dem` must be positive") + logt = np.log10(temperature.value) lower = {axis_temperature: slice(None, -1)} upper = {axis_temperature: slice(1, None)} @@ -137,18 +172,13 @@ def emission_measure( d = na.broadcast_to(dem.outputs, na.shape_broadcasted(temperature, dem.outputs)) - # the width of the interval below and above each temperature, zero - # beyond either end + # Each interval of width step adds step / 3 times the DEM at either end + # to the temperature there, and step / 6 times the DEM at the other end. + # The result is accumulated in place, since a map of DEMs can be large. zero = 0 * step[{axis_temperature: slice(0, 1)}] - below = na.concatenate([zero, step], axis=axis_temperature) - above = na.concatenate([step, zero], axis=axis_temperature) - - # the neighbor below and above each temperature, any value beyond either - # end, since its interval there is zero - first = d[{axis_temperature: slice(0, 1)}] - last = d[{axis_temperature: slice(-1, None)}] - d_below = na.concatenate([first, d[lower]], axis=axis_temperature) - d_above = na.concatenate([d[upper], last], axis=axis_temperature) - - result = (below * (d_below + 2 * d) + above * (2 * d + d_above)) / 6 + width = na.concatenate([step, zero], axis=axis_temperature) + width = width + na.concatenate([zero, step], axis=axis_temperature) + result = d * (width / 3) + result[lower] += d[upper] * (step / 6) + result[upper] += d[lower] * (step / 6) return cast(na.AbstractScalar, result) diff --git a/utu/dem/_tests/test_emission_measure.py b/utu/dem/_tests/test_emission_measure.py index 3808362..dada448 100644 --- a/utu/dem/_tests/test_emission_measure.py +++ b/utu/dem/_tests/test_emission_measure.py @@ -106,6 +106,64 @@ def test_emission_measure_matrix() -> None: assert np.allclose(result, expected, rtol=1e-14) +def test_emission_measure_grid_per_pixel() -> None: + """Each pixel may have its own temperatures, and is weighted by its own.""" + grids = np.stack([_logt_uneven, _logt_uneven + 0.1]) + rng = np.random.default_rng(2) + dem = rng.uniform(1e26, 1e28, grids.shape) + + result = utu.dem.emission_measure( + na.FunctionArray( + inputs=na.ScalarArray(10**grids * u.K, axes=("x", "temperature")), + outputs=na.ScalarArray(dem / u.cm**5, axes=("x", "temperature")), + ), + axis_temperature="temperature", + ) + + expected = na.ScalarArray( + np.stack([d @ _mass(g) for g, d in zip(grids, dem)]) / u.cm**5, + axes=("x", "temperature"), + ) + assert np.allclose(result, expected, rtol=1e-14) + + +def test_emission_measure_uncertain_temperature() -> None: + """ + Uncertain temperatures weight the DEM by the grid of the nominal + temperatures and of each sample. + """ + axis = na.UncertainScalarArray.axis_distribution + grids = np.stack([_logt_uneven + 0.02, _logt_uneven - 0.03]) + dem = np.linspace(1e27, 2e27, _logt_uneven.size) + + result = utu.dem.emission_measure( + na.FunctionArray( + inputs=na.UncertainScalarArray( + nominal=_temperature(_logt_uneven), + distribution=na.ScalarArray( + 10**grids * u.K, + axes=(axis, "temperature"), + ), + ), + outputs=na.ScalarArray(dem / u.cm**5, axes="temperature"), + ), + axis_temperature="temperature", + ) + + expected = na.UncertainScalarArray( + nominal=na.ScalarArray( + dem @ _mass(_logt_uneven) / u.cm**5, + axes="temperature", + ), + distribution=na.ScalarArray( + np.stack([dem @ _mass(g) for g in grids]) / u.cm**5, + axes=(axis, "temperature"), + ), + ) + assert isinstance(result, na.UncertainScalarArray) + assert np.allclose(result, expected, rtol=1e-14) + + def test_emission_measure_uncertain() -> None: """The nominal DEM and each sample of its distribution are weighted alike.""" axis = na.UncertainScalarArray.axis_distribution @@ -166,9 +224,9 @@ def test_emission_measure_plowman() -> None: inputs=response.inputs, outputs=amplitude * np.exp(-np.square((logt - 6.2) / 0.15) / 2) / u.cm**5, ) - intensity = ( - response.outputs * utu.dem.emission_measure(dem_true, "temperature") - ).sum("temperature") + # the intensities of the functions themselves, as the docstring advises, + # which the sum times the step integrates for these smooth functions + intensity = (response.outputs * dem_true.outputs).sum("temperature") * 0.05 # noise, so that the fit stops short of the data rng = np.random.default_rng(0) @@ -205,6 +263,26 @@ def test_emission_measure_plowman() -> None: ), (_dem(_logt[:1], np.ones(1)), "at least two temperatures"), (_dem(_logt[::-1], np.ones(_logt.size)), "must be increasing"), + ( + na.FunctionArray( + inputs=na.ScalarArray( + np.array([0, 1e5, 1e6]) * u.K, axes="temperature" + ), + outputs=na.ScalarArray(np.ones(3) / u.cm**5, axes="temperature"), + ), + "must be positive", + ), + ( + na.FunctionArray( + inputs=na.ScalarArray( + np.array([-1, 1e5, 1e6]) * u.K, axes="temperature" + ), + outputs=na.ScalarArray(np.ones(3) / u.cm**5, axes="temperature"), + ), + "must be positive", + ), + (_dem(_logt, np.ones(1)), "have 1 temperatures, but its inputs have 31"), + (_dem(_logt, np.ones(5)), "have 5 temperatures, but its inputs have 31"), ], ) def test_emission_measure_invalid(dem: na.FunctionArray, match: str) -> None: From 4793426a36beb16f80b657d979dcbe500f18c974 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 16 Oct 2026 12:10:59 -0600 Subject: [PATCH 3/3] Define the mass matrix once, for emission_measure and plowman _mass gives the diagonal and off-diagonal of the mass matrix of the triangle functions, with the arithmetic of the reference in the same order. emission_measure weights the DEM with them, and plowman's _matrices builds its matrix from them, bit for bit the one it built before on even and uneven grids. A change to the basis now reaches both. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017xqpht97bF4FWsACAA4wQa --- utu/dem/_emission_measure.py | 53 ++++++++++++++++++++++++++---------- utu/dem/_plowman.py | 18 +++++++----- 2 files changed, 49 insertions(+), 22 deletions(-) diff --git a/utu/dem/_emission_measure.py b/utu/dem/_emission_measure.py index 4c83618..a4ef2dd 100644 --- a/utu/dem/_emission_measure.py +++ b/utu/dem/_emission_measure.py @@ -162,23 +162,46 @@ def emission_measure( raise ValueError("the temperatures of `dem` must be positive") logt = np.log10(temperature.value) - lower = {axis_temperature: slice(None, -1)} - upper = {axis_temperature: slice(1, None)} - - # the width of each interval between neighboring temperatures - step = logt[upper] - logt[lower] - if not np.all(step > 0): + diagonal, off_diagonal = _mass(logt, axis_temperature) + if not np.all(off_diagonal > 0): raise ValueError("the temperatures of `dem` must be increasing") d = na.broadcast_to(dem.outputs, na.shape_broadcasted(temperature, dem.outputs)) - # Each interval of width step adds step / 3 times the DEM at either end - # to the temperature there, and step / 6 times the DEM at the other end. - # The result is accumulated in place, since a map of DEMs can be large. - zero = 0 * step[{axis_temperature: slice(0, 1)}] - width = na.concatenate([step, zero], axis=axis_temperature) - width = width + na.concatenate([zero, step], axis=axis_temperature) - result = d * (width / 3) - result[lower] += d[upper] * (step / 6) - result[upper] += d[lower] * (step / 6) + # M D, accumulated in place since a map of DEMs can be large: each + # temperature's own term, then the term of its neighbor above and of its + # neighbor below + lower = {axis_temperature: slice(None, -1)} + upper = {axis_temperature: slice(1, None)} + result = d * diagonal + result[lower] += d[upper] * off_diagonal + result[upper] += d[lower] * off_diagonal return cast(na.AbstractScalar, result) + + +def _mass( + logt: na.AbstractScalar, + axis: str, +) -> tuple[na.AbstractScalar, na.AbstractScalar]: + r""" + The diagonal and the off-diagonal of the mass matrix of the triangle + functions on a grid of :math:`\log_{10} T` along `axis`. + + The matrix is tridiagonal and symmetric. This is the one definition of it + in the package: :func:`emission_measure` weights a DEM with it, and + :func:`utu.dem.plowman` fits the intensities with it. Its arithmetic is + that of the reference implementation of :func:`utu.dem.plowman`, in the + same order, since a rounding error in the matrix is carried through every + step of the inversion of every pixel. + """ + step = logt[{axis: slice(1, None)}] - logt[{axis: slice(None, -1)}] + + # the width of the interval above and below each temperature, zero + # beyond either end + zero = 0 * step[{axis: slice(0, 1)}] + above = na.concatenate([step, zero], axis=axis) + below = na.concatenate([zero, step], axis=axis) + + diagonal = (above + below) * 2.0 / 6.0 + off_diagonal = step / 6.0 + return cast(na.AbstractScalar, diagonal), cast(na.AbstractScalar, off_diagonal) diff --git a/utu/dem/_plowman.py b/utu/dem/_plowman.py index 7a6b58a..e7d8751 100644 --- a/utu/dem/_plowman.py +++ b/utu/dem/_plowman.py @@ -6,6 +6,8 @@ import named_arrays as na import numpy as np +from ._emission_measure import _mass + __all__ = [ "plowman", ] @@ -70,20 +72,22 @@ def _matrices( The DEM and the responses are both taken to be piecewise linear between the temperatures of the grid, so the mapping is the response times the - mass matrix of those triangle functions, and the regularization is their + mass matrix of those triangle functions, from :func:`_mass`, which + :func:`utu.dem.emission_measure` shares, and the regularization is their stiffness matrix. Both are tridiagonal, and both are written with the operations of the reference in the same order, since a rounding error in - a matrix is carried through every step of every pixel. That, and the - kernel taking plain arrays, is why this is :mod:`numpy` rather than + a matrix is carried through every step of every pixel. The kernel takes + plain arrays, so the rest of this is :mod:`numpy` rather than :mod:`named_arrays`. """ num_temperature, num_channel = response.shape dt = logt[1:] - logt[:-1] - left = np.concatenate([dt, [0]]) - right = np.concatenate([[0], dt]) - mass = np.diag((left + right) * 2.0) / 6.0 - mass += np.diag(dt, k=1) / 6.0 + np.diag(dt, k=-1) / 6.0 + diagonal, off_diagonal = _mass(na.ScalarArray(logt, axes="t"), axis="t") + diagonal = diagonal.ndarray + off_diagonal = off_diagonal.ndarray + mass = np.diag(diagonal) + mass += np.diag(off_diagonal, k=1) + np.diag(off_diagonal, k=-1) inverse_left = np.concatenate([1.0 / dt, [0]]) inverse_right = np.concatenate([[0], 1.0 / dt])