isospec 0.1.0+2.6.0

/*
 *   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 <iostream>
#include <vector>
#include <functional>
#include <utility>
#include <memory>
#include "aligned_ptr.h"
#include "platform.h"
#include "iso_simd.h"   // ISOSPEC_LPROB_GUARDIANS: how many sentinels each marginal appends
#include "conf.h"
#include "allocator.h"
#include "operators.h"
#include "summator.h"
#include "pod_vector.h"


namespace IsoSpec
{

//! The marginal distribution class (a subisotopologue).
/*!
    This class mostly provides some basic common API for subclasses, but itself is not abstract.
    This class represents the probability distribution generated by one element only -- a subisotopologue.
    For instance, it might be the distribution of C200, that might be part of, say, C200H402.
    It corresponds to the multinomial distribution, where each configuration can also be attributed a precise mass.
    The constructor method perform initial hill-climbing to find the most probable sub-isotopologue (the mode).
*/
class Marginal
{
 private:
    bool disowned;

    //! Delegated-to constructor. Takes ownership of two already-allocated arrays
    //! via unique_ptr so that, should any initialization step throw, the arrays
    //! are released by the unique_ptr temporaries instead of being leaked. See
    //! the leak-on-OOM note for the public/copy constructors.
    Marginal(
        std::unique_ptr<const double[]> _atom_lProbs,
        std::unique_ptr<const double[]> _atom_masses,
        int _isotopeNo,
        int _atomCnt
    );
 protected:
    const unsigned int isotopeNo;       /*!< The number of isotopes of the given element. */
    const unsigned int atomCnt;         /*!< The number of atoms of the given element. */
    const double* const atom_lProbs;    /*!< Table of log-probabilities of all the isotopeNo isotopes. */
    const double* const atom_masses;    /*!< Table of atomic masses of all the isotopeNo isotopes. */
    const double loggamma_nominator;    /*!< The constant nominator that appears in the expressions for the multinomial probabilities. */
    Conf mode_conf;               /*!< A subisotopologue with most probability. If not unique, one of the representatives of that class of subisotopologues. */
    double mode_lprob;            /*!< The log-probability of the mode subisotopologue.*/


 public:
    //! Class constructor.
    /*!
        \param _masses A table of masses of the stable isotopes of the investigated element, e.g. for C10 it is 2: C12 and C13.
        \param _probs A table of natural frequencies of the stable isotopes of the investigated element, see IUPAC at https://iupac.org/isotopesmatter/
        \param _isotopeNo Number of isotopes of a given element.
        \param _atomCnt The number of atoms of the given element, e.g. 10 for C10.
        \return An instance of the Marginal class.
    */
    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;

    //! Copy constructor
    Marginal(const Marginal& other);

    //! Move constructor.
    Marginal(Marginal&& other);

    //! Destructor.
    virtual ~Marginal();

    //! Get the number of isotopes of the investigated element.
    /*!
        \return The integer number of isotopes of the investigated element.
    */
    inline int get_isotopeNo() const { return isotopeNo; }

    inline const double* get_lProbs() const { return atom_lProbs; }
    inline const double* get_masses() const { return atom_masses; }
    inline unsigned int get_atomCnt() const { return atomCnt; }

    //! Get the mass of the lightest subisotopologue.
    /*! This is trivially obtained by considering all atomNo atoms to be the lightest isotope possible.
        \return The mass of the lightiest subisotopologue.
    */
    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;

    //! The heaviest isotope's mass minus the lightest's: the most that moving one atom between two isotopes changes a configuration's mass (0 with no atoms, or one isotope).
    double getAtomMassSpan() 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.
    */
    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.
    */

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

    //! Get the log-probability of the mode subisotopologue.
    /*!
        \return The log-probability of a/the most probable subisotopologue.
    */
    inline double getModeLProb() const { return mode_lprob; }

    //! Same as getModeLProb(). Kept for source compatibility: the mode used to be computed lazily, and this was the unchecked accessor.
    inline double fastGetModeLProb() const { return mode_lprob; }

    //! The the probability of the mode subisotopologue.
    /*!
        \return The probability of a/the most probable subisotopologue.
    */
//    inline double getModeProb() const { return exp(getModeLProb()); }

