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 "fixedEnvelopes.h"
#include <cmath>
#include <limits>
#include <memory>
#include <cassert>
#include <algorithm>
#include <stdexcept>
#include <tuple>
#include <utility>
#include <vector>
#include "isoMath.h"
#include "aligned_ptr.h"

namespace IsoSpec
{

// Copies of an envelope are always made into aligned arrays of our own, even
// when the source adopted foreign, arbitrarily aligned ones. The size is known
// up front, so these never grow and stay on the small backend.
FixedEnvelope::FixedEnvelope(const FixedEnvelope& other) :
_masses(nullptr),
_probs(nullptr),
_confs(nullptr),
_confs_no(other._confs_no),
allDim(other.allDim),
sorted_by_mass(other.sorted_by_mass),
sorted_by_prob(other.sorted_by_prob),
total_prob(other.total_prob),
current_size(other._confs_no),
allDimSizeofInt(other.allDimSizeofInt)
{
    if(other._masses != nullptr)
    {
        _masses_owner.reset(_confs_no);
        _masses = _masses_owner.get();
        if(_confs_no > 0)
            memcpy(_masses, other._masses, sizeof(double) * _confs_no);
    }

    if(other._probs != nullptr)
    {
        _probs_owner.reset(_confs_no);
        _probs = _probs_owner.get();
        if(_confs_no > 0)
            memcpy(_probs, other._probs, sizeof(double) * _confs_no);
    }

    if(other._confs != nullptr)
    {
        const size_t no_ints = _confs_no * static_cast<size_t>(allDim);
        _confs_owner.reset(no_ints);
        _confs = _confs_owner.get();
        if(no_ints > 0)
            memcpy(_confs, other._confs, sizeof(int) * no_ints);
    }
}

FixedEnvelope::FixedEnvelope(FixedEnvelope&& other) :
_masses(other._masses),
_probs(other._probs),
_confs(other._confs),
_confs_no(other._confs_no),
allDim(other.allDim),
sorted_by_mass(other.sorted_by_mass),
sorted_by_prob(other.sorted_by_prob),
total_prob(other.total_prob),
current_size(other.current_size),
allDimSizeofInt(other.allDimSizeofInt),
_masses_owner(std::move(other._masses_owner)),
_probs_owner(std::move(other._probs_owner)),
_confs_owner(std::move(other._confs_owner))
{
other._masses = nullptr;
other._probs  = nullptr;
other._confs  = nullptr;
other._confs_no = 0;
other.total_prob = 0.0;
other.current_size = 0;
}

FixedEnvelope& FixedEnvelope::operator=(const FixedEnvelope& other)
{
    if(this == &other)
        return *this;
    // Copy-and-swap: construct the copy first so that if any allocation inside
    // the copy constructor throws, *this is left untouched.  Only after tmp is
    // fully built do we pilfer its pointers; tmp then destructs holding the old
    // ones, which frees them without any possibility of failure.
    FixedEnvelope tmp(other);
    std::swap(_masses,         tmp._masses);
    std::swap(_probs,          tmp._probs);
    std::swap(_confs,          tmp._confs);
    std::swap(_masses_owner,   tmp._masses_owner);
    std::swap(_probs_owner,    tmp._probs_owner);
    std::swap(_confs_owner,    tmp._confs_owner);
    _confs_no       = tmp._confs_no;
    allDim          = tmp.allDim;
    allDimSizeofInt = tmp.allDimSizeofInt;
    sorted_by_mass  = tmp.sorted_by_mass;
    sorted_by_prob  = tmp.sorted_by_prob;
    total_prob      = tmp.total_prob;
    current_size    = tmp.current_size;
    return *this;
}

FixedEnvelope& FixedEnvelope::operator=(FixedEnvelope&& other)
{
    if(this == &other)
        return *this;
    free_array(_masses_owner, _masses);
    free_array(_probs_owner,  _probs);
    free_array(_confs_owner,  _confs);
    _masses_owner   = std::move(other._masses_owner);
    _probs_owner    = std::move(other._probs_owner);
    _confs_owner    = std::move(other._confs_owner);
    _masses         = other._masses;         other._masses      = nullptr;
    _probs          = other._probs;          other._probs       = nullptr;
    _confs          = other._confs;          other._confs       = nullptr;
    _confs_no       = other._confs_no;       other._confs_no    = 0;
    allDim          = other.allDim;
    allDimSizeofInt = other.allDimSizeofInt;
    sorted_by_mass  = other.sorted_by_mass;
    sorted_by_prob  = other.sorted_by_prob;
    total_prob      = other.total_prob;      other.total_prob   = 0.0;
    current_size    = other.current_size;    other.current_size = 0;
    return *this;
}

FixedEnvelope::FixedEnvelope(double* in_masses, double* in_probs, size_t in_confs_no, bool masses_sorted, bool probs_sorted, double _total_prob) :
_masses(in_masses),
_probs(in_probs),
_confs(nullptr),
_confs_no(in_confs_no),
allDim(0),
sorted_by_mass(masses_sorted),
sorted_by_prob(probs_sorted),
total_prob(_total_prob),
current_size(in_confs_no),
allDimSizeofInt(0)
{}

FixedEnvelope::FixedEnvelope(double* in_masses, double* in_probs, int* in_confs, size_t in_confs_no, int _allDim, bool masses_sorted, bool probs_sorted, double _total_prob) :
_masses(in_masses),
_probs(in_probs),
_confs(in_confs),
_confs_no(in_confs_no),
allDim(_allDim),
sorted_by_mass(masses_sorted),
sorted_by_prob(probs_sorted),
total_prob(_total_prob),
current_size(in_confs_no),
allDimSizeofInt(_allDim * sizeof(int))
{}

FixedEnvelope::FixedEnvelope(aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT>&& in_masses,
                             aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT>&& in_probs,
                             size_t in_confs_no, bool masses_sorted, bool probs_sorted, double _total_prob) :
_masses(in_masses.get()),
_probs(in_probs.get()),
_confs(nullptr),
_confs_no(in_confs_no),
allDim(0),
sorted_by_mass(masses_sorted),
sorted_by_prob(probs_sorted),
total_prob(_total_prob),
current_size(in_confs_no),
allDimSizeofInt(0),
_masses_owner(std::move(in_masses)),
_probs_owner(std::move(in_probs))
{}

FixedEnvelope FixedEnvelope::operator+(const FixedEnvelope& other) const
{
    const size_t new_size = _confs_no + other._confs_no;

    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nprobs;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nmasses;
    nprobs.reset(new_size);
    nmasses.reset(new_size);

    // An empty envelope has null arrays, and memcpy() is declared nonnull: a
    // zero-length copy from it is undefined behaviour, so skip it entirely.
    if(_confs_no > 0)
    {
        memcpy(nprobs.get(),  _probs,  sizeof(double) * _confs_no);
        memcpy(nmasses.get(), _masses, sizeof(double) * _confs_no);
    }

    if(other._confs_no > 0)
    {
        memcpy(nprobs.get()+_confs_no,  other._probs,  sizeof(double) * other._confs_no);
        memcpy(nmasses.get()+_confs_no, other._masses, sizeof(double) * other._confs_no);
    }

    return FixedEnvelope(std::move(nmasses), std::move(nprobs), new_size);
}

namespace
{
// Whether every probability of e times factor stays finite: O(1) for factor <= 1,
// a branch-free pass otherwise. Used where a scaled envelope enters arithmetic.
bool scaled_probs_fit(const double* probs, size_t n, double factor)
{
    if(factor <= 1.0)
        return true;
    bool fits = true;
    for(size_t ii = 0; ii < n; ii++)
        fits &= probs[ii] * factor <= (std::numeric_limits<double>::max)();
    return fits;
}

struct Extremes { double max_prob, min_mass, max_mass; };

Extremes extremes(const double* masses, const double* probs, size_t n)
{
    Extremes ret{0.0, masses[0], masses[0]};
    for(size_t ii = 0; ii < n; ii++)
    {
        ret.max_prob = (std::max)(ret.max_prob, probs[ii]);
        ret.min_mass = (std::min)(ret.min_mass, masses[ii]);
        ret.max_mass = (std::max)(ret.max_mass, masses[ii]);
    }
    return ret;
}
}  // namespace

FixedEnvelope FixedEnvelope::operator*(const FixedEnvelope& other) const
{
    if(other._confs_no != 0 && _confs_no > SIZE_MAX / other._confs_no)
        throw std::bad_alloc();

    // Every product and sum in the loop below must stay finite. Rounding is
    // monotonic, so the largest product is that of the two largest
    // probabilities and the extreme sums are those of the extreme masses:
    // checking those, in O(n + m), clears the whole O(n * m) loop exactly.
    if(_confs_no > 0 && other._confs_no > 0)
    {
        const Extremes a = extremes(_masses, _probs, _confs_no);
        const Extremes b = extremes(other._masses, other._probs, other._confs_no);
        if(!(a.max_prob * b.max_prob <= (std::numeric_limits<double>::max)()) ||
           !(std::fabs(a.min_mass + b.min_mass) <= (std::numeric_limits<double>::max)()) ||
           !(std::fabs(a.max_mass + b.max_mass) <= (std::numeric_limits<double>::max)()))
            throw std::overflow_error("Convolution overflows: a probability product or a mass sum is not finite");
    }

    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nprobs;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nmasses;
    nprobs.reset(_confs_no * other._confs_no);
    nmasses.reset(_confs_no * other._confs_no);

    size_t tgt_idx = 0;

    for(size_t ii = 0; ii < _confs_no; ii++)
        for(size_t jj = 0; jj < other._confs_no; jj++)
        {
            nprobs[tgt_idx]  = _probs[ii]  * other._probs[jj];
            nmasses[tgt_idx] = _masses[ii] + other._masses[jj];
            tgt_idx++;
        }

    return FixedEnvelope(std::move(nmasses), std::move(nprobs), tgt_idx);
}

void FixedEnvelope::sort_by_mass()
{
    if(sorted_by_mass)
        return;

    sort_by(_masses);

    sorted_by_mass = true;
    sorted_by_prob = false;
}


void FixedEnvelope::sort_by_prob()
{
    if(sorted_by_prob)
        return;

    sort_by(_probs);

    sorted_by_prob = true;
    sorted_by_mass = false;
}

template<typename T> void reorder_array(T* arr, size_t* order, size_t size, bool can_destroy = false)
{
    std::unique_ptr<size_t[]> order_owned;
    if(!can_destroy)
    {
        order_owned = std::make_unique<size_t[]>(size);
        memcpy(order_owned.get(), order, sizeof(size_t)*size);
        order = order_owned.get();
    }

    for(size_t ii = 0; ii < size; ii++)
        while(order[ii] != ii)
        {
            std::swap(arr[ii], arr[order[ii]]);
            std::swap(order[order[ii]], order[ii]);
        }
}

void FixedEnvelope::sort_by(double* order)
{
    if(_confs_no <= 1)
        return;

    size_t* indices = new size_t[_confs_no];

    for(size_t ii = 0; ii < _confs_no; ii++)
        indices[ii] = ii;

    std::sort<size_t*>(indices, indices + _confs_no, TableOrder<double>(order));

    size_t* inverse = new size_t[_confs_no];

    for(size_t ii = 0; ii < _confs_no; ii++)
        inverse[indices[ii]] = ii;

    delete[] indices;

    reorder_array(_masses, inverse, _confs_no);
    reorder_array(_probs,  inverse, _confs_no, _confs == nullptr);
    if(_confs != nullptr)
    {
        int* swapspace = new int[allDim];
        for(size_t ii = 0; ii < _confs_no; ii++)
            while(inverse[ii] != ii)
            {
                memcpy(swapspace, &_confs[ii*allDim], allDimSizeofInt);
                memcpy(&_confs[ii*allDim], &_confs[inverse[ii]*allDim], allDimSizeofInt);
                memcpy(&_confs[inverse[ii]*allDim], swapspace, allDimSizeofInt);
                std::swap(inverse[inverse[ii]], inverse[ii]);
            }
        delete[] swapspace;
    }
    delete[] inverse;
}


double FixedEnvelope::get_total_prob()
{
    if(std::isnan(total_prob))
    {
        total_prob = 0.0;
        for(size_t ii = 0; ii < _confs_no; ii++)
            total_prob += _probs[ii];
    }
    return total_prob;
}

void FixedEnvelope::scale(double factor)
{
    // A NaN or negative probability later reaches resample()'s float-to-integer
    // conversions (undefined behaviour) or walks it off the end of the array.
    if(!(factor >= 0.0 && factor <= (std::numeric_limits<double>::max)()))
        throw std::invalid_argument("Scaling factor must be finite and nonnegative");
    // Above 1, a large enough probability overflows to +inf. Checked before
    // anything is written, so a rejected call leaves the envelope unchanged.
    // No probability exceeds the total, so a known total that survives the
    // scaling clears them all in O(1) -- always the case in normalize();
    // otherwise a branch-free pass, which vectorises. The total is known only
    // up to rounding (normalize() records exactly 1.0 for probabilities that
    // may sum to 1 + 2^-52; found by fuzz_envelope), hence the factor-2
    // headroom, far beyond any rounding error.
    if(factor > 1.0 && !(total_prob * factor <= 0.5 * (std::numeric_limits<double>::max)()))
    {
        bool fits = true;
        for(size_t ii = 0; ii < _confs_no; ii++)
            fits &= _probs[ii] * factor <= (std::numeric_limits<double>::max)();
        if(!fits)
            throw std::overflow_error("Scaling overflows a probability");
    }
    for(size_t ii = 0; ii < _confs_no; ii++)
        _probs[ii] *= factor;
    total_prob *= factor;
}

void FixedEnvelope::normalize()
{
    if(_confs_no == 0)
        return;
    double tp = get_total_prob();
    if(tp != 1.0)
    {
        // A zero total has nothing to normalise to; an infinite one (finite
        // probabilities summing past DBL_MAX) would scale everything to zero.
        if(!(tp > 0.0 && tp <= (std::numeric_limits<double>::max)()))
            throw std::invalid_argument("Cannot normalize: total probability is zero or not finite");
        scale(1.0/tp);
        total_prob = 1.0;
    }
}

void FixedEnvelope::shift_mass(double value)
{
    // A NaN mass never satisfies bin()'s <= and hangs it; +inf then -inf makes NaN.
    if(!std::isfinite(value))
        throw std::invalid_argument("Mass shift must be finite");
    // A finite shift can still carry a mass past DBL_MAX, and a +inf mass in
    // one envelope convolved with a -inf one in another makes NaN. Checked
    // before anything is written, as in scale(): in O(1) from the two ends
    // when sorted by mass, otherwise in a pass over all of them.
    const double dmax = (std::numeric_limits<double>::max)();
    bool fits = true;
    if(sorted_by_mass && _confs_no > 0)
    {
        fits = std::fabs(_masses[0] + value) <= dmax && std::fabs(_masses[_confs_no-1] + value) <= dmax;
    }
    else
    {
        for(size_t ii = 0; ii < _confs_no; ii++)
            fits &= std::fabs(_masses[ii] + value) <= dmax;
    }
    if(!fits)
        throw std::overflow_error("Mass shift overflows a mass");
    for(size_t ii = 0; ii < _confs_no; ii++)
        _masses[ii] += value;
}

void FixedEnvelope::resample(size_t samples, double beta_bias, std::mt19937& rng)
{
    if(_confs_no == 0)
        throw std::logic_error("Resample called on an empty spectrum");
    if(!(beta_bias >= 0.0))
        throw std::invalid_argument("beta_bias must be nonnegative");

    // The probabilities are about to be replaced by molecule counts, so any
    // cached sum is stale; invalidate it and let get_total_prob() rescan.
    total_prob = NAN;

    double pprob = 0.0;
    double cprob = 0.0;
    size_t pidx = -1;  // Overflows - but it doesn't matter.

    // Sentinel: prevents the inner while(pprob < cprob) from walking off the end
    // if floating-point rounding leaves the accumulated probability sum slightly
    // below 1.0. Safe to overwrite in-place: rdvariate_binom guards succ_prob>=1.0
    // so infinity is handled correctly, and the sentinel is always cleaned up before
    // return — either by the loop zeroing the slot when it advances pidx to here,
    // or by the memset below when the loop terminates earlier.
    _probs[_confs_no-1] = (std::numeric_limits<double>::max)();

    while(samples > 0)
    {
        pprob += _probs[++pidx];
        _probs[pidx] = 0.0;
        double covered_part = (pprob - cprob) / (1.0 - cprob);
        while(samples * covered_part < beta_bias && samples > 0)
        {
            cprob += rdvariate_beta_1_b(samples, rng) * (1.0 - cprob);
            while(pprob < cprob)
            {
                pprob += _probs[++pidx];
                _probs[pidx] = 0.0;
            }
            _probs[pidx] += 1.0;
            samples--;
            covered_part = (pprob - cprob) / (1.0 - cprob);
        }
        if(samples <= 0)
                break;
        size_t nrtaken = rdvariate_binom(samples, covered_part, rng);
        _probs[pidx] += static_cast<double>(nrtaken);
        samples -= nrtaken;
        cprob = pprob;
    }

    pidx++;
    memset(_probs + pidx, 0, sizeof(double)*(_confs_no - pidx));
}

FixedEnvelope FixedEnvelope::LinearCombination(const std::vector<const FixedEnvelope*>& spectra, const std::vector<double>& intensities)
{
    if(spectra.size() != intensities.size())
        throw std::invalid_argument("Spectra and intensities must have the same length");
    return LinearCombination(spectra.data(), intensities.data(), spectra.size());
}

FixedEnvelope FixedEnvelope::LinearCombination(const FixedEnvelope* const * spectra, const double* intensities, size_t size)
{
    // As in scale().
    for(size_t ii = 0; ii < size; ii++)
        if(!(intensities[ii] >= 0.0 && intensities[ii] <= (std::numeric_limits<double>::max)()))
            throw std::invalid_argument("Intensities must be finite and nonnegative");
    size_t ret_size = 0;
    for(size_t ii = 0; ii < size; ii++)
        ret_size += spectra[ii]->_confs_no;

    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> newprobs;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> newmasses;
    newprobs.reset(ret_size);
    newmasses.reset(ret_size);

    size_t cntr = 0;
    bool fits = true;
    for(size_t ii = 0; ii < size; ii++)
    {
        double mul = intensities[ii];
        for(size_t jj = 0; jj < spectra[ii]->_confs_no; jj++)
        {
            newprobs[jj+cntr] = spectra[ii]->_probs[jj] * mul;
            fits &= newprobs[jj+cntr] <= (std::numeric_limits<double>::max)();
        }
        if(spectra[ii]->_confs_no > 0)
            memcpy(newmasses.get() + cntr, spectra[ii]->_masses, sizeof(double) * spectra[ii]->_confs_no);
        cntr += spectra[ii]->_confs_no;
    }
    if(!fits)
        throw std::overflow_error("Linear combination overflows a probability");
    return FixedEnvelope(std::move(newmasses), std::move(newprobs), cntr);
}

double FixedEnvelope::WassersteinDistance(FixedEnvelope& other)
{
    double ret = 0.0;
    if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
        throw std::logic_error("Spectra must be normalized before computing Wasserstein Distance");

    if(_confs_no == 0 || other._confs_no == 0)
        return 0.0;

    sort_by_mass();
    other.sort_by_mass();

    // The merge reads both envelopes through locals rather than members:
    // gcc reloaded other's array pointers from the object on every step of
    // it, where clang hoisted them out of the loop.
    const double* const masses1 = _masses;
    const double* const probs1 = _probs;
    const size_t n1 = _confs_no;
    const double* const masses2 = other._masses;
    const double* const probs2 = other._probs;
    const size_t n2 = other._confs_no;

    size_t idx_this = 0;
    size_t idx_other = 0;

    double acc_prob = 0.0;
    double last_point = 0.0;


    while(idx_this < n1 && idx_other < n2)
    {
        if(masses1[idx_this] < masses2[idx_other])
        {
            ret += (masses1[idx_this] - last_point) * std::abs(acc_prob);
            acc_prob += probs1[idx_this];
            last_point = masses1[idx_this];
            idx_this++;
        }
        else
        {
            ret += (masses2[idx_other] - last_point) * std::abs(acc_prob);
            acc_prob -= probs2[idx_other];
            last_point = masses2[idx_other];
            idx_other++;
        }
    }

    acc_prob = std::abs(acc_prob);

    while(idx_this < n1)
    {
        ret += (masses1[idx_this] - last_point) * acc_prob;
        acc_prob -= probs1[idx_this];
        last_point = masses1[idx_this];
        idx_this++;
    }

    while(idx_other < n2)
    {
        ret += (masses2[idx_other] - last_point) * acc_prob;
        acc_prob -= probs2[idx_other];
        last_point = masses2[idx_other];
        idx_other++;
    }

    return ret;
}


double FixedEnvelope::OrientedWassersteinDistance(FixedEnvelope& other)
{
    double ret = 0.0;
    if((get_total_prob()*0.999 > other.get_total_prob()) || (other.get_total_prob() > get_total_prob()*1.001))
        throw std::logic_error("Spectra must be normalized before computing Wasserstein Distance");

    if(_confs_no == 0 || other._confs_no == 0)
        return 0.0;

    sort_by_mass();
    other.sort_by_mass();

    // Locals rather than members, as in WassersteinDistance().
    const double* const masses1 = _masses;
    const double* const probs1 = _probs;
    const size_t n1 = _confs_no;
    const double* const masses2 = other._masses;
    const double* const probs2 = other._probs;
    const size_t n2 = other._confs_no;

    size_t idx_this = 0;
    size_t idx_other = 0;

    double acc_prob = 0.0;
    double last_point = 0.0;


    while(idx_this < n1 && idx_other < n2)
    {
        if(masses1[idx_this] < masses2[idx_other])
        {
            ret += (masses1[idx_this] - last_point) * acc_prob;
            acc_prob += probs1[idx_this];
            last_point = masses1[idx_this];
            idx_this++;
        }
        else
        {
            ret += (masses2[idx_other] - last_point) * acc_prob;
            acc_prob -= probs2[idx_other];
            last_point = masses2[idx_other];
            idx_other++;
        }
    }

    // acc_prob is the signed CDF difference (this minus other), so peaks of
    // *this* add to it — in the main loop above and here in the tail alike.
    // (The unsigned version folds the sign into abs() before its tails and so
    // subtracts in both; this one must not.)
    while(idx_this < n1)
    {
        ret += (masses1[idx_this] - last_point) * acc_prob;
        acc_prob += probs1[idx_this];
        last_point = masses1[idx_this];
        idx_this++;
    }

    while(idx_other < n2)
    {
        ret += (masses2[idx_other] - last_point) * acc_prob;
        acc_prob -= probs2[idx_other];
        last_point = masses2[idx_other];
        idx_other++;
    }

    return ret;
}

double FixedEnvelope::AbyssalWassersteinDistance(FixedEnvelope& other, double abyss_depth, double other_scale)
{
    // Infinite depth is meaningful (nothing falls into the abyss); NaN or a
    // negative depth or scale only gives a meaningless distance.
    if(!(abyss_depth >= 0.0) || !(other_scale >= 0.0 && other_scale <= (std::numeric_limits<double>::max)()))
        throw std::invalid_argument("Abyss depth must be nonnegative, and the scale finite and nonnegative");
    // A scaled probability overflowing to inf makes NaN of the flows (inf - inf).
    if(!scaled_probs_fit(other._probs, other._confs_no, other_scale))
        throw std::overflow_error("Scaling the other envelope overflows a probability");
    sort_by_mass();
    other.sort_by_mass();

    // Locals rather than members, as in WassersteinDistance().
    const double* const masses1 = _masses;
    const double* const probs1 = _probs;
    const size_t n1 = _confs_no;
    const double* const masses2 = other._masses;
    const double* const probs2 = other._probs;
    const size_t n2 = other._confs_no;

    std::vector<std::pair<double, double>> carried;

    size_t idx_this = 0;
    size_t idx_other = 0;

    auto finished = [&]() -> bool { return idx_this >= n1 && idx_other >= n2; };
    auto next = [&]() -> std::pair<double, double> {
                            if(idx_this >= n1 || (idx_other < n2 && masses1[idx_this] > masses2[idx_other]))
                            {
                                std::pair<double, double> res = std::pair<double, double>(masses2[idx_other], probs2[idx_other]*other_scale);
                                idx_other++;
                                return res;
                            }
                            else
                            {
                                std::pair<double, double> res = std::pair<double, double>(masses1[idx_this], -probs1[idx_this]);
                                idx_this++;
                                return res;
                            }
                        };
    double accd = 0.0;
    double condemned = 0.0;

    while(!finished())
    {
        auto pair = next();
        double m = pair.first;
        double p = pair.second;
        if(!carried.empty() && carried[0].second * p > 0.0)
        {
            carried.emplace_back(m, p);
            continue;
        }

        while(!carried.empty())
        {
            double cm = carried.back().first;
            double cp = carried.back().second;
            if(m - cm >= abyss_depth)
            {
                for(auto it = carried.cbegin(); it != carried.cend(); it++)
                    condemned += fabs(it->second);
                carried.clear();
                break;
            }
            if((cp+p)*p > 0.0)
            {
                accd += fabs((m-cm)*cp);
                p += cp;
                carried.pop_back();
            }
            else
            {
                // The incoming peak is exhausted against part of the carried
                // one; write the remainder back (cp is a copy of the carried
                // amount, so updating it alone would lose the match and the
                // matched mass would be counted a second time — as transport
                // against a later peak, or as condemned mass at the end).
                accd += fabs((m-cm)*p);
                cp += p;
                if(cp == 0.0)
                    carried.pop_back();
                else
                    carried.back().second = cp;
                p = 0.0;
                break;
            }
        }
        if(p != 0.0)
            carried.emplace_back(m, p);
    }

    for(auto it = carried.cbegin(); it != carried.cend(); it++)
        condemned += fabs(it->second);

    return accd + condemned * abyss_depth * 0.5;
}

std::tuple<double, double, double> FixedEnvelope::WassersteinMatch(FixedEnvelope& other, double flow_distance, double other_scale)
{
    // A NaN or negative flow distance never lets the matching loop advance.
    // Infinite is fine: everything may be matched.
    if(!(flow_distance >= 0.0) || !(other_scale >= 0.0 && other_scale <= (std::numeric_limits<double>::max)()))
        throw std::invalid_argument("Flow distance must be nonnegative, and the scale finite and nonnegative");
    // A scaled probability overflowing to inf makes NaN of the flows (inf - inf).
    if(!scaled_probs_fit(other._probs, other._confs_no, other_scale))
        throw std::overflow_error("Scaling the other envelope overflows a probability");
    if(_confs_no == 0)
        return {0.0, other.get_total_prob() * other_scale, 0.0};

    double unmatched1 = 0.0;
    double unmatched2 = 0.0;
    double massflow = 0.0;

    sort_by_mass();
    other.sort_by_mass();

    // Locals rather than members, as in WassersteinDistance().
    const double* const masses1 = _masses;
    const double* const probs1 = _probs;
    const size_t n1 = _confs_no;
    const double* const masses2 = other._masses;
    const double* const probs2 = other._probs;
    const size_t n2 = other._confs_no;

    // What is left of the current peak on each side. Tracked as remainders, not
    // as amounts used: subtracting the smaller remainder from the larger is
    // never negative in floating point, whereas p - used, with used a rounded
    // running sum, can come out a few ulps below zero (the fuzzer found one).
    size_t idx_this = 0;
    size_t idx_other = 0;
    double rem_this = probs1[0];
    double rem_other = n2 > 0 ? probs2[0]*other_scale : 0.0;

    while(idx_this < n1 && idx_other < n2)
    {
        if(masses1[idx_this] < masses2[idx_other] - flow_distance)
        {
            // Too light to reach anything left in `other`.
            unmatched1 += rem_this;
            if(++idx_this < n1)
                rem_this = probs1[idx_this];
        }
        else if(masses2[idx_other] < masses1[idx_this] - flow_distance)
        {
            unmatched2 += rem_other;
            if(++idx_other < n2)
                rem_other = probs2[idx_other]*other_scale;
        }
        else if(rem_this < rem_other)
        {
            massflow += rem_this;
            rem_other -= rem_this;
            if(++idx_this < n1)
                rem_this = probs1[idx_this];
        }
        else
        {
            massflow += rem_other;
            rem_this -= rem_other;
            if(++idx_other < n2)
                rem_other = probs2[idx_other]*other_scale;
        }
    }

    if(idx_this < n1)
    {
        unmatched1 += rem_this;
        for(idx_this++; idx_this < n1; idx_this++)
            unmatched1 += probs1[idx_this];
    }
    if(idx_other < n2)
    {
        unmatched2 += rem_other;
        for(idx_other++; idx_other < n2; idx_other++)
            unmatched2 += probs2[idx_other]*other_scale;
    }

    return {unmatched1, unmatched2, massflow};
}

FixedEnvelope FixedEnvelope::bin(double bin_width, double middle)
{
    // Zero width has a separate, useful meaning here: coalesce exact masses.
    // Negative/NaN widths can leave the loop below unable to consume a peak.
    if(!std::isfinite(bin_width) || bin_width < 0.0 || !std::isfinite(middle))
        throw std::invalid_argument("Bin width must be finite and nonnegative, and middle finite");
    sort_by_mass();

    FixedEnvelope ret;

    if(_confs_no == 0)
        return ret;

    ret.reallocate_memory<false>(ISOSPEC_INIT_TABLE_SIZE);

    // Bin sums can overflow a probability, and with a tiny width the bin
    // arithmetic can overflow a mass; either breaks what FromArrays checked.
    // Once per bin, on the bin's way out.
    bool fits = true;
    const double dmax = (std::numeric_limits<double>::max)();

    if(bin_width == 0)
    {
        double curr_mass = _masses[0];
        double accd_prob = _probs[0];
        for(size_t ii = 1; ii < _confs_no; ii++)
        {
            if(curr_mass != _masses[ii])
            {
                fits &= accd_prob <= dmax;
                ret.store_conf(curr_mass, accd_prob);
                curr_mass = _masses[ii];
                accd_prob = _probs[ii];
            }
            else
                accd_prob += _probs[ii];
        }
        fits &= accd_prob <= dmax;
        ret.store_conf(curr_mass, accd_prob);
        if(!fits)
            throw std::overflow_error("Binning overflows a probability");
        return ret;
    }

    size_t ii = 0;

    // The same grid, with the middle brought within one width of zero (fmod is
    // exact): a middle far from the masses, say 1e308, otherwise leaves the
    // bin arithmetic below with no precision at the masses' scale.
    middle = std::fmod(middle, bin_width);
    double half_width = 0.5*bin_width;
    double hwmm = half_width-middle;

    while(ii < _confs_no)
    {
        double current_bin_middle = floor((_masses[ii]+hwmm)/bin_width)*bin_width + middle;
        double current_bin_end = current_bin_middle + half_width;
        double bin_prob = 0.0;

        // At least one peak per bin: where masses dwarf the width, rounding can
        // put the computed end below the very peak that defined the bin, and a
        // bin that takes nothing would repeat forever, growing the output.
        do
        {
            bin_prob += _probs[ii];
            ii++;
        }
        while(ii < _confs_no && _masses[ii] <= current_bin_end);  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that is the end of a do{...}while(...)
        fits &= (bin_prob <= dmax) & (std::fabs(current_bin_middle) <= dmax);
        ret.store_conf(current_bin_middle, bin_prob);
    }

    if(!fits)
        throw std::overflow_error("Binning overflows a probability or a bin's mass");
    return ret;
}

FixedEnvelope FixedEnvelope::FromArrays(const double* in_masses, const double* in_probs, size_t confs_no)
{
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nmasses;
    aligned_unique_ptr<double, DOUBLE_SIMD_ALIGNMENT> nprobs;
    nmasses.reset(confs_no);
    nprobs.reset(confs_no);

    // Checked during the copy, branch-free, so that the check rides along
    // with the memory traffic the copy pays anyway.
    bool valid = true;
    for(size_t ii = 0; ii < confs_no; ii++)
    {
        const double m = in_masses[ii];
        const double p = in_probs[ii];
        nmasses[ii] = m;
        nprobs[ii] = p;
        valid &= (std::fabs(m) <= max_abs_mass) & (p >= 0.0) & (p <= (std::numeric_limits<double>::max)());
    }
    if(!valid)
        throw std::invalid_argument("Envelope masses must be finite with |mass| <= 1e100, probabilities finite and nonnegative");

    return FixedEnvelope(std::move(nmasses), std::move(nprobs), confs_no);
}

template<bool tgetConfs> void FixedEnvelope::reallocate_memory(size_t new_size)
{
    // aligned_unique_ptr::realloc() does the overflow-checked size arithmetic
    // and throws std::bad_alloc on failure, so there is nothing to check here.
    grow_array(_masses_owner, _masses, new_size, current_size);
    tmasses = _masses + _confs_no;

    grow_array(_probs_owner, _probs, new_size, current_size);
    tprobs  = _probs  + _confs_no;

    constexpr_if(tgetConfs)
    {
        if(allDim > 0 && new_size > SIZE_MAX / static_cast<size_t>(allDim))
            throw std::bad_alloc();
        const size_t new_ints = new_size * static_cast<size_t>(allDim);
        grow_array(_confs_owner, _confs, new_ints, current_size * static_cast<size_t>(allDim));
        tconfs = _confs + (allDim * _confs_no);
    }
    current_size = new_size;
}

template<bool tgetConfs> void FixedEnvelope::allocate_memory(size_t new_size)
{
    fresh_array(_masses_owner, _masses, new_size);
    tmasses = _masses + _confs_no;

    fresh_array(_probs_owner, _probs, new_size);
    tprobs  = _probs  + _confs_no;

    constexpr_if(tgetConfs)
    {
        if(allDim > 0 && new_size > SIZE_MAX / static_cast<size_t>(allDim))
            throw std::bad_alloc();
        fresh_array(_confs_owner, _confs, new_size * static_cast<size_t>(allDim));
        tconfs = _confs + (allDim * _confs_no);
    }
    current_size = new_size;
}

template<bool tgetConfs> void FixedEnvelope::threshold_init(Iso&& iso, double threshold, bool absolute)
{
    IsoThresholdGenerator generator(std::move(iso), threshold, absolute);

    size_t tab_size = generator.count_confs();
    this->allDim = generator.getAllDim();
    this->allDimSizeofInt = this->allDim * sizeof(int);

    // count_confs() has already fixed the exact output size, so this allocates
    // it once and never grows.
    this->allocate_memory<tgetConfs>(tab_size);

    double* ttmasses = this->_masses;
    double* ttprobs = this->_probs;
    ISOSPEC_MAYBE_UNUSED int* ttconfs;
    constexpr_if(tgetConfs)
        ttconfs = _confs;

    constexpr_if(tgetConfs) {
        while(generator.advanceToNextConfiguration())
        {
            *ttmasses = generator.mass(); ttmasses++;
            *ttprobs = generator.prob(); ttprobs++;
            generator.get_conf_signature(ttconfs); ttconfs += allDim;
        }
    }
    else
    {
        // Batched fill. Each marginal-0 run (one per higher-dimensional carry state)
        // starts at index 0 and descends until it drops below the cutoff; the kernel
        // takes the bulk of a run W at a time and the < W tail is drained scalar.
        // Which build of the kernel runs is decided from the CPU, once -- see
        // isa_kernels.h. When the toolchain has no <experimental/simd> at all the
        // kernel emits nothing and the scalar drain below does the whole job.
        //
        // Convention bridge: the kernel and advanceToNextConfiguration_no_carry()
        // both use "lProbs_ptr points at the last-emitted config" (advance, then
        // emit), matching the generator's initial position of one-before-index-0.
        // carry() instead leaves lProbs_ptr *at* the index-0 config to be emitted.
        // So after each successful carry we emit that index-0 config scalar (it is
        // always above the cutoff when carry() succeeds, exactly as the scalar path
        // relies on), which restores the last-emitted convention and lets the next
        // batch resume cleanly from index 1.
        const isa_fill_run_fn fill_run = isa_kernels().fill_run;
        do {
            const size_t batched = generator.batch_fill_run(fill_run, ttmasses, ttprobs);
            ttmasses += batched;
            ttprobs  += batched;

            while(generator.advanceToNextConfiguration_no_carry())
            {
                *ttmasses = generator.mass(); ttmasses++;
                *ttprobs = generator.prob(); ttprobs++;
            }
            if(!generator.carry())
                break;
            *ttmasses = generator.mass(); ttmasses++;
            *ttprobs = generator.prob(); ttprobs++;
        } while(true);
    }

    // The count_confs pre-pass fixes tab_size to the exact number of above-threshold
    // configurations; the fill above must emit precisely that many. A mismatch means
    // the fill and the pre-pass disagree (e.g. a carry-convention or off-by-one bug),
    // which would leave _confs_no inconsistent with the buffer contents.
    ISOSPEC_IMPOSSIBLE(static_cast<size_t>(ttmasses - this->_masses) != tab_size);

    this->_confs_no = tab_size;
}

template void FixedEnvelope::threshold_init<true>(Iso&& iso, double threshold, bool absolute);
template void FixedEnvelope::threshold_init<false>(Iso&& iso, double threshold, bool absolute);


// Batched like threshold_init()'s fill when no configuration signatures are
// stored. This used not to pay: with the layered generator ordering marginals
// ascending, its runs over marginal 0 held 1.36 configurations on average
// (C520H817N139O147S8 to p=0.99999) and a W=4 batch fired on 0.0% of them.
// With marginal 0 the largest marginal the same case averages 7.8 per run, 91%
// of configurations in runs of 5 or more -- see
// IsoLayeredGenerator::batch_fill_run(). With signatures the fill stays scalar:
// get_conf_signature() needs the generator positioned on each configuration.

// Store the rest of the generator's current marginal-0 run: the bulk W at a time
// by the batched kernel, the < W tail scalar. Leaves the generator past the end
// of the run, where carry() continues from; returns how many it stored.
size_t FixedEnvelope::store_run(IsoLayeredGenerator& generator, isa_fill_run_fn fill_run)
{
    const size_t room_needed = generator.max_run_length();
    if(current_size - _confs_no < room_needed)
        reallocate_memory<false>((std::max)(current_size*2, _confs_no + room_needed));

    size_t stored = generator.batch_fill_run(fill_run, tmasses, tprobs);
    tmasses += stored;
    tprobs  += stored;

    while(generator.advanceToNextConfigurationWithinLayer_no_carry())
    {
        *tmasses = generator.mass(); tmasses++;
        *tprobs  = generator.prob(); tprobs++;
        stored++;
    }

    _confs_no += stored;
    return stored;
}

template<bool tgetConfs> void FixedEnvelope::store_layer(IsoLayeredGenerator& generator)
{
    constexpr_if(tgetConfs)
    {
        while(generator.advanceToNextConfigurationWithinLayer())
            this->template addConfILG<tgetConfs>(generator);
    }
    else
    {
        const isa_fill_run_fn fill_run = isa_kernels().fill_run;
        do
            store_run(generator, fill_run);
        while(generator.carry());  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that's not an empty loop body, that's a do{...}while(...) construct
    }
}

template<bool tgetConfs> bool FixedEnvelope::store_layer_to_prob(IsoLayeredGenerator& generator,
                                                                 double& prob_so_far,
                                                                 double target_total_prob)
{
    constexpr_if(tgetConfs)
    {
        while(generator.advanceToNextConfigurationWithinLayer())
        {
            this->template addConfILG<tgetConfs>(generator);
            prob_so_far += *(tprobs-1);  // The just-stored probability
            if(prob_so_far >= target_total_prob)
                return true;
        }
        return false;
    }
    else
    {
        const isa_fill_run_fn fill_run = isa_kernels().fill_run;
        do
        {
            // An index, not a pointer: store_run() may reallocate the arrays.
            const size_t run_start = _confs_no;
            const size_t stored = store_run(generator, fill_run);
            const double* const run_probs = _probs + run_start;
            // The same additions, in the same order, as the scalar fill makes.
            for(size_t ii = 0; ii < stored; ii++)
            {
                prob_so_far += run_probs[ii];
                if(prob_so_far >= target_total_prob)
                {
                    // Reached partway through the run: drop what was stored past
                    // this configuration, and put the generator back on it, as if
                    // the scalar fill had stopped here. store_run() left it one
                    // past the run's last configuration.
                    const size_t surplus = stored - ii - 1;
                    tmasses -= surplus;
                    tprobs  -= surplus;
                    _confs_no -= surplus;
                    generator.step_back_within_run(surplus + 1);
                    return true;
                }
            }
        }
        while(generator.carry());  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that's not an empty loop body, that's a do{...}while(...) construct
        return false;
    }
}

template<bool tgetConfs> void FixedEnvelope::total_prob_init(Iso&& iso, double target_total_prob, bool optimize)
{
    // NaN passes both range tests below and the layered loop never reaches it.
    if(std::isnan(target_total_prob))
        throw std::invalid_argument("Target total probability must not be NaN");
    if(target_total_prob <= 0.0)
        return;

    if(target_total_prob >= 1.0)
    {
        threshold_init<tgetConfs>(std::move(iso), 0.0, true);
        return;
    }

    IsoLayeredGenerator generator(std::move(iso), default_tab_size, 1000, true, std::min<double>(target_total_prob, 0.9999));

    this->allDim = generator.getAllDim();
    this->allDimSizeofInt = this->allDim*sizeof(int);


    this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);

