From 615726ead4c40f3039fe59965fcfff586b85dac3 Mon Sep 17 00:00:00 2001
From: Danilo Ferreira de Lima <danilo.enoque.ferreira.de.lima@xfel.de>
Date: Fri, 20 Oct 2023 11:42:06 +0200
Subject: [PATCH] Bug fixes in resolution estimate.
---
pes_to_spec/__init__.py | 2 +-
pes_to_spec/model.py | 38 +++++++++++++++++++-------------------
2 files changed, 20 insertions(+), 20 deletions(-)
diff --git a/pes_to_spec/__init__.py b/pes_to_spec/__init__.py
index 02af2a4..d5e0f8e 100644
--- a/pes_to_spec/__init__.py
+++ b/pes_to_spec/__init__.py
@@ -2,4 +2,4 @@
Estimate high-resolution photon spectrometer data from low-resolution non-invasive measurements.
"""
-VERSION = "0.3.7"
+VERSION = "0.3.8"
diff --git a/pes_to_spec/model.py b/pes_to_spec/model.py
index 44d7a14..422a602 100644
--- a/pes_to_spec/model.py
+++ b/pes_to_spec/model.py
@@ -79,14 +79,14 @@ def deconv(y: np.ndarray, yhat: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.
def fit_gaussian(x: np.ndarray, y: np.ndarray) -> lmfit.ModelResult:
"""Fit Gaussian and return the fit result."""
- def gaussian(x, amplitude, centre, sigma):
- return amplitude * np.exp(-0.5 * (x - centre)**2 / (sigma**2))
+ def gaussian(x, amplitude, centre, log_sigma):
+ return amplitude * np.exp(-0.5 * ((x - centre)**2)*np.exp(-2*log_sigma))
gmodel = lmfit.Model(gaussian)
- result = gmodel.fit(y, x=x, centre=0.0, amplitude=np.amax(y), sigma=1.0)
+ result = gmodel.fit(y, x=x, centre=0.0, amplitude=np.amax(y), log_sigma=0.0)
return result
def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
- e_center: Optional[float]=None, e_width: Optional[float]=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, lmfit.ModelResult]:
+ e_center: Optional[float]=None, e_width: Optional[float]=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, lmfit.ModelResult]:
"""
Given the true y and the predicted y, together with the energy axis e,
estimate the impulse response of the system and return the Gaussian fit to it.
@@ -99,7 +99,7 @@ def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
e_center: If given the energy value, for which to probe the resolution.
e_width: The width of the energy neighbourhood to probe if e_center is given.
- Returns: The centered energy axis, the impulse response, the transfer function and the fit result.
+ Returns: The centered energy axis, the impulse response, the transfer function, the grating spectral density and the fit result.
"""
e_range = e[-1] - e[0]
e_axis = np.linspace(-0.5*e_range, 0.5*e_range, len(e))
@@ -112,7 +112,7 @@ def get_resolution(y: np.ndarray, y_hat: np.ndarray, e: np.ndarray,
y_sel = y_sel*f
y_hat_sel = y_hat_sel*f
h, H, S = deconv(y_sel, y_hat_sel)
- return e_axis, h, H, fit_gaussian(e_axis, np.absolute(h))
+ return e_axis, h, H, S, fit_gaussian(e_axis, np.absolute(h))
class PromptNotFoundError(Exception):
"""
@@ -796,7 +796,7 @@ class Model(TransformerMixin, BaseEstimator):
"n_obs",
"pca_threshold",
"high_res_fwhm",
- "resolution_per_energy", "resolution_per_energy_unc"
+ "resolution_per_energy", "resolution_per_energy_unc",
"resolution_energy_bins"
]
@@ -1014,10 +1014,10 @@ class Model(TransformerMixin, BaseEstimator):
e = high_res_photon_energy[0,:] if len(high_res_photon_energy.shape) == 2 else high_res_photon_energy
# get average resolution
- e_axis, h, H, result = get_resolution(y,
- y_hat,
- e
- )
+ e_axis, h, H, S, result = get_resolution(y,
+ y_hat,
+ e
+ )
de = e[1] - e[0]
E = np.fft.fftfreq(len(e), de)
@@ -1063,23 +1063,23 @@ class Model(TransformerMixin, BaseEstimator):
#if self.resolution < 0:
# print("Warning: Resolution calculation failed. The model can still be used, but this is probably a red flag!")
#
- self.resolution = result.best_values["sigma"]*2.355
+ self.resolution = np.exp(result.best_values["log_sigma"])*2.355
print("Resolution:", self.resolution)
# estimate resolution using power spectra
e_min = e.min()
e_max = e.max()
e_probe = np.linspace(e_min, e_max, 5)
- e_width = (e_max - e_min)/(len(e_probe[key])-1)
+ e_width = (e_max - e_min)/(len(e_probe)-1)
width = list()
width_unc = list()
for e_ in e_probe:
- e_axis_centered, h, H, result = get_resolution(y,
- y_hat,
- e_axis,
- e_center=e_,
- e_width=e_width)
- width += [result.best_values["sigma"]*2.355]
+ _, _, _, _, result = get_resolution(y,
+ y_hat,
+ e,
+ e_center=e_,
+ e_width=e_width)
+ width += [np.exp(result.best_values["log_sigma"])*2.355]
width_unc += [np.sqrt(result.covar[2,2])*2.355]
width = np.array(width)
width_unc = np.array(width_unc)
--
GitLab