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 <unordered_map>
#include <queue>
#include <limits>
#include <string>
#include <string_view>
#include <vector>
#include <algorithm>
#include "platform.h"
#include "iso_simd.h"
#include "isa_kernels.h"
#include "dirtyAllocator.h"
#include "summator.h"
#include "operators.h"
#include "marginalTrek++.h"
#include "unimod.h"



namespace IsoSpec
{

// Any input is safe: a malformed formula throws std::invalid_argument (or
// std::out_of_range for a count past int), never undefined behaviour. Checked by
// fuzzing under ASan+UBSan (tests/fuzz/fuzz_formula); keep it that way.
// The parsers take a std::string_view and read exactly its bytes: no NUL
// terminator is needed, and a NUL inside is an invalid character rather than
// the end of the input. The const char* overloads read up to the first NUL.
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 = false);
unsigned int parse_formula(const char* formula,
                           std::vector<double>& isotope_masses,
                           std::vector<double>& isotope_probabilities,
                           int** isotopeNumbers,
                           int** atomCounts,
                           unsigned int* confSize,
                           bool use_nominal_masses = false);

// Shared tokenizing/element-resolution core of parse_formula -- see its doc
// comment in isoSpec++.cpp. Reused by the Unimod composition-delta loader
// (unimod.cpp) so a modification's atomic composition is read through the
// same formula-string parser as everything else, not a second one.
void parse_formula_tokens(std::string_view formula,
                           std::vector<int>& element_first_indexes,
                           std::vector<int>& counts);
void parse_formula_tokens(const char* formula,
                           std::vector<int>& element_first_indexes,
                           std::vector<int>& counts);


//! The Iso class for the calculation of the isotopic distribution.
/*!
    It contains full description of the molecule for which one would like to calculate the isotopic distribution.
*/
class ISOSPEC_EXPORT_SYMBOL Iso {
 private:
    //! Set up the marginal isotopic envelopes, corresponding to subisotopologues.
    /*!
        \param _isotopeMasses A table of masses of isotopes of the elements in the chemical formula,
                              e.g. {12.0, 13.003355, 1.007825, 2.014102} for C100H202.
        \param _isotopeProbabilities A table of isotope frequencies of the elements in the chemical formula,
                                     e.g. {.989212, .010788, .999885, .000115} for C100H202.
    */
    void setupMarginals(const double* _isotopeMasses,
                        const double* _isotopeProbabilities);
    bool            disowned;       /*!< A variable showing if the Iso class was specialized by its child-class. If so, then the description of the molecules has been transfered there and Iso is a carcass class, dead as a dodo, an ex-class if you will. */

 protected:
    int             dimNumber;      /*!< The number of elements in the chemical formula of the molecule. */
    int*            isotopeNumbers; /*!< A table with numbers of isotopes for each element. */
    int*            atomCounts;     /*!< A table with numbers of isotopes for each element. */
    unsigned int    confSize;       /*!< The number of bytes needed to represent the counts of isotopes present in the extended chemical formula. */
    int             allDim;         /*!< The total number of isotopes of elements present in a chemical formula, e.g. for H2O (water) it is 2+3=5. */
    Marginal**      marginals;      /*!< The table of pointers to the distributions of individual subisotopologues. */

    bool doMarginalsNeedSorting() const;

 public:
    Iso();

    //! General constructror.
    /*!
        \param _dimNumber The number of elements in the formula, e.g. for C100H202 it would be 2, as there are only carbon and hydrogen atoms.
        \param _isotopeNumbers A table with numbers of isotopes for each element, e.g. for C100H202 it would be {2, 2}, because both C and H have two stable isotopes.
        \param _atomCounts Number of atoms of each element in the formula, e.g. for C100H202 corresponds to {100, 202}.
        \param _isotopeMasses A table of tables of masses of isotopes of the elements in the chemical formula, e.g. {{12.0, 13.003355}, {1.007825, 2.014102}} for C100H202.
        \param _isotopeProbabilities A table of tables of isotope frequencies of the elements in the chemical formula, e.g. {{.989212, .010788}, {.999885, .000115}} for C100H202.

        Every input is validated; anything the algorithms cannot represent throws
        std::invalid_argument (or std::length_error for sizes): a negative
        _dimNumber or atom count, an element with no isotopes, a probability
        outside (0, 1] (NaN included), a mass that is not finite or exceeds
        1e100 in magnitude, or an element whose probabilities sum further than
        1e-3 from 1. Probabilities within that tolerance (rounded abundances)
        are divided by their sum; within 1e-12 of 1 they are used as given.
    */
    Iso(
        int             _dimNumber,
        const int*      _isotopeNumbers,
        const int*      _atomCounts,
        const double*   _isotopeMasses,
        const double*   _isotopeProbabilities
    );
    Iso(
        int             _dimNumber,
        const int*      _isotopeNumbers,
        const int*      _atomCounts,
        const double* const *  _isotopeMasses,
        const double* const *  _isotopeProbabilities
    );

