XRT DEM Iterative Solver Module#367
Conversation
|
Ensure the function checks for minimum required inputs (filters, observed intensities, etc.) and validates that all observation values are non-zero before proceeding. |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #367 +/- ##
===========================================
- Coverage 87.35% 36.28% -51.07%
===========================================
Files 26 40 +14
Lines 1819 4379 +2560
Branches 106 334 +228
===========================================
Hits 1589 1589
- Misses 184 2744 +2560
Partials 46 46 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
2041e0a to
bd4687b
Compare
| util | ||
| xrt_dem_iterative | ||
| image_correction |
There was a problem hiding this comment.
| util | |
| xrt_dem_iterative | |
| image_correction | |
| xrt_dem_iterative | |
| image_correction | |
| util |
|
|
||
| .. automodule:: xrtpy.response.tools | ||
| :members: | ||
| :undoc-members: |
There was a problem hiding this comment.
🎗️ I want to look up :undoc-members:
| ************************* | ||
| xrtpy.xrt_dem_iterative | ||
| ************************* | ||
|
|
There was a problem hiding this comment.
| ************************* | |
| xrtpy.xrt_dem_iterative | |
| ************************* | |
| *********************** | |
| xrtpy.xrt_dem_iterative | |
| *********************** | |
| @@ -0,0 +1,6 @@ | |||
| ************************* | |||
There was a problem hiding this comment.
🎗️ Does pre-commit have sphinx autoformatting whatever it was? 🤔
| from packaging.version import Version | ||
|
|
||
| # Ensure Sphinx imports the local xrtpy checkout (repo root), not site-packages | ||
| sys.path.insert(0, str(Path(__file__).resolve().parents[1])) |
There was a problem hiding this comment.
🎗️ Was there a new recommended Python way of doing this? 🤔
There was a problem hiding this comment.
I'm wondering if this would be better in a subdirectory maybe? 📂
|
|
||
| Hinode/XRT is well suited for DEM analysis because it observes the corona | ||
| through multiple broadband filters, each sensitive to different temperature | ||
| ranges. By combining these channels, we can infer a DEM(T) that reproduces |
There was a problem hiding this comment.
| ranges. By combining these channels, we can infer a DEM(T) that reproduces | |
| ranges. By combining these channels, we can infer a :math:`\mathrm{DEM}(T)` that reproduces |
| =================================== | ||
| Differential Emission Measure (DEM) | ||
| =================================== | ||
|
|
There was a problem hiding this comment.
| =================================== | |
| Differential Emission Measure (DEM) | |
| =================================== | |
| ============================== | |
| Differential Emission Measures | |
| ============================== | |
| Introduction | ||
| ------------ | ||
|
|
||
| This page describes the XRTpy iterative DEM solver and the inputs and outputs needed to run it. |
There was a problem hiding this comment.
| This page describes the XRTpy iterative DEM solver and the inputs and outputs needed to run it. | |
| This page describes the XRTpy iterative differential emission measure (DEM) solver and the inputs and outputs needed to run it. |
|
|
||
| This page describes the XRTpy iterative DEM solver and the inputs and outputs needed to run it. | ||
|
|
||
| The differential emission measure (DEM) describes how much plasma is present |
There was a problem hiding this comment.
| The differential emission measure (DEM) describes how much plasma is present | |
| The DEM describes how much plasma is present |
| This page describes the XRTpy iterative DEM solver and the inputs and outputs needed to run it. | ||
|
|
||
| The differential emission measure (DEM) describes how much plasma is present | ||
| in the solar corona as a function of temperature. It is a key diagnostic for |
There was a problem hiding this comment.
| in the solar corona as a function of temperature. It is a key diagnostic for | |
| in the solar corona as a function of temperature along a line of sight. DEMs are key diagnostics for |
| 1. Observed channels (filters) | ||
| ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | ||
| * Type: ``list`` of ``str`` | ||
| * Description: Names of the filters used in the observation, for example ``"Al-mesh"`` or ``"Be-thin"``. |
There was a problem hiding this comment.
Either list all of the filters here or link to where the filters are defined? (Maybe a glossary entry?) Linking would reduce duplication of text. 🤔
|
|
||
|
|
| responses = generate_temperature_responses( | ||
| filters, | ||
| "2012-07-10T12:03:20", | ||
| ) | ||
|
|
There was a problem hiding this comment.
🎗️ We should consider whether responses should be included as part of the interface defined in sunkit-dem.
|
|
||
| .. rubric:: Notes | ||
|
|
||
| This mirrors the "regular logT grid" used by the IDL routine ``xrt_dem_iterative2.pro``. |
There was a problem hiding this comment.
| This mirrors the "regular logT grid" used by the IDL routine ``xrt_dem_iterative2.pro``. | |
| This mirrors the "regular logT grid" option used by the IDL routine :file:`xrt_dem_iterative2.pro`. |
Or "keyword"? I can't remember IDL. 😅
| The regular :math:`\log_{10}(T)` grid (dimensionless). | ||
|
|
||
| T : `~astropy.units.Quantity` | ||
| The linear temperature grid in Kelvin. |
There was a problem hiding this comment.
| The linear temperature grid in Kelvin. | |
| The linear temperature grid in kelvin. |
|
|
||
|
|
||
|
|
||
|
|
| Comparison with IDL | ||
| ------------------- |
There was a problem hiding this comment.
I'm wondering about how to qualitatively (and perhaps also quantitatively?) describe the differences that users should expect.
Like...for temperatures where there isn't much plasma, I think there were higher relative errors, if I remember correctly.
| I_i = \int DEM(T)\, R_i(T)\, dT | ||
|
|
||
| where :math:`R_i(T)` is the temperature response function of the filter and | ||
| :math:`DEM(T)` describes the amount of emitting plasma as a function of |
There was a problem hiding this comment.
| :math:`DEM(T)` describes the amount of emitting plasma as a function of | |
| :math:`\mathrm{DEM}(T)` describes the amount of emitting plasma as a function of |
| SolarSoft/IDL routine `xrt_dem_iterative2.pro <https://hesperia.gsfc.nasa.gov/ssw/hinode/xrt/idl/util/xrt_dem_iterative2.pro>`__: | ||
|
|
||
| * Uses a regular :math:`\log_{10}(T)` grid. | ||
| * Represents :math:`\log_{10}(DEM)` at a set of spline knots. |
There was a problem hiding this comment.
| * Represents :math:`\log_{10}(DEM)` at a set of spline knots. | |
| * Represents :math:`\log_{10}(\mathrm{DEM})` at a set of spline knots. |
| The values shown above correspond to the solver defaults and are written | ||
| out here to illustrate which parameters can be tuned. You can adjust these | ||
| to suit your specific analysis needs. This mirrors the flexibility of the | ||
| IDL routine ``xrt_dem_iterative2.pro``. |
There was a problem hiding this comment.
| IDL routine ``xrt_dem_iterative2.pro``. | |
| IDL routine :file:`xrt_dem_iterative2.pro`. |
| monte_carlo_runs=0, | ||
| ) | ||
|
|
||
| x.validate_inputs() |
There was a problem hiding this comment.
Should validate_inputs() be changed to _validate_inputs()? Is it something that needs to be user-facing?
| minimum_bound_temperature=5.5, | ||
| maximum_bound_temperature=7.5, | ||
| logarithmic_temperature_step_size=0.1, |
There was a problem hiding this comment.
| minimum_bound_temperature=5.5, | |
| maximum_bound_temperature=7.5, | |
| logarithmic_temperature_step_size=0.1, | |
| **temperature_kwargs, |
And define temperature_kwargs = {"minimum_bound_temperature": 5.5, ...} outside of the function, so that we can use keyword argument unpacking here?
|
|
||
| x.create_logT_grid() | ||
|
|
||
| assert (x.logT[0], x.logT[-1], x.dlogT) == pytest.approx((5.5, 7.5, 0.1)) |
There was a problem hiding this comment.
| assert (x.logT[0], x.logT[-1], x.dlogT) == pytest.approx((5.5, 7.5, 0.1)) | |
| np.testing.assert_allclose( | |
| [x.logT[0], x.logT[-1], x.dlogT], | |
| [5.5, 7.5, 0.1], | |
| ) |
| dem.validate_inputs() | ||
|
|
||
|
|
||
| def test_create_logT_grid(): |
There was a problem hiding this comment.
This would be a good use case for pytest fixtures to avoid having repeated assertions, but it's probably not worth the effort changing here.
| x.create_logT_grid() | ||
| x._interpolate_responses_to_grid() | ||
|
|
||
| dem, modeled, chi2, result = x._solve_single_dem( |
There was a problem hiding this comment.
🎗️ I'm wondering if this would also be useful to think about adding upstream to the interface in sunkit-dem, especially if many DEM solvers are iterative. 🤔
| temperature_responses=responses, | ||
| ) | ||
| x.solve() | ||
| x.plot_dem() |
There was a problem hiding this comment.
🎗️ Is there a way to do something like a regression test for figures?
| def test_init_rejects_negative_max_iterations(): | ||
| """max_iterations=-100 must raise ValueError. | ||
| Expected: ValueError because max_iterations must be a positive integer. | ||
| """ | ||
| filters = ["Al-poly", "Ti-poly"] | ||
| intensities = np.array([500.0, 800.0]) | ||
| responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") | ||
|
|
||
| with pytest.raises(ValueError, match="max_iterations must be a positive integer"): | ||
| XRTDEMIterative(filters, intensities, responses, max_iterations=-100) | ||
|
|
||
|
|
||
| def test_init_rejects_zero_normalization_factor(): | ||
| """normalization_factor=0 must raise ValueError. | ||
| Expected: ValueError because normalization_factor must be positive. | ||
| """ | ||
| filters = ["Al-poly", "Ti-poly"] | ||
| intensities = np.array([500.0, 800.0]) | ||
| responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") | ||
|
|
||
| with pytest.raises(ValueError, match="normalization_factor must be a positive"): | ||
| XRTDEMIterative(filters, intensities, responses, normalization_factor=0) |
There was a problem hiding this comment.
We might be able to use @pytest.mark.parametrize for some of these. 🤔
| filters = ["Al-poly", "Ti-poly", "Be-thin"] | ||
| intensities = np.array([100.0, 200.0, 300.0], dtype=float) |
There was a problem hiding this comment.
Some of these be defined outside of the test functions to reduce duplication.
2 participants
Summary
This PR introduces
XRTDEMIterative, a Python implementation of the IDLroutine
xrt_dem_iterative2.profor estimating Differential EmissionMeasures (DEMs) from Hinode/XRT multi-filter observations.
What's Included
Core Solver (
dem_solver.py)XRTDEMIterativeclass with full DEM solving pipelinefitting via
lmfitmax(0.03 * intensity, 2 DN/s/pix), matching IDL default behaviormonte_carlo_runsparameterintensities, missing errors)
summary()method for diagnostic outputPlotting (
dem_plotting.py)plot_dem()— base DEM solution visualizationplot_dem_mc()— base DEM + Monte Carlo ensemble overlayTests
response interpolation, spline system, residuals, solver, and
Monte Carlo behavior
Documentation
dem_overview.rst— introduction to DEM analysis with XRTpyxrtpy.xrt_dem_iterativedocs/changelog/0.6.0.rstValidation -
xrt_dem_iterative2.proacross multiplesolar observation cases using Monte Carlo comparison