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>.
 */


#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <queue>
#include <utility>
#include <cstring>
#include <string>
#include <limits>
#include <memory>
#include "platform.h"
#include "marginalTrek++.h"
#include "conf.h"
#include "allocator.h"
#include "operators.h"
#include "summator.h"
#include "element_tables.h"
#include "misc.h"


namespace IsoSpec
{

//! Find one of the most probable subisotopologues.
/*!
    The algorithm uses the hill-climbing algorithm.
    It starts from a subisotopologue close to the mean of the underlying multinomial distribution.
    There might be more than one modes, in case of which this function will return only one of them, close to the mean.

    \param atomCnt

*/
void writeInitialConfiguration(const int atomCnt, const int isotopeNo, const double* lprobs, int* res)
{
    /*!
    Here we perform hill climbing to the mode of the marginal distribution (the subisotopologue distribution).
    We start from the point close to the mean of the underlying multinomial distribution.
    */

    // This approximates the mode (heuristics: the mean is close to the mode).
    for(int i = 0; i < isotopeNo; ++i)
        res[i] = static_cast<int>( atomCnt * exp(lprobs[i]) ) + 1;

    // The number of assigned atoms above.
    int s = 0;

    for(int i = 0; i < isotopeNo; ++i) s += res[i];

    int diff = atomCnt - s;

    // Too little: enlarging fist index.
    if( diff > 0 ){
        res[0] += diff;
    }
    // Too much: trying to redistribute the difference: hopefully the first element is the largest.
    if( diff < 0 ){
        diff = abs(diff);
        int i = 0;

        while( diff > 0){
            int coordDiff = res[i] - diff;

            if( coordDiff >= 0 ){
                res[i] -= diff;
                diff = 0;
            } else {
                res[i] = 0;
                i++;
                diff = abs(coordDiff);
            }
        }
    }

    // What we computed so far will be very close to the mode: hillclimb the rest of the way

    bool modified = true;
    double LP = unnormalized_logProb(res, lprobs, isotopeNo);
    double NLP;

    while(modified)
    {
        modified = false;
        for(int ii = 0; ii < isotopeNo; ii++)
            for(int jj = 0; jj < isotopeNo; jj++)
                if(ii != jj && res[ii] > 0)
                {
                    res[ii]--;
                    res[jj]++;
                    NLP = unnormalized_logProb(res, lprobs, isotopeNo);
                    if(NLP > LP || (NLP == LP && ii > jj))
                    {
                        modified = true;
                        LP = NLP;
                    }
                    else
                    {
                        res[ii]++;
                        res[jj]--;
                    }
                }
    }
}


namespace {

//! The element table's rows sorted by probability, so that getMLogProbs()
//! finds a probability's row by binary search rather than by scanning.
struct ProbabilityIndex {
    double prob[ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES];
    int row[ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES];

    ProbabilityIndex()
    {
        for(int ii = 0; ii < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; ii++)
            row[ii] = ii;
        std::sort(row, row + ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES,
                  [](int a, int b) { return elem_table_probability[a] < elem_table_probability[b]; });
        for(int ii = 0; ii < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES; ii++)
            prob[ii] = elem_table_probability[row[ii]];
    }

