marginalTrek++.h

/*
 *   Copyright (C) 2015-2020 Mateusz Łącki and Michał Startek.
 *
 *   This file is part of IsoSpec.
 *
 *   IsoSpec is free software: you can redistribute it and/or modify
 *   it under the terms of the Simplified ("2-clause") BSD licence.
 *
 *   IsoSpec is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 *   You should have received a copy of the Simplified BSD Licence
 *   along with IsoSpec.  If not, see <https://opensource.org/licenses/BSD-2-Clause>.
 */
 
#pragma once
 
#include <queue>
#include <algorithm>
#include <vector>
#include <functional>
#include <utility>
#include <memory>
#include "conf.h"
#include "allocator.h"
#include "operators.h"
#include "summator.h"
#include "pod_vector.h"
 
 
namespace IsoSpec
{


class Marginal
{
 private:
    bool disowned;
 protected:
    const unsigned int isotopeNo;       
    const unsigned int atomCnt;         
    const double* const atom_lProbs;    
    const double* const atom_masses;    
    const double loggamma_nominator;    
    Conf mode_conf;               
    double mode_lprob;            
 
 
 public:

    Marginal(
        const double* _masses,   // masses size = logProbs size = isotopeNo
        const double* _probs,
        int _isotopeNo,
        int _atomCnt
    );
 
    // Get rid of the C++ generated assignment constructor.
    Marginal& operator= (const Marginal& other) = delete;

    Marginal(const Marginal& other);

    Marginal(Marginal&& other);

    virtual ~Marginal();


    inline int get_isotopeNo() const { return isotopeNo; }
 
    inline const double* get_lProbs() const { return atom_lProbs; }


    double getLightestConfMass() const;


    double getHeaviestConfMass() const;


    double getMonoisotopicConfMass() const;


    size_t getMonoisotopicAtomIndex() const;


 
    inline double getModeMass() { ensureModeConf(); return calc_mass(mode_conf, atom_masses, isotopeNo); }


    inline double getModeLProb() { ensureModeConf(); return mode_lprob; }

    inline double fastGetModeLProb() { return mode_lprob; }


//    inline double getModeProb() const { return exp(getModeLProb()); }

    Conf computeModeConf() const;


    inline double getSmallestLProb() const { return atomCnt * *std::min_element(atom_lProbs, atom_lProbs+isotopeNo); }


    double getAtomAverageMass() const;


    inline double getTheoreticalAverageMass() const { return getAtomAverageMass() * atomCnt; }


 protected:
    ISOSPEC_FORCE_INLINE double unnormalized_logProb(Conf conf) const { double ret = 0.0; for(size_t ii = 0; ii < isotopeNo; ii++) ret += minuslogFactorial(conf[ii]) + conf[ii] * atom_lProbs[ii]; return ret; }
    ISOSPEC_FORCE_INLINE double logProb(Conf conf) const { return loggamma_nominator + unnormalized_logProb(conf); }
 public:
    double variance() const;

    double getLogSizeEstimate(double logEllipsoidRadius) const;
 
    inline void ensureModeConf() { if (mode_conf == nullptr) setupMode(); }
 private:
    void setupMode();
};
 

class MarginalTrek : public Marginal
{
 private:
    int current_count;
    const ConfOrderMarginal orderMarginal;
    std::priority_queue<ProbAndConfPtr, pod_vector<ProbAndConfPtr> > pq;
    pod_vector<unsafe_pod_vector<ProbAndConfPtr> > fringe;
    Allocator<int> allocator;
    pod_vector<double> _conf_lprobs;
    pod_vector<double> _conf_masses;
    pod_vector<int*> _confs;
 
    const double min_lprob;
    size_t current_bucket;
    size_t initialized_until;

    bool add_next_conf();
 
    size_t bucket_no(double lprob) { return static_cast<size_t>( (mode_lprob - lprob) * 100.0 ); }
 
 public:

