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 "isoSpec++.h"
#include <cmath>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <cstring>
#include <unordered_map>
#include <queue>
#include <utility>
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <string>
#include <cstdint>
#include <limits>
#include <memory>
#include <cassert>
#include <cctype>
#include <charconv>
#include <string_view>
#include <system_error>
#include "platform.h"
#include "conf.h"
#include "dirtyAllocator.h"
#include "operators.h"
#include "summator.h"
#include "marginalTrek++.h"
#include "misc.h"
#include "element_tables.h"
#include "element_lookup.h"
#include "composition.h"
#include "fasta.h"



namespace IsoSpec
{

Iso::Iso() :
disowned(false),
dimNumber(0),
isotopeNumbers(nullptr),
atomCounts(nullptr),
confSize(0),
allDim(0),
marginals(nullptr)
{}


// The shape every raw-array constructor is handed, checked before anything is
// sized from it: a negative element count, an element with no isotopes, or an
// isotope total past int would otherwise size the arrays below (and allDim)
// from garbage. The values themselves are checked per element, in Marginal.
static int checked_dim_number(int dimNumber, const int* isotopeNumbers)
{
    if(dimNumber < 0)
        throw std::invalid_argument("The number of elements must not be negative");
    std::int64_t total = 0;
    for(int ii = 0; ii < dimNumber; ++ii)
    {
        if(isotopeNumbers[ii] < 1)
            throw std::invalid_argument("Every element must have at least one isotope");
        total += isotopeNumbers[ii];
    }
    if(total > std::numeric_limits<int>::max())
        throw std::length_error("Too many isotopes in total");
    return dimNumber;
}

Iso::Iso(
    int             _dimNumber,
    const int*      _isotopeNumbers,
    const int*      _atomCounts,
    const double* const *  _isotopeMasses,
    const double* const *  _isotopeProbabilities
) :
disowned(false),
dimNumber(checked_dim_number(_dimNumber, _isotopeNumbers)),
isotopeNumbers(nullptr),
atomCounts(nullptr),
confSize(_dimNumber * sizeof(int)),
allDim(0),
marginals(nullptr)
{
    // Owned here until construction succeeds: a constructor that throws runs
    // no destructor, so raw members allocated before the throw would leak.
    std::unique_ptr<int[]> isotope_numbers(array_copy<int>(_isotopeNumbers, dimNumber));
    std::unique_ptr<int[]> atom_counts(array_copy<int>(_atomCounts, dimNumber));
    isotopeNumbers = isotope_numbers.get();
    atomCounts = atom_counts.get();

    for(int ii = 0; ii < dimNumber; ++ii)
        allDim += isotopeNumbers[ii];

    std::unique_ptr<double[]> masses(new double[allDim]);
    std::unique_ptr<double[]> probs(new double[allDim]);
    size_t idx = 0;

    for(int ii = 0; ii < dimNumber; ++ii)
        for(int jj = 0; jj < isotopeNumbers[ii]; ++jj)
        {
            masses[idx] = _isotopeMasses[ii][jj];
            probs[idx]  = _isotopeProbabilities[ii][jj];
            ++idx;
        }

    allDim = 0;  // setupMarginals will recalculate it, assuming it's set to 0

    setupMarginals(masses.get(), probs.get());
    isotope_numbers.release();
    atom_counts.release();
}

Iso::Iso(
    int             _dimNumber,
    const int*      _isotopeNumbers,
    const int*      _atomCounts,
    const double*   _isotopeMasses,
    const double*   _isotopeProbabilities
) :
disowned(false),
dimNumber(checked_dim_number(_dimNumber, _isotopeNumbers)),
isotopeNumbers(nullptr),
atomCounts(nullptr),
confSize(_dimNumber * sizeof(int)),
allDim(0),
marginals(nullptr)
{
    // As above: owned here until construction succeeds.
    std::unique_ptr<int[]> isotope_numbers(array_copy<int>(_isotopeNumbers, dimNumber));
    std::unique_ptr<int[]> atom_counts(array_copy<int>(_atomCounts, dimNumber));
    isotopeNumbers = isotope_numbers.get();
    atomCounts = atom_counts.get();

    setupMarginals(_isotopeMasses, _isotopeProbabilities);
    isotope_numbers.release();
    atom_counts.release();
}

Iso::Iso(Iso&& other) :
disowned(other.disowned),
dimNumber(other.dimNumber),
isotopeNumbers(other.isotopeNumbers),
atomCounts(other.atomCounts),
confSize(other.confSize),
allDim(other.allDim),
marginals(other.marginals)
{
    other.disowned = true;
}


Iso::Iso(const Iso& other, bool fullcopy) :
disowned(!fullcopy),
dimNumber(other.dimNumber),
isotopeNumbers(other.isotopeNumbers),
atomCounts(other.atomCounts),
confSize(other.confSize),
allDim(other.allDim),
marginals(other.marginals)
{
    if(fullcopy)
    {
        // Every copy owned here until all of them exist, then handed over:
        // the release() calls cannot throw.
        std::unique_ptr<int[]> isotope_numbers(array_copy<int>(other.isotopeNumbers, dimNumber));
        std::unique_ptr<int[]> atom_counts(array_copy<int>(other.atomCounts, dimNumber));
        std::unique_ptr<Marginal*[]> marginal_table(new Marginal*[dimNumber]);
        std::vector<std::unique_ptr<Marginal>> copies;
        copies.reserve(dimNumber);
        for(int ii = 0; ii < dimNumber; ii++)
            copies.push_back(std::make_unique<Marginal>(*other.marginals[ii]));

        for(int ii = 0; ii < dimNumber; ii++)
            marginal_table[ii] = copies[ii].release();
        isotopeNumbers = isotope_numbers.release();
        atomCounts = atom_counts.release();
        marginals = marginal_table.release();
    }
}

Iso Iso::FromFASTA(const char* fasta, bool use_nominal_masses, bool add_water)
{
    return FromFASTA(std::string_view(fasta), use_nominal_masses, add_water);
}

Iso Iso::FromFASTA(std::string_view fasta, bool use_nominal_masses, bool add_water)
{
    int atomCounts[6];

    parse_fasta(fasta, atomCounts);

    if(add_water)
    {
        atomCounts[1] += 2;
        atomCounts[3] += 1;
    }

    const int dimNr = atomCounts[5] > 0 ? 6 : 5;

    return Iso(dimNr, aa_isotope_numbers, atomCounts, use_nominal_masses ? aa_elem_nominal_masses : aa_elem_masses, aa_elem_probabilities);
}

inline void Iso::setupMarginals(const double* _isotopeMasses, const double* _isotopeProbabilities)
{
    if (marginals == nullptr)
    {
        int ii = 0;
        marginals = new Marginal*[dimNumber];
        try
        {
            while(ii < dimNumber)
            {
                marginals[ii] = new Marginal(
                        &_isotopeMasses[allDim],
                        &_isotopeProbabilities[allDim],
                        isotopeNumbers[ii],
                        atomCounts[ii]
                    );
                allDim += isotopeNumbers[ii];
                ii++;
            }
        }
        catch(...)
        {
            ii--;
            while(ii >= 0)
            {
                delete marginals[ii];
                ii--;
            }
            delete[] marginals;
            marginals = nullptr;
            throw;
        }
    }
}

Iso::~Iso()
{
    if(!disowned)
    {
    if (marginals != nullptr)
        dealloc_table(marginals, dimNumber);
    delete[] isotopeNumbers;
    delete[] atomCounts;
    }
}

bool Iso::doMarginalsNeedSorting() const
{
    int nontrivial_marginals = 0;
    for(int ii = 0; ii < dimNumber; ii++)
    {
        if(marginals[ii]->get_isotopeNo() > 1)
            nontrivial_marginals++;
        if(nontrivial_marginals > 1)
            return true;
    }
    return false;
}

double Iso::getLightestPeakMass() const
{
    double mass = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        mass += marginals[ii]->getLightestConfMass();
    return mass;
}

double Iso::getLightestPeakLProb() const
{
    double lprob = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        lprob += marginals[ii]->getLightestConfLProb();
    return lprob;
}

void Iso::getLightestPeakSignature(int* space) const
{
    for (int ii = 0; ii < dimNumber; ii++)
    {
        size_t lightest_idx = marginals[ii]->getLightestAtomIndex();
        for(size_t jj = 0; jj < static_cast<size_t>(isotopeNumbers[ii]); jj++)
            space[jj] = (jj == lightest_idx ? atomCounts[ii] : 0);
        space += isotopeNumbers[ii];
    }
}

double Iso::getHeaviestPeakMass() const
{
    double mass = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        mass += marginals[ii]->getHeaviestConfMass();
    return mass;
}

double Iso::getMassRoundingBound() const
{
    // A marginal's configuration mass is calc_mass(), k = isotopeNumbers[ii]
    // products and sums, off by at most ~k ulps of atomCnt * max |atom mass|;
    // the generators then sum dimNumber of those, and the extremes are summed
    // the same way from one rounded product each. Both errors are relative to
    // T = sum over marginals of atomCnt * max |atom mass| (which is what the
    // max below computes, since the extremes are atomCnt times the extreme
    // atom masses), so (max k + 2 dimNumber) ulps of T cover their sum; the +4
    // covers T's own rounding and the caller's subtraction of the bound.
    double total = 0.0;
    int max_isotopes = 0;
    for (int ii = 0; ii < dimNumber; ii++)
    {
        total += (std::max)(std::fabs(marginals[ii]->getLightestConfMass()),
                            std::fabs(marginals[ii]->getHeaviestConfMass()));
        max_isotopes = (std::max)(max_isotopes, isotopeNumbers[ii]);
    }
    return (max_isotopes + 2.0 * dimNumber + 4.0) * std::numeric_limits<double>::epsilon() * total;
}

double Iso::getMaxAtomMassSpan() const
{
    double span = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        span = (std::max)(span, marginals[ii]->getAtomMassSpan());
    return span;
}

double Iso::getHeaviestPeakLProb() const
{
    double lprob = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        lprob += marginals[ii]->getHeaviestConfLProb();
    return lprob;
}

void Iso::getHeaviestPeakSignature(int* space) const
{
    for (int ii = 0; ii < dimNumber; ii++)
    {
        size_t heaviest_idx = marginals[ii]->getHeaviestAtomIndex();
        for(size_t jj = 0; jj < static_cast<size_t>(isotopeNumbers[ii]); jj++)
            space[jj] = (jj == heaviest_idx ? atomCounts[ii] : 0);
        space += isotopeNumbers[ii];
    }
}

double Iso::getMonoisotopicPeakMass() const
{
    double mass = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        mass += marginals[ii]->getMonoisotopicConfMass();
    return mass;
}

double Iso::getMonoisotopicPeakLProb() const
{
    double lprob = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        lprob += marginals[ii]->getMonoisotopicConfLProb();
    return lprob;
}

void Iso::getMonoisotopicPeakSignature(int* space) const
{
    for (int ii = 0; ii < dimNumber; ii++)
    {
        size_t monoisotopic_idx = marginals[ii]->getMonoisotopicAtomIndex();
        for(size_t jj = 0; jj < static_cast<size_t>(isotopeNumbers[ii]); jj++)
            space[jj] = (jj == monoisotopic_idx ? atomCounts[ii] : 0);
        space += isotopeNumbers[ii];
    }
}

double Iso::getUnlikeliestPeakLProb() const
{
    double lprob = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        lprob += marginals[ii]->getSmallestLProb();
    return lprob;
}

double Iso::getModeMass() const
{
    double mass = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        mass += marginals[ii]->getModeMass();
    return mass;
}

double Iso::getModeLProb() const
{
    double ret = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        ret += marginals[ii]->getModeLProb();
    return ret;
}

double Iso::getTheoreticalAverageMass() const
{
    double mass = 0.0;
    for (int ii = 0; ii < dimNumber; ii++)
        mass += marginals[ii]->getTheoreticalAverageMass();
    return mass;
}

double Iso::variance() const
{
    double ret = 0.0;
    for(int ii = 0; ii < dimNumber; ii++)
        ret += marginals[ii]->variance();
    return ret;
}


Iso::Iso(const char* formula, bool use_nominal_masses) : Iso(std::string_view(formula), use_nominal_masses) {}

Iso::Iso(std::string_view formula, bool use_nominal_masses) :
disowned(false),
allDim(0),
marginals(nullptr)
{
    std::vector<double> isotope_masses;
    std::vector<double> isotope_probabilities;

    dimNumber = parse_formula(formula, isotope_masses, isotope_probabilities, &isotopeNumbers, &atomCounts, &confSize, use_nominal_masses);

    try{
        setupMarginals(isotope_masses.data(), isotope_probabilities.data());
    }
    catch(...)
    {
        // parse_formula allocated these; without this catch they leak whenever
        // setupMarginals throws (e.g. an atom count over the lgamma-table limit).
        delete[] isotopeNumbers;
        delete[] atomCounts;
        isotopeNumbers = nullptr;
        atomCounts = nullptr;
        throw;
    }
}


void Iso::addElement(int atomCount, int noIsotopes, const double* isotopeMasses, const double* isotopeProbabilities)
{
    if(noIsotopes > std::numeric_limits<int>::max() - allDim)
        throw std::length_error("Too many isotopes in total");
    // Owned here until the table holds it. An append that throws leaves its
    // array as it was; one that succeeded before it leaves an array one longer
    // than dimNumber says, which is harmless.
    std::unique_ptr<Marginal> m(new Marginal(isotopeMasses, isotopeProbabilities, noIsotopes, atomCount));
    realloc_append<int>(&isotopeNumbers, noIsotopes, dimNumber);
    realloc_append<int>(&atomCounts, atomCount, dimNumber);
    realloc_append<Marginal*>(&marginals, m.get(), dimNumber);
    m.release();
    dimNumber++;
    confSize += sizeof(int);
    allDim += noIsotopes;
}

void Iso::saveMarginalLogSizeEstimates(double* priorities, double target_total_prob) const
{
    /*
     * We shall now use Gaussian approximations of the marginal multinomial distributions to estimate
     * how many configurations we shall need to visit from each marginal. This should be approximately
     * proportional to the volume of the optimal P-ellipsoid of the marginal, which, in turn is defined
     * by the quantile function of the chi-square distribution plus some modifications.
     *
     * We're dropping the constant factor and the (monotonic) exp() transform - these will be used as keys
     * for sorting, so only the ordering is important.
     */

    double K = allDim - dimNumber;

    double log_R2 = log(InverseChiSquareCDF2(K, target_total_prob));

    for(int ii = 0; ii < dimNumber; ii++)
        priorities[ii] = marginals[ii]->getLogSizeEstimate(log_R2);
}

// Shared core of parse_formula: tokenizes `formula` (IsoSpec's native
// concatenated grammar, e.g. "H3C2N1O1", signed counts allowed, e.g.
// "H-2O-1") into (element_table_first_index, signed_count) pairs -- one pair
// per element occurrence, in the order it appears in the string. Same
// grammar, same validation, same element resolution parse_formula has always
// used; also reused as-is by the Unimod composition-delta loader
// (unimod.cpp), so both read a formula string through one code path instead
// of two.
void parse_formula_tokens(const char* formula, std::vector<int>& element_first_indexes, std::vector<int>& counts)
{
    parse_formula_tokens(std::string_view(formula), element_first_indexes, counts);
}

void parse_formula_tokens(std::string_view formula_view, std::vector<int>& element_first_indexes, std::vector<int>& counts)
{
    // Every index below is bounds-checked against slen: the view need not be
    // NUL-terminated, and a NUL inside it fails the character check.
    const char* const formula = formula_view.data();
    const size_t slen = formula_view.size();
    std::vector<std::pair<const char*, size_t> > elements;

    if(slen == 0)
        throw std::invalid_argument("Invalid formula: can't be empty");

    if(!isdigit(static_cast<unsigned char>(formula[slen-1])))
        throw std::invalid_argument("Invalid formula: every element must be followed by a number - write H2O1 and not H2O for water");

    for(size_t ii = 0; ii < slen; ii++)
        if(!isdigit(static_cast<unsigned char>(formula[ii])) && !isalpha(static_cast<unsigned char>(formula[ii])) && formula[ii] != '-')
            throw std::invalid_argument("Invalid formula: contains invalid (non-digit, non-alpha, non-'-') character");

    size_t position = 0;

    while(position < slen)
    {
        size_t elem_end = position;
        while(elem_end < slen && isalpha(static_cast<unsigned char>(formula[elem_end])))
            elem_end++;
        // An element's count may carry a leading '-' (a signed delta, e.g. a
        // Unimod modification's atom-count change) -- consumed here so the
        // digit run below still starts exactly at the first digit.
        size_t digit_end = elem_end;
        if(digit_end < slen && formula[digit_end] == '-')
            digit_end++;
        size_t digits_start = digit_end;
        while(digit_end < slen && isdigit(static_cast<unsigned char>(formula[digit_end])))
            digit_end++;
        if(digit_end == digits_start)
            throw std::invalid_argument("Invalid formula: every element must be followed by a number - write H2O1 and not H2O for water");
        elements.emplace_back(&formula[position], elem_end-position);
        // from_chars reads only [elem_end, digit_end), unlike std::stoi, which
        // would need a terminator; out_of_range as stoi threw it.
        int count;
        const std::from_chars_result parsed = std::from_chars(formula + elem_end, formula + digit_end, count);
        if(parsed.ec == std::errc::result_out_of_range)
            throw std::out_of_range("Invalid formula: atom count out of range");
        counts.push_back(count);
        position = digit_end;
    }

    for (unsigned int i = 0; i < elements.size(); i++)
    {
        int idx = find_element_table_first_index(elements[i].first, elements[i].second);
        if(idx < 0)
            throw std::invalid_argument("Invalid formula");
        element_first_indexes.push_back(idx);
    }
}

ElementComposition parse_formula_composition(std::string_view formula)
{
    // parse_formula_tokens' two output vectors are exactly an
    // ElementComposition's two members.
    ElementComposition composition;
    parse_formula_tokens(formula, composition.element_first_index, composition.count);

    // An element written twice ("C1C1") would otherwise become two independent marginals of the same element: the
    // right distribution, but an extra dimension the generators pay for, and not what a formula means. Rejected
    // here rather than in parse_formula_tokens, which the Unimod composition-delta loader shares. IsoSpecPy's
    // ParseFormula rejects it too.
    const std::vector<int>& first_indexes = composition.element_first_index;
    for(size_t ii = 1; ii < first_indexes.size(); ii++)
        if(std::find(first_indexes.begin(), first_indexes.begin() + ii, first_indexes[ii]) != first_indexes.begin() + ii)
            throw std::invalid_argument(std::string("Invalid formula: element ") + elem_table_symbol[first_indexes[ii]] +
                                        " appears more than once - write C2H6O1 and not C1H3C1H3O1");
    return composition;
}

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)
{
    return parse_formula(std::string_view(formula), isotope_masses, isotope_probabilities, isotopeNumbers, atomCounts, confSize, use_nominal_masses);
}