    size_t last_switch = 0;
    double prob_at_last_switch = 0.0;
    double prob_so_far = 0.0;
    double layer_delta;

    const double sum_above = log1p(-target_total_prob) - 2.3025850929940455;  // log(0.1);

    do
    {  // Store confs until we accumulate more prob than needed - and, if optimizing,
       // store also the rest of the last layer
        if(this->template store_layer_to_prob<tgetConfs>(generator, prob_so_far, target_total_prob))
        {
            if(!optimize)
                return;
            this->template store_layer<tgetConfs>(generator);
            break;
        }

        last_switch = this->_confs_no;
        prob_at_last_switch = prob_so_far;

        layer_delta = sum_above - log1p(-prob_so_far);
        layer_delta = (std::max)((std::min)(layer_delta, -0.1), ISOSPEC_TOTALPROB_LAYER_MAXSTEP);
    } while(generator.nextLayer(layer_delta));

    if(!optimize || prob_so_far <= target_total_prob)
        return;

    // Right. We have extra configurations and we have been asked to produce an optimal p-set, so
    // now we shall trim unneeded configurations, using an algorithm dubbed "quicktrim"
    // - similar to the quickselect algorithm, except that we use the cumulative sum of elements
    // left of pivot to decide whether to go left or right, instead of the positional index.
    // We'll be sorting by the prob array, permuting the other ones in parallel.
    //
    // The partition is Hoare-style, with the left side's sum accumulated inside
    // the partition scans themselves. This pass is memory-bound, so what counts
    // is passes over the data and stores into it: the fused sum keeps it at one
    // pass per level (a separate summing pass gives back most of the win), and
    // Hoare swaps only misplaced *pairs* where the previous Lomuto partition
    // swapped once per every element above the pivot. Measured on the harvested
    // real trim inputs: 1.65x faster than Lomuto on Piledriver, 1.3x on Alder
    // Lake, never slower, including L2-resident segment sizes.

