Pearson4 model now functional, amplitude and fwhm added

Author: lellid

src/ramanchada2/fitting_functions/models.py

@@ -1,7 +1,12 @@
 from lmfit import Model
+from lmfit import Parameter
 from lmfit.models import guess_from_peak, update_param_vals
-from pearsonivamplitudeparametrizationhpw import PearsonIVAmplitudeParametrizationHPW
-from voigtareaparametrizationnu import VoigtAreaParametrizationNu
+from ramanchada2.fitting_functions.pearsonivamplitudeparametrizationhpw import (
+    PearsonIVAmplitudeParametrizationHPW,
+)
+from ramanchada2.fitting_functions.voigtareaparametrizationnu import (
+    VoigtAreaParametrizationNu,
+)
 import numpy as np
 
 
@@ -11,25 +16,80 @@ class PearsonIVParametrizationHPWModel(Model):
     `sigma` (:math:`\sigma`), `expon` (:math:`m`) and `skew` (:math:`\nu`).
     """
 
-    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
-                 **kwargs):
-        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
-                       'independent_vars': independent_vars})
+    def __init__(self, independent_vars=["x"], prefix="", nan_policy="raise", **kwargs):
+        kwargs.update(
+            {
+                "prefix": prefix,
+                "nan_policy": nan_policy,
+                "independent_vars": independent_vars,
+            }
+        )
         super().__init__(PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm, **kwargs)
         self._set_paramhints_prefix()
 
     def _set_paramhints_prefix(self):
-        self.set_param_hint('expon', value=1.5, min=0.5 + 1 / 1024.0, max=1000)
-        self.set_param_hint('skew', value=0.0, min=-1000, max=1000)
+        self.set_param_hint("expon", value=1.5, min=0.5 + 1 / 1024.0, max=1000)
+        self.set_param_hint("skew", value=0.0, min=-1000, max=1000)
 
     def guess(self, data, x, negative=False, **kwargs):
         """Estimate initial model parameter values from data."""
         pars = guess_from_peak(self, data, x, negative)
-        pars[f'{self.prefix}height'].set(value=pars[f'{self.prefix}amplitude'] / pars[f'{self.prefix}sigma'])
-        pars[f'{self.prefix}expon'].set(value=1)
-        pars[f'{self.prefix}skew'].set(value=0.0)
         return update_param_vals(pars, self.prefix, **kwargs)
 
+    def make_params(self, verbose=False, **kwargs):
+        pars = super().make_params(verbose=verbose, **kwargs)
+        pars[f"{self.prefix}height"].set(
+            value=kwargs[f"{self.prefix}amplitude"] / kwargs[f"{self.prefix}sigma"]
+        )
+        pars[f"{self.prefix}expon"].set(value=1.0)
+        pars[f"{self.prefix}skew"].set(value=0.0)
+        return pars
+
+    def fit(
+        self,
+        data,
+        params=None,
+        weights=None,
+        method="leastsq",
+        iter_cb=None,
+        scale_covar=True,
+        verbose=False,
+        fit_kws=None,
+        nan_policy=None,
+        calc_covar=True,
+        max_nfev=None,
+        **kwargs,
+    ):
+        """overwrite fit in order to amend area and fwhm to the parameters"""
+        result = super().fit(
+            data,
+            params=params,
+            weights=weights,
+            method=method,
+            iter_cb=iter_cb,
+            scale_covar=scale_covar,
+            verbose=verbose,
+            fit_kws=fit_kws,
+            nan_policy=nan_policy,
+            calc_covar=calc_covar,
+            max_nfev=max_nfev,
+            **kwargs,
+        )
+        pahf = PearsonIVAmplitudeParametrizationHPW.GetPositionAreaHeightFWHMFromPeakParameters(
+            result.params[f"{self.prefix}height"],
+            result.params[f"{self.prefix}center"],
+            result.params[f"{self.prefix}sigma"],
+            result.params[f"{self.prefix}expon"],
+            result.params[f"{self.prefix}skew"],
+            result.covar
+        )
+        p1 = Parameter(f"{self.prefix}amplitude", value=pahf.Area)
+        p1.stderr = pahf.AreaStdDev
+        p2 = Parameter(f"{self.prefix}fwhm", value=pahf.FWHM)
+        p2.stderr = pahf.FWHMStdDev
+        result.params.add_many(p1, p2)
+        return result
+
 
 class VoigtAreaParametrizationNuModel(Model):
     r"""A model based on a Voigt distribution function.