    //! Constructor from a chemical formula: exactly the bytes of the view (see parse_formula).
    Iso(std::string_view formula, bool use_nominal_masses = false);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    //! Constructor from a chemical formula as a C string.
    Iso(const char* formula, bool use_nominal_masses = false);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    //! Constructor from C++ std::string chemical formula, all of it.
    inline Iso(const std::string& formula, bool use_nominal_masses = false) : Iso(std::string_view(formula), use_nominal_masses) {}  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    //! Constructor (named) from aminoacid FASTA sequence as C string.
    /*!
        \param fasta An aminoacid FASTA sequence. May be upper/lower/mixed case, may contain selenocystein (U) or xleucine (J).
                     Other characters, including FASTA codes of indeterminate chemical formula (X, *, -, B, ...) are silently ignored.
                     That means "AEDA", "AE-DA", "EAXXDA", "AE DA" will all result in the same chemical formula.
                     Subisotopologues will be in order: CHNOS, possibly with Se added at an end if present.
        \use_nominal_masses Whether to use nucleon number instead of the real mass of each isotope during calculations.
        \add_water Whether the chain should have the terminating -H and -OH groups at the N and C terminus, respectively.
    */
    static Iso FromFASTA(const char* fasta, bool use_nominal_masses = false, bool add_water = true);

    //! As above, reading exactly the bytes of the view: a NUL is one more ignored character.
    static Iso FromFASTA(std::string_view fasta, bool use_nominal_masses = false, bool add_water = true);

    //! Constructor (named) from aminoacid FASTA sequence as C++ std::string. See above for details.
    static inline Iso FromFASTA(const std::string& fasta, bool use_nominal_masses = false, bool add_water = true) { return FromFASTA(std::string_view(fasta), use_nominal_masses, add_water); }

    //! Constructor (named) from a peptide sequence recognizing [UNIMOD:<id>]
    //! modification brackets (see fasta_mods.h's parse_fasta_with_mods for
    //! the exact notation/placement). Despite this constructor's own name
    //! being consistent with FromFASTA above, "peptide_sequence" -- IsoSpecPy's
    //! own honest name for what this parameter actually is -- would be the
    //! better name if backward compatibility with FromFASTA's name were not
    //! a concern; see docs/ai/unimod.md.
    //! \param mods The Unimod composition-delta table to resolve bracket ids
    //!             against -- embedded_unimod_table() (the compile-time
    //!             default) unless the caller has a specific override table
    //!             (e.g. from unimod_table_for_path()).
    static Iso FromFASTAWithMods(const char* sequence, const UnimodTable& mods, bool use_nominal_masses = false, bool add_water = true);
    static Iso FromFASTAWithMods(std::string_view sequence, const UnimodTable& mods, bool use_nominal_masses = false, bool add_water = true);

    //! Convenience overload: unimod_db_path == nullptr or "" uses the
    //! compile-time embedded default table; otherwise resolves and caches
    //! the table at that path (unimod_table_for_path -- a single cached
    //! (path, table) pair, not a map, since a process realistically uses 0
    //! or 1 distinct override paths in its lifetime).
    static Iso FromFASTAWithMods(const char* sequence, bool use_nominal_masses = false, bool add_water = true, const char* unimod_db_path = nullptr);
    static Iso FromFASTAWithMods(std::string_view sequence, bool use_nominal_masses = false, bool add_water = true, const char* unimod_db_path = nullptr);

    //! The move constructor.
    Iso(Iso&& other);

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

    //! The copy constructor.
    /*!
        \param other The other instance of the Iso class.
        \param fullcopy If false, copy only the number of atoms in the formula, the size of the configuration, the total number of isotopes, and the probability of the mode isotopologue.
    */
    Iso(const Iso& other, bool fullcopy);

    //! Destructor.
    virtual ~Iso();

    //! Get the mass of the lightest peak in the isotopic distribution.
    double getLightestPeakMass() const;

    //! Get the log-probability of the lightest peak in the isotopic distribution.
    double getLightestPeakLProb() const;

    /*!
        Write the counts of isotopes in the lightest peak into target memory location. It must be large enough to accomodate it.
        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
    */
    void getLightestPeakSignature(int* space) const;

    //! Get the mass of the heaviest peak in the isotopic distribution.
    double getHeaviestPeakMass() const;

    //! A bound on the rounding error of every configuration mass the generators compute, and of getLightestPeakMass() and getHeaviestPeakMass().
    /*!
        The exact masses all lie between the exact lightest and heaviest ones, so
        every computed configuration mass lies in [getLightestPeakMass() - bound,
        getHeaviestPeakMass() + bound]. These are summed in different orders (the
        generators add the marginals from the last one down, the extremes from
        the first one up), so without the margin a configuration can come out
        lighter than the computed lightest peak.
    */
    double getMassRoundingBound() const;