    //! Some row holding exactly this probability, or -1.
    int find(double p) const
    {
        const double* it = std::lower_bound(prob, prob + ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES, p);
        if(it != prob + ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES && *it == p)
            return row[it - prob];
        return -1;
    }
};

const ProbabilityIndex& probability_index()
{
    static const ProbabilityIndex instance;
    return instance;
}

}  // namespace

double* getMLogProbs(const double* probs, int isoNo)
{
    /*!
    Here we order the processor to round the numbers up rather than down.
    Rounding down could result in the algorithm falling in an infinite loop
    because of the numerical instability of summing.
    */
    if(isoNo < 1)
        throw std::invalid_argument("An element must have at least one isotope");
    // Written as !(in range) so that a NaN, for which every comparison is
    // false, is rejected too.
    double sum = 0.0;
    for(int ii = 0; ii < isoNo; ii++)
    {
        if(!(probs[ii] > 0.0 && probs[ii] <= 1.0))
            throw std::invalid_argument("All isotope probabilities p must fulfill: 0.0 < p <= 1.0");
        sum += probs[ii];
    }

    // An element's probabilities must sum to 1. Far off, the molecule's
    // probabilities are meaningless and numerically dangerous: a sum near 0
    // underflows every configuration to 0 (and normalising that gives NaN), a
    // sum above 1 overflows large molecules' to infinity ((1.0001)^(10^7) is
    // e^1000). So: rejected beyond kMaxProbabilitySumError -- loose enough for
    // abundances rounded to four decimals -- and divided by the sum beyond
    // kProbabilitySumRoundoff, below which the element table's rows (within
    // 2.2e-16 of 1) and anything else exact up to rounding pass untouched.
    constexpr double kMaxProbabilitySumError = 1e-3;
    constexpr double kProbabilitySumRoundoff = 1e-12;
    if(!(std::fabs(sum - 1.0) <= kMaxProbabilitySumError))
        throw std::invalid_argument("The isotope probabilities of an element must sum to 1 (within 1e-3)");
    std::unique_ptr<double[]> normalized;
    if(std::fabs(sum - 1.0) > kProbabilitySumRoundoff)
    {
        normalized.reset(new double[isoNo]);
        for(int ii = 0; ii < isoNo; ii++)
            normalized[ii] = probs[ii] / sum;
        probs = normalized.get();
    }

    double* ret = new double[isoNo];

    // A probability found in the element table takes the table's log (more
    // precise than libm's in about a third of the rows), anything else gets
    // log(). Rows of equal probability all carry the same log, so which
    // matching row is found does not matter. Isotopes coming from a formula
    // are consecutive table rows, so the row after the previous isotope's is
    // tried first, and only a miss there goes to the binary search.
    //
    // This used to be a linear scan of the whole table per isotope, with log()
    // computed up front whether needed or not, and it was the largest single
    // cost of parsing a formula -- and compiler-dependent: clang vectorised the
    // early-exit scan, gcc did not, which left gcc's formula parsing 1.5x
    // slower than clang's.
    const ProbabilityIndex& index = probability_index();
    int next_row = -1;
    for(int i = 0; i < isoNo; i++)
    {
        int found;
        if(next_row >= 0 && next_row < ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES && elem_table_probability[next_row] == probs[i])
            found = next_row;
        else
            found = index.find(probs[i]);
        ret[i] = found >= 0 ? elem_table_log_probability[found] : log(probs[i]);
        next_row = found >= 0 ? found + 1 : -1;
    }
    return ret;
}

double get_loggamma_nominator(int x)
{
    // calculate log gamma of the nominator calculated in the binomial exression.
    return -minuslogFactorial(x);
}

const double* checked_masses(const double* masses, int isoNo)
{
    if(isoNo < 1)
        throw std::invalid_argument("An element must have at least one isotope");
    for(int ii = 0; ii < isoNo; ii++)
        if(!(std::fabs(masses[ii]) <= max_abs_mass))
            throw std::invalid_argument("Isotope masses must be finite, with |mass| <= 1e100");
    return masses;
}

int verify_atom_cnt(int atomCnt)
{
    if(atomCnt < 0)
        throw std::invalid_argument("Atom counts must not be negative");
    #if !ISOSPEC_BUILDING_OPENMS
    if(ISOSPEC_G_FACT_TABLE_SIZE-1 <= atomCnt)
        throw std::length_error("Subisotopologue too large, size limit (that is, the maximum number of atoms of a single element in a molecule) is: " + std::to_string(ISOSPEC_G_FACT_TABLE_SIZE-1));
    #endif
    return atomCnt;
}

// Delegated-to constructor: the two arrays arrive already allocated and owned
// by unique_ptrs, so any throw during the rest of construction (e.g. from
// verify_atom_cnt) frees them via the unique_ptr parameters rather than leaking.
Marginal::Marginal(
    std::unique_ptr<const double[]> _atom_lProbs,
    std::unique_ptr<const double[]> _atom_masses,
    int _isotopeNo,
    int _atomCnt
) :
disowned(false),
isotopeNo(_isotopeNo),
atomCnt(verify_atom_cnt(_atomCnt)),
atom_lProbs(_atom_lProbs.release()),
atom_masses(_atom_masses.release()),
loggamma_nominator(get_loggamma_nominator(_atomCnt)),
mode_conf(nullptr)
// Deliberately not initializing mode_lprob: every constructor delegating here
// sets both mode members in its own body (see the public one below).
{}

Marginal::Marginal(
    const double* _masses,
    const double* _probs,
    int _isotopeNo,
    int _atomCnt
) :
// The two allocating expressions are wrapped in unique_ptr temporaries: if the
// second one throws, the first is freed during cleanup of the full-expression
// instead of being leaked (constructor bodies/init-lists don't free raw members
// on throw). The target constructor then takes ownership.
Marginal(std::unique_ptr<const double[]>(getMLogProbs(_probs, _isotopeNo)),
         std::unique_ptr<const double[]>(array_copy<double>(checked_masses(_masses, _isotopeNo), _isotopeNo)),
         _isotopeNo, _atomCnt)
{
    // The mode is computed here, eagerly, rather than on first request. A lazy
    // fill would be a write hidden behind const accessors (Iso::getModeLProb()
    // const reaches it through the marginals pointer), i.e. a data race for
    // two threads reading one shared Iso. Nearly every user of a marginal
    // needs the mode anyway -- every generator starts from it. Done in this
    // body, not in the delegated-to constructor, so that a throw from it runs
    // the destructor and frees the two arrays.
    mode_conf = computeModeConf();
    mode_lprob = logProb(mode_conf);
}

Marginal::Marginal(const Marginal& other) :
// Same exception-safety reasoning as the public constructor above. Because this
// delegates to the target constructor, once that completes the object is fully
// constructed, so if the mode_conf copy below throws the destructor runs and
// frees atom_lProbs/atom_masses too.
Marginal(std::unique_ptr<const double[]>(array_copy<double>(other.atom_lProbs, other.isotopeNo)),
         std::unique_ptr<const double[]>(array_copy<double>(other.atom_masses, other.isotopeNo)),
         other.isotopeNo, other.atomCnt)
{
    mode_conf = array_copy<int>(other.mode_conf, isotopeNo);
    mode_lprob = other.mode_lprob;
}


// the move-constructor: used in the specialization of the marginal.
Marginal::Marginal(Marginal&& other) :
disowned(other.disowned),
isotopeNo(other.isotopeNo),
atomCnt(other.atomCnt),
atom_lProbs(other.atom_lProbs),
atom_masses(other.atom_masses),
loggamma_nominator(other.loggamma_nominator),
mode_conf(other.mode_conf),
mode_lprob(other.mode_lprob)
{
    other.disowned = true;
}

Marginal::~Marginal()
{
    if(!disowned)
    {
        delete[] atom_masses;
        delete[] atom_lProbs;
        delete[] mode_conf;
    }
}


Conf Marginal::computeModeConf() const
{
    Conf res = new int[isotopeNo];
    writeInitialConfiguration(atomCnt, isotopeNo, atom_lProbs, res);
    return res;
}


double Marginal::getLightestConfMass() const
{
    double ret_mass = std::numeric_limits<double>::infinity();
    for(unsigned int ii = 0; ii < isotopeNo; ii++)
        if( ret_mass > atom_masses[ii] )
            ret_mass = atom_masses[ii];
    return ret_mass*atomCnt;
}

double Marginal::getLightestConfLProb() const
{
    std::unique_ptr<int[]> lightest_conf = std::make_unique<int[]>(isotopeNo);
    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());
}

