From 4a5a544f8b0361eec8170338d68b609b56296d0c Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Thu, 12 Sep 2024 11:57:17 +0100 Subject: [PATCH 01/21] First draft --- sunkit_spex/models/physical/albedo.py | 81 +++++++++++++++++++ .../models/physical/tests/test_albedo.py | 13 +++ 2 files changed, 94 insertions(+) create mode 100644 sunkit_spex/models/physical/albedo.py create mode 100644 sunkit_spex/models/physical/tests/test_albedo.py diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py new file mode 100644 index 00000000..78905103 --- /dev/null +++ b/sunkit_spex/models/physical/albedo.py @@ -0,0 +1,81 @@ +import numpy as np +from scipy.interpolate import RegularGridInterpolator +from scipy.io import readsav + +import astropy.units as u +from astropy.units import Quantity + +from sunpy.data import cache + + +@u.quantity_input +def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotropy=1): + r""" + Add albedo correction to input spectrum + + Parameters + ---------- + spec : + Input spectrum + energy : + Energy edges associated with the spectrum + theta : + Angle between Sun-observer line and X-ray source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, anisotropy=1 (default) the source is isotropic + + Example + ------- + >>> import astropy.units as u + >>> import numpy as np + >>> from sunkit_spex.models.physical import albedo + >>> e = np.linspace(5, 500, 1000) + >>> e_c = e[:-1] + np.diff(e) + >>> s = 125*e_c**-3 + >>> corrected = albedo(s, energy, theta=45*u.deg) + """ + base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" + mu = np.cos(theta) + + # what bout 0 and 1 i assume so close to 05 and 95 that it doesn't matter + # load precomputed green matrices + if 0.5 <= mu <= 0.95: + low = 5 * np.floor(mu * 20) + high = 5 * np.ceil(mu * 20) + low_name = f"green_compton_mu{low:03.0f}.dat" + high_name = f"green_compton_mu{high:03.0f}.dat" + low_file = cache.download(base_url + low_name) + high_file = cache.download(base_url + high_name) + green = readsav(low_file) + albedo_low = green["p"].albedo[0] + green_high = readsav(high_file) + albedo_high = green_high["p"].albedo[0] + # why 20? + albedo = albedo_low + (albedo_high - albedo_low) * (mu - (np.floor(mu * 20)) / 20) + + elif mu < 0.5: + file = "green_compton_mu005.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + elif mu > 0.95: + file = "green_compton_mu095.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + + albedo = albedo.T + energy_grid_edges = green["p"].edges[0] * u.keV + energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) + + interp = RegularGridInterpolator((energy_grid_centers.to_value(u.keV), energy_grid_centers.to_value(u.keV)), albedo) + + de = np.diff(energy) + engery_centers = energy[:-1] + de / 2 + + X, Y = np.meshgrid(engery_centers.to_value(u.keV), engery_centers.to_value(u.keV)) + albedo_interp = interp((X, Y)) + + albedo_interp = albedo_interp * de.value / anisotropy + + return spec + spec @ albedo_interp.T diff --git a/sunkit_spex/models/physical/tests/test_albedo.py b/sunkit_spex/models/physical/tests/test_albedo.py new file mode 100644 index 00000000..c36106af --- /dev/null +++ b/sunkit_spex/models/physical/tests/test_albedo.py @@ -0,0 +1,13 @@ +import numpy as np + +import astropy.units as u + +from sunkit_spex.models.physical.albedo import albedo + + +def test_albedo(): + e = np.linspace(5, 500, 250) + e_c = e[:-1] + np.diff(e) + s = 125 * e_c**-3 + out = albedo(s, e * u.keV, 45 * u.deg) + print(out) From 087ea163b508999f2362015677a263b2493b3529 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Thu, 19 Sep 2024 16:11:56 +0100 Subject: [PATCH 02/21] eh not sure --- sunkit_spex/models/physical/albedo.py | 13 +++++++++---- sunkit_spex/models/physical/tests/test_albedo.py | 2 +- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 78905103..a5c532f1 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -13,6 +13,11 @@ def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotr r""" Add albedo correction to input spectrum + Correct input model spectrum for the component reflected by the solar atmosphere following [Kontar20006]_ using + precomputed green matrices from SSW. + + .. [Kontar20006] https://doi.org/10.1051/0004-6361:20053672 + Parameters ---------- spec : @@ -32,12 +37,12 @@ def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotr >>> e = np.linspace(5, 500, 1000) >>> e_c = e[:-1] + np.diff(e) >>> s = 125*e_c**-3 - >>> corrected = albedo(s, energy, theta=45*u.deg) + >>> corrected = albedo(s, e, theta=45*u.deg) """ base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" mu = np.cos(theta) - # what bout 0 and 1 i assume so close to 05 and 95 that it doesn't matter + # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter # load precomputed green matrices if 0.5 <= mu <= 0.95: low = 5 * np.floor(mu * 20) @@ -71,9 +76,9 @@ def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotr interp = RegularGridInterpolator((energy_grid_centers.to_value(u.keV), energy_grid_centers.to_value(u.keV)), albedo) de = np.diff(energy) - engery_centers = energy[:-1] + de / 2 + energy_centers = energy[:-1] + de / 2 - X, Y = np.meshgrid(engery_centers.to_value(u.keV), engery_centers.to_value(u.keV)) + X, Y = np.meshgrid(energy_centers.to_value(u.keV), energy_centers.to_value(u.keV)) albedo_interp = interp((X, Y)) albedo_interp = albedo_interp * de.value / anisotropy diff --git a/sunkit_spex/models/physical/tests/test_albedo.py b/sunkit_spex/models/physical/tests/test_albedo.py index c36106af..58911aca 100644 --- a/sunkit_spex/models/physical/tests/test_albedo.py +++ b/sunkit_spex/models/physical/tests/test_albedo.py @@ -6,7 +6,7 @@ def test_albedo(): - e = np.linspace(5, 500, 250) + e = np.linspace(4, 600, 597) e_c = e[:-1] + np.diff(e) s = 125 * e_c**-3 out = albedo(s, e * u.keV, 45 * u.deg) From f7fc55d37b26580af7e0aaf24adf183e22b82d80 Mon Sep 17 00:00:00 2001 From: natsat0919 Date: Thu, 19 Dec 2024 19:38:38 +0000 Subject: [PATCH 03/21] Speed up albedo correction, add AlbedoModel --- sunkit_spex/models/physical/albedo.py | 153 +++++++++++++----- .../models/physical/tests/test_albedo.py | 7 +- 2 files changed, 121 insertions(+), 39 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index a5c532f1..63db9b50 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -1,47 +1,62 @@ +import warnings +import functools + import numpy as np from scipy.interpolate import RegularGridInterpolator from scipy.io import readsav -import astropy.units as u -from astropy.units import Quantity +# from astropy.units import Quantity +from astropy.modeling import Fittable1DModel, Parameter from sunpy.data import cache +__all__ = ["AlbedoModel"] -@u.quantity_input -def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotropy=1): - r""" - Add albedo correction to input spectrum - Correct input model spectrum for the component reflected by the solar atmosphere following [Kontar20006]_ using - precomputed green matrices from SSW. +class AlbedoModel(Fittable1DModel): + def __init__(self, energy, anisotropy, theta): + self.energy = energy + self.anisotropy = Parameter( + default=anisotropy, description="The anisotropy used for albedo correction", fixed=True + ) + self.theta = Parameter( + default=theta, description="Angle between the flare and the telescope in radians", fixed=True + ) + super().__init__() - .. [Kontar20006] https://doi.org/10.1051/0004-6361:20053672 + def evaluate(self, count_model): + """Corrects the composite count model for albedo. To be used as: ct_model = (ph_model|srm)|albedo""" + albedo_matrix_T = albedo(self.energy, self.theta.value, self.anisotropy.value) - Parameters - ---------- - spec : - Input spectrum - energy : - Energy edges associated with the spectrum - theta : - Angle between Sun-observer line and X-ray source - anisotropy : - Ratio of the flux in observer direction to the flux downwards, anisotropy=1 (default) the source is isotropic + return count_model + count_model @ albedo_matrix_T - Example - ------- - >>> import astropy.units as u - >>> import numpy as np - >>> from sunkit_spex.models.physical import albedo - >>> e = np.linspace(5, 500, 1000) - >>> e_c = e[:-1] + np.diff(e) - >>> s = 125*e_c**-3 - >>> corrected = albedo(s, e, theta=45*u.deg) - """ - base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" + +# Wrapper for chaching function with numpy array args +def np_cache(function): + @functools.cache + def cached_wrapper(*args, **kwargs): + args = [np.array(a) if isinstance(a, tuple) else a for a in args] + kwargs = {k: np.array(v) if isinstance(v, tuple) else v for k, v in kwargs.items()} + + return function(*args, **kwargs) + + @functools.wraps(function) + def wrapper(*args, **kwargs): + args = [tuple(a) if isinstance(a, np.ndarray) else a for a in args] + kwargs = {k: tuple(v) if isinstance(v, np.ndarray) else v for k, v in kwargs.items()} + return cached_wrapper(*args, **kwargs) + + wrapper.cache_info = cached_wrapper.cache_info + wrapper.cache_clear = cached_wrapper.cache_clear + + return wrapper + + +@functools.cache +def load_green_matrices(theta): mu = np.cos(theta) + base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter # load precomputed green matrices if 0.5 <= mu <= 0.95: @@ -70,17 +85,83 @@ def albedo(spec, energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotr albedo = green["p"].albedo[0] albedo = albedo.T - energy_grid_edges = green["p"].edges[0] * u.keV + + # energy_grid_edges = green["p"].edges[0] * u.keV + energy_grid_edges = green["p"].edges[0] + energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) - interp = RegularGridInterpolator((energy_grid_centers.to_value(u.keV), energy_grid_centers.to_value(u.keV)), albedo) + # interp_green_matrix = RegularGridInterpolator((energy_grid_centers.to_value(u.keV), energy_grid_centers.to_value(u.keV)), albedo) + interp_green_matrix = RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) + + return interp_green_matrix + +@np_cache +def calc_albedo_matrix(energy, interp_green_matrix, anisotropy): + # Re-introduce units + # energy = energy * u.keV de = np.diff(energy) energy_centers = energy[:-1] + de / 2 - X, Y = np.meshgrid(energy_centers.to_value(u.keV), energy_centers.to_value(u.keV)) - albedo_interp = interp((X, Y)) + # X, Y = np.meshgrid(energy_centers.to_value(u.keV), energy_centers.to_value(u.keV)) + X, Y = np.meshgrid(energy_centers, energy_centers) + + albedo_interp = interp_green_matrix((X, Y)) + + # Scale by anisotropy + # albedo_interp = albedo_interp * de.value / anisotropy + albedo_interp = albedo_interp * de / anisotropy + + # Take a transpose + albedo_interp_T = albedo_interp.T + + return albedo_interp_T + + +# @u.quantity_input +# def albedo(energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotropy=1): +def albedo(energy, theta, anisotropy=1): + r""" + Add albedo correction to input spectrum + + Correct input model spectrum for the component reflected by the solar atmosphere following [Kontar20006]_ using + precomputed green matrices from SSW. + + .. [Kontar20006] https://doi.org/10.1051/0004-6361:20053672 + + Parameters + ---------- + spec : + Input count spectrum + energy : + Energy edges associated with the spectrum + theta : + Angle between Sun-observer line and X-ray source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, anisotropy=1 (default) the source is isotropic + + Example + ------- + >>> import astropy.units as u + >>> import numpy as np + >>> from sunkit_spex.models.physical.albedo import albedo + >>> e = np.linspace(5, 500, 1000) + >>> e_c = e[:-1] + np.diff(e) + >>> s = 125*e_c**-3 + >>> albedo_matrix = albedo(e, theta=0.2) + >>> corrected = s + s@albedo_matrix + """ + + # Add check for energy range restricted by the Green matrices + if energy[0] < 3 or energy[-1] > 600: + warnings.warn("\nCount energies are required to be >= 3keV and <=600 keV ") + + # Green matrix for a given theta + interp_green = load_green_matrices(theta) - albedo_interp = albedo_interp * de.value / anisotropy + # Transpose of a matrix used for albedo correction, need to strip energy of units otherwise problem with caching + # albedo_matrix_T = calc_albedo_matrix(energy.value, interp_green, anisotropy) + albedo_matrix_T = calc_albedo_matrix(energy, interp_green, anisotropy) - return spec + spec @ albedo_interp.T + return albedo_matrix_T diff --git a/sunkit_spex/models/physical/tests/test_albedo.py b/sunkit_spex/models/physical/tests/test_albedo.py index 58911aca..5d33964c 100644 --- a/sunkit_spex/models/physical/tests/test_albedo.py +++ b/sunkit_spex/models/physical/tests/test_albedo.py @@ -1,7 +1,5 @@ import numpy as np -import astropy.units as u - from sunkit_spex.models.physical.albedo import albedo @@ -9,5 +7,8 @@ def test_albedo(): e = np.linspace(4, 600, 597) e_c = e[:-1] + np.diff(e) s = 125 * e_c**-3 - out = albedo(s, e * u.keV, 45 * u.deg) + # e = e *u.keV + theta = 0.2 # u.rad + albedo_matrix = albedo(e, theta=theta) + out = s + s @ albedo_matrix print(out) From 503ff07c6b64c02c16bbe9916124cf3d35c21eca Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Thu, 26 Sep 2024 19:15:07 +0100 Subject: [PATCH 04/21] WIP --- sunkit_spex/models/physical/albedo.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 63db9b50..dfe4516b 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -1,6 +1,8 @@ import warnings import functools +from functools import lru_cache + import numpy as np from scipy.interpolate import RegularGridInterpolator from scipy.io import readsav @@ -128,7 +130,7 @@ def albedo(energy, theta, anisotropy=1): Correct input model spectrum for the component reflected by the solar atmosphere following [Kontar20006]_ using precomputed green matrices from SSW. - .. [Kontar20006] https://doi.org/10.1051/0004-6361:20053672 + .. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672 Parameters ---------- From b5c88953ae646551b6484e3ad06d26c22b0edaca Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 20 Dec 2024 16:25:17 +0000 Subject: [PATCH 05/21] Tidyup, tests and docs --- docs/reference/fitting.rst | 10 + docs/reference/index.rst | 17 +- docs/reference/legacy.rst | 13 ++ docs/reference/models.rst | 8 + sunkit_spex/models/physical/albedo.py | 200 +++++++++++------- .../models/physical/tests/test_albedo.py | 45 +++- 6 files changed, 192 insertions(+), 101 deletions(-) create mode 100644 docs/reference/fitting.rst create mode 100644 docs/reference/legacy.rst create mode 100644 docs/reference/models.rst diff --git a/docs/reference/fitting.rst b/docs/reference/fitting.rst new file mode 100644 index 00000000..ba24555e --- /dev/null +++ b/docs/reference/fitting.rst @@ -0,0 +1,10 @@ +Fitting (`sunkit_spex.fitting`) +******************************* + +``sunkit_spex.fitting`` contains the sunkit-spex fitting infrastructure including statistics, optimisers and fitters. + + +.. automodapi:: sunkit_spex.fitting.objective_functions +.. automodapi:: sunkit_spex.fitting.optimizer_tools +.. automodapi:: sunkit_spex.fitting.statistics +.. automodapi:: sunkit_spex.fitting diff --git a/docs/reference/index.rst b/docs/reference/index.rst index a77f7489..8537a8a3 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -4,13 +4,10 @@ Reference Software and API. -.. automodapi:: sunkit_spex -.. automodapi:: sunkit_spex.legacy.thermal -.. automodapi:: sunkit_spex.legacy.integrate -.. automodapi:: sunkit_spex.legacy.io -.. automodapi:: sunkit_spex.legacy.emission -.. automodapi:: sunkit_spex.legacy.constants -.. automodapi:: sunkit_spex.legacy.fitting -.. automodapi:: sunkit_spex.legacy.fitting.io -.. automodapi:: sunkit_spex.legacy.photon_power_law -.. automodapi:: sunkit_spex.extern.rhessi + +.. toctree:: + :maxdepth: 2 + + fitting/index + models/index + legacy/index diff --git a/docs/reference/legacy.rst b/docs/reference/legacy.rst new file mode 100644 index 00000000..71c14f23 --- /dev/null +++ b/docs/reference/legacy.rst @@ -0,0 +1,13 @@ +Legacy (`sunkit_spex.legacy`) +***************************** + + +.. automodapi:: sunkit_spex.legacy +.. automodapi:: sunkit_spex.legacy.thermal +.. automodapi:: sunkit_spex.legacy.integrate +.. automodapi:: sunkit_spex.legacy.io +.. automodapi:: sunkit_spex.legacy.emission +.. automodapi:: sunkit_spex.legacy.constants +.. automodapi:: sunkit_spex.legacy.fitting_legacy +.. automodapi:: sunkit_spex.legacy.fitting_legacy.io +.. automodapi:: sunkit_spex.legacy.photon_power_law diff --git a/docs/reference/models.rst b/docs/reference/models.rst new file mode 100644 index 00000000..7d71a80c --- /dev/null +++ b/docs/reference/models.rst @@ -0,0 +1,8 @@ +Models (`sunkit_spex.models`) +***************************** + +``sunkit_spex.models`` contains a number of functional and physical models. + + +.. automodapi:: sunkit_spex.models +.. automodapi:: sunkit_spex.models.physical diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index dfe4516b..4c7c2880 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -1,61 +1,93 @@ -import warnings -import functools - from functools import lru_cache import numpy as np +from numpy.typing import NDArray from scipy.interpolate import RegularGridInterpolator from scipy.io import readsav -# from astropy.units import Quantity -from astropy.modeling import Fittable1DModel, Parameter +import astropy.units as u +from astropy.modeling import FittableModel, Parameter +from astropy.units import Quantity from sunpy.data import cache -__all__ = ["AlbedoModel"] +__all__ = ["Albedo", "get_albedo_matrix"] -class AlbedoModel(Fittable1DModel): - def __init__(self, energy, anisotropy, theta): - self.energy = energy - self.anisotropy = Parameter( - default=anisotropy, description="The anisotropy used for albedo correction", fixed=True - ) - self.theta = Parameter( - default=theta, description="Angle between the flare and the telescope in radians", fixed=True - ) - super().__init__() +class Albedo(FittableModel): + r""" + Aldedo model which adds albdeo correction to input spectrum. - def evaluate(self, count_model): - """Corrects the composite count model for albedo. To be used as: ct_model = (ph_model|srm)|albedo""" - albedo_matrix_T = albedo(self.energy, self.theta.value, self.anisotropy.value) + Following [Kontar2006]_ using precomputed green matrices distributed as part of SSW_. - return count_model + count_model @ albedo_matrix_T + .. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672 + .. [SSW] https://www.lmsal.com/solarsoft/ + Parameters + ========== + energy_edges : + Energy edges associated with input spectrum + theta : + Angle between Sun-observer line and X-ray source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source + + Examples + ======== + >>> import astropy.units as u + >>> import numpy as np + >>> from astropy.modeling.powerlaws import PowerLaw1D + >>> from sunkit_spex.models.physical.albedo import Albedo + >>> e_edges = np.linspace(10, 300, 10) * u.keV + >>> e_centers = e_edges[0:-1] + (0.5 * np.diff(e_edges)) + >>> source = PowerLaw1D(amplitude=100*u.ph, x_0=10*u.keV, alpha=4) + >>> observed = source | Albedo(energy_edges=e_edges) + >>> observed(e_centers) + + """ + + theta = Parameter( + name="theta", + default=0, + unit=u.deg, + min=-90, + max=90, + description="Angle between the observer and the source", + fixed=True, + ) + anisotropy = Parameter(default=1, description="The anisotropy used for albedo correction", fixed=True) + + n_inputs = 1 + n_outputs = 1 -# Wrapper for chaching function with numpy array args -def np_cache(function): - @functools.cache - def cached_wrapper(*args, **kwargs): - args = [np.array(a) if isinstance(a, tuple) else a for a in args] - kwargs = {k: np.array(v) if isinstance(v, tuple) else v for k, v in kwargs.items()} + def __init__(self, *args, **kwargs): + self.energy_edges = kwargs.pop("energy_edges") + return super().__init__(*args, **kwargs) - return function(*args, **kwargs) + def evaluate(self, spectrum, theta, anisotropy): + albedo_matrix_T = get_albedo_matrix(self.energy_edges, theta, anisotropy) + return spectrum + spectrum @ albedo_matrix_T - @functools.wraps(function) - def wrapper(*args, **kwargs): - args = [tuple(a) if isinstance(a, np.ndarray) else a for a in args] - kwargs = {k: tuple(v) if isinstance(v, np.ndarray) else v for k, v in kwargs.items()} - return cached_wrapper(*args, **kwargs) - wrapper.cache_info = cached_wrapper.cache_info - wrapper.cache_clear = cached_wrapper.cache_clear +@lru_cache +def _get_green_matrix(theta: float) -> RegularGridInterpolator: + r""" + Get greens matrix for given angle. - return wrapper + Interpolates pre-computed green matrices for fixed angles. The resulting greens matrix is then loaded into an + interpolator for later energy interpolation. + Parameters + ========== + theta : float + Angle between the observer and the source -@functools.cache -def load_green_matrices(theta): + Returns + ======= + Greens matrix interpolator + """ mu = np.cos(theta) base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" @@ -88,32 +120,41 @@ def load_green_matrices(theta): albedo = albedo.T - # energy_grid_edges = green["p"].edges[0] * u.keV + # By construction in keV energy_grid_edges = green["p"].edges[0] - energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) - # interp_green_matrix = RegularGridInterpolator((energy_grid_centers.to_value(u.keV), energy_grid_centers.to_value(u.keV)), albedo) - interp_green_matrix = RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) + green_matrix_interpolator = RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) + + return green_matrix_interpolator - return interp_green_matrix +@lru_cache +def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotropy: float) -> NDArray: + r""" + Calculate green matrix for given energies and angle. -@np_cache -def calc_albedo_matrix(energy, interp_green_matrix, anisotropy): - # Re-introduce units - # energy = energy * u.keV - de = np.diff(energy) - energy_centers = energy[:-1] + de / 2 + Interpolates precomputed green matrices for given energies and angle. + + Parameters + ========== + energy_edges : + Energy edges associated with the spectrum + theta : + Angle between the observer and the source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source + """ + albedo_interpolator = _get_green_matrix(theta) + de = np.diff(energy_edges) + energy_centers = energy_edges[:-1] + de / 2 - # X, Y = np.meshgrid(energy_centers.to_value(u.keV), energy_centers.to_value(u.keV)) X, Y = np.meshgrid(energy_centers, energy_centers) - albedo_interp = interp_green_matrix((X, Y)) + albedo_interp = albedo_interpolator((X, Y)) # Scale by anisotropy - # albedo_interp = albedo_interp * de.value / anisotropy - albedo_interp = albedo_interp * de / anisotropy + albedo_interp = (albedo_interp * de) / anisotropy # Take a transpose albedo_interp_T = albedo_interp.T @@ -121,49 +162,46 @@ def calc_albedo_matrix(energy, interp_green_matrix, anisotropy): return albedo_interp_T -# @u.quantity_input -# def albedo(energy: Quantity[u.keV], theta: Quantity[u.deg, 0, 90], anisotropy=1): -def albedo(energy, theta, anisotropy=1): +@u.quantity_input +def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], anisotropy=1): r""" - Add albedo correction to input spectrum + Get albedo correction matrix. - Correct input model spectrum for the component reflected by the solar atmosphere following [Kontar20006]_ using - precomputed green matrices from SSW. + Matrix used to correct a photon spectrum for the component reflected by the solar atmosphere following + [Kontar2006]_ using precomputed green matrices distributed as part of SSW_. .. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672 + .. [SSW] https://www.lmsal.com/solarsoft/ Parameters ---------- - spec : - Input count spectrum - energy : + energy_edges : Energy edges associated with the spectrum theta : Angle between Sun-observer line and X-ray source anisotropy : - Ratio of the flux in observer direction to the flux downwards, anisotropy=1 (default) the source is isotropic + Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source Example ------- >>> import astropy.units as u >>> import numpy as np - >>> from sunkit_spex.models.physical.albedo import albedo - >>> e = np.linspace(5, 500, 1000) - >>> e_c = e[:-1] + np.diff(e) - >>> s = 125*e_c**-3 - >>> albedo_matrix = albedo(e, theta=0.2) - >>> corrected = s + s@albedo_matrix + >>> from sunkit_spex.models.physical.albedo import get_albedo_matrix + >>> e = np.linspace(5, 500, 5)*u.keV + >>> albedo_matrix = get_albedo_matrix(e,theta=45*u.deg) + >>> albedo_matrix + array([[3.80274484e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], + [5.10487362e-01, 8.35309813e-07, 0.00000000e+00, 0.00000000e+00], + [3.61059918e-01, 2.48711099e-01, 2.50744411e-09, 0.00000000e+00], + [3.09323903e-01, 2.66485260e-01, 1.23563372e-01, 1.81846722e-10]]) """ - - # Add check for energy range restricted by the Green matrices - if energy[0] < 3 or energy[-1] > 600: - warnings.warn("\nCount energies are required to be >= 3keV and <=600 keV ") - - # Green matrix for a given theta - interp_green = load_green_matrices(theta) - - # Transpose of a matrix used for albedo correction, need to strip energy of units otherwise problem with caching - # albedo_matrix_T = calc_albedo_matrix(energy.value, interp_green, anisotropy) - albedo_matrix_T = calc_albedo_matrix(energy, interp_green, anisotropy) - - return albedo_matrix_T + if energy_edges[0].to_value(u.keV) < 3 or energy_edges[-1].to_value(u.keV) > 600: + raise ValueError("Supported energy range 3 <= E <= 600 keV") + theta = np.array(theta).squeeze() << theta.unit + if np.abs(theta) > 90 * u.deg: + raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") + anisotropy = np.array(anisotropy).squeeze() + albedo_matrix = _calculate_albedo_matrix( + tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item() + ) + return albedo_matrix diff --git a/sunkit_spex/models/physical/tests/test_albedo.py b/sunkit_spex/models/physical/tests/test_albedo.py index 5d33964c..542d9ada 100644 --- a/sunkit_spex/models/physical/tests/test_albedo.py +++ b/sunkit_spex/models/physical/tests/test_albedo.py @@ -1,14 +1,39 @@ import numpy as np +import pytest -from sunkit_spex.models.physical.albedo import albedo +import astropy.units as u +from astropy.modeling.powerlaws import PowerLaw1D +from astropy.units import UnitsError +from sunkit_spex.models.physical.albedo import Albedo, get_albedo_matrix -def test_albedo(): - e = np.linspace(4, 600, 597) - e_c = e[:-1] + np.diff(e) - s = 125 * e_c**-3 - # e = e *u.keV - theta = 0.2 # u.rad - albedo_matrix = albedo(e, theta=theta) - out = s + s @ albedo_matrix - print(out) + +def test_get_albedo_matrix(): + e = np.linspace(4, 600, 597) * u.keV + theta = 0 * u.deg + albedo_matrix = get_albedo_matrix(e, theta=theta) + assert albedo_matrix[0, 0] == 0.006154127884656191 + assert albedo_matrix[298, 298] == 5.956079300274577e-24 + assert albedo_matrix[-1, -1] == 2.3302891436400413e-27 + + +def test_get_albedo_matrix_bad_energy(): + e = [1, 4, 10, 20] * u.keV + with pytest.raises(ValueError, match=r"Supported energy range 3.*"): + get_albedo_matrix(e, theta=0 * u.deg) + + +def test_get_albedo_matrix_bad_angle(): + e = [4, 10, 20] * u.keV + with pytest.raises(UnitsError, match=r".*Argument 'theta' to function 'get_albedo_matrix'.*'deg'."): + get_albedo_matrix(e, theta=91 * u.m) + with pytest.raises(ValueError, match=r"Theta must be between -90 and 90 degrees.*"): + get_albedo_matrix(e, theta=-91 * u.deg) + + +def test_albedo_model(): + e_edges = np.linspace(10, 300, 10) * u.keV + e_centers = e_edges[0:-1] + (0.5 * np.diff(e_edges)) + source = PowerLaw1D(amplitude=100 * u.ph, x_0=10 * u.keV, alpha=4) + observed = source | Albedo(energy_edges=e_edges) + observed(e_centers) From 6ee7756e235f745c83e290ad8deada7a3e3a16ad Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 27 Dec 2024 13:39:37 +0000 Subject: [PATCH 06/21] Add model and tidy up code --- changelog/161.feature.rst | 1 + docs/conf.py | 1 + docs/index.rst | 2 +- docs/reference/index.rst | 8 ++- docs/reference/legacy.rst | 6 +- docs/reference/models.rst | 4 ++ sunkit_spex/models/physical/albedo.py | 80 +++++++++++++++------------ 7 files changed, 61 insertions(+), 41 deletions(-) create mode 100644 changelog/161.feature.rst diff --git a/changelog/161.feature.rst b/changelog/161.feature.rst new file mode 100644 index 00000000..37333b52 --- /dev/null +++ b/changelog/161.feature.rst @@ -0,0 +1 @@ +Add function `~sunkit_spex.models.physical.albedo.get_albedo_matrix` to calculate Albedo correction for given input spectrum and add model `~sunkit_spex.models.physical.albedo.Albedo` to correct modeled photon spectrum for albedo. diff --git a/docs/conf.py b/docs/conf.py index 07261173..26b1cd61 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -48,6 +48,7 @@ "sphinx_automodapi.smart_resolver", "sphinx_changelog", "sphinx_gallery.gen_gallery", + "matplotlib.sphinxext.plot_directive" ] # Add any paths that contain templates here, relative to this directory. diff --git a/docs/index.rst b/docs/index.rst index 007f3393..82e91fdf 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -16,7 +16,7 @@ will help you know where to look for certain things: * :doc:`Topic guides ` discuss key topics and concepts at a fairly high level and provide useful background information and explanation. -* :doc:`Reference guides ` contain technical reference for APIs and +* :doc:`Reference ` contain technical reference for APIs and other aspects of sunkit-spex machinery. They describe how it works and how to use it but assume that you have a basic understanding of key concepts. diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 8537a8a3..6c60dd2d 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -4,10 +4,12 @@ Reference Software and API. +.. automodapi:: sunkit_spex .. toctree:: :maxdepth: 2 - fitting/index - models/index - legacy/index + + fitting + models + legacy diff --git a/docs/reference/legacy.rst b/docs/reference/legacy.rst index 71c14f23..733fc562 100644 --- a/docs/reference/legacy.rst +++ b/docs/reference/legacy.rst @@ -1,6 +1,8 @@ Legacy (`sunkit_spex.legacy`) ***************************** +.. warning:: + The legacy module contains legacy code which will no longer be maintained and will be removed in the near future. .. automodapi:: sunkit_spex.legacy .. automodapi:: sunkit_spex.legacy.thermal @@ -8,6 +10,6 @@ Legacy (`sunkit_spex.legacy`) .. automodapi:: sunkit_spex.legacy.io .. automodapi:: sunkit_spex.legacy.emission .. automodapi:: sunkit_spex.legacy.constants -.. automodapi:: sunkit_spex.legacy.fitting_legacy -.. automodapi:: sunkit_spex.legacy.fitting_legacy.io +.. automodapi:: sunkit_spex.legacy.fitting +.. automodapi:: sunkit_spex.legacy.fitting.io .. automodapi:: sunkit_spex.legacy.photon_power_law diff --git a/docs/reference/models.rst b/docs/reference/models.rst index 7d71a80c..c588800e 100644 --- a/docs/reference/models.rst +++ b/docs/reference/models.rst @@ -5,4 +5,8 @@ Models (`sunkit_spex.models`) .. automodapi:: sunkit_spex.models +.. automodapi:: sunkit_spex.models.models +.. automodapi:: sunkit_spex.models.instrument_response .. automodapi:: sunkit_spex.models.physical +.. automodapi:: sunkit_spex.models.physical.thermal +.. automodapi:: sunkit_spex.models.physical.albedo diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 4c7c2880..f60d68e3 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -18,7 +18,7 @@ class Albedo(FittableModel): r""" Aldedo model which adds albdeo correction to input spectrum. - Following [Kontar2006]_ using precomputed green matrices distributed as part of SSW_. + Following [Kontar2006]_ using precomputed green matrices distributed as part of [SSW]_. .. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672 .. [SSW] https://www.lmsal.com/solarsoft/ @@ -34,20 +34,43 @@ class Albedo(FittableModel): Examples ======== - >>> import astropy.units as u - >>> import numpy as np - >>> from astropy.modeling.powerlaws import PowerLaw1D - >>> from sunkit_spex.models.physical.albedo import Albedo - >>> e_edges = np.linspace(10, 300, 10) * u.keV - >>> e_centers = e_edges[0:-1] + (0.5 * np.diff(e_edges)) - >>> source = PowerLaw1D(amplitude=100*u.ph, x_0=10*u.keV, alpha=4) - >>> observed = source | Albedo(energy_edges=e_edges) - >>> observed(e_centers) - + .. plot:: + :include-source: + + import astropy.units as u + import numpy as np + import matplotlib.pyplot as plt + + from astropy.modeling.powerlaws import PowerLaw1D + from astropy.visualization import quantity_support + + from sunkit_spex.models.physical.albedo import Albedo + + e_edges = np.linspace(5, 550, 600) * u.keV + e_centers = e_edges[0:-1] + (0.5 * np.diff(e_edges)) + source = PowerLaw1D(amplitude=1*u.ph/(u.cm*u.s), x_0=5*u.keV, alpha=3) + albedo = Albedo(energy_edges=e_edges) + observed = source | albedo + + with quantity_support(): + plt.figure() + plt.plot(e_centers, source(e_centers), 'k', label='Source') + for i, t in enumerate([0, 45, 90]*u.deg): + albedo.theta = t + plt.plot(e_centers, observed(e_centers), '--', label=f'Observed, theta={t}', color=f'C{i+1}') + plt.plot(e_centers, observed(e_centers) - source(e_centers), ':', + label=f'Reflected, theta={t}', color=f'C{i+1}') + + plt.ylim(1e-6, 1) + plt.xlim(5, 550) + plt.loglog() + plt.legend() + plt.show() + """ + n_inputs = 1 + n_outputs = 1 theta = Parameter( name="theta", default=0, @@ -55,20 +78,17 @@ class Albedo(FittableModel): min=-90, max=90, description="Angle between the observer and the source", - fixed=True, + fixed=False, ) anisotropy = Parameter(default=1, description="The anisotropy used for albedo correction", fixed=True) - n_inputs = 1 - n_outputs = 1 - def __init__(self, *args, **kwargs): self.energy_edges = kwargs.pop("energy_edges") - return super().__init__(*args, **kwargs) + super().__init__(*args, **kwargs) def evaluate(self, spectrum, theta, anisotropy): - albedo_matrix_T = get_albedo_matrix(self.energy_edges, theta, anisotropy) - return spectrum + spectrum @ albedo_matrix_T + albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) + return spectrum + spectrum @ albedo_matrix @lru_cache @@ -124,9 +144,7 @@ def _get_green_matrix(theta: float) -> RegularGridInterpolator: energy_grid_edges = green["p"].edges[0] energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) - green_matrix_interpolator = RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) - - return green_matrix_interpolator + return RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) @lru_cache @@ -157,9 +175,7 @@ def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotrop albedo_interp = (albedo_interp * de) / anisotropy # Take a transpose - albedo_interp_T = albedo_interp.T - - return albedo_interp_T + return albedo_interp.T @u.quantity_input @@ -167,11 +183,8 @@ def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], ani r""" Get albedo correction matrix. - Matrix used to correct a photon spectrum for the component reflected by the solar atmosphere following - [Kontar2006]_ using precomputed green matrices distributed as part of SSW_. - - .. [Kontar2006] https://doi.org/10.1051/0004-6361:20053672 - .. [SSW] https://www.lmsal.com/solarsoft/ + Matrix used to correct a photon spectrum for the component reflected by the solar atmosphere following interpolated + to given angle and energy indices. Parameters ---------- @@ -201,7 +214,4 @@ def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], ani if np.abs(theta) > 90 * u.deg: raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") anisotropy = np.array(anisotropy).squeeze() - albedo_matrix = _calculate_albedo_matrix( - tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item() - ) - return albedo_matrix + return _calculate_albedo_matrix(tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item()) From 9e0137c9f16f2650f2d54c97c4e4e74e4ca097c4 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 27 Dec 2024 13:52:55 +0000 Subject: [PATCH 07/21] Fix test issue --- sunkit_spex/legacy/fitting/fitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 6545d752..60cc09ab 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -2441,7 +2441,7 @@ def _calc_hessian(self, mod_without_x, fparams, step=0.01, _abs_step=None): replacement_steps = np.mean(self._minimize_solution.final_simplex[0][:, zero_steps], axis=0) * step except AttributeError: replacement_steps = 0 - steps[zero_steps] = replacement_steps if replacement_steps > 0 else step + steps[zero_steps] = replacement_steps if replacement_steps.size > 0 else step else: steps = _abs_step From b15fe408826ffcdd6aaf2170ca29c2f770708572 Mon Sep 17 00:00:00 2001 From: Jake Mitchell Date: Fri, 10 Jan 2025 10:20:35 +0100 Subject: [PATCH 08/21] Unit handling now works with Scipy fitting A very minor edit to the evaluate step in the Albedo class which allows unit handling to work with Scipy fitting. This was previously not working at the point of calling get_albedo_matrix, simply adding in self.theta here solves the problem. --- sunkit_spex/models/physical/albedo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index f60d68e3..ac866097 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -87,7 +87,7 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) def evaluate(self, spectrum, theta, anisotropy): - albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) + albedo_matrix = get_albedo_matrix(self.energy_edges, self.theta, anisotropy) return spectrum + spectrum @ albedo_matrix From 7cd2333a93bc9d0009bcd437bf7712b8d8d64c7f Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Mon, 13 Jan 2025 16:40:03 +0100 Subject: [PATCH 09/21] Added in unit handling for Astropy The model can now be applied to a spectrum which can be fit by the astropy fitting routines with units handled correctly. However there is still an issue with fitting theta as a free parameter, practically this isnt an issue, however should in principle work. --- sunkit_spex/models/physical/albedo.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index ac866097..b16c255c 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -82,14 +82,22 @@ class Albedo(FittableModel): ) anisotropy = Parameter(default=1, description="The anisotropy used for albedo correction", fixed=True) + _input_units_allow_dimensionless = True + def __init__(self, *args, **kwargs): self.energy_edges = kwargs.pop("energy_edges") super().__init__(*args, **kwargs) def evaluate(self, spectrum, theta, anisotropy): - albedo_matrix = get_albedo_matrix(self.energy_edges, self.theta, anisotropy) + if hasattr(theta, "unit"): + albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) + else: + albedo_matrix = get_albedo_matrix(self.energy_edges, theta*u.deg, anisotropy) return spectrum + spectrum @ albedo_matrix + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return {"theta": u.deg} + @lru_cache def _get_green_matrix(theta: float) -> RegularGridInterpolator: @@ -179,7 +187,7 @@ def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotrop @u.quantity_input -def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], anisotropy=1): +def get_albedo_matrix(energy_edges: Quantity[u.keV], theta:Quantity[u.deg], anisotropy=1): r""" Get albedo correction matrix. @@ -211,6 +219,7 @@ def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], ani if energy_edges[0].to_value(u.keV) < 3 or energy_edges[-1].to_value(u.keV) > 600: raise ValueError("Supported energy range 3 <= E <= 600 keV") theta = np.array(theta).squeeze() << theta.unit + # theta = np.array(theta)*u.deg if np.abs(theta) > 90 * u.deg: raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") anisotropy = np.array(anisotropy).squeeze() From bce28e6e1cc2723ffa13b29e4dcb94f02f86d307 Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Tue, 14 Jan 2025 15:19:43 +0100 Subject: [PATCH 10/21] Astropy fitting is now fully functional The code has been restrcutred such that Astropy fitting is now fully functional and all required funtionality is contained within the class rather than externally. --- sunkit_spex/models/physical/albedo.py | 218 +++++++++++--------------- 1 file changed, 91 insertions(+), 127 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index b16c255c..3b31dfd3 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -11,7 +11,7 @@ from sunpy.data import cache -__all__ = ["Albedo", "get_albedo_matrix"] +__all__ = ["Albedo"] class Albedo(FittableModel): @@ -84,143 +84,107 @@ class Albedo(FittableModel): _input_units_allow_dimensionless = True - def __init__(self, *args, **kwargs): - self.energy_edges = kwargs.pop("energy_edges") - super().__init__(*args, **kwargs) - - def evaluate(self, spectrum, theta, anisotropy): - if hasattr(theta, "unit"): - albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) - else: - albedo_matrix = get_albedo_matrix(self.energy_edges, theta*u.deg, anisotropy) - return spectrum + spectrum @ albedo_matrix - - def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"theta": u.deg} - + def __init__(self, theta=u.Quantity(theta.default, theta.unit), + anisotropy=anisotropy.default, + **kwargs): -@lru_cache -def _get_green_matrix(theta: float) -> RegularGridInterpolator: - r""" - Get greens matrix for given angle. + self.energy_edges = kwargs.pop("energy_edges") - Interpolates pre-computed green matrices for fixed angles. The resulting greens matrix is then loaded into an - interpolator for later energy interpolation. + self._get_green_matrix(theta) - Parameters - ========== - theta : float - Angle between the observer and the source + super().__init__(theta=theta, + anisotropy=anisotropy, + **kwargs) - Returns - ======= - Greens matrix interpolator - """ - mu = np.cos(theta) - - base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" - # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter - # load precomputed green matrices - if 0.5 <= mu <= 0.95: - low = 5 * np.floor(mu * 20) - high = 5 * np.ceil(mu * 20) - low_name = f"green_compton_mu{low:03.0f}.dat" - high_name = f"green_compton_mu{high:03.0f}.dat" - low_file = cache.download(base_url + low_name) - high_file = cache.download(base_url + high_name) - green = readsav(low_file) - albedo_low = green["p"].albedo[0] - green_high = readsav(high_file) - albedo_high = green_high["p"].albedo[0] - # why 20? - albedo = albedo_low + (albedo_high - albedo_low) * (mu - (np.floor(mu * 20)) / 20) - - elif mu < 0.5: - file = "green_compton_mu005.dat" - file = cache.download(base_url + file) - green = readsav(file) - albedo = green["p"].albedo[0] - elif mu > 0.95: - file = "green_compton_mu095.dat" - file = cache.download(base_url + file) - green = readsav(file) - albedo = green["p"].albedo[0] - - albedo = albedo.T - - # By construction in keV - energy_grid_edges = green["p"].edges[0] - energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) - - return RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) - - -@lru_cache -def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotropy: float) -> NDArray: - r""" - Calculate green matrix for given energies and angle. - - Interpolates precomputed green matrices for given energies and angle. + def evaluate(self, spectrum, theta, anisotropy): - Parameters - ========== - energy_edges : - Energy edges associated with the spectrum - theta : - Angle between the observer and the source - anisotropy : - Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source - """ - albedo_interpolator = _get_green_matrix(theta) - de = np.diff(energy_edges) - energy_centers = energy_edges[:-1] + de / 2 + if hasattr(theta, "unit"): + theta = Quantity(theta.value,theta.unit) + else: + theta = theta*u.deg - X, Y = np.meshgrid(energy_centers, energy_centers) + if self.energy_edges[0].to_value(u.keV) < 3 or self.energy_edges[-1].to_value(u.keV) > 600: + raise ValueError("Supported energy range 3 <= E <= 600 keV") + theta = np.array(theta).squeeze() << theta.unit + if np.abs(theta) > 90 * u.deg: + raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") + anisotropy = np.array(anisotropy).squeeze() + + energy_edges = tuple(self.energy_edges.to_value(u.keV)) + theta = theta.to_value(u.deg) + anisotropy = anisotropy.item() - albedo_interp = albedo_interpolator((X, Y)) + albedo_interpolator = self._get_green_matrix(theta) + de = np.diff(energy_edges) + energy_centers = energy_edges[:-1] + de / 2 - # Scale by anisotropy - albedo_interp = (albedo_interp * de) / anisotropy + X, Y = np.meshgrid(energy_centers, energy_centers) - # Take a transpose - return albedo_interp.T + albedo_interp = albedo_interpolator((X, Y)) + # Scale by anisotropy + albedo_interp = (albedo_interp * de) / anisotropy -@u.quantity_input -def get_albedo_matrix(energy_edges: Quantity[u.keV], theta:Quantity[u.deg], anisotropy=1): - r""" - Get albedo correction matrix. + albedo_matrix = albedo_interp.T - Matrix used to correct a photon spectrum for the component reflected by the solar atmosphere following interpolated - to given angle and energy indices. + return spectrum + spectrum @ albedo_matrix + + + @lru_cache + def _get_green_matrix(self, theta: float) -> RegularGridInterpolator: + r""" + Get greens matrix for given angle. + + Interpolates pre-computed green matrices for fixed angles. The resulting greens matrix is then loaded into an + interpolator for later energy interpolation. + + Parameters + ========== + theta : float + Angle between the observer and the source + + Returns + ======= + Greens matrix interpolator + """ + mu = np.cos(theta) + + base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" + # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter + # load precomputed green matrices + if 0.5 <= mu <= 0.95: + low = 5 * np.floor(mu * 20) + high = 5 * np.ceil(mu * 20) + low_name = f"green_compton_mu{low:03.0f}.dat" + high_name = f"green_compton_mu{high:03.0f}.dat" + low_file = cache.download(base_url + low_name) + high_file = cache.download(base_url + high_name) + green = readsav(low_file) + albedo_low = green["p"].albedo[0] + green_high = readsav(high_file) + albedo_high = green_high["p"].albedo[0] + # why 20? + albedo = albedo_low + (albedo_high - albedo_low) * (mu - (np.floor(mu * 20)) / 20) + + elif mu < 0.5: + file = "green_compton_mu005.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + elif mu > 0.95: + file = "green_compton_mu095.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + + albedo = albedo.T + + # By construction in keV + energy_grid_edges = green["p"].edges[0] + energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) + + return RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) - Parameters - ---------- - energy_edges : - Energy edges associated with the spectrum - theta : - Angle between Sun-observer line and X-ray source - anisotropy : - Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return {"theta": u.deg} - Example - ------- - >>> import astropy.units as u - >>> import numpy as np - >>> from sunkit_spex.models.physical.albedo import get_albedo_matrix - >>> e = np.linspace(5, 500, 5)*u.keV - >>> albedo_matrix = get_albedo_matrix(e,theta=45*u.deg) - >>> albedo_matrix - array([[3.80274484e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], - [5.10487362e-01, 8.35309813e-07, 0.00000000e+00, 0.00000000e+00], - [3.61059918e-01, 2.48711099e-01, 2.50744411e-09, 0.00000000e+00], - [3.09323903e-01, 2.66485260e-01, 1.23563372e-01, 1.81846722e-10]]) - """ - if energy_edges[0].to_value(u.keV) < 3 or energy_edges[-1].to_value(u.keV) > 600: - raise ValueError("Supported energy range 3 <= E <= 600 keV") - theta = np.array(theta).squeeze() << theta.unit - # theta = np.array(theta)*u.deg - if np.abs(theta) > 90 * u.deg: - raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") - anisotropy = np.array(anisotropy).squeeze() - return _calculate_albedo_matrix(tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item()) From da1407c0dd67bcadd7d1437f97009749c81ce716 Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Mon, 20 Jan 2025 15:11:36 +0100 Subject: [PATCH 11/21] Functions moved out of class, fitting works. The functions have been moved out of the class which now functions as a wrapper. The lru_cache method is applied to both the _calculate_albedo_matrix and _get_green_matrix functions, the outputs when fitting are identical with or without the decorators, and the code runs much quicker this way. --- sunkit_spex/models/physical/albedo.py | 212 ++++++++++++++++---------- 1 file changed, 129 insertions(+), 83 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 3b31dfd3..87993acf 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -11,7 +11,7 @@ from sunpy.data import cache -__all__ = ["Albedo"] +__all__ = ["Albedo","get_albedo_matrix"] class Albedo(FittableModel): @@ -84,107 +84,153 @@ class Albedo(FittableModel): _input_units_allow_dimensionless = True - def __init__(self, theta=u.Quantity(theta.default, theta.unit), - anisotropy=anisotropy.default, + def __init__(self, *args, **kwargs): self.energy_edges = kwargs.pop("energy_edges") - self._get_green_matrix(theta) - - super().__init__(theta=theta, - anisotropy=anisotropy, + super().__init__(*args, **kwargs) + + def evaluate(self, spectrum, theta, anisotropy): if hasattr(theta, "unit"): theta = Quantity(theta.value,theta.unit) else: theta = theta*u.deg + + albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) - if self.energy_edges[0].to_value(u.keV) < 3 or self.energy_edges[-1].to_value(u.keV) > 600: - raise ValueError("Supported energy range 3 <= E <= 600 keV") - theta = np.array(theta).squeeze() << theta.unit - if np.abs(theta) > 90 * u.deg: - raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") - anisotropy = np.array(anisotropy).squeeze() - - energy_edges = tuple(self.energy_edges.to_value(u.keV)) - theta = theta.to_value(u.deg) - anisotropy = anisotropy.item() + return spectrum + spectrum @ albedo_matrix - albedo_interpolator = self._get_green_matrix(theta) - de = np.diff(energy_edges) - energy_centers = energy_edges[:-1] + de / 2 + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return {"theta": u.deg} - X, Y = np.meshgrid(energy_centers, energy_centers) - albedo_interp = albedo_interpolator((X, Y)) - # Scale by anisotropy - albedo_interp = (albedo_interp * de) / anisotropy +@lru_cache +def _get_green_matrix(theta: float) -> RegularGridInterpolator: + r""" + Get greens matrix for given angle. - albedo_matrix = albedo_interp.T + Interpolates pre-computed green matrices for fixed angles. The resulting greens matrix is then loaded into an + interpolator for later energy interpolation. - return spectrum + spectrum @ albedo_matrix - + Parameters + ========== + theta : float + Angle between the observer and the source - @lru_cache - def _get_green_matrix(self, theta: float) -> RegularGridInterpolator: - r""" - Get greens matrix for given angle. - - Interpolates pre-computed green matrices for fixed angles. The resulting greens matrix is then loaded into an - interpolator for later energy interpolation. - - Parameters - ========== - theta : float - Angle between the observer and the source - - Returns - ======= - Greens matrix interpolator - """ - mu = np.cos(theta) - - base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" - # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter - # load precomputed green matrices - if 0.5 <= mu <= 0.95: - low = 5 * np.floor(mu * 20) - high = 5 * np.ceil(mu * 20) - low_name = f"green_compton_mu{low:03.0f}.dat" - high_name = f"green_compton_mu{high:03.0f}.dat" - low_file = cache.download(base_url + low_name) - high_file = cache.download(base_url + high_name) - green = readsav(low_file) - albedo_low = green["p"].albedo[0] - green_high = readsav(high_file) - albedo_high = green_high["p"].albedo[0] - # why 20? - albedo = albedo_low + (albedo_high - albedo_low) * (mu - (np.floor(mu * 20)) / 20) - - elif mu < 0.5: - file = "green_compton_mu005.dat" - file = cache.download(base_url + file) - green = readsav(file) - albedo = green["p"].albedo[0] - elif mu > 0.95: - file = "green_compton_mu095.dat" - file = cache.download(base_url + file) - green = readsav(file) - albedo = green["p"].albedo[0] - - albedo = albedo.T - - # By construction in keV - energy_grid_edges = green["p"].edges[0] - energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) - - return RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) + Returns + ======= + Greens matrix interpolator + """ + mu = np.cos(theta) + + base_url = "https://soho.nascom.nasa.gov/solarsoft/packages/xray/dbase/albedo/" + # what about 0 and 1 assume so close to 05 and 95 that it doesn't matter + # load precomputed green matrices + if 0.5 <= mu <= 0.95: + low = 5 * np.floor(mu * 20) + high = 5 * np.ceil(mu * 20) + low_name = f"green_compton_mu{low:03.0f}.dat" + high_name = f"green_compton_mu{high:03.0f}.dat" + low_file = cache.download(base_url + low_name) + high_file = cache.download(base_url + high_name) + green = readsav(low_file) + albedo_low = green["p"].albedo[0] + green_high = readsav(high_file) + albedo_high = green_high["p"].albedo[0] + # why 20? + albedo = albedo_low + (albedo_high - albedo_low) * (mu - (np.floor(mu * 20)) / 20) + + elif mu < 0.5: + file = "green_compton_mu005.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + elif mu > 0.95: + file = "green_compton_mu095.dat" + file = cache.download(base_url + file) + green = readsav(file) + albedo = green["p"].albedo[0] + + albedo = albedo.T + + # By construction in keV + energy_grid_edges = green["p"].edges[0] + energy_grid_centers = energy_grid_edges[:, 0] + (np.diff(energy_grid_edges, axis=1) / 2).reshape(-1) + + return RegularGridInterpolator((energy_grid_centers, energy_grid_centers), albedo) + + +@lru_cache +def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotropy: float) -> NDArray: + r""" + Calculate green matrix for given energies and angle. - def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): - return {"theta": u.deg} + Interpolates precomputed green matrices for given energies and angle. + Parameters + ========== + energy_edges : + Energy edges associated with the spectrum + theta : + Angle between the observer and the source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source + """ + albedo_interpolator = _get_green_matrix(theta) + de = np.diff(energy_edges) + energy_centers = energy_edges[:-1] + de / 2 + + X, Y = np.meshgrid(energy_centers, energy_centers) + + albedo_interp = albedo_interpolator((X, Y)) + + # Scale by anisotropy + albedo_interp = (albedo_interp * de) / anisotropy + + # Take a transpose + return albedo_interp.T + + +def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], anisotropy=1): + r""" + Get albedo correction matrix. + + Matrix used to correct a photon spectrum for the component reflected by the solar atmosphere following interpolated + to given angle and energy indices. + + Parameters + ---------- + energy_edges : + Energy edges associated with the spectrum + theta : + Angle between Sun-observer line and X-ray source + anisotropy : + Ratio of the flux in observer direction to the flux downwards, 1 for an isotropic source + + Example + ------- + >>> import astropy.units as u + >>> import numpy as np + >>> from sunkit_spex.models.physical.albedo import get_albedo_matrix + >>> e = np.linspace(5, 500, 5)*u.keV + >>> albedo_matrix = get_albedo_matrix(e,theta=45*u.deg) + >>> albedo_matrix + array([[3.80274484e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], + [5.10487362e-01, 8.35309813e-07, 0.00000000e+00, 0.00000000e+00], + [3.61059918e-01, 2.48711099e-01, 2.50744411e-09, 0.00000000e+00], + [3.09323903e-01, 2.66485260e-01, 1.23563372e-01, 1.81846722e-10]]) + """ + if energy_edges[0].to_value(u.keV) < 3 or energy_edges[-1].to_value(u.keV) > 600: + raise ValueError("Supported energy range 3 <= E <= 600 keV") + theta = np.array(theta).squeeze() << theta.unit + if np.abs(theta) > 90 * u.deg: + raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") + anisotropy = np.array(anisotropy).squeeze() + + return _calculate_albedo_matrix(tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item()) \ No newline at end of file From 725001c56072f779e0f3ef4a2388d2966e22ebb3 Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Tue, 21 Jan 2025 16:23:23 +0100 Subject: [PATCH 12/21] Update sunkit_spex/models/physical/albedo.py Fixes errors in the test process. Co-authored-by: Shane Maloney --- sunkit_spex/models/physical/albedo.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 87993acf..d6ba884c 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -197,6 +197,7 @@ def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotrop return albedo_interp.T +@quantity_input def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], anisotropy=1): r""" Get albedo correction matrix. From ac154bc05dee450d103aec2fd12cbf467351eac6 Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Wed, 22 Jan 2025 10:59:46 +0100 Subject: [PATCH 13/21] Update sunkit_spex/models/physical/albedo.py Co-authored-by: Shane Maloney --- sunkit_spex/models/physical/albedo.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index d6ba884c..db418265 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -96,9 +96,7 @@ def __init__(self, *args, def evaluate(self, spectrum, theta, anisotropy): - if hasattr(theta, "unit"): - theta = Quantity(theta.value,theta.unit) - else: + if not isinstance(theta, Quantity) theta = theta*u.deg albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) From b9299a89beedae4d77a7f51981581e362df6abbf Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 22 Jan 2025 10:05:17 +0000 Subject: [PATCH 14/21] Update sunkit_spex/models/physical/albedo.py --- sunkit_spex/models/physical/albedo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index db418265..376bfd60 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -96,7 +96,7 @@ def __init__(self, *args, def evaluate(self, spectrum, theta, anisotropy): - if not isinstance(theta, Quantity) + if not isinstance(theta, Quantity): theta = theta*u.deg albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) From 737a7ff531f18c37dee5106861cd6d67957af472 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 22 Jan 2025 10:09:02 +0000 Subject: [PATCH 15/21] Update sunkit_spex/models/physical/albedo.py --- sunkit_spex/models/physical/albedo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 376bfd60..fe586294 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -195,7 +195,7 @@ def _calculate_albedo_matrix(energy_edges: tuple[float], theta: float, anisotrop return albedo_interp.T -@quantity_input +@u.quantity_input def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], anisotropy=1): r""" Get albedo correction matrix. From a5ef1e23d979b4c77b0a69fb49f300cb9699ffad Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 22 Jan 2025 10:11:23 +0000 Subject: [PATCH 16/21] Ruff fixes --- sunkit_spex/models/physical/albedo.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index fe586294..f51b865b 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -11,7 +11,7 @@ from sunpy.data import cache -__all__ = ["Albedo","get_albedo_matrix"] +__all__ = ["Albedo", "get_albedo_matrix"] class Albedo(FittableModel): @@ -84,21 +84,15 @@ class Albedo(FittableModel): _input_units_allow_dimensionless = True - def __init__(self, *args, - **kwargs): - + def __init__(self, *args, **kwargs): self.energy_edges = kwargs.pop("energy_edges") - super().__init__(*args, - **kwargs) - - + super().__init__(*args, **kwargs) def evaluate(self, spectrum, theta, anisotropy): - if not isinstance(theta, Quantity): - theta = theta*u.deg - + theta = theta * u.deg + albedo_matrix = get_albedo_matrix(self.energy_edges, theta, anisotropy) return spectrum + spectrum @ albedo_matrix @@ -107,7 +101,6 @@ def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): return {"theta": u.deg} - @lru_cache def _get_green_matrix(theta: float) -> RegularGridInterpolator: r""" @@ -232,4 +225,4 @@ def get_albedo_matrix(energy_edges: Quantity[u.keV], theta: Quantity[u.deg], ani raise ValueError(f"Theta must be between -90 and 90 degrees: {theta}.") anisotropy = np.array(anisotropy).squeeze() - return _calculate_albedo_matrix(tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item()) \ No newline at end of file + return _calculate_albedo_matrix(tuple(energy_edges.to_value(u.keV)), theta.to_value(u.deg), anisotropy.item()) From fcdbc09c4be83870ad4da1f2f3e839e08bbee67e Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Wed, 22 Jan 2025 11:56:12 +0000 Subject: [PATCH 17/21] Update spectrum example for NDCube docstring change --- sunkit_spex/spectrum/spectrum.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index 6bab8128..6c958d1d 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -158,7 +158,7 @@ class Spectrum(NDCube): Date: Wed, 22 Jan 2025 12:45:12 +0000 Subject: [PATCH 18/21] Fix doctest issue disable on older version of NDCube --- pytest.ini | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pytest.ini b/pytest.ini index c9c1945e..0a6723e8 100644 --- a/pytest.ini +++ b/pytest.ini @@ -14,6 +14,8 @@ norecursedirs = .history sunkit_spex/extern doctest_plus = enabled +doctest_subpackage_requires = + * = ndcube>2.3.0 doctest_optionflags = NORMALIZE_WHITESPACE FLOAT_CMP From a93d551952a251213f53741b0a0713eeb9476427 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Thu, 23 Jan 2025 16:27:17 +0000 Subject: [PATCH 19/21] Only skip doctest on specific class --- pytest.ini | 2 -- sunkit_spex/spectrum/spectrum.py | 3 +++ 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pytest.ini b/pytest.ini index 0a6723e8..c9c1945e 100644 --- a/pytest.ini +++ b/pytest.ini @@ -14,8 +14,6 @@ norecursedirs = .history sunkit_spex/extern doctest_plus = enabled -doctest_subpackage_requires = - * = ndcube>2.3.0 doctest_optionflags = NORMALIZE_WHITESPACE FLOAT_CMP diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index 6c958d1d..9f322732 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -11,6 +11,9 @@ __all__ = ["gwcs_from_array", "SpectralAxis", "Spectrum"] +__doctest_requires__ = {"Spectrum": ["ndcube>=2.3"]} + + def gwcs_from_array(array): """ Create a new WCS from provided tabular data. This defaults to being From 2c4069ecef15cc5783183689fa7e75df62c17f86 Mon Sep 17 00:00:00 2001 From: jajmitchell Date: Thu, 30 Jan 2025 16:36:41 +0100 Subject: [PATCH 20/21] Update sunkit_spex/models/physical/albedo.py Co-authored-by: Shane Maloney --- sunkit_spex/models/physical/albedo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index f51b865b..99803370 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -80,7 +80,7 @@ class Albedo(FittableModel): description="Angle between the observer and the source", fixed=False, ) - anisotropy = Parameter(default=1, description="The anisotropy used for albedo correction", fixed=True) + anisotropy = Parameter(name="anisotropy", default=1, description="The anisotropy used for albedo correction", fixed=True) _input_units_allow_dimensionless = True From ea984152cc3e19877d8c6de50b4680296c8ebd9b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 30 Jan 2025 15:44:30 +0000 Subject: [PATCH 21/21] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- sunkit_spex/models/physical/albedo.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/models/physical/albedo.py b/sunkit_spex/models/physical/albedo.py index 99803370..faca26e6 100644 --- a/sunkit_spex/models/physical/albedo.py +++ b/sunkit_spex/models/physical/albedo.py @@ -80,7 +80,9 @@ class Albedo(FittableModel): description="Angle between the observer and the source", fixed=False, ) - anisotropy = Parameter(name="anisotropy", default=1, description="The anisotropy used for albedo correction", fixed=True) + anisotropy = Parameter( + name="anisotropy", default=1, description="The anisotropy used for albedo correction", fixed=True + ) _input_units_allow_dimensionless = True