    //! The largest mass change from moving one atom between two isotopes of an element: the largest Marginal::getAtomMassSpan().
    /*!
        Configurations above a probability cutoff are connected by such moves
        (their log-probability is M-concave), so the masses of a complete set of
        them leave no gap wider than this.
    */
    double getMaxAtomMassSpan() const;

    //! Get the log-probability of the heaviest peak in the isotopic distribution.
    double getHeaviestPeakLProb() const;

    /*!
        Write the counts of isotopes in the heaviest peak into target memory location. It must be large enough to accomodate it.
        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
    */
    void getHeaviestPeakSignature(int* space) const;

    /*!
        Get the mass of the monoisotopic peak in the isotopic distribution. Monoisotopc molecule is defined as
        consisting only of the most frequent isotopes of each element. These are often (but not always) the
        lightest ones, making this often (but again, not always) equal to getLightestPeakMass()
    */
    double getMonoisotopicPeakMass() const;

    /*!
        Get the log-probability of the monoisotopic peak in the isotopic distribution. Monoisotopc molecule is defined as
        consisting only of the most frequent isotopes of each element. These are often (but not always) the
        lightest ones. Making this often (but again, not always) equal to getLightestPeakLProb()
    */
    double getMonoisotopicPeakLProb() const;

    /*!
        Write the counts of isotopes in the monoisotopic peak into target memory location. It must be large enough to accomodate it.
        \param space An array where counts of isotopes shall be written. Must be as big as the overall number of isotopes.
    */
    void getMonoisotopicPeakSignature(int* space) const;

    //! Get the log-probability of the mode-configuration (if there are many modes, they share this value).
    double getModeLProb() const;

    //! Get the logprobability of the least probable subisotopologue.
    double getUnlikeliestPeakLProb() const;

    //! Get the mass of the mode-configuration (if there are many modes, it is undefined which one will be selected).
    double getModeMass() const;

    //! Get the theoretical average mass of the molecule.
    double getTheoreticalAverageMass() const;

    //! Get the theoretical variance of the distribution.
    double variance() const;

    //! Get the standard deviation of the theoretical distribution.
    double stddev() const { return sqrt(variance()); }

    //! Get the number of elements in the chemical formula of the molecule.
    inline int getDimNumber() const { return dimNumber; }

    //! Get the total number of isotopes of elements present in a chemical formula.
    inline int getAllDim() const { return allDim; }

    //! The number of isotopes of each element: getDimNumber() entries.
    inline const int* getIsotopeNumbers() const { return isotopeNumbers; }

    //! The number of atoms of each element: getDimNumber() entries.
    inline const int* getAtomCounts() const { return atomCounts; }

    //! The isotope masses of element `element` (0 <= element < getDimNumber()):
    //! getIsotopeNumbers()[element] entries, as the molecule was built with them.
    inline const double* getElementIsotopeMasses(int element) const { return marginals[element]->get_masses(); }

    //! The isotope log-probabilities of element `element`, renormalised to sum
    //! to 1 (see the constructor): getIsotopeNumbers()[element] entries.
    inline const double* getElementIsotopeLProbs(int element) const { return marginals[element]->get_lProbs(); }

    //! Add an element to the molecule. Note: this method can only be used BEFORE Iso is used to construct an IsoGenerator instance.
    void addElement(int atomCount, int noIsotopes, const double* isotopeMasses, const double* isotopeProbabilities);

    //! Save estimates of logarithms of target sizes of marginals using Gaussian approximation into argument array. Array priorities must have length equal to dimNumber.
    void saveMarginalLogSizeEstimates(double* priorities, double target_total_prob) const;
};


//! The generator of isotopologues.
/*!
    This class provides the common interface for all isotopic generators.
*/
class ISOSPEC_EXPORT_SYMBOL IsoGenerator : public Iso
{
 public:
    const double mode_lprob;

 protected:
    double* partialLProbs;  /*!< The prefix sum of the log-probabilities of the current isotopologue. */
    double* partialMasses;  /*!< The prefix sum of the masses of the current isotopologue. */
    double* partialProbs;   /*!< The prefix product of the probabilities of the current isotopologue. */

 public:
    //! Advance to the next, not yet visited, most probable isotopologue.
    /*!
        \return Return false if it is not possible to advance.
    */
    virtual bool advanceToNextConfiguration() = 0;

    //! Get the log-probability of the current isotopologue.
    /*!
        \return The log-probability of the current isotopologue.
    */
    virtual double lprob() const { return partialLProbs[0]; }