    //! Computes and returns the mode configuration, a isotopeNo-large array that the caller is responsible for delete[]-ing.
    Conf computeModeConf() const;

    //! The the log-probability of the lightest subisotopologue.
    /*!
        \return The logarithm of the  smallest non-zero probability of a subisotopologue.
    */
    inline double getSmallestLProb() const { return atomCnt * *std::min_element(atom_lProbs, atom_lProbs+isotopeNo); }

    //! The average mass of a single atom.
    /*!
        \return The average mass of a single atom.
     */
    double getAtomAverageMass() const;

    //! The theoretical average mass of the molecule.
    /*!
        \return The theoretical average mass of the molecule.
     */
    inline double getTheoreticalAverageMass() const { return getAtomAverageMass() * atomCnt; }

    //! Calculate the log-probability of a given subisotopologue.
    /*!
        \param conf A subisotopologue (a table of integers describing subsequent isotope-counts).
        \return The log-probability of the input subisotopologue.
    */
 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:
    //! Calculate the variance of the theoretical distribution describing the subisotopologue
    double variance() const;

    //! Return estimated logarithm of size of the marginal at a given ellipsoid radius
    double getLogSizeEstimate(double logEllipsoidRadius) const;

    //! No-op, kept for source compatibility: the mode is now computed by the constructor.
    inline void ensureModeConf() const {}
};


//! The marginal distribution class (a subisotopologue).
class MarginalTrek : public Marginal
{
 private:
    int current_count;
    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;

    size_t current_bucket;
    size_t initialized_until;

    //! Proceed to the next configuration and memoize it (as it will be surely needed).
    bool add_next_conf();

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

 public:
    //! Move constructor: specializes the Marginal class.
    /*!
        \param tabSize The size of the table used to store configurations in the allocator.
        \param hashSize The size of the hash table used to store visited subisotopologues.
    */
    MarginalTrek(
        Marginal&& m,
        int tabSize = default_tab_size,
        int hashSize = 1000
    );  // NOLINT(runtime/explicit) - Constructor deliberately left usable as a conversion.

    MarginalTrek(const MarginalTrek& other) = delete;
    MarginalTrek& operator=(const MarginalTrek& other) = delete;

    //! Check if the table of computed subisotopologues does not have to be extended.
    /*!
        This function checks if the idx-th most probable subisotopologue was memoized and if not, computes it and memoizes it.

        \param idx The number of the idx-th most probable subisotopologue.
        \return Returns false if it the provided idx exceeds the total number of subisotopologues.
    */
    inline bool probeConfigurationIdx(int idx)
    {
        while(current_count <= idx)
            if(!add_next_conf())
                return false;
        return true;
    }

    //! Get the log-probability of the mode subisotopologue.
    /*!
        \return The log-probability of a/the most probable subisotopologue.
    */
    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();
};


//! Precalculated Marginal class
/*!
    This class serves to calculate a set of isotopologues that
    is defined by the minimal probability threshold.

    This works faster than if you did not know the threshold.
    If you have no idea about the threshold, you would need to call us,
    to change encode the layered version of the marginal.
*/
class PrecalculatedMarginal : public Marginal
{
 protected:
    pod_vector<Conf> configurations;
    Conf* confs;
    unsigned int no_confs;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> masses;
    pod_vector<double> lProbs;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> probs;
    Allocator<int> allocator;
 public:
    //! The move constructor (disowns the Marginal).
    /*!
        This constructor memoizes all subisotopologues with log-probability above the provided threshold lCutOff
        \param Marginal An instance of the Marginal class this class is about to disown.
        \param lCutOff The lower limit on the log-probability of the precomputed subisotopologues.
        \param sort Should the subisotopologues be stored with descending probability ?
        \return An instance of the PrecalculatedMarginal class.
    */
    PrecalculatedMarginal(
        Marginal&& m,
        double lCutOff,
        bool sort = true,
        int tabSize = default_tab_size,
        int hashSize = 1000
    );

    PrecalculatedMarginal(const PrecalculatedMarginal& other) = delete;
    PrecalculatedMarginal& operator=(const PrecalculatedMarginal& other) = delete;

    //! Destructor.
    virtual ~PrecalculatedMarginal();

