Modelling

class AutoPeakDetection:
    """Automatic detection of peaks in a spectrum."""

    def __init__(
        self,
        x: NDArray[np.float64],
        data: NDArray[np.float64],
        args: Dict[str, Any],
    ) -> None:
        """Initialize the AutoPeakDetection class.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            data (NDArray[np.float64]): `y`-values of the data as 1d-array.
            args (Dict[str, Any]): The input file arguments as a dictionary with
                 additional information beyond the command line arguments.
        """
        self.x = x
        self.data = data
        self._args = args["autopeak"]

    @staticmethod
    def check_key_exists(
        key: str, args: Dict[str, Any], value: Union[float, Tuple[Any, Any]]
    ) -> Any:
        """Check if a key exists in a dictionary.

        Please check for the reference key also [scipy.signal.find_peaks][1].

        [1]:
            https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html


        Args:
            key (str): Reference key of `scipy.signal.find_peaks`.
            args (Dict[str, Any]): Reference values of `scipy.signal.find_peaks`, if not
                 defined will be set to estimated default values.
            value (Union[float, Tuple[float,float]]): Default value for the reference
                 key.

        Returns:
            Any: The reference value for `scipy.signal.find_peaks`.
        """
        return args.get(key, value)

    @property
    def estimate_height(self) -> Tuple[float, float]:
        r"""Estimate the initial height based on an inverse noise ratio of a signal.

        !!! info "About the estimation of the height"

            The lower end of the height is the inverse noise ratio of the `data`, and
            upper limit is the maximum value of the `data`. The noise ratio of the
            `data` is based on the original implementation by `SciPy`:

            ```python
            def signaltonoise(a, axis=0, ddof=0):
                a = np.asanyarray(a)
                m = a.mean(axis)
                sd = a.std(axis=axis, ddof=ddof)
                return np.where(sd == 0, 0, m / sd)
            ```

        Returns:
            Tuple[float, float]: Tuple of the inverse signal to noise ratio and
                 the maximum value of the `data`.
        """
        return 1 - self.data.mean() / self.data.std(), self.data.max()

    @property
    def estimate_threshold(self) -> Tuple[float, float]:
        """Estimate the threshold value for the peak detection.

        Returns:
            Tuple[float, float]: Minimum and maximum value of the spectrum `data`,
                 respectively, `intensity`.
        """
        return self.data.min(), self.data.max()

    @property
    def estimate_distance(self) -> float:
        """Estimate the initial distance between peaks.

        Returns:
            float: Estimated distance between peaks.
        """
        min_step = np.diff(self.x).min()
        return max(min_step, 1.0)

    @property
    def estimate_prominence(self) -> Tuple[float, float]:
        """Estimate the prominence of a peak.

        !!! info "About the estimation of the prominence"

            The prominence is the difference between the height of the peak and the
            bottom. To get a estimate of the prominence, the height of the peak is
            calculated by maximum value of the `data` and the bottom is calculated by
            the harmonic mean of the `data`.

        Returns:
            Tuple[float, float]: Tuple of the harmonic-mean and maximum value of `data`.
        """
        try:
            return hmean(self.data), self.data.max()
        except ValueError as exc:
            print(f"{exc}: Using standard arithmetic mean of NumPy.\n")
        return self.data.mean(), self.data.max()

    @property
    def estimated_width(self) -> Tuple[float, float]:
        """Estimate the width of a peak.

        !!! info "About the estimation of the width"

            The width of a peak is estimated for a lower and an upper end. For the lower
            end, the minimum stepsize is used. For the upper end, the stepsize between
            the half maximum and the minimum value of the `data` is used as the width.

        Returns:
            Tuple[float, float]: Estimated width lower and uper end of the peaks.
        """
        return (
            np.diff(self.x).min(),
            np.abs(self.x[self.data.argmax()] - self.x[self.data.argmin()]) / 2,
        )

    @property
    def estimated_rel_height(self) -> float:
        """Estimate the relative height of a peak.

        !!! info "About the estimation of the relative height"

            The relative height of a peak is approximated by the difference of the
            harmonic mean value of the `data` and the minimum value of the `data`
            divided by the factor of `4`. In case of negative ratios, the value will be
            set to `Zero`.

        Returns:
            float: Estimated relative height of a peak.
        """
        try:
            rel_height = (hmean(self.data) - self.data.min()) / 4
        except ValueError as exc:
            print(f"{exc}: Using standard arithmetic mean of NumPy.\n")
            rel_height = (self.data.mean() - self.data.min()) / 4
        return rel_height if rel_height > 0 else 0.0

    @property
    def estimated_wlen(self) -> float:
        r"""Estimate the window length for the peak detection.

        !!! info "About the estimation of the window length"

            The window length is the length of the window for the peak detection is
            defined to be 1% of the length of the `data`, consequently the len of the
            `data` is divided by 100. In case of a window length smaller than 1, the
            window length will be set to numerical value of 1, which is defined by
            `1 + 1e-9`.

        Returns:
            float: Estimated window length is set to the numeric value of > 1.
        """
        wlen = self.data.size / 100
        return wlen if wlen > 1.0 else 1 + 1e-9

    @property
    def estimated_plateau_size(self) -> Tuple[float, float]:
        """Estimate the plateau size for the peak detection.

        Returns:
            Tuple[float, float]: Estimated plateau size is set to `zero` for the lower
                 end and the maximum value of the `x` for the upper end.
        """
        return 0.0, self.x.max()

    def initialize_peak_detection(self) -> None:
        """Initialize the peak detection.

        !!! note "Initialize the peak detection"

            This method is used to initialize the peak detection. The initialization can
            be activated by setting the `initialize` attribute to `True`, which will
            automatically estimate the default parameters for the peak detection. In
            case of the `initialize` attribute is defined as dictionary, the proposed
            values are taken from the dictionary if th

        Raise:
            TypeError: If the `initialize` attribute is not of type `bool` or `dict`.
        """
        if isinstance(self._args, bool):
            self.default_values()
        elif isinstance(self._args, dict):
            ReferenceKeys().detection_check(self._args)
            self.height = self.check_key_exists(
                key="height", args=self._args, value=self.estimate_height
            )
            self.threshold = self.check_key_exists(
                key="threshold", args=self._args, value=self.estimate_threshold
            )
            self.distance = self.check_key_exists(
                key="distance", args=self._args, value=self.estimate_distance
            )
            self.prominence = self.check_key_exists(
                key="prominence", args=self._args, value=self.estimate_prominence
            )
            self.width = self.check_key_exists(
                key="width", args=self._args, value=self.estimated_width
            )
            self.wlen = self.check_key_exists(
                key="wlen", args=self._args, value=self.estimated_wlen
            )
            self.rel_height = self.check_key_exists(
                key="rel_height", args=self._args, value=self.estimated_rel_height
            )
            self.plateau_size = self.check_key_exists(
                key="plateau_size", args=self._args, value=0.0
            )
        else:
            raise TypeError(
                f"The type of the `args` is not supported: {type(self._args)}"
            )

    def default_values(self) -> None:
        """Set the default values for the peak detection."""
        self.height = self.estimate_height
        self.threshold = self.estimate_threshold
        self.distance = self.estimate_distance
        self.prominence = self.estimate_prominence
        self.width = self.estimated_width
        self.wlen = self.estimated_wlen
        self.rel_height = self.estimated_rel_height
        self.plateau_size = 0

    def __autodetect__(self) -> Any:
        """Return peak positions and properties."""
        return find_peaks(
            self.data,
            height=self.height,
            threshold=self.threshold,
            distance=self.distance,
            prominence=self.prominence,
            width=self.width,
            wlen=self.wlen,
            rel_height=self.rel_height,
            plateau_size=self.plateau_size,
        )