    int* conf_swapspace = nullptr;
    constexpr_if(tgetConfs)
        conf_swapspace = reinterpret_cast<int*>(malloc(this->allDimSizeofInt));

    size_t start = last_switch;
    size_t end = this->_confs_no;
    double sum_to_start = prob_at_last_switch;

    PivotSequence pivots;
    while(start < end)
    {
        size_t len = end - start;
        if(len == 1)
        {
            // A single element cannot be partitioned (both Hoare sides must be
            // nonempty): it is kept iff the prefix sum still needs it.
            if(sum_to_start >= target_total_prob)
                end = start;
            break;
        }
        size_t pivot = pivots.next(len) + start;
        const double pprob = this->_probs[pivot];

        // Descending Hoare partition of [start, end) around the value pprob:
        // afterwards [start, split) all >= pprob and [split, end) all <= pprob.
        // left_sum accumulates the values that end up on the left: everything
        // the i-scan walks past stays there, and each swap moves the j-side
        // value into a left slot.
        std::ptrdiff_t ii = static_cast<std::ptrdiff_t>(start) - 1;
        std::ptrdiff_t jj = static_cast<std::ptrdiff_t>(end);
        double left_sum = 0.0;
        while(true)
        {
            while(true)
            {
                ii++;
                if(!(this->_probs[ii] > pprob))
                    break;
                left_sum += this->_probs[ii];
            }
            do { jj--; } while(this->_probs[jj] < pprob);
            if(ii >= jj)
                break;
            swap<tgetConfs>(ii, jj, conf_swapspace);
            left_sum += this->_probs[ii];
        }
        // The scans cross with at most one element between them (it then equals
        // the pivot). jj == ii is the case where that element sits on the left
        // without the i-scan ever having counted it.
        if(jj == ii)
            left_sum += this->_probs[ii];

        size_t split = static_cast<size_t>(jj) + 1;
        if(split == end)
        {
            // Right side came out empty; then _probs[end-1] == pprob (the j-scan
            // stopped on it immediately), so splitting just before it is a valid
            // partition too, and keeps both sides strictly shrinking.
            split = end - 1;
            left_sum -= this->_probs[end-1];
        }

        // Selection part
        double new_csum = sum_to_start + left_sum;
        if(new_csum < target_total_prob)
        {
            sum_to_start = new_csum;
            start = split;
        }
        else
            end = split;
    }