double Marginal::getAtomMassSpan() const
{
    if(atomCnt == 0 || isotopeNo < 2)
        return 0.0;
    double lightest = atom_masses[0], heaviest = atom_masses[0];
    for(unsigned int ii = 1; ii < isotopeNo; ii++)
    {
        lightest = (std::min)(lightest, atom_masses[ii]);
        heaviest = (std::max)(heaviest, atom_masses[ii]);
    }
    return heaviest - lightest;
}

double Marginal::getHeaviestConfMass() const
{
    double ret_mass = 0.0;
    for(unsigned int ii = 0; ii < isotopeNo; ii++)
        if( ret_mass < atom_masses[ii] )
            ret_mass = atom_masses[ii];
    return ret_mass*atomCnt;
}

double Marginal::getHeaviestConfLProb() const
{
    std::unique_ptr<int[]> heaviest_conf = std::make_unique<int[]>(isotopeNo);
    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());
}

size_t Marginal::getMonoisotopicAtomIndex() const
{
    double found_prob = -std::numeric_limits<double>::infinity();
    size_t found_idx = 0;
    for(unsigned int ii = 0; ii < isotopeNo; ii++)
        if( found_prob < atom_lProbs[ii] )
        {
            found_prob = atom_lProbs[ii];
            found_idx = ii;
        }
    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();
    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;
    for(unsigned int ii = 0; ii < isotopeNo; ii++)
        ret += exp(atom_lProbs[ii]) * atom_masses[ii];
    return ret;
}