    //! Get the mass of the current isotopologue.
    /*!
        \return The mass of the current isotopologue.
    */
    virtual double mass()  const { return partialMasses[0]; }

    //! Get the probability of the current isotopologue.
    /*!
        \return The probability of the current isotopologue.
    */
    virtual double prob() const { return partialProbs[0]; }

    //! Write the signature of configuration into target memory location. It must be large enough to accomodate it.
    virtual void get_conf_signature(int* space) const = 0;

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

    //! Destructor.
    virtual ~IsoGenerator();
};



//! The generator of isotopologues sorted by their probability of occurrence.
/*!
    The subsequent isotopologues are generated with diminishing probability, starting from the mode.
    This algorithm take O(N*log(N)) to compute the N isotopologues because of using the Priority Queue data structure.
    Obtaining the N isotopologues can be achieved in O(N) if they are not required to be spit out in the descending order.
*/
template<typename MarginalType>
class ISOSPEC_EXPORT_SYMBOL IsoOrderedGeneratorTemplate: public IsoGenerator
{
 private:
    MarginalType**              marginalResults;                    /*!< Table of pointers to marginal distributions of subisotopologues. */
    std::priority_queue<void*, pod_vector<void*>, ConfOrder> pq;   /*!< The priority queue used to generate isotopologues ordered by descending probability. */
    void*                       topConf;                            /*!< Most probable configuration. */
    DirtyAllocator              allocator;                          /*!< Structure used for alocating memory for isotopologues. */
    const pod_vector<double>**  logProbs;                           /*!< Obtained log-probabilities. */
    const pod_vector<double>**  masses;                             /*!< Obtained masses. */
    double                      currentLProb;                       /*!< The log-probability of the current isotopologue. */
    double                      currentMass;                        /*!< The mass of the current isotopologue. */
    double                      currentProb;                        /*!< The probability of the current isotopologue. */
    int                         ccount;

 public:
    IsoOrderedGeneratorTemplate(const IsoOrderedGeneratorTemplate& other) = delete;
    IsoOrderedGeneratorTemplate& operator=(const IsoOrderedGeneratorTemplate& other) = delete;

    bool advanceToNextConfiguration() override final;

    //! Save the counts of isotopes in the space.
    /*!
        \param space An array where counts of isotopes shall be written.
                     Must be as big as the overall number of isotopes.
    */
    inline void get_conf_signature(int* space) const override final
    {
        if constexpr (std::is_same<MarginalType, MarginalTrek>::value)
        {
            int* c = getConf(topConf);

            if (ccount >= 0)
                c[ccount]--;

            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResults[ii]->confs()[c[ii]], isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }

            if (ccount >= 0)
                c[ccount]++;
        }
        else
            throw std::runtime_error("IsoOrderedGeneratorTemplate::get_conf_signature() called on a non-MarginalTrek generator. This is not supported yet.");
    };

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

    //! Destructor.
    virtual ~IsoOrderedGeneratorTemplate();

    inline void get_conf_by_indexes(int* space)
    {
        if constexpr (std::is_same<MarginalType, SingleAtomMarginal<false>>::value)
        {
            if(dimNumber == 0)
                return;

            int* c = getConf(topConf);
            space[0] = (std::max)(c[0]-1, 0);

            for(int ii = 1; ii < dimNumber; ii++)
                space[ii] = c[ii];
        }
    }
};

using IsoOrderedGenerator = IsoOrderedGeneratorTemplate<MarginalTrek>;


//! The generator of isotopologues above a given threshold value.
/*!
    Attention: the calculated configurations are only partially ordeded and the user should not assume they will be ordered.
    This algorithm computes N isotopologues in O(N).
    It is a considerable advantage w.r.t. the IsoOrderedGenerator.
*/
class ISOSPEC_EXPORT_SYMBOL IsoThresholdGenerator: public IsoGenerator
{
 private:
    int*                    counter;            /*!< An array storing the position of an isotopologue in terms of the subisotopologues ordered by decreasing probability. */
    double*                 maxConfsLPSum;
    const double            Lcutoff;            /*!< The logarithm of the lower bound on the calculated probabilities. */
    PrecalculatedMarginal** marginalResults;
    PrecalculatedMarginal** marginalResultsUnsorted;
    int* marginalOrder;

    const double* lProbs_ptr;
    const double* lProbs_ptr_start;
    double* partialLProbs_second;
    double partialLProbs_second_val, lcfmsv;
    bool empty;

 public:
    IsoThresholdGenerator(const IsoThresholdGenerator& other) = delete;
    IsoThresholdGenerator& operator=(const IsoThresholdGenerator& other) = delete;

    inline void get_conf_signature(int* space) const override final
    {
        counter[0] = lProbs_ptr - lProbs_ptr_start;
        if(marginalOrder != nullptr)
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                int jj = marginalOrder[ii];
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
        else
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
    };

