Skip to content

Commit

Permalink
make RhessiLoader compatible with legacy fitter args (#181)
Browse files Browse the repository at this point in the history
* make RhessiLoader compatible with legacy fitter args

* RHESSI lightcurves: display proper energy ranges

* add time legacy properties to RHESSI loader

* add back old bg subtract property

* changelog

* spectrum: ndcube print change (for docs to pass)

* remove select_times method
  • Loading branch information
settwi authored Feb 1, 2025
1 parent 2125410 commit 44a7611
Show file tree
Hide file tree
Showing 4 changed files with 85 additions and 56 deletions.
1 change: 1 addition & 0 deletions changelog/181.bugfix.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Add back legacy parameters/functionality relating to `~sunkit_spex.extern.stix.RhessiLoader`.
14 changes: 0 additions & 14 deletions examples/fitting_RHESSI_spectra.py
Original file line number Diff line number Diff line change
Expand Up @@ -116,20 +116,6 @@
plt.show()
plt.rcParams["font.size"] = default_text

#####################################################
#
# An alternative to set the event and background times like above is by using the `select_time` method. E.g.,
#
# .. code-block::python
#
# ress_spec.loaded_spec_data["spectrum1"].select_time(start="2002-10-05T10:38:32", end="2002-10-05T10:40:32", background=True)
#
# ress_spec.loaded_spec_data["spectrum1"].select_time(start="2002-10-05T10:41:20", end="2002-10-05T10:42:24")
#
# Both and end and a start time needs to be defined for the background, whereas the event time is assumed to
# commence/finish at the first/final data time if the start/end time is not given.
#

#####################################################
#
# Now let's get going with a model and explicitly stating a fit statistic
Expand Down
124 changes: 83 additions & 41 deletions sunkit_spex/extern/rhessi.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,15 @@ class RhessiLoader(instruments.InstrumentBlueprint):
[1] https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html
"""

def __init__(self, spectrum_fn, srm_fn, **kwargs):
def __init__(self, spectrum_fn, srm_fn):
"""
Spectrum and SRM files are both required: attenuator state change times
are in the spectrum file,
and the state determines which SRM will be used.
"""
self._construction_string = f"RhessiLoader(spectrum_fn={spectrum_fn}, " f"srm_fn={srm_fn}," f"**{kwargs})"
self._construction_string = f"RhessiLoader(spectrum_fn={spectrum_fn}, " f"srm_fn={srm_fn})"
self._systematic_error = 0
self.load_prepare_spectrum_srm(spectrum_fn, srm_fn)
self._start_background_time, self._end_background_time = None, None

@property
def systematic_error(self):
Expand Down Expand Up @@ -96,9 +95,11 @@ def load_prepare_spectrum_srm(self, spectrum_fn, srm_fn):
}
self._update_event_data_with_times()
self._original_data = copy.deepcopy(self._loaded_spec_data)
# Set these times to sensible defaults
self.update_background_times(self._start_event_time, self._end_event_time)

@property
def subtract_background(self):
def subtract_background(self) -> bool:
"""
States whether the the data is event minus background or not.
"""
Expand Down Expand Up @@ -411,10 +412,7 @@ def _rebin_ts(self, times, clump_bins):
def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):
"""Creates a RHESSI lightcurve.
Helps the user see the RHESSI time profile. The defined event time (defined either through `start_event_time`,
`end_event_time` setters, or `select_time(...)` method) is shown with a purple shaded region and if a background
time (defined either through `start_background_time`, `end_background_time` setters, or
`select_time(...,background=True)` method) is defined then it is shown with an orange shaded region.
The background and event times are highlighted based on what the user has set.
Parameters
----------
Expand Down Expand Up @@ -445,10 +443,9 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):
elif len(np.shape(energy_ranges)) == 1:
energy_ranges = [energy_ranges]
elif len(np.shape(energy_ranges)) == 0:
print(
raise ValueError(
"The `energy_ranges` input should be a range of two energy values (e.g., [4,8] meaning 4-8 keV inclusive) or a list of these ranges."
)
return

ax = axes or plt.gca()
font_size = plt.rcParams["font.size"]
Expand All @@ -462,11 +459,11 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):

# plot each energy range
time_binning = (flat_times[1:] - flat_times[:-1]).to_value("s")
left_edges, right_edges = self._loaded_spec_data["count_channel_bins"].T
for er in energy_ranges:
i = np.where(
(self._loaded_spec_data["count_channel_bins"][:, 0] >= er[0])
& (self._loaded_spec_data["count_channel_bins"][:, -1] <= er[-1])
)
# Restrict the energy range to the data energy bins
ea, eb = max(er[0], left_edges[0]), min(er[1], right_edges[-1])
i = np.where((left_edges >= ea) & (right_edges <= eb))
e_range_cts = np.sum(self._spectrum["counts"][:, i].reshape((len(time_binning), -1)), axis=1)

if rebin_time > 1:
Expand All @@ -475,7 +472,7 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):

