Implement helper functions for getting the logprobs of stuff

Author: michalsta

src/IsoSpec++/cwrapper.cpp

@@ -59,16 +59,31 @@ double getLightestPeakMassIso(void* iso)
     return reinterpret_cast<Iso*>(iso)->getLightestPeakMass();
 }
 
+double getLightestPeakLProbIso(void* iso)
+{
+    return reinterpret_cast<Iso*>(iso)->getLightestPeakLProb();
+}
+
 double getHeaviestPeakMassIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getHeaviestPeakMass();
 }
 
+double getHeaviestPeakLProbIso(void* iso)
+{
+    return reinterpret_cast<Iso*>(iso)->getHeaviestPeakLProb();
+}
+
 double getMonoisotopicPeakMassIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getMonoisotopicPeakMass();
 }
 
+double getMonoisotopicPeakLProbIso(void* iso)
+{
+    return reinterpret_cast<Iso*>(iso)->getMonoisotopicPeakLProb();
+}
+
 double getModeLProbIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getModeLProb();

src/IsoSpec++/cwrapper.h

@@ -40,8 +40,11 @@ ISOSPEC_C_API void * setupIso(int             dimNumber,
 ISOSPEC_C_API void * isoFromFasta(const char* fasta, bool use_nominal_masses, bool add_water);
 
 ISOSPEC_C_API double getLightestPeakMassIso(void* iso);
+ISOSPEC_C_API double getLightestPeakLProbIso(void* iso);
 ISOSPEC_C_API double getHeaviestPeakMassIso(void* iso);
+ISOSPEC_C_API double getHeaviestPeakLProbIso(void* iso);
 ISOSPEC_C_API double getMonoisotopicPeakMassIso(void* iso);
+ISOSPEC_C_API double getMonoisotopicPeakLProbIso(void* iso);
 ISOSPEC_C_API double getModeLProbIso(void* iso);
 ISOSPEC_C_API double getModeMassIso(void* iso);
 ISOSPEC_C_API double getTheoreticalAverageMassIso(void* iso);

src/IsoSpec++/isoSpec++.cpp

@@ -250,6 +250,14 @@ double Iso::getLightestPeakMass() const
     return mass;
 }
 
+double Iso::getLightestPeakLProb() const
+{
+    double lprob = 0.0;
+    for (int ii = 0; ii < dimNumber; ii++)
+        lprob += marginals[ii]->getLightestConfLProb();
+    return lprob;
+}
+
 double Iso::getHeaviestPeakMass() const
 {
     double mass = 0.0;
@@ -258,6 +266,14 @@ double Iso::getHeaviestPeakMass() const
     return mass;
 }
 
+double Iso::getHeaviestPeakLProb() const
+{
+    double lprob = 0.0;
+    for (int ii = 0; ii < dimNumber; ii++)
+        lprob += marginals[ii]->getHeaviestConfLProb();
+    return lprob;
+}
+
 double Iso::getMonoisotopicPeakMass() const
 {
     double mass = 0.0;
@@ -266,6 +282,14 @@ double Iso::getMonoisotopicPeakMass() const
     return mass;
 }
 
