Change parameter name from sigma to w

Author: lellid

src/ramanchada2/fitting_functions/models.py

@@ -1,6 +1,7 @@
 from lmfit import Model
 from lmfit import Parameter
 from lmfit.models import guess_from_peak, update_param_vals
+import numpy as np
 from ramanchada2.fitting_functions.pearsonivamplitudeparametrizationhpw import (
     PearsonIVAmplitudeParametrizationHPW,
 )
@@ -12,7 +13,7 @@
 class PearsonIVParametrizationHPWModel(Model):
     r"""A model based on a Pearson IV distribution.
     The model has five parameters: `height` (:math:`a`), `position`,
-    `sigma` (:math:`\sigma`), `expon` (:math:`m`) and `skew` (:math:`\nu`).
+    `w` (:math:`w`), `expon` (:math:`m`) and `skew` (:math:`\nu`).
     """
 
     def __init__(self, independent_vars=["x"], prefix="", nan_policy="raise", **kwargs):
@@ -33,15 +34,24 @@ def _set_paramhints_prefix(self):
     def guess(self, data, x, negative=False, **kwargs):
         """Estimate initial model parameter values from data."""
         pars = guess_from_peak(self, data, x, negative)
+        del pars['sigma']  # sigma is no longer needed
         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}w"].set(kwargs[f"{self.prefix}sigma"])
         pars[f"{self.prefix}expon"].set(value=1.0)
         pars[f"{self.prefix}skew"].set(value=0.0)
+
+        # Here a parameter sigma is added again, but this is merely
+        # to not make guess_from_peak crashing
+        # sigma is deleted from the dictionary after the call to 'guess_from_peak'
+        par = Parameter(name='sigma')
+        pars.add(par)
+
         return pars
 
     def fit(
@@ -77,7 +87,7 @@ def fit(
         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}w"],
             result.params[f"{self.prefix}expon"],
             result.params[f"{self.prefix}skew"],
             result.covar
@@ -112,16 +122,30 @@ def __init__(self, independent_vars=["x"], prefix="", nan_policy="raise", **kwar
         self._set_paramhints_prefix()
 
     def _set_paramhints_prefix(self):
-        self.set_param_hint("sigma", min=1E-100)
+        self.set_param_hint("w", min=1E-100)
         self.set_param_hint("nu", min=0.0, max=1.0)
         self.set_param_hint("gamma", expr='sigma*(1-nu)')
 
     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)
+        del pars['sigma']  # sigma is no longer needed
         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}w"].set(kwargs[f"{self.prefix}sigma"])
+
+        # Here a parameter sigma is added again, but this is merely
+        # to avoid crashing in the call to  'guess_from_peak'
+        # sigma is deleted from the dictionary after the call to 'guess_from_peak'
+        par = Parameter(name='sigma')
+        pars.add(par)
+
+        return pars
+
     def fit(
         self,
         data,
@@ -155,13 +179,20 @@ def fit(
         pahf = VoigtAreaParametrizationNu.GetPositionAreaHeightFWHMFromPeakParameters(
             result.params[f"{self.prefix}amplitude"],
             result.params[f"{self.prefix}center"],
-            result.params[f"{self.prefix}sigma"],
+            result.params[f"{self.prefix}w"],
             result.params[f"{self.prefix}nu"],
             result.covar
         )
         p1 = Parameter(f"{self.prefix}height", value=pahf.Height)
         p1.stderr = pahf.HeightStdDev
         p2 = Parameter(f"{self.prefix}fwhm", value=pahf.FWHM)
         p2.stderr = pahf.FWHMStdDev
-        result.params.add_many(p1, p2)
+        p3 = Parameter(f"{self.prefix}sigma",
+                       value=result.params[f"{self.prefix}w"] * np.sqrt(result.params[f"{self.prefix}nu"] / np.log(4)))
+        p3.stderr = 0.0  # TBD
+        p4 = Parameter(f"{self.prefix}gamma",
+                       value=result.params[f"{self.prefix}w"] * (1 - result.params[f"{self.prefix}nu"]))
+        p4.stderr = 0.0  # TBD
+
+        result.params.add_many(p1, p2, p3, p4)
         return result

src/ramanchada2/fitting_functions/pearsonivamplitudeparametrizationhpw.py

@@ -71,12 +71,12 @@ def GetParameterNamesForPeak(self, indexOfPeak):
             f"v{indexOfPeak}",
         ]
 
-    def GetYOfOneTerm(x, height=1.0, center=0.0, sigma=1.0, expon=1.0, skew=0.0):
+    def GetYOfOneTerm(x, height=1.0, center=0.0, w=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
+            / w
             - skew
         )
         return height * np.exp(

src/ramanchada2/fitting_functions/voigtareaparametrizationnu.py

@@ -95,10 +95,10 @@ def GetParameterNamesForPeak(self, indexOfPeak):
             f"nu{indexOfPeak}",
         ]
 
-    def GetYOfOneTerm(x, amplitude=1.0, center=0.0, sigma=1.0, nu=1.0):
+    def GetYOfOneTerm(x, amplitude=1.0, center=0.0, w=1.0, nu=1.0):
         """Returns the y-value of one peak in dependence on x and the peak parameters."""
         return amplitude * voigt_profile(
-            x - center, sigma * np.sqrt(nu) * OneBySqrtLog4, sigma * (1 - nu)
+            x - center, w * np.sqrt(nu) * OneBySqrtLog4, w * (1 - nu)
         )
 
     def func(self, pars, x, data=None):

tests/fitting_functions/test_pearsonivamplitudeparametrizationhpw.py

@@ -472,7 +472,7 @@ def test_fit_with_model_oneterm():
     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["w"], w, 7)
     np.testing.assert_almost_equal(result.params["expon"], m, 7)
     np.testing.assert_almost_equal(result.params["skew"], v, 7)
     np.testing.assert_almost_equal(result.params["amplitude"], 92.774482626349894596, 7)

tests/fitting_functions/test_voigtareaparametrizationnu.py

@@ -811,7 +811,9 @@ def test_fit_with_model_oneterm():
     result = model.fit(data, x=x, params=params)
     np.testing.assert_almost_equal(result.params["amplitude"], a, 3)
     np.testing.assert_almost_equal(result.params["center"], pos, 3)
-    np.testing.assert_almost_equal(result.params["sigma"], w, 3)
+    np.testing.assert_almost_equal(result.params["w"], w, 3)
     np.testing.assert_almost_equal(result.params["nu"], nu, 3)
+    np.testing.assert_almost_equal(result.params["sigma"], sigma, 2)
+    np.testing.assert_almost_equal(result.params["gamma"], gamma, 2)
     np.testing.assert_almost_equal(result.params["height"], 0.542599801404569082663570, 3)
     np.testing.assert_almost_equal(result.params["fwhm"], 10.111072472153491175127, 3)