diff --git a/src/ramanchada2/fitting_functions/__init__.py b/src/ramanchada2/fitting_functions/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/ramanchada2/fitting_functions/models.py b/src/ramanchada2/fitting_functions/models.py new file mode 100644 index 00000000..d1a6c330 --- /dev/null +++ b/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 diff --git a/src/ramanchada2/fitting_functions/models_lmfit_like.py b/src/ramanchada2/fitting_functions/models_lmfit_like.py new file mode 100644 index 00000000..8dd8974b --- /dev/null +++ b/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 diff --git a/src/ramanchada2/fitting_functions/pearson4_hpw.py b/src/ramanchada2/fitting_functions/pearson4_hpw.py new file mode 100644 index 00000000..b1fd06f2 --- /dev/null +++ b/src/ramanchada2/fitting_functions/pearson4_hpw.py @@ -0,0 +1,159 @@ +import numpy as np +from scipy.special import betaln, gammaln, loggamma + + +def perarson4_hpw_profile(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.""" + np.sqrt((2**(1/expon)-1) * (1+skew**2)) * (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 + * ( + np.log((1 + skew * skew) / (1 + arg * arg)) + - 2 * skew * (np.arctan(arg) + np.arctan(skew)) + ) + ) + + +def pearson4_hpw_area(amp, w, m, v): + lnprefactor = ( + betaln(m - 0.5, 0.5) + + 2 * gammaln(m) + - 2 * np.real(loggamma(m + m * v * 1.0j)) + ) + tmm1 = np.power(2, 1 / m) - 1 + lnbody = ( + -2 * m * v * np.arctan(v) + + (m - 0.5) * np.log(1 + v * v) + - 0.5 * np.log(tmm1) + ) + return amp * w * np.exp(lnprefactor + lnbody) + + +def pearson4_hpw_FWHM_approx(w, m, v): + """ + Gets the approximate full-width-half-maximum of a peak + in dependence on the parameters 'w', 'm' and 'v'. + The approximation has an accuracy of 39%. + """ + sq = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) + ws = w / sq + ms = m + vs = 2 * m * v + return ( + ws + * np.sqrt(np.power(2, 1 / ms) - 1) + * (np.pi / np.arctan2(np.exp(1) * ms, np.abs(vs))) + ) + + +def pearson4_hpw_LeftWHM(w, m, v): + sq = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) + return _pearson4_original_HFHM(w / sq, m, 2 * m * v, False) + + +def pearson4_hpw_RightWHM(w, m, v): + sq = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) + return _pearson4_original_HFHM(w / sq, m, 2 * m * v, True) + + +def pearson4_hpw_FWHM(w, m, v): + """ + Gets the full-width-half-maximum of a peak + in dependence on the parameters 'w', 'm' and 'v'. + """ + return pearson4_hpw_LeftWHM(w, m, v) + pearson4_hpw_RightWHM(w, m, v) + + +def _pearson4_original_HFHM(w, m, v, rightSide): + """ + Gets the half-width-half-maximum of the side of a peak + in dependence on the parameters 'w', 'm' and 'v'. + ATTENTION: Here, the parameters 'w', 'm' and 'v' + are from the original definition of PearsonIV, not + the parametrized definition used in this class! + If 'rightSide' is true, the HWHM of the right side + of the peak is returned, else that of the left side. + """ + if not m > 0: + return np.NaN + + w = np.abs(w) + sign = 1 if rightSide else -1 + z0 = -v / (2 * m) + + def funcsimp(z, m, v): + return ( + np.log(2) + + v * (np.arctan(z0) - np.arctan(z)) + + m + * ( + np.log(1 + z0 * z0) + - ( + 2 * np.log(np.abs(z)) + if np.abs(z) > 1e100 + else np.log(1 + z * z) + ) + ) + ) + + def dervsimp(z, m, v): + return ( + (-v - 2 * m * z) / (1 + z * z) + if np.abs(z) < 1e100 + else (-v / z - 2 * m) / z + ) + + # go forward in exponentially increasing steps, until the amplitude falls below ymaxHalf, + # in order to bracked the solution + zNear = z0 + zFar = z0 + d = 1.0 + + while np.isfinite(d): + zFar = z0 + d * sign + y = funcsimp(zFar, m, v) + if y < 0: + break + else: + zNear = zFar + d *= 2 + if zNear > zFar: + (zNear, zFar) = (zFar, zNear) + + # use Newton-Raphson to refine the result + z = 0.5 * (zNear + zFar) # starting value + for i in range(40): + funcVal = funcsimp(z, m, v) + if rightSide: + if funcVal > 0 and z > zNear: + zNear = z + if funcVal < 0 and z < zFar: + zFar = z + else: # leftSide + if funcVal < 0 and z > zNear: + zNear = z + if funcVal > 0 and z < zFar: + zFar = z + + dz = funcVal / dervsimp(z, m, v) + znext = z - dz + if znext <= zNear: + znext = (z + zNear) / 2 + elif znext >= zFar: + znext = (z + zFar) / 2 + + if z == znext: + break + + z = znext + + if np.abs(dz) < 5e-15 * np.abs(z): + break + return np.abs((z - z0) * w) diff --git a/src/ramanchada2/fitting_functions/pearsonivamplitudeparametrizationhpw.py b/src/ramanchada2/fitting_functions/pearsonivamplitudeparametrizationhpw.py new file mode 100644 index 00000000..cb26209d --- /dev/null +++ b/src/ramanchada2/fitting_functions/pearsonivamplitudeparametrizationhpw.py @@ -0,0 +1,549 @@ +#!/usr/bin/env python3 + +import numpy as np +from ramanchada2.fitting_functions.voigtareaparametrizationnu import ( + PositionAreaHeightFWHM, +) +from scipy.special import betaln, gammaln, loggamma + +Log2 = 0.69314718055994530941723212145818 # Math.Log(2) + + +def SafeSqrt(x): + """Returns 0 if x is negative, otherwise the square root of x""" + return 0.0 if x < 0 else np.sqrt(x) + + +class PearsonIVAmplitudeParametrizationHPW: + """ + PearsonIV peak fitting function with multiple peaks and polynomial baseline + + A special parametrization is used for the PearsonIV function: + Each of the peaks has 5 parameters: a, pos, w, m, n + Parameter a designates the height of the peak (maximum y-value) + Parameter 'pos' designates the position of the peak (x-value at maximum) + Parameter 'w' is approximately (38%) the HWHM value of the peak + Parameter 'm' determines the Gauss'ness of the peak (1=Lorentian, the higher, + the more gaussian) + Parameter v determines the skewness + + Note: the higher the value of 'm', the less influence the skewness parameter 'v' has + that means that a skewed Gaussian peak can not be modeled with this function + """ + + # Number of peaks + numberOfTerms = 1 + + # Order of the baseline polynomial (-1=no baseline, 0=constant, 1=linear, etc.) + orderOfBaselinePolynomial = -1 + + def __init__(self, numberOfTerms=1, orderOfBackgroundPolynomial=-1): + """ + Initializes an instance of the class with a given number of + peak terms, and the order of the baseline polynomial + """ + self.numberOfTerms = numberOfTerms + self.orderOfBaselinePolynomial = orderOfBackgroundPolynomial + + def WithNumberOfTerms(self, numberOfTerms): + """Returns a new instance of the class with the given number of peak terms.""" + return PearsonIVAmplitudeParametrizationHPW( + numberOfTerms, self.orderOfBaselinePolynomial + ) + + def WithOrderOfBaselinePolynomial(self, orderOfBaselinePolynomial): + """Returns a new instance of the class with the given order of the baseline polynomial.""" + return PearsonIVAmplitudeParametrizationHPW( + self.numberOfTerms, orderOfBaselinePolynomial + ) + + def GetNumberOfParametersPerPeak(self): + """Returns the number of parameters per peak term.""" + return 5 + + def GetParameterNamesForPeak(self, indexOfPeak): + """Returns the parameter names for a given peak index.""" + return [ + f"a{indexOfPeak}", + f"pos{indexOfPeak}", + f"w{indexOfPeak}", + f"m{indexOfPeak}", + f"v{indexOfPeak}", + ] + + 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) + / w + - skew + ) + return height * np.exp( + expon + * ( + np.log((1 + skew * skew) / (1 + arg * arg)) + - 2 * skew * (np.arctan(arg) + np.arctan(skew)) + ) + ) + + def func(self, pars, x, data=None): + """Returns the y-values of the fitting function.""" + sum = np.zeros(len(x)) + for i in range(self.orderOfBaselinePolynomial, -1, -1): + sum *= x + sum += pars[f"b{i}"] + + for i in range(self.numberOfTerms): + a, xc, w, m, v = ( + pars[f"a{i}"], + pars[f"pos{i}"], + pars[f"w{i}"], + pars[f"m{i}"], + pars[f"v{i}"], + ) + arg = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) * (x - xc) / w - v + sum += a * np.exp( + m + * ( + np.log((1 + v * v) / (1 + arg * arg)) + - 2 * v * (np.arctan(arg) + np.arctan(v)) + ) + ) + + if data is None: + return sum + return sum - data + + def dfunc(self, pars, x, data=None): + """Returns the derivatives of the fitting function w.r.t. the parameters.""" + result = [] + + # first the baseline derivatives + if self.orderOfBaselinePolynomial >= 0: + xn = np.ones(len(x)) + for i in range(self.orderOfBaselinePolynomial + 1): + if pars[f"b{i}"].vary: + result.append(np.copy(xn)) + xn *= x + + # second the derivatives of the peak terms + for i in range(self.numberOfTerms): + height, pos, w, m, v = ( + pars[f"a{i}"], + pars[f"pos{i}"], + pars[f"w{i}"], + pars[f"m{i}"], + pars[f"v{i}"], + ) + if not (height.vary or pos.vary or w.vary or m.vary or v.vary): + continue + + twoToOneByM_1 = np.power(2, 1 / m) - 1 + ww = w / np.sqrt(twoToOneByM_1 * (1 + v * v)) + z = (x - pos) / ww - v + log1v2_1z2 = np.log((1 + v * v) / (1 + z * z)) + atanZ_p_V = 2 * v * (np.arctan(z) + np.arctan(v)) + body = np.exp(m * (log1v2_1z2 - atanZ_p_V)) + dbodydz = -body * (2 * m * (z + v)) / (1 + z * z) + + dfdheight = np.full(len(x), np.NaN) + dfdpos = np.full(len(x), np.NaN) + dfdw = np.full(len(x), np.NaN) + dfdm = np.full(len(x), np.NaN) + dfdv = np.full(len(x), np.NaN) + + dfdheight = body + dfdpos = -height * dbodydz / ww + dfdw = -height * dbodydz * (z + v) / w + dfdm = ( + height + * body + * ( + -atanZ_p_V + + log1v2_1z2 + + Log2 + * np.square(z + v) + * np.power(2, 1 / m) + / (twoToOneByM_1 * (1 + z * z) * m) + ) + ) + dfdv = ( + -height + * 2 + * m + * body + * ( + (z + v) * (z * v - 1) / ((1 + z * z) * (1 + v * v)) + + np.arctan(z) + + np.arctan(v) + ) + ) + + if height.vary: + result.append(dfdheight) + if pos.vary: + result.append(dfdpos) + if w.vary: + result.append(dfdw) + if m.vary: + result.append(dfdm) + if v.vary: + result.append(dfdv) + + return np.array(result) + + def GetInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + height, position, fullWidth, relativeHeight=0.5 + ): + """ + Returns an array of initial parameters for peak fitting, + given the height, the position, the full width and the relative + height at which the full width was measured. + """ + if not (relativeHeight > 0 and relativeHeight < 1): + raise ValueError("RelativeHeight should be in the open interval (0, 1)") + + # we evaluate the parameters here for a pure Lorentzian (m=1, v=0) + w = np.abs(0.5 * fullWidth * np.sqrt(relativeHeight / (1 - relativeHeight))) + return [height, position, w, 1, 0] # Parameters for the Lorentz limit + + def AddInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + self, paras, indexOfPeak, height, position, fullWidth, relativeHeight=0.5 + ): + """ + Add to an existing 'Parameters' instance the initial parameters for a peak. + + Parameters + ---------- + paras: Parameters + Existing instance of Parameters, to which the new parameters will be added. + indexOfPeak: + The index of the peak whose parameters are added. + height: + Approximate height of the peak. + pos: + Approximate position of the peak + fullWidth: + Approximate full width of the peak + relativeHeight: + Relative height, at which full width was measured (usually 0.5) + """ + apwmv = PearsonIVAmplitudeParametrizationHPW.GetInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + height, position, fullWidth, relativeHeight + ) + paras.add(f"a{indexOfPeak}", apwmv[0]) + paras.add(f"pos{indexOfPeak}", apwmv[1]) + paras.add(f"w{indexOfPeak}", apwmv[2]) + paras.add(f"m{indexOfPeak}", apwmv[3]) + paras.add(f"v{indexOfPeak}", apwmv[4]) + + def AddInitialParametersForBaselinePolynomial(self, paras): + """Add the required parameters for the baseline polynomial""" + for i in range(self.orderOfBaselinePolynomial + 1): + paras.add(f"b{i}", 0.0) + + def GetFWHMApproximation(w, m, v): + """ + Gets the approximate full-width-half-maximum of a peak + in dependence on the parameters 'w', 'm' and 'v'. + The approximation has an accuracy of 39%. + """ + sq = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) + ws = w / sq + ms = m + vs = 2 * m * v + return ( + ws + * np.sqrt(np.power(2, 1 / ms) - 1) + * (np.pi / np.arctan2(np.exp(1) * ms, np.abs(vs))) + ) + + def GetFWHM(w, m, v): + """ + Gets the full-width-half-maximum of a peak + in dependence on the parameters 'w', 'm' and 'v'. + """ + return PearsonIVAmplitudeParametrizationHPW.GetHWHM( + w, m, v, True + ) + PearsonIVAmplitudeParametrizationHPW.GetHWHM(w, m, v, False) + + def GetHWHM(w, m, v, rightSide): + """ + Gets the half-width-half-maximum of the side of a peak + in dependence on the parameters 'w', 'm' and 'v'. + If 'rightSide' is true, the HWHM of the right side + of the peak is returned, else that of the left side. + """ + sq = np.sqrt((np.power(2, 1 / m) - 1) * (1 + v * v)) + return PearsonIVAmplitudeParametrizationHPW.GetHWHMOfPearson4Original( + w / sq, m, 2 * m * v, rightSide + ) + + def GetHWHMOfPearson4Original(w, m, v, rightSide): + """ + Gets the half-width-half-maximum of the side of a peak + in dependence on the parameters 'w', 'm' and 'v'. + ATTENTION: Here, the parameters 'w', 'm' and 'v' + are from the original definition of PearsonIV, not + the parametrized definition used in this class! + If 'rightSide' is true, the HWHM of the right side + of the peak is returned, else that of the left side. + """ + if not m > 0: + return np.NaN + + w = np.abs(w) + sign = 1 if rightSide else -1 + z0 = -v / (2 * m) + + def funcsimp(z, m, v): + return ( + np.log(2) + + v * (np.arctan(z0) - np.arctan(z)) + + m + * ( + np.log(1 + z0 * z0) + - ( + 2 * np.log(np.abs(z)) + if np.abs(z) > 1e100 + else np.log(1 + z * z) + ) + ) + ) + + def dervsimp(z, m, v): + return ( + (-v - 2 * m * z) / (1 + z * z) + if np.abs(z) < 1e100 + else (-v / z - 2 * m) / z + ) + + # go forward in exponentially increasing steps, until the amplitude falls below ymaxHalf, + # in order to bracked the solution + zNear = z0 + zFar = z0 + d = 1.0 + + while np.isfinite(d): + zFar = z0 + d * sign + y = funcsimp(zFar, m, v) + if y < 0: + break + else: + zNear = zFar + d *= 2 + if zNear > zFar: + (zNear, zFar) = (zFar, zNear) + + # use Newton-Raphson to refine the result + z = 0.5 * (zNear + zFar) # starting value + for i in range(40): + funcVal = funcsimp(z, m, v) + if rightSide: + if funcVal > 0 and z > zNear: + zNear = z + if funcVal < 0 and z < zFar: + zFar = z + else: # leftSide + if funcVal < 0 and z > zNear: + zNear = z + if funcVal > 0 and z < zFar: + zFar = z + + dz = funcVal / dervsimp(z, m, v) + znext = z - dz + if znext <= zNear: + znext = (z + zNear) / 2 + elif znext >= zFar: + znext = (z + zFar) / 2 + + if z == znext: + break + + z = znext + + if np.abs(dz) < 5e-15 * np.abs(z): + break + return np.abs((z - z0) * w) + + def GetArea(amp, w, m, v): + lnprefactor = ( + betaln(m - 0.5, 0.5) + + 2 * gammaln(m) + - 2 * np.real(loggamma(m + m * v * 1.0j)) + ) + tmm1 = np.power(2, 1 / m) - 1 + lnbody = ( + -2 * m * v * np.arctan(v) + + (m - 0.5) * np.log(1 + v * v) + - 0.5 * np.log(tmm1) + ) + return amp * w * np.exp(lnprefactor + lnbody) + + def GetPositionAreaHeightFWHMFromSinglePeakParameters( + self, parameters, indexOfPeak, 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 + ---------- + parameters: Parameters + Existing instance of Parameters, which contains the result of the fit. + indexOfPeak: int + The index of the peak into consideration. + cv: np.array + Covariance matrix of the fit. + """ + amp = parameters[f"a{indexOfPeak}"] + loc = parameters[f"pos{indexOfPeak}"] + 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 + fwhm = PearsonIVAmplitudeParametrizationHPW.GetHWHM( + wm, m, vm, True + ) + PearsonIVAmplitudeParametrizationHPW.GetHWHM(wm, m, vm, False) + + posStdDev = 0 + areaStdDev = 0 + heightStdDev = 0 + fwhmStdDev = 0 + + if cv is not None: + deriv = np.zeros(5) + resVec = np.zeros(5) + + # PositionVariance + posStdDev = SafeSqrt(cv[1, 1]) + + # Height variance + heightStdDev = SafeSqrt(cv[0, 0]) + + # Area variance + deriv[0] = PearsonIVAmplitudeParametrizationHPW.GetArea(1, wm, m, vm) + deriv[1] = 0 + absDelta = wm * 1e-5 + deriv[2] = ( + PearsonIVAmplitudeParametrizationHPW.GetArea(amp, wm + absDelta, m, vm) + - PearsonIVAmplitudeParametrizationHPW.GetArea( + amp, wm - absDelta, m, vm + ) + ) / (2 * absDelta) + absDelta = m * 1e-5 + deriv[3] = ( + PearsonIVAmplitudeParametrizationHPW.GetArea(amp, wm, m + absDelta, vm) + - PearsonIVAmplitudeParametrizationHPW.GetArea( + amp, wm, m - absDelta, vm + ) + ) / (2 * absDelta) + absDelta = 1e-5 if vm == 0 else np.abs(vm * 1e-5) + deriv[4] = ( + PearsonIVAmplitudeParametrizationHPW.GetArea(amp, wm, m, vm + absDelta) + - PearsonIVAmplitudeParametrizationHPW.GetArea( + amp, wm, m, vm - absDelta + ) + ) / (2 * absDelta) + resVec = cv.dot(deriv) + areaStdDev = SafeSqrt(np.dot(deriv, resVec)) + + # FWHM variance + deriv[0] = 0 + deriv[1] = 0 + absDelta = wm * 1e-5 + deriv[2] = ( + PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm + absDelta, m, vm) + - PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm - absDelta, m, vm) + ) / (2 * absDelta) + absDelta = m * 1e-5 + deriv[3] = ( + PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm, m + absDelta, vm) + - PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm, m - absDelta, vm) + ) / (2 * absDelta) + absDelta = 1e-5 if vm == 0 else np.abs(vm * 1e-5) + deriv[4] = ( + PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm, m, vm + absDelta) + - PearsonIVAmplitudeParametrizationHPW.GetFWHM(wm, m, vm - absDelta) + ) / (2 * absDelta) + resVec = cv.dot(deriv) + fwhmStdDev = SafeSqrt(np.dot(deriv, resVec)) + + return PositionAreaHeightFWHM( + pos, posStdDev, area, areaStdDev, height, heightStdDev, fwhm, fwhmStdDev + ) + + def SetParameterBoundariesForPositivePeaks( + self, + params, + minimalPosition=None, + maximalPosition=None, + minimalFWHM=None, + maximalFWHM=None, + ): + """ + Amends the parameters in an existing Parameters instance + with fit boundaries. + + Parameters + ---------- + params: Parameters + Existing instance of Parameters, which already contain all parameters. + minimalPosition: float (or None) + The minimal allowed position of a peak. + maximalPosition: float (or None) + The maximal allowed position of a peak. + minimalFWHM: float (or None) + The minimal allowed FWHM value of a peak. + maximalFWHM: float (or None) + The maximal allowed FWHM value of a peak. + + Note: the minimal height is set to 0 in order to disallow negative peaks. + the maximal height of a peak is not limited + """ + DefaultMinWidth = 1.4908919308538355e-81 # Math.Pow(double.Epsilon, 0.25); + DefaultMaxWidth = 1.157920892373162e77 # Math.Pow(double.MaxValue, 0.25); + + for i in range(self.numberOfTerms): + params[f"a{i}"].min = 0 # minimal amplitude is 0 + if minimalPosition is not None: + params[f"pos{i}"].min = minimalPosition + if maximalPosition is not None: + params[f"pos{i}"].max = maximalPosition + if minimalFWHM is not None: + params[f"w{i}"].min = minimalFWHM / 2.0 + else: + params[f"w{i}"].min = DefaultMinWidth + if maximalFWHM is not None: + params[f"w{i}"].max = maximalFWHM / 2.0 + else: + params[f"w{i}"].max = DefaultMaxWidth + params[f"m{i}"].min = 1 / 2.0 + 1 / 1024.0 + params[f"m{i}"].max = 1000 + params[f"v{i}"].min = -1000 + params[f"v{i}"].max = 1000 diff --git a/src/ramanchada2/fitting_functions/problem.ipynb b/src/ramanchada2/fitting_functions/problem.ipynb new file mode 100644 index 00000000..2ef65f23 --- /dev/null +++ b/src/ramanchada2/fitting_functions/problem.ipynb @@ -0,0 +1,204 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "7b91f3b1-d08b-4d84-a241-819d11fa9301", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import ramanchada2 as rc2\n", + "from ramanchada2.fitting_functions import models_lmfit_like as rc2models_lm\n", + "from ramanchada2.fitting_functions import models as rc2models\n", + "from ramanchada2.io.simulated import read_simulated_lines\n", + "from typing import Dict\n", + "import numpy as np\n", + "import pandas as pd\n", + "import lmfit\n", + "import matplotlib.pyplot as plt\n", + "from scipy import fft, signal, optimize, sparse\n", + "from scipy.sparse import diags\n", + "def plt_fig_ax(nrows=1, ncols=1, figsize=(20,5), **kwargs):\n", + " return plt.subplots(nrows=nrows, ncols=ncols, figsize=figsize, **kwargs)\n", + "from ramanchada2.auxiliary.spectra.datasets2 import get_filenames" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc744b17-3cf9-4b39-8379-2aac75674e0c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "x = np.linspace(0, 200, 1000, endpoint=False)\n", + "mod = lmfit.models.Pearson4Model()\n", + "\n", + "par = mod.make_params()\n", + "par['amplitude'].set(value=1000)\n", + "par['center'].set(value=100)\n", + "par['sigma'].set(value=5)\n", + "par['skew'].set(value=10)\n", + "par['expon'].set(value=4)\n", + "arr = mod.eval(params=par, x=x)\n", + "arr += np.random.uniform(-3,3, size=len(arr))\n", + "plt.plot(x,arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69ab3608-6277-44c5-acaa-ed1084373e02", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "a=rc2models_lm.Pearson4HPWModel()\n", + "\n", + "\n", + "params = a.guess(arr, x=x)\n", + "\n", + "fit_res = a.fit(data=arr, x=x, params=params)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fit_res.eval())\n", + "\n", + "fit_res.params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "436ca200-f24d-4fde-87c1-f5e713e9b3b9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "b = rc2models.PearsonIVParametrizationHPWModel()\n", + "par = b.guess(arr, x=x)\n", + "fitres = b.fit(arr, x=x, params=par)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fitres.eval(x=x))\n", + "\n", + "fitres.params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a7515a4-d0a8-4a5c-a834-0774fd202891", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "b = lmfit.models.Pearson4Model()\n", + "par = b.guess(arr, x=x)\n", + "fitres = b.fit(arr, x=x, params=par)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fitres.eval(x=x))\n", + "\n", + "fitres.params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "392c9f61-8413-49ae-a8fe-24ca8c8137e8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "x = np.linspace(0, 200, 1000, endpoint=False)\n", + "mod = lmfit.models.VoigtModel()\n", + "\n", + "par = mod.make_params()\n", + "par['amplitude'].set(value=1000)\n", + "par['center'].set(value=100)\n", + "par['sigma'].set(value=3)\n", + "par['gamma'].set(value=5, vary=True)\n", + "arr = mod.eval(params=par, x=x)\n", + "arr += np.random.uniform(-3,3, size=len(arr))\n", + "plt.plot(x,arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d85562f2-b98d-481e-8b6f-e8ab4aa9de53", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "a = rc2models_lm.VoigtNuModel()\n", + "par = a.guess(arr, x=x)\n", + "fitres = a.fit(arr, x=x, params=par)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fitres.eval(x=x))\n", + "fitres.params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17640903-7d74-4a9c-81be-97b31ccedc36", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "a = rc2models.VoigtAreaParametrizationNuModel()\n", + "par = a.guess(arr, x=x)\n", + "fitres = a.fit(arr, x=x, params=par)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fitres.eval(x=x))\n", + "fitres.params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbfbbb88-8005-4f29-8fec-ad31c107e4a1", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "a = lmfit.models.VoigtModel()\n", + "par = a.guess(arr, x=x)\n", + "par['gamma'].set(value=par['sigma'].value, vary=True)\n", + "fitres = a.fit(arr, x=x, params=par)\n", + "plt.plot(x, arr)\n", + "plt.plot(x, fitres.eval(x=x))\n", + "fitres.params" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ramanchada2", + "language": "python", + "name": "ramanchada2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/ramanchada2/fitting_functions/voigt_nu.py b/src/ramanchada2/fitting_functions/voigt_nu.py new file mode 100644 index 00000000..cf1e55e1 --- /dev/null +++ b/src/ramanchada2/fitting_functions/voigt_nu.py @@ -0,0 +1,143 @@ +import numpy as np +from scipy.special import erfc, voigt_profile + +DBL_EPSILON = np.finfo(float).eps +OneBySqrtLog4 = 0.84932180028801904272150283410289 # 1/sqrt(log(4)) +Log4 = 1.3862943611198906188344642429164 # Log(4) +Sqrt2 = 1.4142135623730950488016887242097 # Math.Sqrt(2) +Sqrt2Pi = 2.5066282746310005024157652848110 # Math.Sqrt(2*Math.Pi) +SqrtLog2 = 0.832554611157697756353165 # Math.Sqrt(Math.Log(2)) +Sqrt2PiByLog4 = 2.12893403886245235863054 # Math.Sqrt(2*Math.Pi/Math.Log(4)) +SqrtLog4 = 1.1774100225154746910 # Math.Sqrt(Math.Log(4)) +SqrtPi = 1.77245385090551602729817 # Math.Sqrt(Math.PI) +Log2 = 0.693147180559945309417232 # Math.Log(2) +Log2P32 = 0.577082881386139784459964 # Math.Pow(Log2, 1.5) +SqrtPiLog2 = 1.47566462663560588938882 # Math.Sqrt(Math.PI*Math.Log(2)) +C2_FWHM = 0.21669 # Approximation constant for FWHM of Voigt +C1_FWHM = 1 - np.sqrt(C2_FWHM) # Second constant for FWHM of Voigt + + +def voigt_nu_profile(x, amplitude=1.0, center=0.0, sigma=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) + ) + + +def voigt_nu_FWHM_approx(*, w, nu): + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + return 2*voigt_nu_HWHM_approx(sigma=sigma, gamma=gamma) + + +def voigt_nu_FWHM(*, w, nu): + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + return 2*voigt_nu_HWHM(sigma=sigma, gamma=gamma) + + +def voigt_nu_HWHM_approx(*, sigma, gamma): + """ + Gets the approximate half-width-half-maximum of a peak + in dependence on the original Voigt function parameters 'sigma' and 'gamma'. + The approximation has an accuracy of 0.0216%. + """ + C2 = 0.86676 # for relative error always < 0.000216 _and_ asymptotic correct behaviour + C2S = 2 - np.sqrt(C2) + return 0.5 * ( + C2S * gamma + + np.sqrt(C2 * gamma * gamma + 4 * 2 * np.log(2) * sigma * sigma) + ) + + +def voigt_nu_HWHM(*, sigma, gamma): + """ + Gets the half-width-half-maximum of a peak + in dependence on the original Voigt function parameters 'sigma' and 'gamma'. + """ + side = 0 + if sigma == 0 and gamma == 0: + return 0 + if np.isnan(sigma) or np.isnan(gamma): + return np.NaN + + # Reduce order of magnitude to prevent overflow + prefac = 1.0 + s = np.abs(sigma) + g = np.abs(gamma) + while s > 1e100 or g > 1e100: + prefac *= 1e30 + s /= 1e30 + g /= 1e30 + + # Increase order of magnitude to prevent underflow + while s < 1e-100 and g < 1e-100: + prefac /= 1e30 + s *= 1e30 + g *= 1e30 + + HM = voigt_profile(0.0, s, g) / 2 + + # Choose initial points a,b that bracket the expected root + c = voigt_nu_HWHM_approx(sigma=s, gamma=g) + a = c * 0.995 + b = c * 1.005 + del_a = voigt_profile(a, s, g) - HM + del_b = voigt_profile(b, s, g) - HM + + # Iteration using regula falsi (Illinois variant). + # Empirically, this takes <5 iterations to converge to FLT_EPSILON + # and <10 iterations to converge to DBL_EPSILON. + # We have never seen convergence worse than k = 15. + for k in range(30): + if np.abs(del_a - del_b) < 2 * DBL_EPSILON * HM: + return prefac * (a + b) / 2 + c = (b * del_a - a * del_b) / (del_a - del_b) + if np.abs(b - a) < 2 * DBL_EPSILON * np.abs(b + a): + return prefac * c + del_c = voigt_profile(c, s, g) - HM + + if del_b * del_c > 0: + b = c + del_b = del_c + if side < 0: + del_a /= 2 + side = -1 + elif del_a * del_c > 0: + a = c + del_a = del_c + if side > 0: + del_b /= 2 + side = 1 + else: + return prefac * c + + +def voigt_nu_height(*, amplitude, sigma, nu): + area = amplitude + w = sigma + + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + + if nu < (1 / 20.0): # approximation by series expansion is needed in the Lorentz limit + expErfcTermBySqrtNu = ( + ( + ( + ((4.21492418972454262863068) * nu + 0.781454843465470036843478) + * nu + + 0.269020594012510850443784 + ) + * nu + + 0.188831848731661377618477 + ) + * nu + ) + 0.677660751603104996618259 + else: # normal case + expTerm = np.exp(0.5 * np.square(gamma / sigma)) + erfcTerm = erfc(-gamma / (np.sqrt(2) * sigma)) + expErfcTermBySqrtNu = expTerm * (2 - erfcTerm) / np.sqrt(nu) + + bodyheight = expErfcTermBySqrtNu / (w * OneBySqrtLog4 * Sqrt2Pi) + height = area * bodyheight + return height diff --git a/src/ramanchada2/fitting_functions/voigtareaparametrizationnu.py b/src/ramanchada2/fitting_functions/voigtareaparametrizationnu.py new file mode 100644 index 00000000..f06523e2 --- /dev/null +++ b/src/ramanchada2/fitting_functions/voigtareaparametrizationnu.py @@ -0,0 +1,585 @@ +#!/usr/bin/env python3 + +from decimal import InvalidOperation +import numpy as np +from scipy.special import voigt_profile, wofz, erfc +from collections import namedtuple + + +DBL_EPSILON = np.finfo(float).eps +OneBySqrtLog4 = 0.84932180028801904272150283410289 # 1/sqrt(log(4)) +Log4 = 1.3862943611198906188344642429164 # Log(4) +Sqrt2 = 1.4142135623730950488016887242097 # Math.Sqrt(2) +Sqrt2Pi = 2.5066282746310005024157652848110 # Math.Sqrt(2*Math.Pi) +SqrtLog2 = 0.832554611157697756353165 # Math.Sqrt(Math.Log(2)) +Sqrt2PiByLog4 = 2.12893403886245235863054 # Math.Sqrt(2*Math.Pi/Math.Log(4)) +SqrtLog4 = 1.1774100225154746910 # Math.Sqrt(Math.Log(4)) +SqrtPi = 1.77245385090551602729817 # Math.Sqrt(Math.PI) +Log2 = 0.693147180559945309417232 # Math.Log(2) +Log2P32 = 0.577082881386139784459964 # Math.Pow(Log2, 1.5) +SqrtPiLog2 = 1.47566462663560588938882 # Math.Sqrt(Math.PI*Math.Log(2)) +C2_FWHM = 0.21669 # Approximation constant for FWHM of Voigt +C1_FWHM = 1 - np.sqrt(C2_FWHM) # Second constant for FWHM of Voigt + + +def SafeSqrt(x): + """Returns 0 if x is negative, otherwise the square root of x""" + return 0.0 if x < 0 else np.sqrt(x) + + +# Gets the maximal relative error of the VoigtHalfWidthHalfMaximumApproximation(sigma, gamma) function +VoigtHalfWidthHalfMaximumApproximationMaximalRelativeError = 0.000216 + +PositionAreaHeightFWHM = namedtuple( + "PositionAreaHeightFWHM", + [ + "Position", + "PositionStdDev", + "Area", + "AreaStdDev", + "Height", + "HeightStdDev", + "FWHM", + "FWHMStdDev", + ], +) + + +class VoigtAreaParametrizationNu: + """ + Voigt peak fitting function with multiple peaks and polynomial baseline + + A special parametrization is used for the Voigt function: + Each of the peaks has 4 parameters: area, pos, w, nu + Parameter 'area' designates the area under the peak + Parameter 'pos' designates the position of the peak (x-value at maximum) + Parameter 'w' is approximately (2%) the HWHM value of the peak + Parameter 'nu' (0..1) determines the Gauss'ness of the peak (0=Lorentzian, + 1=Gaussian) + + Note: the parameters 'sigma' and 'gamma' of the original Voigt function are + calculated as follows: + sigma = w * Sqrt(nu / log(4)) + gamma = w * (1 - nu) + """ + + numberOfTerms = 1 + orderOfBaselinePolynomial = -1 + + def __init__(self, numberOfTerms=1, orderOfBackgroundPolynomial=-1): + """ + Initializes an instance of the class with a given number of + peak terms, and the order of the baseline polynomial + """ + self.numberOfTerms = numberOfTerms + self.orderOfBaselinePolynomial = orderOfBackgroundPolynomial + + def WithNumberOfTerms(self, numberOfTerms): + """Returns a new instance of the class with the given number of peak terms.""" + return VoigtAreaParametrizationNu(numberOfTerms, self.orderOfBaselinePolynomial) + + def WithOrderOfBaselinePolynomial(self, orderOfBaselinePolynomial): + """Returns a new instance of the class with the given order of the baseline polynomial.""" + return VoigtAreaParametrizationNu(self.numberOfTerms, orderOfBaselinePolynomial) + + def GetNumberOfParametersPerPeak(self): + """Returns the number of parameters per peak term.""" + return 4 + + def GetParameterNamesForPeak(self, indexOfPeak): + """Returns the parameter names for a given peak index.""" + return [ + f"area{indexOfPeak}", + f"pos{indexOfPeak}", + f"w{indexOfPeak}", + f"nu{indexOfPeak}", + ] + + 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, w * np.sqrt(nu) * OneBySqrtLog4, w * (1 - nu) + ) + + def func(self, pars, x, data=None): + """Returns the y-values of the fitting function.""" + sum = np.zeros(len(x)) + for i in range(self.orderOfBaselinePolynomial, -1, -1): + sum *= x + sum += pars[f"b{i}"] + + for i in range(self.numberOfTerms): + a, xc, w, nu = ( + pars[f"area{i}"], + pars[f"pos{i}"], + pars[f"w{i}"], + pars[f"nu{i}"], + ) + sum += a * voigt_profile( + x - xc, w * np.sqrt(nu) * OneBySqrtLog4, w * (1 - nu) + ) + + if data is None: + return sum + return sum - data + + def dfunc(self, pars, x, data=None): + """Returns the derivatives of the fitting function w.r.t. the parameters.""" + result = [] + + # first the baseline derivatives + if self.orderOfBaselinePolynomial >= 0: + xn = np.ones(len(x)) + for i in range(self.orderOfBaselinePolynomial + 1): + if pars[f"b{i}"].vary: + result.append(np.copy(xn)) + xn *= x + + # second the derivatives of the peak terms + for i in range(self.numberOfTerms): + area, xc, w, nu = ( + pars[f"area{i}"], + pars[f"pos{i}"], + pars[f"w{i}"], + pars[f"nu{i}"], + ) + if not (area.vary or xc.vary or w.vary or nu.vary): + continue + + arg = x - xc + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + + dfdarea = np.full(len(x), np.NaN) + dfdpos = np.full(len(x), np.NaN) + dfdw = np.full(len(x), np.NaN) + dfdnu = np.full(len(x), np.NaN) + + if w > 0 and nu >= 0: + if nu < 1e-4: # approximately this is a Lorentzian + arg /= w + onePlusArg2 = 1 + arg * arg + body = 1 / (np.pi * w * onePlusArg2) + dfdarea = body + dfdpos = area * body * 2 * arg / (w * onePlusArg2) + dfdw = -area * body * (1 - arg * arg) / (w * onePlusArg2) + dfdnu = ( + area + * body + * (arg * arg * (3 - arg * arg * Log4) + Log4 - 1) + / (onePlusArg2 * onePlusArg2 * Log4) + ) + elif ( + nu <= 1 + ): # general case including nu==1 (which means gamma==0, i.e. pure Gaussian) + arg_gammaj = np.empty(len(arg), dtype=np.complex128) + arg_gammaj.real = arg + arg_gammaj.imag = gamma # arg_gammaj is now complex(arg, gamma) + z = arg_gammaj / (Sqrt2 * sigma) + wOfZ = wofz(z) + body = wOfZ / (Sqrt2Pi * sigma) + dbodydz = (Sqrt2 * 1.0j - z * wOfZ * Sqrt2Pi) / ( + np.pi * sigma + ) # Derivative of wOfZBySqrt2PiSigma w.r.t. z + dfdarea = np.real(body) # Derivative w.r.t. amplitude + dfdpos = ( + -area * np.real(dbodydz) / (Sqrt2 * sigma) + ) # Derivative w.r.t. position + dfdw = ( + -area / w * np.real(dbodydz * arg / (Sqrt2 * sigma) + body) + ) # Derivative w.r.t. w + argwonepnu = np.empty(len(arg), dtype=np.complex128) + argwonepnu.real = arg + argwonepnu.imag = w * (1 + nu) # complex(arg, w * (1 + nu)) + dfdnu = ( + -area + / (2 * nu) + * np.real((dbodydz * argwonepnu / (Sqrt2 * sigma) + body)) + ) # Derivative w.r.t. nu + + if area.vary: + result.append(dfdarea) + if xc.vary: + result.append(dfdpos) + if w.vary: + result.append(dfdw) + if nu.vary: + result.append(dfdnu) + + return np.array(result) + + def GetInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + height, position, fullWidth, relativeHeight=0.5 + ): + """ + Returns an array of initial parameters for peak fitting, + given the height, the position, the full width and the relative + height at which the full width was measured. + """ + useLorenzLimit = False + + if not (relativeHeight > 0 and relativeHeight < 1): + raise ValueError("RelativeHeight should be in the open interval (0, 1)") + + if useLorenzLimit: + # we calculate at the Lorentz limit (nu==0) + w = 0.5 * fullWidth * np.sqrt(relativeHeight / (1 - relativeHeight)) + area = height * w * np.pi + return [area, position, w, 0] + else: + # use Gaussian limit + w = 0.5 * fullWidth * SqrtLog2 / np.sqrt(-np.log(relativeHeight)) + area = height * w * Sqrt2PiByLog4 + return [area, position, w, 1] + + def AddInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + self, paras, indexOfPeak, height, position, fullWidth, relativeHeight=0.5 + ): + """ + Add to an existing 'Parameters' instance the initial parameters for a peak. + + Parameters + ---------- + paras: Parameters + Existing instance of Parameters, to which the new parameters will be added. + indexOfPeak: + The index of the peak whose parameters are added. + height: + Approximate height of the peak. + pos: + Approximate position of the peak + fullWidth: + Approximate full width of the peak + relativeHeight: + Relative height, at which full width was measured (usually 0.5) + """ + apwnu = VoigtAreaParametrizationNu.GetInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + height, position, fullWidth, relativeHeight + ) + paras.add(f"area{indexOfPeak}", apwnu[0]) + paras.add(f"pos{indexOfPeak}", apwnu[1]) + paras.add(f"w{indexOfPeak}", apwnu[2]) + paras.add(f"nu{indexOfPeak}", apwnu[3]) + + def AddInitialParametersForBaselinePolynomial(self, paras): + """Add the required parameters for the baseline polynomial""" + for i in range(self.orderOfBaselinePolynomial + 1): + paras.add(f"b{i}", 0.0) + + def VoigtHalfWidthHalfMaximumOfSigmaGammaApproximation(sigma, gamma): + """ + Gets the approximate half-width-half-maximum of a peak + in dependence on the original Voigt function parameters 'sigma' and 'gamma'. + The approximation has an accuracy of 0.0216%. + """ + C2 = 0.86676 # for relative error always < 0.000216 _and_ asymptotic correct behaviour + C2S = 2 - np.sqrt(C2) + return 0.5 * ( + C2S * gamma + + np.sqrt(C2 * gamma * gamma + 4 * 2 * np.log(2) * sigma * sigma) + ) + + def VoigtHalfWidthHalfMaximumOfWNuApproximation(w, nu): + """ + Gets the approximate half-width-half-maximum of a peak + in dependence on the parameters 'w' and 'nu'. + The approximation has an accuracy of 0.0216%. + """ + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + return VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGammaApproximation( + sigma, gamma + ) + + def VoigtHalfWidthHalfMaximumOfSigmaGamma(sigma, gamma): + """ + Gets the half-width-half-maximum of a peak + in dependence on the original Voigt function parameters 'sigma' and 'gamma'. + """ + side = 0 + if sigma == 0 and gamma == 0: + return 0 + if np.isnan(sigma) or np.isnan(gamma): + return np.NaN + + # Reduce order of magnitude to prevent overflow + prefac = 1.0 + s = np.abs(sigma) + g = np.abs(gamma) + while s > 1e100 or g > 1e100: + prefac *= 1e30 + s /= 1e30 + g /= 1e30 + + # Increase order of magnitude to prevent underflow + while s < 1e-100 and g < 1e-100: + prefac /= 1e30 + s *= 1e30 + g *= 1e30 + + HM = voigt_profile(0.0, s, g) / 2 + + # Choose initial points a,b that bracket the expected root + c = VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGammaApproximation( + s, g + ) + a = c * 0.995 + b = c * 1.005 + del_a = voigt_profile(a, s, g) - HM + del_b = voigt_profile(b, s, g) - HM + + # Iteration using regula falsi (Illinois variant). + # Empirically, this takes <5 iterations to converge to FLT_EPSILON + # and <10 iterations to converge to DBL_EPSILON. + # We have never seen convergence worse than k = 15. + for k in range(30): + if np.abs(del_a - del_b) < 2 * DBL_EPSILON * HM: + return prefac * (a + b) / 2 + c = (b * del_a - a * del_b) / (del_a - del_b) + if np.abs(b - a) < 2 * DBL_EPSILON * np.abs(b + a): + return prefac * c + del_c = voigt_profile(c, s, g) - HM + + if del_b * del_c > 0: + b = c + del_b = del_c + if side < 0: + del_a /= 2 + side = -1 + elif del_a * del_c > 0: + a = c + del_a = del_c + if side > 0: + del_b /= 2 + side = 1 + else: + return prefac * c + + raise InvalidOperation # One should never arrive here + + def VoigtHalfWidthHalfMaximumOfWNu(w, nu): + """ + Gets the half-width-half-maximum of a peak + in dependence on the parameters 'w' and 'nu'. + """ + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + return VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGamma( + sigma, gamma + ) + + def GetPositionAreaHeightFWHMFromSinglePeakParameters( + self, parameters, indexOfPeak, 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 + ---------- + parameters: Parameters + Existing instance of Parameters, which contains the result of the fit. + indexOfPeak: int + The index of the peak into consideration. + cv: np.array + Covariance matrix of the fit. + """ + area = parameters[f"area{indexOfPeak}"] + pos = parameters[f"pos{indexOfPeak}"] + w = parameters[f"w{indexOfPeak}"] + nu = parameters[f"nu{indexOfPeak}"] + + if cv is not None: + cv = cv[ + indexOfPeak * 4: (indexOfPeak + 1) * 4, + indexOfPeak * 4: (indexOfPeak + 1) * 4, + ] + return VoigtAreaParametrizationNu.GetPositionAreaHeightFWHMFromPeakParameters(area, pos, w, nu, cv) + + def GetPositionAreaHeightFWHMFromPeakParameters( + area, pos, w, nu, 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 + ---------- + parameters: Parameters + Existing instance of Parameters, which contains the result of the fit. + indexOfPeak: int + The index of the peak into consideration. + cv: np.array + Covariance matrix of the fit. + """ + sigma = w * np.sqrt(nu) * OneBySqrtLog4 + gamma = w * (1 - nu) + + if cv is not None: + areaStdDev = SafeSqrt(cv[0, 0]) + posStdDev = SafeSqrt(cv[1, 1]) + else: + areaStdDev = 0.0 + posStdDev = 0.0 + + height = np.NaN + fwhm = np.NaN + heightStdDev = 0.0 + fwhmStdDev = 0.0 + + # expErfcTerm is normally: Math.Exp(0.5 * RMath.Pow2(gamma / sigma)) * Erfc(gamma / (Math.Sqrt(2) * sigma)) + # but for large gamma/sigma, we need a approximation, because the exp term becomes too large + expErfcTermBySqrtNu = np.NaN + + # for gamma > 20*sigma we need an approximation of the expErfcTerm, + # since the expTerm will get too large and the Erfc term too small + # we use a series expansion + + if nu < ( + 1 / 20.0 + ): # approximation by series expansion is needed in the Lorentz limit + expErfcTermBySqrtNu = ( + ( + ( + ((4.21492418972454262863068) * nu + 0.781454843465470036843478) + * nu + + 0.269020594012510850443784 + ) + * nu + + 0.188831848731661377618477 + ) + * nu + ) + 0.677660751603104996618259 + else: # normal case + expTerm = np.exp(0.5 * np.square(gamma / sigma)) + erfcTerm = erfc(-gamma / (np.sqrt(2) * sigma)) + expErfcTermBySqrtNu = expTerm * (2 - erfcTerm) / np.sqrt(nu) + + bodyheight = expErfcTermBySqrtNu / (w * OneBySqrtLog4 * Sqrt2Pi) + height = area * bodyheight + + # fwhmSqrtTerm = np.sqrt(C2_FWHM * gamma * gamma + 8 * sigma * sigma * np.log(2)); + # fwhm = C1_FWHM * gamma + fwhmSqrtTerm; // this is the approximation formula + fwhm = 2 * VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGamma( + sigma, gamma + ) # we use the exact formula for fwhm + + if cv is not None: + dHeightByDArea = bodyheight + dHeightByDW = -bodyheight * area / w + dHeightByDNu = np.NaN + + if nu < (1 / 20.0): + # Series expansion for Lorentzian limit + dHeightByDNu = ( + ( + ( + ( + ( + ( + 88758.4100339208444038841 * nu + - 7108.49617246574864439307 + ) + * nu + + 641.286392761620402559762 + ) + * nu + - 66.1018151520036002556775 + ) + * nu + + 7.91931382144031445585169 + ) + * nu + - 1.10119171735779477287440 + ) + * nu + + 0.252727974753276908699454 + ) * nu + 0.0886978390521480859157182 + else: + # General case + dHeightByDNu = ( + (1 + nu) * Log2 + - (expErfcTermBySqrtNu) + * ((1 - nu * nu) * SqrtPi * Log2P32 + 0.5 * nu * SqrtPiLog2) + ) / (nu * nu * np.pi) + + dHeightByDNu *= area / w + + heightStdDev = SafeSqrt( + cv[0, 0] * np.square(dHeightByDArea) + + cv[2, 2] * np.square(dHeightByDW) + + cv[3, 3] * np.square(dHeightByDNu) + + dHeightByDNu * (cv[0, 3] + cv[3, 0]) * dHeightByDArea + + dHeightByDW * (cv[0, 2] + cv[2, 0]) * dHeightByDArea + + dHeightByDNu * (cv[2, 3] + cv[3, 2]) * dHeightByDW + ) + + dFwhmByDW = 2 * ( + 1 + + np.sqrt(C2_FWHM) * (-1 + nu) + - nu + + np.sqrt(C2_FWHM * np.square(1 - nu) + nu) + ) + + dFwhmByDNu = ( + 2 + * w + * ( + -1 + + np.sqrt(C2_FWHM) + + (1 + 2 * C2_FWHM * (-1 + nu)) + / (2 * np.sqrt(C2_FWHM * np.square(1 - nu) + nu)) + ) + ) + + fwhmStdDev = SafeSqrt( + cv[3, 3] * np.square(dFwhmByDNu) + + dFwhmByDW * ((cv[2, 3] + cv[3, 2]) * dFwhmByDNu + cv[2, 2] * dFwhmByDW) + ) + + return PositionAreaHeightFWHM( + pos, posStdDev, area, areaStdDev, height, heightStdDev, fwhm, fwhmStdDev + ) + + def SetParameterBoundariesForPositivePeaks( + self, + params, + minimalPosition=None, + maximalPosition=None, + minimalFWHM=None, + maximalFWHM=None, + ): + """ + Amends the parameters in an existing Parameters instance + with fit boundaries. + + Parameters + ---------- + params: Parameters + Existing instance of Parameters, which already contain all parameters. + minimalPosition: float (or None) + The minimal allowed position of a peak. + maximalPosition: float (or None) + The maximal allowed position of a peak. + minimalFWHM: float (or None) + The minimal allowed FWHM value of a peak. + maximalFWHM: float (or None) + The maximal allowed FWHM value of a peak. + + Note: the minimal height is set to 0 in order to disallow negative peaks. + the maximal height of a peak is not limited + """ + for i in range(self.numberOfTerms): + params[f"area{i}"].min = 0 # minimal area is 0 + if minimalPosition is not None: + params[f"pos{i}"].min = minimalPosition + if maximalPosition is not None: + params[f"pos{i}"].max = maximalPosition + if minimalFWHM is not None: + params[f"w{i}"].min = minimalFWHM / 2.0 + else: + params[f"w{i}"].min = np.sqrt(np.sqrt(DBL_EPSILON)) + if maximalFWHM is not None: + params[f"w{i}"].max = maximalFWHM / 2.0 + params[f"nu{i}"].min = 0 + params[f"nu{i}"].max = 1 diff --git a/src/ramanchada2/spectrum/peaks/find_peaks_ByIncrementalPeakAddition.py b/src/ramanchada2/spectrum/peaks/find_peaks_ByIncrementalPeakAddition.py new file mode 100644 index 00000000..1ecbd7e9 --- /dev/null +++ b/src/ramanchada2/spectrum/peaks/find_peaks_ByIncrementalPeakAddition.py @@ -0,0 +1,207 @@ +import numpy as np +from scipy import signal +import lmfit + + +def IsInRange(x, idx): + return idx >= 0 and idx <= len(x) - 1 + + +def InterpolateLinear(index, array, extendToSides=False): + if not (index >= 0): + if extendToSides: + return array[0] + else: + raise AttributeError + elif not (index <= len(array) - 1): + if extendToSides: + return array[-1] + else: + raise AttributeError + elif index == len(array) - 1: + return array[-1] + + idx = int(np.floor(index)) + fraction = index - idx + if fraction == 0: + return array[idx] # return left value + elif fraction == 1: + return array[idx + 1] # return right value + else: + return (1 - fraction) * array[idx] + fraction * array[idx + 1] + + +def GetWidthValue(x, leftIdx, middleIdx, rightIdx): + if IsInRange(x, leftIdx) and IsInRange(x, rightIdx): + return InterpolateLinear(rightIdx, x) - InterpolateLinear(leftIdx, x) + elif IsInRange(x, middleIdx) and IsInRange(x, leftIdx): + return 2 * (InterpolateLinear(middleIdx, x) - InterpolateLinear(leftIdx, x)) + elif IsInRange(x, rightIdx) and IsInRange(x, middleIdx): + return 2 * (InterpolateLinear(rightIdx, x) - InterpolateLinear(middleIdx, x)) + else: + # Try to interpolate over the full range of x + return (x[-1] - x[0]) * (rightIdx - leftIdx) / (len(x) - 1.0) + + +class FindPeaksByIncrementalPeakAddition: + + # Order of the baseline polynomial (-1=no baseline, 0=constant, 1=linear, etc.) + orderOfBaselinePolynomial = -1 + + # Maximum number of peak that will be fitted + maximumNumberOfPeaks = 50 + + minimalRelativeHeight = 2.5e-3 + + minimalSignalToNoiseRatio = 8 + + fitWidthScalingFactor = 2 + + prunePeaksSumChiSquareFactor = 0.1 + + minimalFWHMValue = None + + isMinimalFWHMValueInXUnits = True + + def Execute(self, x, y, fitFuncTemplate): + # estimate the properties of the x and y array + + minimalX = np.min(x) + maximalX = np.max(x) + minimalY = np.min(y) + maximalY = np.max(y) + minIncrement = np.min(np.abs(x[1:] - x[:-1])) + spanX = maximalX - minimalX + spanY = maximalY - minimalY + + # estimate the noise level + noiseArray = np.abs(y[1:-1] - 0.5 * y[2:] - 0.5 * y[:-2]) + np.sort(noiseArray) + noiseLevel = ( + noiseArray[int(len(noiseArray) / 2)] * 1.22 + ) # take the 50% percentile ans noise level + + fitFunction = fitFuncTemplate.WithNumberOfTerms( + 0 + ).WithOrderOfBaselinePolynomial(self.orderOfBaselinePolynomial) + + prohibitedPeaks = {} + + yRest = np.copy(y) + + minimalPeakHeightForSearching = max( + np.abs(self.minimalRelativeHeight * spanY), + np.abs(noiseLevel * self.minimalSignalToNoiseRatio), + ) + + minimalFWHM = ( + minIncrement + if self.minimalFWHMValue is None + else ( + self.minimalFWHMValue + if self.isMinimalFWHMValueInXUnits and self.minimalFWHMValue > 0 + else max(1, self.minimalFWHMValue) * minIncrement + ) + ) + maximalFWHM = spanX + + parameters = lmfit.Parameters() + + # before adding peaks, fit the baseline + fitFunction.AddInitialParametersForBaselinePolynomial(parameters) + min2 = lmfit.Minimizer( + fitFunction.func, parameters, fcn_args=(x,), fcn_kws={"data": y} + ) + out2 = min2.leastsq(Dfun=fitFunction.dfunc, col_deriv=1) + parameters = out2.params + yRest = y - fitFunction.func(parameters, x) + + for numberOfTerms in range(1, self.maximumNumberOfPeaks + 1): + (peakIndices, peakProperties) = signal.find_peaks( + yRest, height=minimalPeakHeightForSearching, width=0.0, rel_height=0.5 + ) + if len(peakIndices) == 0: + break + + peakHeights = peakProperties["peak_heights"] + peakWidths = signal.peak_widths(yRest, peakIndices, 0.5) + # Sort the three arrays by peak height + sortKey = np.argsort(peakHeights) + peakHeights = peakHeights[sortKey] + peakWidths = peakWidths[0][sortKey] + peakIndices = peakIndices[sortKey] + + thisPeakIndex = -1 + for j in range(len(peakHeights) - 1, -1, -1): + if ( + not (peakIndices[j] in prohibitedPeaks) + or prohibitedPeaks[peakIndices[j]] == 0 + ): + thisPeakIndex = j + break + + if thisPeakIndex < 0 or peakHeights[thisPeakIndex] == 0: + break # no more peaks to fit + + thisPeakPosition = peakIndices[thisPeakIndex] + + yMax = y[thisPeakPosition] + fwhm = 0.5 * np.abs( + GetWidthValue( + x, + thisPeakPosition - 0.5 * peakWidths[thisPeakIndex], + thisPeakPosition, + thisPeakPosition + 0.5 * peakWidths[thisPeakIndex], + ) + ) + fitFunction = fitFunction.WithNumberOfTerms(numberOfTerms) + # at first, set the previous parameters to fixed + for para in parameters: + parameters[para].vary = False + fitFunction.AddInitialParametersFromHeightPositionAndWidthAtRelativeHeight( + parameters, numberOfTerms - 1, yMax, x[thisPeakPosition], fwhm, 0.5 + ) + fitFunction.SetParameterBoundariesForPositivePeaks( + parameters, minimalX, maximalX, minimalFWHM, maximalFWHM + ) + min1 = lmfit.Minimizer( + fitFunction.func, parameters, fcn_args=(x,), fcn_kws={"data": y} + ) + out1 = min1.leastsq(Dfun=fitFunction.dfunc, col_deriv=1) + parameters = out1.params + + # if the amplitude of the new peak is zero, then this peak has a peculiar shape, + # we set this peak on a list of prohibited peaks + if ( + parameters[fitFunction.GetParameterNamesForPeak(numberOfTerms - 1)[0]] + == 0 + ): + if thisPeakPosition in prohibitedPeaks: + numberOfHits = prohibitedPeaks[thisPeakPosition] + prohibitedPeaks[thisPeakPosition] = 2 * (numberOfHits + 1) + else: + prohibitedPeaks[thisPeakPosition] = 2 + + numberOfTerms -= 1 + continue + + # relax the prohibition status of the keys by one + for key in prohibitedPeaks: + prohibitedPeaks[key] = max(0, prohibitedPeaks[key] - 1) + + # at second, set all parameters to vary + for para in parameters: + parameters[para].vary = True + + # now fit again, with all parameters to be able to vary + min2 = lmfit.Minimizer( + fitFunction.func, parameters, fcn_args=(x,), fcn_kws={"data": y} + ) + out2 = min2.leastsq(Dfun=fitFunction.dfunc, col_deriv=1) + parameters = out2.params + + # calculate the yRest as the difference between y and the fit + yRest = y - fitFunction.func(parameters, x) + # the parameters now contains all peaks + # lets see if some of them can be pruned, because they do not contribute much to the improvement of the fit + return out2 diff --git a/tests/fitting_functions/test_pearsonivamplitudeparametrizationhpw.py b/tests/fitting_functions/test_pearsonivamplitudeparametrizationhpw.py new file mode 100644 index 00000000..025bddbb --- /dev/null +++ b/tests/fitting_functions/test_pearsonivamplitudeparametrizationhpw.py @@ -0,0 +1,479 @@ +from ramanchada2.fitting_functions.pearsonivamplitudeparametrizationhpw import ( + PearsonIVAmplitudeParametrizationHPW, + ) +from ramanchada2.fitting_functions.models import PearsonIVParametrizationHPWModel +import lmfit +import numpy as np + + +def almostequal(desired, actual, abstol, reltol): + delta = np.abs(abstol) + np.abs(reltol * desired) + return np.abs(desired - actual) <= delta + + +def test_FunctionValuesBasics(): + # amp, pos, w, m, v + testData = [ + (3, 7, 5, 11, 13), + (3, 7, 1 / 5.0, 1 / 1000.0, 1000), + (3, 7, 1 / 5.0, 1 / 1000.0, -1000), + (3, 7, 1 / 5.0, 1 / 1000.0, 0), + (3, 7, 1 / 5.0, 1000.0, 1000), + (3, 7, 1 / 5.0, 1000.0, -1000), + (3, 7, 1 / 5.0, 1000.0, 0), + ] + + func = PearsonIVAmplitudeParametrizationHPW(1, -1) + for (amp, pos, w, m, v) in testData: + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + + # Y at x=pos should be equal to amp + X = np.array([pos]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + assert almostequal(amp, y, 1e-10, 1e-10) + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + # Y at x= pos - w /1000 should be less than amp + X = np.array([pos - w / 1000.0]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + assert amp > y + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + # Y at x= pos + w /1000 should be less than amp + X = np.array([pos + w / 1000.0]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + assert amp > y + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + # Y at x= pos -3 w should be less than amp/2 + X = np.array([pos - 3 * w]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + assert amp / 2 > y + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + # Y at x= pos + 3 w should be less than amp/2 + X = np.array([pos + 3 * w]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + assert amp / 2 > y + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + +def test_FunctionValuesInDependenceOnX(): + testData = [ + (-1e9, 3, 7, 5, 11, 13, 2.7867315365921926531e-160), + (-19, 3, 7, 5, 11, 13, 0.0071554075200805234436), + (-1, 3, 7, 5, 11, 13, 0.94609302487197904495), + (0, 3, 7, 5, 11, 13, 1.1845028216303614895), + (6, 3, 7, 5, 11, 13, 2.9207184451185306856), + (7, 3, 7, 5, 11, 13, 3.0000000000000000000), + (8, 3, 7, 5, 11, 13, 2.9093699544835131623), + (23, 3, 7, 5, 11, 13, 5.7924047779440271306e-25), + (28, 3, 7, 5, 11, 13, 7.9007462868389145363e-260), + ] + + func = PearsonIVAmplitudeParametrizationHPW(1, -1) + for (x, amp, pos, w, m, v, yexpected) in testData: + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + + X = np.array([x]) + y = PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(X[0], amp, pos, w, m, v) + if yexpected < 1e-100: + assert almostequal(yexpected, y, 0, 1e-7) + else: + assert almostequal(yexpected, y, 0, 1e-13) + + Y = func.func(p, X) + assert almostequal(y, Y[0], 1e-15, 1e-15) + + +def test_DerivativesWrtParameters(): + ff = PearsonIVAmplitudeParametrizationHPW(1, -1) + + # General case + amp = 17 + pos = 7 + w = 3 + m = 5 + v = 7 + + expectedFunctionValue = 10.456864176583425249219583024250 + expectedDerivativeWrtAmplitude = 0.61510965744608383818938723672062 + expectedDerivativeWrtPosition = 6.2023944308469781499266683603126 + expectedDerivativeWrtW = 4.1349296205646520999511122402084 + expectedDerivativeWrtM = 0.31213582841409496793210643128803 + expectedDerivativeWrtV = -0.0080620681798224796387712054749779 + + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + X = np.array([9]) + FV = ff.func(p, X) + DY = ff.dfunc(p, X) + + assert almostequal(expectedFunctionValue, FV[0], 0, 1e-12) + assert almostequal(expectedDerivativeWrtAmplitude, DY[0, 0], 0, 1e-12) + assert almostequal(expectedDerivativeWrtPosition, DY[1, 0], 0, 1e-12) + assert almostequal(expectedDerivativeWrtW, DY[2, 0], 0, 1e-12) + assert almostequal(expectedDerivativeWrtM, DY[3, 0], 0, 1e-12) + assert almostequal(expectedDerivativeWrtV, DY[4, 0], 0, 1e-11) + + +def test_AreaFwhm(): + testData = [ + (1, 7, 5, 11, 13), + (3, 7, 5, 1000, -1000), + (3, 7, 5, 1000, 0), + (3, 7, 5, 1000, 1000), + (3, 7, 5, 513 / 1024.0, -1000), + (3, 7, 5, 513 / 1024.0, 0), + (3, 7, 5, 513 / 1024.0, 1000), + ] + + expectedAreas = [ + 11.017024491957193080, + 31.945780266160134464, + 31.940456418616702875, + 31.945780266160134464, + 12084.703859367018981, + 8896.1196272036681069, + 12084.703859367018981, + ] + + expectedFwhm = [ + 10.069306246733044125, + 10.000770127169960848, + 10.000000000000000000, + 10.000770127169960848, + 11.364892411874591570, + 10.000000000000000000, + 11.364892411874591570, + ] + + func = PearsonIVAmplitudeParametrizationHPW(1, -1) + + for i in range(len(testData)): + (amp, pos, w, m, v) = testData[i] + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + ( + position, + posErr, + area, + areaErr, + height, + heightErr, + fwhm, + fwhmErr, + ) = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0) + + assert almostequal(pos, position, 0, 1e-14) + assert almostequal(amp, height, 0, 1e-14) + assert almostequal(expectedAreas[i], area, 0, 1e-8) + assert almostequal(expectedFwhm[i], fwhm, 0, 1e-8) + + # with the first test set, we test additionally + (amp, pos, w, m, v) = testData[0] + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + ( + position, + posErr, + area, + areaErr, + height, + heightErr, + fwhm, + fwhmErr, + ) = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0) + + # (i) when we vary the amplitude, area and height should change accordingly, pos and fwhm stay constant + for newAmp in [-11, 0, 2, 17, 33, 1023, -11]: + p = lmfit.Parameters() + p.add("a0", newAmp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + ( + newposition, + _, + newarea, + _, + newheight, + _, + newfwhm, + _, + ) = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0) + assert almostequal(position, newposition, 0, 1e-14) + assert almostequal(height * newAmp, newheight, 0, 1e-14) + assert almostequal(area * newAmp, newarea, 0, 1e-14) + assert almostequal(fwhm, newfwhm, 0, 1e-14) + + # (ii) when we vary the position, position should change accordingly, everything else is the same + for newPos in [-11, 0, 2, 17, 33, 1023, -11]: + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", newPos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + ( + newposition, + _, + newarea, + _, + newheight, + _, + newfwhm, + _, + ) = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0) + assert almostequal(newPos, newposition, 0, 1e-14) + assert almostequal(height, newheight, 0, 1e-14) + assert almostequal(area, newarea, 0, 1e-14) + assert almostequal(fwhm, newfwhm, 0, 1e-14) + + # (iii) when we vary the width, area and fwhm should change accordingly, everything else is the same + for newW in [1 / 1000.0, 2, 17, 33, 1023]: + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", newW) + p.add("m0", m) + p.add("v0", v) + ( + newposition, + _, + newarea, + _, + newheight, + _, + newfwhm, + _, + ) = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0) + assert almostequal(position, newposition, 0, 1e-14) + assert almostequal(height, newheight, 0, 1e-14) + assert almostequal(area * (newW / w), newarea, 0, 1e-14) + assert almostequal(fwhm * (newW / w), newfwhm, 0, 1e-14) + + +def test_FwhmApproximation(): + claimedMaxApproximationError = 0.19 # 19% error + w = 7 + for idx_v in range(-1, 16): + v = 0 if idx_v < 0 else np.power(10, (idx_v - 7.5) / 2.5) + for idx_m in range(0, 16): + m = np.power(10, (idx_m - 7.5) / 2.5) + fwhmP = PearsonIVAmplitudeParametrizationHPW.GetFWHM(w, m, v) + fwhmN = PearsonIVAmplitudeParametrizationHPW.GetFWHM(w, m, -v) + + fwhmApproxP = PearsonIVAmplitudeParametrizationHPW.GetFWHMApproximation( + w, m, v + ) + fwhmApproxN = PearsonIVAmplitudeParametrizationHPW.GetFWHMApproximation( + w, m, -v + ) + + # FWHM should be independent of whether v is positive or negative + assert almostequal(fwhmP, fwhmN, 1e-8, 1e-8) + + assert almostequal(fwhmP, fwhmApproxP, 0, claimedMaxApproximationError) + assert almostequal(fwhmN, fwhmApproxN, 0, claimedMaxApproximationError) + + +def test_DerivativesOfAreaFwhm(): + # see PearsonIV (amplitudeModified).nb + amp = 3 + pos = 7 + w = 5 + m = 11 + v = 13 + + ymaxDerivs = [1, 0, 0, 0, 0] + xmaxDerivs = [0, 1, 0, 0, 0] + areaDerivs = [ + 11.0170244919571930798, + 0, + 6.6102146951743158479, + -0.106359067713399883296, + 0.00045052750734765666701, + ] + fwhmDerivs = [ + 0, + 0, + 2.0138612493466088250930420306, + -0.0062732720000624223737773164, + 0.0000628521667333096080982645, + ] + + p = lmfit.Parameters() + p.add("a0", amp) + p.add("pos0", pos) + p.add("w0", w) + p.add("m0", m) + p.add("v0", v) + + func = PearsonIVAmplitudeParametrizationHPW(1, -1) + + for i in range(0, 5): + cov = np.zeros((5, 5)) + cov[i, i] = 1 + result = func.GetPositionAreaHeightFWHMFromSinglePeakParameters(p, 0, cov) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-4) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-4) + + +def test_fit_without_derivatives_oneterm(): + ff = PearsonIVAmplitudeParametrizationHPW(1, -1) + 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) + params = lmfit.Parameters() + params.add("a0", 10) + params.add("pos0", 49) + params.add("w0", 6) + params.add("m0", 1) + params.add("v0", 0) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq() + np.testing.assert_almost_equal(out1.params["a0"], a, 12) + np.testing.assert_almost_equal(out1.params["pos0"], pos, 12) + np.testing.assert_almost_equal(out1.params["w0"], w, 12) + np.testing.assert_almost_equal(out1.params["m0"], m, 12) + np.testing.assert_almost_equal(out1.params["v0"], v, 12) + + +def test_fit_with_derivatives_oneterm(): + ff = PearsonIVAmplitudeParametrizationHPW(1, -1) + 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) + params = lmfit.Parameters() + # important: the parameters must be added in exactly the order + # as you can see here (area0, pos0, w0, nu0, ... areaN, posN, wN, nuN, b0, b1, ... bM) + params.add("a0", 10) + params.add("pos0", 49) + params.add("w0", 6) + params.add("m0", 1) + params.add("v0", 0) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq(Dfun=ff.dfunc, col_deriv=1) + np.testing.assert_almost_equal(out1.params["a0"], a, 12) + np.testing.assert_almost_equal(out1.params["pos0"], pos, 12) + np.testing.assert_almost_equal(out1.params["w0"], w, 12) + np.testing.assert_almost_equal(out1.params["m0"], m, 12) + np.testing.assert_almost_equal(out1.params["v0"], v, 12) + + +def test_fit_with_derivatives_twoterms_linearbaseline(): + ff = PearsonIVAmplitudeParametrizationHPW(2, 1) + x = np.linspace(0, 99, 100) + a0 = 7 + pos0 = 25 + w0 = 5 + m0 = 16 / 11.0 + v0 = 1 / 4.0 + a1 = 6 + pos1 = 75 + w1 = 6 + m1 = 21 / 11.0 + v1 = 1 / 8.0 + b0 = 100 + b1 = 0.25 + data = ( + PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(x, a0, pos0, w0, m0, v0) + + PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(x, a1, pos1, w1, m1, v1) + + b0 + + b1 * x + ) + + # important: the parameters must be added in exactly the order + # as you can see here (b0, b1, ... bM, area0, pos0, w0, nu0, ... areaN, posN, wN, nuN, ) + params = lmfit.Parameters() + params.add("b0", 99) + params.add("b1", 0.24) + params.add("a0", 10) + params.add("pos0", 24) + params.add("w0", 6) + params.add("m0", 1) + params.add("v0", 0) + params.add("a1", 10) + params.add("pos1", 74) + params.add("w1", 7) + params.add("m1", 1) + params.add("v1", 0) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq(Dfun=ff.dfunc, col_deriv=1) + cv = out1.covar + print(np.shape(cv)) + + np.testing.assert_almost_equal(out1.params["b0"], b0, 12) + np.testing.assert_almost_equal(out1.params["b1"], b1, 12) + np.testing.assert_almost_equal(out1.params["a0"], a0, 12) + np.testing.assert_almost_equal(out1.params["pos0"], pos0, 12) + np.testing.assert_almost_equal(out1.params["w0"], w0, 12) + np.testing.assert_almost_equal(out1.params["m0"], m0, 12) + np.testing.assert_almost_equal(out1.params["v0"], v0, 12) + np.testing.assert_almost_equal(out1.params["a1"], a1, 12) + np.testing.assert_almost_equal(out1.params["pos1"], pos1, 12) + 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["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) + np.testing.assert_almost_equal(result.params["fwhm"], 10.02952431724565395, 7) diff --git a/tests/fitting_functions/test_voigtareaparametrizationnu.py b/tests/fitting_functions/test_voigtareaparametrizationnu.py new file mode 100644 index 00000000..d26b6aac --- /dev/null +++ b/tests/fitting_functions/test_voigtareaparametrizationnu.py @@ -0,0 +1,819 @@ +from ramanchada2.fitting_functions import voigtareaparametrizationnu +from ramanchada2.fitting_functions.models import VoigtAreaParametrizationNuModel +from ramanchada2.fitting_functions.voigtareaparametrizationnu import ( + VoigtAreaParametrizationNu, +) +import lmfit +import numpy as np +from scipy.special import voigt_profile + +# The voigt HWHM with sigma=1, and gamma = (1+d)/(1-d), d = [-1, 1/64, 1) +voigt_hwhm_sigmaone_dvar = [ + 1.177410022515474691, + 1.1816092691813030468, + 1.1858873739520897034, + 1.1902465142050570863, + 1.194688946189857397, + 1.1992170085518593889, + 1.2038331260421231307, + 1.2085398134254873464, + 1.2133396795989858683, + 1.2182354319336629671, + 1.2232298808537770661, + 1.2283259446683743323, + 1.2335266546712841194, + 1.2388351605267440826, + 1.2442547359591114981, + 1.2497887847664671425, + 1.2554408471793780475, + 1.2612146065876654677, + 1.2671138966597353836, + 1.2731427088808827742, + 1.2793052005389909397, + 1.2856057031882278137, + 1.2920487316237085011, + 1.2986389934026648102, + 1.3053813989504578017, + 1.3122810722928098214, + 1.3193433624589418804, + 1.3265738556039068772, + 1.3339783879023381003, + 1.3415630592701179674, + 1.3493342479751497996, + 1.3572986262035252986, + 1.3654631766529664197, + 1.3738352102315315844, + 1.3824223849462672922, + 1.3912327260738180846, + 1.4002746477130484284, + 1.4095569758285553194, + 1.4190889729036451368, + 1.4288803643320075071, + 1.4389413666890491483, + 1.4492827180367713744, + 1.4599157104303203962, + 1.4708522248100608452, + 1.4821047684803902488, + 1.4936865153957175192, + 1.5056113494952888279, + 1.5178939113521050052, + 1.5305496484273140468, + 1.5435948692504962536, + 1.5570468018785470874, + 1.5709236570218130826, + 1.5852446962662159034, + 1.6000303058648426117, + 1.6153020766224977944, + 1.6310828904527072039, + 1.6473970142494402595, + 1.6642702017863107582, + 1.681729804435296653, + 1.69980489158633592, + 1.7185263817499517539, + 1.7379271854390185723, + 1.7580423610548628274, + 1.7789092851493945634, + 1.8005678386015786904, + 1.823060610436444369, + 1.8464331212317292609, + 1.8707340683055766275, + 1.8960155951636819839, + 1.9223335880121482867, + 1.9497480025204815182, + 1.9783232244565282241, + 2.008128468322387214, + 2.0392382187102560495, + 2.0717327197852989941, + 2.1056985191077205408, + 2.141229072951108682, + 2.1784254213865807416, + 2.2173969427163275031, + 2.2582621983975355419, + 2.3011498814497513326, + 2.3461998835490209547, + 2.393564498659336335, + 2.4434097842340786489, + 2.4959171048598266576, + 2.551284887865954588, + 2.6097306260799226329, + 2.6714931698159537872, + 2.7368353586570783274, + 2.8060470540283742602, + 2.8794486464792013645, + 2.9573951276637198642, + 3.0402808371050941667, + 3.1285450190960804891, + 3.2226783570437941673, + 3.3232306932306941103, + 3.4308201940528437269, + 3.5461442879894212769, + 3.6699927908828172788, + 3.8032637475035756812, + 3.9469826695032671633, + 4.1023260513255730736, + 4.2706503168149816862, + 4.4535277180208407742, + 4.6527912147146390367, + 4.8705910685478533488, + 5.1094668797693460142, + 5.3724402144314312228, + 5.6631350290659021695, + 5.9859361342534139213, + 6.3462004902284115273, + 6.7505430891683662786, + 7.2072300562145225172, + 7.7267290051348977324, + 8.3224952766020927734, + 9.0121210735938730413, + 9.8190591847972300101, + 10.775286943841163681, + 11.925568590252918077, + 13.334559234025861601, + 15.0992368270188772, + 17.371990147302994501, + 20.406794809689181251, + 24.660800535533617726, + 31.04829937949178548, + 41.702630469368096869, + 63.023799037287215586, + 127.01180974246945707, +] + +# The voigt HWHM with gamma=1, and sigma = (1+d)/(1-d), d = [-1, 1/64, 1) +voigt_hwhm_gammaone_dvar = [ + 1, + 1.0000929900981847014, + 1.0003777624966224696, + 1.0008631312648343249, + 1.0015580444997350155, + 1.0024715664851064116, + 1.0036128594929105533, + 1.0049911655464542273, + 1.0066157884679251467, + 1.0084960765229643227, + 1.0106414059536371252, + 1.0130611656602803461, + 1.0157647432548858631, + 1.0187615126671334742, + 1.0220608234423622704, + 1.0256719918320660707, + 1.0296042937449317882, + 1.033866959602362403, + 1.0384691711282855227, + 1.0434200600992189404, + 1.0487287090862781795, + 1.0544041542342061278, + 1.0604553901408076633, + 1.0668913769200059716, + 1.0737210495495320854, + 1.0809533296167089258, + 1.0885971395802894494, + 1.0966614196612918117, + 1.1051551474609148057, + 1.1140873603798353005, + 1.1234671808824950853, + 1.133303844615175872, + 1.143606731350947909, + 1.1543853987011884804, + 1.1656496185052021456, + 1.1774094157888474959, + 1.189675110171558587, + 1.2024573595995344503, + 1.2157672062912183539, + 1.2296161247989505185, + 1.2440160721168537852, + 1.2589795397983230491, + 1.2745196080855436114, + 1.2906500020968946739, + 1.3073851501646711064, + 1.3247402444642595085, + 1.3427313041259691636, + 1.3613752410716541342, + 1.3806899288698507996, + 1.4006942749554334374, + 1.4214082966130304507, + 1.4428532011781249068, + 1.4650514709665480456, + 1.4880269535027891822, + 1.5118049576811641308, + 1.5364123565625216811, + 1.5618776975840789442, + 1.5882313210425649123, + 1.6155054878026846865, + 1.6437345172857499843, + 1.6729549369091311623, + 1.7032056442782115564, + 1.7345280835813214269, + 1.7669664378076306961, + 1.8005678386015786904, + 1.8353825957890578829, + 1.8714644488648539776, + 1.9088708430231843335, + 1.9476632326499453211, + 1.9879074155840199742, + 2.0296739019046683743, + 2.0730383215232993655, + 2.1180818754635660479, + 2.1648918364190477434, + 2.2135621050012006813, + 2.2641938290540225637, + 2.3168960945429309357, + 2.3717866978564628894, + 2.4289930109305334564, + 2.4886529524650857967, + 2.5509160807121900781, + 2.6159448259472447961, + 2.6839158838829061716, + 2.7550217950632123133, + 2.8294727398261995659, + 2.9074985839268644614, + 2.989351216595417954, + 3.0753072309560758432, + 3.1656710067159081262, + 3.2607782673217607212, + 3.3610001989823174293, + 3.4667482378486090287, + 3.5784796552666793171, + 3.6967041007104309105, + 3.8219912995573272195, + 3.9549801506665303188, + 4.0963895299588992592, + 4.247031185217482386, + 4.4078252100521560122, + 4.579818719577299268, + 4.7642085282226360724, + 4.9623688672590590593, + 5.1758854988773873769, + 5.406598017846376464, + 5.6566527287853171409, + 5.9285693177078176118, + 6.2253257069142318691, + 6.5504671543400179081, + 6.9082480829105510743, + 7.3038186983166432837, + 7.7434738129206051218, + 8.2349894900724039365, + 8.7880859302556463323, + 9.4150755119073858072, + 10.131788564238479342, + 10.958926420044274577, + 11.924090995155746488, + 13.064921411472708807, + 14.43411259407456938, + 16.107779600011765899, + 18.200095193984788025, + 20.89047391778342413, + 24.477940229523170324, + 29.500738410375740968, + 37.035357331885579308, + 49.593604758544045262, + 74.710904558981651312, + 150.06437718602548694, +] + + +def almostequal(desired, actual, abstol, reltol): + delta = np.abs(abstol) + np.abs(reltol * desired) + return np.abs(desired - actual) <= delta + + +def test_empty_func(): + ff = voigtareaparametrizationnu.VoigtAreaParametrizationNu(0, -1) + p = lmfit.Parameters() + x = np.array([3]) + y = ff.func(p, x) + assert y[0] == 0, "empty function must return 0" + + +def test_constantbaseline(): + ff = voigtareaparametrizationnu.VoigtAreaParametrizationNu(0, 0) + p = lmfit.Parameters() + p.add("b0", 66) + x = np.array([3]) + y = ff.func(p, x) + assert y[0] == 66, "const background must return parameter b0" + + +def test_linearbaseline(): + ff = voigtareaparametrizationnu.VoigtAreaParametrizationNu(0, 1) + p = lmfit.Parameters() + p.add("b0", 66) + p.add("b1", 7) + x = np.array([3, 5]) + y = ff.func(p, x) + assert y[0] == 66 + 7 * 3, "linear background must return parameter b0 + b1 * x" + assert y[1] == 66 + 7 * 5, "linear background must return parameter b0 + b1 * x" + + +def test_oneterm_plus_linearbaseline(): + ff = VoigtAreaParametrizationNu(1, 1) + p = lmfit.Parameters() + p.add("b0", 4) + p.add("b1", -1) + p.add("area0", 17) + p.add("pos0", 1) + p.add("w0", 3) + p.add("nu0", 5.0 / 7.0) + x = np.array([3]) + y = ff.func(p, x) + assert ( + y[0] == 1.72519653836887579809275 + 4 - 1 * 3 + ), "test with one term and linear background" + + +def test_derivatives_generalcase(): + ff = VoigtAreaParametrizationNu(1, 1) + area = 17 + position = 7 + w = 3 + nu = 5 / 11.0 + b0 = 6 + b1 = -1 + expectedFunctionValue = 1.53969430505912493288384 + expectedDerivativeWrtArea = 0.0905702532387720548755199 + expectedDerivativeWrtPosition = 0.469766282659418172420053 + expectedDerivativeWrtW = -0.200053913246762862681244 + expectedDerivativeWrtNu = 0.677517478495227285639865 + p = lmfit.Parameters() + p.add("b0", b0) + p.add("b1", b1) + p.add("area0", area) + p.add("pos0", position) + p.add("w0", w) + p.add("nu0", nu) + x = np.array([9, 9]) + y = ff.func(p, x) + np.testing.assert_almost_equal(y[0], expectedFunctionValue + b0 + b1 * x[0], 12) + np.testing.assert_almost_equal(y[1], expectedFunctionValue + b0 + b1 * x[1], 12) + dy = ff.dfunc(p, x) + np.testing.assert_almost_equal(dy[0, 0], 1, decimal=12) + np.testing.assert_almost_equal(dy[1, 0], x[0], decimal=12) + np.testing.assert_almost_equal(dy[2, 0], expectedDerivativeWrtArea, decimal=12) + np.testing.assert_almost_equal(dy[3, 0], expectedDerivativeWrtPosition, decimal=12) + np.testing.assert_almost_equal(dy[4, 0], expectedDerivativeWrtW, decimal=12) + np.testing.assert_almost_equal(dy[5, 0], expectedDerivativeWrtNu, decimal=12) + + +def test_derivatives_lorentzlimit(): + ff = VoigtAreaParametrizationNu(1, -1) + area = 17 + position = 7 + w = 3 + nu = 1 / 16383.0 + expectedFunctionValue = 1.248792269568032230171214 + expectedDerivativeWrtArea = 0.0734561275808747703548694 + expectedDerivativeWrtPosition = 0.384232051961498798779317 + expectedDerivativeWrtW = -0.160096688317291166158049 + expectedDerivativeWrtNu = 0.624202576968383211902258 + p = lmfit.Parameters() + p.add("area0", area) + p.add("pos0", position) + p.add("w0", w) + p.add("nu0", nu) + x = np.array([9, 9]) + y = ff.func(p, x) + np.testing.assert_almost_equal(y[0], expectedFunctionValue, 12) + np.testing.assert_almost_equal(y[1], expectedFunctionValue, 12) + dy = ff.dfunc(p, x) + np.testing.assert_almost_equal(dy[0, 0], expectedDerivativeWrtArea, decimal=12) + np.testing.assert_almost_equal(dy[1, 0], expectedDerivativeWrtPosition, decimal=12) + np.testing.assert_almost_equal(dy[2, 0], expectedDerivativeWrtW, decimal=12) + np.testing.assert_almost_equal(dy[3, 0], expectedDerivativeWrtNu, decimal=12) + + +def test_derivatives_gausslimit(): + ff = VoigtAreaParametrizationNu(1, -1) + area = 17 + position = 7 + w = 3 + nu = 1 + expectedFunctionValue = 1.95602477676540327048287 + expectedDerivativeWrtArea = 0.1150602809862001923813453 + expectedDerivativeWrtPosition = 0.602583581831260309934069 + expectedDerivativeWrtW = -0.250285871034294216871578 + expectedDerivativeWrtNu = 0.865084621539303204435295 + p = lmfit.Parameters() + p.add("area0", area) + p.add("pos0", position) + p.add("w0", w) + p.add("nu0", nu) + x = np.array([9, 9]) + y = ff.func(p, x) + np.testing.assert_almost_equal(y[0], expectedFunctionValue, 12) + np.testing.assert_almost_equal(y[1], expectedFunctionValue, 12) + dy = ff.dfunc(p, x) + np.testing.assert_almost_equal(dy[0, 0], expectedDerivativeWrtArea, decimal=12) + np.testing.assert_almost_equal(dy[1, 0], expectedDerivativeWrtPosition, decimal=12) + np.testing.assert_almost_equal(dy[2, 0], expectedDerivativeWrtW, decimal=12) + np.testing.assert_almost_equal(dy[3, 0], expectedDerivativeWrtNu, decimal=12) + + +def test_fit_without_derivatives_oneterm(): + ff = VoigtAreaParametrizationNu(1, -1) + x = np.linspace(0, 99, 100) + w = 5 + nu = 6 / 11.0 + sigma = w * np.sqrt(nu) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma = w * (1 - nu) + data = 7 * voigt_profile(x - 50, sigma, gamma) + params = lmfit.Parameters() + params.add("area0", 10) + params.add("pos0", 49) + params.add("w0", 6) + params.add("nu0", 0.5) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq() + np.testing.assert_almost_equal(out1.params["area0"], 7, 12) + np.testing.assert_almost_equal(out1.params["pos0"], 50, 12) + np.testing.assert_almost_equal(out1.params["w0"], 5, 12) + np.testing.assert_almost_equal(out1.params["nu0"], 6 / 11.0, 12) + + +def test_fit_with_derivatives_oneterm(): + ff = voigtareaparametrizationnu.VoigtAreaParametrizationNu(1, -1) + x = np.linspace(0, 99, 100) + area = 7 + pos = 50 + w = 5 + nu = 6 / 11.0 + sigma = w * np.sqrt(nu) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma = w * (1 - nu) + data = area * voigt_profile(x - pos, sigma, gamma) + params = lmfit.Parameters() + + # important: the parameters must be added in exactly the order + # as you can see here (area0, pos0, w0, nu0, ... areaN, posN, wN, nuN, b0, b1, ... bM) + params.add("area0", 10) + params.add("pos0", 49) + params.add("w0", 6) + params.add("nu0", 0.5) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq(Dfun=ff.dfunc, col_deriv=1) + np.testing.assert_almost_equal(out1.params["area0"], area, 12) + np.testing.assert_almost_equal(out1.params["pos0"], pos, 12) + np.testing.assert_almost_equal(out1.params["w0"], w, 12) + np.testing.assert_almost_equal(out1.params["nu0"], nu, 12) + + +def test_fit_with_derivatives_twoterms_linearbaseline(): + ff = voigtareaparametrizationnu.VoigtAreaParametrizationNu(2, 1) + x = np.linspace(0, 99, 100) + area0 = 7 + pos0 = 25 + w0 = 3 + nu0 = 6 / 11.0 + area1 = 6 + pos1 = 75 + w1 = 4 + nu1 = 8 / 11.0 + b0 = 100 + b1 = 0.25 + sigma0 = w0 * np.sqrt(nu0) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma0 = w0 * (1 - nu0) + sigma1 = w1 * np.sqrt(nu1) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma1 = w1 * (1 - nu1) + data = ( + area0 * voigt_profile(x - pos0, sigma0, gamma0) + + area1 * voigt_profile(x - pos1, sigma1, gamma1) + + b0 + + b1 * x + ) + + # important: the parameters must be added in exactly the order + # as you can see here (b0, b1, ... bM, area0, pos0, w0, nu0, ... areaN, posN, wN, nuN, ) + params = lmfit.Parameters() + params.add("b0", 99) + params.add("b1", 0.24) + params.add("area0", 10) + params.add("pos0", 24) + params.add("w0", 3.5) + params.add("nu0", 0.5) + params.add("area1", 10) + params.add("pos1", 74) + params.add("w1", 4.5) + params.add("nu1", 0.5) + + min1 = lmfit.Minimizer(ff.func, params, fcn_args=(x,), fcn_kws={"data": data}) + out1 = min1.leastsq(Dfun=ff.dfunc, col_deriv=1) + cv = out1.covar + print(np.shape(cv)) + + np.testing.assert_almost_equal(out1.params["b0"], b0, 12) + np.testing.assert_almost_equal(out1.params["b1"], b1, 12) + np.testing.assert_almost_equal(out1.params["area0"], area0, 12) + np.testing.assert_almost_equal(out1.params["pos0"], pos0, 12) + np.testing.assert_almost_equal(out1.params["w0"], w0, 12) + np.testing.assert_almost_equal(out1.params["nu0"], nu0, 12) + np.testing.assert_almost_equal(out1.params["area1"], area1, 12) + np.testing.assert_almost_equal(out1.params["pos1"], pos1, 12) + np.testing.assert_almost_equal(out1.params["w1"], w1, 12) + np.testing.assert_almost_equal(out1.params["nu1"], nu1, 12) + + +def test_voigt_hwhm_exact_vargamma(): + sigma = 1 + for i in range(128): + d = -1 + i / 64.0 + gamma = (1 + d) / (1 - d) + hwhm = VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGamma( + sigma, gamma + ) + np.testing.assert_almost_equal(voigt_hwhm_sigmaone_dvar[i], hwhm, 0, 14) + i += 1 + + +def test_voigt_hwhm_Exact_varsigma(): + gamma = 1 + for i in range(128): + d = -1 + i / 64.0 + sigma = (1 + d) / (1 - d) + hwhm = VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGamma( + sigma, gamma + ) + np.testing.assert_almost_equal(voigt_hwhm_gammaone_dvar[i], hwhm, 0, 1e-14) + i += 1 + + +def test_voigt_hwhm_approx_vargamma(): + sigma = 1 + for i in range(128): + d = -1 + i / 64.0 + gamma = (1 + d) / (1 - d) + hwhm = VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGammaApproximation( + sigma, gamma + ) + assert ( + np.abs(1 - voigt_hwhm_sigmaone_dvar[i] / hwhm) + < voigtareaparametrizationnu.VoigtHalfWidthHalfMaximumApproximationMaximalRelativeError + ) + i += 1 + + +def test_voigt_hwhm_approx_varsigma(): + gamma = 1 + for i in range(128): + d = -1 + i / 64.0 + sigma = (1 + d) / (1 - d) + hwhm = VoigtAreaParametrizationNu.VoigtHalfWidthHalfMaximumOfSigmaGammaApproximation( + sigma, gamma + ) + assert ( + np.abs(1 - hwhm / voigt_hwhm_gammaone_dvar[i]) + < voigtareaparametrizationnu.VoigtHalfWidthHalfMaximumApproximationMaximalRelativeError + ) + i += 1 + + +def test_derivatives_SecondaryParameters_ExactlyGaussianLimit(): + # see VoigtArea-Derivatives-ParametrizationSqrtNuLog4.nb + area = 17 + pos = 7 + w = 3 + nu = 1 + ff = VoigtAreaParametrizationNu(1, -1) + + # position, area, height and FWHM + pahf = [pos, area, 2.66173895631567877902163, 6] + # derivatives + ymaxDerivs = [ + 0.156572879783275222295390, + 0, + -0.887246318771892926340544, + 1.16966732357221200231569, + ] + xmaxDerivs = [0, 1, 0, 0] + areaDerivs = [1, 0, 0, 0] + fwhmDerivs = [0, 0, 2, -0.207001611171248813604190] + + params = lmfit.Parameters() + params.add("area0", area) + params.add("pos0", pos) + params.add("w0", w) + params.add("nu0", nu) + + for i in range(4): + cov = np.zeros((4, 4)) + cov[i, i] = 1 + result = ff.GetPositionAreaHeightFWHMFromSinglePeakParameters(params, 0, cov) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-7) + + assert almostequal(pahf[0], result.Position, 1e-13, 1e-10) + assert almostequal(pahf[1], result.Area, 1e-13, 1e-10) + assert almostequal(pahf[2], result.Height, 1e-13, 1e-10) + assert almostequal(pahf[3], result.FWHM, 1e-13, 1e-10) + + +def test_derivatives_SecondaryParameters_NearlyGaussianLimit(): + # see VoigtArea-Derivatives-ParametrizationSqrtNuLog4.nb + area = 17 + pos = 7 + w = 3 + nu = 1 - 1 / 32767.0 + ff = VoigtAreaParametrizationNu(1, -1) + + # position, area, height and FWHM + pahf = [pos, area, 2.66170326013146224987024, 6.00000595967162241141443] + # derivatives + ymaxDerivs = [ + 0.156570780007733073521779, + 0, + -0.887234420043820749956747, + 1.16964641292668431368984, + ] + xmaxDerivs = [0, 1, 0, 0] + areaDerivs = [1, 0, 0, 0] + fwhmDerivs = [0, 0, 2.00000210576251930356528, -0.206995511602332887432860] + + params = lmfit.Parameters() + params.add("area0", area) + params.add("pos0", pos) + params.add("w0", w) + params.add("nu0", nu) + + for i in range(4): + cov = np.zeros((4, 4)) + cov[i, i] = 1 + result = ff.GetPositionAreaHeightFWHMFromSinglePeakParameters(params, 0, cov) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-7) + + assert almostequal(pahf[0], result.Position, 1e-13, 1e-10) + assert almostequal(pahf[1], result.Area, 1e-13, 1e-10) + assert almostequal(pahf[2], result.Height, 1e-13, 1e-10) + assert almostequal(pahf[3], result.FWHM, 1e-13, 1e-10) + + +def test_derivatives_SecondaryParameters_ExactlyLorentzianLimit(): + # see VoigtArea-Derivatives-ParametrizationSqrtNuLog4.nb + area = 17 + pos = 7 + w = 3 + nu = 0 + ff = VoigtAreaParametrizationNu(1, -1) + + # position, area, height and FWHM + pahf = [pos, area, 1.80375602170814713871402, 6] + # derivatives + ymaxDerivs = [ + 0.106103295394596890512589, + 0, + -0.601252007236049046238005, + 0.502621087962172486855736, + ] + xmaxDerivs = [0, 1, 0, 0] + areaDerivs = [1, 0, 0, 0] + fwhmDerivs = [0, 0, 2, 0.444686854097446089796046] + + params = lmfit.Parameters() + params.add("area0", area) + params.add("pos0", pos) + params.add("w0", w) + params.add("nu0", nu) + + for i in range(4): + cov = np.zeros((4, 4)) + cov[i, i] = 1 + result = ff.GetPositionAreaHeightFWHMFromSinglePeakParameters(params, 0, cov) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-7) + + assert almostequal(pahf[0], result.Position, 1e-13, 1e-10) + assert almostequal(pahf[1], result.Area, 1e-13, 1e-10) + assert almostequal(pahf[2], result.Height, 1e-13, 1e-10) + assert almostequal(pahf[3], result.FWHM, 1e-13, 1e-10) + + +def test_derivatives_SecondaryParameters_NearlyLorentzianLimit(): + # see VoigtArea-Derivatives-ParametrizationSqrtNuLog4.nb + area = 17 + pos = 7 + w = 3 + nu = 1 / 32767.0 + ff = VoigtAreaParametrizationNu(1, -1) + + # position, area, height and FWHM + pahf = [pos, area, 1.80377136162144979963231, 6.00001501741722156713951] + # derivatives + ymaxDerivs = [ + 0.106104197742438223507783, + 0, + -0.601257120540483266544103, + 0.502664788477749705785319, + ] + xmaxDerivs = [0, 1, 0, 0] + areaDerivs = [1, 0, 0, 0] + fwhmDerivs = [0, 0, 2.00000452341909985950189, 0.444626388979414132950005] + + params = lmfit.Parameters() + params.add("area0", area) + params.add("pos0", pos) + params.add("w0", w) + params.add("nu0", nu) + + for i in range(4): + cov = np.zeros((4, 4)) + cov[i, i] = 1 + result = ff.GetPositionAreaHeightFWHMFromSinglePeakParameters(params, 0, cov) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-7) + + assert almostequal(pahf[0], result.Position, 1e-13, 1e-10) + assert almostequal(pahf[1], result.Area, 1e-13, 1e-10) + assert almostequal(pahf[2], result.Height, 1e-13, 1e-10) + assert almostequal(pahf[3], result.FWHM, 1e-13, 1e-10) + + +def test_derivatives_SecondaryParameters_GeneralCase(): + # see VoigtArea-Derivatives-ParametrizationSqrtNuLog4.nb + area = 17 + pos = 7 + w = 3 + nu = 5 / 7.0 + ff = VoigtAreaParametrizationNu(1, -1) + + # position, area, height and FWHM + pahf = [pos, area, 2.35429674800053710973966, 6.04835222264692750740263] + # derivatives + ymaxDerivs = [ + 0.138488044000031594690569, + 0, + -0.784765582666845703246555, + 0.986294916351210818615044, + ] + xmaxDerivs = [0, 1, 0, 0] + areaDerivs = [1, 0, 0, 0] + fwhmDerivs = [0, 0, 2.01653911742859190892489, -0.134690859032882012668096] + + params = lmfit.Parameters() + params.add("area0", area) + params.add("pos0", pos) + params.add("w0", w) + params.add("nu0", nu) + + for i in range(4): + cov = np.zeros((4, 4)) + cov[i, i] = 1 + result = ff.GetPositionAreaHeightFWHMFromSinglePeakParameters(params, 0, cov) + assert almostequal(np.abs(areaDerivs[i]), result.AreaStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(xmaxDerivs[i]), result.PositionStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(ymaxDerivs[i]), result.HeightStdDev, 1e-13, 1e-7) + assert almostequal(np.abs(fwhmDerivs[i]), result.FWHMStdDev, 1e-13, 1e-7) + + assert almostequal(pahf[0], result.Position, 1e-13, 1e-10) + assert almostequal(pahf[1], result.Area, 1e-13, 1e-10) + assert almostequal(pahf[2], result.Height, 1e-13, 1e-10) + assert almostequal(pahf[3], result.FWHM, 1e-13, 1e-10) + + +def test_parameter_boundaries(): + # add directly, including max + a = lmfit.Parameters() + a.add("foo", value=0, max=1) + assert a["foo"].max == 1 + + # add parameter, and then modify + b = lmfit.Parameters() + b.add("foo", 1) + b["foo"].max = 66 + assert b["foo"].max == 66 + + # add parameter, and then modify + c = lmfit.Parameters() + d = lmfit.Parameter("foo", 34) + d.max = 77 + c.add(d) + assert c["foo"] == 34 + assert c["foo"].max == 77 + + +def test_fit_with_model_oneterm(): + x = np.linspace(0, 99, 100) + a = 7 + pos = 50 + w = 5 + nu = 6 / 11.0 + sigma = w * np.sqrt(nu) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma = w * (1 - nu) + data = a * voigt_profile(x - pos, sigma, gamma) + model = VoigtAreaParametrizationNuModel(x=x) + params = model.guess(data, x) + 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["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) diff --git a/tests/spectrum/peaks/test_PeakFittingByIncrementalPeakAddition.py b/tests/spectrum/peaks/test_PeakFittingByIncrementalPeakAddition.py new file mode 100644 index 00000000..891eab01 --- /dev/null +++ b/tests/spectrum/peaks/test_PeakFittingByIncrementalPeakAddition.py @@ -0,0 +1,93 @@ +from ramanchada2.fitting_functions import voigtareaparametrizationnu +from ramanchada2.fitting_functions.pearsonivamplitudeparametrizationhpw import ( + PearsonIVAmplitudeParametrizationHPW, +) +from ramanchada2.spectrum.peaks.find_peaks_ByIncrementalPeakAddition import ( + FindPeaksByIncrementalPeakAddition, +) +import numpy as np +from scipy.special import voigt_profile + + +def test_twotermsvoigt_linearbaseline(): + x = np.linspace(0, 99, 100) + area0 = 7 + pos0 = 25 + w0 = 3 + nu0 = 6 / 11.0 + area1 = 6 + pos1 = 75 + w1 = 4 + nu1 = 8 / 11.0 + b0 = 0.5 + b1 = 1 / 2048.0 + sigma0 = w0 * np.sqrt(nu0) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma0 = w0 * (1 - nu0) + sigma1 = w1 * np.sqrt(nu1) * voigtareaparametrizationnu.OneBySqrtLog4 + gamma1 = w1 * (1 - nu1) + y = ( + area0 * voigt_profile(x - pos0, sigma0, gamma0) + + area1 * voigt_profile(x - pos1, sigma1, gamma1) + + b0 + + b1 * x + ) + + fit = FindPeaksByIncrementalPeakAddition() + fit.maximumNumberOfPeaks = 2 + fit.orderOfBaselinePolynomial = 1 + + out = fit.Execute(x, y, voigtareaparametrizationnu.VoigtAreaParametrizationNu()) + params = out.params + np.testing.assert_almost_equal(params["b0"], b0, decimal=8) + np.testing.assert_almost_equal(params["b1"], b1, decimal=8) + np.testing.assert_almost_equal(params["area0"], area0, decimal=8) + np.testing.assert_almost_equal(params["pos0"], pos0, decimal=8) + np.testing.assert_almost_equal(params["w0"], w0, decimal=8) + np.testing.assert_almost_equal(params["nu0"], nu0, decimal=8) + np.testing.assert_almost_equal(params["area1"], area1, decimal=8) + np.testing.assert_almost_equal(params["pos1"], pos1, decimal=8) + np.testing.assert_almost_equal(params["w1"], w1, decimal=8) + np.testing.assert_almost_equal(params["nu1"], nu1, decimal=8) + + +def test_twotermspearsoniv_linearbaseline(): + x = np.linspace(0, 99, 100) + a0 = 7 + pos0 = 25 + w0 = 5 + m0 = 16 / 11.0 + v0 = 1 / 4.0 + a1 = 6 + pos1 = 75 + w1 = 6 + m1 = 21 / 11.0 + v1 = 1 / 8.0 + b0 = 100 + b1 = 0.25 + y = ( + PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(x, a0, pos0, w0, m0, v0) + + PearsonIVAmplitudeParametrizationHPW.GetYOfOneTerm(x, a1, pos1, w1, m1, v1) + + b0 + + b1 * x + ) + + fit = FindPeaksByIncrementalPeakAddition() + fit.maximumNumberOfPeaks = 2 + fit.orderOfBaselinePolynomial = 1 + + out = fit.Execute(x, y, PearsonIVAmplitudeParametrizationHPW()) + assert len(out.params) == 12 + # TODO: fit does not converge in the moment + # params = out.params + # np.testing.assert_almost_equal(params["b0"], b0, decimal=8) + # np.testing.assert_almost_equal(params["b1"], b1, decimal=8) + # np.testing.assert_almost_equal(params["a0"], a0, decimal=8) + # np.testing.assert_almost_equal(params["pos0"], pos0, decimal=8) + # np.testing.assert_almost_equal(params["w0"], w0, decimal=8) + # np.testing.assert_almost_equal(params["m0"], m0, decimal=8) + # np.testing.assert_almost_equal(params["v0"], v0, decimal=8) + # np.testing.assert_almost_equal(params["a1"], a1, decimal=8) + # np.testing.assert_almost_equal(params["pos1"], pos1, decimal=8) + # np.testing.assert_almost_equal(params["w1"], w1, decimal=8) + # np.testing.assert_almost_equal(params["m1"], m1, decimal=8) + # np.testing.assert_almost_equal(params["v1"], v1, decimal=8)