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