@@ -41,22 +101,26 @@ class VoigtAreaParametrizationNuModel(Model):
     For more information, see: https://en.wikipedia.org/wiki/Voigt_profile
     """
 
-    def __init__(self, independent_vars=['x'], prefix='', nan_policy='raise',
-                 **kwargs):
-        kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
-                       'independent_vars': independent_vars})
+    def __init__(self, independent_vars=["x"], prefix="", nan_policy="raise", **kwargs):
+        kwargs.update(
+            {
+                "prefix": prefix,
+                "nan_policy": nan_policy,
+                "independent_vars": independent_vars,
+            }
+        )
         super().__init__(VoigtAreaParametrizationNu.GetYOfOneTerm, **kwargs)
         self._set_paramhints_prefix()
 
     def _set_paramhints_prefix(self):
-        self.set_param_hint('sigma', min=np.tiny)
-        self.set_param_hint('nu', min=0, max=1)
+        self.set_param_hint("sigma", min=np.tiny)
+        self.set_param_hint("nu", min=0, max=1)
 
-        fexpr = ("2.0*{pre:s}sigma")
-        self.set_param_hint('fwhm', expr=fexpr.format(pre=self.prefix))
+        fexpr = "2.0*{pre:s}sigma"
+        self.set_param_hint("fwhm", expr=fexpr.format(pre=self.prefix))
 
     def guess(self, data, x, negative=False, **kwargs):
         """Estimate initial model parameter values from data."""
         pars = guess_from_peak(self, data, x, negative)
-        pars[f'{self.prefix}nu'].set(value=1.0)
+        pars[f"{self.prefix}nu"].set(value=1.0)
         return update_param_vals(pars, self.prefix, **kwargs)

src/ramanchada2/fitting_functions/pearsonivamplitudeparametrizationhpw.py

@@ -73,7 +73,12 @@ def GetParameterNamesForPeak(self, indexOfPeak):
 
     def GetYOfOneTerm(x, height=1.0, center=0.0, sigma=1.0, expon=1.0, skew=0.0):
         """Returns the y-value of one peak in dependence on x and the peak parameters."""
-        arg = np.sqrt((np.power(2, 1 / expon) - 1) * (1 + skew * skew)) * (x - center) / sigma - skew
+        arg = (
+            np.sqrt((np.power(2, 1 / expon) - 1) * (1 + skew * skew))
+            * (x - center)
+            / sigma
+            - skew
+        )
         return height * np.exp(
             expon
             * (
@@ -397,12 +402,28 @@ def GetPositionAreaHeightFWHMFromSinglePeakParameters(
         wm = parameters[f"w{indexOfPeak}"]
         m = parameters[f"m{indexOfPeak}"]
         vm = parameters[f"v{indexOfPeak}"]
+
         if cv is not None:
             cv = cv[
                 indexOfPeak * 5: (indexOfPeak + 1) * 5,
                 indexOfPeak * 5: (indexOfPeak + 1) * 5,
             ]
+        return PearsonIVAmplitudeParametrizationHPW.GetPositionAreaHeightFWHMFromPeakParameters(
+            amp, loc, wm, m, vm, cv
+        )
 
+    def GetPositionAreaHeightFWHMFromPeakParameters(amp, loc, wm, m, vm, cv=None):
+        """
+        Get position, area, height, and FWHM of one peak from the peak parameters.
+        If the covariance matrix is given, the corresponding errors are calculated.
+
+        Parameters
+        ----------
+        amp, loc, wm, m, vm:
+            Peak parameters height, center, sigma, m, v
+        cv: np.array
+            Covariance matrix of the fit (5 x 5 matrix)
+        """
         area = PearsonIVAmplitudeParametrizationHPW.GetArea(amp, wm, m, vm)
         pos = loc
         height = amp

tests/fitting_functions/test_pearsonivamplitudeparametrizationhpw.py

@@ -1,6 +1,7 @@
 from ramanchada2.fitting_functions.pearsonivamplitudeparametrizationhpw import (
     PearsonIVAmplitudeParametrizationHPW,
-)
+    )
+from ramanchada2.fitting_functions.models import PearsonIVParametrizationHPWModel
 import lmfit
 import numpy as np
 
@@ -456,3 +457,21 @@ def test_fit_with_derivatives_twoterms_linearbaseline():
     np.testing.assert_almost_equal(out1.params["w1"], w1, 12)
     np.testing.assert_almost_equal(out1.params["m1"], m1, 12)
     np.testing.assert_almost_equal(out1.params["v1"], v1, 12)
+
+
+def test_fit_with_model_oneterm():
+    x = np.linspace(0, 99, 100)
+    a = 7
+    pos = 50
+    w = 5
+    m = 16 / 11.0
+    v = 1 / 4.0
+    data = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(x, a, pos, w, m, v)
+    model = PearsonIVParametrizationHPWModel(x=x)
+    params = model.guess(data, x)
+    result = model.fit(data, x=x, params=params)
+    np.testing.assert_almost_equal(result.params["height"], a, 7)
+    np.testing.assert_almost_equal(result.params["center"], pos, 7)
+    np.testing.assert_almost_equal(result.params["sigma"], w, 7)
+    np.testing.assert_almost_equal(result.params["expon"], m, 7)
+    np.testing.assert_almost_equal(result.params["skew"], v, 7)