double Marginal::variance() const
{
    double ret = 0.0;
    double mean = getAtomAverageMass();
    for(size_t ii = 0; ii < isotopeNo; ii++)
    {
        double msq = atom_masses[ii] - mean;
        ret += exp(atom_lProbs[ii]) * msq * msq;
    }
    return ret * atomCnt;
}

double Marginal::getLogSizeEstimate(double logEllipsoidRadius) const
{
    if(isotopeNo <= 1)
        return -std::numeric_limits<double>::infinity();

    const double i = static_cast<double>(isotopeNo);
    const double k = i - 1.0;
    const double n = static_cast<double>(atomCnt);

    double sum_lprobs = 0.0;
    for(int jj = 0; jj < i; jj++)
        sum_lprobs += atom_lProbs[jj];

    // lgamma_half() takes twice its argument: (i+1)/2 is a half-integer for even i.
    const std::int64_t two_i = 2 * static_cast<std::int64_t>(isotopeNo);
    const std::int64_t two_n = 2 * static_cast<std::int64_t>(atomCnt);
    const double lgamma_i = lgamma_half(two_i);
    double log_V_simplex = k * log(n) - lgamma_i;
    double log_N_simplex = lgamma_half(two_n + two_i) - lgamma_half(two_n + 2) - lgamma_i;
    double log_V_ellipsoid = (k * (log(n) + logpi + logEllipsoidRadius) + sum_lprobs) * 0.5 - lgamma_half(static_cast<std::int64_t>(isotopeNo) + 1);

    return log_N_simplex + log_V_ellipsoid - log_V_simplex;
}


// this is roughly an equivalent of IsoSpec-Threshold-Generator
MarginalTrek::MarginalTrek(
    Marginal&& m,
    int tabSize,
    int
) :
Marginal(std::move(m)),
current_count(0),
pq(),
allocator(isotopeNo, tabSize)
{
    int* initialConf = allocator.makeCopy(mode_conf);

    pq.push({mode_lprob, initialConf});

    current_count = 0;

    fringe.resize_and_wipe(1);

    current_bucket = 0;
    initialized_until = 1;

    add_next_conf();
}