unsigned int parse_formula(std::string_view formula, std::vector<double>& isotope_masses, std::vector<double>& isotope_probabilities, int** isotopeNumbers, int** atomCounts, unsigned int* confSize, bool use_nominal_masses)
{
    // The formula path joins the peptide-sequence path here and shares
    // everything downstream of it.
    ElementComposition composition = parse_formula_composition(formula);

    // parse_formula_tokens accepts signed counts because the Unimod
    // composition-delta loader needs them ("H-2O-1" is a perfectly good
    // modification delta). A molecule is not a delta, and this is the point
    // where these tokens become one -- see composition.h for why a negative
    // count reaching Iso's constructor is undefined behaviour rather than
    // merely a wrong answer. Before signed counts were introduced this was
    // unreachable: '-' was rejected outright by the character check in the
    // tokenizer. Keep it unreachable.
    reject_negative_counts(composition);

    ExpandedComposition expanded;
    expand_composition_into(composition, expanded, use_nominal_masses);

    const unsigned int dimNumber = static_cast<unsigned int>(expanded.atomCounts.size());

    // Appended, not assigned: these are out-parameters of a public function
    // and have always added to whatever the caller passed in.
    isotope_masses.insert(isotope_masses.end(),
                          expanded.isotope_masses.begin(), expanded.isotope_masses.end());
    isotope_probabilities.insert(isotope_probabilities.end(),
                                 expanded.isotope_probabilities.begin(), expanded.isotope_probabilities.end());

    // Both allocated before either is handed out, so that neither leaks if the
    // other throws.
    std::unique_ptr<int[]> isotope_numbers(array_copy<int>(expanded.isotopeNumbers.data(), dimNumber));
    std::unique_ptr<int[]> atom_counts(array_copy<int>(expanded.atomCounts.data(), dimNumber));
    *isotopeNumbers = isotope_numbers.release();
    *atomCounts = atom_counts.release();
    *confSize = dimNumber * sizeof(int);

    return dimNumber;
}


