corems.ms_peak.calc.MSPeakCalc

   1__author__ = "Yuri E. Corilo"
   2__date__ = "Jun 04, 2019"
   3
   4import warnings
   5
   6import pyswarm
   7from lmfit import models
   8from numpy import (
   9    ceil,
  10    exp,
  11    flip,
  12    floor,
  13    linspace,
  14    log,
  15    nan,
  16    pi,
  17    poly1d,
  18    polyfit,
  19    rint,
  20    sqrt,
  21    square,
  22)
  23
  24try:
  25    from numpy import trapezoid
  26except ImportError:  # numpy < 2.0
  27    from numpy import trapz as trapezoid
  28
  29from corems.encapsulation.constant import Atoms
  30from corems.encapsulation.factory.parameters import MSParameters
  31
  32
  33class MSPeakCalculation:
  34    """Class to perform calculations on MSPeak objects.
  35
  36    This class provides methods to perform various calculations on MSPeak objects, such as calculating Kendrick Mass Defect (KMD) and Kendrick Mass (KM), calculating peak area, and fitting peak lineshape using different models.
  37
  38    Parameters
  39    ----------
  40    None
  41
  42    Attributes
  43    ----------
  44    _ms_parent : MSParent
  45        The parent MSParent object associated with the MSPeakCalculation object.
  46    mz_exp : float
  47        The experimental m/z value of the peak.
  48    peak_left_index : int
  49        The start scan index of the peak.
  50    peak_right_index : int
  51        The final scan index of the peak.
  52    resolving_power : float
  53        The resolving power of the peak.
  54
  55    Methods
  56    -------
  57    * _calc_kmd(dict_base).
  58        Calculate the Kendrick Mass Defect (KMD) and Kendrick Mass (KM) for a given base formula.
  59    * calc_area().
  60        Calculate the peak area using numpy's trapezoidal fit.
  61    * fit_peak(mz_extend=6, delta_rp=0, model='Gaussian').
  62        Perform lineshape analysis on a peak using lmfit module.
  63    * voigt_pso(w, r, yoff, width, loc, a).
  64        Calculate the Voigt function for particle swarm optimization (PSO) fitting.
  65    * objective_pso(x, w, u).
  66        Calculate the objective function for PSO fitting.
  67    * minimize_pso(lower, upper, w, u).
  68        Minimize the objective function using the particle swarm optimization algorithm.
  69    * fit_peak_pso(mz_extend=6, upsample_multiplier=5).
  70        Perform lineshape analysis on a peak using particle swarm optimization (PSO) fitting.
  71    * voigt(oversample_multiplier=1, delta_rp=0, mz_overlay=1).
  72        [Legacy] Perform voigt lineshape analysis on a peak.
  73    * pseudovoigt(oversample_multiplier=1, delta_rp=0, mz_overlay=1, fraction=0.5).
  74        [Legacy] Perform pseudovoigt lineshape analysis on a peak.
  75    * lorentz(oversample_multiplier=1, delta_rp=0, mz_overlay=1).
  76        [Legacy] Perform lorentz lineshape analysis on a peak.
  77    * gaussian(oversample_multiplier=1, delta_rp=0, mz_overlay=1).
  78        [Legacy] Perform gaussian lineshape analysis on a peak.
  79    * get_mz_domain(oversample_multiplier, mz_overlay).
  80        [Legacy] Resample/interpolate datapoints for lineshape analysis.
  81    * number_possible_assignments().
  82        Return the number of possible molecular formula assignments for the peak.
  83    * molecular_formula_lowest_error().
  84        Return the molecular formula with the smallest absolute mz error.
  85    * molecular_formula_highest_prob_score().
  86        Return the molecular formula with the highest confidence score.
  87    * molecular_formula_earth_filter(lowest_error=True).
  88        Filter molecular formula using the 'Earth' filter.
  89    * molecular_formula_water_filter(lowest_error=True).
  90        Filter molecular formula using the 'Water' filter.
  91    * molecular_formula_air_filter(lowest_error=True).
  92        Filter molecular formula using the 'Air' filter.
  93    * cia_score_S_P_error().
  94        Compound Identification Algorithm SP Error - Assignment Filter.
  95    * cia_score_N_S_P_error().
  96        Compound Identification Algorithm NSP Error - Assignment Filter.
  97
  98    """
  99
 100    def _calc_kmd(self, dict_base):
 101        """Calculate the Kendrick Mass Defect (KMD) and Kendrick Mass (KM) for a given base formula
 102
 103        Parameters
 104        ----------
 105        dict_base : dict
 106            dictionary with the base formula to be used in the calculation
 107            Default is CH2, e.g.
 108                dict_base = {"C": 1, "H": 2}
 109        """
 110
 111        if self._ms_parent:
 112            # msPeak obj does have a ms object parent
 113            kendrick_rounding_method = (
 114                self._ms_parent.mspeaks_settings.kendrick_rounding_method
 115            )  # rounding method can be one of floor, ceil or round
 116            # msPeak obj does not have a ms object parent
 117        else:
 118            kendrick_rounding_method = MSParameters.ms_peak.kendrick_rounding_method
 119
 120        mass = 0
 121        for atom in dict_base.keys():
 122            mass += Atoms.atomic_masses.get(atom) * dict_base.get(atom)
 123
 124        kendrick_mass = (int(mass) / mass) * self.mz_exp
 125
 126        if kendrick_rounding_method == "ceil":
 127            nominal_km = ceil(kendrick_mass)
 128
 129        elif kendrick_rounding_method == "round":
 130            nominal_km = rint(kendrick_mass)
 131
 132        elif kendrick_rounding_method == "floor":
 133            nominal_km = floor(kendrick_mass)
 134
 135        else:
 136            raise Exception(
 137                "%s method was not implemented, please refer to corems.ms_peak.calc.MSPeakCalc Class"
 138                % kendrick_rounding_method
 139            )
 140
 141        kmd = nominal_km - kendrick_mass
 142
 143        # kmd = (nominal_km - km) * 1
 144        # kmd = round(kmd,0)
 145
 146        return kmd, kendrick_mass, nominal_km
 147
 148    def calc_area(self):
 149        """Calculate the peak area using numpy's trapezoidal fit
 150
 151        uses provided mz_domain to accurately integrate areas independent of digital resolution
 152
 153        Returns
 154        -------
 155        float
 156            peak area
 157        """
 158        if self.peak_right_index > self.peak_left_index:
 159            yy = self._ms_parent.abundance_profile[
 160                self.peak_left_index : self.peak_right_index
 161            ]
 162            xx = self._ms_parent.mz_exp_profile[
 163                self.peak_left_index : self.peak_right_index
 164            ]
 165            # check if the axis is high to low m/z or not. if its MSFromFreq its high mz first, if its from Profile, its low mz first
 166            if xx[0] > xx[-1]:
 167                xx = flip(xx)
 168                yy = flip(yy)
 169            return float(trapezoid(yy, xx))
 170
 171        else:
 172            warnings.warn(
 173                "Peak Area Calculation for m/z {} has failed".format(self.mz_exp)
 174            )
 175            return nan
 176
 177    def fit_peak(self, mz_extend=6, delta_rp=0, model="Gaussian"):
 178        """Lineshape analysis on a peak using lmfit module.
 179
 180        Model and fit peak lineshape by defined function - using lmfit module
 181        Does not oversample/resample/interpolate data points
 182        Better to go back to time domain and perform more zero filling - if possible.
 183
 184        Parameters
 185        ----------
 186        mz_extend : int
 187            extra points left and right of peak definition to include in fitting
 188        delta_rp : float
 189            delta resolving power to add to resolving power
 190        model : str
 191            Type of lineshape model to use.
 192            Models allowed: Gaussian, Lorentz, Voigt
 193
 194        Returns
 195        -----
 196        mz_domain : ndarray
 197            x-axis domain for fit
 198        fit_peak : lmfit object
 199            fit results object from lmfit module
 200
 201        Notes
 202        -----
 203        Returns the calculated mz domain, initial defined abundance profile, and the fit peak results object from lmfit module
 204        mz_extend here extends the x-axis domain so that we have sufficient points either side of the apex to fit.
 205        Takes about 10ms per peak
 206        """
 207        start_index = (
 208            self.peak_left_index - mz_extend if not self.peak_left_index == 0 else 0
 209        )
 210        final_index = (
 211            self.peak_right_index + mz_extend
 212            if not self.peak_right_index == len(self._ms_parent.mz_exp_profile)
 213            else self.peak_right_index
 214        )
 215
 216        # check if MSPeak contains the resolving power info
 217        if self.resolving_power:
 218            # full width half maximum distance
 219            self.fwhm = self.mz_exp / (self.resolving_power + delta_rp)
 220
 221            mz_domain = self._ms_parent.mz_exp_profile[start_index:final_index]
 222            abundance_domain = self._ms_parent.abundance_profile[
 223                start_index:final_index
 224            ]
 225
 226            if model == "Gaussian":
 227                # stardard deviation
 228                sigma = self.fwhm / (2 * sqrt(2 * log(2)))
 229                amplitude = (sqrt(2 * pi) * sigma) * self.abundance
 230                model = models.GaussianModel()
 231                params = model.make_params(
 232                    center=self.mz_exp, amplitude=amplitude, sigma=sigma
 233                )
 234
 235            elif model == "Lorentz":
 236                # stardard deviation
 237                sigma = self.fwhm / 2
 238                amplitude = sigma * pi * self.abundance
 239                model = models.LorentzianModel()
 240                params = model.make_params(
 241                    center=self.mz_exp, amplitude=amplitude, sigma=sigma
 242                )
 243
 244            elif model == "Voigt":
 245                # stardard deviation
 246                sigma = self.fwhm / 3.6013
 247                amplitude = (sqrt(2 * pi) * sigma) * self.abundance
 248                model = models.VoigtModel()
 249                params = model.make_params(
 250                    center=self.mz_exp, amplitude=amplitude, sigma=sigma, gamma=sigma
 251                )
 252            else:
 253                raise LookupError("model lineshape not known or defined")
 254
 255            # calc_abundance = model.eval(params=params, x=mz_domain) #Same as initial fit, returned in fit_peak object
 256            fit_peak = model.fit(abundance_domain, params=params, x=mz_domain)
 257            return mz_domain, fit_peak
 258
 259        else:
 260            raise LookupError(
 261                "resolving power is not defined, try to use set_max_resolving_power()"
 262            )
 263
 264    def voigt_pso(self, w, r, yoff, width, loc, a):
 265        """Voigt function for particle swarm optimisation (PSO) fitting
 266
 267        From https://github.com/pnnl/nmrfit/blob/master/nmrfit/equations.py.
 268        Calculates a Voigt function over w based on the relevant properties of the distribution.
 269
 270        Parameters
 271        ----------
 272        w : ndarray
 273            Array over which the Voigt function will be evaluated.
 274        r : float
 275            Ratio between the Guassian and Lorentzian functions.
 276        yoff : float
 277            Y-offset of the Voigt function.
 278        width : float
 279            The width of the Voigt function.
 280        loc : float
 281            Center of the Voigt function.
 282        a : float
 283            Area of the Voigt function.
 284        Returns
 285        -------
 286        V : ndarray
 287            Array defining the Voigt function over w.
 288
 289        References
 290        ----------
 291        1. https://github.com/pnnl/nmrfit
 292
 293        Notes
 294        -----
 295        Particle swarm optimisation (PSO) fitting function can be significantly more computationally expensive than lmfit, with more parameters to optimise.
 296
 297        """
 298        # Lorentzian component
 299        L = (2 / (pi * width)) * 1 / (1 + ((w - loc) / (0.5 * width)) ** 2)
 300
 301        # Gaussian component
 302        G = (
 303            (2 / width)
 304            * sqrt(log(2) / pi)
 305            * exp(-(((w - loc) / (width / (2 * sqrt(log(2))))) ** 2))
 306        )
 307
 308        # Voigt body
 309        V = (yoff + a) * (r * L + (1 - r) * G)
 310
 311        return V
 312
 313    def objective_pso(self, x, w, u):
 314        """Objective function for particle swarm optimisation (PSO) fitting
 315
 316        The objective function used to fit supplied data.  Evaluates sum of squared differences between the fit and the data.
 317
 318        Parameters
 319        ----------
 320        x : list of floats
 321            Parameter vector.
 322        w : ndarray
 323            Array of frequency data.
 324        u : ndarray
 325            Array of data to be fit.
 326
 327        Returns
 328        -------
 329        rmse : float
 330            Root mean square error between the data and fit.
 331
 332        References
 333        ----------
 334        1. https://github.com/pnnl/nmrfit
 335
 336        """
 337        # global parameters
 338        r, width, loc, a = x
 339        yoff = 0
 340
 341        # calculate fit for V
 342        V_fit = self.voigt_pso(w, r, yoff, width, loc, a)
 343
 344        # real component RMSE
 345        rmse = sqrt(square((u - V_fit)).mean(axis=None))
 346
 347        # return the total RMSE
 348        return rmse
 349
 350    def minimize_pso(self, lower, upper, w, u):
 351        """Minimization function for particle swarm optimisation (PSO) fitting
 352
 353        Minimizes the objective function using the particle swarm optimization algorithm.
 354        Minimization function based on defined parameters
 355
 356
 357        Parameters
 358        ----------
 359        lower : list of floats
 360            Lower bounds for the parameters.
 361        upper : list of floats
 362            Upper bounds for the parameters.
 363        w : ndarray
 364            Array of frequency data.
 365        u : ndarray
 366            Array of data to be fit.
 367
 368        Notes
 369        -----
 370        Particle swarm optimisation (PSO) fitting function can be significantly more computationally expensive than lmfit, with more parameters to optimise.
 371        Current parameters take ~2 seconds per peak.
 372
 373
 374        References
 375        ----------
 376        1. https://github.com/pnnl/nmrfit
 377
 378        """
 379        # TODO - allow support to pass swarmsize, maxiter, omega, phip, phig parameters.
 380        # TODO - Refactor PSO fitting into its own class?
 381
 382        xopt, fopt = pyswarm.pso(
 383            self.objective_pso,
 384            lower,
 385            upper,
 386            args=(w, u),
 387            swarmsize=1000,
 388            maxiter=5000,
 389            omega=-0.2134,
 390            phip=-0.3344,
 391            phig=2.3259,
 392        )
 393        return xopt, fopt
 394
 395    def fit_peak_pso(self, mz_extend: int = 6, upsample_multiplier: int = 5):
 396        """Lineshape analysis on a peak using particle swarm optimisation (PSO) fitting
 397
 398        Function to fit a Voigt peakshape using particle swarm optimisation (PSO).
 399        Should return better results than lmfit, but much more computationally expensive
 400
 401        Parameters
 402        ----------
 403        mz_extend : int, optional
 404            extra points left and right of peak definition to include in fitting. Defaults to 6.
 405        upsample_multiplier : int, optional
 406            factor to increase x-axis points by for simulation of fitted lineshape function. Defaults to 5.
 407
 408        Returns
 409        -------
 410        xopt : array
 411            variables describing the voigt function.
 412            G/L ratio, width (fwhm), apex (x-axis), area.
 413            y-axis offset is fixed at 0
 414        fopt : float
 415            objective score (rmse)
 416        psfit : array
 417            recalculated y values based on function and optimised fit
 418        psfit_hdp : tuple of arrays
 419            0 - linspace x-axis upsampled grid
 420            1 - recalculated y values based on function and upsampled x-axis grid
 421            Does not change results, but aids in visualisation of the 'true' voigt lineshape
 422
 423        Notes
 424        -----
 425        Particle swarm optimisation (PSO) fitting function can be significantly more computationally expensive than lmfit, with more parameters to optimise.
 426        """
 427        # TODO - Add ability to pass pso args (i.e. swarm size, maxiter, omega, phig, etc)
 428        # TODO: fix xopt. Magnitude mode data through CoreMS/Bruker starts at 0 but is noise centered well above 0.
 429        # Thermo data is noise reduced by also noise subtracted, so starts at 0
 430        # Absorption mode/phased data will have positive and negative components and may not be baseline corrected
 431
 432        start_index = (
 433            self.peak_left_index - mz_extend if not self.peak_left_index == 0 else 0
 434        )
 435        final_index = (
 436            self.peak_right_index + mz_extend
 437            if not self.peak_right_index == len(self._ms_parent.mz_exp_profile)
 438            else self.peak_right_index
 439        )
 440
 441        # check if MSPeak contains the resolving power info
 442        if self.resolving_power:
 443            # full width half maximum distance
 444            self.fwhm = self.mz_exp / (self.resolving_power)
 445
 446            mz_domain = self._ms_parent.mz_exp_profile[start_index:final_index]
 447            abundance_domain = self._ms_parent.abundance_profile[
 448                start_index:final_index
 449            ]
 450            lower = [0, self.fwhm * 0.8, (self.mz_exp - 0.0005), 0]
 451            upper = [
 452                1,
 453                self.fwhm * 1.2,
 454                (self.mz_exp + 0.0005),
 455                self.abundance / self.signal_to_noise,
 456            ]
 457            xopt, fopt = self.minimize_pso(lower, upper, mz_domain, abundance_domain)
 458
 459            psfit = self.voigt_pso(mz_domain, xopt[0], 0, xopt[1], xopt[2], xopt[3])
 460            psfit_hdp_x = linspace(
 461                min(mz_domain), max(mz_domain), num=len(mz_domain) * upsample_multiplier
 462            )
 463            psfit_hdp = self.voigt_pso(
 464                psfit_hdp_x, xopt[0], 0, xopt[1], xopt[2], xopt[3]
 465            )
 466            return xopt, fopt, psfit, (psfit_hdp_x, psfit_hdp)
 467        else:
 468            raise LookupError(
 469                "resolving power is not defined, try to use set_max_resolving_power()"
 470            )
 471
 472    def voigt(self, oversample_multiplier=1, delta_rp=0, mz_overlay=1):
 473        """[Legacy] Voigt lineshape analysis function
 474        Legacy function for voigt lineshape analysis
 475
 476        Parameters
 477        ----------
 478        oversample_multiplier : int
 479            factor to increase x-axis points by for simulation of fitted lineshape function
 480        delta_rp : float
 481            delta resolving power to add to resolving power
 482        mz_overlay : int
 483            extra points left and right of peak definition to include in fitting
 484
 485        Returns
 486        -------
 487        mz_domain : ndarray
 488            x-axis domain for fit
 489        calc_abundance : ndarray
 490            calculated abundance profile based on voigt function
 491        """
 492
 493        if self.resolving_power:
 494            # full width half maximum distance
 495            self.fwhm = self.mz_exp / (
 496                self.resolving_power + delta_rp
 497            )  # self.resolving_power)
 498
 499            # stardart deviation
 500            sigma = self.fwhm / 3.6013
 501
 502            # half width baseline distance
 503
 504            # mz_domain = linspace(self.mz_exp - hw_base_distance,
 505            #                     self.mz_exp + hw_base_distance, datapoint)
 506            mz_domain = self.get_mz_domain(oversample_multiplier, mz_overlay)
 507
 508            # gaussian_pdf = lambda x0, x, s: (1/ math.sqrt(2*math.pi*math.pow(s,2))) * math.exp(-1 * math.pow(x-x0,2) / 2*math.pow(s,2) )
 509
 510            # TODO derive amplitude
 511            amplitude = (sqrt(2 * pi) * sigma) * self.abundance
 512
 513            model = models.VoigtModel()
 514
 515            params = model.make_params(
 516                center=self.mz_exp, amplitude=amplitude, sigma=sigma, gamma=sigma
 517            )
 518
 519            calc_abundance = model.eval(params=params, x=mz_domain)
 520
 521            return mz_domain, calc_abundance
 522
 523        else:
 524            raise LookupError(
 525                "resolving power is not defined, try to use set_max_resolving_power()"
 526            )
 527
 528    def pseudovoigt(
 529        self, oversample_multiplier=1, delta_rp=0, mz_overlay=1, fraction=0.5
 530    ):
 531        """[Legacy] pseudovoigt lineshape function
 532
 533        Legacy function for pseudovoigt lineshape analysis.
 534        Note - Code may not be functional currently.
 535
 536        Parameters
 537        ----------
 538        oversample_multiplier : int, optional
 539            factor to increase x-axis points by for simulation of fitted lineshape function. Defaults to 1.
 540        delta_rp : float, optional
 541            delta resolving power to add to resolving power. Defaults to 0.
 542        mz_overlay : int, optional
 543            extra points left and right of peak definition to include in fitting. Defaults to 1.
 544        fraction : float, optional
 545            fraction of gaussian component in pseudovoigt function. Defaults to 0.5.
 546
 547        """
 548        if self.resolving_power:
 549            # full width half maximum distance
 550            self.fwhm = self.mz_exp / (
 551                self.resolving_power + delta_rp
 552            )  # self.resolving_power)
 553
 554            # stardart deviation
 555            sigma = self.fwhm / 2
 556
 557            # half width baseline distance
 558
 559            # mz_domain = linspace(self.mz_exp - hw_base_distance,
 560            #                     self.mz_exp + hw_base_distance, datapoint)
 561            mz_domain = self.get_mz_domain(oversample_multiplier, mz_overlay)
 562
 563            # gaussian_pdf = lambda x0, x, s: (1/ math.sqrt(2*math.pi*math.pow(s,2))) * math.exp(-1 * math.pow(x-x0,2) / 2*math.pow(s,2) )
 564            model = models.PseudoVoigtModel()
 565
 566            # TODO derive amplitude
 567            gamma = sigma
 568
 569            amplitude = (sqrt(2 * pi) * sigma) * self.abundance
 570            amplitude = (sqrt(pi / log(2)) * (pi * sigma * self.abundance)) / (
 571                (pi * (1 - gamma)) + (sqrt(pi * log(2)) * gamma)
 572            )
 573
 574            params = model.make_params(center=self.mz_exp, sigma=sigma)
 575
 576            calc_abundance = model.eval(params=params, x=mz_domain)
 577
 578            return mz_domain, calc_abundance
 579
 580        else:
 581            raise LookupError(
 582                "resolving power is not defined, try to use set_max_resolving_power()"
 583            )
 584
 585    def lorentz(self, oversample_multiplier=1, delta_rp=0, mz_overlay=1):
 586        """[Legacy] Lorentz lineshape analysis function
 587
 588        Legacy function for lorentz lineshape analysis
 589
 590        Parameters
 591        ----------
 592        oversample_multiplier : int
 593            factor to increase x-axis points by for simulation of fitted lineshape function
 594        delta_rp : float
 595            delta resolving power to add to resolving power
 596        mz_overlay : int
 597            extra points left and right of peak definition to include in fitting
 598
 599        Returns
 600        -------
 601        mz_domain : ndarray
 602            x-axis domain for fit
 603        calc_abundance : ndarray
 604            calculated abundance profile based on lorentz function
 605
 606        """
 607        if self.resolving_power:
 608            # full width half maximum distance
 609            self.fwhm = self.mz_exp / (
 610                self.resolving_power + delta_rp
 611            )  # self.resolving_power)
 612
 613            # stardart deviation
 614            sigma = self.fwhm / 2
 615
 616            # half width baseline distance
 617            hw_base_distance = 8 * sigma
 618
 619            # mz_domain = linspace(self.mz_exp - hw_base_distance,
 620            #                     self.mz_exp + hw_base_distance, datapoint)
 621
 622            mz_domain = self.get_mz_domain(oversample_multiplier, mz_overlay)
 623            # gaussian_pdf = lambda x0, x, s: (1/ math.sqrt(2*math.pi*math.pow(s,2))) * math.exp(-1 * math.pow(x-x0,2) / 2*math.pow(s,2) )
 624            model = models.LorentzianModel()
 625
 626            amplitude = sigma * pi * self.abundance
 627
 628            params = model.make_params(
 629                center=self.mz_exp, amplitude=amplitude, sigma=sigma
 630            )
 631
 632            calc_abundance = model.eval(params=params, x=mz_domain)
 633
 634            return mz_domain, calc_abundance
 635
 636        else:
 637            raise LookupError(
 638                "resolving power is not defined, try to use set_max_resolving_power()"
 639            )
 640
 641    def gaussian(self, oversample_multiplier=1, delta_rp=0, mz_overlay=1):
 642        """[Legacy] Gaussian lineshape analysis function
 643        Legacy gaussian lineshape analysis function
 644
 645        Parameters
 646        ----------
 647        oversample_multiplier : int
 648            factor to increase x-axis points by for simulation of fitted lineshape function
 649        delta_rp : float
 650            delta resolving power to add to resolving power
 651        mz_overlay : int
 652            extra points left and right of peak definition to include in fitting
 653
 654        Returns
 655        -------
 656        mz_domain : ndarray
 657            x-axis domain for fit
 658        calc_abundance : ndarray
 659            calculated abundance profile based on gaussian function
 660
 661
 662        """
 663
 664        # check if MSPeak contains the resolving power info
 665        if self.resolving_power:
 666            # full width half maximum distance
 667            self.fwhm = self.mz_exp / (
 668                self.resolving_power + delta_rp
 669            )  # self.resolving_power)
 670
 671            # stardart deviation
 672            sigma = self.fwhm / (2 * sqrt(2 * log(2)))
 673
 674            # half width baseline distance
 675            # hw_base_distance = (3.2 * s)
 676
 677            # match_loz_factor = 3
 678
 679            # n_d = hw_base_distance * match_loz_factor
 680
 681            # mz_domain = linspace(
 682            #    self.mz_exp - n_d, self.mz_exp + n_d, datapoint)
 683
 684            mz_domain = self.get_mz_domain(oversample_multiplier, mz_overlay)
 685
 686            # gaussian_pdf = lambda x0, x, s: (1/ math.sqrt(2*math.pi*math.pow(s,2))) * math.exp(-1 * math.pow(x-x0,2) / 2*math.pow(s,2) )
 687
 688            # calc_abundance = norm.pdf(mz_domain, self.mz_exp, s)
 689
 690            model = models.GaussianModel()
 691
 692            amplitude = (sqrt(2 * pi) * sigma) * self.abundance
 693
 694            params = model.make_params(
 695                center=self.mz_exp, amplitude=amplitude, sigma=sigma
 696            )
 697
 698            calc_abundance = model.eval(params=params, x=mz_domain)
 699
 700            return mz_domain, calc_abundance
 701
 702        else:
 703            raise LookupError(
 704                "resolving power is not defined, try to use set_max_resolving_power()"
 705            )
 706
 707    def get_mz_domain(self, oversample_multiplier, mz_overlay):
 708        """[Legacy] function to resample/interpolate datapoints for lineshape analysis
 709
 710        This code is used for the legacy line fitting functions and not recommended.
 711        Legacy function to support expanding mz domain for legacy lineshape functions
 712
 713        Parameters
 714        ----------
 715        oversample_multiplier : int
 716            factor to increase x-axis points by for simulation of fitted lineshape function
 717        mz_overlay : int
 718            extra points left and right of peak definition to include in fitting
 719
 720        Returns
 721        -------
 722        mz_domain : ndarray
 723            x-axis domain for fit
 724
 725        """
 726        start_index = (
 727            self.peak_left_index - mz_overlay if not self.peak_left_index == 0 else 0
 728        )
 729        final_index = (
 730            self.peak_right_index + mz_overlay
 731            if not self.peak_right_index == len(self._ms_parent.mz_exp_profile)
 732            else self.peak_right_index
 733        )
 734
 735        if oversample_multiplier == 1:
 736            mz_domain = self._ms_parent.mz_exp_profile[start_index:final_index]
 737
 738        else:
 739            # we assume a linear correlation for m/z and datapoits
 740            # which is only true if the m/z range in narrow (within 1 m/z unit)
 741            # this is not true for a wide m/z range
 742
 743            indexes = range(start_index, final_index + 1)
 744            mz = self._ms_parent.mz_exp_profile[indexes]
 745            pol = poly1d(polyfit(indexes, mz, 1))
 746            oversampled_indexes = linspace(
 747                start_index,
 748                final_index,
 749                (final_index - start_index) * oversample_multiplier,
 750            )
 751            mz_domain = pol(oversampled_indexes)
 752
 753        return mz_domain
 754
 755    @property
 756    def number_possible_assignments(
 757        self,
 758    ):
 759        return len(self.molecular_formulas)
 760
 761    def molecular_formula_lowest_error(self):
 762        """Return the molecular formula with the smallest absolute mz error"""
 763
 764        return min(self.molecular_formulas, key=lambda m: abs(m.mz_error))
 765
 766    def molecular_formula_highest_prob_score(self):
 767        """Return the molecular formula with the highest confidence score score"""
 768
 769        return max(self.molecular_formulas, key=lambda m: abs(m.confidence_score))
 770
 771    def molecular_formula_earth_filter(self, lowest_error=True):
 772        """Filter molecular formula using the 'Earth' filter
 773
 774        This function applies the Formularity-esque 'Earth' filter to possible molecular formula assignments.
 775        Earth Filter:
 776            O > 0 AND N <= 3 AND P <= 2 AND 3P <= O
 777
 778        If the lowest_error method is also used, it will return the single formula annotation with the smallest absolute error which also fits the Earth filter.
 779        Otherwise, it will return all Earth-filter compliant formulas.
 780
 781        Parameters
 782        ----------
 783        lowest_error : bool, optional.
 784            Return only the lowest error formula which also fits the Earth filter.
 785            If False, return all Earth-filter compliant formulas. Default is True.
 786
 787        Returns
 788        -------
 789        list
 790            List of molecular formula objects which fit the Earth filter
 791
 792        References
 793        ----------
 794        1. Nikola Tolic et al., "Formularity: Software for Automated Formula Assignment of Natural and Other Organic Matter from Ultrahigh-Resolution Mass Spectra"
 795            Anal. Chem. 2017, 89, 23, 12659–12665
 796            doi: 10.1021/acs.analchem.7b03318
 797        """
 798
 799        candidates = list(
 800            filter(
 801                lambda mf: (
 802                    mf.get("O") > 0
 803                    and mf.get("N") <= 3
 804                    and mf.get("P") <= 2
 805                    and (3 * mf.get("P")) <= mf.get("O")
 806                ),
 807                self.molecular_formulas,
 808            )
 809        )
 810        if len(candidates) > 0:
 811            if lowest_error:
 812                return min(candidates, key=lambda m: abs(m.mz_error))
 813            else:
 814                return candidates
 815        else:
 816            return candidates
 817
 818    def molecular_formula_water_filter(self, lowest_error=True):
 819        """Filter molecular formula using the 'Water' filter
 820
 821        This function applies the Formularity-esque 'Water' filter to possible molecular formula assignments.
 822        Water Filter:
 823            O > 0 AND N <= 3 AND S <= 2 AND P <= 2
 824
 825        If the lowest_error method is also used, it will return the single formula annotation with the smallest absolute error which also fits the Water filter.
 826        Otherwise, it will return all Water-filter compliant formulas.
 827
 828        Parameters
 829        ----------
 830        lowest_error : bool, optional
 831            Return only the lowest error formula which also fits the Water filter.
 832            If False, return all Water-filter compliant formulas. Defaults to 2
 833
 834        Returns
 835        -------
 836        list
 837            List of molecular formula objects which fit the Water filter
 838
 839        References
 840        ----------
 841        1. Nikola Tolic et al., "Formularity: Software for Automated Formula Assignment of Natural and Other Organic Matter from Ultrahigh-Resolution Mass Spectra"
 842            Anal. Chem. 2017, 89, 23, 12659–12665
 843            doi: 10.1021/acs.analchem.7b03318
 844        """
 845
 846        candidates = list(
 847            filter(
 848                lambda mf: (
 849                    mf.get("O") > 0
 850                    and mf.get("N") <= 3
 851                    and mf.get("S") <= 2
 852                    and mf.get("P") <= 2
 853                ),
 854                self.molecular_formulas,
 855            )
 856        )
 857        if len(candidates) > 0:
 858            if lowest_error:
 859                return min(candidates, key=lambda m: abs(m.mz_error))
 860            else:
 861                return candidates
 862        else:
 863            return candidates
 864
 865    def molecular_formula_air_filter(self, lowest_error=True):
 866        """Filter molecular formula using the 'Air' filter
 867
 868        This function applies the Formularity-esque 'Air' filter to possible molecular formula assignments.
 869        Air Filter:
 870            O > 0 AND N <= 3 AND S <= 1 AND P = 0 AND 3(S+N) <= O
 871
 872        If the lowest_error method is also used, it will return the single formula annotation with the smallest absolute error which also fits the Air filter.
 873        Otherwise, it will return all Air-filter compliant formulas.
 874
 875        Parameters
 876        ----------
 877        lowest_error : bool, optional
 878            Return only the lowest error formula which also fits the Air filter.
 879            If False, return all Air-filter compliant formulas. Defaults to True.
 880
 881        Returns
 882        -------
 883        list
 884            List of molecular formula objects which fit the Air filter
 885
 886        References
 887        ----------
 888        1. Nikola Tolic et al., "Formularity: Software for Automated Formula Assignment of Natural and Other Organic Matter from Ultrahigh-Resolution Mass Spectra"
 889            Anal. Chem. 2017, 89, 23, 12659–12665
 890            doi: 10.1021/acs.analchem.7b03318
 891        """
 892
 893        candidates = list(
 894            filter(
 895                lambda mf: (
 896                    mf.get("O") > 0
 897                    and mf.get("N") <= 2
 898                    and mf.get("S") <= 1
 899                    and mf.get("P") == 0
 900                    and 3 * (mf.get("S") + mf.get("N")) <= mf.get("O")
 901                ),
 902                self.molecular_formulas,
 903            )
 904        )
 905
 906        if len(candidates) > 0:
 907            if lowest_error:
 908                return min(candidates, key=lambda m: abs(m.mz_error))
 909            else:
 910                return candidates
 911        else:
 912            return candidates
 913
 914    def cia_score_S_P_error(self):
 915        """Compound Identification Algorithm SP Error - Assignment Filter
 916
 917        This function applies the Compound Identification Algorithm (CIA) SP Error filter to possible molecular formula assignments.
 918
 919        It takes the molecular formula with the lowest S+P count, and returns the formula with the lowest absolute error from this subset.
 920
 921        Returns
 922        -------
 923        MolecularFormula
 924            A single molecular formula which fits the rules of the CIA SP Error filter
 925
 926
 927        References
 928        ----------
 929        1. Elizabeth B. Kujawinski and Mark D. Behn, "Automated Analysis of Electrospray Ionization Fourier Transform Ion Cyclotron Resonance Mass Spectra of Natural Organic Matter"
 930            Anal. Chem. 2006, 78, 13, 4363–4373
 931            doi: 10.1021/ac0600306
 932        """
 933        # case EFormulaScore.HAcap:
 934
 935        lowest_S_P_mf = min(
 936            self.molecular_formulas, key=lambda mf: mf.get("S") + mf.get("P")
 937        )
 938        lowest_S_P_count = lowest_S_P_mf.get("S") + lowest_S_P_mf.get("P")
 939
 940        list_same_s_p = list(
 941            filter(
 942                lambda mf: mf.get("S") + mf.get("P") == lowest_S_P_count,
 943                self.molecular_formulas,
 944            )
 945        )
 946
 947        # check if list is not empty
 948        if list_same_s_p:
 949            return min(list_same_s_p, key=lambda m: abs(m.mz_error))
 950
 951        else:
 952            return lowest_S_P_mf
 953
 954    def cia_score_N_S_P_error(self):
 955        """Compound Identification Algorithm NSP Error - Assignment Filter
 956
 957        This function applies the Compound Identification Algorithm (CIA) NSP Error filter to possible molecular formula assignments.
 958
 959        It takes the molecular formula with the lowest N+S+P count, and returns the formula with the lowest absolute error from this subset.
 960
 961        Returns
 962        -------
 963        MolecularFormula
 964            A single molecular formula which fits the rules of the CIA NSP Error filter
 965
 966        References
 967        ----------
 968        1. Elizabeth B. Kujawinski and Mark D. Behn, "Automated Analysis of Electrospray Ionization Fourier Transform Ion Cyclotron Resonance Mass Spectra of Natural Organic Matter"
 969            Anal. Chem. 2006, 78, 13, 4363–4373
 970            doi: 10.1021/ac0600306
 971
 972        Raises
 973        -------
 974        Exception
 975            If no molecular formula are associated with mass spectrum peak.
 976        """
 977        # case EFormulaScore.HAcap:
 978        if self.molecular_formulas:
 979            lowest_N_S_P_mf = min(
 980                self.molecular_formulas,
 981                key=lambda mf: mf.get("N") + mf.get("S") + mf.get("P"),
 982            )
 983            lowest_N_S_P_count = (
 984                lowest_N_S_P_mf.get("N")
 985                + lowest_N_S_P_mf.get("S")
 986                + lowest_N_S_P_mf.get("P")
 987            )
 988
 989            list_same_N_S_P = list(
 990                filter(
 991                    lambda mf: (
 992                        mf.get("N") + mf.get("S") + mf.get("P") == lowest_N_S_P_count
 993                    ),
 994                    self.molecular_formulas,
 995                )
 996            )
 997
 998            if list_same_N_S_P:
 999                SP_filtered_list = list(
1000                    filter(
1001                        lambda mf: (mf.get("S") <= 3) and (mf.get("P") <= 1),
1002                        list_same_N_S_P,
1003                    )
1004                )
1005
1006                if SP_filtered_list:
1007                    return min(SP_filtered_list, key=lambda m: abs(m.mz_error))
1008
1009                else:
1010                    return min(list_same_N_S_P, key=lambda m: abs(m.mz_error))
1011
1012            else:
1013                return lowest_N_S_P_mf
1014        else:
1015            raise Exception(
1016                "No molecular formula associated with the mass spectrum peak at m/z: %.6f"
1017                % self.mz_exp
1018            )