    //! The move-constructor.
    /*!
        \param iso An instance of the Iso class.
        \param _threshold The threshold value.
        \param _absolute If true, the _threshold is interpreted as the absolute minimal peak height for the isotopologues.
                         If false, the _threshold is the fraction of the heighest peak's probability.
        \param tabSize The size of the extension of the table with configurations.
        \param hashSize The size of the hash-table used to store subisotopologues and check if they have been already calculated.
    */
    IsoThresholdGenerator(Iso&& iso, double _threshold, bool _absolute = true, int _tabSize = default_tab_size, int _hashSize = 1000, bool reorder_marginals = true);

    ~IsoThresholdGenerator();

    // Perform highly aggressive inling as this function is often called as while(advanceToNextConfiguration()) {}
    // which leads to an extremely tight loop and some compilers miss this (potentially due to the length of the function).
    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        lProbs_ptr++;

        if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
            return true;

        // If we reached this point, a carry is needed

        return carry();
    }

    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration_no_carry()
    {
        lProbs_ptr++;

        if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
            return true;

        // If we reached this point, a carry is needed

        return false;
    }

    // Forced inline, and with no out-of-line call left on any of its paths
    // (see terminate_search_inline()): this is what makes the caller's
    // while(advanceToNextConfiguration()) loop call-free. With a call anywhere
    // in that loop -- even on the exhaustion exit -- clang kept the caller's
    // own accumulators in stack slots across the whole loop, a store-forwarding
    // round trip per configuration: measured 1.7-2.1x slower than gcc, which
    // inlined carry() unasked.
    ISOSPEC_FORCE_INLINE bool carry()
    {
        int idx = 0;
        lProbs_ptr = lProbs_ptr_start;

        int * cntr_ptr = counter;

        while(idx < dimNumber-1)
        {
            // counter[idx] = 0;
            *cntr_ptr = 0;
            idx++;
            cntr_ptr++;
            // counter[idx]++;
            (*cntr_ptr)++;
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            if(partialLProbs[idx] + maxConfsLPSum[idx-1] >= Lcutoff)
            {
                partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
                partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
                recalc(idx-1);
                return true;
            }
        }

        terminate_search_inline();
        return false;
    }

    //! Hand the current marginal-0 run to a batched kernel, which consumes as
    //! much of it as fills whole vectors and returns how many configurations it
    //! emitted; the caller's scalar loop drains the rest.
    /*! The kernel is passed in rather than called directly because which build
        of it runs is decided at run time, from the CPU -- see isa_kernels.h.
        Everything crossing that boundary is a plain pointer or double, so the
        kernel can live in a translation unit compiled for another ISA. */
    ISOSPEC_FORCE_INLINE size_t batch_fill_run(isa_fill_run_fn kernel, double* out_masses, double* out_probs)
    {
        return kernel(&lProbs_ptr, lProbs_ptr_start,
                      marginalResults[0]->get_masses_ptr(),
                      marginalResults[0]->get_probs_ptr(),
                      lcfmsv, partialMasses[1], partialProbs[1],
                      out_masses, out_probs);
    }

    //! As above, but scattering into a dense bin accumulator.
    ISOSPEC_FORCE_INLINE size_t batch_bin_run(isa_bin_run_fn kernel, double* acc, size_t no_bins,
                                              double hwmm, double inv_bin_width)
    {
        return kernel(&lProbs_ptr, lProbs_ptr_start,
                      marginalResults[0]->get_masses_ptr(),
                      marginalResults[0]->get_probs_ptr(),
                      lcfmsv, partialMasses[1], partialProbs[1],
                      acc, hwmm, inv_bin_width, static_cast<double>(no_bins));
    }


    ISOSPEC_FORCE_INLINE double lprob() const override final { return partialLProbs_second_val + (*(lProbs_ptr)); }
    ISOSPEC_FORCE_INLINE double mass()  const override final { return partialMasses[1] + marginalResults[0]->get_mass(lProbs_ptr - lProbs_ptr_start); }
    ISOSPEC_FORCE_INLINE double prob()  const override final { return partialProbs[1] * marginalResults[0]->get_prob(lProbs_ptr - lProbs_ptr_start); }

    //! Block the subsequent search of isotopologues.
    void terminate_search();

    /*! Reset the generator to the beginning of the sequence. Allows it to be reused, eg. to go through the conf space once, calculate
        the amount of space needed to store configurations, then to allocate that memory, and go through it again, this time saving
        configurations (and *is* in fact faster than allocating a std::vector and depending on it to grow as needed. This is cheaper
        than throwing away the generator and making a new one too: marginal distributions don't need to be recalculated. */
    void reset();

    /*! Count the number of configurations in the distribution. This can be used to pre-allocate enough memory to store it (e.g.
     * std::vector's reserve() method - this is faster than depending on the vector's dynamic resizing, even though it means that
     * the configuration space is walked through twice. This method has to be called before the first call to advanceToNextConfiguration
     * and has undefined results (incl. segfaults) otherwise. */
    size_t count_confs();

 private:
    //! Recalculate the current partial log-probabilities, masses, and probabilities.
    ISOSPEC_FORCE_INLINE void recalc(int idx)
    {
        for(; idx > 0; idx--)
        {
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
            partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
        }
        // Through a local: partialLProbs[0] is a double store that may alias the
        // member, so using the member after it forces a reload of what was just
        // stored -- a store-forwarding round trip on the carry path.
        const double second_val = *partialLProbs_second;
        partialLProbs_second_val = second_val;
        partialLProbs[0] = second_val + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = Lcutoff - second_val;
    }

    //! terminate_search(), inlinable, for carry(); terminate_search() itself
    //! stays out of line so that the exported symbol does not go away.
    ISOSPEC_FORCE_INLINE void terminate_search_inline()
    {
        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;
    }

    ISOSPEC_FORCE_INLINE void short_recalc(int idx)
    {
        for(; idx > 0; idx--)
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
        // Through a local: partialLProbs[0] is a double store that may alias the
        // member, so using the member after it forces a reload of what was just
        // stored -- a store-forwarding round trip on the carry path.
        const double second_val = *partialLProbs_second;
        partialLProbs_second_val = second_val;
        partialLProbs[0] = second_val + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = Lcutoff - second_val;
    }
};