/*
 * ----------------------------------------------------------------------------------------------------------
 */



IsoGenerator::IsoGenerator(Iso&& iso, bool alloc_partials) :
    Iso(std::move(iso)),
    mode_lprob(getModeLProb()),
    partialLProbs(alloc_partials ? new double[dimNumber+1] : nullptr),
    partialMasses(alloc_partials ? new double[dimNumber+1] : nullptr),
    partialProbs(alloc_partials ? new double[dimNumber+1] : nullptr)
{
    if(alloc_partials)
    {
        partialLProbs[dimNumber] = 0.0;
        partialMasses[dimNumber] = 0.0;
        partialProbs[dimNumber] = 1.0;
    }
}


IsoGenerator::~IsoGenerator()
{
    if(partialLProbs != nullptr)
        delete[] partialLProbs;
    if(partialMasses != nullptr)
        delete[] partialMasses;
    if(partialProbs != nullptr)
        delete[] partialProbs;
}



/*
 * ----------------------------------------------------------------------------------------------------------
 */

static const double minsqrt = -1.3407796239501852e+154;  // == constexpr(-sqrt(std::numeric_limits<double>::max()));

IsoThresholdGenerator::IsoThresholdGenerator(Iso&& iso, double _threshold, bool _absolute, int tabSize, int hashSize, bool reorder_marginals)
: IsoGenerator(std::move(iso)),
Lcutoff(_threshold <= 0.0 ? minsqrt : (_absolute ? log(_threshold) : log(_threshold) + mode_lprob))
{
    // Reject empty Iso early: the maxConfsLPSum allocation below would otherwise
    // request new double[(size_t)(-1)] == new double[SIZE_MAX] and throw a
    // misleading bad_alloc.
    if (dimNumber == 0)
        throw std::invalid_argument("IsoThresholdGenerator requires a non-empty Iso");
    // NaN would silently give an empty result; anything <= 0 means everything.
    if (std::isnan(_threshold))
        throw std::invalid_argument("Threshold must not be NaN");

    counter = new int[dimNumber];
    maxConfsLPSum = new double[dimNumber-1];
    marginalResultsUnsorted = new PrecalculatedMarginal*[dimNumber];

    empty = false;

    const bool marginalsNeedSorting = doMarginalsNeedSorting();

    for(int ii = 0; ii < dimNumber; ii++)
    {
        counter[ii] = 0;
        marginalResultsUnsorted[ii] = new PrecalculatedMarginal(std::move(*(marginals[ii])),
                                                        Lcutoff - mode_lprob + marginals[ii]->fastGetModeLProb(),
                                                        marginalsNeedSorting,
                                                        tabSize,
                                                        hashSize);

        if(!marginalResultsUnsorted[ii]->inRange(0))
            empty = true;
    }

    if(reorder_marginals && dimNumber > 1)
    {
        OrderMarginalsBySizeDecresing<PrecalculatedMarginal> comparator(marginalResultsUnsorted);
        int* tmpMarginalOrder = new int[dimNumber];

        for(int ii = 0; ii < dimNumber; ii++)
            tmpMarginalOrder[ii] = ii;

        std::sort(tmpMarginalOrder, tmpMarginalOrder + dimNumber, comparator);
        marginalResults = new PrecalculatedMarginal*[dimNumber];

        for(int ii = 0; ii < dimNumber; ii++)
            marginalResults[ii] = marginalResultsUnsorted[tmpMarginalOrder[ii]];

        marginalOrder = new int[dimNumber];
        for(int ii = 0; ii < dimNumber; ii++)
            marginalOrder[tmpMarginalOrder[ii]] = ii;

        delete[] tmpMarginalOrder;
    }
    else
    {
        marginalResults = marginalResultsUnsorted;
        marginalOrder = nullptr;
    }

    lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();

    if(dimNumber > 1)
        maxConfsLPSum[0] = marginalResults[0]->fastGetModeLProb();

    for(int ii = 1; ii < dimNumber-1; ii++)
        maxConfsLPSum[ii] = maxConfsLPSum[ii-1] + marginalResults[ii]->fastGetModeLProb();

    lProbs_ptr = lProbs_ptr_start;

    partialLProbs_second = partialLProbs;
    partialLProbs_second++;

    if(!empty)
    {
        recalc(dimNumber-1);
        counter[0]--;
        lProbs_ptr--;
    }
    else
    {
        terminate_search();
        lcfmsv = std::numeric_limits<double>::infinity();
    }
}

