ParticleSet.h

Go to the documentation of this file.

#pragma once

 
#include <concepts>
#include <cstddef>
#include <expected>
#include <functional>
#include <optional>
#include <random>
#include <span>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>
 
#include <Eigen/Cholesky>
 
#include <EmrysFilter/Core/Allocation.h>
#include <EmrysFilter/Core/Execution.h>
#include <EmrysFilter/Core/Forwarding.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Multiply.h>
#include <EmrysFilter/Core/Print.h>
#include <EmrysFilter/Core/Resampling.h>
#include <EmrysFilter/Core/Solve.h>
#include <EmrysFilter/Core/Stamp.h>
#include <EmrysFilter/Core/Weights.h>
#include <EmrysFilter/Distributions/Gaussian.h>
#include <EmrysFilter/Distributions/Moments.h>
#include <EmrysFilter/Distributions/SampleSet.h>
#include <EmrysFilter/Distributions/Statistics.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
 
namespace Emrys
{
 
template <Space Space_, RealScalar Scalar_, typename Allocator_>
class ParticleSet;
 
template <Space Space_, RealScalar Scalar_, typename Allocator_>
class KernelDensity;
 
template <typename Source, typename Selection, typename SubSpace>
class ParticleSetMarginal;
 
namespace internal
{

struct ParticleSetAccess
{
    template <typename Set, typename Derived>
    static void
    setNormalizedLogWeights(Set & set, Eigen::DenseBase<Derived> const & logWeights)
    {
        set.setLogWeights(logWeights);
        set.m_logTotalWeight = typename Set::Scalar(0);
    }

    template <typename Set, typename Rule, typename OnKept>
    [[nodiscard]] static std::expected<typename Set::Scalar, WeightError>
    prune(Set & set, Rule const & rule, OnKept const & onKept)
    {
        return set.pruneParticles(rule, onKept);
    }
};
 
}  // namespace internal
 
namespace internal
{

template <typename MomentsType, ExecutionPolicy P, typename Space, typename StatesDerived, typename WeightsDerived>
[[nodiscard]] typename MomentsType::MatrixType
silvermanBandwidth(
    P && policy,
    Space const & space,
    Eigen::MatrixBase<StatesDerived> const & states,
    Eigen::ArrayBase<WeightsDerived> const & logWeights
)
{
    using Scalar = typename MomentsType::Scalar;
    Scalar const logTotal = logSumExp<Scalar>(policy, logWeights.size(), [&](Eigen::Index i) {
        return logWeights.derived()(i);
    });
    Scalar const d = Scalar(differenceDimension(space));
    Scalar const factor = Eigen::numext::pow(
        Scalar(4) / (ess(policy, logWeights, logTotal) * (d + Scalar(2))),
        Scalar(2) / (d + Scalar(4))
    );
    return factor * weightedMomentsIn<MomentsType>(policy, space, states, logWeights, logTotal).covarianceLower();
}

template <typename MomentsType, typename Space, typename StatesDerived, typename WeightsDerived>
[[nodiscard]] typename MomentsType::MatrixType
silvermanBandwidth(
    Space const & space,
    Eigen::MatrixBase<StatesDerived> const & states,
    Eigen::ArrayBase<WeightsDerived> const & logWeights
)
{
    return silvermanBandwidth<MomentsType>(execution::sequential, space, states, logWeights);
}

template <
    typename Source,
    typename BandwidthDerived,
    typename Target,
    typename Given,
    typename Derived,
    typename TargetSpace,
    typename GivenSpace>
[[nodiscard]] auto
conditionalParticlesIn(
    Source const & source,
    Eigen::MatrixBase<BandwidthDerived> const & bandwidth,
    Target const & target,
    Given const & given,
    Eigen::MatrixBase<Derived> const & values,
    TargetSpace const & targetSpace,
    GivenSpace const & givenSpace
)
{
    using Scalar = typename Source::Scalar;
    EMRYS_FILTER_ASSERT(distinctIndices(target.values, given.values) && "target and given indices must be distinct");
    EMRYS_FILTER_ASSERT(values.size() == given.values.size());
    constexpr int TargetSize = valueDimensionAtCompileTime<DifferenceSpace<TargetSpace>>;
    constexpr int TargetMaxSize = maxValueDimensionAtCompileTime<DifferenceSpace<TargetSpace>>;
    constexpr int GivenSize = valueDimensionAtCompileTime<DifferenceSpace<GivenSpace>>;
    constexpr int GivenMaxSize = maxValueDimensionAtCompileTime<DifferenceSpace<GivenSpace>>;
    using GivenMatrix = MatrixOf<Scalar, GivenSize, GivenSize, GivenMaxSize, GivenMaxSize>;
    using CrossMatrix = MatrixOf<Scalar, TargetSize, GivenSize, TargetMaxSize, GivenMaxSize>;
    using GivenOffsets = SampleSet<DifferenceSpace<GivenSpace>, Scalar>;
    using TargetStates = SampleSet<TargetSpace, Scalar>;
 
    using Result = ParticleSet<TargetSpace, Scalar, typename Source::allocator_type>;
 
    auto const llt =
        Eigen::LLT<GivenMatrix>{GivenMatrix{lowerSymmetricBlock(bandwidth, given.differences, given.differences)}};
    if (llt.info() != Eigen::Success)
    {
        return std::expected<Result, GaussianError>{std::unexpected{GaussianError::NotPositiveDefinite}};
    }
    CrossMatrix const cross = lowerSymmetricBlock(bandwidth, target.differences, given.differences);
    CrossMatrix const gain = transposedSolve<CrossMatrix>(llt, cross.transpose());
 
    GivenOffsets offsets = Emrys::minus(givenSpace, source.states()(given.values, Eigen::placeholders::all), values);
    GivenOffsets solved = offsets;
    solveInPlace(llt, solved);
    Eigen::Array<Scalar, 1, Eigen::Dynamic> const squaredDistance = (offsets.array() * solved.array()).colwise().sum();
    Scalar const log2Pi = Emrys::log2Pi<Scalar>();
    Scalar const logNormalizer =
        Scalar(-0.5)
        * (Scalar(given.differences.size()) * log2Pi + Scalar(2) * llt.matrixLLT().diagonal().array().log().sum());
 
    TargetStates const centers =
        Emrys::plus(targetSpace, source.states()(target.values, Eigen::placeholders::all), -(gain * offsets));
    auto result = Result{
        targetSpace,
        centers,
        Eigen::Array<Scalar, Eigen::Dynamic, 1>{
            source.logWeights() + logNormalizer - Scalar(0.5) * squaredDistance.transpose()
        },
        source.get_allocator()
    };
    result.setBandwidth(
        lowerSymmetricBlock(bandwidth, target.differences, target.differences) - multiply(gain, cross.transpose())
    );
    return std::expected<Result, GaussianError>{std::move(result)};
}

template <typename GivenPath, typename Source, typename BandwidthDerived, typename Derived>
[[nodiscard]] auto
tagConditionalParticles(
    Source const & source,
    Eigen::MatrixBase<BandwidthDerived> const & bandwidth,
    Eigen::MatrixBase<Derived> const & values
)
{
    return tagConditional<GivenPath>(
        source.space(),
        [&](auto const & target, auto const & given, auto const & targetSpace, auto const & givenSpace) {
            return conditionalParticlesIn(source, bandwidth, target, given, values, targetSpace, givenSpace);
        }
    );
}
 
}  // namespace internal

template <Space Space_, RealScalar Scalar_ = double, typename Allocator_ = DefaultAllocator>
class ParticleSet
{
public:
 