    //! Is there a subisotopologue with a given number?
    /*!
        \return Returns true if idx does not exceed the number of pre-computed configurations.
    */
    inline bool inRange(unsigned int idx) const { return idx < no_confs; }

    //! Get the log-probability of the idx-th subisotopologue.
    /*!
        \param idx The number of the considered subisotopologue.
        \return The log-probability of the idx-th subisotopologue.
    */
    inline const double& get_lProb(int idx) const { return lProbs[idx]; }

    //! Get the probability of the idx-th subisotopologue.
    /*!
        \param idx The number of the considered subisotopologue.
        \return The probability of the idx-th subisotopologue.
    */
    inline const double& get_prob(int idx) const { return probs[idx]; }

    //! Get the mass of the idx-th subisotopologue.
    /*!
        \param idx The number of the considered subisotopologue.
        \return The mass of the idx-th subisotopologue.
    */
    inline const double& get_mass(int idx) const { return masses[idx]; }

    //! Get the table of the log-probabilities of subisotopologues.
    /*!
        \return Pointer to the first element in the table storing log-probabilities of subisotopologues.
    */
    inline const double* get_lProbs_ptr() const { return lProbs.data(); }

    //! Get the table of the masses of subisotopologues.
    /*!
        \return Pointer to the first element in the table storing masses of subisotopologues.
    */
    inline const double* get_masses_ptr() const { return masses.get(); }

    //! Get the table of the probabilities of subisotopologues.
    inline const double* get_probs_ptr() const { return probs.get(); }

    //! Get the counts of isotopes that define the subisotopologue.
    /*!
        \param idx The number of the considered subisotopologue.
        \return The counts of isotopes that define the subisotopologue.
    */
    inline const Conf& get_conf(int idx) const { return confs[idx]; }

    //! Get the number of precomputed subisotopologues.
    /*!
        \return The number of precomputed subisotopologues.
    */
    inline unsigned int get_no_confs() const { return no_confs; }

    //! Get the log-probability of the mode subisotopologue.
    /*!
        \return The log-probability of a/the most probable subisotopologue.
    */
    inline double getModeLProb() const { return mode_lprob; }
};