bool MarginalTrek::add_next_conf()
{
    /*!
    Add next configuration.
    If visited all, return false.
    */
    if(pq.empty())
    {
        current_bucket++;
        while(current_bucket < initialized_until && fringe[current_bucket].empty())
        {
//            std::cout << "EMPTY bucket, id: " << current_bucket << std::endl;
            current_bucket++;
        }

//        std::cout << "Entering bucket, size: " << fringe[current_bucket].size() << std::endl;

        if(current_bucket >= initialized_until)
            return false;

    //    std::cout << "Fringe size at pop: " << fringe[current_bucket].size() << std::endl;
        pq = std::priority_queue<ProbAndConfPtr, pod_vector<ProbAndConfPtr> >(std::less<ProbAndConfPtr>(), pod_vector<ProbAndConfPtr>(std::move(fringe[current_bucket])));
    };

    double logprob = pq.top().first;
    Conf topConf = pq.top().second;

    pq.pop();
    ++current_count;

    _confs.push_back(topConf);

    _conf_masses.push_back(calc_mass(topConf, atom_masses, isotopeNo));
    _conf_lprobs.push_back(logprob);

    for( unsigned int j = 0; j < isotopeNo; ++j )
    {
        if( topConf[j] > mode_conf[j])
            continue;

        if( topConf[j] > 0 )
        {
            for( unsigned int i = 0; i < isotopeNo; ++i )
            {
                if( topConf[i] < mode_conf[i] )
                    continue;
                // Growing index different than decreasing one AND
                // Remain on simplex condition.
                if( i != j ){
                    Conf acceptedCandidate = allocator.makeCopy(topConf);

                    ++acceptedCandidate[i];
                    --acceptedCandidate[j];

                    double new_prob = logProb(acceptedCandidate);
                    size_t bucket_nr = bucket_no(new_prob);

                    if(bucket_nr >= initialized_until)
                    {
                    //    std::cout << "Extending to: " << bucket_nr << std::endl;
                        initialized_until = bucket_nr+1;
                        fringe.resize_and_wipe(initialized_until);
                    }

                    ISOSPEC_IMPOSSIBLE(bucket_nr < current_bucket);
                    if(bucket_nr == current_bucket)
                        pq.push({new_prob, acceptedCandidate});
                    else
                        fringe[bucket_nr].push_back({new_prob, acceptedCandidate});
                }

                if( topConf[i] > mode_conf[i] )
                    break;
            }
        }
        if( topConf[j] < mode_conf[j] )
            break;
    }

    return true;
}


MarginalTrek::~MarginalTrek()
{
    const size_t fringe_size = fringe.size();
    for(size_t ii = 0; ii < fringe_size; ii++)
        fringe[ii].clear();
}



PrecalculatedMarginal::PrecalculatedMarginal(Marginal&& m,
    double lCutOff,
    bool sort,
    int tabSize,
    int
) : Marginal(std::move(m)),
allocator(isotopeNo, tabSize)
{
    Conf currentConf = allocator.makeCopy(mode_conf);
    if(logProb(currentConf) >= lCutOff)
    {
        configurations.push_back(currentConf);
        lProbs.push_back(mode_lprob);
    }

    unsigned int idx = 0;

    std::unique_ptr<double[]> prob_partials(new double[isotopeNo]);
    std::unique_ptr<double[]> prob_part_acc(new double[isotopeNo+1]);
    prob_part_acc[0] = loggamma_nominator;

    while(idx < configurations.size())
    {
        currentConf = configurations[idx];
        idx++;

        for(size_t ii = 0; ii < isotopeNo; ii++)
            prob_partials[ii] = minuslogFactorial(currentConf[ii]) + currentConf[ii] * atom_lProbs[ii];

        for(unsigned int ii = 0; ii < isotopeNo; ii++ )
        {
            if(currentConf[ii] > mode_conf[ii])
                continue;

            if(currentConf[ii] != 0)
            {
                double prev_partial_ii = prob_partials[ii];
                currentConf[ii]--;
                prob_partials[ii] = minuslogFactorial(currentConf[ii]) + currentConf[ii] * atom_lProbs[ii];

                for(unsigned int jj = 0; jj < isotopeNo; jj++ )
                {
                    prob_part_acc[jj+1] = prob_part_acc[jj] + prob_partials[jj];

                    if(currentConf[jj] < mode_conf[jj])
                        continue;

                    if( ii != jj )
                    {
                        double logp = prob_part_acc[jj] + minuslogFactorial(1+currentConf[jj]) + (1+currentConf[jj]) * atom_lProbs[jj];
                        for(size_t kk = jj+1; kk < isotopeNo; kk++)
                            logp += prob_partials[kk];

                        if (logp >= lCutOff)
                        {
                            auto tmp = allocator.makeCopy(currentConf);
                            tmp[jj]++;
                            configurations.push_back(tmp);
                            lProbs.push_back(logp);
                        }
                    }

                    if (currentConf[jj] > mode_conf[jj])
                        break;
                }
                currentConf[ii]++;
                prob_partials[ii] = prev_partial_ii;
            }

            if(currentConf[ii] < mode_conf[ii])
                break;
        }
    }

    no_confs = configurations.size();
    confs  = configurations.data();

    if(sort && no_confs > 0)
    {
            std::unique_ptr<size_t[]> order_arr(get_inverse_order(lProbs.data(), no_confs));
            impose_order(order_arr.get(), no_confs, lProbs.data(), confs);
    }

    probs.reset(no_confs);
    masses.reset(no_confs);


    for(unsigned int ii = 0; ii < no_confs; ii++)
    {
        probs[ii] = exp(lProbs[ii]);
        masses[ii] = calc_mass(confs[ii], atom_masses, isotopeNo);
    }

    for(size_t ii = 0; ii < ISOSPEC_LPROB_GUARDIANS; ii++)
        lProbs.push_back(-std::numeric_limits<double>::infinity());
}


