isoSpec++.h

 
#pragma once
 
#include <unordered_map>
#include <queue>
#include <limits>
#include <string>
#include <vector>
#include <algorithm>
#include "platform.h"
#include "dirtyAllocator.h"
#include "summator.h"
#include "operators.h"
#include "marginalTrek++.h"
 
 
 
namespace IsoSpec
{
 
// This function is NOT guaranteed to be secure against malicious input. It should be used only for debugging.
unsigned int parse_formula(const char* formula,
                           std::vector<double>& isotope_masses,
                           std::vector<double>& isotope_probabilities,
                           int** isotopeNumbers,
                           int** atomCounts,
                           unsigned int* confSize,
                           bool use_nominal_masses = false);
 


class ISOSPEC_EXPORT_SYMBOL Iso {
 private:

    void setupMarginals(const double* _isotopeMasses,
                        const double* _isotopeProbabilities);
    bool            disowned;       
 
 protected:
    int             dimNumber;      
    int*            isotopeNumbers; 
    int*            atomCounts;     
    unsigned int    confSize;       
    int             allDim;         
    Marginal**      marginals;      
 
    bool doMarginalsNeedSorting() const;
 
 public:
    Iso();


    Iso(
        int             _dimNumber,
        const int*      _isotopeNumbers,
        const int*      _atomCounts,
        const double*   _isotopeMasses,
        const double*   _isotopeProbabilities
    );
    Iso(
        int             _dimNumber,
        const int*      _isotopeNumbers,
        const int*      _atomCounts,
        const double* const *  _isotopeMasses,
        const double* const *  _isotopeProbabilities
    );

    Iso(const char* formula, bool use_nominal_masses = false);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    inline Iso(const std::string& formula, bool use_nominal_masses = false) : Iso(formula.c_str(), use_nominal_masses) {}  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion


    static Iso FromFASTA(const char* fasta, bool use_nominal_masses = false, bool add_water = true);

    static inline Iso FromFASTA(const std::string& fasta, bool use_nominal_masses = false, bool add_water = true) { return FromFASTA(fasta.c_str(), use_nominal_masses, add_water); }

    Iso(Iso&& other);
 
    /* We're not exactly following standard copy and assign semantics with Iso objects, so delete the default assign constructor just in case, so noone tries to use it. Copy ctor declared below. */
    Iso& operator=(const Iso& other) = delete;


    Iso(const Iso& other, bool fullcopy);

    virtual ~Iso();

    double getLightestPeakMass() const;

    double getHeaviestPeakMass() const;

    double getMonoisotopicPeakMass() const;

    double getModeLProb() const;

    double getUnlikeliestPeakLProb() const;

    double getModeMass() const;

    double getTheoreticalAverageMass() const;

    double variance() const;

    double stddev() const { return sqrt(variance()); }

    inline int getDimNumber() const { return dimNumber; }

    inline int getAllDim() const { return allDim; }

    void addElement(int atomCount, int noIsotopes, const double* isotopeMasses, const double* isotopeProbabilities);

    void saveMarginalLogSizeEstimates(double* priorities, double target_total_prob) const;
};
 


class ISOSPEC_EXPORT_SYMBOL IsoGenerator : public Iso
{
 public:
    const double mode_lprob;
 
 protected:
    double* partialLProbs;  
    double* partialMasses;  
    double* partialProbs;   
 
 public:

    virtual bool advanceToNextConfiguration() = 0;


    virtual double lprob() const { return partialLProbs[0]; }


    virtual double mass()  const { return partialMasses[0]; }


    virtual double prob() const { return partialProbs[0]; }

    virtual void get_conf_signature(int* space) const = 0;

    IsoGenerator(Iso&& iso, bool alloc_partials = true);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    virtual ~IsoGenerator();
};
 
 


template<typename MarginalType>
class ISOSPEC_EXPORT_SYMBOL IsoOrderedGeneratorTemplate: public IsoGenerator
{
 private:
    MarginalType**              marginalResults;                    
    std::priority_queue<void*, pod_vector<void*>, ConfOrder> pq;   
    void*                       topConf;                            
    DirtyAllocator              allocator;                          
    const pod_vector<double>**  logProbs;                           
    const pod_vector<double>**  masses;                             
    double                      currentLProb;                       
    double                      currentMass;                        
    double                      currentProb;                        
    int                         ccount;
 