    MarginalTrek(
        Marginal&& m,
        int tabSize = 1000,
        int hashSize = 1000
    );  // NOLINT(runtime/explicit) - Constructor deliberately left usable as a conversion.
 
    MarginalTrek(const MarginalTrek& other) = delete;
    MarginalTrek& operator=(const MarginalTrek& other) = delete;


    inline bool probeConfigurationIdx(int idx)
    {
        while(current_count <= idx)
            if(!add_next_conf())
                return false;
        return true;
    }


    inline double getModeLProb() const { return mode_lprob; }
 
 
    inline const pod_vector<double>& conf_lprobs() const { return _conf_lprobs; }
    inline const pod_vector<double>& conf_masses() const { return _conf_masses; }
    inline const pod_vector<Conf>& confs() const { return _confs; }
 
 
    virtual ~MarginalTrek();
};
 


class PrecalculatedMarginal : public Marginal
{
 protected:
    pod_vector<Conf> configurations;
    Conf* confs;
    unsigned int no_confs;
    double* masses;
    pod_vector<double> lProbs;
    double* probs;
    Allocator<int> allocator;
 public:

    PrecalculatedMarginal(
        Marginal&& m,
        double lCutOff,
        bool sort = true,
        int tabSize = 1000,
        int hashSize = 1000
    );
 
    PrecalculatedMarginal(const PrecalculatedMarginal& other) = delete;
    PrecalculatedMarginal& operator=(const PrecalculatedMarginal& other) = delete;

    virtual ~PrecalculatedMarginal();


    inline bool inRange(unsigned int idx) const { return idx < no_confs; }


    inline const double& get_lProb(int idx) const { return lProbs[idx]; }


    inline const double& get_prob(int idx) const { return probs[idx]; }


    inline const double& get_mass(int idx) const { return masses[idx]; }


    inline const double* get_lProbs_ptr() const { return lProbs.data(); }


    inline const double* get_masses_ptr() const { return masses; }
 


    inline const Conf& get_conf(int idx) const { return confs[idx]; }


    inline unsigned int get_no_confs() const { return no_confs; }


    inline double getModeLProb() const { return mode_lprob; }
};
 
 


class LayeredMarginal : public Marginal
{
 private:
    double current_threshold;
    pod_vector<Conf> configurations;
    pod_vector<Conf> fringe;
    pod_vector<double> fringe_unn_lprobs;
    Allocator<int> allocator;
    const ConfEqual equalizer;
    const KeyHasher keyHasher;
    pod_vector<double> lProbs;
    pod_vector<double> probs;
    pod_vector<double> masses;
    double* guarded_lProbs;
 
 public:

    LayeredMarginal(Marginal&& m, int tabSize = 1000, int hashSize = 1000);  // NOLINT(runtime/explicit) - constructor deliberately left usable as a conversion
 
    LayeredMarginal(const LayeredMarginal& other) = delete;
    LayeredMarginal& operator=(const LayeredMarginal& other) = delete;


    bool extend(double new_threshold, bool do_sort = true);

    inline double get_lProb(int idx) const { return guarded_lProbs[idx]; }  // access to idx == -1 is valid and gives a guardian of +inf

    inline double get_prob(int idx) const { return probs[idx]; }

    inline double get_mass(int idx) const { return masses[idx]; }

    inline const double* get_lProbs_ptr() const { return lProbs.data()+1; }

    inline const Conf& get_conf(int idx) const { return configurations[idx]; }

    inline unsigned int get_no_confs() const { return configurations.size(); }

    double get_min_mass() const;

    double get_max_mass() const;


    inline double getModeLProb() const { return mode_lprob; }
    const pod_vector<double>& conf_lprobs() const { return lProbs; }
    const pod_vector<double>& conf_masses() const { return masses; }
 
};
 
 
 
 
 
 
 
 
template <bool add_guards>
class SingleAtomMarginal : public Marginal
{
 private:
    double current_threshold;
    pod_vector<double> lProbs;
    pod_vector<double> probs;
    pod_vector<double> masses;
    pod_vector<size_t> original_indexes;
    double* guarded_lProbs;
    int extended_to_idx;
 
