Implement retrieval of special confs

Author: michalsta

src/IsoSpec++/cwrapper.cpp

@@ -64,6 +64,11 @@ double getLightestPeakLProbIso(void* iso)
     return reinterpret_cast<Iso*>(iso)->getLightestPeakLProb();
 }
 
+void getLightestPeakSignature(void* iso, int* space)
+{
+    reinterpret_cast<Iso*>(iso)->getLightestPeakSignature(space);
+}
+
 double getHeaviestPeakMassIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getHeaviestPeakMass();
@@ -74,6 +79,11 @@ double getHeaviestPeakLProbIso(void* iso)
     return reinterpret_cast<Iso*>(iso)->getHeaviestPeakLProb();
 }
 
+void getHeaviestPeakSignature(void* iso, int* space)
+{
+    reinterpret_cast<Iso*>(iso)->getHeaviestPeakSignature(space);
+}
+
 double getMonoisotopicPeakMassIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getMonoisotopicPeakMass();
@@ -84,6 +94,11 @@ double getMonoisotopicPeakLProbIso(void* iso)
     return reinterpret_cast<Iso*>(iso)->getMonoisotopicPeakLProb();
 }
 
+void getMonoisotopicPeakSignature(void* iso, int* space)
+{
+    reinterpret_cast<Iso*>(iso)->getMonoisotopicPeakSignature(space);
+}
+
 double getModeLProbIso(void* iso)
 {
     return reinterpret_cast<Iso*>(iso)->getModeLProb();

src/IsoSpec++/cwrapper.h

@@ -41,10 +41,13 @@ ISOSPEC_C_API void * isoFromFasta(const char* fasta, bool use_nominal_masses, bo
 
 ISOSPEC_C_API double getLightestPeakMassIso(void* iso);
 ISOSPEC_C_API double getLightestPeakLProbIso(void* iso);
+ISOSPEC_C_API void getLightestPeakSignature(void* iso, int* space);
 ISOSPEC_C_API double getHeaviestPeakMassIso(void* iso);
 ISOSPEC_C_API double getHeaviestPeakLProbIso(void* iso);
+ISOSPEC_C_API void getHeaviestPeakSignature(void* iso, int* space);
 ISOSPEC_C_API double getMonoisotopicPeakMassIso(void* iso);
 ISOSPEC_C_API double getMonoisotopicPeakLProbIso(void* iso);
+ISOSPEC_C_API void getMonoisotopicPeakSignature(void* iso, int* space);
 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

@@ -258,6 +258,17 @@ double Iso::getLightestPeakLProb() const
     return lprob;
 }
 
+void Iso::getLightestPeakSignature(int* space) const
+{
+    for (int ii = 0; ii < dimNumber; ii++)
+    {
+        size_t lightest_idx = marginals[ii]->getLightestAtomIndex();
+        for(int jj = 0; jj < isotopeNumbers[ii]; jj++)
+            space[jj] = (jj == lightest_idx ? atomCounts[ii] : 0);
+        space += isotopeNumbers[ii];
+    }
+}
+
 double Iso::getHeaviestPeakMass() const
 {
     double mass = 0.0;
@@ -274,6 +285,17 @@ double Iso::getHeaviestPeakLProb() const
     return lprob;
 }
 
+void Iso::getHeaviestPeakSignature(int* space) const
+{
+    for (int ii = 0; ii < dimNumber; ii++)
+    {
+        size_t heaviest_idx = marginals[ii]->getHeaviestAtomIndex();
+        for(int jj = 0; jj < isotopeNumbers[ii]; jj++)
+            space[jj] = (jj == heaviest_idx ? atomCounts[ii] : 0);
+        space += isotopeNumbers[ii];
+    }
+}
+
 double Iso::getMonoisotopicPeakMass() const
 {
     double mass = 0.0;
@@ -290,6 +312,17 @@ double Iso::getMonoisotopicPeakLProb() const
     return lprob;
 }
 
+void Iso::getMonoisotopicPeakSignature(int* space) const
+{
+    for (int ii = 0; ii < dimNumber; ii++)
+    {
+        size_t monoisotopic_idx = marginals[ii]->getMonoisotopicAtomIndex();
+        for(int jj = 0; jj < isotopeNumbers[ii]; jj++)
+            space[jj] = (jj == monoisotopic_idx ? atomCounts[ii] : 0);
+        space += isotopeNumbers[ii];
+    }
+}
+
 double Iso::getUnlikeliestPeakLProb() const
 {
     double lprob = 0.0;

src/IsoSpec++/isoSpec++.h

@@ -138,12 +138,24 @@ class ISOSPEC_EXPORT_SYMBOL Iso {
     //! Get the log-probability of the lightest peak in the isotopic distribution.
     double getLightestPeakLProb() const;
 
+    /*!
+        Write the counts of isotopes in the lightest peak into target memory location. It must be large enough to accomodate it.
+        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
+    */
+    void getLightestPeakSignature(int* space) 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;
 
+    /*!
+        Write the counts of isotopes in the heaviest peak into target memory location. It must be large enough to accomodate it.
+        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
+    */
+    void getHeaviestPeakSignature(int* space) 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
@@ -158,6 +170,12 @@ class ISOSPEC_EXPORT_SYMBOL Iso {
     */
     double getMonoisotopicPeakLProb() const;
 
+    /*!
+        Write the counts of isotopes in the monoisotopic peak into target memory location. It must be large enough to accomodate it.
+        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
+    */
+    void getMonoisotopicPeakSignature(int* space) 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

@@ -259,14 +259,7 @@ double Marginal::getLightestConfMass() const
 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;
-        }
+    size_t lightest_idx = getLightestAtomIndex();
     for(unsigned int ii = 0; ii < isotopeNo; ii++)
         lightest_conf[ii] = (ii == lightest_idx ? atomCnt : 0);
     return logProb(lightest_conf.get());
@@ -284,14 +277,7 @@ double Marginal::getHeaviestConfMass() const
 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;
-        }
+    size_t heaviest_idx = getHeaviestAtomIndex();
     for(unsigned int ii = 0; ii < isotopeNo; ii++)
         heaviest_conf[ii] = (ii == heaviest_idx ? atomCnt : 0);
     return logProb(heaviest_conf.get());
@@ -310,6 +296,32 @@ size_t Marginal::getMonoisotopicAtomIndex() const
     return found_idx;
 }
 
+size_t Marginal::getLightestAtomIndex() const
+{
+    double found_mass = std::numeric_limits<double>::infinity();
+    size_t found_idx = 0;
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        if( found_mass > atom_masses[ii] )
+        {
+            found_mass = atom_masses[ii];
+            found_idx = ii;
+        }
+    return found_idx;
+}
+
+size_t Marginal::getHeaviestAtomIndex() const
+{
+    double found_mass = -std::numeric_limits<double>::infinity();
+    size_t found_idx = 0;
+    for(unsigned int ii = 0; ii < isotopeNo; ii++)
+        if( found_mass < atom_masses[ii] )
+        {
+            found_mass = atom_masses[ii];
+            found_idx = ii;
+        }
+    return found_idx;
+}
+
 double Marginal::getMonoisotopicConfMass() const
 {
     size_t idx = getMonoisotopicAtomIndex();

src/IsoSpec++/marginalTrek++.h

@@ -124,6 +124,18 @@ class Marginal
     */
     size_t getMonoisotopicAtomIndex() const;
 
+    //! Get the index of the lightest atom.
+    /*!
+        \return The index of the lightest isotope of the element.
+    */
+    size_t getLightestAtomIndex() const;
+
+    //! Get the index of the heaviest atom.
+    /*!
+        \return The index of the heaviest isotope of the element.
+    */
+    size_t getHeaviestAtomIndex() const;
+    
     //! The the mass of the mode subisotopologue.
     /*!
         \return The mass of one of the most probable subisotopologues.

src/IsoSpecPy/IsoSpecPy.py

@@ -202,6 +202,8 @@ def __init__(self, formula="",
             offsets.append(tuple(newl))
         self.offsets = tuple(offsets)
 
+        self.conf_space = isoFFI.ffi.new("int[" + str(sum(self.isotopeNumbers)) + "]")
+
         self.iso = self.ffi.setupIso(self.dimNumber, self.isotopeNumbers,
                                      self.atomCounts,
                                      [i/charge for s in self.isotopeMasses for i in s],
@@ -223,6 +225,11 @@ def getLightestPeakLProb(self):
         """Get the log probability of the lightest peak in the isotopic distribution."""
         return self.ffi.getLightestPeakLProbIso(self.iso)
 
+    def getLightestPeakConf(self):
+        """Get the isotopic configuration of the lightest peak in the isotopic distribution."""
+        self.ffi.getLightestPeakSignature(self.iso, self.conf_space)
+        return self.parse_conf(self.conf_space)
+
     def getHeaviestPeakMass(self):
         """Get the heaviest peak in the isotopic distribution."""
         return self.ffi.getHeaviestPeakMassIso(self.iso)
@@ -231,6 +238,11 @@ def getHeaviestPeakLProb(self):
         """Get the log probability of the heaviest peak in the isotopic distribution."""
         return self.ffi.getHeaviestPeakLProbIso(self.iso)
 
+    def getHeaviestPeakConf(self):
+        """Get the isotopic configuration of the heaviest peak in the isotopic distribution."""
+        self.ffi.getHeaviestPeakSignature(self.iso, self.conf_space)
+        return self.parse_conf(self.conf_space)
+
     def getMonoisotopicPeakMass(self):
         """Get the monoisotopic mass of the peak."""
         return self.ffi.getMonoisotopicPeakMassIso(self.iso)
@@ -239,6 +251,11 @@ def getMonoisotopicPeakLProb(self):
         """Get the log probability of the monoisotopic peak in the isotopic distribution."""
         return self.ffi.getMonoisotopicPeakLProbIso(self.iso)
 
+    def getMonoisotopicPeakConf(self):
+        """Get the isotopic configuration of the monoisotopic peak in the isotopic distribution."""
+        self.ffi.getMonoisotopicPeakSignature(self.iso, self.conf_space)
+        return self.parse_conf(self.conf_space)
+
     def getModeLProb(self):
         """Get the log probability of the most probable peak(s) in the isotopic distribution."""
         return self.ffi.getModeLProbIso(self.iso)
@@ -700,7 +717,6 @@ def __init__(self, formula="", get_confs=False, **kwargs):
         """
         self.cgen = None
         super(IsoGenerator, self).__init__(formula=formula, get_confs=get_confs, **kwargs)
-        self.conf_space = isoFFI.ffi.new("int[" + str(sum(self.isotopeNumbers)) + "]")
         self.firstuse = True
 
     def __iter__(self):

src/IsoSpecPy/isoFFI.py

@@ -23,10 +23,13 @@ def __init__(self):
 
         double getLightestPeakMassIso(void* iso);
         double getLightestPeakLProbIso(void* iso);
+        void getLightestPeakSignature(void* iso, int* space);
         double getHeaviestPeakMassIso(void* iso);
         double getHeaviestPeakLProbIso(void* iso);
+        void getHeaviestPeakSignature(void* iso, int* space);
         double getMonoisotopicPeakMassIso(void* iso);
         double getMonoisotopicPeakLProbIso(void* iso);
+        void getMonoisotopicPeakSignature(void* iso, int* space);
         double getModeLProbIso(void* iso);
         double getModeMassIso(void* iso);
         double getTheoreticalAverageMassIso(void* iso);

tests/Python/test_iface.py

@@ -173,11 +173,14 @@ def test_lightest_peak():
     iso = IsoSpecPy.Iso(formula=formula)
     lightest_mass = iso.getLightestPeakMass()
     lightest_lprob = iso.getLightestPeakLProb()
-    iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula)
+    lightest_conf = iso.getLightestPeakConf()
+    iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula, get_confs=True)
     masses = list(iso_threshold.masses)
     probs = list(iso_threshold.probs)
+    confs = list(iso_threshold.confs)
     min_index = masses.index(min(masses))
-    print(lightest_mass, lightest_lprob, end=" ")
+    print(lightest_conf, lightest_mass, lightest_lprob, end=" ")
+    assert lightest_conf == confs[min_index]
     assert lightest_mass == masses[min_index]
     assert math.isclose(lightest_lprob, math.log(probs[min_index]), rel_tol=1e-9)
     print("OK!")
@@ -188,11 +191,14 @@ def test_heaviest_peak():
     iso = IsoSpecPy.Iso(formula=formula)
     heaviest_mass = iso.getHeaviestPeakMass()
     heaviest_lprob = iso.getHeaviestPeakLProb()
-    iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula)
+    heaviest_conf = iso.getHeaviestPeakConf()
+    iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula, get_confs=True)
     masses = list(iso_threshold.masses)
     probs = list(iso_threshold.probs)
+    confs = list(iso_threshold.confs)
     max_index = masses.index(max(masses))
-    print(heaviest_mass, heaviest_lprob, end=" ")
+    print(heaviest_conf, heaviest_mass, heaviest_lprob, end=" ")
+    assert heaviest_conf == confs[max_index]
     assert heaviest_mass == masses[max_index]
     assert math.isclose(heaviest_lprob, math.log(probs[max_index]), rel_tol=1e-9)
     print("OK!")
@@ -203,13 +209,16 @@ def test_monoisotopic_peak():
     iso = IsoSpecPy.Iso(formula=formula)
     monoisotopic_mass = iso.getMonoisotopicPeakMass()
     monoisotopic_lprob = iso.getMonoisotopicPeakLProb()
+    monoisotopic_conf = iso.getMonoisotopicPeakConf()
     iso_threshold = IsoSpecPy.IsoThreshold(0.0, formula=formula, get_confs=True)
     masses = list(iso_threshold.masses)
     probs = list(iso_threshold.probs)
     confs = list(iso_threshold.confs)
-    monoisotopic_conf = ((10, 0), (0, 1000), (10, 0), (0,0,0,0,0,0,0,1,0,0))
+    monoisotopic_conf2 = ((10, 0), (0, 1000), (10, 0), (0,0,0,0,0,0,0,1,0,0))
     monoisotopic_peak_idx = confs.index(monoisotopic_conf)
-    print(monoisotopic_mass, monoisotopic_lprob, end=" ")
+    print(monoisotopic_conf, monoisotopic_mass, monoisotopic_lprob, end=" ")
+    assert monoisotopic_conf == confs[monoisotopic_peak_idx]
+    assert monoisotopic_conf == monoisotopic_conf2
     assert monoisotopic_mass == masses[monoisotopic_peak_idx]
     assert math.isclose(monoisotopic_lprob, math.log(probs[monoisotopic_peak_idx]), rel_tol=1e-9)
     print("OK!")