//! LayeredMarginal class
/*!
    An extendable version of the PrecalculatedMarginal, where you can extend the threshold at will.
*/
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;
    pod_vector<double> lProbs;
    pod_vector<double> probs;
    pod_vector<double> masses;
    double* guarded_lProbs;

 public:
    //! Move constructor: specializes the Marginal class.
    /*!
        \param tabSize The size of the table used to store configurations in the allocator.
        \param hashSize The size of the hash table used to store visited subisotopologues.
    */
    LayeredMarginal(Marginal&& m, int tabSize = default_tab_size, int hashSize = 1000);  // NOLINT(runtime/explicit) - constructor deliberately left usable as a conversion

    LayeredMarginal(const LayeredMarginal& other) = delete;
    LayeredMarginal& operator=(const LayeredMarginal& other) = delete;

    //! Extend the set of computed subisotopologues to those above the new threshold.
    /*!
        \param new_threshold The new log-probability limiting the subisotopologues from below.
        \return Returns false, if there are no fringe-subisotopologues (subisotopologues that were neighbours of the previously calculated subisotopologues, with log-probability below the previous threshold).
    */
    bool extend(double new_threshold, bool do_sort = true);

    //! get the log-probability of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_lProb.
    inline double get_lProb(int idx) const { return guarded_lProbs[idx]; }  // access to idx == -1 is valid and gives a guardian of +inf

    //! get the probability of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_eProb.
    inline double get_prob(int idx) const { return probs[idx]; }

    //! get the mass of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_mass.
    inline double get_mass(int idx) const { return masses[idx]; }

    //! Marginal masses and probabilities as raw arrays, indexed as get_lProbs_ptr() is -- for the batched kernels
    //! (isa_kernels.h). Invalidated, like get_lProbs_ptr(), by calls to extend().
    inline const double* get_masses_ptr() const { return masses.data(); }
    inline const double* get_probs_ptr() const { return probs.data(); }

    //! get the pointer to lProbs array. Accessing index -1 is legal and returns a guardian of +inf (trailing guardians of -inf follow the data). Warning: The pointer gets invalidated on calls to extend()
    inline const double* get_lProbs_ptr() const { return lProbs.data()+1; }

    //! get the counts of isotopes that define the subisotopologue, see details in @ref PrecalculatedMarginal::get_conf.
    inline const Conf& get_conf(int idx) const { return configurations[idx]; }

    //! Get the number of precomputed subisotopologues, see details in @ref PrecalculatedMarginal::get_no_confs.
    inline unsigned int get_no_confs() const { return configurations.size(); }

    //! Get the minimal mass in current layer
    double get_min_mass() const;

    //! Get the maximal mass in current layer
    double get_max_mass() const;

    //! Get the log-probability of the mode subisotopologue.
    /*!
        \return The log-probability of a/the most probable subisotopologue.
    */
    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:
    //! Move constructor: specializes the Marginal class.
    /*!
    */
    SingleAtomMarginal(Marginal&& m, int tabSize = default_tab_size, int hashSize = 1000);  // NOLINT(runtime/explicit) - constructor deliberately left usable as a conversion

    SingleAtomMarginal(const SingleAtomMarginal& other) = delete;
    SingleAtomMarginal& operator=(const SingleAtomMarginal& other) = delete;

    //! Extend the set of computed subisotopologues to those above the new threshold.
    /*!
        \param new_threshold The new log-probability limiting the subisotopologues from below.
        \return Returns false, if there are no fringe-subisotopologues (subisotopologues that were neighbours of the previously calculated subisotopologues, with log-probability below the previous threshold).
    */
    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());
    }

    //! get the log-probability of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_lProb.
    inline double get_lProb(int idx) const { return guarded_lProbs[idx]; }  // access to idx == -1 is valid and gives a guardian of +inf

    //! get the probability of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_eProb.
    inline double get_prob(int idx) const { return probs[idx]; }

    //! get the mass of the idx-th subisotopologue, see details in @ref PrecalculatedMarginal::get_mass.
    inline double get_mass(int idx) const { return masses[idx]; }

    //! Marginal masses and probabilities as raw arrays, indexed as get_lProbs_ptr() is -- for the batched kernels
    //! (isa_kernels.h). Invalidated, like get_lProbs_ptr(), by calls to extend().
    inline const double* get_masses_ptr() const { return masses.data(); }
    inline const double* get_probs_ptr() const { return probs.data(); }

    //! get the pointer to lProbs array. Accessing index -1 is legal and returns a guardian of +inf (trailing guardians of -inf follow the data). Warning: The pointer gets invalidated on calls to extend()
    inline const double* get_lProbs_ptr() const { return lProbs.data()+1; }

    //! get the counts of isotopes that define the subisotopologue, see details in @ref PrecalculatedMarginal::get_conf.
    inline const Conf& get_conf([[maybe_unused]] int idx) const { throw std::logic_error("SingleAtomMarginal.get_conf: not implemented"); /*return configurations[idx];*/ }

    //! Get the number of precomputed subisotopologues, see details in @ref PrecalculatedMarginal::get_no_confs.
    inline unsigned int get_no_confs() const {
        if constexpr(add_guards)
            return extended_to_idx;
        else
            return static_cast<unsigned int>(original_indexes.size());}

    //! Get the minimal mass in current layer
    inline double get_min_mass() const { throw std::logic_error("SingleAtomMarginal.get_min_mass: not implemented"); }

    //! Get the maximal mass in current layer
    double get_max_mass() const { throw std::logic_error("SingleAtomMarginal.get_max_mass: not implemented"); }

    //! Get the log-probability of the mode subisotopologue.
    /*!
        \return The log-probability of a/the most probable subisotopologue.
    */
    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::get_original_position");
        #endif
        return original_indexes[idx];
    }
};




template <typename T>
class LoggingMarginal : public Marginal
{
 private:
    std::unique_ptr<T> real_marginal;
 public:
    // The cast moves the Marginal base only: T itself is not move-constructible
    // (a bare std::move(m) would select T's deleted copy constructor), but every
    // marginal specialization constructs from a Marginal&&.
    explicit LoggingMarginal(T&& m) : Marginal(m), real_marginal(std::make_unique<T>(static_cast<Marginal&&>(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 { const Conf& 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