 public:

    SingleAtomMarginal(Marginal&& m, int tabSize = 1000, int hashSize = 1000);  // NOLINT(runtime/explicit) - constructor deliberately left usable as a conversion
 
    SingleAtomMarginal(const SingleAtomMarginal& other) = delete;
    SingleAtomMarginal& operator=(const SingleAtomMarginal& other) = delete;


    bool extend(double new_threshold, [[maybe_unused]] bool do_sort = true) {
        static_assert(add_guards, "SingleAtomMarginal::extend: add_guards must be true");
        current_threshold = new_threshold;
        bool extended = false;
        while(guarded_lProbs[extended_to_idx] >= current_threshold)
        {
            extended_to_idx++;
            extended = true;
        }
        return extended || (guarded_lProbs[extended_to_idx] != -std::numeric_limits<double>::infinity());
    };

    inline double get_lProb(int idx) const { return guarded_lProbs[idx]; }  // access to idx == -1 is valid and gives a guardian of +inf

    inline double get_prob(int idx) const { return probs[idx]; }

    inline double get_mass(int idx) const { return masses[idx]; }

    inline const double* get_lProbs_ptr() const { return lProbs.data()+1; }

    inline const Conf& get_conf([[maybe_unused]] int idx) const { throw std::logic_error("SingleAtomMarginal.get_conf: not implemented"); /*return configurations[idx];*/ }

    inline unsigned int get_no_confs() const {
        if constexpr(add_guards)
            return extended_to_idx;
        else
            return static_cast<unsigned int>(original_indexes.size());}

    inline double get_min_mass() const { throw std::logic_error("SingleAtomMarginal.get_min_mass: not implemented"); };

    double get_max_mass() const { throw std::logic_error("SingleAtomMarginal.get_max_mass: not implemented"); };


    inline double getModeLProb() const { return mode_lprob; }
 
    inline bool probeConfigurationIdx(int idx)
    {
        return idx < static_cast<int>(original_indexes.size());
    }
    const pod_vector<double>& conf_lprobs() const { return lProbs; }
    const pod_vector<double>& conf_masses() const { return masses; }
 
    int get_original_position(int idx) const
    {
        #ifdef ISOSPEC_DEBUG
        if (idx < 0 || static_cast<size_t>(idx) >= original_indexes.size())
            throw std::out_of_range("Index out of range in SingleAtomMarginal::et_original_position");
        #endif
        return original_indexes[idx];
    }
};
 
 
 
 
template <typename T>
class LoggingMarginal : public Marginal
{
 private:
    std::unique_ptr<T> real_marginal;
 public:
 
 
    LoggingMarginal(T&& m) : Marginal(m), real_marginal(std::make_unique(std::move(m))) {}
    LoggingMarginal(T&& m, int tabSize, int hashSize)
        : Marginal(m), real_marginal(std::make_unique<T>(std::move(m), tabSize, hashSize)) {}
 
    LoggingMarginal(const LoggingMarginal& other) = delete;
    LoggingMarginal& operator=(const LoggingMarginal& other) = delete;
 
