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