 public:
    IsoOrderedGeneratorTemplate(const IsoOrderedGeneratorTemplate& other) = delete;
    IsoOrderedGeneratorTemplate& operator=(const IsoOrderedGeneratorTemplate& other) = delete;
 
    bool advanceToNextConfiguration() override final;


    inline void get_conf_signature(int* space) const override final
    {
        if constexpr (std::is_same<MarginalType, MarginalTrek>::value)
        {
            int* c = getConf(topConf);
 
            if (ccount >= 0)
                c[ccount]--;
 
            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResults[ii]->confs()[c[ii]], isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
 
            if (ccount >= 0)
                c[ccount]++;
        }
        else
            throw std::runtime_error("IsoOrderedGeneratorTemplate::get_conf_signature() called on a non-MarginalTrek generator. This is not supported yet.");
    };

    IsoOrderedGeneratorTemplate(Iso&& iso, int _tabSize  = 1000, int _hashSize = 1000);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    virtual ~IsoOrderedGeneratorTemplate();
 
    inline void get_conf_by_indexes(int* space)
    {
        if constexpr (std::is_same<MarginalType, SingleAtomMarginal<false>>::value)
        {
            if(dimNumber == 0)
                return;
 
            int* c = getConf(topConf);
            space[0] = std::max(c[0]-1, 0);
 
            for(int ii = 1; ii < dimNumber; ii++)
                space[ii] = c[ii];
        }
    }
};
 
using IsoOrderedGenerator = IsoOrderedGeneratorTemplate<MarginalTrek>;
 


class ISOSPEC_EXPORT_SYMBOL IsoThresholdGenerator: public IsoGenerator
{
 private:
    int*                    counter;            
    double*                 maxConfsLPSum;
    const double            Lcutoff;            
    PrecalculatedMarginal** marginalResults;
    PrecalculatedMarginal** marginalResultsUnsorted;
    int* marginalOrder;
 
    const double* lProbs_ptr;
    const double* lProbs_ptr_start;
    double* partialLProbs_second;
    double partialLProbs_second_val, lcfmsv;
    bool empty;
 
 public:
    IsoThresholdGenerator(const IsoThresholdGenerator& other) = delete;
    IsoThresholdGenerator& operator=(const IsoThresholdGenerator& other) = delete;
 
    inline void get_conf_signature(int* space) const override final
    {
        counter[0] = lProbs_ptr - lProbs_ptr_start;
        if(marginalOrder != nullptr)
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                int jj = marginalOrder[ii];
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
        else
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
    };


    IsoThresholdGenerator(Iso&& iso, double _threshold, bool _absolute = true, int _tabSize = 1000, int _hashSize = 1000, bool reorder_marginals = true);
 
    ~IsoThresholdGenerator();
 
    // Perform highly aggressive inling as this function is often called as while(advanceToNextConfiguration()) {}
    // which leads to an extremely tight loop and some compilers miss this (potentially due to the length of the function).
    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        lProbs_ptr++;
 
        if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
        {
            return true;
        }
 
        // If we reached this point, a carry is needed
 
        int idx = 0;
        lProbs_ptr = lProbs_ptr_start;
 
        int * cntr_ptr = counter;
 
        while(idx < dimNumber-1)
        {
            // counter[idx] = 0;
            *cntr_ptr = 0;
            idx++;
            cntr_ptr++;
            // counter[idx]++;
            (*cntr_ptr)++;
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            if(partialLProbs[idx] + maxConfsLPSum[idx-1] >= Lcutoff)
            {
                partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
                partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
                recalc(idx-1);
                return true;
            }
        }
 
        terminate_search();
        return false;
    }
 
 
    ISOSPEC_FORCE_INLINE double lprob() const override final { return partialLProbs_second_val + (*(lProbs_ptr)); }
    ISOSPEC_FORCE_INLINE double mass()  const override final { return partialMasses[1] + marginalResults[0]->get_mass(lProbs_ptr - lProbs_ptr_start); }
    ISOSPEC_FORCE_INLINE double prob()  const override final { return partialProbs[1] * marginalResults[0]->get_prob(lProbs_ptr - lProbs_ptr_start); }