PrecalculatedMarginal::~PrecalculatedMarginal()
{}







LayeredMarginal::LayeredMarginal(Marginal&& m, int tabSize, int)
: Marginal(std::move(m)), current_threshold(1.0), allocator(isotopeNo, tabSize)
{
    fringe.push_back(mode_conf);
    lProbs.push_back(std::numeric_limits<double>::infinity());
    fringe_unn_lprobs.push_back(unnormalized_logProb(mode_conf));
    for(size_t ii = 0; ii < ISOSPEC_LPROB_GUARDIANS; ii++)
        lProbs.push_back(-std::numeric_limits<double>::infinity());
    guarded_lProbs = lProbs.data()+1;
}

bool LayeredMarginal::extend(double new_threshold, bool do_sort)
{
    new_threshold -= loggamma_nominator;
    if(fringe.empty())
        return false;

    for(size_t ii = 0; ii < ISOSPEC_LPROB_GUARDIANS; ii++)
        lProbs.pop_back();  // Remove the trailing -inf guardians; restored at the end

    pod_vector<Conf> new_fringe;
    pod_vector<double> new_fringe_unn_lprobs;

    while(!fringe.empty())
    {
        Conf currentConf = fringe.back();
        fringe.pop_back();

        double opc = fringe_unn_lprobs.back();

        fringe_unn_lprobs.pop_back();
        if(opc < new_threshold)
        {
            new_fringe.push_back(currentConf);
            new_fringe_unn_lprobs.push_back(opc);
        }

        else
        {
            configurations.push_back(currentConf);
            lProbs.push_back(opc+loggamma_nominator);
            for(unsigned int ii = 0; ii < isotopeNo; ii++ )
            {
                if(currentConf[ii] > mode_conf[ii])
                    continue;

                if(currentConf[ii] > 0)
                {
                    currentConf[ii]--;
                    for(unsigned int jj = 0; jj < isotopeNo; jj++ )
                    {
                        if(currentConf[jj] < mode_conf[jj])
                            continue;

                        if( ii != jj )
                        {
                            Conf nc = allocator.makeCopy(currentConf);
                            nc[jj]++;

                            double lpc = unnormalized_logProb(nc);
                            if(lpc >= new_threshold)
                            {
                                fringe.push_back(nc);
                                fringe_unn_lprobs.push_back(lpc);
                            }
                            else
                            {
                                new_fringe.push_back(nc);
                                new_fringe_unn_lprobs.push_back(lpc);
                            }
                        }

                        if(currentConf[jj] > mode_conf[jj])
                            break;
                    }
                    currentConf[ii]++;
                }

                if(currentConf[ii] < mode_conf[ii])
                    break;
            }
        }
    }

    current_threshold = new_threshold;
    fringe.swap(new_fringe);
    fringe_unn_lprobs.swap(new_fringe_unn_lprobs);

    if(do_sort)
    {
        size_t to_sort_size = configurations.size() - probs.size();
        if(to_sort_size > 0)
        {
            std::unique_ptr<size_t[]> order_arr(get_inverse_order(lProbs.data()+1+probs.size(), to_sort_size));
            double* P = lProbs.data()+1+probs.size();
            Conf* C = configurations.data()+probs.size();
            size_t* O = order_arr.get();
            impose_order(O, to_sort_size, P, C);
        }
    }

    if(probs.capacity() * 2 < configurations.size() + 2)
    {
        // Reserve space for new values
        probs.reserve(configurations.size());
        masses.reserve(configurations.size());
    }  // Otherwise we're growing slowly enough that standard reallocations on push_back work better - we waste some extra memory
       // but don't reallocate on every call

//    printVector(lProbs);
    for(unsigned int ii = probs.size(); ii < configurations.size(); ii++)
    {
        probs.push_back(exp(lProbs[ii+1]));
        masses.push_back(calc_mass(configurations[ii], atom_masses, isotopeNo));
    }

    for(size_t ii = 0; ii < ISOSPEC_LPROB_GUARDIANS; ii++)
        lProbs.push_back(-std::numeric_limits<double>::infinity());  // Restore guardians

    guarded_lProbs = lProbs.data()+1;  // Vector might have reallocated its backing storage

    return true;
}