void IsoThresholdGenerator::terminate_search()
{
    terminate_search_inline();
}

// Find, scanning downward (towards lProbs_ptr_start) from p, the last position
// whose value is >= lcfmsv. lProbs descends in value as the index rises, so the
// qualifying positions are exactly [start, B] for some boundary B; count_confs
// only scans on lines whose carry was accepted, which guarantees
// *start >= lcfmsv, so B exists. Galloping (probe at distance 1, 2, 4, ...)
// then binary search inside the final bracket: O(log d) in the distance d
// actually moved, instead of a linear scan's O(d).
static inline const double* boundary_scan_down(const double* p, const double* start, double lcfmsv)
{
    if(*p >= lcfmsv)
        return p;
    const double* hi = p;      // invariant: *hi < lcfmsv
    size_t step = 1;
    while(true)
    {
        const double* q = (static_cast<size_t>(hi - start) >= step) ? hi - step : start;
        if(*q >= lcfmsv)
        {
            // Boundary within [q, hi): binary search, kept branchless. Its
            // direction is a data-dependent coin flip, so a branch here
            // mispredicts about half the time; gcc made a cmov of it unasked,
            // clang only with ISOSPEC_UNPREDICTABLE. The window [lo, lo+n)
            // keeps *lo >= lcfmsv, and everything from lo+n on is < lcfmsv.
            const double* lo = q;
            size_t n = static_cast<size_t>(hi - lo);
            while(n > 1)
            {
                const size_t half = n / 2;
                lo = ISOSPEC_UNPREDICTABLE(lo[half] >= lcfmsv) ? lo + half : lo;
                n -= half;
            }
            return lo;
        }
        ISOSPEC_IMPOSSIBLE(q == start);
        hi = q;
        step *= 2;
    }
}