    void terminate_search();

    void reset();

    size_t count_confs();
 
 private:
    ISOSPEC_FORCE_INLINE void recalc(int idx)
    {
        for(; idx > 0; idx--)
        {
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
            partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
        }
        partialLProbs_second_val = *partialLProbs_second;
        partialLProbs[0] = *partialLProbs_second + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = Lcutoff - partialLProbs_second_val;
    }
 
    ISOSPEC_FORCE_INLINE void short_recalc(int idx)
    {
        for(; idx > 0; idx--)
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
        partialLProbs_second_val = *partialLProbs_second;
        partialLProbs[0] = *partialLProbs_second + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = Lcutoff - partialLProbs_second_val;
    }
};
 
 
 
 
template<typename MarginalType>
class ISOSPEC_EXPORT_SYMBOL IsoLayeredGeneratorTemplate : public IsoGenerator
{
 private:
    int*                    counter;            
    double*                 maxConfsLPSum;
    double currentLThreshold, lastLThreshold;
    MarginalType** marginalResults;
    MarginalType** marginalResultsUnsorted;
    int* marginalOrder;
 
    const double* lProbs_ptr;
    const double* lProbs_ptr_start;
    const double** resetPositions;
    double* partialLProbs_second;
    double partialLProbs_second_val, lcfmsv, last_lcfmsv;
    bool marginalsNeedSorting;
 
 
 public:
    IsoLayeredGeneratorTemplate(const IsoLayeredGeneratorTemplate& other) = delete;
    IsoLayeredGeneratorTemplate& operator=(const IsoLayeredGeneratorTemplate& other) = delete;
 
    inline void get_conf_signature(int* space) const override final
    {
        counter[0] = lProbs_ptr - lProbs_ptr_start;
        if(marginalOrder != nullptr)
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                int jj = marginalOrder[ii];
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
        else
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
    };
 
    inline double get_currentLThreshold() const { return currentLThreshold; }
 
    IsoLayeredGeneratorTemplate(Iso&& iso, int _tabSize = 1000, int _hashSize = 1000, bool reorder_marginals = true, double t_prob_hint = 0.99);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion
 
    ~IsoLayeredGeneratorTemplate();
 
    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        do
        {
            if(advanceToNextConfigurationWithinLayer())
                return true;
        } while(IsoLayeredGeneratorTemplate<MarginalType>::nextLayer(-2.0));
        return false;
    }
 
    ISOSPEC_FORCE_INLINE bool advanceToNextConfigurationWithinLayer()
    {
        do{
            lProbs_ptr++;
 
            if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
                return true;
        }
        while(carry());  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that's not an empty loop body, that's a do{...}while(...) construct
        return false;
    }
 
    ISOSPEC_FORCE_INLINE double lprob() const override final { return partialLProbs_second_val + (*(lProbs_ptr)); };
    ISOSPEC_FORCE_INLINE double mass()  const override final { return partialMasses[1] + marginalResults[0]->get_mass(lProbs_ptr - lProbs_ptr_start); };
    ISOSPEC_FORCE_INLINE double prob()  const override final { return partialProbs[1] * marginalResults[0]->get_prob(lProbs_ptr - lProbs_ptr_start); };

    void terminate_search();
 

    ISOSPEC_FORCE_INLINE void recalc(int idx)
    {
        for(; idx > 0; idx--)
        {
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
            partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
        }
        partialLProbs_second_val = *partialLProbs_second;
        partialLProbs[0] = partialLProbs_second_val + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = currentLThreshold - partialLProbs_second_val;
        last_lcfmsv = lastLThreshold - partialLProbs_second_val;
    }
 
    bool nextLayer(double offset);
 
    void get_conf_by_indexes(int* space) const
    {
        if constexpr (std::is_same<MarginalType, SingleAtomMarginal<true>>::value)
        {
            counter[0] = lProbs_ptr - lProbs_ptr_start;
            if(marginalOrder != nullptr)
            {
                for(int ii = 0; ii < dimNumber; ii++)
                {
                    int jj = marginalOrder[ii];
                    space[ii] = marginalResultsUnsorted[ii]->get_original_position(counter[jj]);
                }
            }
            else
            {
                for(int ii = 0; ii < dimNumber; ii++)
                    space[ii] = marginalResultsUnsorted[ii]->get_original_position(counter[ii]);
            }
        }
        else
            throw std::runtime_error("IsoLayeredGeneratorTemplate::get_conf_by_indexes() called on a non-SingleAtomMarginal generator. This is not supported yet.");
    }
 
 private:
    bool carry();
};
using IsoLayeredGenerator = IsoLayeredGeneratorTemplate<LayeredMarginal>;
 