    inline double getModeLProb() const { auto ret = real_marginal->getModeLProb(); std::cout << "LoggingMarginal::getModeLProb: " << ret << std::endl; return ret; }
    inline double get_lProb(int idx) const { auto ret = real_marginal->get_lProb(idx); std::cout << "LoggingMarginal::get_lProb: " << idx << " " << ret << std::endl; return ret; }
    inline double get_prob(int idx) const { auto ret = real_marginal->get_prob(idx); std::cout << "LoggingMarginal::get_prob: " << idx << " " << ret << std::endl; return ret; }
    inline double get_mass(int idx) const { auto ret = real_marginal->get_mass(idx); std::cout << "LoggingMarginal::get_mass: " << idx << " " << ret << std::endl; return ret; }
    inline const double* get_lProbs_ptr() const { auto ret = real_marginal->get_lProbs_ptr(); std::cout << "LoggingMarginal::get_lProbs_ptr: "; printArray<double>(ret, real_marginal->get_no_confs()); return ret; }
    inline const Conf& get_conf(int idx) const { auto ret = real_marginal->get_conf(idx); std::cout << "LoggingMarginal::get_conf: " << idx << std::endl; return ret; }
    inline unsigned int get_no_confs() const { auto ret = real_marginal->get_no_confs(); std::cout << "LoggingMarginal::get_no_confs: " << ret << std::endl; return ret; }
    inline bool probeConfigurationIdx(int idx) { auto ret = real_marginal->probeConfigurationIdx(idx); std::cout << "LoggingMarginal::probeConfigurationIdx: " << idx << " " << ret << std::endl; return ret; }
    inline double get_min_mass() const { auto ret = real_marginal->get_min_mass(); std::cout << "LoggingMarginal::get_min_mass: " << ret << std::endl; return ret; }
    inline double get_max_mass() const { auto ret = real_marginal->get_max_mass(); std::cout << "LoggingMarginal::get_max_mass: " << ret << std::endl; return ret; }
    inline double getModeMass() const { auto ret = real_marginal->getModeMass(); std::cout << "LoggingMarginal::getModeMass: " << ret << std::endl; return ret; }
    inline double getLightestConfMass() const { auto ret = real_marginal->getLightestConfMass(); std::cout << "LoggingMarginal::getLightestConfMass: " << ret << std::endl; return ret; }
    inline double getHeaviestConfMass() const { auto ret = real_marginal->getHeaviestConfMass(); std::cout << "LoggingMarginal::getHeaviestConfMass: " << ret << std::endl; return ret; }
    inline double getMonoisotopicConfMass() const { auto ret = real_marginal->getMonoisotopicConfMass(); std::cout << "LoggingMarginal::getMonoisotopicConfMass: " << ret << std::endl; return ret; }
    inline double getAtomAverageMass() const { auto ret = real_marginal->getAtomAverageMass(); std::cout << "LoggingMarginal::getAtomAverageMass: " << ret << std::endl; return ret; }
    inline double variance() const { auto ret = real_marginal->variance(); std::cout << "LoggingMarginal::variance: " << ret << std::endl; return ret; }
    inline double getTheoreticalAverageMass() const { auto ret = real_marginal->getTheoreticalAverageMass(); std::cout << "LoggingMarginal::getTheoreticalAverageMass: " << ret << std::endl; return ret; }
    inline double getLogSizeEstimate(double logEllipsoidRadius) const
    {
        auto ret = real_marginal->getLogSizeEstimate(logEllipsoidRadius);
        std::cout << "LoggingMarginal::getLogSizeEstimate: " << logEllipsoidRadius << " " << ret << std::endl;
        return ret;
    }
    inline void ensureModeConf() { real_marginal->ensureModeConf(); std::cout << "LoggingMarginal::ensureModeConf" << std::endl; }
    inline const pod_vector<double>& conf_lprobs() const { return real_marginal->conf_lprobs(); }
    inline const pod_vector<double>& conf_masses() const { return real_marginal->conf_masses(); }
    inline int get_original_position(int idx) const { return real_marginal->get_original_position(idx); }
    inline int get_isotopeNo() const { return real_marginal->get_isotopeNo(); }
    inline const double* get_lProbs() const { return real_marginal->get_lProbs(); }
    inline const double* get_masses() const { return real_marginal->get_masses(); }
    inline const double* get_atom_lProbs() const { return real_marginal->get_lProbs(); }
    inline const double* get_atom_masses() const { return real_marginal->get_masses(); }
    inline bool extend(double new_threshold, bool do_sort = true)
    {
        auto ret = real_marginal->extend(new_threshold, do_sort);
        std::cout << "LoggingMarginal::extend: " << new_threshold << " " << ret << std::endl;
        return ret;
    }
};
 
 
}  // namespace IsoSpec