size_t IsoThresholdGenerator::count_confs()
{
    if(empty)
        return 0;

    if(dimNumber == 1)
        return marginalResults[0]->get_no_confs();

    const double* lProbs_ptr_l = marginalResults[0]->get_lProbs_ptr() + marginalResults[0]->get_no_confs();

    std::unique_ptr<const double* []> lProbs_restarts(new const double*[dimNumber]);

    for(int ii = 0; ii < dimNumber; ii++)
        lProbs_restarts[ii] = lProbs_ptr_l;

    size_t count = 0;

    lProbs_ptr_l = boundary_scan_down(lProbs_ptr_l, lProbs_ptr_start, lcfmsv);

    while(true)
    {
        count += lProbs_ptr_l - lProbs_ptr_start + 1;

        int idx = 0;
        int * cntr_ptr = counter;

        while(idx < dimNumber - 1)
        {
            *cntr_ptr = 0;
            idx++;
            cntr_ptr++;
            (*cntr_ptr)++;

            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            if(partialLProbs[idx] + maxConfsLPSum[idx-1] >= Lcutoff)
            {
                short_recalc(idx-1);
                // Level-1 carries continue from the previous line's boundary
                // (which lProbs_ptr_l still holds): within a sweep over
                // marginal 1 the cutoff never drops, so the boundary never
                // rises, and walking it incrementally costs its total movement
                // instead of a re-descent from the sweep top on every line.
                // The monotonicity holds in every reachable case: whenever
                // marginal 1 can produce a multi-line sweep and marginal 0 has
                // more than one configuration, at least two marginals are
                // nontrivial and all are therefore sorted; with a single-
                // configuration marginal 0 the boundary of every accepted line
                // is that lone configuration regardless of sorting. Carries at
                // higher levels can move the boundary back up, and use the
                // saved per-level restart as before.
                lProbs_ptr_l = boundary_scan_down(idx == 1 ? lProbs_ptr_l : lProbs_restarts[idx],
                                                  lProbs_ptr_start, lcfmsv);
                for(idx--; idx > 0; idx--)
                    lProbs_restarts[idx] = lProbs_ptr_l;
                break;
            }
        }
        if(idx == dimNumber - 1)
        {
            reset();
            return count;
        }
    }
}