    constexpr_if(tgetConfs)
        free(conf_swapspace);

    // Set before the shrink below, so that the tmasses/tprobs/tconfs cursors it
    // recomputes land inside the smaller buffer rather than past its end.
    this->_confs_no = end;

    if(end <= current_size/2)
        // Overhead in memory of 2x or more, shrink to fit
        this->template reallocate_memory<tgetConfs>(end);
}

template void FixedEnvelope::total_prob_init<true>(Iso&& iso, double target_total_prob, bool optimize);
template void FixedEnvelope::total_prob_init<false>(Iso&& iso, double target_total_prob, bool optimize);

template<bool tgetConfs> void FixedEnvelope::stochastic_init(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias, std::mt19937& rng)
{
    IsoStochasticGenerator generator(std::move(iso), _no_molecules, _precision, _beta_bias, rng);

    this->allDim = generator.getAllDim();
    this->allDimSizeofInt = this->allDim * sizeof(int);

    this->reallocate_memory<tgetConfs>(ISOSPEC_INIT_TABLE_SIZE);

    while(generator.advanceToNextConfiguration())
        addConfILG<tgetConfs, IsoStochasticGenerator>(generator);
}

template void FixedEnvelope::stochastic_init<true>(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias, std::mt19937& rng);
template void FixedEnvelope::stochastic_init<false>(Iso&& iso, size_t _no_molecules, double _precision, double _beta_bias, std::mt19937& rng);