e_range_ctr, e_range_ctr_err = e_range_cts / time_binning, np.sqrt(e_range_cts) / time_binning
lc = e_range_ctr
p = ax.stairs(lc, flat_times.datetime, label=f"{er[0]}$-${er[-1]} keV")
p = ax.stairs(lc, flat_times.datetime, label=f"{ea}$-${eb} keV")
ax.errorbar(
x=time_mids.datetime, y=e_range_ctr, yerr=e_range_ctr_err, c=p.get_edgecolor(), ls=""
) # error bar in middle of the bin
Expand Down Expand Up @@ -515,10 +512,7 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1):
def spectrogram(self, axes=None, rebin_time=1, rebin_energy=1, **kwargs):
"""Creates a RHESSI spectrogram.
Helps the user see the RHESSI time and energy evolution. The defined event time (defined either through
`start_event_time`, `end_event_time` setters, or `select_time(...)` method) is shown with a violet line
and if a background time (defined either through `start_background_time`, `end_background_time` setters,
or `select_time(...,background=True)` method) is defined then it is shown with an orange line.
The background and event times are highlighted based on what the user has set.
Parameters
----------
Expand Down Expand Up @@ -604,32 +598,65 @@ def spectrogram(self, axes=None, rebin_time=1, rebin_energy=1, **kwargs):

return ax

def select_time(self, start=None, end=None, background=False):
"""Provides method to set start and end time of the event or background in one line.
# "retroactive" properties for setting the times like in older versions,
# for legacy compatibility.
@property
def start_event_time(self):
time_depr_warn('event')
return self._start_event_time

Parameters
----------
start, end : str, `astropy.Time`, None
String to be given to astropy's Time, `astropy.Time` is used directly, None sets the
start/end event time to be the first time of the data. None doesn't add, or will remove,
any background data in `_loaded_spec_data["extras"]` if background is set to True.
Default: None
@start_event_time.setter
def start_event_time(self, new_time: atime.Time):
time_depr_warn('event')
self.update_event_times(new_time, self._end_event_time)

background : bool
Determines whether the start and end times are for the event (False) or background
time (True).
Default: False
@property
def end_event_time(self):
time_depr_warn('event')
return self._end_event_time

Returns
-------
None.
"""
self.__warn = False if (start is not None) and (end is not None) else True
@end_event_time.setter
def end_event_time(self, new_time: atime.Time):
time_depr_warn('event')
self.update_event_times(self._start_event_time, new_time)

if background:
self.start_background_time, self.end_background_time = start, end
else:
self.start_event_time, self.end_event_time = start, end
@property
def start_background_time(self):
time_depr_warn('background')
return self._start_background_time

@start_background_time.setter
def start_background_time(self, new_time: atime.Time):
time_depr_warn('background')
self.update_background_times(new_time, self._end_background_time)

@property
def end_background_time(self):
time_depr_warn('background')
return self._end_background_time

@end_background_time.setter
def end_background_time(self, new_time: atime.Time):
time_depr_warn('background')
self.update_background_times(self._start_background_time, new_time)

# "retroactive" properties for specifying background subtraction,
# like in older versions.
@property
def data2data_minus_background(self) -> bool:
warnings.warn(
"Please use the `subtract_background` property instead.",
stacklevel=2
)
return self.subtract_background

@data2data_minus_background.setter
def data2data_minus_background(self, do_subtract: bool):
warnings.warn(
"Please use the `subtract_background` property instead.",
stacklevel=2
)
self.subtract_background = do_subtract


def load_spectrum(spec_fn: str):
Expand Down Expand Up @@ -801,3 +828,18 @@ def load_srm(srm_file: str):
}

return dict(channel_bins=channel_bins, photon_bins=photon_bins, srm_options=ret_srms)


class LegacyRhessiLoader(RhessiLoader):
'''Allow legacy `Fitter` to load in RHESSI data with the expected args format'''
def __init__(self, spec_fn, **kw):
srm_fn = kw.pop('srm_file')
super().__init__(spec_fn, srm_fn)


def time_depr_warn(partial_name: str):
warnings.warn(
"This time-related property is deprecated. "
f"In the future, please use `RhessiLoader.update_{partial_name}_times`",
stacklevel=2
)
2 changes: 1 addition & 1 deletion sunkit_spex/legacy/fitting/data_loader.py
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ def __init__(
f"srm_file={srm_file},srm_custom={srm_custom},custom_channel_bins={custom_channel_bins}, **{kwargs})"
)

self.instrument_loaders = {"NuSTAR": inst.NustarLoader, "RHESSI": rhessi.RhessiLoader, "Solar Orbiter": stix.STIXLoader}
self.instrument_loaders = {"NuSTAR": inst.NustarLoader, "RHESSI": rhessi.LegacyRhessiLoader, "Solar Orbiter": stix.STIXLoader}

pha_file, arf_file, rmf_file, srm_file, srm_custom, custom_channel_bins, instruments = self._sort_files(
pha_file=pha_file,
Expand Down

0 comments on commit 44a7611

Please sign in to comment.