@dataclass(frozen=True)
class Constants:
    r"""Mathematical constants for the curve models.

    !!! info "Constants"

        1. Natural logarithm of 2

            $$
            ln2 = \log{2}
            $$

        2. Square root of 2 times pi

            $$
            sq2pi = \sqrt{2 \pi}
            $$

        3. Square root of pi

            $$
            sqpi = \sqrt{ \pi}
            $$

        4. Square root of 2

            $$
            sq2 = \sqrt{2}
            $$

        5. Full width at half maximum to sigma for Gaussian

            $$
            fwhmg2sig = \frac{1}{ 2 \sqrt{2\log{2}}}
            $$

        6. Full width at half maximum to sigma for Lorentzian

            $$
            fwhml2sig = \frac{1}{2}
            $$

        7. Full width at half maximum to sigma for Voigt according to the article by
            Olivero and Longbothum[^1], check also
            [XPSLibary website](https://xpslibrary.com/voigt-peak-shape/).

            $$
            fwhm_{\text{Voigt}} \approx 0.5346 \cdot fwhm_{\text{Gaussian}} +
              \sqrt{ 0.2166 fwhm_{\text{Lorentzian}}^2  + fwhm_{\text{Gaussian}}^2 }

            $$

            In case of equal FWHM for Gaussian and Lorentzian, the Voigt FWHM can be
            defined as:

            $$
            fwhm_{\text{Voigt}} \approx 1.0692 + 2 \sqrt{0.2166 + 2 \ln{2}} \cdot \sigma
            $$

            $$
            fwhmv2sig = \frac{1}{fwhm_{\text{Voigt}}}
            $$

        [^1]:
            J.J. Olivero, R.L. Longbothum,
            _Empirical fits to the Voigt line width: A brief review_,
            **Journal of Quantitative Spectroscopy and Radiative Transfer**,
            Volume 17, Issue 2, 1977, Pages 233-236, ISSN 0022-4073,
            https://doi.org/10.1016/0022-4073(77)90161-3.
    """

    ln2 = log(2.0)
    sq2pi = sqrt(2.0 * pi)
    sqpi = sqrt(pi)
    sq2 = sqrt(2.0)
    fwhmg2sig = 1 / (2.0 * sqrt(2.0 * log(2.0)))
    fwhml2sig = 1 / 2.0
    fwhmv2sig = 1 / (2 * 0.5346 + 2 * sqrt(0.2166 + log(2) * 2))