double FixedEnvelope::empiric_average_mass()
{
    double ret = 0.0;
    for(size_t ii = 0; ii < _confs_no; ii++)
    {
        ret += _masses[ii] * _probs[ii];
    }
    return ret / get_total_prob();
}

double FixedEnvelope::empiric_variance()
{
    // Single sweep replacing the previous implementation, which swept the arrays
    // up to three times (empiric_average_mass(), the centered-square loop, and
    // get_total_prob()). Masses are shifted by a reference m0 (the first mass)
    // before accumulating, so the running sums stay O(spread) rather than
    // O(mass^2); this keeps the final subtraction between same-magnitude
    // quantities and avoids the catastrophic cancellation that a raw
    // Sum(m^2 p) - avg^2 would suffer for the large absolute masses seen here.
    const double m0 = (_confs_no > 0) ? _masses[0] : 0.0;
    double sum_p   = 0.0;  // Sum p
    double sum_dp  = 0.0;  // Sum p (m - m0)
    double sum_d2p = 0.0;  // Sum p (m - m0)^2
    for(size_t ii = 0; ii < _confs_no; ii++)
    {
        const double p  = _probs[ii];
        const double d  = _masses[ii] - m0;
        const double dp = d * p;
        sum_p   += p;
        sum_dp  += dp;
        sum_d2p += d * dp;
    }

    // Match the memoised get_total_prob(): reuse the cached value if present,
    // otherwise the total we just computed (and cache it, as get_total_prob would).
    if(std::isnan(total_prob))
        total_prob = sum_p;

    // avg - m0, then Var = Sum_ii p_ii ((m_ii - m0) - (avg - m0))^2 / total,
    // expanded so the loop above is the only pass over the arrays.
    const double avg_d = (sum_dp + m0 * sum_p) / total_prob - m0;
    return (sum_d2p - 2.0 * avg_d * sum_dp + avg_d * avg_d * sum_p) / total_prob;
}