double LayeredMarginal::get_min_mass() const
{
    double ret = std::numeric_limits<double>::infinity();
    for(pod_vector<double>::const_iterator it = masses.cbegin(); it != masses.cend(); ++it)
        if(*it < ret)
            ret = *it;
    return ret;
}


double LayeredMarginal::get_max_mass() const
{
    double ret = -std::numeric_limits<double>::infinity();
    for(pod_vector<double>::const_iterator it = masses.cbegin(); it != masses.cend(); ++it)
        if(*it > ret)
            ret = *it;
    return ret;
}

/* ===============================================================   */

template<bool add_guards>
SingleAtomMarginal<add_guards>::SingleAtomMarginal(Marginal&& m, int, int)
: Marginal(std::move(m)), current_threshold(1.0), extended_to_idx(0)
{
    original_indexes.resize(isotopeNo);
    for(size_t ii = 0; ii < isotopeNo; ++ii)
        original_indexes[ii] = ii;

    std::sort(original_indexes.begin(), original_indexes.end(), [&](int a, int b) {
        return atom_lProbs[a] > atom_lProbs[b];
    });


    masses.reserve(isotopeNo);
    probs.reserve(isotopeNo);

    if constexpr (add_guards)
    {
        lProbs.reserve(isotopeNo+1+ISOSPEC_LPROB_GUARDIANS);
        lProbs.push_back(std::numeric_limits<double>::infinity());
    }
    else
        lProbs.reserve(isotopeNo);

    for(size_t idx : original_indexes)
    {
        lProbs.push_back(atom_lProbs[idx]);
        probs.push_back(exp(lProbs.back()));
        masses.push_back(atom_masses[idx]);
    }

    if constexpr (add_guards)
    {
        for(size_t ii = 0; ii < ISOSPEC_LPROB_GUARDIANS; ii++)
            lProbs.push_back(-std::numeric_limits<double>::infinity());
        guarded_lProbs = lProbs.data()+1;
    }
    else
        guarded_lProbs = lProbs.data();
}

template class SingleAtomMarginal<true>;
template SingleAtomMarginal<false>::SingleAtomMarginal(IsoSpec::Marginal&&, int, int);

}  // namespace IsoSpec