    using MomentsType = Moments<Space_, Scalar_>;                             
    using Scalar = typename MomentsType::Scalar;                              
    using Index = typename MomentsType::Index;                                
    using VectorType = typename MomentsType::VectorType;                      
    using DifferenceVectorType = typename MomentsType::DifferenceVectorType;  
    using MatrixType = typename MomentsType::MatrixType;  

    using SampleSetType = SampleSet<Space_, Scalar_>;

    using StatesType = ConstSampleSetMap<Space_, Scalar_>;

    using ConstStateType = Eigen::Block<StatesType const, SampleSetType::RowsAtCompileTime, 1, true>;
 
    using WeightArrayType = Eigen::Array<Scalar, Eigen::Dynamic, 1>;  

    using LogWeightsType = Eigen::Map<WeightArrayType const>;
 
    using SpaceType = Space_;                                        
    using allocator_type = Allocator_;                               
    using DensityType = KernelDensity<Space_, Scalar_, Allocator_>;  

    static constexpr int ValueDimensionAtCompileTime = MomentsType::ValueDimensionAtCompileTime;

    static constexpr int MaxValueDimensionAtCompileTime = MomentsType::MaxValueDimensionAtCompileTime;

    static constexpr int DifferenceDimensionAtCompileTime = MomentsType::DifferenceDimensionAtCompileTime;

    static constexpr int MaxDifferenceDimensionAtCompileTime = MomentsType::MaxDifferenceDimensionAtCompileTime;
 
    static constexpr bool IsParams = true;  

    template <bool ReturnRemovedWeight>
    using RemovalResult = std::conditional_t<ReturnRemovedWeight, std::expected<Scalar, WeightError>, void>;

    class MutableParticle
    {
    public:
 
        MutableParticle(MutableParticle const &) = delete;  

        MutableParticle &
        operator=(MutableParticle const &) = delete;

        [[nodiscard]] Index
        index() const
        {
            return m_index;
        }

        [[nodiscard]] Scalar
        logWeight() const
        {
            return m_set.logWeight(m_index);
        }

        [[nodiscard]] ConstStateType
        state() const
        {
            return m_set.state(m_index);
        }

        template <typename Derived>
        void
        setState(Eigen::MatrixBase<Derived> const & state)
        {
            EMRYS_FILTER_ASSERT(state.size() == m_set.valueDimension());
            m_set.statesMap().col(m_index) = state;
            normalizeValues(m_set.m_space, m_set.statesMap().col(m_index));
        }
 
    private:
 
        friend ParticleSet;  

        MutableParticle(ParticleSet & set, Index index)
            : m_set{set}
            , m_index{index}
        {}
 
        ParticleSet & m_set;  
        Index m_index;        
    };
 
    ParticleSet() = default;

    explicit ParticleSet(Allocator_ const & allocator)
        : m_states(StateAllocator(allocator))
        , m_logWeights(StateAllocator(allocator))
        , m_stateScratch(StateAllocator(allocator))
        , m_weightScratch(StateAllocator(allocator))
        , m_resamplerScratch(StateAllocator(allocator))
        , m_ancestors(IndexAllocator(allocator))
    {
        m_rows = Emrys::valueDimension(m_space);
    }

    explicit ParticleSet(Capacity capacity, Allocator_ const & allocator = {})
        requires FixedSizeSpace<Space_>
        : ParticleSet{allocator}
    {
        setCapacity(capacity.count);
    }

    explicit ParticleSet(Index valueDimension, Allocator_ const & allocator = {})
        requires internal::SizedByLength<Space_>
        : ParticleSet{allocator}
    {
        m_rows = valueDimension;
        fitSpace();
    }

    template <typename StatesDerived>
        requires internal::SizedByLength<Space_>
    explicit ParticleSet(Eigen::MatrixBase<StatesDerived> const & states, Allocator_ const & allocator = {})
        : ParticleSet{allocator}
    {
        assign(states);
    }

    template <typename StatesDerived, typename WeightsDerived>
        requires internal::SizedByLength<Space_>
    ParticleSet(
        Eigen::MatrixBase<StatesDerived> const & states,
        Eigen::DenseBase<WeightsDerived> const & logWeights,
        Allocator_ const & allocator = {}
    )
        : ParticleSet{allocator}
    {
        assign(states, logWeights);
    }

    explicit ParticleSet(Space_ space, Allocator_ const & allocator = {})
        : ParticleSet{allocator}
    {
        m_space = std::move(space);
        m_rows = Emrys::valueDimension(m_space);
    }

    ParticleSet(Space_ space, Capacity capacity, Allocator_ const & allocator = {})
        : ParticleSet{std::move(space), allocator}
    {
        setCapacity(capacity.count);
    }

    ParticleSet(Space_ space, Index valueDimension, Allocator_ const & allocator = {})
        : ParticleSet{allocator}
    {
        m_space = std::move(space);
        m_rows = valueDimension;
        fitSpace();
    }

    template <typename StatesDerived>
    ParticleSet(Space_ space, Eigen::MatrixBase<StatesDerived> const & states, Allocator_ const & allocator = {})
        : ParticleSet{allocator}
    {
        m_space = std::move(space);
        assign(states);
    }