// Map a mass to its bin index.  Multiplying by the precomputed reciprocal avoids
// a per-configuration division on this hot path.  The slot range in Binned()
// MUST come from this exact same expression (as must bin_slot() and the vector
// one in the bin_run kernel): every step of it -- the rounded add, the rounded
// multiply by a positive number, floor -- is monotonic, so a mass within the
// range Binned() bounds lands on a slot within the range it allocates.
// (Reciprocal vs true division can disagree by <=1 ULP, which at an exact bin
// boundary could move a peak to an adjacent bin, but never out of the
// allocated range.)
static ISOSPEC_FORCE_INLINE double bin_index_fp(double mass, double hwmm, double inv_bin_width)
{
    return floor((mass + hwmm) * inv_bin_width);
}

// The accumulator slot of a mass, for the scatter.  Binned() shifts hwmm so that
// slots, not absolute bin indices, come out of this directly, and makes sure the
// lowest mass maps to bin_index_fp() >= 0, so every (mass + hwmm) * inv_bin_width
// is >= 0, by monotonicity: there truncation, which the conversion does anyway,
// equals floor.  Without the floor() the scalar fills make no call per
// configuration -- below SSE4.1, the wheel's baseline build of this file, it is
// one to libm -- which had them keep the running sums in stack slots across it.
static ISOSPEC_FORCE_INLINE std::ptrdiff_t bin_slot(size_t no_bins, double mass, double hwmm, double inv_bin_width)
{
    const std::ptrdiff_t slot = static_cast<std::ptrdiff_t>((mass + hwmm) * inv_bin_width);
    ISOSPEC_IMPOSSIBLE(slot < 0 || static_cast<size_t>(slot) >= no_bins);
    return slot;
}