void IsoThresholdGenerator::reset()
{
    if(empty)
    {
        terminate_search();
        return;
    }

    partialLProbs[dimNumber] = 0.0;

    memset(counter, 0, sizeof(int)*dimNumber);
    recalc(dimNumber-1);
    counter[0]--;

    lProbs_ptr = lProbs_ptr_start - 1;
}

IsoThresholdGenerator::~IsoThresholdGenerator()
{
    delete[] counter;
    delete[] maxConfsLPSum;
    if (marginalResultsUnsorted != marginalResults)
        delete[] marginalResultsUnsorted;
    dealloc_table(marginalResults, dimNumber);
    if(marginalOrder != nullptr)
        delete[] marginalOrder;
}


/*
 * ------------------------------------------------------------------------------------------------------------------------
 */

template <typename MarginalType>
IsoLayeredGeneratorTemplate<MarginalType>::IsoLayeredGeneratorTemplate(Iso&& iso, int tabSize, int hashSize, bool reorder_marginals, double t_prob_hint)
: IsoGenerator(std::move(iso))
{
    // Reject empty Iso early; see comment in IsoThresholdGenerator's ctor.
    if (dimNumber == 0)
        throw std::invalid_argument("IsoLayeredGenerator requires a non-empty Iso");

    counter = new int[dimNumber];
    maxConfsLPSum = new double[dimNumber-1];
    currentLThreshold = nextafter(mode_lprob, -std::numeric_limits<double>::infinity());
    lastLThreshold = (std::numeric_limits<double>::min)();
    marginalResultsUnsorted = new MarginalType*[dimNumber];
    resetPositions = new const double*[dimNumber];
    marginalsNeedSorting = doMarginalsNeedSorting();

    memset(counter, 0, sizeof(int)*dimNumber);

    for(int ii = 0; ii < dimNumber; ii++)
        marginalResultsUnsorted[ii] = new MarginalType(std::move(*(marginals[ii])), tabSize, hashSize);

    if(reorder_marginals && dimNumber > 1)
    {
        double* marginal_priorities = new double[dimNumber];

        saveMarginalLogSizeEstimates(marginal_priorities, t_prob_hint);

        int* tmpMarginalOrder = new int[dimNumber];

        for(int ii = 0; ii < dimNumber; ii++)
            tmpMarginalOrder[ii] = ii;

        // Largest marginal innermost, as in IsoThresholdGenerator: long rows of marginal 0 mean few carries.
        // (Sorting ascending - smallest, or single-isotope, marginal innermost - gave about one configuration per
        // row per layer, so nearly every configuration paid for a carry.) Single-isotope marginals have priority
        // -inf and so end up outermost.
        std::sort(tmpMarginalOrder, tmpMarginalOrder + dimNumber, [&](int a, int b) { return marginal_priorities[a] > marginal_priorities[b]; });
        marginalResults = new MarginalType*[dimNumber];

        for(int ii = 0; ii < dimNumber; ii++)
            marginalResults[ii] = marginalResultsUnsorted[tmpMarginalOrder[ii]];

        marginalOrder = new int[dimNumber];
        for(int ii = 0; ii < dimNumber; ii++)
            marginalOrder[tmpMarginalOrder[ii]] = ii;

        delete[] tmpMarginalOrder;
        delete[] marginal_priorities;
    }
    else
    {
        marginalResults = marginalResultsUnsorted;
        marginalOrder = nullptr;
    }

    lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();

    if(dimNumber > 1)
        maxConfsLPSum[0] = marginalResults[0]->fastGetModeLProb();

    for(int ii = 1; ii < dimNumber-1; ii++)
        maxConfsLPSum[ii] = maxConfsLPSum[ii-1] + marginalResults[ii]->fastGetModeLProb();

    lProbs_ptr = lProbs_ptr_start;

    partialLProbs_second = partialLProbs;
    partialLProbs_second++;

    counter[0]--;
    lProbs_ptr--;
    lastLThreshold = 10.0;
    IsoLayeredGeneratorTemplate<MarginalType>::nextLayer(-0.00001);
}