    template <typename StatesDerived, typename WeightsDerived>
    ParticleSet(
        Space_ space,
        Eigen::MatrixBase<StatesDerived> const & states,
        Eigen::DenseBase<WeightsDerived> const & logWeights,
        Allocator_ const & allocator = {}
    )
        : ParticleSet{allocator}
    {
        m_space = std::move(space);
        assign(states, logWeights);
    }

    ParticleSet(ParticleSet const & other)
        : ParticleSet{other, other.get_allocator()}
    {}

    ParticleSet(ParticleSet const & other, Allocator_ const & allocator)
        : ParticleSet{allocator}
    {
        m_capacity = other.m_capacity;
        reserveRows(other.m_rows);
        copyFrom(other);
    }

    ParticleSet(ParticleSet && other) noexcept
        : m_rows{other.m_rows}
        , m_states{std::move(other.m_states)}
        , m_logWeights{std::move(other.m_logWeights)}
        , m_space{std::move(other.m_space)}
        , m_bandwidth{std::move(other.m_bandwidth)}
        , m_capacity{other.m_capacity}
        , m_stateScratch{std::move(other.m_stateScratch)}
        , m_weightScratch{std::move(other.m_weightScratch)}
        , m_resamplerScratch{std::move(other.m_resamplerScratch)}
        , m_ancestors{std::move(other.m_ancestors)}
        , m_resampleCount{other.m_resampleCount}
        , m_logTotalWeight{other.m_logTotalWeight}
    {
        other.touchShape();
    }

    ParticleSet(ParticleSet && other, Allocator_ const & allocator)
        : m_rows{other.m_rows}
        , m_states{std::move(other.m_states), StateAllocator(allocator)}
        , m_logWeights{std::move(other.m_logWeights), StateAllocator(allocator)}
        , m_space{std::move(other.m_space)}
        , m_bandwidth{std::move(other.m_bandwidth)}
        , m_capacity{other.m_capacity}
        , m_stateScratch{std::move(other.m_stateScratch), StateAllocator(allocator)}
        , m_weightScratch{std::move(other.m_weightScratch), StateAllocator(allocator)}
        , m_resamplerScratch{std::move(other.m_resamplerScratch), StateAllocator(allocator)}
        , m_ancestors{std::move(other.m_ancestors), IndexAllocator(allocator)}
        , m_resampleCount{other.m_resampleCount}
        , m_logTotalWeight{other.m_logTotalWeight}
    {
        other.touchShape();
    }

    ParticleSet &
    operator=(ParticleSet const & other)
    {
        touchShape();
        if (this != &other)
        {
            copyFrom(other);
        }
        return *this;
    }

    ParticleSet &
    operator=(ParticleSet && other) noexcept
    {
        touchShape();
        other.touchShape();
        if (this == &other)
        {
            return *this;
        }
        if constexpr (std::allocator_traits<Allocator_>::propagate_on_container_move_assignment::value)
        {
            takeStorage(other);
        }
        else
        {
            if (get_allocator() == other.get_allocator())
            {
                takeStorage(other);
            }
            else
            {
                copyFrom(other);
                m_ancestors.assign(other.m_ancestors.begin(), other.m_ancestors.end());
                m_resampleCount = other.m_resampleCount;
            }
        }
        return *this;
    }

    [[nodiscard]] allocator_type
    get_allocator() const
    {
        return allocator_type(m_states.get_allocator());
    }

    [[nodiscard]] Index
    size() const
    {
        return Index(m_logWeights.size());
    }

    [[nodiscard]] bool
    empty() const
    {
        return size() == 0;
    }

    [[nodiscard]] std::optional<Index>
    capacity() const
    {
        return m_capacity;
    }

    void
    setCapacity(Index capacity)
    {
        EMRYS_FILTER_ASSERT(capacity >= size() && "the capacity is smaller than the set");
        m_capacity = capacity;
        reserveRows(m_rows);
    }

    [[nodiscard]] Index
    valueDimension() const
    {
        return m_rows;
    }

    [[nodiscard]] Index
    differenceDimension() const
    {
        return Emrys::differenceDimension(m_space);
    }

    [[nodiscard]] Space_ const &
    space() const
    {
        return m_space;
    }

    [[nodiscard]] StatesType
    states() const
    {
        return StatesType{m_states.data(), m_rows, size()};
    }

    [[nodiscard]] LogWeightsType
    logWeights() const
    {
        return LogWeightsType{m_logWeights.data(), size()};
    }

    [[nodiscard]] ConstStateType
    state(Index i) const
    {
        StatesType const all = states();
        return ConstStateType{all, i};
    }

    [[nodiscard]] Scalar
    logWeight(Index i) const
    {
        assertIndex(i);
        return m_logWeights[std::size_t(i)];
    }

    [[nodiscard]] std::optional<MatrixType> const &
    bandwidth() const
    {
        return m_bandwidth;
    }

    [[nodiscard]] std::expected<MatrixType, WeightError>
    kdeBandwidth() const
    {
        return kdeBandwidth(execution::sequential);
    }

    template <ExecutionPolicy P>
    [[nodiscard]] std::expected<MatrixType, WeightError>
    kdeBandwidth(P && policy) const
    {
        if (m_bandwidth)
        {
            return *m_bandwidth;
        }
        return silvermanBandwidth(std::forward<P>(policy), *this);
    }

    template <MatrixExpression Bandwidth>
    void
    setBandwidth(Bandwidth && bandwidth)
    {
        m_stamp.touchValues();
        EMRYS_FILTER_ASSERT(bandwidth.rows() == differenceDimension() && bandwidth.cols() == differenceDimension());
        m_bandwidth = MatrixType{std::forward<Bandwidth>(bandwidth)};
    }

    void
    clearBandwidth()
    {
        m_stamp.touchValues();
        m_bandwidth.reset();
    }

    void
    resizeParticles(Index n)
        requires internal::SizedByLength<Space_>
    {
        touchShape();
        internal::fitValueDimension(m_space, n);
        m_rows = n;
        reserveRows(m_rows);
        m_states.resize(std::size_t(m_rows * size()));
        m_bandwidth.reset();
    }

