diff --git a/docs/about_xrt.rst b/docs/about_xrt.rst index c9848a872..575c6901b 100644 --- a/docs/about_xrt.rst +++ b/docs/about_xrt.rst @@ -8,8 +8,8 @@ Hinode ====== .. image:: _static/images/hinode_satellite.png - :alt: Hinode Satellite - :align: center + :alt: Hinode Satellite + :align: center Hinode is a joint mission involving the space agencies of Japan, the United States, Europe, and the United Kingdom. It is depicted in the *illustration shown above*. @@ -25,9 +25,9 @@ The X-Ray Telescope =================== .. image:: _static/images/XRT_composite_image_full_disk_14February2015.png - :alt: XRT Composite Image - :align: center - :scale: 50% + :alt: XRT Composite Image + :align: center + :scale: 50% The X-Ray Telescope (XRT), depicted as a long linear black tube on the Hinode spacecraft is a crucial instrument for observing the solar corona's most intense regions, with temperatures ranging from 1,000,000 to 10,000,000 Kelvin. The image below is a synoptic composite from February 14, 2015, created using the Al-Mesh/Be-Thin/Al-Med filters. @@ -36,7 +36,7 @@ For a comprehensive overview of XRT's mission and capabilities, please visit the .. tip:: - Visit the `XRT Picture of the Week`_ and the `Hinode-XRT YouTube`_ page for captivating visual content showcasing the XRT's solar observations. + Visit the `XRT Picture of the Week`_ and the `Hinode-XRT YouTube`_ page for captivating visual content showcasing the XRT's solar observations. XRT uses two sequentially positioned filter wheels, as shown in the diagram below, where each wheel houses a variety of filters. By rotating these wheels, scientists can select different filters to study the Sun in different wavelengths, thereby enhancing the resolution and quality of solar images. @@ -73,8 +73,8 @@ The existing filters are structured as follows: The process is the same for all XRT filter channels. .. image:: _static/images/XRT_filter_wheels_Sun_View_Diagram.png - :alt: Diagram of the XRT Filter Wheels - :align: center + :alt: Diagram of the XRT Filter Wheels + :align: center Data Products ************* @@ -99,8 +99,8 @@ The XRT software was originally created in the Interactive Data Language (IDL). .. note:: - Please note that the `SolarSoft XRT Analysis Guide`_ does not serve as a guide for using XRTpy. - It focuses solely on the analysis of XRT data using the IDL software. + Please note that the `SolarSoft XRT Analysis Guide`_ does not serve as a guide for using XRTpy. + It focuses solely on the analysis of XRT data using the IDL software. .. _Hinode-XRT YouTube: https://www.youtube.com/user/xrtpow .. _Interactive Data Language: https://www.l3harrisgeospatial.com/Software-Technology/IDL diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 9c26a4eb2..59a01cde0 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -171,3 +171,23 @@ @ARTICLE{velasquez:2024 doi = {10.21105/joss.06396}, url = {https://doi.org/10.21105/joss.06396} } + +@ARTICLE{weber:2004, + author = {{Weber}, M.~A. and {DeLuca}, E.~E. and {Golub}, L. and {Sette}, A.~L.}, + title = "{Temperature diagnostics with multichannel imaging telescopes}", + journal = "{IAU Symposium 223: Multi-Wavelength Investigations of Solar Activity}", + year = 2004, + pages = {321--328}, + publisher = {Cambridge University Press}, + doi = {10.1017/S1743921304006088} +} + +@ARTICLE{golub:2004, + author = {{Golub}, L. and {DeLuca}, E.~E. and {Sette}, A. and {Weber}, M.}, + title = "{DEM analysis with the X-Ray Telescope (XRT) for Solar-B}", + journal = "{ASP Conference Series, The Solar-B Mission and the Outlook for Space-Based Solar Physics}", + year = 2004, + volume = {325}, + pages = {217--222}, + publisher = {Astronomical Society of the Pacific} +} diff --git a/docs/conf.py b/docs/conf.py index 941c35134..c3681e6e0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -70,7 +70,7 @@ ogp_image = "https://raw.githubusercontent.com/HinodeXRT/xrtpy/main/docs/_static/images/xrtpy_logo.png" ogp_use_first_image = True ogp_description_length = 160 -ogp_custom_meta_tags = ('',) +#ogp_custom_meta_tags = ('',) # Suppress warnings about overriding directives as we overload some of the # doctest extensions. @@ -230,6 +230,7 @@ "feedback_communication*": [], "contributing*": [], "code_of_conduct*": [], + "dem_overview*": [], } # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, diff --git a/docs/dem_overview.rst b/docs/dem_overview.rst new file mode 100644 index 000000000..1c9ba5eaf --- /dev/null +++ b/docs/dem_overview.rst @@ -0,0 +1,264 @@ +.. _xrtpy-dem-overview: + +=================================== +Differential Emission Measure (DEM) +=================================== + +.. contents:: + :local: + :depth: 2 + +Introduction +------------ + +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 +understanding coronal heating, solar flares, and the thermal structure of +active regions. + +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 temperature distribution +DEM(T) that explains the observed X-ray intensities. + + +DEM in XRTpy +------------ +XRTpy provides a Python implementation of the iterative spline fitting method +originally available in IDL as `xrt_dem_iterative2.pro `__. +The core solver is implemented in :class:`xrtpy.xrt_dem_iterative.dem_solver.XRTDEMIterative`. + +Conceptually, the solver: + 1. Builds a regular grid in log10(T) between user-specified bounds. + 2. Interpolates the filter temperature responses onto that grid. + 3. Represents log10(DEM) as a spline in log10(T). + 4. Uses least-squares fitting (via ``lmfit``) to adjust the spline values so that the modeled filter intensities match the observed intensities. + 5. Optionally performs Monte Carlo runs by perturbing the observed intensities with their errors and re-solving the DEM many times to estimate uncertainties. + +This approach mirrors the structure and behavior of the IDL routine while providing +a modern, fully open-source implementation in Python. + +Required inputs +--------------- +The DEM workflow requires three main input pieces: + +1. Observed channels (filters) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +* Type: ``str`` or ``list`` of ``str`` +* Description: Names of the filters used in the observation, for example ``"Al-mesh"`` or ``"Be-thin"``. +* These must correspond to valid XRT filters and must match the provided temperature responses one-to-one. + +2. Observed intensities +~~~~~~~~~~~~~~~~~~~~~~~ +* Type: array-like +* Units: DN/s (normalized per pixel) +* Description: Measured intensities in each filter channel. +* Length must match the number of filters. + + +3. Temperature response functions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +* Type: ``list`` of :class:`xrtpy.response.TemperatureResponseFundamental` +* Units: DN s\ :sup:`-1` pix\ :sup:`-1` cm\ :sup:`5` +* Description: Instrument response as a function of temperature for each filter, matching the order of the filters. +* Can be generated using :func:`xrtpy.response.tools.generate_temperature_responses`. + + + + +Example +------- +A simple example with two filters: + + +.. code-block:: python + + from xrtpy.response.tools import generate_temperature_responses + + filters = ["Al-poly", "Ti-poly"] + responses = generate_temperature_responses( + filters, + "2012-10-27T00:00:00", + ) + + +Overview of the XRTDEMIterative API +----------------------------------- +The main entry point is :class:`xrtpy.xrt_dem_iterative.dem_solver.XRTDEMIterative`. + + + +Constructor +~~~~~~~~~~~ +.. code-block:: python + + from xrtpy.xrt_dem_iterative import XRTDEMIterative + + dem_solver = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + monte_carlo_runs=0, + max_iterations=2000, + ) + + # Solve for the DEM + dem = dem_solver.solve() # returns the DEM array, also stored in dem_solver.dem + + # Plot the DEM + dem_solver.plot_dem() + + +Enabling Monte Carlo error estimates +------------------------------------ +To estimate uncertainties, you can enable Monte Carlo iterations. The solver +will perturb the observed intensities by their errors and re-solve the DEM +for each realization. + +.. code-block:: python + + N_mc = 100 # number of Monte Carlo runs + + dem_solver = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + monte_carlo_runs=N_mc, + max_iterations=2000, + ) + + dem_solver.solve() + + # Monte Carlo DEM Plot + dem_solver.plot_dem_mc() # base DEM plus Monte Carlo curves + +The arrays ``dem_solver.mc_dem``, ``dem_solver.mc_chisq``, +``dem_solver.mc_base_obs``, and ``dem_solver.mc_mod_obs`` are then available +for custom analysis. + +Comparison with IDL +------------------- +The Python solver is designed to closely follow the logic of the +SolarSoft/IDL routine `xrt_dem_iterative2.pro `__: + + + + +* Uses a regular log10(T) grid. +* Represents log10(DEM) at a set of spline knots. +* Uses a least-squares algorithm to minimize chi-square. +* Supports Monte Carlo noise realizations for uncertainty estimation. + + +Small numerical differences can arise due to: + +* Different interpolation choices (for example, cubic splines from SciPy). +* Differences in optimization libraries (lmfit versus IDL MPFIT). +* Floating-point rounding and platform-specific details. + +Within these limits, the Python implementation is intended to produce +results that are consistent with the IDL tool. + + + +Mathematical background +----------------------- +The DEM inversion problem is ill posed. For each filter channel i, +the observed intensity :math:`I_i` is related to the DEM through: + +.. math:: + + I_i = \int DEM(T) \, R_i(T) \, dT + +where :math:`R_i(T)` is the temperature response function for the filter, +and :math:`DEM(T)` is the unknown thermal distribution. + +Since the number of temperature bins typically exceeds the number of +observed channels, the inversion does not have a unique solution. The +XRTpy solver uses a forward-fitting approach: + +1. Assume a parametric form for log10(DEM(T)) using spline knots. +2. Compute model intensities: + + .. math:: + + I_i^{model} = \sum_j DEM(T_j)\, R_i(T_j)\, T_j\, \Delta(\ln T) + +3. Adjust the spline values to minimize: + + .. math:: + + \chi^2 = \sum_i \left[ + \frac{I_i^{model} - I_i^{obs}}{\sigma_i} + \right]^2 + +Here :math:`\sigma_i` are the observational uncertainties. Smoothness and +the low number of spline knots help regularize the solution. + + +Monte Carlo iterations perturb the observed intensities: + +.. math:: + + I_i^{(k)} = I_i^{obs} + \mathcal{N}(0, \sigma_i) + +and re-fit the DEM for each realization k. The spread in the resulting DEM +curves provides an estimate of the uncertainty in DEM(T). + + +Extended example with options +----------------------------- +Below is an extended example showing more constructor options explicitly. +These values match current defaults but are written out here for clarity. + +.. code-block:: python + + from xrtpy.response.tools import generate_temperature_responses + from xrtpy.xrt_dem_iterative import XRTDEMIterative + + filters = ["Al-poly", "Ti-poly", "Be-thin", "C-poly"] + intensities = [2500.0, 1800.0, 900.0, 450.0] # DN/s + observation_date="2012-10-27T00:00:00" + + responses = generate_temperature_responses( + filters, + observation_date, + ) + + dem_solver = XRTDEMIterative( + observed_channel=filters, # Filter names + observed_intensities=intensities, # Observed intensity values + temperature_responses=responses, # Instrument responses + + # Optional configuration: + intensity_errors=None, # Obs. uncertainties - default: auto-estimated (3%) + minimum_bound_temperature=5.5, # Minimum log T (default: 5.5) + maximum_bound_temperature=8.0, # Maximum log T (default: 8.0) + logarithmic_temperature_step_size=0.1, # Bin width in log T (default: 0.1) + monte_carlo_runs=100, # # of Monte Carlo runs (default: none) + max_iterations=2000, # Solver max iterations (default: 2000) + normalization_factor=1e21, # Normalization saling factor (default: 1e21) + ) + + dem_solver.solve() + dem_solver.plot_dem_mc() + +.. note:: + The values shown above correspond to existing defaults in the solver, + but they are written out here to illustrate what can be tuned. + You can adjust these to best suit your analysis needs. + This mirrors the flexibility of the IDL routine ``xrt_dem_iterative2.pro``. + + +.. Acknowledgement +.. --------------- +.. *Development of the DEM solver in XRTpy has been supported in part by +.. a NASA Heliophysics Tools and Methods (HTM) program grant (ROSES-2025, +.. element B.20). This effort reflects the ongoing transition of DEM +.. capabilities from legacy IDL routines into modern, open-source Python +.. tools for the solar physics community.* + +References +---------- +- Golub, L., et al. (2004), Solar Physics, 243, 63. :cite:p:`golub:2004` +- Weber, M. A., et al. (2004), Astrophysical Journal, 605, 528. :cite:p:`weber:2004` diff --git a/docs/gallery/data_processing/README.txt b/docs/gallery/data_processing/README.txt new file mode 100644 index 000000000..f34326aac --- /dev/null +++ b/docs/gallery/data_processing/README.txt @@ -0,0 +1,11 @@ +Advanced Data Processing +========================= + +This section covers advanced image correction and diagnostic tools. + +Examples here show how to: +- Deconvolve images to correct for blurring +- Identify and remove light leaks +- Estimate temperature and emission measure using multi-filter observations + +These are core XRTpy functionalities for high-level solar data analysis. \ No newline at end of file diff --git a/docs/gallery/data_processing/__init__.py b/docs/gallery/data_processing/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/docs/gallery/data_processing/deconvolving.py b/docs/gallery/data_processing/deconvolving.py new file mode 100644 index 000000000..503a36e03 --- /dev/null +++ b/docs/gallery/data_processing/deconvolving.py @@ -0,0 +1,45 @@ +""" +======================= +Deconvolving XRT Images +======================= + +This example demonstrates deconvolvoing X-Ray Telescope (XRT) images using the +`xrtpy.image_correction.deconvolve` function in XRTpy. +""" + +import matplotlib.pyplot as plt +import sunpy.map +from sunpy.net import Fido +from sunpy.net import attrs as a + +from xrtpy.image_correction import deconvolve + +############################################################################## +# We will search for XRT data from the Virtual Solar Observatory (VSO) and fetch the first result. + +result = Fido.search( + a.Time("2012-06-05 21:58:39", "2012-06-05 21:59:00"), a.Instrument("xrt") +) +data_file = Fido.fetch(result[0]) + +############################################################################## +# Typically most deconvolve routines use the Richardson-Lucy deconvolution algorithm. + +xrt_map = sunpy.map.Map(data_file) +deconv_map = deconvolve(xrt_map) + +############################################################################## +# To see the effects of the deconvolution we plot both the before and after images. + +fig = plt.figure(figsize=(15, 10)) + +ax = fig.add_subplot(121, projection=xrt_map) +xrt_map.plot(axes=ax, title="Original") +ax1 = fig.add_subplot(122, projection=deconv_map) +deconv_map.plot(axes=ax1, title="Deconvolved") + +ax1.coords[1].set_ticks_visible(False) +ax1.coords[1].set_ticklabel_visible(False) +fig.tight_layout() + +plt.show() diff --git a/docs/gallery/data_processing/remove_lightleak.py b/docs/gallery/data_processing/remove_lightleak.py new file mode 100644 index 000000000..51e5c6f4f --- /dev/null +++ b/docs/gallery/data_processing/remove_lightleak.py @@ -0,0 +1,68 @@ +""" +=================== +Removing Light Leak +=================== + +In this example, we show how to remove the light leak (visible stray light) +from XRT synoptic composite images. +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import sunpy.map +from astropy.utils.data import get_pkg_data_path +from astropy.visualization import ImageNormalize, SqrtStretch + +from xrtpy.image_correction import remove_lightleak + +############################################################################## +# This example will be using XRT synoptic data from the first day of summer of 2015. +# This is stored in the ``example_data`` directory of the `xrtpy` package. + +directory = get_pkg_data_path("data/example_data", package="xrtpy.image_correction") +data_file = Path(directory) / "comp_XRT20150621_055911.7.fits" +xrt_map = sunpy.map.Map(data_file) + +############################################################################## +# Removing the light leak from the composite image is done using the `xrtpy.image_correction.remove_lightleak` function. + +lightleak_map = remove_lightleak(xrt_map) + +############################################################################## +# Finally, we plot the original and light leak subtracted images side by side. + +fig = plt.figure(figsize=(12, 6)) + +ax = fig.add_subplot(121, projection=xrt_map) +xrt_map.plot( + axes=ax, + title="Original", + norm=ImageNormalize(vmin=0, vmax=7e3, stretch=SqrtStretch()), +) +ax1 = fig.add_subplot(122, projection=lightleak_map) +lightleak_map.plot( + axes=ax1, + title="Light Leak Subtracted", + norm=ImageNormalize(vmin=0, vmax=7e3, stretch=SqrtStretch()), +) + +ax1.coords[1].set_ticks_visible(False) +ax1.coords[1].set_ticklabel_visible(False) +fig.tight_layout() + +############################################################################## +# They look almost identical, but the light leak has been removed from the second image. +# To confirm this we can plot the difference between the two images. + +diff_data = xrt_map.data - lightleak_map.data +# For this image, the difference is very small. +print(diff_data.min(), diff_data.max()) + +fig = plt.figure() +ax = fig.add_subplot(111) +ax.set_title("Lightleak Difference") +im = ax.imshow(diff_data, origin="lower") +fig.colorbar(im) + +plt.show() diff --git a/docs/gallery/data_processing/temperature_from_filter_ratios.py b/docs/gallery/data_processing/temperature_from_filter_ratios.py new file mode 100644 index 000000000..04d1aa9f0 --- /dev/null +++ b/docs/gallery/data_processing/temperature_from_filter_ratios.py @@ -0,0 +1,70 @@ +""" +================================================ +Calculating the temperature and emission measure +================================================ + +In this example, we will showcase how to use the filter method to calculate +the temperature and emission measure of the X-ray Telescope (XRT) on Hinode. +""" + +import matplotlib.pyplot as plt +import sunpy.map +from astropy.visualization import ImageNormalize, LogStretch +from sunpy.net import Fido +from sunpy.net import attrs as a + +from xrtpy.response import temperature_from_filter_ratio + +############################################################################## +# To start, we will get XRT data via ``sunpy``. +# +# It is important to use images that same size and with the smallest time separation. +# Note that not all filter ratios produce good results. + +query = Fido.search( + a.Time("2011-01-28 01:31:55", "2011-01-28 01:32:05"), a.Instrument("xrt") +) +data_files = Fido.fetch(query) +xrt_map_1 = sunpy.map.Map(data_files[0]) +xrt_map_2 = sunpy.map.Map(data_files[1]) + +############################################################################## +# The `xrtpy.response.temperature_from_filter_ratio` function has several options, mirroring +# the IDL routine xrt_teem.pro in SolarSoft in most respects.A simple call with +# no extra parameters calculates the temperature and (volume) emission measure +# for the two images without any binning or masking of the data. + +T_EM = temperature_from_filter_ratio(xrt_map_1, xrt_map_2) + +############################################################################## +# The output is a namedtuple with attributes ``Tmap``, ``EMmap``, ``Terrmap``, and ``EMerrmap``. +# As with the SolarSoft IDL routine xrt_teem.pro, the output images are logs of the quantities. +# +# ``Tmap`` is the electron temperature, ``EMmap`` is the volume emission measure, ``Terrmap`` +# is a measure of the uncertainties in the temperature determined for each pixel and ``EMerrmap`` +# is the same for the emission measure. + +T_e = T_EM.Tmap + +fig = plt.figure() + +ax = plt.subplot(projection=T_e) +T_e.plot( + title="Derived Temperature", + cmap="inferno", + norm=ImageNormalize(vmin=6, vmax=7, stretch=LogStretch(10)), +) +T_e.draw_limb() +T_e.draw_grid() +plt.colorbar(label="T (K)") +plt.tight_layout() + +plt.show() + +############################################################################## +# If you want to do the same for Level 2 synoptic composite images, you have to use +# `~.make_exposure_map` to generate the exposure maps for the composite images. +# This is then passed to `xrtpy.response.temperature_from_filter_ratio` as the +# ``expmap1`` and ``expmap2`` arguments. +# Otherwise without accounting for the different exposure time per pixel, +# the temperature and emission measure will be incorrect. diff --git a/docs/gallery/image_filtering/README.txt b/docs/gallery/image_filtering/README.txt new file mode 100644 index 000000000..8b6076c10 --- /dev/null +++ b/docs/gallery/image_filtering/README.txt @@ -0,0 +1,6 @@ +Image Filtering & Visualization +=============================== + +This section includes examples that demonstrate how to sort, filter, and visualize XRT image data. + +The tools shown here help users prepare XRT observations for further scientific analysis. \ No newline at end of file diff --git a/docs/gallery/image_filtering/__init__.py b/docs/gallery/image_filtering/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/docs/gallery/image_filtering/sorting_data.py b/docs/gallery/image_filtering/sorting_data.py new file mode 100644 index 000000000..ed4fe0b98 --- /dev/null +++ b/docs/gallery/image_filtering/sorting_data.py @@ -0,0 +1,92 @@ +""" +========================= +Filtering and Visualizing +========================= + +This example provides a simple overview of filtering and visualizing XRT data. +""" + +import astropy.units as u +import matplotlib.pyplot as plt +import sunpy.map +from astropy.visualization import ImageNormalize, SqrtStretch +from sunpy.net import Fido +from sunpy.net import attrs as a + +############################################################################## +# To start we will download a range of XRT data from the Virtual Solar Observatory (VSO). +# The goal is to acquire a large set of files we can sort through and visualize. + +query = Fido.search( + a.Time("2021-05-21 18:51:00", "2021-05-22 00:00:00"), a.Instrument("xrt") +) +print(query) + +############################################################################## +# This query will return a large number of files, this is due to the fact we do not +# specify any additional filters. We can filter the data by specifying additional +# attributes in the query. +# +# For wavelength, we use a range that focuses the data to return only Al-Poly filter images. +# This will cut the results down in half. + +query = Fido.search( + a.Time("2021-05-21 20:51:00", "2021-05-22 00:00:00"), + a.Instrument("xrt"), + a.Wavelength(4 * u.nm, 5 * u.nm), +) +print(query) + +############################################################################## +# Now we will download the data. +# As this is still over 60 files, this process can take some time. + +xrt_files = Fido.fetch(query) + +############################################################################## +# We can now load the data into a `~sunpy.map.MapSequence` and create a animation. + +xrt_seq = sunpy.map.Map(xrt_files, sequence=True) + +fig = plt.figure() +ax = fig.add_subplot(projection=xrt_seq.maps[0]) +ani = xrt_seq.plot( + axes=ax, norm=ImageNormalize(vmin=0, vmax=5e3, stretch=SqrtStretch()) +) + +############################################################################## +# You might notice that there is a jump in the sequence. +# The size of the data and the pointing changes. +# We can exclude these images by filtering the data further. + +xrt_seq_filtered_shape = sunpy.map.Map( + [m for m in xrt_seq if m.data.shape == (384, 384)], sequence=True +) + +fig = plt.figure() +ax = fig.add_subplot(projection=xrt_seq.maps[0]) +ani = xrt_seq_filtered_shape.plot( + axes=ax, norm=ImageNormalize(vmin=0, vmax=5e3, stretch=SqrtStretch()) +) + +############################################################################## +# In fact, `sunpy.map.Map` provides many attributes that can be used to filter the data. +# This provides a lot of flexibility in how you can filter the data for your science objective. +# +# For example, we can filter the data by the exposure time or the detector. + +xrt_seq_filtered_exp_time = sunpy.map.Map( + [m for m in xrt_seq_filtered_shape if m.exposure_time < 0.1 * u.s], sequence=True +) + +fig = plt.figure() +ax = fig.add_subplot(projection=xrt_seq.maps[0]) +ani = xrt_seq_filtered_exp_time.plot( + axes=ax, norm=ImageNormalize(vmin=0, vmax=5e3, stretch=SqrtStretch()) +) + +############################################################################## +# If you want to save this animation to a file, you can use the ``save`` method. +# For more information on how to use this method, `see the matplotlib documentation `__. + +plt.show() diff --git a/docs/gallery/instrument_response/README.txt b/docs/gallery/instrument_response/README.txt new file mode 100644 index 000000000..c22eab767 --- /dev/null +++ b/docs/gallery/instrument_response/README.txt @@ -0,0 +1,10 @@ +Instrument Response & Calibration +================================= + +These examples illustrate how to evaluate and interpret the XRT instrument's characteristics. + +You will learn to work with: +- The instrument's effective area across different filters +- Temperature response functions for solar plasma + +This section is helpful for understanding how the instrument "sees" the Sun at different temperatures. \ No newline at end of file diff --git a/docs/gallery/instrument_response/__init__.py b/docs/gallery/instrument_response/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/docs/gallery/instrument_response/channels.py b/docs/gallery/instrument_response/channels.py new file mode 100644 index 000000000..485fbfaba --- /dev/null +++ b/docs/gallery/instrument_response/channels.py @@ -0,0 +1,103 @@ +""" +============================= +Exploring XRT's Configuration +============================= + +This example explores the X-Ray Telescope (XRT) instrument properties +using XRTpy's `xrtpy.response.Channel`. It provides convenient methods and attributes +to access and analyze various aspects of the XRT instrument configuration. +""" + +import matplotlib.pyplot as plt + +import xrtpy + +############################################################################## +# We begin by defining a filter channel by its common abbreviation. +# In this example we will be exploring the titanium-on-polyimide filter. +# For detailed information about various filter channels and their characteristics, you can refer to :ref:`xrtpy-about-xrt-filters`. +# +# To explore the properties and characteristics of a defined filter channel, we will create a +# `xrtpy.response.Channel`. By passing in the filter name as an input, we can work +# with the properties associated with the titanium-on-polyimide filter. + +channel = xrtpy.response.Channel("Ti-poly") + +############################################################################## +# Now that we have created our channel, we can delve into the XRT instrument and its properties. +# We will start by examining basic information about the XRT instrument. + +print("Selected filter:", channel.name) +print("\nObservatory:", channel.observatory) +print("Instrument:", channel.instrument) + +############################################################################## +# It is important to note that most instrument properties of XRT remain the same +# regardless of the specific filter being used. This means that many characteristics +# and specifications of the XRT instrument, such as its dimensions, +# field of view, and detector properties, are independent of the selected filter. +# +# We can explore various characteristics of the the Charge-Coupled-Device (CCD) +# camera camera, such as its quantum efficiency and pixel size to list a few. + +print(channel.ccd.ccd_name) +print("\nPixel size: ", channel.ccd.ccd_pixel_size) +print("Full well: ", channel.ccd.ccd_full_well) +print("Gain left: ", channel.ccd.ccd_gain_left) +print("Gain right: ", channel.ccd.ccd_gain_right) +print("eV pre electron: ", channel.ccd.ccd_energy_per_electron) + +############################################################################## +# We can explore the XRT entrance filter properties utilizing ``entrancefilter``. + +print(channel.entrancefilter.entrancefilter_name) +print("Material: ", channel.entrancefilter.entrancefilter_material) +print("Thickness: ", channel.entrancefilter.entrancefilter_thickness) +print("Density: ", channel.entrancefilter.entrancefilter_density) + +############################################################################## +# XRT data is recorded through nine X-ray filters, which are implemented using two filter wheels. +# +# By utilizing the ``channel.filter_#`` notation, where ``#`` represents filter wheel 1 or 2, +# we can explore detailed information about the selected XRT channel filter. +# +# It's worth noting that sometimes the other filter will yield the result "Open," as it's not use. +# For more comprehensive information about the XRT filters, you can refer to :ref:`xrtpy-about-xrt-filters`. + +print("Filter Wheel:", channel.filter_2.filter_name) +print("\nFilter material:", channel.filter_2.filter_material) +print("Thickness: ", channel.filter_2.filter_thickness) +print("Density: ", channel.filter_2.filter_density) + +############################################################################## +# We can explore geometry factors in the XRT using ``geometry``. + +print(channel.geometry.geometry_name) +print("\nFocal length:", channel.geometry.geometry_focal_len) +print("Aperture Area:", channel.geometry.geometry_aperture_area) + +############################################################################## +# The XRT is equipped with two mirrors and We can access the properties of these +# mirrors using the ``channel_mirror_#`` notation, where ``#`` represents the +# first or second mirror surface. + +print(channel.mirror_1.mirror_name) +print("Material: ", channel.mirror_1.mirror_material) +print("Density: ", channel.mirror_1.mirror_density) +print("Graze_angle: ", channel.mirror_1.mirror_graze_angle) + +############################################################################## +# Finally we can explore the XRT transmission properties + +plt.figure() + +plt.plot(channel.wavelength, channel.transmission, label=f"{channel.name}") +plt.title(f"{channel.name} filter") +plt.xlabel(r"$\lambda$ [Å]") +plt.ylabel(r"Transmittance") +# The full range goes up to 400 Å, but we will limit it to 80 Å for better visualization +plt.xlim(0, 80) +plt.grid(color="lightgrey") +plt.tight_layout() + +plt.show() diff --git a/docs/gallery/instrument_response/effective_area.py b/docs/gallery/instrument_response/effective_area.py new file mode 100644 index 000000000..f2052cfd4 --- /dev/null +++ b/docs/gallery/instrument_response/effective_area.py @@ -0,0 +1,75 @@ +""" +======================= +Effective Area Analysis +======================= + +In this example, we will explore the effective areas for different XRT filter channels. +Understanding the effective areas is important for accurately interpreting and quantifying the data. +""" + +import matplotlib.pyplot as plt + +import xrtpy + +############################################################################## +# Let us begin by defining a filter channel using its abbreviation. +# For example, if we want to explore the effective area for an aluminum-on-polyimide filter +# channel, we need to specify the relevant abbreviation. + +xrt_filter = "Al-poly" + +############################################################################## +# `~.EffectiveAreaFundamental` allows us to accurately determine the effective area +# based on the specified filter channel, date, and time. + +date_time = "2023-09-22T22:59:59" +eaf = xrtpy.response.EffectiveAreaFundamental(xrt_filter, date_time) + +############################################################################## +# To actually calculate the effective area function we can call :meth:`~xrtpy.response.EffectiveAreaFundamental.effective_area`. + +effective_area = eaf.effective_area() +print("Effective Area:\n", effective_area) + +############################################################################## +# Differences overtime arise from an increase of the contamination layer on the +# CCD which blocks some of the X-rays thus reducing the effective area. +# For detailed information about the calculation of the XRT CCD contaminant layer thickness, +# you can refer to +# `Montana State University `__. +# +# Additional information is provided by +# `Narukage et. al. (2011) `__. + +relative_launch_date_time = "2006-09-22T22:59:59" +eaf_launch = xrtpy.response.EffectiveAreaFundamental( + xrt_filter, relative_launch_date_time +) +launch_effective_area = eaf_launch.effective_area() + +############################################################################## +# Finally, we can plot how the effective area has changed over time. + +plt.figure() + +plt.plot( + eaf.wavelength, + effective_area, + label=f"{date_time}", +) +plt.plot( + eaf.wavelength, + launch_effective_area, + label=f"{relative_launch_date_time}", +) + +plt.title("XRT Effective Area - Al-Poly") +plt.xlabel("Wavelength (Å)") +plt.ylabel("Effective Area ($cm^{2}$)") +plt.legend() +plt.xlim(0, 60) + +plt.grid(color="lightgrey") +plt.tight_layout() + +plt.show() diff --git a/docs/gallery/instrument_response/temperature_response.py b/docs/gallery/instrument_response/temperature_response.py new file mode 100644 index 000000000..031d5153c --- /dev/null +++ b/docs/gallery/instrument_response/temperature_response.py @@ -0,0 +1,88 @@ +""" +==================== +Temperature Response +==================== + +In this example, we will explore the temperature response of the filters on XRT. +The temperature response provides important information on how XRT responds to +the different temperatures of X-ray emissions. +""" + +import matplotlib.pyplot as plt +import numpy as np + +import xrtpy + +############################################################################## +# A filter channel is defined by its common abbreviation, which represents +# a specific type of filter used to modify the X-ray radiation observed. +# In this example, we will explore the carbon-on-polyimide filter (abbreviated as "C-poly"). + +xrt_filter = "C-poly" + +############################################################################## +# `~.TemperatureResponseFundamental` provides the functions and properties for +# calculating the temperature response. + +date_time = "2023-09-22T21:59:59" +tpf = xrtpy.response.TemperatureResponseFundamental( + xrt_filter, date_time, abundance_model="Photospheric" +) + +############################################################################## +# To calculate the temperature response,we can do the following: + +temperature_response = tpf.temperature_response() +print("Temperature Response:\n", temperature_response) + +############################################################################## +# We will now visualize the temperature response function using CHIANTI. +# These temperatures are of the plasma and are independent of the channel filter. +# +# We use the log of the these temperatures, to enhance the visibility of the +# variations at lower temperatures. + +chianti_temperature = np.log10(tpf.CHIANTI_temperature.to_value()) + +############################################################################## +# Differences overtime arise from an increase of the contamination layer on the +# CCD which blocks some of the X-rays thus reducing the effective area. +# For detailed information about the calculation of the XRT CCD contaminant layer thickness, +# you can refer to +# `Montana State University `__. +# +# Additional information is provided by +# `Narukage et. al. (2011) `__. + + +launch_datetime = "2006-09-22T23:59:59" + +launch_temperature_response = xrtpy.response.TemperatureResponseFundamental( + xrt_filter, launch_datetime, abundance_model="Photospheric" +).temperature_response() + +############################################################################## +# Now we can plot the temperature response versus the log of the CHIANTI temperature +# and compare the results for the launch date and the chosen date. + +plt.figure() + +plt.plot( + chianti_temperature, + np.log10(temperature_response.value), + label=f"{date_time}", +) +plt.plot( + chianti_temperature, + np.log10(launch_temperature_response.value), + label=f"{launch_datetime}", + color="red", +) + +plt.title("XRT Temperature Response") +plt.xlabel("Log(T) ($K$)") +plt.ylabel("$DN$ $cm^5$ $ s^-1$ $pix^-1$") +plt.legend() +plt.grid() + +plt.show() diff --git a/docs/getting_started.rst b/docs/getting_started.rst index 8f3cae3cf..d6af2ed45 100644 --- a/docs/getting_started.rst +++ b/docs/getting_started.rst @@ -7,7 +7,7 @@ on board the Hinode spacecraft. This page is intended for new users of XRTpy. Fo please refer to the `SolarSoft XRT Analysis Guide`_. XRTpy Objects -************* +------------- XRTpy currently provides access to the following core classes: @@ -24,41 +24,41 @@ It also includes functionality to: Visit our Examples page for step-by-step Jupyter notebook guides on how to use each feature. Channel -------- +^^^^^^^ The ``Channel`` class describes the configuration of a specific XRT filter channel. It includes details for the Charge-Coupled Device (CCD), Entrance Filter, Focal Plane Filters, Geometry, and Mirrors. Effective Area --------------- +^^^^^^^^^^^^^^ XRTpy calculates the effective area for each XRT filter channel, accounting for time-dependent contamination on the CCD. For more details, refer to the `SolarSoft XRT Analysis Guide`_. Temperature Response --------------------- +^^^^^^^^^^^^^^^^^^^^ XRTpy calculates the temperature response of XRT filter channels using the CHIANTI_ atomic database (version 10.0) and coronal abundances (:cite:t:`feldman:1992`). This produces a response function as a function of temperature, using an assumed emission model (:cite:t:`narukage:2011`, :cite:t:`narukage:2014`). Deriving Temperature and Emission Measure ------------------------------------------ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The ``temperature_from_filter_ratio`` function allows you to derive plasma temperature and emission measure from a pair of XRT images using the filter-ratio method. This mirrors the logic in the SolarSoft IDL routine of the same name. A usage example is available in the Examples section. Image Deconvolution with the PSF --------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The ``deconvolve`` function applies image deconvolution using the instrument's Point Spread Function (PSF) to sharpen XRT images. This is especially useful for recovering detail around bright or sharp solar structures. Light Leak Correction ---------------------- +^^^^^^^^^^^^^^^^^^^^^ The ``remove_lightleak`` function subtracts visible stray light from XRT synoptic composite images. This correction improves the quality of long-term coronal evolution studies. See our Examples section for how to use this function. Abundance Model Options ------------------------ +^^^^^^^^^^^^^^^^^^^^^^^ By default, XRTpy uses CHIANTI coronal abundances (:cite:t:`feldman:1992`). You may also choose: @@ -72,7 +72,9 @@ To use a different abundance model: from xrtpy.response import TemperatureResponseFundamental TemperatureResponseFundamental( - "Al-poly", "2022-07-04T23:43:12", abundance_model="hybrid" + "Al-poly", + "2022-07-04T23:43:12", + abundance_model="hybrid" ) You may also pass the ``abundance_model`` keyword to ``temperature_from_filter_ratio``. @@ -83,12 +85,12 @@ You may also pass the ``abundance_model`` keyword to ``temperature_from_filter_r Tools -***** +----- The ``xrtpy.response.tools`` module includes helpful utility functions to streamline workflows. It includes the following: Generate Temperature Responses ------------------------------- +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Use the ``generate_temperature_responses`` tool to compute the temperature response for one or more filters — including combinations like ``"Al-poly/Ti-poly"`` — with a single command. @@ -101,15 +103,15 @@ This function returns a list of ``TemperatureResponseFundamental`` objects, one from xrtpy.response.tools import generate_temperature_responses responses = generate_temperature_responses( - filters=["Al-poly", "Be-thick", "Al-poly/Ti-poly"], - obs_date="2020-07-04T00:00:00", - abundance="Hybrid", + filters=["Al-poly", "Be-thick", "Al-poly/Ti-poly"], + obs_date="2020-07-04T00:00:00", + abundance="Hybrid", ) for resp in responses: - print(f"Filter: {resp.filter_name}") - print(f" Temperatures: {resp.temperature[:3]}") - print(f" Response: {resp.response[:3]}") + print(f"Filter: {resp.filter_name}") + print(f" Temperatures: {resp.temperature[:3]}") + print(f" Response: {resp.response[:3]}") **Example Output:** @@ -117,16 +119,16 @@ This function returns a list of ``TemperatureResponseFundamental`` objects, one .. code-block:: text Filter: Al-poly - Temperatures: [100000. 112201.9 125892.516] K - Response: [8.34e-31 1.07e-30 1.53e-30] cm5 DN / (pix s) + Temperatures: [100000. 112201.9 125892.516] K + Response: [8.34e-31 1.07e-30 1.53e-30] cm5 DN / (pix s) Filter: Be-thick - Temperatures: [100000. 112201.9 125892.516] K - Response: [0.00e+00 1.73e-94 2.43e-84] cm5 DN / (pix s) + Temperatures: [100000. 112201.9 125892.516] K + Response: [0.00e+00 1.73e-94 2.43e-84] cm5 DN / (pix s) Filter: Al-poly/Ti-poly - Temperatures: [100000. 112201.9 125892.516] K - Response: [5.34e-34 7.24e-34 1.11e-33] cm5 DN / (pix s) + Temperatures: [100000. 112201.9 125892.516] K + Response: [5.34e-34 7.24e-34 1.11e-33] cm5 DN / (pix s) Each response object has the following attributes: @@ -138,7 +140,7 @@ This tool is useful on its own, but it also serves as a foundation for upcoming Data Products -************* +------------- XRT data products are available through the XRT website. These include: @@ -148,7 +150,7 @@ XRT data products are available through the XRT website. These include: For more information, visit the `XRT data products`_ page. Double Filter Combinations -************************** +-------------------------- XRTpy now supports double filter combinations such as ``"Al-poly/Ti-poly"`` in both the :obj:`EffectiveAreaFundamental ` and :obj:`TemperatureResponseFundamental ` classes. @@ -159,10 +161,10 @@ XRTpy now supports double filter combinations such as ``"Al-poly/Ti-poly"`` in b from xrtpy.response import EffectiveAreaFundamental, TemperatureResponseFundamental eff_area = EffectiveAreaFundamental( - "Al-poly/Ti-poly", "2020-08-17T09:00:00", abundance_model="photospheric" + "Al-poly/Ti-poly", "2020-08-17T09:00:00", abundance_model="photospheric" ) temp_resp = TemperatureResponseFundamental( - "C-poly/Ti-poly", "2025-07-10T12:00:00", abundance_model="coronal" + "C-poly/Ti-poly", "2025-07-10T12:00:00", abundance_model="coronal" ) The following combinations are currently supported: @@ -178,7 +180,7 @@ The following combinations are currently supported: X-Ray Filter Channels -********************* +--------------------- XRT uses two filter wheels to configure the imaging filter channel. Each wheel includes several filters and an open slot: @@ -204,12 +206,6 @@ Filter Wheel 2: Filter names in XRTpy are passed as strings like ``'Ti-poly'``. -References -********** - -Velasquez, J., Murphy, N., Reeves, K. K., Slavin, J., Weber, M., & Barnes, W. (2024). -*XRTpy: A Hinode-X-Ray Telescope Python Package*. JOSS, 9(100), 6396. -https://doi.org/10.21105/joss.06396 .. _CHIANTI: https://www.chiantidatabase.org/chianti_database_history.html .. _SolarSoft XRT Analysis Guide: https://xrt.cfa.harvard.edu/resources/documents/XAG/XAG.pdf diff --git a/docs/glossary.rst b/docs/glossary.rst index 178a371cf..15ffcdd8e 100644 --- a/docs/glossary.rst +++ b/docs/glossary.rst @@ -17,7 +17,13 @@ Glossary A numerical image processing technique used to correct for the blurring caused by the telescope's Point Spread Function (PSF), improving sharpness and visibility of fine structures. DEM - Differential Emission Measure (DEM) — a function that describes the distribution of plasma as a function of temperature along the line of sight. XRTpy will support DEM modeling in future versions. + Differential Emission Measure (DEM) — a function that describes how much plasma is present along the line of sight as a function of temperature. See :ref:`xrtpy-dem-overview` for a detailed overview of DEM theory, + usage, and the solver provided in XRTpy. + + DEM Inversion + The process of determining the temperature distribution of coronal plasma (the DEM) from a small number of filter intensities. Since more temperature bins are used than available filters, the problem is mathematically + underconstrained (“ill posed”), so regularization and smoothing are required to obtain a stable, physical solution. + DN Data Number (DN) — the digital value recorded by the CCD, representing the detected photon flux, usually in DN s\ :sup:`−1`\ . @@ -38,9 +44,20 @@ Glossary Contamination (related to the XRT) Refers to the gradual accumulation of material on the CCD and focal plane filters (FPFs), which reduces instrument throughput. This time-dependent degradation impacts effective area calculations and must be accounted for in data analysis. Refer to Section 2.5.3 *Contamination* in the `SolarSoft XRT Analysis Guide`_ for more information. + Monte Carlo DEM + A set of DEM solutions computed by adding random noise (based on intensity errors) to the observed intensities and re-solving the DEM multiple times. + The spread of these Monte Carlo solutions provides an estimate ofuncertainty in the DEM at each temperature. + PSF Point Spread Function — describes the response of the telescope to a point source of light. In XRTpy, it is used in deconvolution routines to sharpen images. + Response Matrix + A two-dimensional array containing the temperature response of each XRT filter interpolated onto the solver’s regular log10 temperature grid.This matrix connects + the DEM to the modeled filter intensities through the forward model: + + :math:`I_i^{model} = \sum_j DEM(T_j)\, R_i(T_j)\, T_j\, \Delta(\ln T)`. + + Solar Emission Spectra Emission spectra produced by solar plasma across a range of temperatures, calculated using spectral models such as CHIANTI. These spectra are used in temperature response and filter ratio methods diff --git a/docs/index.rst b/docs/index.rst index b2812ead1..bc56f7ed6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -21,6 +21,7 @@ Whether you're conducting research or just beginning to explore the world of X-r about_xrt install getting_started + dem_overview generated/gallery/index reference/index acknowledging_xrtpy @@ -30,3 +31,10 @@ Whether you're conducting research or just beginning to explore the world of X-r contributing code_of_conduct changelog/index + +Published Work +-------------- + +The following paper describes the XRTpy package and its initial release- v0.4.0: + +:cite:p:`velasquez:2024` diff --git a/pyproject.toml b/pyproject.toml index 0a44a9723..29ad1ce98 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,6 +45,7 @@ dependencies = [ "scikit-image>=0.21", "scipy>=1.11.1", "sunpy[map]>=5.1", + "lmfit>=1.2.2", ] optional-dependencies.dev = [ diff --git a/xrtpy/xrt_dem_iterative/__init__.py b/xrtpy/xrt_dem_iterative/__init__.py new file mode 100644 index 000000000..c45214265 --- /dev/null +++ b/xrtpy/xrt_dem_iterative/__init__.py @@ -0,0 +1,10 @@ +""" +XRT DEM Iterative Solver +""" + +from .dem_solver import XRTDEMIterative + +__all__ = [ + "XRTDEMIterative", + ] + diff --git a/xrtpy/xrt_dem_iterative/dem_plotting.py b/xrtpy/xrt_dem_iterative/dem_plotting.py new file mode 100644 index 000000000..a88085b6e --- /dev/null +++ b/xrtpy/xrt_dem_iterative/dem_plotting.py @@ -0,0 +1,145 @@ +__all__ = [ + "plot_dem", + "plot_dem_mc", +] + +import matplotlib.pyplot as plt +import numpy as np + + +def plot_dem(solver): + """ + Plot the base DEM solution in log10 psace (no Monte Carlo needed). + + Parameters + ---------- + solver : XRTDEMIterative + Fully initialized DEM solver. + """ + # Ensure base DEM is available + if not hasattr(solver, "dem"): + solver.solve() + + logT = solver.logT + dem = np.asarray(solver.dem, dtype=float) + + # Avoid log10(0) + dem_safe = np.clip(dem, 1e-40, None) + log_dem = np.log10(dem_safe) + + plt.figure(figsize=(8, 5)) + plt.step( + logT, + log_dem, + where="mid", + color="black", + linewidth=2.0, + label="Base DEM", + ) + + plt.xlabel(r"log$_{10} T$ [K]") + plt.ylabel(r"log$_{10}$ DEM [cm$^{-5}$ K$^{-1}$]") + plt.title("Base DEM Solution") + plt.grid(visible=True, alpha=0.3) + plt.xlim(logT.min(), logT.max()) + + # y-limits based on finite values + finite = np.isfinite(log_dem) + if np.any(finite): + y = np.log10(dem_safe[finite]) + pad = 0.25 * (y.max() - y.min() + 1e-6) + plt.ylim(y.min() - pad, y.max() + pad) + + plt.legend() + plt.tight_layout() + plt.show() + + +def plot_dem_mc( + solver, + mc_color="black", + base_color="#1E90FF", + alpha_mc=0.18, # 0.2 + lw_mc=1.0, + lw_base=2.0, + figsize=(9, 6), +): + """ + Plot DEM with Monte Carlo ensemble (if present). + + - Base DEM: thick colored step curve - Dodger Blue + - MC curves: thin transparent Black + - Automatically chooses limits, even if MC not present + + If no Monte Carlo results are present, this gracefully falls back + to plotting only the base DEM. + """ + + # Ensure base DEM is available + if not hasattr(solver, "dem"): + solver.solve() + + logT = solver.logT + base_dem = np.asarray(solver.dem, dtype=float) + base_safe_dem = np.clip(base_dem, 1e-100, None) # 1e-40 + log_base_dem = np.log10(base_safe_dem) + + # CHecking for Monte Carlo + has_mc = hasattr(solver, "mc_dem") and solver.mc_dem is not None + if has_mc: + mc_dem = np.asarray(solver.mc_dem, dtype=float) # shape (N+1, n_T) + # If someone set monte_carlo_runs=0, mc_dem may be (1, n_T) + N = max(0, mc_dem.shape[0] - 1) + else: + mc_dem = None + N = 0 + + plt.figure(figsize=figsize) + + # Plot Monte Carlo curves (index 1..N) + if has_mc and N > 0: + for i in range(1, N + 1): + dem_i = np.clip(mc_dem[i, :], 1e-100, None) # 1e-40 + plt.step( + logT, + np.log10(dem_i), + where="mid", + color=mc_color, + alpha=alpha_mc, + linewidth=lw_mc, + ) + + # Base DEM + plt.step( + logT, + log_base_dem, + where="mid", + color=base_color, + linewidth=lw_base, + label="Base DEM", + ) + + suffix = f" ({N} MC realizations)" if N > 0 else "" + plt.title("DEM with Monte Carlo" + suffix) + plt.xlabel(r"log$_{10} T$ [K]") + plt.ylabel(r"log$_{10}$ DEM [cm$^{-5}$ K$^{-1}$]") + plt.grid(visible=True, alpha=0.3) + plt.xlim(logT.min(), logT.max()) + + # y-limits based on base DEM and MC if available + logs = [np.log10(log_base_dem)] + + if has_mc and N > 0: + logs.append(np.log10(np.clip(mc_dem[1:, :], 1e-100, None)).ravel()) # 1e-40 + + logs_all = np.concatenate(logs) + + finite = np.isfinite(logs_all) + if np.any(finite): + y = logs_all[finite] + pad = 0.25 * (y.max() - y.min() + 1e-6) + plt.ylim(y.min() - pad, y.max() + pad) + + plt.legend() + plt.tight_layout() + plt.show() diff --git a/xrtpy/xrt_dem_iterative/dem_solver.py b/xrtpy/xrt_dem_iterative/dem_solver.py new file mode 100644 index 000000000..bd29a5fd3 --- /dev/null +++ b/xrtpy/xrt_dem_iterative/dem_solver.py @@ -0,0 +1,1233 @@ +__all__ = [ + "XRTDEMIterative", +] + +import warnings + +import astropy.units as u +import numpy as np +from lmfit import Parameters, minimize +from scipy.interpolate import interp1d + +from xrtpy.util.filters import validate_and_format_filters +from xrtpy.xrt_dem_iterative import dem_plotting + + +class XRTDEMIterative: + """ + Differential Emission Measure (DEM) solver for Hinode/XRT observations. + + This class implements a Python version of the IDL routine + `xrt_dem_iterative2.pro`, using spline-parameterized DEM curves and + iterative least-squares fitting. It supports Monte Carlo error analysis + and closely mirrors the logic of the original IDL algorithm. + + Parameters + ---------- + observed_channel : str or list of str, required + Names of the filters used in the observation (for example, + "Al-mesh", "Be-thin"). Must correspond one-to-one with the + temperature_responses argument. + observed_intensities : array-like, required + Observed intensities for each filter channel. Units are DN/s/pix. + temperature_responses : list, required + List of TemperatureResponseFundamental objects matching the filters. + Units are DN s^-1 pix^-1 cm^5. These can be created using + xrtpy.response.tools.generate_temperature_responses(). + intensity_errors : array-like, optional + Uncertainties in the observed intensities. If None, a default model + is used: max(0.03 * intensity, 2 DN/s/pix). + minimum_bound_temperature : float, optional + Minimum value of the log10(T) grid. Default is 5.5. + maximum_bound_temperature : float, optional + Maximum value of the log10(T) grid. Default is 8.0. + logarithmic_temperature_step_size : float, optional + Step size for the log10(T) grid. Default is 0.1. + monte_carlo_runs : int, optional + Number of Monte Carlo repetitions to perform. Default is 0 (disabled). + max_iterations : int, optional + Maximum number of function evaluations for lmfit. Default is 2000. + normalization_factor : float, optional + Internal scaling factor used during optimization. Default is 1e21. + + Notes + ----- + - All lists (observed_channel, observed_intensities, + temperature_responses) must be the same length. + - The log10(T) range must lie inside the native temperature grid + provided by all filter responses. + - If intensity_errors is not provided, a default model is used to + estimate uncertainties. + + SELFNOTEJOY + Add web-link to IDL script. + """ + + def __init__( + self, + observed_channel, + observed_intensities, + temperature_responses, + intensity_errors=None, + minimum_bound_temperature=5.5, + maximum_bound_temperature=8.0, + logarithmic_temperature_step_size=0.1, + monte_carlo_runs=0, + max_iterations=2000, + normalization_factor=1e21, + ): + """ + Notes + ----- + - All input lists (`observed_channel`, `observed_intensities`, and `temperature_responses`) + must be the same length. Each entry should correspond to one filter. + + - The temperature grid range (`minimum_bound_temperature`, `maximum_bound_temperature`) must lie entirely within the + response temperature ranges for **all** filters provided. + + - If `intensity_errors` is not provided, a model-based error estimate is used: + max(0.03 * observed_intensity, 2 (DN/s/pix)), as in the IDL original. + + - Default XRT filter names include: + {'Al-mesh', 'Al-poly', 'C-poly', 'Ti-poly', 'Be-thin', 'Be-med', 'Al-med', 'Al-thick', 'Be-thick', + 'Al-poly/Al-mesh', 'Al-poly/Ti-poly', 'Al-poly/Al-thick', 'Al-poly/Be-thick'} + """ + + # Validate and store filter names + if observed_channel is None or ( + hasattr(observed_channel, "__len__") and len(observed_channel) == 0 + ): + raise ValueError("`observed_channel` is required and cannot be empty.") + self.observed_channel = validate_and_format_filters(observed_channel) + + # Store intensity and error arrays + if observed_intensities is None or len(observed_intensities) == 0: + raise ValueError("`observed_intensities` is required and cannot be empty.") + self._observed_intensities = np.asarray(observed_intensities, dtype=float) + + if not np.all(np.isfinite(self._observed_intensities)): + raise ValueError("`observed_intensities` must be finite numbers.") + + # Errors + # NEW (defer to validate_inputs) + if intensity_errors is not None: + self._intensity_errors = np.asarray(intensity_errors, dtype=float) + else: + self._intensity_errors = None + + # Store temperature grid parameters + self._logarithmic_temperature_step_size = float( + logarithmic_temperature_step_size + ) + self._minimum_bound_temperature = float(minimum_bound_temperature) + self._maximum_bound_temperature = float(maximum_bound_temperature) + if not (self._minimum_bound_temperature < self._maximum_bound_temperature): + raise ValueError( + "minimum_bound_temperature must be < maximum_bound_temperature." + ) + + n_pts = ( + int( + np.floor( + (self._maximum_bound_temperature - self._minimum_bound_temperature) + / logarithmic_temperature_step_size + + 1e-9 + ) + ) + + 1 + ) + if n_pts < 4: + raise ValueError("Temperature grid must have at least 4 points.") + + # Validate Monte Carlo setting + if isinstance(monte_carlo_runs, bool): + raise TypeError( + "monte_carlo_runs must be a non-negative whole number, not a boolean." + ) + elif ( + isinstance(monte_carlo_runs, int | np.integer) + or isinstance(monte_carlo_runs, float) + and monte_carlo_runs.is_integer() + ): + self._monte_carlo_runs = int(monte_carlo_runs) + else: + raise ValueError( + "monte_carlo_runs must be a non-negative whole number (e.g., 0, 1, 100). " + "Decimal values are not allowed." + ) + + if self._monte_carlo_runs < 0: + raise ValueError("monte_carlo_runs must be ≥ 0.") + + # Validate max_iterations + if not isinstance(max_iterations, int | np.integer) or max_iterations <= 0: + raise ValueError("max_iterations must be a positive integer.") + + self._max_iterations = int(max_iterations) + + # Check logarithmic_temperature_step_size is positive + if self._logarithmic_temperature_step_size <= 0: + raise ValueError( + "logarithmic_temperature_step_size must be a positive scalar." + ) + + # Store temperature response objects + self.responses = temperature_responses + + if temperature_responses is None or len(temperature_responses) == 0: + raise ValueError("`temperature_responses` is required and cannot be empty.") + + # Validate that the temperature grid falls within the responses + for r in self.responses: + logT_grid = np.log10(r.temperature.to_value(u.K)) + if not ( + self._minimum_bound_temperature >= logT_grid.min() + and self._maximum_bound_temperature <= logT_grid.max() + ): + raise ValueError( + f"The specified temperature range [{minimum_bound_temperature}, {maximum_bound_temperature}] is outside the bounds of one or more filter response grids.\n" + "Please ensure the temperature range fits within all responses.\n" + "Hint: Default response range is logT = 5.5 to 8.0. You can view each response's logT range via: [r.temperature for r in responses]" + ) + + # Check consistency between inputs + if not ( + len(self._observed_intensities) + == len(self.responses) + == len(self.observed_channel) + ): + raise ValueError( + f"\nLength mismatch in inputs:\n" + f" Observed intensities: {len(self._observed_intensities)}\n" + f" Responses: {len(self.responses)}\n" + f" Filter channels: {len(self.observed_channel)}\n" + ) + + try: + value = float(normalization_factor) + except (TypeError, ValueError) as err: + raise ValueError(f"Invalid normalization_factor: {err}") from err + + if value <= 0: + raise ValueError("normalization_factor must be a positive number.") + + self._normalization_factor = value + # self._normalization_factor = value + + self._using_estimated_errors = ( + False # track whether default error model has been used + ) + + def validate_inputs(self) -> None: + """ + Validate user-provided inputs. Raises ValueError on any issue. + + Intentionally separate from __init__ so tests and users can construct + the object first, then explicitly validate (matches test expectations). + """ + # 1) observed_channel non-empty + if self.observed_channel is None or len(self.observed_channel) == 0: + raise ValueError("`observed_channel` is required and cannot be empty.") + + # 2) intensities: length & finite + if self._observed_intensities is None or len(self._observed_intensities) == 0: + raise ValueError("`observed_intensities` is required and cannot be empty.") + if not np.all(np.isfinite(self._observed_intensities)): + raise ValueError("`observed_intensities` must be finite numbers.") + + # 3) responses present + if self.responses is None or len(self.responses) == 0: + raise ValueError("`temperature_responses` is required and cannot be empty.") + + # 4) lengths consistent between filters/intensities/responses + if not ( + len(self._observed_intensities) + == len(self.responses) + == len(self.observed_channel) + ): + raise ValueError( + "Length mismatch: intensities, responses, and observed_channel must match." + ) + + # 5) temperature grid sanity + if not (self._minimum_bound_temperature < self._maximum_bound_temperature): + raise ValueError( + "minimum_bound_temperature must be < maximum_bound_temperature." + ) + if self._logarithmic_temperature_step_size <= 0: + raise ValueError( + "logarithmic_temperature_step_size must be a positive scalar." + ) + n_pts = ( + int( + np.floor( + (self._maximum_bound_temperature - self._minimum_bound_temperature) + / self._logarithmic_temperature_step_size + + 1e-9 + ) + ) + + 1 + ) + if n_pts < 4: + raise ValueError("Temperature grid must have at least 4 points.") + + # 6) grid range inside every response + for r in self.responses: + # logT_grid = np√.log10(r.temperature.value) + logT_grid = np.log10(r.temperature.to_value(u.K)) + if not ( + self._minimum_bound_temperature >= logT_grid.min() + and self._maximum_bound_temperature <= logT_grid.max() + ): + raise ValueError( + f"The specified temperature range [{self._minimum_bound_temperature}, {self._maximum_bound_temperature}] " + "is outside the bounds of one or more filter response grids." + ) + + # 7) intensity_errors length & finiteness (only if provided) + if self._intensity_errors is not None: + if self._intensity_errors.shape != self._observed_intensities.shape: + raise ValueError( + "Length of intensity_errors must match observed_intensities." + ) + if not np.all(np.isfinite(self._intensity_errors)) or np.any( + self._intensity_errors < 0 + ): + raise ValueError("`intensity_errors` must be finite and >= 0.") + + # 8 warning + if np.all(self._observed_intensities == 0): + warnings.warn( + "\n\n All observed intensities are zero. DEM solution will yield zero. " + "Object created, but solving will return DEM=0. \n\n", + stacklevel=2, + ) + + # success -> no return value + return None + + def __repr__(self): + return ( + f"" + ) + + # @property #Removed if not used + # def name(self) -> str: + # """ + # The XRT filter channel name, standardized (e.g. "Al-mesh"). + # """ + # return self._name + + @property + def observed_intensities( + self, + ) -> ( + u.Quantity + ): # Add method to account for known values not worth observed_intensities + """ + Observed intensities with physical units. + + Returns + ------- + `~astropy.units.Quantity` + Intensities in DN/s/pix for each filter channel. + Where "pix" means a one-arcsecond, full -resolution XRT pixel. + """ + return self._observed_intensities * (u.DN / u.s) + + @property + def filter_names(self): + """ + List of filter names from the temperature responses. + """ + return [r.filter_name for r in self.responses] + + @property + def response_temperatures(self): + """ + List of temperature grids (K) for each filter response. + """ + return [r.temperature for r in self.responses] + + @property + def response_values(self): + """ + List of response values (DN cm^5 / pix / s) for each filter. + """ + return [r.response for r in self.responses] + + @property + def minimum_bound_temperature(self): + """ + Lower bound of log10 temperature grid. + """ + return self._minimum_bound_temperature + + @property + def maximum_bound_temperature(self): + """ + Upper bound of log10 temperature grid. + """ + return self._maximum_bound_temperature + + @property + def logarithmic_temperature_step_size(self): + """ + Bin width of log10 temperature grid. + """ + return self._logarithmic_temperature_step_size + + @property + def min_observational_error(self): + """ + Default - Minimum absolute observational intensity error applied to DN/s/pix when intensity error is not provided. + """ + return 2 * (u.DN / u.s) + + @property + def relative_error(self): + """ + Relative error used to scale intensity if an error is not provided. + Default is 0.03 (3%). + """ + return 0.03 + + @property + def intensity_errors(self) -> u.Quantity: + """ + Return the intensity uncertainty values. + + If the user supplied intensity_errors, those values are returned. + Otherwise a default model is used: + + sigma = max(0.03 * intensity, 2 DN/s/pix) + + This behavior mirrors the default uncertainty logic of the IDL routine + xrt_dem_iterative2.pro. + + `~astropy.units.Quantity` + Intensity errors in DN/s for each filter. + + For details, see: + https://hesperia.gsfc.nasa.gov/ssw/hinode/xrt/idl/util/xrt_dem_iterative2.pro + """ + if self._intensity_errors is not None: + return self._intensity_errors * (u.DN / u.s) + + if not self._using_estimated_errors: + warnings.warn( + ( + "\n\nNo intensity_errors provided. Using default model: " + "max(relative-error * observed_intensity, min_observational_error)\n" + f"=> relative_error = {self.relative_error}, " + f"min_observational_error = {self.min_observational_error.value} DN/s\n" + "See: https://hesperia.gsfc.nasa.gov/ssw/hinode/xrt/idl/util/xrt_dem_iterative2.pro\n\n" + ), + category=UserWarning, + stacklevel=2, + ) + + self._using_estimated_errors = True + + # Fixed in units + estimated = np.maximum( + (self.relative_error * self._observed_intensities) * (u.DN / u.s), + self.min_observational_error, + ) + return estimated + + @property + def monte_carlo_runs(self) -> int: + """ + Return + ------ + int + Number of Monte Carlo runs to perform (0 = disabled). + """ + return self._monte_carlo_runs + + @property + def normalization_factor(self): + """ + Scaling factor used during DEM optimization to stabilize the spline fit. + Corresponds to `normalization_factor` in IDL (default 1e21). + """ + return self._normalization_factor + + @property + def max_iterations(self): + """ + Maximum number of iterations allowed in the least-squares DEM solver + (e.g., when using `lmfit.minimize`). Default is 2000. + """ + return self._max_iterations + + def create_logT_grid(self): + """ + Construct the regular log10 temperature grid for DEM calculations. + + This builds a regularly spaced grid in log10(temperature), then converts it + to linear temperature for use in the DEM integral. + + Notes + ----- + - IDL's `xrt_dem_iterative2.pro` describes this as the "regular logT grid". + - Two forms of the temperature grid are stored: + * self.logT : log10(T) values (dimensionless) + * self.T : linear temperatures (Kelvin, astropy.units.Quantity) + - The grid is inclusive of both `minimum_bound_temperature` and `maximum_bound_temperature`, with step size `logarithmic_temperature_step_size`. + + Additional attributes created: + - self.dlogT : float + Step size in log10(T) (dimensionless). + - self.dlnT : float + Step size in natural log(T). Useful for IDL-style integrals of the form: + F = int. DEM(T) * R(T) * T d(ln T) + """ + # number of bins including endpoints - if default values are used - end value is 26 + self.n_bins = ( + int( + round( + (self._maximum_bound_temperature - self._minimum_bound_temperature) + / self._logarithmic_temperature_step_size + ) + ) + + 1 + ) + + # This matches the IDL temperature grid exactly. self.logT & self.T. + # inclusive logT grid (IDL-style regular grid) + # Units = 'log K. Runs from minimum_bound_temperature to Max_T with bin-width = DT + # SELFNOTEJOY- Do we need to add units - current holds no units- it's wokring correctly - Should this on the Test as well?- I don't think it- it's noted in IDL but used with units + + # np.linspace over np.arange - simple reproduces that reliably:endpoint included, exact number of bins, and no accumulating floating-point drift - Best match to IDL + self.logT = np.linspace( + self._minimum_bound_temperature, + self._maximum_bound_temperature, + self.n_bins, + ) + + # linear temperature grid in Kelvin + self.T = (10.0**self.logT) * u.K + + self.dlogT = float( + self._logarithmic_temperature_step_size + ) # scalar spacing (dimensionless and natural-log equivalent) + + self.dlnT = ( + np.log(10.0) * self.dlogT + ) # for IDL-style intergral DEM(T) * R(T) * T dlnT - IDL “regular logT grid” + + def _interpolate_responses_to_grid(self): + """ + Interpolate all filter responses onto the solver's regular log10(T) grid. + + This constructs the response matrix used in the DEM forward model. + Each filter's native temperature response (R(T)) is interpolated to the + grid defined by self.logT. Extrapolated values outside the native + response range are set to zero. + + Equivalent to the "Res_Mat" construction in the IDL routine xrt_dem_iterative2.pro. + + Notes + ----- + - Response units (from XRTpy) are DN s^-1 pix^-1 cm^5. + - Output matrix has shape (n_filters, n_temperatures). + - Rows correspond to filters; columns correspond to temperature bins. + + ------- + IDL method of Interpolate emissivity. + Interpolate all filter responses onto the common logT grid and build + the response matrix. + + Equivalent to constructing `Res_Mat` in IDL's `xrt_dem_iterative2.pro` + and in the DEM_Solver PDF documentation. + + Attributes Created + ------------------ + interpolated_responses : list of ndarray + _response_matrix : ndarray + Stacked filter responses on the uniform logT grid. + """ + if not hasattr(self, "logT"): + raise AttributeError( + "Temperature grid missing. Call create_logT_grid() first." + ) + + rows = [] + for T_orig, R_orig in zip( + self.response_temperatures, self.response_values, strict=False + ): # Make sure that R_orig.value is indeed in DN/s/pix per cm^5 + logT_orig = np.log10(T_orig.to_value(u.K)) + + response_vals = R_orig.value # already in correct physical units for XRTpy #NOTEFORJOY- TRIPLE check this + + interp_func = interp1d( + logT_orig, + response_vals, + kind="linear", + bounds_error=False, + fill_value=0.0, + assume_sorted=True, + ) + rows.append(interp_func(self.logT)) + + self.interpolated_responses = rows + self._response_matrix = np.vstack(rows).astype(float) + + # Store the physical unit for clarity + # self._response_unit = u.DN / u.s / u.pix / (u.cm**5) + self._response_unit = (u.DN / u.s / u.pix) * u.cm**5 + + # Quick sanity check + if self._response_matrix.shape != (len(self.responses), self.logT.size): + raise RuntimeError("Interpolated response matrix has unexpected shape.") + + @property + def response_matrix(self): + """ + Response matrix (n_filters x n_temperatures) after interpolation. + + Units: DN s^-1 pix^-1 cm⁵ per emission measure. + + Equivalent to `Res_Mat` in IDL's `xrt_dem_iterative2.pro`. + + Raises + ------ + AttributeError + If `_interpolate_responses_to_grid()` has not been called yet. + """ + if not hasattr(self, "_response_matrix"): + raise AttributeError( + "Response matrix not available. Call _interpolate_responses_to_grid() first." + ) + return self._response_matrix + + def _prepare_scaled_observations(self): + """ + Prepare the scaled observed intensities and uncertainties + exactly as done in the IDL routine xrt_dem_iterative2.pro. + + IDL equivalent: + input1.i_obs = input1.i_obs / solv_factor + input1.i_err = input1.i_err / solv_factor + """ + # Extract values as plain floats (DN/s/pix) + intensities_scaled_raw = ( + self.observed_intensities.value + ) # Might just remove this line and up in the normalization + sigma_intensity_errors_raw = self.intensity_errors.to_value( + u.DN / u.s + ) # Might just remove this line and up in the normalization + + # Apply normalization + self.intensities_scaled = intensities_scaled_raw / self.normalization_factor + self.sigma_scaled_intensity_errors = ( + sigma_intensity_errors_raw / self.normalization_factor + ) + + # Store for solver + self._scaled_prepared = True + + # ====================================================================== + # DEM INITIALIZATION AND SOLVER METHODS + # ====================================================================== + + def _estimate_initial_dem(self, cutoff: float = 1.0 / np.e) -> np.ndarray: + """ + Compute an initial DEM estimate closely following the structure of the IDL routine xrt_dem_iter_estim. + + The IDL code performs a rough DEM estimate by evaluating intensities + relative to response peaks, but xrt_dem_iterative2 ultimately replaces + that estimate with a flat log10(DEM) curve before calling the solver. + + This method repeats the peak-finding logic for diagnostic purposes, but + the final DEM passed into the solver is always: + + log10(DEM(T)) = 0.0 for all temperature bins + + which corresponds to DEM(T) = 1 in arbitrary units. This reproduces the + IDL initial condition exactly. + + Parameters + ---------- + cutoff : float, optional + Fraction of the peak response used to define the usable window + around a channel's emissivity peak. Default is 1/e (approximately + 0.3679). + + Returns + ------- + ndarray + Initial log10(DEM) estimate on self.logT. This is always a flat + array of zeros (IDL-equivalent behavior). + """ + if not hasattr(self, "logT"): + raise AttributeError( + "Temperature grid missing. Call create_logT_grid() first." + ) + if not hasattr(self, "_response_matrix"): + raise AttributeError( + "Response matrix missing. Call _interpolate_responses_to_grid() first." + ) + + # Optional: store the peak-based estimates for diagnostics only. + # These are NOT used to set the initial DEM (IDL overwrites them + # with a flat DEM before handing off to the solver). + t_peaks = [] + log_dem_estimates = [] + + # Loop over each filter/channel with non-zero intensity + for T_orig, R_orig, I_obs in zip( + self.response_temperatures, + self.response_values, + self._observed_intensities, + strict=False, + ): + logT_orig = np.log10(T_orig.to_value(u.K)) + # Make sure the response is in DN s^-1 pix^-1 per EM (cm^-5) + R_vals = R_orig.to_value((u.DN / u.s / u.pix) * u.cm**5) + + if I_obs <= 0 or np.all(R_vals <= 0): + continue # skip unusable channel + + # 1. Peak location (logT) + max_idx = np.argmax(R_vals) + peak_val = R_vals[max_idx] + t_peak_raw = logT_orig[max_idx] + + # Round to nearest grid step in logT, similar to round_off(..., 0.1) + step = self._logarithmic_temperature_step_size + t_peak = np.round(t_peak_raw / step) * step + + # 2. Good window (where R > cutoff * peak) + good = np.where(R_vals > peak_val * cutoff)[0] + if good.size < 1: + continue + + # 3. Compute denominator integral: sum(T * R * dlnT) + T_good = 10.0 ** logT_orig[good] # [K] + R_good = R_vals[good] + # Native spacing in log10(T) + if logT_orig.size > 1: + dlogT_native = np.mean(np.diff(logT_orig)) + else: + # Degenerate case; fall back to solver grid spacing + dlogT_native = step + dlnT_native = np.log(10.0) * dlogT_native + denom = np.sum(T_good * R_good * dlnT_native) + + if denom <= 0: + continue + + # 4. DEM estimate at the peak (for diagnostics only) + dem_peak = I_obs / denom # [cm^-5 K^-1] + if dem_peak <= 0 or not np.isfinite(dem_peak): + continue + + log_dem_est = np.log10(dem_peak) + t_peaks.append(t_peak) + log_dem_estimates.append(log_dem_est) + + # Compact duplicate peak temperatures by averaging (diagnostic only) + if t_peaks: + uniq = {} + for t, val in zip(t_peaks, log_dem_estimates, strict=False): + uniq.setdefault(t, []).append(val) + t_peaks_uniq = np.array(sorted(uniq.keys())) + log_dem_uniq = np.array([np.mean(uniq[t]) for t in t_peaks_uniq]) + # Store raw estimated peaks for debugging/inspection if desired + self._raw_estimated_dem_peaks = (t_peaks_uniq, log_dem_uniq) + else: + self._raw_estimated_dem_peaks = (np.array([]), np.array([])) + + # IDL BEHAVIOR: override with flat initial DEM + # xrt_dem_iter_estim ultimately does: dem = 0.0*findgen(nt) + 1.0 ; Use flat dem for initial guess on a regular logT grid. We mirror that here exactly: + + # est_log_dem_on_grid = np.ones_like(self.logT, dtype=float) * 1.0 NOV20 + # est_log_dem_on_grid = np.ones_like(self.logT, dtype=float) * 0.0 #NOTEFORJOY + est_log_dem_on_grid = np.zeros_like(self.logT) + + # Return the intial first guessed DEM + + # Store for later use by the solver + self._initial_log_dem = est_log_dem_on_grid + + return est_log_dem_on_grid + + def _prepare_spline_system(self): + """ + Prepare the spline-based DEM parameterization. + + This mirrors the IDL routine mp_prep and sets up all arrays needed by + the least-squares solver, including: + + - self.n_spl : number of spline knots + - self.spline_logT : knot positions (evenly spaced in log10(T)) + - self.spline_log_dem : initial values of log10(DEM) at each knot + - self.pm_matrix : response matrix multiplied by T * d(ln T) + - self.weights : all ones (IDL uses a channel weighting mask) + - self.abundances : all ones + + pm_matrix corresponds to: + pm[i, j] = R_i(T_j) * T_j * d(ln T) + + which appears in the forward model: + I_model_i = sum_j DEM(T_j) * pm[i, j] + """ + + # Number of channels + n_line = len(self._observed_intensities) + + # IDL: n_spl = min(n_line - 1, 7) - Make this a keyword in the class so use can tune it? + self.n_spl = min(max(n_line - 1, 1), 7) + + # Weights and abundances (IDL sets all =1) + # Later, should I a use_line mask (IDL ignores lines with i_obs=0), but you can add that when you need it. + self.weights = np.ones(n_line, dtype=float) + self.abundances = np.ones(n_line, dtype=float) + + # pm_matrix = R(T) * T * dlnT (IDL line: emis * 10^t * alog(10^dt)) + # units - DN/s/pix/cm^5 * K * dLnT * DEM == DN/s/PIX + T_linear = self.T.to_value(u.K) + + self.pm_matrix = (self._response_matrix * T_linear * self.dlnT).astype(float) + + # Knot positions are evenly spaced in logT (IDL spl_t) + self.spline_logT = np.linspace(self.logT.min(), self.logT.max(), self.n_spl) + + # Initial spline DEM values: sample from initial logDEM grid + # (IDL spline(est_t, est_dem, spl_t)) + interp_init = interp1d( + self.logT, + self._initial_log_dem, # IDL is flat logDEM = 1.0 + kind="linear", # IDL uses a cubic spline later NOTEFORJOY NOV20 + bounds_error=False, + fill_value="extrapolate", + ) + + self.spline_log_dem = interp_init(self.spline_logT) + + def _build_lmfit_parameters(self): + """ + Build lmfit.Parameters object representing log10(DEM) at the spline knots. + IDL limits each spline DEM parameter to [-20, 0]. + """ + + if not hasattr(self, "spline_log_dem"): + raise RuntimeError("Run _prepare_spline_system() first.") + + params = Parameters() + + for i in range(self.n_spl): + params.add( + f"knot_{i}", + value=float(self.spline_log_dem[i]), + min=-20.0, + max=0.0, + vary=True, + ) + + return params + + def _reconstruct_dem_from_knots(self, params): + """ + Construct DEM(T) on self.logT using spline of log10(DEM) at knot positions. + Uses a natural cubic spline interpolation in log10(DEM) space. + """ + from scipy.interpolate import CubicSpline + + knot_vals = np.array([params[f"knot_{i}"].value for i in range(self.n_spl)]) + + # Or used the code above but switch from linear to kind="cubic" + cs = CubicSpline(self.spline_logT, knot_vals, bc_type="natural") + log_dem = cs(self.logT) + dem = 10.0**log_dem + return dem + + def _residuals(self, params): + """ + Compute residuals for use by the least-squares optimizer. + + Residuals are computed as: + residual_i = (I_model_i - I_observed_i) / sigma_i + + where: + I_model = (pm_matrix @ DEM) * abundances + + Returns + ------- + ndarray + Residuals for each filter channel. + """ + + # 1. DEM(T) + dem = self._reconstruct_dem_from_knots(params) + + # 2. Modeled intensities (IDL: i_mod = (dem ## pm) * abunds) + i_mod = (self.pm_matrix @ dem) * self.abundances + + # 3. Observed (scaled) + y_scaled = self.intensities_scaled # i_obs / solv_factor + sigma_scaled = self.sigma_scaled_intensity_errors + + # 4. Residuals = (i_mod - y_obs) * weights / sigma + residuals = (i_mod - y_scaled) * self.weights / sigma_scaled + + # chi^2 history, mostly for debugging + chi2_val = np.sum(residuals**2) + if not hasattr(self, "_iteration_chi2"): + self._iteration_chi2 = [] + self._iteration_chi2.append(chi2_val) + + return residuals + + def _solve_single_dem(self, observed_intensities_vals: np.ndarray): + """ + This method solves the DEM for one set of intensities only, without + Monte Carlo perturbation. + """ + + nf = self._normalization_factor + + # 1. scaled obs/errors + self.intensities_scaled = observed_intensities_vals / nf + sigma_phys = self.intensity_errors.to_value(u.DN / u.s) + self.sigma_scaled_intensity_errors = sigma_phys / nf + + # 2. trivial nosolve case + if np.all(self.intensities_scaled == 0.0): + dem_model = np.zeros_like(self.logT) + dem_phys = dem_model * nf + modeled_intensities_phys = np.zeros_like(observed_intensities_vals) + return dem_phys, modeled_intensities_phys, 0.0, None + + # 3. initial guess (log10 DEM_model on grid) + init_log_dem = self._estimate_initial_dem() # flat ~ 1.0 in IDL + self._initial_log_dem = init_log_dem + + # 4. spline system using that initial guess + self._prepare_spline_system() + params0 = self._build_lmfit_parameters() # values = initial_log_dem at knots + + # 5. run minimizer + result = minimize( + self._residuals, params0, max_nfev=self._max_iterations + ) # method='leastsq' + + # THIS is the critical part – use *result.params*, not params0 <<< + dem_model = self._reconstruct_dem_from_knots(result.params) # DEM_model(T) + dem_phys = dem_model * nf + + i_mod_scaled = (self.pm_matrix @ dem_model) * self.abundances + modeled_intensities_phys = i_mod_scaled * nf + + resid = self._residuals(result.params) + chisq = float(np.sum(resid**2)) + + return dem_phys, modeled_intensities_phys, chisq, result + + def _run_monte_carlo(self): + """ + Replicates IDL's Monte Carlo loop. + Produces: + - self.mc_dem shape (n_T, N+1) + - self.mc_base_obs shape (n_obs, N+1) + - self.mc_mod_obs shape (n_obs, N+1) + - self.mc_chisq shape (N+1,) + + Note that N+1 rows means: row 0 = base case, rows 1..N = MC. + """ + + n_obs = len(self._observed_intensities) + nT = len(self.logT) + N = self._monte_carlo_runs + + # Prepare arrays + mc_dem = np.zeros((nT, N + 1)) + mc_base = np.zeros((n_obs, N + 1)) + mc_mod = np.zeros((n_obs, N + 1)) + mc_chi = np.zeros(N + 1) + + # Base run first (IDL puts real data in column 0) + dem = self.dem # already scaled back by normalization + mc_dem[:, 0] = dem + mc_base[:, 0] = self._observed_intensities # unscaled + mc_mod[:, 0] = self.modeled_intensities # unscaled + mc_chi[0] = self.current_chi2 + + # --- Run the MC loops --- + rng = np.random.default_rng() # like systime(1) + + for ii in range(1, N + 1): + # Step 1: Perturbed intensities (scaled) + perturbed = ( + self.intensities_scaled + + rng.normal(size=n_obs) * self.sigma_scaled_intensity_errors + ) + perturbed = np.clip(perturbed, 0, None) + + # Store unscaled in mc_base + mc_base[:, ii] = perturbed * self.normalization_factor + + # If all zero → nosolve=True + if np.all(perturbed == 0): + mc_dem[:, ii] = 0.0 + mc_mod[:, ii] = 0.0 + mc_chi[ii] = 0.0 + continue + + # Step 2: assign perturbed intensities + self.intensities_scaled = perturbed + + # Step 3: Rebuild spline system + self._prepare_spline_system() + + # Step 4: Solve via lmfit + params = self._build_lmfit_parameters() + out = minimize(self._residuals, params, max_nfev=self.max_iterations) + + # Step 5: Reconstruct DEM + dem = self._reconstruct_dem_from_knots(out.params) + dem_scaled = dem * self.normalization_factor # unscale + mc_dem[:, ii] = dem_scaled + + # Step 6: Compute modeled intensities + modeled = (self.pm_matrix @ dem) * self.abundances + mc_mod[:, ii] = modeled * self.normalization_factor + + # Step 7: Compute chi-square + resid = self._residuals(out.params) + mc_chi[ii] = np.sum(resid**2) + + # store results + self.mc_dem = mc_dem + self.mc_base_obs = mc_base + self.mc_mod_obs = mc_mod + self.mc_chisq = mc_chi + + def solve(self): + """ + High-level DEM solver. + + Python analogue of IDL's xrt_dem_iterative2.pro: + + 1. Validate inputs. + 2. Build the logT grid and interpolate temperature responses. + 3. Solve ONE base DEM using the original (unperturbed) intensities. + 4. If Monte Carlo is requested (monte_carlo_runs > 0), perform N + perturbed solves by adding Gaussian noise to the base intensities. + 5. Store all outputs on the instance for later analysis/plotting. + + After calling solve(), the following attributes are defined: + + Base solution + ------------- + logT : ndarray (n_T,) + log10 temperature grid [K]. + dem : ndarray (n_T,) + Best-fit DEM(T) in physical units [cm^-5 K^-1]. + chisq : float + Chi-square of the base fit (sum of squared residuals). + modeled_intensities : ndarray (n_channels,) + Best-fit modeled intensities in [DN s^-1 pix^-1]. + _base_fit_result : lmfit.MinimizerResult + Full lmfit result object for diagnostics. + + Monte Carlo products + -------------------- + mc_dem : ndarray (N+1, n_T) + DEM curves for base (row 0) and each Monte Carlo run (rows 1..N), + in physical units [cm^-5 K^-1]. + mc_chisq : ndarray (N+1,) + Chi-square values for base (index 0) and each MC run. + mc_base_obs : ndarray (N+1, n_channels) + Observed intensities [DN s^-1 pix^-1] for base + each MC run. + Row 0 = original observation; rows 1..N = perturbed. + mc_mod_obs : ndarray (N+1, n_channels) + Modeled intensities [DN s^-1 pix^-1] corresponding to mc_dem. + """ + + # Validate inputs (IDL: argument checks near top) + self.validate_inputs() + + # 1) Build logT grid and response matrix - IDL: regular logT grid + interpolated emissivities + self.create_logT_grid() + self._interpolate_responses_to_grid() + + # Base observed intensities in physical units [DN/s/pix] + base_obs_phys = np.asarray(self._observed_intensities, dtype=float) + + # 2) Solve BASE DEM (unperturbed intensities) Corresponds to the first call to xrt_dem_iter_nowidget in IDL. + + dem_base, mod_base, chisq_base, base_result = self._solve_single_dem( + observed_intensities_vals=base_obs_phys + ) + + # Store base solution + self.logT_solution = self.logT.copy() # alias + self.dem = dem_base # [cm^-5 K^-1] + self.chisq = chisq_base # chi-square + self.modeled_intensities = mod_base # [DN/s/pix] + self._base_fit_result = base_result + + # 3) Allocate Monte Carlo arrays (IDL: base_obs, dem_out, chisq, mod_obs) + + n_T = self.logT.size + n_ch = base_obs_phys.size + N = self.monte_carlo_runs + + self.mc_dem = np.zeros((N + 1, n_T), dtype=float) + self.mc_chisq = np.zeros((N + 1,), dtype=float) + self.mc_base_obs = np.zeros((N + 1, n_ch), dtype=float) + self.mc_mod_obs = np.zeros((N + 1, n_ch), dtype=float) + + # Row 0 = base solution (unperturbed) + self.mc_dem[0, :] = dem_base + self.mc_chisq[0] = chisq_base + self.mc_base_obs[0, :] = base_obs_phys + self.mc_mod_obs[0, :] = mod_base + + # 4) Monte Carlo loop + if N > 0: + rng = np.random.default_rng() # like IDL's systime(1) seeding + + # Intensity errors in physical units [DN/s/pix] + sigma_phys = self.intensity_errors.to_value(u.DN / u.s) + + for ii in range(1, N + 1): + # Lightweight progress indicator + if ii % max(1, N // 20) == 0: + print(f" - Monte Carlo run {ii}/{N}") + + # 4a) Perturb intensities: I' = I + N(0, sigma), clipped at 0 + noise = rng.normal(loc=0.0, scale=sigma_phys, size=base_obs_phys.shape) + obs_pert = base_obs_phys + noise + obs_pert = np.maximum(obs_pert, 0.0) # IDL: >0 to avoid negatives + + # 4b) Solve DEM for this perturbed realization + dem_i, mod_i, chisq_i, _ = self._solve_single_dem( + observed_intensities_vals=obs_pert + ) + + # 4c) Store Monte Carlo results + self.mc_dem[ii, :] = dem_i + self.mc_chisq[ii] = chisq_i + self.mc_base_obs[ii, :] = obs_pert + self.mc_mod_obs[ii, :] = mod_i + + return self.dem + + def summary(self): + """ + Print a detailed, diagnostic summary of the DEM solver state. + + This provides: + - Input observation details + - Temperature grid configuration + - Response matrix status + - Spline system configuration + - Base DEM fit results + - Monte Carlo statistics (if available) + - Available plotting helpers + """ + + print("\n" + "=" * 76) + print(" XRTpy DEM Iterative — Solver Summary") + print("=" * 76) + + print("\nINPUT DATA") + print("-" * 70) + print(f" Filters: {self.filter_names}") + print( + f" Observed Intensities: {np.array(self._observed_intensities)} DN/s/pix" + ) + print(f" Number of channels: {len(self._observed_intensities)}") + + # Error model + if self._intensity_errors is not None: + print(" Intensity Errors: User-provided") + else: + print(" Intensity Errors: Auto-estimated (3% of I, min=2 DN/s)") + + print(f" Error values (DN/s): {self.intensity_errors.to_value('DN/s')}\n") + + print("\nTEMPERATURE GRID") + print("-" * 70) + if hasattr(self, "logT"): + print(f" logT range: {self.logT[0]:.2f} to {self.logT[-1]:.2f}") + print(f" Number of temperature bins: {len(self.logT)}") + print(f" logT (grid spacing): {self.dlogT:.3f}") + print(f" lnT (natural log spacing): {self.dlnT:.3f}") + else: + print(" Grid has not been constructed (call solve()).") + + print("\nRESPONSE MATRIX") + print("-" * 70) + if hasattr(self, "_response_matrix"): + print( + f" Matrix shape: {self._response_matrix.shape} (filters x T bins)" + ) + print(f" Response units: {self._response_unit}") + else: + print(" Response matrix not constructed.") + + print("\nSOLVER CONFIGURATION") + print("-" * 70) + print(f" Normalization factor: {self.normalization_factor:.2e}") + print(f" Max iterations: {self.max_iterations}") + print(f" Monte Carlo runs: {self.monte_carlo_runs}") + if hasattr(self, "n_spl"): + print(f" Number of spline knots: {self.n_spl}") + print(f" Knot positions (logT): {getattr(self, 'spline_logT', 'N/A')}") + else: + print(" Spline system not prepared yet.") + + print("\nINITIAL DEM GUESS") + print("-" * 70) + if hasattr(self, "_initial_log_dem"): + print(" Initial DEM assumption: flat log10(DEM) (IDL-style)") + print(f" First 5 bins (log10): {self._initial_log_dem[:5]}") + else: + print(" Initial DEM has not been estimated.") + + print("\nBASE DEM SOLUTION") + print("-" * 70) + if hasattr(self, "dem"): + print(f" DEM shape: {self.dem.shape}") + print(f" First 5 DEM bins: {self.dem[:5]}") + print(f" log10(DEM) first 5: {np.log10(self.dem[:5] + 1e-99)}") + print(f" Chi-square: {self.chisq:.4e}") + print(f" Modeled intensities: {self.modeled_intensities}") + else: + print(" No DEM solution computed yet (call solve()).") + + print("\nMONTE CARLO ENSEMBLE") + print("-" * 70) + if hasattr(self, "mc_dem"): + N = self.mc_dem.shape[0] - 1 + print(f" MC realizations: {N}") + if N > 0: + median = np.median(self.mc_dem[1:], axis=0) + p16, p84 = np.percentile(self.mc_dem[1:], [16, 84], axis=0) + print(" MC DEM statistics (first T-bin):") + print(f" Median (first 5): {median[:5]}") + print( + f" 1x bounds (log10): " + f"{np.log10(p16[0] + 1e-99):.2f} – {np.log10(p84[0] + 1e-99):.2f}" + ) + else: + print(" MC array allocated but N=0 (no Monte Carlo).") + else: + print(" No Monte Carlo results available.") + + print("\nPLOTTING HELPERS") + print("-" * 76) + print(" • plot_dem() – Base DEM only") + print(" • plot_dem_mc() – Base DEM + MC ensemble") + print("\n" + "=" * 76 + "\n") + + +XRTDEMIterative.plot_dem = dem_plotting.plot_dem +XRTDEMIterative.plot_dem_mc = dem_plotting.plot_dem_mc diff --git a/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/obs_20090730_DEM_MC100_IDL_2026.sav b/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/obs_20090730_DEM_MC100_IDL_2026.sav new file mode 100644 index 000000000..0115d0ed6 Binary files /dev/null and b/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/obs_20090730_DEM_MC100_IDL_2026.sav differ diff --git a/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/xrt_IDL_DEM_2012_10_27_MC100.sav b/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/xrt_IDL_DEM_2012_10_27_MC100.sav new file mode 100644 index 000000000..eac9fc1a1 Binary files /dev/null and b/xrtpy/xrt_dem_iterative/test/IDL_DEM_testing_sav_files/xrt_IDL_DEM_2012_10_27_MC100.sav differ diff --git a/xrtpy/xrt_dem_iterative/test/__init__.py b/xrtpy/xrt_dem_iterative/test/__init__.py new file mode 100644 index 000000000..2bd939b8a --- /dev/null +++ b/xrtpy/xrt_dem_iterative/test/__init__.py @@ -0,0 +1,8 @@ +""" +Test suite for xrtpy.xrt_dem_iterative and related DEM tools. + +This file marks the 'tests' directory as a Python package and allows +test modules to import from xrtpy and related packages cleanly. + +Use this space to define shared fixtures or helper utilities for tests if needed. +""" diff --git a/xrtpy/xrt_dem_iterative/test/test_dem_input_validation.py b/xrtpy/xrt_dem_iterative/test/test_dem_input_validation.py new file mode 100644 index 000000000..23de3ef5b --- /dev/null +++ b/xrtpy/xrt_dem_iterative/test/test_dem_input_validation.py @@ -0,0 +1,595 @@ +from pathlib import Path + +import astropy.units as u +import numpy as np +import pytest +from lmfit import Parameters +from scipy.io import readsav + +from xrtpy.response.channel import Channel +from xrtpy.response.tools import generate_temperature_responses +from xrtpy.xrt_dem_iterative import XRTDEMIterative + +channel_names = [ + "Al-mesh", + "Al-poly", + "C-poly", + "Ti-poly", + "Be-thin", + "Be-med", + "Al-med", + "Al-thick", + "Be-thick", + "Al-poly/Al-mesh", + "Al-poly/Ti-poly", + "Al-poly/Al-thick", + "Al-poly/Be-thick", + "C-poly/Ti-poly", +] + + +@pytest.mark.parametrize("channel_name", channel_names) +def test_channel_name(channel_name): + channel = Channel(channel_name) + assert channel.name == channel_name + + +def test_dem_validate_inputs_basic(): + filters = ["Al-poly", "Ti-poly"] + intensities = np.array([2500.0, 1800.0]) + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=8.0, + logarithmic_temperature_step_size=0.1, + monte_carlo_runs=0, + ) + + # Should NOT raise any error + x.validate_inputs() + + +def test_dem_temperature_grid(): + filters = ["Al-poly"] + intensities = np.array([1500.0]) + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=7.5, + logarithmic_temperature_step_size=0.1, + ) + + x.create_logT_grid() + + assert np.isclose(x.logT[0], 5.5) + assert np.isclose(x.logT[-1], 7.5) + assert len(x.logT) == 21 # (7.5-5.5)/0.1 + 1 = 21 + assert np.isclose(x.dlogT, 0.1) + assert np.isclose(x.dlnT, np.log(10) * 0.1) + + +def test_validate_inputs_good_case(): + filters = ["Be-thin", "Be-med"] + i_obs = [100.0, 200.0] + resp = generate_temperature_responses(filters, "2007-07-10") + dem = XRTDEMIterative(filters, i_obs, resp) + dem.validate_inputs() # Should NOT raise + + +def test_validate_inputs_mismatched_errors(): + filters = ["Be-thin", "Be-med"] + i_obs = [100.0, 200.0] + i_err = [100.0] # Wrong length - should be two error/ uncertainties + resp = generate_temperature_responses(filters, "2007-07-10") + dem = XRTDEMIterative(filters, i_obs, resp, intensity_errors=i_err) + with pytest.raises(ValueError, match="intensity_errors must match"): + dem.validate_inputs() + + +def test_create_logT_grid(): + + filters = ["Al-poly"] + intensities = np.array([150.0]) + responses = generate_temperature_responses(filters, "2018-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=7.5, + logarithmic_temperature_step_size=0.1, + ) + + x.create_logT_grid() + + # 1 — Correct start and end + assert x.logT[0] == pytest.approx(5.5) + assert x.logT[-1] == pytest.approx(7.5) + + # 2 — Correct number of bins: (7.5 - 5.5)/0.1 + 1 = 21 + assert len(x.logT) == 21 + assert x.n_bins == 21 + + # 3 — Correct spacing (uniform) + diffs = np.diff(x.logT) + assert np.allclose(diffs, 0.1, atol=1e-12) + + # 4 — dlogT and dlnT correct + assert x.dlogT == pytest.approx(0.1) + assert x.dlnT == pytest.approx(np.log(10) * 0.1) + + # 5 — T = 10**logT + assert np.allclose(x.T.to_value(u.K), 10**x.logT) + + # 6 — logT strictly increasing + assert np.all(np.diff(x.logT) > 0) + + +def test_estimate_initial_dem(): + filters = ["Al-poly", "Ti-poly"] + intensities = np.array([150.0, 230.0]) + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=7.5, + ) + + # Step 2: Create temperature grid & response matrix + x.create_logT_grid() + x._interpolate_responses_to_grid() + + # Step 3: Compute initial DEM + est = x._estimate_initial_dem() + + # TEST 1: Correct length + assert len(est) == len(x.logT) + + # TEST 2: All values should be exactly 0.0 ( Python implementation overrides with flat logDEM = 0) + assert np.allclose(est, 0.0) + + # TEST 3: Internal storage _initial_log_dem should match + assert np.allclose(x._initial_log_dem, est) + + # TEST 4: Returned DEM should be finite + assert np.all(np.isfinite(est)) + + +def test_prepare_spline_system(): + """ + 1. _prepare_spline_system runs without errors + 2. n_spl computed correctly + 3. spline_logT shape and monotonicity + 4. spline_log_dem has correct values + 5. pm_matrix has correct shape + 6. weights and abundances are all ones + """ + filters = ["Al-poly", "Ti-poly", "Be-thin"] + intensities = np.array([1000.0, 2000.0, 1500.0]) + responses = generate_temperature_responses(filters, "2012-10-27T12:30:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + ) + + # logT, response matrix, initial DEM + x.create_logT_grid() + x._interpolate_responses_to_grid() + x._estimate_initial_dem() # sets _initial_log_dem + + # Prepare spline system + x._prepare_spline_system() + + # TEST 1 — n_spl formula: n_channels=3 > n_spl=2 + assert x.n_spl == 2 + + # TEST 2 — spline_logT: correct shape and increasing + assert len(x.spline_logT) == 2 + assert np.all(np.diff(x.spline_logT) > 0) + + # TEST 3 — spline_log_dem should be zeros (flat initial DEM) + assert len(x.spline_log_dem) == 2 + assert np.allclose(x.spline_log_dem, 0.0) + + # TEST 4 — pm_matrix has correct shape + assert x.pm_matrix.shape == (3, len(x.logT)) + + # pm_matrix should be >= 0 + assert np.all(x.pm_matrix >= 0) + + # TEST 5 — weights and abundances + assert np.all(x.weights == 1.0) + assert np.all(x.abundances == 1.0) + + +def test_residuals_simple_case(): + """ + Create a fully synthetic DEM / response case so the forward model has a predictable value. + This isolates and tests the math inside `_residuals`. + """ + filters = ["Dummy"] + intensities = np.array([10.0]) # I_obs + responses = generate_temperature_responses(["Al-poly"], "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=6.5, + logarithmic_temperature_step_size=0.1, + ) + + # STEP 2 — Temperature grid + x.create_logT_grid() + N = len(x.logT) + + # Synthetic pm_matrix = constant 2 everywhere + x.pm_matrix = np.ones((1, N)) * 2.0 + + # STEP 3 — Construct synthetic spline state + x.spline_logT = np.array([x.logT[0], x.logT[-1]]) + x.spline_log_dem = np.array([0.0, 0.0]) # log10(DEM)=0 → DEM=1 + x.n_spl = 2 + + params = Parameters() + params.add("knot_0", value=0.0, min=-20, max=0) + params.add("knot_1", value=0.0, min=-20, max=0) + + # STEP 4 synthetic errors + x.intensities_scaled = np.array([10.0]) + x.sigma_scaled_intensity_errors = np.array([1.0]) + + # MISSING IN ORIGINAL TEST: Need abundances and weights (normally set in _prepare_spline_system) + x.abundances = np.ones(1) + x.weights = np.ones(1) + + # STEP 5 — Compute residuals + residuals = x._residuals(params) + + # Expected residual: + # DEM(T)=1 > pm(T)=2 > I_model = 2*N + # residual = (2N - I_obs) / sigma + expected_I_model = 2.0 * N + expected_residual = expected_I_model - 10.0 # sigma = 1 + + assert residuals.shape == (1,) + assert np.isfinite(residuals[0]) + assert np.isclose(residuals[0], expected_residual) + + +def test_solve_single_dem_zero_case(): + """ + If all observed intensities are zero, the solver must return: + - DEM = all zeros or value of normalization_factor + - modeled intensities = all zeros + - chi sqr = 0 + - result = None + This run will output a warning - expected due to no intensity_errors provided. + """ + # filterwarnings = ignore:No intensity_errors provided + + filters = ["Al-poly", "Ti-poly"] + intensities = np.array([0.0, 0.0]) # all zero → triggers nosolve + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + ) + + # STEP 2 — Build grid + response matrix + x.create_logT_grid() + x._interpolate_responses_to_grid() + + # STEP 3 — Call solver on zero intensities + dem, modeled, chi2, result = x._solve_single_dem( + observed_intensities_vals=intensities + ) + + # STEP 4 — Assertions matching IDL behavior -DEM must be zero everywhere + assert np.all(dem == 0.0) + + # modeled intensities must be all zero + assert np.all(modeled == 0.0) + + # chi² must be zero + assert chi2 == 0.0 + + # result object must be None (no lmfit minimization) + assert result is None + + # Shape correctness + assert dem.shape == x.logT.shape + assert modeled.shape == intensities.shape + + +def test_monte_carlo_different_realizations(): + """ + Monte Carlo DEM runs should produce DEMs that differ from the base DEM + when observational noise causes perturbed intensities to change. + + This test verifies: + • Monte Carlo output arrays have correct shapes + • Each perturbed DEM differs from the base DEM + • Perturbed intensities differ from the base intensities + • No invalid values appear (no NaNs, no negatives) + """ + + filters = ["Al-poly", "Ti-poly", "Be-thin"] + intensities = np.array([3000.0, 1500.0, 800.0], dtype=float) + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + N = 5 # small Monte Carlo batch for fast testing + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + minimum_bound_temperature=5.5, + maximum_bound_temperature=7.5, + logarithmic_temperature_step_size=0.5, + monte_carlo_runs=N, + ) + + # STEP 2: Run full DEM + Monte Carlo solver + x.solve() # this computes base and MC DEMs + + # STEP 3: Basic shape checks + n_T = len(x.logT) + n_obs = len(filters) + + # mc_dem shape: (N+1, n_T) + assert x.mc_dem.shape == (N + 1, n_T) + + # mc_base_obs shape: (N+1, n_obs) + assert x.mc_base_obs.shape == (N + 1, n_obs) + + # mc_mod_obs shape: (N+1, n_obs) + assert x.mc_mod_obs.shape == (N + 1, n_obs) + + # mc_chisq shape: (N+1,) + assert x.mc_chisq.shape == (N + 1,) + + # STEP 4: DEMs should *not* all be identical + base_dem = x.mc_dem[0] + + # At least one MC DEM must differ from base DEM + diffs = [ + not np.allclose(base_dem, x.mc_dem[ii], rtol=1e-5, atol=1e-8) + for ii in range(1, N + 1) + ] + + assert any( + diffs + ), "Monte Carlo realizations did not change the DEM; noise not applied?" + + # STEP 5: Perturbed observed intensities should differ + base_obs = x.mc_base_obs[0] + + obs_diffs = [not np.allclose(base_obs, x.mc_base_obs[ii]) for ii in range(1, N + 1)] + assert any( + obs_diffs + ), "Monte Carlo observed intensities identical; noise not applied?" + + # STEP 6: No invalid numbers after fitting + assert np.all(np.isfinite(x.mc_dem)) + assert np.all(x.mc_dem >= 0.0) + assert np.all(np.isfinite(x.mc_mod_obs)) + assert np.all(np.isfinite(x.mc_chisq)) + + +def test_reconstruct_dem_from_knots(): + """ + Test whether DEM reconstruction from spline knots behaves predictably. + + We construct: + • a synthetic temperature grid + • synthetic knots at logT endpoints + • synthetic log10(DEM) values at knots + And verify that: + • reconstructed logDEM matches knot values at endpoints + • DEM(T) is smooth in between (no NaNs or jumps) + • the DEM shape increases when knot values increase + """ + + filters = ["Dummy"] + intensities = np.array([1000.0]) + responses = generate_temperature_responses(["Al-poly"], "2012-10-27T00:00:00") + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + ) + + x.create_logT_grid() + + # STEP 2 — Define synthetic spline knot positions - Use ONLY endpoints (simplest nontrivial spline) + x.spline_logT = np.array([x.logT[0], x.logT[-1]]) + x.n_spl = 2 + + # STEP 3 — Create synthetic lmfit Parameters for knots - Case: log10(DEM) goes from -2 at low T to +1 at high T + params = Parameters() + params.add("knot_0", value=-2.0) # DEM = 1e-2 + params.add("knot_1", value=+1.0) # DEM = 1e+1 + + # Also required by reconstruct: store initial logDEM values + x.spline_log_dem = np.array([-2.0, 1.0]) + + # STEP 4 — Reconstruct DEM + dem = x._reconstruct_dem_from_knots(params) + + # STEP 5 — Assertions + # # Shape matches temperature grid + assert dem.shape == x.logT.shape + + # No NaNs or negatives + assert np.all(np.isfinite(dem)) + assert np.all(dem >= 0.0) + + # Endpoint matches exactly 10^knot_value + assert np.isclose(dem[0], 10 ** (-2.0), rtol=1e-6) + assert np.isclose(dem[-1], 10 ** (1.0), rtol=1e-6) + + # The DEM should increase monotonically between these endpoints + # (Cubic spline + monotonic knots → smooth monotonic increase) + assert ( + dem[0] < dem[len(dem) // 2] < dem[-1] + ), "Reconstructed DEM is not increasing between knots" + + +def test_full_pipeline_end_to_end(): + """ + Full DEM solving pipeline using real XRT filter responses. + Verifies: + • solve() runs without errors + • base DEM exists, finite, and positive + • modeled intensities computed + • chi-square finite + • Monte Carlo arrays created correctly (when N>0) + • No NaNs, no negative DEM, no shape mismatches + """ + + filters = ["Al-poly", "Ti-poly", "Be-thin", "C-poly"] + intensities = np.array([2500.0, 1800.0, 900.0, 450.0], dtype=float) + responses = generate_temperature_responses(filters, "2012-10-27T00:00:00") + N = 3 + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + monte_carlo_runs=N, + ) + + x.solve() + + # STEP 3: Base DEM checks + assert hasattr(x, "dem") + assert x.dem.shape == x.logT.shape + assert np.all(np.isfinite(x.dem)) + assert np.all(x.dem >= 0.0) + + # Peak must not be zero everywhere + assert np.max(x.dem) > 0.0 + + # STEP 4: Modeled intensities + assert hasattr(x, "modeled_intensities") + assert x.modeled_intensities.shape == intensities.shape + assert np.all(np.isfinite(x.modeled_intensities)) + + # STEP 5: Chi-square + assert hasattr(x, "chisq") + assert np.isfinite(x.chisq) + + # STEP 6: Monte Carlo arrays exist and are valid + assert x.mc_dem.shape == (N + 1, len(x.logT)) + assert x.mc_base_obs.shape == (N + 1, len(filters)) + assert x.mc_mod_obs.shape == (N + 1, len(filters)) + assert x.mc_chisq.shape == (N + 1,) + + # MC fields finite + assert np.all(np.isfinite(x.mc_dem)) + assert np.all(x.mc_dem >= 0.0) + assert np.all(np.isfinite(x.mc_base_obs)) + assert np.all(np.isfinite(x.mc_mod_obs)) + assert np.all(np.isfinite(x.mc_chisq)) + + # STEP 7: At least one MC DEM must differ from base DEM + base_dem = x.mc_dem[0] + different = [ + not np.allclose(base_dem, x.mc_dem[i], rtol=1e-5, atol=1e-8) + for i in range(1, N + 1) + ] + assert any(different), "Monte Carlo DEMs identical to base DEM — noise not applied?" + + # STEP 8: Modes must differ for perturbed cases + base_mod = x.mc_mod_obs[0] + different_mod = [ + not np.allclose(base_mod, x.mc_mod_obs[i]) for i in range(1, N + 1) + ] + assert any(different_mod) + +# ----------------------------------- TEST Against IDL ------------------------------------------ + + +def test_compare_with_idl_dem(): + """ + Compare xrtpy DEM solver output with reference IDL DEM stored in a .sav file. + The comparison is tolerant (allclose), since Python spline fitting and + MPFIT/LMFit differences mean we cannot expect bit-exact equality. + """ + + TEST_DIR = Path(__file__).parent + data_path = ( + TEST_DIR / "IDL_DEM_testing_sav_files" / "obs_20090730_DEM_MC100_IDL_2026.sav"# "xrt_IDL_DEM_2012_10_27_MC100.sav" + ) + data = readsav(data_path) + logT_idl = data["logt"] + + # XRTpy - 2012-10-27 Data Set + # filters = ["Ti-poly", "Be-thin", "Al-poly", "C-poly"] + # intensities = np.array([311680.2, 135815.0, 2351258.9, 2352987.7]) + # date = "2012-10-27T16:27:46" + + filters = ["Al-mesh", "Ti-poly", "Al-poly", "Be-thin"] + intensities = [178.482 ,44.919,132.193,3.149] # DN/s + observation_date="2009-07-30T00:38" + + responses = generate_temperature_responses(filters, observation_date) + + x = XRTDEMIterative( + observed_channel=filters, + observed_intensities=intensities, + temperature_responses=responses, + monte_carlo_runs=0, # IDL comparison → no MC + ) + + # STEP 3 — Run Python DEM + x.solve() # linear DEM + logT_python = x.logT + + # STEP 4 — Compare logT grids - IDL usually uses exactly the same grid, but check with tolerance + assert np.allclose( + logT_python, logT_idl, atol=1e-6 + ), "Temperature grids differ significantly between IDL and Python." + + # #log_dem_python = np.log10(np.maximum(dem_python, 1e-99)) + # #dem_idl = data["dem"] + # #log_dem_idl = np.log10(np.maximum(dem_idl, 1e-99)) + # # STEP 5 — Compare DEM values + # # Numerical tolerances: + # # - 0.3 dex tolerance → factor of ~2 + # # - Acceptable because: + # # • Spline fits differ slightly between MPFIT and LMFit + # # • IDL interpolation behavior differs from CubicSpline + # # • Floating-point roundoff differences + # tol_dex = 0.9 + + # diff = np.abs(log_dem_python - log_dem_idl) + + # assert np.all(diff < tol_dex), ( + # "DEM shape diverges too far from IDL reference.\n" + # f"Max difference {diff.max():.3f} dex (allowed {tol_dex})." + # ) + + # # Ensure the overall trends match: peak in same region + # peak_idl = np.argmax(log_dem_idl) + # peak_python = np.argmax(log_dem_python) + # assert abs(peak_idl - peak_python) <= 1, \ + # "DEM peak location differs more than 1 temperature bin."