From 49e8a25b52327fa6abacf194ae2d97d4f45b96ca Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Tue, 14 Jul 2026 22:58:06 +0100 Subject: [PATCH 01/11] fix to albedo bumps Clunky way of doing it but means make_model input is still photons/s/cm^2 and internal calc only changes if also provide photon_energy_widths. --- sunkit_spex/legacy/fitting/fitter.py | 41 ++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 1343b467..d51e9373 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -1474,14 +1474,21 @@ def _counts_model(self, **kwargs): # assign the spectrum parameter values to their model param counterpart sep_params = dict(zip(self._orig_params, ordered_kwarg_values)) - # calculate the [photon s^-1 cm^-2] + # calculate the [photon s^-1 cm^-2] -> keep this as before m = ( self._model(**sep_params, energies=kwargs["photon_channels"][s]) * kwargs["photon_channel_widths"][s] ) # np.diff(kwargs["photon_channels"][s]).flatten() # remove energy bin dependence # fold the photon model through the SRM to create the count rate model, [photon s^-1 cm^-2] * [count photon^-1 cm^2] = [count s^-1] cts_model, cts_albedo_exess = make_model( - energies=kwargs["photon_channels"][s], photon_model=m, parameters=None, srm=kwargs["total_responses"][s], albedo_corr=self.albedo_corr, albedo_angle =self.albedo_angle, albedo_anisotropy=self.albedo_anisotropy + energies=kwargs["photon_channels"][s], + photon_model=m, + parameters=None, + photon_energy_widths=kwargs["photon_channel_widths"][s], + srm=kwargs["total_responses"][s], + albedo_corr=self.albedo_corr, + albedo_angle=self.albedo_angle, + albedo_anisotropy=self.albedo_anisotropy, ) if "scaled_background_spectrum" + str(s + 1) in self._scaled_backgrounds: @@ -2603,8 +2610,9 @@ def _calc_counts_model(self, photon_model, parameters=None, spectrum="spectrum1" if for_plotting: cts_model, _ = make_model( energies=photon_channel_bins, - photon_model=m * np.diff(photon_channel_bins).flatten(), + photon_model=m * np.diff(photon_channel_bins).flatten(), # keep this as before [photon s^-1 cm^-2] parameters=None, + photon_energy_widths=np.diff(photon_channel_bins).flatten(), srm=srm, albedo_corr=False, albedo_angle=self.albedo_angle, @@ -2614,8 +2622,9 @@ def _calc_counts_model(self, photon_model, parameters=None, spectrum="spectrum1" else: cts_model, _ = make_model( energies=photon_channel_bins, - photon_model=m * np.diff(photon_channel_bins).flatten(), + photon_model=m * np.diff(photon_channel_bins).flatten(), # keep this as before [photon s^-1 cm^-2] parameters=None, + photon_energy_widths=np.diff(photon_channel_bins).flatten(), srm=srm, albedo_corr=self.albedo_corr, albedo_angle=self.albedo_angle, @@ -5655,7 +5664,16 @@ def imports(): return _imps -def make_model(energies=None, photon_model=None, parameters=None, srm=None, albedo_corr=False, albedo_angle=0, albedo_anisotropy=1): +def make_model( + energies=None, + photon_model=None, + parameters=None, + photon_energy_widths=None, + srm=None, + albedo_corr=False, + albedo_angle=0, + albedo_anisotropy=1, +): """Takes a photon model array ( or function if you provide the pinputs with parameters), the spectral response matrix and returns a model count spectrum. Parameters @@ -5665,13 +5683,18 @@ def make_model(energies=None, photon_model=None, parameters=None, srm=None, albe Default : None photon_model : function/array/list - Array -OR- function representing the photon model (if it's a function, provide the parameters of the function as a list, e.g. parameters = [energies, const, power]). + Array [photon s^-1 cm^-2] -OR- function representing the photon model (if it's a function, provide the parameters of the function as a list, e.g. parameters = [energies, const, power]). Default : None parameters : list List representing the inputs a photon model function, if a function is provided, excludeing the energies the spectrum is over. Default : None + photon_energy_widths : array/list + Optional photon-energy bin widths. If provided, the photon model is divided by these widths, the SRM multiplied by them before + folding through the SRM. + Default : None + srm : matrix/array Spectral response matrix. Default : None @@ -5687,6 +5710,12 @@ def make_model(energies=None, photon_model=None, parameters=None, srm=None, albe else: photon_spec = photon_model(energies, *parameters) + # if photon_energy_widths provided remove from photon spec and apply to the srm instead + # since the photon_spec has been provided + if photon_energy_widths is not None: + photon_spec = photon_spec / photon_energy_widths + srm = photon_energy_widths[:, None] * srm + if albedo_corr: photon_spec, albedo_excess_phot = albedo(photon_spec, energies, albedo_angle, anisotropy=albedo_anisotropy) albedo_excess_count = np.matmul(albedo_excess_phot, srm) From 3937c189eca5d1a04458e7aba23180a6737e5517 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Wed, 15 Jul 2026 15:05:22 +0100 Subject: [PATCH 02/11] Slight sped up to make_model --- sunkit_spex/legacy/fitting/fitter.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index d51e9373..659a4aab 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -5704,25 +5704,30 @@ def make_model( A model count spectrum. """ + def _fold_through_srm(spec, matrix, widths=None): + """Fold photon spectrum through SRM with optional SRM-row width scaling.""" + if widths is None: + return np.matmul(spec, matrix) + # Equivalent to spec @ (widths[:, None] * matrix) without allocating that dense temporary. + return np.einsum("i,ij,i->j", spec, matrix, widths) + # if parameters is None then assume the photon_model input is already a spectrum to test, else make the model spectrum from the function and parameters if type(parameters) == type(None): photon_spec = photon_model else: photon_spec = photon_model(energies, *parameters) - # if photon_energy_widths provided remove from photon spec and apply to the srm instead - # since the photon_spec has been provided + # Albedo correction is applied to differential flux, so divide by bin widths first. if photon_energy_widths is not None: photon_spec = photon_spec / photon_energy_widths - srm = photon_energy_widths[:, None] * srm if albedo_corr: - photon_spec, albedo_excess_phot = albedo(photon_spec, energies, albedo_angle, anisotropy=albedo_anisotropy) - albedo_excess_count = np.matmul(albedo_excess_phot, srm) + photon_spec, albedo_excess_phot = albedo(photon_spec, energies, albedo_angle, anisotropy=albedo_anisotropy) + albedo_excess_count = _fold_through_srm(albedo_excess_phot, srm, photon_energy_widths) else: albedo_excess_count = np.array([]) - model_cts_spectrum = np.matmul(photon_spec, srm) + model_cts_spectrum = _fold_through_srm(photon_spec, srm, photon_energy_widths) return model_cts_spectrum, albedo_excess_count From 744cbb441ae317d2f9674666806fd59c8b90892c Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Wed, 15 Jul 2026 15:43:08 +0100 Subject: [PATCH 03/11] Update fitter.py --- 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 659a4aab..9612f46c 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -5708,7 +5708,7 @@ def _fold_through_srm(spec, matrix, widths=None): """Fold photon spectrum through SRM with optional SRM-row width scaling.""" if widths is None: return np.matmul(spec, matrix) - # Equivalent to spec @ (widths[:, None] * matrix) without allocating that dense temporary. + # Equivalent to spec @ (widths[:, None] * matrix) without allocating full temporary matrix return np.einsum("i,ij,i->j", spec, matrix, widths) # if parameters is None then assume the photon_model input is already a spectrum to test, else make the model spectrum from the function and parameters From 3a9a46ec2b61828652692d0b0456565e5c62690a Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Wed, 15 Jul 2026 22:34:12 +0100 Subject: [PATCH 04/11] tweaking possible speed up --- sunkit_spex/legacy/fitting/fitter.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 9612f46c..5fcc1cd0 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -5709,7 +5709,9 @@ def _fold_through_srm(spec, matrix, widths=None): if widths is None: return np.matmul(spec, matrix) # Equivalent to spec @ (widths[:, None] * matrix) without allocating full temporary matrix - return np.einsum("i,ij,i->j", spec, matrix, widths) + return np.einsum("i,ij,i->j", spec, matrix, widths, optimize=True) + # Above really faster, what about just + # return np.matmul(spec, widths[:, None] * matrix) # if parameters is None then assume the photon_model input is already a spectrum to test, else make the model spectrum from the function and parameters if type(parameters) == type(None): From 74b06b8ad6dc00e6b76b1184ee3f1152b01a0ea8 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Wed, 15 Jul 2026 23:14:39 +0100 Subject: [PATCH 05/11] tweaked albedo plotting --- sunkit_spex/legacy/fitting/fitter.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 5fcc1cd0..d0059e3a 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3855,7 +3855,18 @@ def _plot_1spec( res.plot(energy_channels_res, residuals, color="k", alpha=0.8) # , drawstyle='steps-mid' if self.albedo_corr and albedo_excess_count.size > 0: - axs.plot(energy_channels, albedo_excess_count, color="grey") + axs.plot(energy_channels, albedo_excess_count, color="plum",alpha=1.0) + rainbow_text_lines( + (0.1, 0.99), + strings=["Albedo"], + colors=["plum"], + xycoords="axes fraction", + verticalalignment="top", + horizontalalignment="left", + ax=axes, + alpha=1.0, + fontsize="small", + ) if self._latest_fit_run == "mcmc": _rebin_info = ( @@ -3888,6 +3899,10 @@ def _plot_1spec( ) = count_rates, energy_channel_error, count_rate_errors # other log_plotting_info["residuals"], log_plotting_info["fitting_range"] = residuals, fitting_range + if self.albedo_corr and np.size(albedo_excess_count) > 0: + log_plotting_info["albedo_excess_count"] = np.array(albedo_excess_count).copy() + else: + log_plotting_info["albedo_excess_count"] = np.array([]) # do we have submodels if hasattr(self, "_corresponding_submod_inputs"): log_plotting_info["submodels"] = spec_submods[0] From 5fe360d9ccadac717887f45936087c30f583d52b Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Thu, 16 Jul 2026 00:01:00 +0100 Subject: [PATCH 06/11] improved plotting param format --- sunkit_spex/legacy/fitting/fitter.py | 36 ++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index d0059e3a..12928cb4 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3542,16 +3542,38 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): _xycoords = "axes fraction" for p in plot_params: par_spec = p.split("_spectrum") + pname = par_spec[0] error = self.params["Error", p] + value = self.params["Value", p] + if pname.startswith("EM"): + error = np.asarray(error) * 1e46 + value = value * 1e46 + if pname.startswith("total_eflux"): + error = np.asarray(error) * 1e35 + value = value * 1e35 + value_fmt = ".2e" + error_fmt = ".2e" + if pname.startswith(("T", "index", "e_c")): + value_fmt = ".2f" + error_fmt = ".2f" if np.all(error) == np.all([0, 0]): - param_str += [par_spec[0] + ": {0:.2e}".format(self.params["Value", p]) + "\n"] + param_str += [pname + f": {value:{value_fmt}}" + "\n"] else: - param_str += [ - par_spec[0] - + ": {0:.2e}".format(self.params["Value", p]) - + f"$^{{+{error[1]:.2e}}}_{{-{error[0]:.2e}}}$" - + "\n" - ] # str(round(self.params["Value", p], 2)) + param_str += [pname + f": {value:{value_fmt}}" + + f"$^{{+{error[1]:{error_fmt}}}}_{{-{error[0]:{error_fmt}}}}$" + + "\n" + ] + # par_spec = p.split("_spectrum") + # error = self.params["Error", p] + # if np.all(error) == np.all([0, 0]): + # param_str += [par_spec[0] + ": {0:.2e}".format(self.params["Value", p]) + "\n"] + # else: + # param_str += [ + # par_spec[0] + # + ": {0:.2e}".format(self.params["Value", p]) + # + f"$^{{+{error[1]:.2e}}}_{{-{error[0]:.2e}}}$" + # + "\n" + # ] # str(round(self.params["Value", p], 2)) # join param strings correctly and annotate the axes depending on whether they should be coloured with sub-models if present self._annotate_params(axs, param_str, submod_param_cols, _xycoords) From 629701c77de9b7cf818691b599906405dc4ec6a6 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Thu, 16 Jul 2026 00:10:08 +0100 Subject: [PATCH 07/11] Update fitter.py --- 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 12928cb4..6a55c815 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3553,7 +3553,7 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): value = value * 1e35 value_fmt = ".2e" error_fmt = ".2e" - if pname.startswith(("T", "index", "e_c")): + if pname.startswith(("T", "index", "e_c", "C")): value_fmt = ".2f" error_fmt = ".2f" if np.all(error) == np.all([0, 0]): From 1d02aa9cc6336871365684dc9da7eb5d2340dfe0 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Thu, 16 Jul 2026 11:20:38 +0100 Subject: [PATCH 08/11] further tidy up of plot params --- sunkit_spex/legacy/fitting/fitter.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 6a55c815..28c7d025 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3545,20 +3545,27 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): pname = par_spec[0] error = self.params["Error", p] value = self.params["Value", p] + scl=1 + # Might not have an error value so just change scale factor + # and apply this at the end when making the nice string if pname.startswith("EM"): - error = np.asarray(error) * 1e46 - value = value * 1e46 + scl=1e46 if pname.startswith("total_eflux"): - error = np.asarray(error) * 1e35 - value = value * 1e35 + scl=1e35 + if pname.startswith("plasma_d"): + scl=1e10 value_fmt = ".2e" error_fmt = ".2e" - if pname.startswith(("T", "index", "e_c", "C")): + # Only have f_vth, thick_fn and thick_warm as the models from photon_models_for_fitting.py + if pname.startswith(("T", "index", "e_c", "C", "loop_temp","length")): value_fmt = ".2f" error_fmt = ".2f" if np.all(error) == np.all([0, 0]): + value = value * scl param_str += [pname + f": {value:{value_fmt}}" + "\n"] else: + value = np.asarray(value) * scl + error = np.asarray(error) * scl param_str += [pname + f": {value:{value_fmt}}" + f"$^{{+{error[1]:{error_fmt}}}}_{{-{error[0]:{error_fmt}}}}$" + "\n" From 33c71222c328d88178d4056f02558241b8d9ec4a Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Thu, 16 Jul 2026 11:36:53 +0100 Subject: [PATCH 09/11] added units to plot params --- sunkit_spex/legacy/fitting/fitter.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 28c7d025..1c2aa04d 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3546,14 +3546,24 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): error = self.params["Error", p] value = self.params["Value", p] scl=1 + punts="" # Might not have an error value so just change scale factor # and apply this at the end when making the nice string if pname.startswith("EM"): scl=1e46 + punts=" cm$^{-3}$" if pname.startswith("total_eflux"): scl=1e35 + punts=" e$^-$s$^{-1}$" if pname.startswith("plasma_d"): scl=1e10 + punts=" cm$^{-3}$" + if pname.startswith(("T", "loop_temp")): + punts=" MK" + if pname.startswith("e_c"): + punts=" keV" + if pname.startswith("length"): + punts=" Mm" value_fmt = ".2e" error_fmt = ".2e" # Only have f_vth, thick_fn and thick_warm as the models from photon_models_for_fitting.py @@ -3562,12 +3572,12 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): error_fmt = ".2f" if np.all(error) == np.all([0, 0]): value = value * scl - param_str += [pname + f": {value:{value_fmt}}" + "\n"] + param_str += [pname + f": {value:{value_fmt}}"+ punts + "\n"] else: value = np.asarray(value) * scl error = np.asarray(error) * scl param_str += [pname + f": {value:{value_fmt}}" - + f"$^{{+{error[1]:{error_fmt}}}}_{{-{error[0]:{error_fmt}}}}$" + + f"$^{{+{error[1]:{error_fmt}}}}_{{-{error[0]:{error_fmt}}}}$" + punts + "\n" ] # par_spec = p.split("_spectrum") From 1cd32bc755eb66aadb8f0c319d9e4a3c5aa91fb4 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Thu, 16 Jul 2026 23:53:45 +0100 Subject: [PATCH 10/11] fixed plotting thick warm params --- sunkit_spex/legacy/fitting/fitter.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sunkit_spex/legacy/fitting/fitter.py b/sunkit_spex/legacy/fitting/fitter.py index 1c2aa04d..cea9cac2 100644 --- a/sunkit_spex/legacy/fitting/fitter.py +++ b/sunkit_spex/legacy/fitting/fitter.py @@ -3552,7 +3552,7 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): if pname.startswith("EM"): scl=1e46 punts=" cm$^{-3}$" - if pname.startswith("total_eflux"): + if pname.startswith(("total_eflux", "tot_eflux")): scl=1e35 punts=" e$^-$s$^{-1}$" if pname.startswith("plasma_d"): @@ -3560,14 +3560,14 @@ def _plot_params(self, axs, res, plot_params, submod_param_cols): punts=" cm$^{-3}$" if pname.startswith(("T", "loop_temp")): punts=" MK" - if pname.startswith("e_c"): + if pname.startswith(("e_c", "ec")): punts=" keV" if pname.startswith("length"): punts=" Mm" value_fmt = ".2e" error_fmt = ".2e" # Only have f_vth, thick_fn and thick_warm as the models from photon_models_for_fitting.py - if pname.startswith(("T", "index", "e_c", "C", "loop_temp","length")): + if pname.startswith(("T", "index", "e_c", "C", "ind", "ec", "loop_temp","length")): value_fmt = ".2f" error_fmt = ".2f" if np.all(error) == np.all([0, 0]): From bd9fc78b763fdbdb58d946944dd8ecb5e9fd44b9 Mon Sep 17 00:00:00 2001 From: Iain Hannah Date: Fri, 17 Jul 2026 00:10:43 +0100 Subject: [PATCH 11/11] Create 286.bugfix.rst --- changelog/286.bugfix.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changelog/286.bugfix.rst diff --git a/changelog/286.bugfix.rst b/changelog/286.bugfix.rst new file mode 100644 index 00000000..5c48ccb3 --- /dev/null +++ b/changelog/286.bugfix.rst @@ -0,0 +1 @@ +Fixed a bug in legacy albedo correction that causes bumps at energy bin width changes, by applying photon bin width scaling to SRM instead of albedo matrix in :func:`sunkit_spex.legacy.fitting.fitter.make_model`.