    template <typename Path, typename Derived>
    void
    insert(Index element, Eigen::MatrixBase<Derived> const & states)
    {
        touchShape();
        EMRYS_FILTER_ASSERT(states.cols() == size() && "insert one column per particle");
        auto const rows = internal::growBlock<Path>(m_space, element, states.rows());
        Index const newRows = m_rows + states.rows();
        reserveRows(newRows);
        m_stateScratch.resize(std::size_t(newRows * size()));
        auto grown = SampleSetMap<Space_, Scalar_>{m_stateScratch.data(), newRows, size()};
        grown.topRows(rows.start) = this->states().topRows(rows.start);
        grown.middleRows(rows.start, rows.count) = states;
        grown.bottomRows(m_rows - rows.start) = this->states().bottomRows(m_rows - rows.start);
        m_states.swap(m_stateScratch);
        m_rows = newRows;
        normalizeValues(m_space, statesMap());
        m_bandwidth.reset();
    }

    template <typename Path>
    void
    erase(Index element, Index count = 1)
    {
        touchShape();
        auto const rows = internal::shrinkBlock<Path>(m_space, element, count);
        Index const newRows = m_rows - rows.count;
        m_stateScratch.resize(std::size_t(newRows * size()));
        auto shrunk = SampleSetMap<Space_, Scalar_>{m_stateScratch.data(), newRows, size()};
        shrunk.topRows(rows.start) = states().topRows(rows.start);
        shrunk.bottomRows(newRows - rows.start) = states().bottomRows(newRows - rows.start);
        m_states.swap(m_stateScratch);
        m_rows = newRows;
        if (m_bandwidth)
        {
            internal::eraseSymmetric(*m_bandwidth, rows.differenceStart, rows.differenceCount);
        }
    }

    void
    clear()
    {
        touchShape();
        m_states.clear();
        m_logWeights.clear();
    }

    template <typename StatesDerived>
    void
    assign(Eigen::MatrixBase<StatesDerived> const & states)
    {
        touchShape();
        copyStates(states);
        logWeightsMap().setConstant(-Eigen::numext::log(Scalar(size())));
        fitSpace();
    }

    template <typename StatesDerived, typename WeightsDerived>
    void
    assign(Eigen::MatrixBase<StatesDerived> const & states, Eigen::DenseBase<WeightsDerived> const & logWeights)
    {
        touchShape();
        EMRYS_FILTER_ASSERT(logWeights.size() == states.cols());
        copyStates(states);
        logWeightsMap() = logWeights.derived().array();
        fitSpace();
    }

    template <typename Derived>
    Index
    add(Scalar logWeight, Eigen::MatrixBase<Derived> const & state)
    {
        touchShape();
        EMRYS_FILTER_ASSERT(state.size() == valueDimension());
        // Copied first: growing the storage would move a state read from this set.
        VectorType const value = state;
        Index const i = size();
        setCount(i + 1);
        statesMap().col(i) = value;
        normalizeValues(m_space, statesMap().col(i));
        m_logWeights[std::size_t(i)] = logWeight;
        return i;
    }

    template <typename OtherAllocator>
    void
    append(ParticleSet<Space_, Scalar_, OtherAllocator> const & other)
    {
        touchShape();
        EMRYS_FILTER_ASSERT(other.space() == m_space && "cannot append a set in a different space");
        Index const n = size();
        Index const m = other.size();
        // Read after growing, which may move this set's storage when other is this set.
        auto const * source = &other;
        setCount(n + m);
        statesMap().rightCols(m) = source->states().leftCols(m);
        logWeightsMap().tail(m) = source->logWeights().head(m);
    }

    void
    remove(Index i)
    {
        touchShape();
        assertIndex(i);
        Index const last = size() - 1;
        if (i != last)
        {
            statesMap().col(i) = states().col(last);
            m_logWeights[std::size_t(i)] = m_logWeights[std::size_t(last)];
        }
        setCount(last);
    }

    template <typename Derived>
    void
    setState(Index i, Eigen::MatrixBase<Derived> const & state)
    {
        m_stamp.touchValues();
        assertIndex(i);
        EMRYS_FILTER_ASSERT(state.size() == valueDimension());
        statesMap().col(i) = state;
        normalizeValues(m_space, statesMap().col(i));
    }

    template <typename Derived>
    void
    setStates(Eigen::MatrixBase<Derived> const & states)
    {
        m_stamp.touchValues();
        EMRYS_FILTER_ASSERT(states.rows() == valueDimension() && states.cols() == size());
        statesMap() = states;
        normalizeValues(m_space, statesMap());
    }

    void
    setLogWeight(Index i, Scalar logWeight)
    {
        touchWeights();
        assertIndex(i);
        m_logWeights[std::size_t(i)] = logWeight;
    }

    template <typename Derived>
    void
    setLogWeights(Eigen::DenseBase<Derived> const & logWeights)
    {
        touchWeights();
        EMRYS_FILTER_ASSERT(logWeights.size() == size());
        logWeightsMap() = logWeights.derived().array();
    }

    void
    scaleWeight(Index i, Scalar logScale)
    {
        touchWeights();
        assertIndex(i);
        m_logWeights[std::size_t(i)] += logScale;
    }

    void
    scaleWeights(Scalar logScale)
    {
        touchWeights();
        logWeightsMap() += logScale;
    }