template<typename IsoType>
class IsoStochasticGeneratorTemplate : public IsoGenerator
{
    IsoType ILG;
    size_t to_sample_left;
    const double precision;
    const double beta_bias;
    double confs_prob;
    double chasing_prob;
    size_t current_count;
    std::mt19937& rdvariate_gen;  
 
 public:
    IsoStochasticGeneratorTemplate(Iso&& iso, size_t no_molecules, double precision = 0.9999, double beta_bias = 5.0, std::mt19937& rdvariate_gen = random_gen);
 
    ISOSPEC_FORCE_INLINE size_t count() const { return current_count; }
 
    ISOSPEC_FORCE_INLINE double mass() const override final { return ILG.mass(); }
 
    ISOSPEC_FORCE_INLINE double prob() const override final { return static_cast<double>(count()); }
 
    ISOSPEC_FORCE_INLINE double lprob() const override final { return log(prob()); }
 
    ISOSPEC_FORCE_INLINE void get_conf_signature(int* space) const override final { ILG.get_conf_signature(space); }
 
    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        /* This function will be used mainly in very small, tight loops, therefore it makes sense to
         * aggressively inline it, despite its seemingly large body.
         */
        while(true)
        {
            double curr_conf_prob_left, current_prob;
 
            if(to_sample_left <= 0)
                return false;
 
            if(confs_prob < chasing_prob)
            {
                // Beta was last
                current_count = 1;
                to_sample_left--;
                if(!ILG.advanceToNextConfiguration())
                    return false;
                current_prob = ILG.prob();
                confs_prob += current_prob;
                while(confs_prob <= chasing_prob)
                {
                    if(!ILG.advanceToNextConfiguration())
                        return false;
                    current_prob = ILG.prob();
                    confs_prob += current_prob;
                }
                if(to_sample_left <= 0)
                    return true;
                curr_conf_prob_left = confs_prob - chasing_prob;
            }
            else
            {
                // Binomial was last
                current_count = 0;
                if(!ILG.advanceToNextConfiguration())
                    return false;
                current_prob = ILG.prob();
                confs_prob += current_prob;
                curr_conf_prob_left = current_prob;
            }
 
            double prob_left_to_1 = precision - chasing_prob;
            double expected_confs = curr_conf_prob_left * to_sample_left / prob_left_to_1;
 
            if(expected_confs <= beta_bias)
            {
                // Beta mode: we keep making beta jumps until we leave the current configuration
                chasing_prob += rdvariate_beta_1_b(to_sample_left, rdvariate_gen) * prob_left_to_1;
                while(chasing_prob <= confs_prob)
                {
                    current_count++;
                    to_sample_left--;
                    if(to_sample_left == 0)
                        return true;
                    prob_left_to_1 = precision - chasing_prob;
                    chasing_prob += rdvariate_beta_1_b(to_sample_left, rdvariate_gen) * prob_left_to_1;
                }
                if(current_count > 0)
                    return true;
            }
            else
            {
                // Binomial mode: a single binomial step
                size_t rbin = rdvariate_binom(to_sample_left, curr_conf_prob_left/prob_left_to_1, rdvariate_gen);
                current_count += rbin;
                to_sample_left -= rbin;
                chasing_prob = confs_prob;
                if(current_count > 0)
                    return true;
            }
        };
    }
 
    ISOSPEC_FORCE_INLINE void get_indexes(int* space)
    {
        ILG.get_conf_by_indexes(space);
    }
};
 
using IsoStochasticGenerator = IsoStochasticGeneratorTemplate<IsoLayeredGenerator>;
 
}  // namespace IsoSpec