+double Iso::getMonoisotopicPeakLProb() const
+{
+    double lprob = 0.0;
+    for (int ii = 0; ii < dimNumber; ii++)
+        lprob += marginals[ii]->getMonoisotopicConfLProb();
+    return lprob;
+}
+
 double Iso::getUnlikeliestPeakLProb() const
 {
     double lprob = 0.0;

src/IsoSpec++/isoSpec++.h

@@ -135,16 +135,29 @@ class ISOSPEC_EXPORT_SYMBOL Iso {
     //! Get the mass of the lightest peak in the isotopic distribution.
     double getLightestPeakMass() const;
 
+    //! Get the log-probability of the lightest peak in the isotopic distribution.
+    double getLightestPeakLProb() const;
+
     //! Get the mass of the heaviest peak in the isotopic distribution.
     double getHeaviestPeakMass() const;
 
+    //! Get the log-probability of the heaviest peak in the isotopic distribution.
+    double getHeaviestPeakLProb() const;
+
     /*!
         Get the mass of the monoisotopic peak in the isotopic distribution. Monoisotopc molecule is defined as
         consisting only of the most frequent isotopes of each element. These are often (but not always) the
         lightest ones, making this often (but again, not always) equal to getLightestPeakMass()
     */
     double getMonoisotopicPeakMass() const;
 
+    /*!
+        Get the log-probability of the monoisotopic peak in the isotopic distribution. Monoisotopc molecule is defined as
+        consisting only of the most frequent isotopes of each element. These are often (but not always) the
+        lightest ones. Making this often (but again, not always) equal to getLightestPeakLProb()
+    */
+    double getMonoisotopicPeakLProb() const;
+
     //! Get the log-probability of the mode-configuration (if there are many modes, they share this value).
     double getModeLProb() const;
 

src/IsoSpec++/marginalTrek++.cpp

@@ -256,6 +256,22 @@ double Marginal::getLightestConfMass() const
     return ret_mass*atomCnt;
 }
 
+double Marginal::getLightestConfLProb() const
+{
+    std::unique_ptr<int[]> lightest_conf = std::make_unique<int[]>(isotopeNo);
+    size_t lightest_idx = 0;
+    double min_mass = std::numeric_limits<double>::infinity();
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        if( min_mass > atom_masses[ii] )
+        {
+            min_mass = atom_masses[ii];
+            lightest_idx = ii;
+        }
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        lightest_conf[ii] = (ii == lightest_idx ? atomCnt : 0);
+    return logProb(lightest_conf.get());
+}
+
 double Marginal::getHeaviestConfMass() const
 {
     double ret_mass = 0.0;
@@ -265,6 +281,22 @@ double Marginal::getHeaviestConfMass() const
     return ret_mass*atomCnt;
 }
 
+double Marginal::getHeaviestConfLProb() const
+{
+    std::unique_ptr<int[]> heaviest_conf = std::make_unique<int[]>(isotopeNo);
+    size_t heaviest_idx = 0;
+    double max_mass = -std::numeric_limits<double>::infinity();
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        if( max_mass < atom_masses[ii] )
+        {
+            max_mass = atom_masses[ii];
+            heaviest_idx = ii;
+        }
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        heaviest_conf[ii] = (ii == heaviest_idx ? atomCnt : 0);
+    return logProb(heaviest_conf.get());
+}
+
 size_t Marginal::getMonoisotopicAtomIndex() const
 {
     double found_prob = -std::numeric_limits<double>::infinity();
@@ -284,6 +316,15 @@ double Marginal::getMonoisotopicConfMass() const
     return atom_masses[idx]*atomCnt;
 }
 
+double Marginal::getMonoisotopicConfLProb() const
+{
+    size_t idx = getMonoisotopicAtomIndex();
+    std::unique_ptr<int[]> monoisotopic_conf = std::make_unique<int[]>(isotopeNo);
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        monoisotopic_conf[ii] = (ii == idx ? atomCnt : 0);
+    return logProb(monoisotopic_conf.get());
+}
+
 double Marginal::getAtomAverageMass() const
 {
     double ret = 0.0;

src/IsoSpec++/marginalTrek++.h

@@ -96,19 +96,28 @@ class Marginal
     */
     double getLightestConfMass() const;
 
+    //! Get the log-probability of the lightest subisotopologue.
+    double getLightestConfLProb() const;
+
     //! Get the mass of the heaviest subisotopologue.
     /*! This is trivially obtained by considering all atomNo atoms to be the heaviest isotope possible.
         \return The mass of the heaviest subisotopologue.
     */
     double getHeaviestConfMass() const;
 
+    //! Get the log-probability of the heaviest subisotopologue.
+    double getHeaviestConfLProb() const;
+
     //! Get the mass of the monoisotopic subisotopologue.
     /*! The monoisotopic subisotopologue is defined as the molecule consiting only
         of the most likely isotope. This is frequently the lightest subisotopologue,
         making this frequently (but not always) equal to getLightestConfMass()
     */
     double getMonoisotopicConfMass() const;
 
+    //! Get the log-probability of the monoisotopic subisotopologue.
+    double getMonoisotopicConfLProb() const;
+
     //! Get the index of the monoisotopic (most probable) atom.
     /*!
         \return The index of the most probable isotope of the element.

src/IsoSpecPy/IsoSpecPy.py

@@ -219,14 +219,26 @@ def getLightestPeakMass(self):
         """Get the lightest peak in the isotopic distribution."""
         return self.ffi.getLightestPeakMassIso(self.iso)
 
+    def getLightestPeakLProb(self):
+        """Get the log probability of the lightest peak in the isotopic distribution."""
+        return self.ffi.getLightestPeakLProbIso(self.iso)
+
     def getHeaviestPeakMass(self):
         """Get the heaviest peak in the isotopic distribution."""
         return self.ffi.getHeaviestPeakMassIso(self.iso)
 
+    def getHeaviestPeakLProb(self):
+        """Get the log probability of the heaviest peak in the isotopic distribution."""
+        return self.ffi.getHeaviestPeakLProbIso(self.iso)
+
     def getMonoisotopicPeakMass(self):
         """Get the monoisotopic mass of the peak."""
         return self.ffi.getMonoisotopicPeakMassIso(self.iso)
 
+    def getMonoisotopicPeakLProb(self):
+        """Get the log probability of the monoisotopic peak in the isotopic distribution."""
+        return self.ffi.getMonoisotopicPeakLProbIso(self.iso)
+
     def getModeLProb(self):
         """Get the log probability of the most probable peak(s) in the isotopic distribution."""
         return self.ffi.getModeLProbIso(self.iso)

src/IsoSpecPy/isoFFI.py

@@ -3,13 +3,16 @@
 import sys
 import glob
 from pathlib import Path
+
 if False:
     import IsoSpecCppPy
 
+
 class IsoFFI:
     def __init__(self):
         self.ffi = cffi.FFI()
-        self.ffi.cdef('''
+        self.ffi.cdef(
+            """
         void * setupIso(int dimNumber,
                 const int* isotopeNumbers,
                 const int* atomCounts,
@@ -19,8 +22,11 @@ def __init__(self):
         void * isoFromFasta(const char* fasta, bool use_nominal_masses, bool add_water);
 
         double getLightestPeakMassIso(void* iso);
+        double getLightestPeakLProbIso(void* iso);
         double getHeaviestPeakMassIso(void* iso);
+        double getHeaviestPeakLProbIso(void* iso);
         double getMonoisotopicPeakMassIso(void* iso);
+        double getMonoisotopicPeakLProbIso(void* iso);
         double getModeLProbIso(void* iso);
         double getModeMassIso(void* iso);
         double getTheoreticalAverageMassIso(void* iso);
@@ -157,39 +163,47 @@ def __init__(self):
         extern const char* elem_table_element[NUMBER_OF_ISOTOPIC_ENTRIES];
         extern const char* elem_table_symbol[NUMBER_OF_ISOTOPIC_ENTRIES];
         extern const bool elem_table_Radioactive[NUMBER_OF_ISOTOPIC_ENTRIES];
-                        ''');
+                        """
+        )
 
         mod_dir = Path(__file__).resolve().parent
 
-        if (mod_dir.parent / 'setup.py').exists():
-            raise ImportError('''Attempted to load IsoSpecPy module from its build directory. This usually
+        if (mod_dir.parent / "setup.py").exists():
+            raise ImportError(
+                """Attempted to load IsoSpecPy module from its build directory. This usually
 won't work and is generally a Bad Idea. Please cd somewhere else, or, if you're really
-sure you want to do that, edit the source and disable this check.''')
+sure you want to do that, edit the source and disable this check."""
+            )
 
-        libnames  = ['IsoSpecCppPy*', 'IsoSpec++*']
-        libprefix = ['', 'lib', 'Lib']
-        extension = ['.so', '.dylib', '.dll']
+        libnames = ["IsoSpecCppPy*", "IsoSpec++*"]
+        libprefix = ["", "lib", "Lib"]
+        extension = [".so", ".dylib", ".dll"]
         try:
-            if platform.system() == 'Linux':
-                extension = ['.so', 'pyd']
-            elif platform.system() == 'Windows':
-                extension = ['.dll', '.pyd']
+            if platform.system() == "Linux":
+                extension = [".so", "pyd"]
+            elif platform.system() == "Windows":
+                extension = [".dll", ".pyd"]
         except:
             pass
 
-        prebuilt =  ['', 'prebuilt-']
+        prebuilt = ["", "prebuilt-"]
 
         def cprod(ll1, ll2):
             res = []
             for l1 in ll1:
                 for l2 in ll2:
-                    res.append(l1+l2)
+                    res.append(l1 + l2)
             return res
 
         paths_to_check = cprod(prebuilt, cprod(libprefix, cprod(libnames, extension)))
         dpc = []
 
-        for dirpath in [mod_dir, mod_dir.parent, mod_dir.parent / 'bin', mod_dir.parent / 'lib']:
+        for dirpath in [
+            mod_dir,
+            mod_dir.parent,
+            mod_dir.parent / "bin",
+            mod_dir.parent / "lib",
+        ]:
             dpc.extend([dirpath / p for p in paths_to_check])
 
         paths_to_check = dpc
@@ -200,7 +214,10 @@ def cprod(ll1, ll2):
         paths_to_check = expanded
         try:
             import importlib
-            paths_to_check.insert(0, Path(importlib.util.find_spec("IsoSpecCppPy").origin))
+
+            paths_to_check.insert(
+                0, Path(importlib.util.find_spec("IsoSpecCppPy").origin)
+            )
         except (ImportError, AttributeError):
             pass
 
@@ -213,11 +230,21 @@ def cprod(ll1, ll2):
                 self.libpath = libpath
                 break
             except (IndexError, OSError) as e:
-                errmsg = "Load libIsoSpec++.so, tried: " + libpath + '\n' + "Got error: " + str(type(e)) + ": " + str(e)
+                errmsg = (
+                    "Load libIsoSpec++.so, tried: "
+                    + libpath
+                    + "\n"
+                    + "Got error: "
+                    + str(type(e))
+                    + ": "
+                    + str(e)
+                )
                 errors.append(errmsg)
 
         if self.clib == None:
-            raise ImportError("Cannot find or load the C++ part of the library\n" + '\n'.join(errors))
+            raise ImportError(
+                "Cannot find or load the C++ part of the library\n" + "\n".join(errors)
+            )
 
 
 isoFFI = IsoFFI()  # This is done while including the module

tests/Python/test_iface.py

@@ -167,6 +167,21 @@ def test_get_monoisotopic_mass():
     assert math.isclose(68885.198515667, mono_mass, rel_tol=1e-9)
     print("OK!")
 
+def test_lightest_peak():
+    print("Checking lightest peak...", end=" ")
+    formula = "C10B10H10Sn1"
+    iso = IsoSpecPy.Iso(formula=formula)
+    lightest_mass = iso.getLightestPeakMass()
+    lightest_lprob = iso.getLightestPeakLProb()
+    iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula)
+    masses = list(iso_threshold.masses)
+    probs = list(iso_threshold.probs)
+    min_index = masses.index(min(masses))
+    print(lightest_mass, lightest_lprob, end=" ")
+    assert lightest_mass == masses[min_index]
+    assert math.isclose(lightest_lprob, math.log(probs[min_index]), rel_tol=1e-9)
+    print("OK!")
+
 if __name__ == "__main__":
     test_wasserstein_distance()
     test_normalization()
@@ -181,3 +196,4 @@ def test_get_monoisotopic_mass():
     test_empiric_variance()
     test_empiric_stddev()
     test_get_monoisotopic_mass()
+    test_lightest_peak()