class DistributionModels:
    """Distribution models for the fit.

    !!! note "About distribution models"

        `DistributionModels` are wrapper functions for the distribution models. The
        overall goal is to extract from the best parameters the single contributions in
        the model. The superposition of the single contributions is the final model.

    !!! note "About the cumulative distribution"

        The cumulative distribution is the sum of the single contributions. The
        cumulative distribution is the model that is fitted to the data. In contrast to
        the single contributions, the cumulative distribution is not normalized and
        therefore the amplitude of the single contributions is not directly comparable
        to the amplitude of the cumulative distribution. Also, the cumulative
        distributions are consquently using the `fwhm` parameter instead of the
        `sigma` parameter.
    """

    @staticmethod
    def gaussian(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhmg: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional Gaussian distribution.

        $$
        {\displaystyle g(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}\exp
        (  -{\frac {1}{2}}{\frac {(x-\mu )^{2}}{\sigma ^{2}}} ) }
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Gaussian distribution.
                 Defaults to 1.0.
            center (float, optional): Center of the Gaussian distribution.
                 Defaults to 0.0.
            fwhmg (float, optional): Full width at half maximum (FWHM) of the Gaussian
                distribution. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Gaussian distribution of `x` given.
        """
        sigma = fwhmg * Constants.fwhmg2sig
        return np.array(amplitude / (Constants.sq2pi * sigma)) * np.exp(
            -((1.0 * x - center) ** 2) / (2 * sigma**2)
        )

    @staticmethod
    def lorentzian(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhml: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional Lorentzian distribution.

        $$
        f(x;x_{0},\gamma )={\frac  {1}{\pi \gamma
        [ 1+ ( {\frac  {x-x_{0}}{\gamma }})^{2} ]
        }} ={1 \over \pi \gamma } [ {\gamma ^{2} \over (x-x_{0})^{2}+\gamma ^{2}} ]
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Lorentzian distribution.
                Defaults to 1.0.
            center (float, optional): Center of the Lorentzian distribution. Defaults to
                0.0.
            fwhml (float, optional): Full width at half maximum (FWHM) of the Lorentzian
                distribution. Defaults to 1.0.

        Returns:
            Union[NDArray[np.float64], float]: Lorentzian distribution of `x` given.
        """
        sigma = fwhml * Constants.fwhml2sig
        return np.array(amplitude / (1 + ((1.0 * x - center) / sigma) ** 2)) / (
            pi * sigma
        )

    @staticmethod
    def voigt(
        x: NDArray[np.float64],
        center: float = 0.0,
        fwhmv: float = 1.0,
        gamma: Optional[float] = None,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional Voigt distribution.

        $$
        {\displaystyle V(x;\sigma ,\gamma )\equiv
        \int_{-\infty }^{\infty }G(x';\sigma )
        L(x-x';\gamma )\,dx'}
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            center (float, optional): Center of the Voigt distribution. Defaults to 0.0.
            fwhmv (float, optional): Full width at half maximum (FWHM) of the Lorentzian
                distribution. Defaults to 1.0.
            gamma (float, optional): Scaling factor of the complex part of the
                [Faddeeva Function](https://en.wikipedia.org/wiki/Faddeeva_function).
                Defaults to None.

        Returns:
            NDArray[np.float64]: Voigt distribution of `x` given.
        """
        sigma = fwhmv * Constants.fwhmv2sig
        if gamma is None:
            gamma = sigma
        z = (x - center + 1j * gamma) / (sigma * Constants.sq2)
        return np.array(wofz(z).real / (sigma * Constants.sq2pi))

    @staticmethod
    def pseudovoigt(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhmg: float = 1.0,
        fwhml: float = 1.0,
    ) -> NDArray[np.float64]:
        """Return a 1-dimensional Pseudo-Voigt distribution.

        !!! note "See also:"

            J. Appl. Cryst. (2000). 33, 1311-1316
            https://doi.org/10.1107/S0021889800010219

        Args:
            x (NDArray[np.float64]):  `x`-values of the data.
            amplitude (float, optional): Amplitude of the Pseudo-Voigt distribution.
                Defaults to 1.0.
            center (float, optional): Center of the Pseudo-Voigt distribution.
                Defaults to 0.0.
            fwhmg (float, optional): Full width half maximum of the Gaussian
                distribution in the Pseudo-Voigt distribution. Defaults to 1.0.
            fwhml (float, optional): Full width half maximum of the Lorentzian
                distribution in the Pseudo-Voigt distribution. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Pseudo-Voigt distribution of `x` given.
        """
        f = np.power(
            fwhmg**5
            + 2.69269 * fwhmg**4 * fwhml
            + 2.42843 * fwhmg**3 * fwhml**2
            + 4.47163 * fwhmg**2 * fwhml**3
            + 0.07842 * fwhmg * fwhml**4
            + fwhml**5,
            0.2,
        )
        n = (
            1.36603 * (fwhml / f)
            - 0.47719 * (fwhml / f) ** 2
            + 0.11116 * (fwhml / f) ** 3
        )
        return np.array(
            n
            * DistributionModels.lorentzian(
                x=x, amplitude=amplitude, center=center, fwhml=fwhml
            )
            + (1 - n)
            * DistributionModels.gaussian(
                x=x, amplitude=amplitude, center=center, fwhmg=fwhmg
            )
        )

    @staticmethod
    def exponential(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        decay: float = 1.0,
        intercept: float = 0.0,
    ) -> NDArray[np.float64]:
        """Return a 1-dimensional exponential decay.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the exponential function.
                 Defaults to 1.0.
            decay (float, optional): Decay of the exponential function. Defaults to 1.0.
            intercept (float, optional): Intercept of the exponential function.
                 Defaults to 0.0.

        Returns:
            NDArray[np.float64]: Exponential decay of `x` given.
        """
        return np.array(amplitude * np.exp(-x / decay) + intercept)

    @staticmethod
    def power(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        exponent: float = 1.0,
        intercept: float = 0.0,
    ) -> NDArray[np.float64]:
        """Return a 1-dimensional power function.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the power function. Defaults to
                1.0.
            exponent (float, optional): Exponent of the power function. Defaults to 1.0.
            intercept (float, optional): Intercept of the power function. Defaults to
                0.0.

        Returns:
            NDArray[np.float64]: power function of `x` given.
        """
        return np.array(amplitude * np.power(x, exponent) + intercept)

    @staticmethod
    def linear(
        x: NDArray[np.float64],
        slope: float = 1.0,
        intercept: float = 0.0,
    ) -> NDArray[np.float64]:
        """Return a 1-dimensional linear function.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            slope (float, optional): Slope of the linear function. Defaults to 1.0.
            intercept (float, optional): Intercept of the linear function.
                 Defaults to 0.0.

        Returns:
            NDArray[np.float64]: Linear function of `x` given.
        """
        return np.array(slope * x + intercept)

    @staticmethod
    def constant(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
    ) -> NDArray[np.float64]:
        """Return a 1-dimensional constant value.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the constant. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Constant value of `x` given.
        """
        return np.array(np.linspace(amplitude, amplitude, len(x)))

    @staticmethod
    def _norm(
        x: NDArray[np.float64], center: float, sigma: float
    ) -> NDArray[np.float64]:
        """Normalize the data for step functions.

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            center (float): Center of the step function.
            sigma (float): Sigma of the step function.

        Returns:
            NDArray[np.float64]: Normalized data.
        """
        if abs(sigma) < 1.0e-13:
            sigma = 1.0e-13
        return np.subtract(x, center) / sigma

    @staticmethod
    def erf(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        sigma: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional error function.

        $$
        f(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the error function.
                    Defaults to 1.0.
            center (float, optional): Center of the error function. Defaults to 0.0.
            sigma (float, optional): Sigma of the error function. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Error function of `x` given.
        """
        return np.array(
            amplitude * 0.5 * (1 + erf(DistributionModels._norm(x, center, sigma)))
        )

    @staticmethod
    def heaviside(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        sigma: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional Heaviside step function.

        $$
        f(x) = \begin{cases}
        0 & x < 0 \\
        0.5 & x = 0 \\
        1 & x > 0
        \end{cases}
        $$
        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Heaviside step function.
                    Defaults to 1.0.
            center (float, optional): Center of the Heaviside step function.
                 Defaults to 0.0.
            sigma (float, optional): Sigma of the Heaviside step function.
                 Defaults to 1.0.


        Returns:
            NDArray[np.float64]: Heaviside step function of `x` given.
        """
        return np.array(
            amplitude * 0.5 * (1 + np.sign(DistributionModels._norm(x, center, sigma)))
        )

    @staticmethod
    def atan(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        sigma: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional arctan step function.

        $$
        f(x) = \frac{1}{\pi} \arctan(\frac{x - c}{s})
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the arctan step function.
                    Defaults to 1.0.
            center (float, optional): Center of the arctan step function.
                 Defaults to 0.0.
            sigma (float, optional): Sigma of the arctan step function.
                 Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Arctan step function of `x` given.
        """
        return np.array(
            amplitude
            * 0.5
            * (1 + np.arctan(DistributionModels._norm(x, center, sigma)) / pi)
        )

    @staticmethod
    def log(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        sigma: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional logarithmic step function.

        $$
        f(x) = \frac{1}{1 + e^{-\frac{x - c}{s}}}
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the logarithmic step function.
                    Defaults to 1.0.
            center (float, optional): Center of the logarithmic step function.
                 Defaults to 0.0.
            sigma (float, optional): Sigma of the logarithmic step function.
                 Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Logarithmic step function of `x` given.
        """
        return np.array(
            amplitude
            * 0.5
            * (1 + np.log(DistributionModels._norm(x, center, sigma)) / pi)
        )

    @staticmethod
    def cgaussian(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhmg: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional cumulative Gaussian function.

        $$
        f(x) = \frac{1}{2} \left[1 + erf\left(\frac{x - c}{s \sqrt{2}}\right)\right]
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Gaussian function. Defaults to
                1.0.
            center (float, optional): Center of the Gaussian function. Defaults to 0.0.
            fwhmg (float, optional): Full width at half maximum of the Gaussian
                 function. Defaults to 1.0.


        Returns:
            NDArray[np.float64]: Cumulative Gaussian function of `x` given.
        """
        sigma = fwhmg * Constants.fwhmg2sig
        return np.array(
            amplitude * 0.5 * (1 + erf((x - center) / (sigma * np.sqrt(2.0))))
        )

    @staticmethod
    def clorentzian(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhml: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional cumulative Lorentzian function.

        $$
        f(x) = \frac{1}{\pi} \arctan\left(\frac{x - c}{s}\right) + \frac{1}{2}
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Lorentzian function.
                    Defaults to 1.0.
            center (float, optional): Center of the Lorentzian function.
                 Defaults to 0.0.
            fwhml (float, optional): Full width at half maximum of the Lorentzian
                function. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Cumulative Lorentzian function of `x` given.
        """
        sigma = fwhml * Constants.fwhml2sig
        return np.array(amplitude * (np.arctan((x - center) / sigma) / pi) + 0.5)

    @staticmethod
    def cvoigt(
        x: NDArray[np.float64],
        amplitude: float = 1.0,
        center: float = 0.0,
        fwhmv: float = 1.0,
        gamma: float = 1.0,
    ) -> NDArray[np.float64]:
        r"""Return a 1-dimensional cumulative Voigt function.

        $$
        f(x) = \frac{1}{2} \left[1 + erf\left(\frac{x - c}{s \sqrt{2}}\right)\right]
        $$

        Args:
            x (NDArray[np.float64]): `x`-values of the data.
            amplitude (float, optional): Amplitude of the Voigt function. Defaults to
                1.0.
            center (float, optional): Center of the Voigt function. Defaults to 0.0.
            fwhmv (float, optional): Full width at half maximum of the Voigt function.
                Defaults to 1.0.
            gamma (float, optional): Gamma of the Voigt function. Defaults to 1.0.

        Returns:
            NDArray[np.float64]: Cumulative Voigt function of `x` given.
        """
        sigma = fwhmv * Constants.fwhmv2sig
        return np.array(
            amplitude
            * 0.5
            * (1 + erf((x - center) / (sigma * np.sqrt(2.0))))
            * np.exp(-(((x - center) / gamma) ** 2))
        )