template <typename MarginalType>
bool IsoLayeredGeneratorTemplate<MarginalType>::nextLayer(double offset)
{
    size_t first_mrg_size = marginalResults[0]->get_no_confs();

    if(lastLThreshold < getUnlikeliestPeakLProb())
        return false;

    lastLThreshold = currentLThreshold;
    currentLThreshold += offset;

    for(int ii = 0; ii < dimNumber; ii++)
    {
        marginalResults[ii]->extend(currentLThreshold - mode_lprob + marginalResults[ii]->fastGetModeLProb(), marginalsNeedSorting);
        counter[ii] = 0;
    }

    lProbs_ptr_start = marginalResults[0]->get_lProbs_ptr();  // vector relocation might have happened

    lProbs_ptr = lProbs_ptr_start + first_mrg_size - 1;

    for(int ii = 0; ii < dimNumber; ii++)
        resetPositions[ii] = lProbs_ptr;

    recalc(dimNumber-1);

    return true;
}


template<typename MarginalType>
void IsoLayeredGeneratorTemplate<MarginalType>::terminate_search()
{
    for(int ii = 0; ii < dimNumber; ii++)
    {
        counter[ii] = marginalResults[ii]->get_no_confs()-1;
        partialLProbs[ii] = -std::numeric_limits<double>::infinity();
    }
    partialLProbs[dimNumber] = -std::numeric_limits<double>::infinity();
    lProbs_ptr = lProbs_ptr_start + marginalResults[0]->get_no_confs()-1;
}

template<typename MarginalType>
IsoLayeredGeneratorTemplate<MarginalType>::~IsoLayeredGeneratorTemplate()
{
    delete[] counter;
    delete[] maxConfsLPSum;
    delete[] resetPositions;
    if (marginalResultsUnsorted != marginalResults)
        delete[] marginalResultsUnsorted;
    dealloc_table(marginalResults, dimNumber);
    if(marginalOrder != nullptr)
      delete[] marginalOrder;
}

template class IsoLayeredGeneratorTemplate<LayeredMarginal>;
// template class IsoLayeredGeneratorTemplate<PrecalculatedMarginal>;
// template class IsoLayeredGeneratorTemplate<MarginalTrek>;
template class IsoLayeredGeneratorTemplate<SingleAtomMarginal<true>>;

/*
 * ------------------------------------------------------------------------------------------------------------------------
 */