// The seeding contract for both fillers below: nonzero_slot just needs to be the
// slot of *some* populated bin, to seed the outward compaction scan in Binned().
// The scan walks both ways from the seed until it hits a gap wider than any the
// configurations can leave (see Binned()), so *any* nonzero bin works equally
// well (it need not be the smallest, leftmost, or the peak) -- we take whichever
// configuration the generator yields first, which is the cheapest choice and
// lands on the mode.  The one requirement is that the seed be nonzero: an empty
// seed lying beyond that gap from the support would let the scan stop before
// ever reaching it.  Hence zero-probability configurations are skipped.  (The
// seed also fixes the order of the output bins.)

// Scatter a generator's entire output into the dense bin accumulator `acc` (used
// for the target>=1 full-enumeration path).  Returns false if empty.
static bool fill_bins_full(IsoThresholdGenerator& generator,
                           double* acc,
                           size_t no_bins,
                           double hwmm,
                           double inv_bin_width,
                           std::ptrdiff_t& nonzero_slot)
{
    bool non_empty;
    while((non_empty = generator.advanceToNextConfiguration()) && generator.prob() == 0.0)
    {}

    if(!non_empty)
        return false;

    nonzero_slot = bin_slot(no_bins, generator.mass(), hwmm, inv_bin_width);
    acc[nonzero_slot] = generator.prob();

    // Same batching as threshold_init()'s fill, including the carry convention:
    // the seeding loop above left the pointer on the configuration it emitted,
    // and IsoThresholdGenerator::carry() likewise leaves it on the index-0
    // configuration of the next run, so that one is scattered by hand before
    // batching resumes at index 1.
    const isa_bin_run_fn bin_run = isa_kernels().bin_run;
    do
    {
        generator.batch_bin_run(bin_run, acc, no_bins, hwmm, inv_bin_width);
        while(generator.advanceToNextConfiguration_no_carry())
            acc[bin_slot(no_bins, generator.mass(), hwmm, inv_bin_width)] += generator.prob();
        if(!generator.carry())
            break;
        acc[bin_slot(no_bins, generator.mass(), hwmm, inv_bin_width)] += generator.prob();
    } while(true);

    return true;
}

// Scatter a layered generator's output into `acc`, stopping once target_total_prob
// is reached.  Drives the layers with a step adaptively tuned to the remaining
// probability (as total_prob_init does), rather than advanceToNextConfiguration()'s
// fixed -2.0 nat step, to cut marginal over-expansion and layer-transition
// (extend()) overhead.  Returns false if the distribution is empty.
static bool fill_bins_to_prob(IsoLayeredGenerator& generator,
                              double* acc,
                              size_t no_bins,
                              double hwmm,
                              double inv_bin_width,
                              double target_total_prob,
                              std::ptrdiff_t& nonzero_slot)
{
    double prob_so_far = 0.0;
    double layer_delta;
    std::ptrdiff_t seed = -1;  // the first configuration's slot, once there is one
    const double sum_above = log1p(-target_total_prob) - 2.3025850929940455;  // log(0.1)

    // Each marginal-0 run is batch-filled into scratch arrays (as store_layer()
    // does into the envelope; see the note above it), then scattered from there
    // configuration by configuration, exactly as the unbatched loop did -- the
    // bin_run kernel cannot stop partway through a run once the target is met.
    const isa_fill_run_fn fill_run = isa_kernels().fill_run;
    std::vector<double> run_masses, run_probs;

    do
    {
        // Marginal 0 grows with every layer, and so does the longest run.
        if(run_masses.size() < generator.max_run_length())
        {
            run_masses.resize(generator.max_run_length());
            run_probs.resize(generator.max_run_length());
        }

        do
        {
            size_t stored = generator.batch_fill_run(fill_run, run_masses.data(), run_probs.data());
            while(generator.advanceToNextConfigurationWithinLayer_no_carry())
            {
                run_masses[stored] = generator.mass();
                run_probs[stored] = generator.prob();
                stored++;
            }

            for(size_t ii = 0; ii < stored; ii++)
            {
                double prob = run_probs[ii];
                if(prob == 0.0)
                    continue;

                std::ptrdiff_t slot = bin_slot(no_bins, run_masses[ii], hwmm, inv_bin_width);
                acc[slot] += prob;
                if(seed < 0)
                    seed = slot;

                prob_so_far += prob;
                if(prob_so_far >= target_total_prob)
                {
                    nonzero_slot = seed;
                    return true;
                }
            }
        }
        while(generator.carry());  // NOLINT(whitespace/empty_loop_body) - cpplint bug, that's not an empty loop body, that's a do{...}while(...) construct

        layer_delta = sum_above - log1p(-prob_so_far);
        // Cap the widest step; see ISOSPEC_BINNED_LAYER_MAXSTEP in fixedEnvelopes.h.
        layer_delta = (std::max)((std::min)(layer_delta, -0.1), ISOSPEC_BINNED_LAYER_MAXSTEP);
    } while(generator.nextLayer(layer_delta));

    nonzero_slot = seed;
    return seed >= 0;
}