    template <typename Derived>
    void
    scaleWeights(Eigen::DenseBase<Derived> const & logScales)
    {
        touchWeights();
        EMRYS_FILTER_ASSERT(logScales.size() == size());
        logWeightsMap() += logScales.derived().array();
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    logTotalWeight() const
    {
        return logTotalWeight(execution::sequential);
    }

    template <ExecutionPolicy P>
    [[nodiscard]] std::expected<Scalar, WeightError>
    logTotalWeight(P && policy) const
    {
        if (m_logTotalWeight)
        {
            return *m_logTotalWeight;
        }
        return internal::logTotalWeight<Scalar>(
            std::forward<P>(policy),
            std::span<Scalar const>{m_logWeights.data(), m_logWeights.size()}
        );
    }

    [[nodiscard]] std::expected<void, WeightError>
    normalize()
    {
        return normalize(execution::sequential);
    }

    template <ExecutionPolicy P>
    [[nodiscard]] std::expected<void, WeightError>
    normalize(P && policy)
    {
        touchWeights();
        auto const total = logTotalWeight(std::forward<P>(policy));
        if (!total)
        {
            return std::unexpected{total.error()};
        }
        scaleWeights(-*total);
        m_logTotalWeight = Scalar(0);
        return {};
    }

    [[nodiscard]] std::expected<WeightArrayType, WeightError>
    normalizedWeights() const
    {
        return logTotalWeight().transform([this](Scalar total) {
            return WeightArrayType{(logWeights() - total).exp()};
        });
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    ess() const
    {
        return ess(execution::sequential);
    }

    template <ExecutionPolicy P>
    [[nodiscard]] std::expected<Scalar, WeightError>
    ess(P && policy) const
    {
        return logTotalWeight(policy).transform([&](Scalar total) {
            return internal::ess(policy, logWeights(), total);
        });
    }

    [[nodiscard]] std::expected<MomentsType, WeightError>
    moments() const
    {
        return moments(execution::sequential);
    }

    template <ExecutionPolicy P>
    [[nodiscard]] std::expected<MomentsType, WeightError>
    moments(P && policy) const
    {
        return logTotalWeight(policy).transform([&](Scalar total) {
            return internal::weightedMomentsIn<MomentsType>(policy, m_space, states(), logWeights(), total);
        });
    }

    template <typename F>
    void
    forEach(F && f)
    {
        forEach(execution::sequential, std::forward<F>(f));
    }

    template <ExecutionPolicy P, typename F>
    void
    forEach(P && policy, F && f)
    {
        static_assert(
            std::is_void_v<std::invoke_result_t<F &, MutableParticle &>>,
            "forEach() takes a callback that returns nothing; one that returns a log scale factor goes to reweight()"
        );
        visitParticles(std::forward<P>(policy), [&](MutableParticle & particle, Index) {
            std::invoke(f, particle);
        });
    }

    template <typename F>
    void
    reweight(F && f)
    {
        reweight(execution::sequential, std::forward<F>(f));
    }

    template <ExecutionPolicy P, typename F>
    void
    reweight(P && policy, F && f)
    {
        static_assert(
            std::floating_point<std::remove_cvref_t<std::invoke_result_t<F &, MutableParticle &>>>,
            "reweight() takes a callback that returns a log scale factor, a floating-point number"
        );
        visitParticles(std::forward<P>(policy), [&](MutableParticle & particle, Index i) {
            m_logWeights[std::size_t(i)] += Scalar(std::invoke(f, particle));
        });
    }

    template <bool ReturnRemovedWeight = false, typename Pred>
        requires(!ReturnRemovedWeight && !ExecutionPolicy<Pred>)
    void
    removeIf(Pred pred)
    {
        removeParticlesIf<false>(execution::sequential, std::move(pred));
    }

    template <bool ReturnRemovedWeight, typename Pred>
        requires(ReturnRemovedWeight && !ExecutionPolicy<Pred>)
    [[nodiscard]] std::expected<Scalar, WeightError>
    removeIf(Pred pred)
    {
        return removeParticlesIf<true>(execution::sequential, std::move(pred));
    }

    template <bool ReturnRemovedWeight = false, ExecutionPolicy P, typename Pred>
        requires(!ReturnRemovedWeight)
    void
    removeIf(P && policy, Pred pred)
    {
        removeParticlesIf<false>(std::forward<P>(policy), std::move(pred));
    }

    template <bool ReturnRemovedWeight, ExecutionPolicy P, typename Pred>
        requires ReturnRemovedWeight
    [[nodiscard]] std::expected<Scalar, WeightError>
    removeIf(P && policy, Pred pred)
    {
        return removeParticlesIf<true>(std::forward<P>(policy), std::move(pred));
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    pruneBelow(Scalar logThreshold)
    {
        return pruneParticles(internal::PruneBelowRule<Scalar>{logThreshold}, [](Index) {});
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    pruneTail(Scalar logFraction)
    {
        return pruneParticles(internal::PruneTailRule<Scalar>{logFraction}, [](Index) {});
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    pruneToSize(Index maxSize)
    {
        return pruneParticles(internal::PruneToSizeRule{maxSize}, [](Index) {});
    }

    template <std::uniform_random_bit_generator G, typename R = resamplers::Systematic>
        requires Resampler<R, Scalar, G>
    [[nodiscard]] std::expected<void, WeightError>
    resample(G & rng, Index newSize, R const & resampler = {})
    {
        return resample(execution::sequential, rng, newSize, resampler);
    }

    template <std::uniform_random_bit_generator G, typename R = resamplers::Systematic>
        requires Resampler<R, Scalar, G>
    [[nodiscard]] std::expected<void, WeightError>
    resample(G & rng, R const & resampler = {})
    {
        return resample(execution::sequential, rng, size(), resampler);
    }

    template <ExecutionPolicy P, std::uniform_random_bit_generator G, typename R = resamplers::Systematic>
        requires Resampler<R, Scalar, G>
    [[nodiscard]] std::expected<void, WeightError>
    resample(P && policy, G & rng, Index newSize, R const & resampler = {})
    {
        touchShape();
        EMRYS_FILTER_ASSERT(newSize > 0);
        assertCapacity(newSize);
        auto const total = logTotalWeight(policy);
        if (!total)
        {
            return std::unexpected{total.error()};
        }
 
        m_weightScratch.resize(std::size_t(size()));
        auto linearWeights = Eigen::Map<WeightArrayType>{m_weightScratch.data(), size()};
        internal::forEachChunk(policy, size(), [&](Index begin, Index length) {
            linearWeights.segment(begin, length) = (logWeights().segment(begin, length) - *total).exp();
        });
        m_ancestors.resize(std::size_t(newSize));
        auto const weights = std::span<Scalar const>{m_weightScratch};
        auto const ancestors = std::span<Index>{m_ancestors};
        if constexpr (ResamplerUnder<R, P, Scalar, G> && !internal::isSequential<P>)
        {
            resampler(policy, weights, ancestors, rng);
        }
        else if constexpr (ResamplerWithScratch<R, Scalar, G>)
        {
            m_resamplerScratch.resize(std::size_t(size() + newSize));
            resampler(weights, ancestors, rng, std::span<Scalar>{m_resamplerScratch});
        }
        else
        {
            resampler(weights, ancestors, rng);
        }
 
        m_stateScratch.resize(std::size_t(m_rows * newSize));
        internal::gather(
            policy,
            std::span<Index const>{ancestors},
            states(),
            SampleSetMap<Space_, Scalar_>{m_stateScratch.data(), m_rows, newSize}
        );
        m_states.swap(m_stateScratch);
        m_logWeights.resize(std::size_t(newSize));
        internal::fill(std::forward<P>(policy), std::span<Scalar>{m_logWeights}, -Eigen::numext::log(Scalar(newSize)));
        m_logTotalWeight = Scalar(0);
        m_bandwidth.reset();
        ++m_resampleCount;
        return {};
    }

    template <ExecutionPolicy P, std::uniform_random_bit_generator G, typename R = resamplers::Systematic>
        requires Resampler<R, Scalar, G>
    [[nodiscard]] std::expected<void, WeightError>
    resample(P && policy, G & rng, R const & resampler = {})
    {
        return resample(std::forward<P>(policy), rng, size(), resampler);
    }

    [[nodiscard]] std::span<Index const>
    ancestors() const
    {
        return {m_ancestors.data(), m_ancestors.size()};
    }

    [[nodiscard]] Index
    resampleCount() const
    {
        return m_resampleCount;
    }

    template <typename Path>
    [[nodiscard]] auto
    marginal() const &
    {
        return internal::tagMarginal<ParticleSetMarginal, Path>(*this);
    }

    template <typename Path>
    void
    marginal() const && = delete;

    template <typename Path>
    [[nodiscard]] auto
    marginal(Index element) const &
    {
        return internal::elementMarginal<ParticleSetMarginal, Path>(*this, element);
    }

    template <typename Path>
    void
    marginal(Index element) const && = delete;

    template <typename GivenPath, typename Derived>
    [[nodiscard]] auto
    conditional(Eigen::MatrixBase<Derived> const & values) const
    {
        auto const bandwidth = kdeBandwidth();
        using Result = decltype(internal::tagConditionalParticles<GivenPath>(*this, *bandwidth, values));
        if (!bandwidth)
        {
            return Result{std::unexpected{GaussianError::NonFinite}};
        }
        return internal::tagConditionalParticles<GivenPath>(*this, *bandwidth, values);
    }

    [[nodiscard]] std::expected<DensityType, NumericError>
    toDensity() const &
    {
        auto const total = logTotalWeight();
        if (!total)
        {
            return std::unexpected{total.error()};
        }
        // The weights are valid, so the bandwidth is found.
        return kernelDensity(*total, *kdeBandwidth());
    }

    void
    toDensity() const && = delete;

    template <MatrixExpression Bandwidth>
    [[nodiscard]] std::expected<DensityType, NumericError>
    toDensity(Bandwidth && bandwidth) const &
    {
        auto const total = logTotalWeight();
        if (!total)
        {
            return std::unexpected{total.error()};
        }
        return kernelDensity(*total, std::forward<Bandwidth>(bandwidth));
    }

    template <MatrixExpression Bandwidth>
    void
    toDensity(Bandwidth && bandwidth) const && = delete;
 
private:

    template <typename>
    friend class internal::ShapeCheck;

    template <typename>
    friend class internal::ValuesCheck;

    friend struct internal::ParticleSetAccess;

    template <bool ReturnRemovedWeight, ExecutionPolicy P, typename Pred>
    RemovalResult<ReturnRemovedWeight>
    removeParticlesIf(P && policy, Pred pred)
    {
        touchShape();
        [[maybe_unused]] Scalar logTotal = 0;
        if constexpr (ReturnRemovedWeight)
        {
            auto const total = logTotalWeight(policy);
            if (!total)
            {
                return std::unexpected{total.error()};
            }
            logTotal = *total;
        }
        auto const guard = internal::CallbackGuard{m_stamp};
        auto const checked = [&](Scalar logWeight, ConstStateType const & state) {
            bool const remove = static_cast<bool>(std::invoke(pred, logWeight, state));
            guard.check();
            return remove;
        };
        [[maybe_unused]] Scalar removed = 0;
        if constexpr (internal::isSequential<P>)
        {
            removed = compact([&](Index i) {
                return checked(m_logWeights[std::size_t(i)], state(i));
            });
        }
        else
        {
            removed = compact(policy, checked, ReturnRemovedWeight);
        }
        if constexpr (ReturnRemovedWeight)
        {
            return removed - logTotal;
        }
    }

    template <typename Rule, typename OnKept>
    [[nodiscard]] std::expected<Scalar, WeightError>
    pruneParticles(Rule const & rule, OnKept const & onKept)
    {
        auto const total = logTotalWeight();
        if (!total)
        {
            return std::unexpected{total.error()};
        }
        touchShape();
        auto cut = rule(
            size(),
            [&](Index i) {
                return m_logWeights[std::size_t(i)];
            },
            *total,
            m_weightScratch
        );
        Scalar const removed = compact([&](Index i) {
            bool const remove = cut.removes(m_logWeights[std::size_t(i)]);
            if (!remove)
            {
                onKept(i);
            }
            return remove;
        });
        // The kept particles include one of positive weight, so their total is finite.
        static_cast<void>(normalize());
        return removed - *total;
    }

    template <ExecutionPolicy P, typename Visit>
    void
    visitParticles(P && policy, Visit const & visit)
    {
        touchWeights();
        auto const guard = internal::CallbackGuard{m_stamp};
        internal::forEachIndex(std::forward<P>(policy), size(), [&](Index i) {
            auto particle = MutableParticle{*this, i};
            visit(particle, i);
            guard.check();
        });
    }

    void
    touchWeights()
    {
        m_stamp.touchValues();
        m_logTotalWeight.reset();
    }

    void
    touchShape()
    {
        m_stamp.touchShape();
        m_logTotalWeight.reset();
    }
 
    using StateAllocator = internal::ReBound<Allocator_, Scalar>;  
    using IndexAllocator = internal::ReBound<Allocator_, Index>;   

    void
    assertCapacity([[maybe_unused]] Index count) const
    {
        EMRYS_FILTER_ASSERT(
            (!m_capacity || count <= *m_capacity) && "the set would hold more particles than its capacity"
        );
    }

    [[nodiscard]] SampleSetMap<Space_, Scalar_>
    statesMap()
    {
        return {m_states.data(), m_rows, size()};
    }

    [[nodiscard]] Eigen::Map<WeightArrayType>
    logWeightsMap()
    {
        return {m_logWeights.data(), size()};
    }

    void
    setCount(Index count)
    {
        assertCapacity(count);
        m_states.resize(std::size_t(m_rows * count));
        m_logWeights.resize(std::size_t(count));
    }

    void
    reserveRows(Index rows)
    {
        if (!m_capacity)
        {
            return;
        }
        Index const bound = maxValueDimensionAtCompileTime<Space_> == Eigen::Dynamic
                              ? rows
                              : Index(maxValueDimensionAtCompileTime<Space_>);
        auto const count = std::size_t(*m_capacity);
        m_states.reserve(std::size_t(bound) * count);
        m_stateScratch.reserve(std::size_t(bound) * count);
        m_logWeights.reserve(count);
        m_weightScratch.reserve(count);
        m_resamplerScratch.reserve(2 * count);
        m_ancestors.reserve(count);
    }

    template <typename StatesDerived>
    void
    copyStates(Eigen::MatrixBase<StatesDerived> const & states)
    {
        if (states.rows() != m_rows)
        {
            m_bandwidth.reset();
        }
        m_rows = states.rows();
        reserveRows(m_rows);
        setCount(states.cols());
        statesMap() = states;
    }

    void
    copyFrom(ParticleSet const & other)
    {
        m_rows = other.m_rows;
        m_space = other.m_space;
        m_bandwidth = other.m_bandwidth;
        setCount(other.size());
        m_states.assign(other.m_states.begin(), other.m_states.end());
        m_logWeights.assign(other.m_logWeights.begin(), other.m_logWeights.end());
        m_logTotalWeight = other.m_logTotalWeight;
        m_ancestors.clear();
        m_resampleCount = 0;
    }

    void
    takeStorage(ParticleSet & other)
    {
        m_rows = other.m_rows;
        m_states = std::move(other.m_states);
        m_logWeights = std::move(other.m_logWeights);
        m_logTotalWeight = other.m_logTotalWeight;
        m_space = std::move(other.m_space);
        m_bandwidth = std::move(other.m_bandwidth);
        m_capacity = other.m_capacity;
        m_stateScratch = std::move(other.m_stateScratch);
        m_weightScratch = std::move(other.m_weightScratch);
        m_resamplerScratch = std::move(other.m_resamplerScratch);
        m_ancestors = std::move(other.m_ancestors);
        m_resampleCount = other.m_resampleCount;
    }

    template <MatrixExpression Bandwidth>
    [[nodiscard]] std::expected<DensityType, NumericError>
    kernelDensity(Scalar logTotalWeight, Bandwidth && bandwidth) const
    {
        auto kde = DensityType::create(*this, logTotalWeight, std::forward<Bandwidth>(bandwidth));
        if (!kde)
        {
            return std::unexpected{kde.error()};
        }
        return std::move(*kde);
    }

    void
    assertIndex([[maybe_unused]] Index i) const
    {
        EMRYS_FILTER_ASSERT(i >= 0 && i < size() && "particle index out of range");
    }

    void
    fitSpace()
    {
        internal::fitValueDimension(m_space, valueDimension());
        normalizeValues(m_space, statesMap());
    }

    template <typename RemoveAtIndex>
    Scalar
    compact(RemoveAtIndex removeAtIndex)
    {
        auto removed = LogSumExp<Scalar>{};
        Index out = 0;
        for (Index i = 0; i < size(); ++i)
        {
            if (removeAtIndex(i))
            {
                removed.add(m_logWeights[std::size_t(i)]);
                continue;
            }
            if (out != i)
            {
                statesMap().col(out) = states().col(i);
                m_logWeights[std::size_t(out)] = m_logWeights[std::size_t(i)];
            }
            ++out;
        }
        setCount(out);
        return removed.value();
    }

    template <ExecutionPolicy P, typename Pred>
    Scalar
    compact(P && policy, Pred const & pred, bool sumRemoved)
    {
        Index const count = size();
        auto remove = internal::scratchVector<char>(policy, static_cast<std::size_t>(count));
        internal::forEachIndex(policy, count, [&](Index i) {
            remove[std::size_t(i)] = static_cast<bool>(std::invoke(pred, m_logWeights[std::size_t(i)], state(i)));
        });
        Scalar const removed = sumRemoved ? logSumExp<Scalar>(
                                                policy,
                                                count,
                                                [&](Index i) {
                                                    return remove[std::size_t(i)]
                                                             ? m_logWeights[std::size_t(i)]
                                                             : -Eigen::NumTraits<Scalar>::infinity();
                                                }
                                            )
                                          : -Eigen::NumTraits<Scalar>::infinity();
        m_stateScratch.resize(m_states.size());
        m_weightScratch.resize(m_logWeights.size());
        auto kept = SampleSetMap<Space_, Scalar_>{m_stateScratch.data(), m_rows, count};
        Index const keptCount = internal::copyIfIndex(
            std::forward<P>(policy),
            count,
            [&](Index i) {
                return remove[std::size_t(i)] == 0;
            },
            [&](Index i, Index j) {
                kept.col(j) = states().col(i);
                m_weightScratch[std::size_t(j)] = m_logWeights[std::size_t(i)];
            }
        );
        m_states.swap(m_stateScratch);
        m_logWeights.swap(m_weightScratch);
        setCount(keptCount);
        return removed;
    }
 
    Index m_rows = valueDimensionAtCompileTime<Space_> == Eigen::Dynamic
                     ? 0
                     : valueDimensionAtCompileTime<Space_>;  
    internal::Buffer<Scalar, Allocator_> m_states;           
    internal::Buffer<Scalar, Allocator_> m_logWeights;       
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Space_ m_space;           
    std::optional<MatrixType> m_bandwidth;                   
    std::optional<Index> m_capacity;                         
 
    // Scratch, kept between calls so that steady-size resampling and resizing don't allocate.
    internal::Buffer<Scalar, Allocator_> m_stateScratch;  
 
    internal::Buffer<Scalar, Allocator_> m_weightScratch;     
    internal::Buffer<Scalar, Allocator_> m_resamplerScratch;  
    internal::Buffer<Index, Allocator_> m_ancestors;          
    Index m_resampleCount = 0;                                
    std::optional<Scalar> m_logTotalWeight;                   
    internal::Stamp m_stamp;                                  
};

template <typename Source, typename Selection, typename SubSpace>
class ParticleSetMarginal
{
public:
 
    using Scalar = typename Source::Scalar;  
    using Index = Eigen::Index;              

    static constexpr int ValueDimensionAtCompileTime = valueDimensionAtCompileTime<SubSpace>;

    static constexpr int MaxValueDimensionAtCompileTime = maxValueDimensionAtCompileTime<SubSpace>;

    static constexpr int DifferenceDimensionAtCompileTime = valueDimensionAtCompileTime<DifferenceSpace<SubSpace>>;

    static constexpr int MaxDifferenceDimensionAtCompileTime =
        maxValueDimensionAtCompileTime<DifferenceSpace<SubSpace>>;
 
    using SpaceType = SubSpace;                                        
    using allocator_type = typename Source::allocator_type;            
    using ParamsType = ParticleSet<SubSpace, Scalar, allocator_type>;  
    using MatrixType = typename ParamsType::MatrixType;                
    using WeightArrayType = typename ParamsType::WeightArrayType;      
    using LogWeightsType = typename Source::LogWeightsType;            

    using IndexVectorType = internal::IndexedSize<
        decltype(Selection::values),
        Source::ValueDimensionAtCompileTime,
        Source::MaxValueDimensionAtCompileTime>;

    using DifferenceIndexVectorType = internal::IndexedSize<
        decltype(Selection::differences),
        Source::DifferenceDimensionAtCompileTime,
        Source::MaxDifferenceDimensionAtCompileTime>;

    using StateType = decltype(std::declval<Source const &>().state(Index{})(std::declval<IndexVectorType const &>()));

    struct Element
    {
        Scalar logWeight;  
        StateType state;   
    };
 
    using value_type = Element;                                             
    using const_iterator = internal::ElementIterator<ParticleSetMarginal>;  
    using iterator = const_iterator;  

    ParticleSetMarginal(Source const & source, Selection const & selection, SubSpace space)
        : m_holder{source}
        , m_indices{internal::resolveIndices<
              Source::ValueDimensionAtCompileTime,
              Source::MaxValueDimensionAtCompileTime>(source.valueDimension(), selection.values)}
        , m_differenceIndices{internal::resolveIndices<
              Source::DifferenceDimensionAtCompileTime,
              Source::MaxDifferenceDimensionAtCompileTime>(source.differenceDimension(), selection.differences)}
        , m_space{std::move(space)}
    {
        EMRYS_FILTER_ASSERT(
            internal::distinctIndices(m_indices) && internal::distinctIndices(m_differenceIndices)
            && "marginal indices must be distinct"
        );
    }

    [[nodiscard]] SubSpace const &
    space() const
    {
        return m_space;
    }

    [[nodiscard]] allocator_type
    get_allocator() const
    {
        return source().get_allocator();
    }

    [[nodiscard]] Index
    size() const
    {
        return source().size();
    }

    [[nodiscard]] Index
    valueDimension() const
    {
        return m_indices.size();
    }

    [[nodiscard]] Index
    differenceDimension() const
    {
        return m_differenceIndices.size();
    }

    [[nodiscard]] IndexVectorType const &
    indices() const
    {
        return m_indices;
    }

    [[nodiscard]] DifferenceIndexVectorType const &
    differenceIndices() const
    {
        return m_differenceIndices;
    }

    [[nodiscard]] auto
    states() const
    {
        return source().states()(m_indices, Eigen::placeholders::all);
    }

    [[nodiscard]] auto
    state(Index i) const
    {
        return source().state(i)(m_indices);
    }

    [[nodiscard]] LogWeightsType
    logWeights() const
    {
        return source().logWeights();
    }

    [[nodiscard]] Scalar
    logWeight(Index i) const
    {
        return source().logWeight(i);
    }

    [[nodiscard]] Element
    operator[](Index i) const
    {
        return Element{logWeight(i), state(i)};
    }

    [[nodiscard]] const_iterator
    begin() const
    {
        return const_iterator{*this, 0};
    }

    [[nodiscard]] const_iterator
    end() const
    {
        return const_iterator{*this, size()};
    }

    [[nodiscard]] std::optional<MatrixType>
    bandwidth() const
    {
        auto const & full = source().bandwidth();
        if (!full)
        {
            return std::nullopt;
        }
        return MatrixType{internal::lowerSymmetricBlock(*full, m_differenceIndices, m_differenceIndices)};
    }

    [[nodiscard]] ParamsType
    copy() const
    {
        return copy(get_allocator());
    }

    [[nodiscard]] ParamsType
    copy(allocator_type const & allocator) const
    {
        auto result = ParamsType{m_space, states(), logWeights(), allocator};
        if (auto const h = bandwidth())
        {
            result.setBandwidth(*h);
        }
        return result;
    }
 
private:

    [[nodiscard]] Source const &
    source() const
    {
        return m_holder.get();
    }
 
    internal::HolderFor<Source> m_holder;           
    IndexVectorType m_indices;                      
    DifferenceIndexVectorType m_differenceIndices;  
    SubSpace m_space;                               
};

template <typename Scalar, typename Space, typename Allocator>
[[nodiscard]] std::expected<typename ParticleSet<Space, Scalar, Allocator>::MatrixType, WeightError>
silvermanBandwidth(ParticleSet<Space, Scalar, Allocator> const & particles)
{
    return silvermanBandwidth(execution::sequential, particles);
}

template <ExecutionPolicy P, typename Scalar, typename Space, typename Allocator>
[[nodiscard]] std::expected<typename ParticleSet<Space, Scalar, Allocator>::MatrixType, WeightError>
silvermanBandwidth(P && policy, ParticleSet<Space, Scalar, Allocator> const & particles)
{
    using MomentsType = typename ParticleSet<Space, Scalar, Allocator>::MomentsType;
    return particles.logTotalWeight(policy).transform([&](Scalar) {
        return internal::silvermanBandwidth<MomentsType>(
            policy,
            particles.space(),
            particles.states(),
            particles.logWeights()
        );
    });
}
 
}  // namespace Emrys
 
// KernelDensity completes ParticleSet::toDensity(); it includes this header first.
#include <EmrysFilter/Distributions/KernelDensity.h>
 
#ifndef EMRYS_FILTER_NO_IO
 
namespace Emrys
{
 
namespace internal
{

template <typename Space_, typename Scalar_, typename Allocator_>
void
describe(std::ostream & stream, ParticleSet<Space_, Scalar_, Allocator_> const & particles, Indent indent)
{
    stream << indent << "particle set of ";
    printCount(stream, particles.size(), "particle");
    printDimension(stream, particles.valueDimension(), particles.differenceDimension());
    if (auto const ess = particles.ess())
    {
        stream << ", effective sample size " << *ess;
    }
    printMoments(stream, indent.deeper(), particles.moments());
}
 
}  // namespace internal
 
}  // namespace Emrys
 
#endif  // EMRYS_FILTER_NO_IO