template<typename MarginalType>
IsoOrderedGeneratorTemplate<MarginalType>::IsoOrderedGeneratorTemplate(Iso&& iso, int _tabSize, int _hashSize) :
IsoGenerator(std::move(iso), false), allocator(dimNumber, _tabSize)
{
    // Reject empty Iso early; see comment in IsoThresholdGenerator's ctor.
    if (dimNumber == 0)
        throw std::invalid_argument("IsoOrderedGenerator requires a non-empty Iso");

    partialLProbs = &currentLProb;
    partialMasses = &currentMass;
    partialProbs = &currentProb;

    marginalResults = new MarginalType*[dimNumber];

    for(int i = 0; i < dimNumber; i++)
        marginalResults[i] = new MarginalType(std::move(*(marginals[i])), _tabSize, _hashSize);

    logProbs        = new const pod_vector<double>*[dimNumber];
    masses          = new const pod_vector<double>*[dimNumber];

    for(int i = 0; i < dimNumber; i++)
    {
        masses[i] = &marginalResults[i]->conf_masses();
        logProbs[i] = &marginalResults[i]->conf_lprobs();
    }

    topConf = allocator.newConf();
    memset(
            reinterpret_cast<char*>(topConf) + sizeof(double),
            0,
            sizeof(int)*dimNumber
    );

    *(reinterpret_cast<double*>(topConf)) =
    combinedSum(
                getConf(topConf),
                logProbs,
                dimNumber
    );

    pq.push(topConf);
}

template<typename MarginalType>
IsoOrderedGeneratorTemplate<MarginalType>::~IsoOrderedGeneratorTemplate()
{
    dealloc_table<MarginalType*>(marginalResults, dimNumber);
    delete[] logProbs;
    delete[] masses;
    partialLProbs = nullptr;
    partialMasses = nullptr;
    partialProbs = nullptr;
}

template<typename MarginalType>
bool IsoOrderedGeneratorTemplate<MarginalType>::advanceToNextConfiguration()
{
    if(pq.size() < 1)
        return false;


    topConf = pq.top();
    pq.pop();

    int* topConfIsoCounts = getConf(topConf);

    currentLProb = *(reinterpret_cast<double*>(topConf));
    currentMass = combinedSum( topConfIsoCounts, masses, dimNumber );
    currentProb = exp(currentLProb);

    ccount = -1;
    for(int j = 0; j < dimNumber; ++j)
    {
        if(marginalResults[j]->probeConfigurationIdx(topConfIsoCounts[j] + 1))
        {
            if(ccount == -1)
            {
                topConfIsoCounts[j]++;
                *(reinterpret_cast<double*>(topConf)) = combinedSum(topConfIsoCounts, logProbs, dimNumber);
                pq.push(topConf);
                topConfIsoCounts[j]--;
                ccount = j;
            }
            else
            {
                void* acceptedCandidate = allocator.newConf();
                int* acceptedCandidateIsoCounts = getConf(acceptedCandidate);
                memcpy(acceptedCandidateIsoCounts, topConfIsoCounts, confSize);

                acceptedCandidateIsoCounts[j]++;

                *(reinterpret_cast<double*>(acceptedCandidate)) = combinedSum(acceptedCandidateIsoCounts, logProbs, dimNumber);

                pq.push(acceptedCandidate);
            }
        }
        if(topConfIsoCounts[j] > 0)
            break;
    }
    if(ccount >=0)
        topConfIsoCounts[ccount]++;

    return true;
}

template class IsoOrderedGeneratorTemplate<MarginalTrek>;
template class IsoOrderedGeneratorTemplate<SingleAtomMarginal<false>>;


/*
 * ---------------------------------------------------------------------------------------------------
 */

template<typename IsoType>
IsoStochasticGeneratorTemplate<IsoType>::IsoStochasticGeneratorTemplate(Iso&& iso, size_t no_molecules, double _precision, double _beta_bias, std::mt19937& _rng) :
IsoGenerator(std::move(iso)),
ILG(std::move(*this)),
to_sample_left(no_molecules),
precision(_precision),
beta_bias(_beta_bias),
confs_prob(0.0),
chasing_prob(0.0),
// The default argument is random_gen as seen by the constructing thread;
// recognise it and defer to the drawing thread's own instead (see user_rng).
user_rng(&_rng == &random_gen ? nullptr : &_rng)
{
    // Precision is the probability mass sampling stops at: NaN reaches the
    // binomial sampler's float-to-integer conversion (undefined behaviour),
    // above 1 is never reached and the sampler loops forever.
    if(!(precision > 0.0 && precision <= 1.0))
        throw std::invalid_argument("Precision must be in (0, 1]");
    if(!(beta_bias >= 0.0))
        throw std::invalid_argument("beta_bias must be nonnegative");
}

template class IsoStochasticGeneratorTemplate<IsoLayeredGeneratorTemplate<LayeredMarginal>>;
template class IsoStochasticGeneratorTemplate<IsoLayeredGeneratorTemplate<SingleAtomMarginal<true>>>;
template class IsoStochasticGeneratorTemplate<IsoOrderedGeneratorTemplate<MarginalTrek>>;
template class IsoStochasticGeneratorTemplate<IsoOrderedGeneratorTemplate<SingleAtomMarginal<false>>>;
// template class IsoStochasticGeneratorTemplate<IsoThresholdGenerator>;

/*
 * ---------------------------------------------------------------------------------------------------
 */





}  // namespace IsoSpec