template<typename MarginalType>
class ISOSPEC_EXPORT_SYMBOL IsoLayeredGeneratorTemplate : public IsoGenerator
{
 private:
    int*                    counter;            /*!< An array storing the position of an isotopologue in terms of the subisotopologues ordered by decreasing probability. */
    double*                 maxConfsLPSum;
    double currentLThreshold, lastLThreshold;
    MarginalType** marginalResults;
    MarginalType** marginalResultsUnsorted;
    int* marginalOrder;

    const double* lProbs_ptr;
    const double* lProbs_ptr_start;
    const double** resetPositions;
    double* partialLProbs_second;
    double partialLProbs_second_val, lcfmsv, last_lcfmsv;
    bool marginalsNeedSorting;


 public:
    IsoLayeredGeneratorTemplate(const IsoLayeredGeneratorTemplate& other) = delete;
    IsoLayeredGeneratorTemplate& operator=(const IsoLayeredGeneratorTemplate& other) = delete;

    inline void get_conf_signature(int* space) const override final
    {
        counter[0] = lProbs_ptr - lProbs_ptr_start;
        if(marginalOrder != nullptr)
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                int jj = marginalOrder[ii];
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[jj]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
        else
        {
            for(int ii = 0; ii < dimNumber; ii++)
            {
                memcpy(space, marginalResultsUnsorted[ii]->get_conf(counter[ii]), isotopeNumbers[ii]*sizeof(int));
                space += isotopeNumbers[ii];
            }
        }
    };

    inline double get_currentLThreshold() const { return currentLThreshold; }

    IsoLayeredGeneratorTemplate(Iso&& iso, int _tabSize = default_tab_size, int _hashSize = 1000, bool reorder_marginals = true, double t_prob_hint = 0.99);  // NOLINT(runtime/explicit) - constructor deliberately left to be used as a conversion

