Fit functions Voigt and PearsonIV using special parametrization

src/ramanchada2/fitting_functions/models.py

@@ -0,0 +1,198 @@
+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,
+)
+from ramanchada2.fitting_functions.voigtareaparametrizationnu import (
+    VoigtAreaParametrizationNu,
+)
+
+
+class PearsonIVParametrizationHPWModel(Model):
+    r"""A model based on a Pearson IV distribution.
+    The model has five parameters: `height` (:math:`a`), `position`,
+    `w` (:math:`w`), `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,
+            }
+        )
+        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)
+
+    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(
+        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}w"],
+            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.
+    The model has four Parameters: `amplitude`, `center`, `sigma`, and
+    `nu`. The parameter `nu` is constrained to the range [0,1]. In addition,
+    parameters `fwhm` are included as constraints to report
+    full width at half maximum.
+
+    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,
+            }
+        )
+        super().__init__(VoigtAreaParametrizationNu.GetYOfOneTerm, **kwargs)
+        self._set_paramhints_prefix()
+
+    def _set_paramhints_prefix(self):
+        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,
+        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 = VoigtAreaParametrizationNu.GetPositionAreaHeightFWHMFromPeakParameters(
+            result.params[f"{self.prefix}amplitude"],
+            result.params[f"{self.prefix}center"],
+            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
+        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/models_lmfit_like.py

@@ -0,0 +1,70 @@
+import numpy as np
+from lmfit import Model
+from lmfit.models import guess_from_peak, update_param_vals
+
+from .pearson4_hpw import (pearson4_hpw_area, pearson4_hpw_FWHM,
+                           pearson4_hpw_FWHM_approx, pearson4_hpw_LeftWHM,
+                           pearson4_hpw_RightWHM, perarson4_hpw_profile)
+from .voigt_nu import (voigt_nu_FWHM, voigt_nu_FWHM_approx, voigt_nu_height,
+                       voigt_nu_profile)
+
+
+class VoigtNuModel(Model):
+    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__(voigt_nu_profile, **kwargs)
+        self._set_paramhints_prefix()
+
+    def _set_paramhints_prefix(self):
+        self.set_param_hint('w', min=1E-100)
+        self.set_param_hint('nu', min=0.0, max=1.0)
+        self.set_param_hint('sigma', expr='w * sqrt(nu) / sqrt(log(4))')
+        self.set_param_hint('gamma', expr='sigma * (1-nu)')
+        self.set_param_hint('fwhm', expr='voigt_nu_FWHM(w=sigma, nu=nu)')
+        self.set_param_hint('fwhm_approx', expr='voigt_nu_FWHM_approx(w=sigma, nu=nu)')
+        self.set_param_hint('fwhm_gamma_sigma', expr='1.0692*gamma+sqrt(0.8664*gamma**2+5.545083*sigma**2)')
+        self.set_param_hint('height', expr='voigt_nu_height(amplitude=amplitude, sigma=sigma, nu=nu)')
+
+    def guess(self, data, x, negative=False, **kwargs):
+        """Estimate initial model parameter values from data."""
+        pars = self.make_params(amplitude=np.max(data), center=x[np.argmax(data)], nu=.5, w=10)
+        return update_param_vals(pars, self.prefix, **kwargs)
+
+    def make_params(self, *args, **kwargs):
+        params = super().make_params(*args, **kwargs)
+        params._asteval.symtable['voigt_nu_height'] = voigt_nu_height
+        params._asteval.symtable['voigt_nu_FWHM'] = voigt_nu_FWHM
+        params._asteval.symtable['voigt_nu_FWHM_approx'] = voigt_nu_FWHM_approx
+        return params
+
+
+class Pearson4HPWModel(Model):
+    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__(perarson4_hpw_profile, **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('amplitude', expr='pearson4_hpw_area(height, sigma, expon, skew)')
+        self.set_param_hint('fwhm', expr='pearson4_hpw_FWHM(sigma, expon, skew)')
+        self.set_param_hint('fwhm_approx', expr='pearson4_hpw_FWHM_approx(sigma, expon, skew)')
+
+    def guess(self, data, x, negative=False, **kwargs):
+        """Estimate initial model parameter values from data."""
+        pars = guess_from_peak(self, data, x, negative)
+        return update_param_vals(pars, self.prefix, **kwargs)
+
+    def make_params(self, *args, **kwargs):
+        params = super().make_params(*args, **kwargs)
+        params._asteval.symtable['pearson4_hpw_area'] = pearson4_hpw_area
+        params._asteval.symtable['pearson4_hpw_FWHM_approx'] = pearson4_hpw_FWHM_approx
+        params._asteval.symtable['pearson4_hpw_LeftWHM'] = pearson4_hpw_LeftWHM
+        params._asteval.symtable['pearson4_hpw_RightWHM'] = pearson4_hpw_RightWHM
+        params._asteval.symtable['pearson4_hpw_FWHM'] = pearson4_hpw_FWHM
+        return params