FixedEnvelope FixedEnvelope::Binned(Iso&& iso, double target_total_prob, double bin_width, double bin_middle)
{
    // Unlike bin(), direct binning requires a genuine grid, so zero width is
    // invalid too. Validate even an empty requested distribution consistently.
    if(!std::isfinite(bin_width) || bin_width <= 0.0 || !std::isfinite(bin_middle))
        throw std::invalid_argument("Bin width must be finite and positive, and middle finite");

    if(std::isnan(target_total_prob))
        throw std::invalid_argument("Target total probability must not be NaN");

    FixedEnvelope ret;

    if(target_total_prob <= 0.0)
        return ret;

    // The same grid, with the middle brought within one width of zero (fmod is
    // exact, and the identity for any |middle| < width): a middle far from the
    // masses would otherwise leave mass + hwmm below with no precision at the
    // masses' scale.
    bin_middle = std::fmod(bin_middle, bin_width);
    const double half_width = 0.5*bin_width;
    const double hwmm = half_width-bin_middle;
    const double inv_bin_width = 1.0/bin_width;

    // Every mass a generator can compute lies in [lo, hi]: a little outside
    // the computed extremes, as they are summed in another order (see
    // getMassRoundingBound()).
    const double rounding = iso.getMassRoundingBound();
    const double lo = iso.getLightestPeakMass() - rounding;
    const double hi = iso.getHeaviestPeakMass() + rounding;
    // Read now: iso is moved into the generator below.
    const double max_step = iso.getMaxAtomMassSpan();

    // Bin indices (signed: masses below bin_middle - half_width have negative
    // ones) beyond 2^48 in magnitude would mean bins narrower than 1/16 of an
    // ulp of the masses: nothing meaningful to compute, and no bound on the
    // rounding of the shift below.  This also bounds the allocation.
    const double first_index = bin_index_fp(lo, hwmm, inv_bin_width);
    const double last_index = bin_index_fp(hi, hwmm, inv_bin_width);
    if(!(std::fabs(first_index) <= 0x1p48 && std::fabs(last_index) <= 0x1p48))
        throw std::length_error("Bins too fine for the precision of the masses");

    // Shift the grid's origin to bin idx_min, so that bin_slot() yields
    // accumulator slots directly, with no subtraction per configuration.
    // idx_min * bin_width rounds, which moves the bin edges by up to 1/16 of a
    // bin at the 2^48 limit (an ulp of the mass, and less at smaller indices):
    // a peak that close to an edge may land in the neighbouring bin, as one
    // already may through rounding the reciprocal.  The same rounding leaves
    // the slot of lo in [1, 3] (true value in [2, 3), give or take 1/4), so
    // spare bins below and one above keep the scatter in bounds should the
    // vector and scalar forms of the index ever disagree by one.  Slots are
    // allocated from the shifted map itself, monotonic like the unshifted one.
    const std::ptrdiff_t idx_min = static_cast<std::ptrdiff_t>(first_index) - 2;
    const double hwmm_slot = hwmm - static_cast<double>(idx_min) * bin_width;
    const double first_slot = bin_index_fp(lo, hwmm_slot, inv_bin_width);
    const double last_slot = bin_index_fp(hi, hwmm_slot, inv_bin_width);
    // Bins so wide that idx_min * bin_width overflows (a few times 1e308; the
    // indices are small then) leave an infinite origin.  Otherwise last_slot
    // is at most 2^49 + 3.
    if(!(std::isfinite(hwmm_slot) && last_slot <= 0x1p50))
        throw std::length_error("Bins too wide for the range of doubles");
    if(!(first_slot >= 0.0))
        throw std::logic_error("Binned: the shifted bin grid starts below its first slot");
    const size_t no_bins = static_cast<size_t>(last_slot) + 2;

    double* acc;
# if ISOSPEC_GOT_MMAN
    acc = reinterpret_cast<double*>(mmap(nullptr, sizeof(double)*no_bins, PROT_READ | PROT_WRITE, MAP_ANONYMOUS | MAP_PRIVATE, -1, 0));
    if(acc == MAP_FAILED)
        throw std::bad_alloc();
#else
    // Platforms without mmap use a zero-initialized heap accumulator.
    acc = reinterpret_cast<double*>(calloc(no_bins, sizeof(double)));
    if(acc == nullptr)
        throw std::bad_alloc();
#endif
    // Keep the original mapping alive through generator construction and output
    // growth, both of which can throw. acc is offset below for absolute indexing;
    // the owner must always retain the actual base address for cleanup.
    const auto release_acc = [no_bins](double* base) noexcept {
#if ISOSPEC_GOT_MMAN
        munmap(base, sizeof(double)*no_bins);
#else
        free(base);
#endif
    };
    std::unique_ptr<double, decltype(release_acc)> acc_owner(acc, release_acc);

#if ISOSPEC_GOT_MMAN && defined(MADV_HUGEPAGE)
    // For a large accumulator the scatter (acc[bin_idx] += prob, at essentially
    // random indices) thrashes the DTLB with 4 KiB pages.  Advise transparent
    // huge pages so the touched region is backed by 2 MiB pages.  Gated by size:
    // for a small array a single on-fault 2 MiB page would just waste memory and
    // add zeroing latency with no TLB benefit.  Overridable at build time.
    // Measured (Alder Lake + Opteron 6380, THP active): ~2-9% on fine-binned large
    // molecules, no regression on coarse/small; needs a genuinely THP-capable host
    // (no effect where the kernel won't hand out anon huge pages).
#  ifndef ISOSPEC_BINNED_HUGEPAGE_MIN_BYTES
#    define ISOSPEC_BINNED_HUGEPAGE_MIN_BYTES (size_t(64) << 20)  // 64 MiB
#  endif
    if(sizeof(double)*no_bins >= ISOSPEC_BINNED_HUGEPAGE_MIN_BYTES)
        madvise(acc, sizeof(double)*no_bins, MADV_HUGEPAGE);
#endif

    std::ptrdiff_t nonzero_slot = 0;
    bool non_empty;

    if(target_total_prob >= 1.0)
    {
        // The whole distribution is requested: a threshold generator with a zero
        // cutoff enumerates everything and is cheaper than the layered generator
        // descending layer-by-layer to the least-probable peak.  Mirrors the
        // >= 1.0 fast path in total_prob_init.
        IsoThresholdGenerator generator(std::move(iso), 0.0, true);
        non_empty = fill_bins_full(
                        generator, acc, no_bins, hwmm_slot, inv_bin_width, nonzero_slot);
    }
    else
    {
        // Pass the requested probability as the layer-sizing hint (as
        // total_prob_init does) so the marginal ordering is tuned to the amount
        // of mass actually needed; the default 0.99 hint over- or under-shoots.
        IsoLayeredGenerator generator(std::move(iso), default_tab_size, 1000, true,
                                      std::min<double>(target_total_prob, 0.9999));
        non_empty = fill_bins_to_prob(
                        generator, acc, no_bins, hwmm_slot, inv_bin_width, target_total_prob, nonzero_slot);
    }

    if(non_empty)
    {
        // How wide a gap the populated bins can leave.  Every set of
        // configurations above a probability cutoff is connected by moves of
        // one atom between two isotopes (log-probability is M-concave in the
        // isotope counts, so its superlevel sets have no holes and no islands),
        // and such a move changes the mass by at most max_step: the masses
        // leave no gap wider than that, however improbable the configurations
        // left out.  It is twice that here, for configurations whose computed
        // log-probability rounds across the cutoff, plus both masses' rounding,
        // plus two bins for where the masses fall within theirs.  (Not covered:
        // the partial last layer of fill_bins_to_prob, an arbitrary subset of
        // one band of probabilities; a configuration from it more than this
        // beyond the rest would be missed.  None found so far.)  Artificial
        // elements may have isotopes far apart, and the 10 Da assumed here
        // before lost whole peaks of them.  Capped at the accumulator's size,
        // which the scans cannot pass anyway, before converting.
        const double gap_bins = (2.0 * max_step + 2.0 * rounding) * inv_bin_width + 2.0;
        const size_t max_gap = gap_bins < static_cast<double>(no_bins) ? static_cast<size_t>(gap_bins) : no_bins;

        ret.reallocate_memory<false>(ISOSPEC_INIT_TABLE_SIZE);

        // Each scan stops after max_gap consecutive empty bins, or at the end
        // of the accumulator: its bound moves out with every populated bin.
        // One compare per bin: with a separate count of empty bins gcc kept
        // both loop conditions as flags, twice the instructions per bin, and
        // the 2000 empty bins every spectrum at 0.01 Da ended with (at 10 Da
        // each way) were 11% of a small peptide's whole Binned().
        const std::ptrdiff_t distance = static_cast<std::ptrdiff_t>(max_gap);
        const std::ptrdiff_t bins = static_cast<std::ptrdiff_t>(no_bins);

        std::ptrdiff_t stop = (std::max)(nonzero_slot - distance, std::ptrdiff_t(0));
        for(std::ptrdiff_t ii = nonzero_slot; ii >= stop; ii--)
            if(acc[ii] > 0.0)
            {
                stop = (std::max)(ii - distance, std::ptrdiff_t(0));
                ret.store_conf(static_cast<double>(ii + idx_min)*bin_width + bin_middle, acc[ii]);
            }

        std::ptrdiff_t end = (std::min)(nonzero_slot + distance + 1, bins);
        for(std::ptrdiff_t ii = nonzero_slot + 1; ii < end; ii++)
            if(acc[ii] > 0.0)
            {
                end = (std::min)(ii + distance + 1, bins);
                ret.store_conf(static_cast<double>(ii + idx_min)*bin_width + bin_middle, acc[ii]);
            }
    }

    return ret;
}

}  // namespace IsoSpec