    ~IsoLayeredGeneratorTemplate();

    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        do
        {
            if(advanceToNextConfigurationWithinLayer())
                return true;
        } while(IsoLayeredGeneratorTemplate<MarginalType>::nextLayer(-2.0));
        return false;
    }

    ISOSPEC_FORCE_INLINE bool advanceToNextConfigurationWithinLayer()
    {
        do{
            lProbs_ptr++;

            if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
                return true;
        }
        while(carry());  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that's not an empty loop body, that's a do{...}while(...) construct
        return false;
    }

    //! advanceToNextConfigurationWithinLayer() without the carry: returns false at the end of the current
    //! marginal-0 run, leaving the generator where carry() continues from.
    ISOSPEC_FORCE_INLINE bool advanceToNextConfigurationWithinLayer_no_carry()
    {
        lProbs_ptr++;

        if(ISOSPEC_LIKELY(*lProbs_ptr >= lcfmsv))
            return true;

        return false;
    }

    //! Move on to the next marginal-0 run within the current layer; false once the layer is exhausted.
    /*! Unlike IsoThresholdGenerator::carry(), leaves the generator one *before* the run's first configuration,
        i.e. in the "last-emitted" convention batch_fill_run() and advanceToNextConfigurationWithinLayer_no_carry()
        continue from -- no configuration needs emitting in between. */
    ISOSPEC_FORCE_INLINE bool carry()
    {
        // If we reached this point, a carry is needed

        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] >= currentLThreshold)
            {
                partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
                partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
                recalc(idx-1);
                // Walk back to the end of the confs emitted in previous layers, starting from the boundary found for the
                // previous line: marginals are sorted by decreasing probability, so that line was at least as probable
                // and its boundary is never before this one. It is 1-2 steps away, where the boundary of the block's
                // first line (resetPositions[idx]) was 8-19 steps away on large molecules - with marginal 0 the largest,
                // lines are long and the band boundary drifts far within a block. Unlike count_confs()'s continuation
                // this needs no extra select: it only changes which slot is read.
                lProbs_ptr = resetPositions[idx-1];

                while(*lProbs_ptr <= last_lcfmsv)
                    lProbs_ptr--;

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

                return true;
            }
        }

        return false;
    }

    //! Hand the rest of the current marginal-0 run to a batched kernel; see IsoThresholdGenerator::batch_fill_run().
    /*! With marginal 0 the largest marginal (see the constructor), runs within a layer are long enough for this to
        pay: 7.8 configurations per run on C520H817N139O147S8 to p=0.99999, 91% of configurations in runs of 5 or
        more, where the ascending order this generator used to have gave 1.36 and batching never fired. The
        kernel's read-ahead stays inside marginal 0's lProbs array: LayeredMarginal keeps
        ISOSPEC_LPROB_GUARDIANS sentinels past its end, as every marginal does. */
    ISOSPEC_FORCE_INLINE size_t batch_fill_run(isa_fill_run_fn kernel, double* out_masses, double* out_probs)
    {
        return kernel(&lProbs_ptr, lProbs_ptr_start,
                      marginalResults[0]->get_masses_ptr(),
                      marginalResults[0]->get_probs_ptr(),
                      lcfmsv, partialMasses[1], partialProbs[1],
                      out_masses, out_probs);
    }

    //! The longest a marginal-0 run can be in the current layer: what batch_fill_run() and the scalar drain after
    //! it can emit between two carries, at most.
    ISOSPEC_FORCE_INLINE size_t max_run_length() const { return marginalResults[0]->get_no_confs(); }

    //! Move the position back by n configurations within the current marginal-0 run.
    /*! For callers that batch-emit a run and then find they only needed part of it: positioned on the last
        configuration they keep, the generator continues exactly as if it had stopped there. */
    ISOSPEC_FORCE_INLINE void step_back_within_run(size_t n) { lProbs_ptr -= n; }

    ISOSPEC_FORCE_INLINE double lprob() const override final { return partialLProbs_second_val + (*(lProbs_ptr)); };
    ISOSPEC_FORCE_INLINE double mass()  const override final { return partialMasses[1] + marginalResults[0]->get_mass(lProbs_ptr - lProbs_ptr_start); };
    ISOSPEC_FORCE_INLINE double prob()  const override final { return partialProbs[1] * marginalResults[0]->get_prob(lProbs_ptr - lProbs_ptr_start); };

    //! Block the subsequent search of isotopologues.
    void terminate_search();


    //! Recalculate the current partial log-probabilities, masses, and probabilities.
    ISOSPEC_FORCE_INLINE void recalc(int idx)
    {
        for(; idx > 0; idx--)
        {
            partialLProbs[idx] = partialLProbs[idx+1] + marginalResults[idx]->get_lProb(counter[idx]);
            partialMasses[idx] = partialMasses[idx+1] + marginalResults[idx]->get_mass(counter[idx]);
            partialProbs[idx] = partialProbs[idx+1] * marginalResults[idx]->get_prob(counter[idx]);
        }
        // Through a local, as in IsoThresholdGenerator::recalc(). Here the reload
        // was worse than a round trip: clang paired the two subtractions below
        // into one subpd and fed it with a 16-byte load straddling the 8-byte
        // store of the member just made, which cannot be store-forwarded and
        // stalls until the store retires -- measured as the single hottest
        // instruction of carry().
        const double second_val = *partialLProbs_second;
        partialLProbs_second_val = second_val;
        partialLProbs[0] = second_val + marginalResults[0]->get_lProb(counter[0]);
        lcfmsv = currentLThreshold - second_val;
        last_lcfmsv = lastLThreshold - second_val;
    }

    bool nextLayer(double offset);

    void get_conf_by_indexes(int* space) const
    {
        if constexpr (std::is_same<MarginalType, SingleAtomMarginal<true>>::value)
        {
            counter[0] = lProbs_ptr - lProbs_ptr_start;
            if(marginalOrder != nullptr)
            {
                for(int ii = 0; ii < dimNumber; ii++)
                {
                    int jj = marginalOrder[ii];
                    space[ii] = marginalResultsUnsorted[ii]->get_original_position(counter[jj]);
                }
            }
            else
            {
                for(int ii = 0; ii < dimNumber; ii++)
                    space[ii] = marginalResultsUnsorted[ii]->get_original_position(counter[ii]);
            }
        }
        else
            throw std::runtime_error("IsoLayeredGeneratorTemplate::get_conf_by_indexes() called on a non-SingleAtomMarginal generator. This is not supported yet.");
    }
};
using IsoLayeredGenerator = IsoLayeredGeneratorTemplate<LayeredMarginal>;

template<typename IsoType>
class IsoStochasticGeneratorTemplate : public IsoGenerator
{
    IsoType ILG;
    size_t to_sample_left;
    const double precision;
    const double beta_bias;
    double confs_prob;
    double chasing_prob;
    size_t current_count;
    // The caller's own generator, or nullptr for the library's default: the
    // thread_local random_gen of whichever thread is *drawing*, looked up at
    // each draw. Holding a reference to the default instead would pin it to
    // the constructing thread -- raced on if the generator is used from
    // another thread while that one keeps drawing, and dangling once it exits.
    std::mt19937* const user_rng;
    ISOSPEC_FORCE_INLINE std::mt19937& rng() const { return user_rng != nullptr ? *user_rng : random_gen; }

 public:
    IsoStochasticGeneratorTemplate(Iso&& iso, size_t no_molecules, double precision = 0.9999, double beta_bias = 5.0, std::mt19937& rdvariate_gen = random_gen);

    ISOSPEC_FORCE_INLINE size_t count() const { return current_count; }

    ISOSPEC_FORCE_INLINE double mass() const override final { return ILG.mass(); }

    ISOSPEC_FORCE_INLINE double prob() const override final { return static_cast<double>(count()); }

    ISOSPEC_FORCE_INLINE double lprob() const override final { return log(prob()); }

    ISOSPEC_FORCE_INLINE void get_conf_signature(int* space) const override final { ILG.get_conf_signature(space); }

    ISOSPEC_FORCE_INLINE bool advanceToNextConfiguration() override final
    {
        /* This function will be used mainly in very small, tight loops, therefore it makes sense to
         * aggressively inline it, despite its seemingly large body.
         */
        // Looked up once per call, not per draw: in a shared library each
        // thread_local access is a __tls_get_addr call. Against the old
        // reference pinned at construction, sampling 10^6 molecules through
        // the C ABI measured +4% resolving at every draw, +3% resolving here
        // (x86-64-v3, gcc, -mtls-dialect=gnu as in the wheels) -- the price of
        // a generator that may change threads.
        std::mt19937& gen = rng();
        while(true)
        {
            double curr_conf_prob_left, current_prob;

            if(to_sample_left <= 0)
                return false;

            if(confs_prob < chasing_prob)
            {
                // Beta was last
                current_count = 1;
                to_sample_left--;
                if(!ILG.advanceToNextConfiguration())
                    return false;
                current_prob = ILG.prob();
                confs_prob += current_prob;
                while(confs_prob <= chasing_prob)
                {
                    if(!ILG.advanceToNextConfiguration())
                        return false;
                    current_prob = ILG.prob();
                    confs_prob += current_prob;
                }
                if(to_sample_left <= 0)
                    return true;
                curr_conf_prob_left = confs_prob - chasing_prob;
            }
            else
            {
                // Binomial was last
                current_count = 0;
                if(!ILG.advanceToNextConfiguration())
                    return false;
                current_prob = ILG.prob();
                confs_prob += current_prob;
                curr_conf_prob_left = current_prob;
            }

            double prob_left_to_1 = precision - chasing_prob;
            double expected_confs = curr_conf_prob_left * to_sample_left / prob_left_to_1;

            if(expected_confs <= beta_bias)
            {
                // Beta mode: we keep making beta jumps until we leave the current configuration
                chasing_prob += rdvariate_beta_1_b(to_sample_left, gen) * prob_left_to_1;
                while(chasing_prob <= confs_prob)
                {
                    current_count++;
                    to_sample_left--;
                    if(to_sample_left == 0)
                        return true;
                    prob_left_to_1 = precision - chasing_prob;
                    chasing_prob += rdvariate_beta_1_b(to_sample_left, gen) * prob_left_to_1;
                }
                if(current_count > 0)
                    return true;
            }
            else
            {
                // Binomial mode: a single binomial step
                size_t rbin = rdvariate_binom(to_sample_left, curr_conf_prob_left/prob_left_to_1, gen);
                current_count += rbin;
                to_sample_left -= rbin;
                chasing_prob = confs_prob;
                if(current_count > 0)
                    return true;
            }
        };
    }

    ISOSPEC_FORCE_INLINE void get_indexes(int* space)
    {
        ILG.get_conf_by_indexes(space);
    }
};

using IsoStochasticGenerator = IsoStochasticGeneratorTemplate<IsoLayeredGenerator>;

}  // namespace IsoSpec