MixtureFilter.h

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <atomic>
#include <concepts>
#include <cstddef>
#include <expected>
#include <functional>
#include <numeric>
#include <optional>
#include <ranges>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Allocation.h>
#include <EmrysFilter/Core/Execution.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Print.h>
#include <EmrysFilter/Core/Weights.h>
#include <EmrysFilter/Distributions/Component.h>
#include <EmrysFilter/Distributions/Gaussian.h>
#include <EmrysFilter/Distributions/Mixture.h>
#include <EmrysFilter/Filters/FilterUpkeep.h>
#include <EmrysFilter/Filters/ImmFilter.h>
#include <EmrysFilter/Filters/KalmanFilter.h>
#include <EmrysFilter/Filters/Recursion.h>
#include <EmrysFilter/Models/Measurement.h>
#include <EmrysFilter/Models/ProcessModel.h>
#include <EmrysFilter/Models/Reduction.h>
#include <EmrysFilter/Models/Transform.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Space.h>
 
namespace Emrys
{

template <typename ComponentUpdate, typename Scalar = double, typename Allocator = DefaultAllocator>
struct MixtureFilterUpdate
{
    using ComponentUpdateType = ComponentUpdate;                          
    using ComponentsType = internal::Buffer<ComponentUpdate, Allocator>;  
 
    MixtureFilterUpdate() = default;  

    explicit MixtureFilterUpdate(Allocator const & allocator)
        : components(internal::ReBound<Allocator, ComponentUpdate>(allocator))
        , m_posteriorEnds(internal::ReBound<Allocator, Eigen::Index>(allocator))
    {}

    [[nodiscard]] Allocator
    get_allocator() const
    {
        return Allocator(components.get_allocator());
    }

    [[nodiscard]] Scalar
    logLikelihood() const
    {
        return m_logLikelihood;
    }

    [[nodiscard]] std::ranges::iota_view<Eigen::Index, Eigen::Index>
    posteriorIndices(Eigen::Index i) const
    {
        EMRYS_FILTER_ASSERT(i >= 0 && i < Eigen::Index(m_posteriorEnds.size()) && "prior component index out of range");
        Eigen::Index const first = i == 0 ? 0 : m_posteriorEnds[std::size_t(i - 1)];
        return std::views::iota(first, m_posteriorEnds[std::size_t(i)]);
    }
 
    ComponentsType components;  
    bool accepted = false;      
 
private:
 
    friend struct internal::UpdateAccess;
 
    Scalar m_logLikelihood = 0;  

    internal::Buffer<Eigen::Index, Allocator> m_posteriorEnds;
};

template <Space Space_, RealScalar Scalar_ = double, typename Allocator_ = DefaultAllocator>
struct GaussianMixturePredictedJoint
{
    using Scalar = Scalar_;                                        
    using SpaceType = Space_;                                      
    using JointType = PredictedJoint<Space_, Scalar_>;             
    using JointsType = internal::Buffer<JointType, Allocator_>;    
    using LogWeightsType = internal::Buffer<Scalar_, Allocator_>;  
    using LineageType = MixtureLineage<Scalar_, Allocator_>;       
 
    LogWeightsType logWeights;  
    JointsType joints;          
    LineageType lineage;        
};

template <
    Component Component_,
    typename Model_,
    typename Method_ = DefaultMethodFor<Model_>,
    typename Allocator_ = DefaultAllocator,
    RecursionOver<Component_> Recursion_ = internal::RecursionForType<Component_, Model_, Method_, Allocator_>,
    typename Upkeep_ = upkeep::None>
    requires std::same_as<typename Recursion_::ModelType, Model_>
          && std::same_as<typename Recursion_::MethodType, Method_>
class MixtureFilter
    : public internal::HeldUpkeep<MixtureFilter<Component_, Model_, Method_, Allocator_, Recursion_, Upkeep_>, Upkeep_>
{
    using UpkeepBase = internal::HeldUpkeep<MixtureFilter, Upkeep_>;  
 
    static_assert(
        internal::predictsBy<Model_, Method_>,
        "a filter with a fixed-step model predicts by the model's method"
    );
    static_assert(internal::checkPredictionMethod<Method_>());
 
public:
 
    using Scalar = typename Params<Component_>::Scalar;        
    using SpaceType = typename Params<Component_>::SpaceType;  
    using ModelType = Model_;                                  
    using MethodType = Method_;                                
    using StateType = MixtureParams<Component_, Allocator_>;   
    using ComponentType = Params<Component_>;                  
    using allocator_type = Allocator_;                         

    using RecursionType = Recursion_;

    using JointType = GaussianMixturePredictedJoint<SpaceType, Scalar, Allocator_>;

    template <Measurement M>
    using ResultType = MixtureFilterUpdate<
        internal::RebindUpdateAllocator<typename RecursionType::template UpdateType<M>::value_type, Allocator_>,
        Scalar,
        Allocator_>;

    template <Measurement M>
    using UpdateType = std::expected<ResultType<M>, NumericError>;

    MixtureFilter(
        Model_ model,
        StateType state,
        MethodFor<Model_, Method_> method = {},
        std::type_identity_t<Upkeep_> upkeep = {}
    )
        : UpkeepBase{std::move(upkeep)}
        , m_state{std::move(state)}
        , m_recursion{makeRecursion(std::move(model), std::move(method), m_state)}
        , m_posterior{m_state.space(), m_state.get_allocator()}
    {
        checkState();
    }

    MixtureFilter(Model_ model, StateType state, std::type_identity_t<Upkeep_> upkeep)
        requires(!std::same_as<Upkeep_, Emrys::upkeep::None>)
        : MixtureFilter{std::move(model), std::move(state), MethodFor<Model_, Method_>{}, std::move(upkeep)}
    {}

    MixtureFilter(Recursion_ recursion, StateType state, std::type_identity_t<Upkeep_> upkeep = {})
        : UpkeepBase{std::move(upkeep)}
        , m_state{std::move(state)}
        , m_recursion{std::move(recursion)}
        , m_posterior{m_state.space(), m_state.get_allocator()}
    {
        checkState();
    }

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

    [[nodiscard]] StateType const &
    state() const
    {
        return m_state;
    }

    [[nodiscard]] StateType &
    state()
    {
        return m_state;
    }

    void
    setState(StateType state)
    {
        EMRYS_FILTER_ASSERT(!state.empty() && "the mixture must have a component");
        m_state = std::move(state);
        if (m_posterior.get_allocator() != m_state.get_allocator() || m_posterior.space() != m_state.space()
            || m_posterior.capacity() != m_state.capacity())
        {
            // The assignment carried these in, so the scratch is made again from them.
            m_posterior = StateType{m_state.space(), m_state.get_allocator()};
        }
        if constexpr (requires { m_recursion.resetScratch(componentAllocator()); })
        {
            if (m_recursion.get_allocator() != componentAllocator())
            {
                m_recursion.resetScratch(componentAllocator());
            }
        }
        checkState();
    }

    [[nodiscard]] RecursionType const &
    recursion() const
    {
        return m_recursion;
    }

    [[nodiscard]] Model_ const &
    model() const
    {
        return m_recursion.model();
    }

    [[nodiscard]] Model_ &
    model()
    {
        return m_recursion.model();
    }

    [[nodiscard]] Method_ const &
    method() const
    {
        return m_recursion.method();
    }

    template <typename... Inputs>
        requires(!isFixedStep<Model_>)
    [[nodiscard]] auto
    predict(TimeStep<Scalar> auto const & dt, Inputs const &... inputs)
    {
        return predict(execution::sequential, dt, inputs...);
    }

    template <ExecutionPolicy P, typename... Inputs>
        requires(!isFixedStep<Model_>)
    [[nodiscard]] auto
    predict(P && policy, TimeStep<Scalar> auto const & dt, Inputs const &... inputs)
    {
        assertRunsUnder<P>();
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the mixture must have a component");
        Scalar const step = internal::toScalarTime<Scalar>(dt);
        return predictComponents(std::forward<P>(policy), [&](ComponentType & component) {
            return m_recursion.predict(component, step, inputs...);
        });
    }

    template <typename... Inputs>
        requires(isFixedStep<Model_> && !internal::leadsWithPolicy<Inputs...> && !internal::leadsWithSteps<Inputs...>)
    [[nodiscard]] auto
    predict(Inputs const &... inputs)
    {
        internal::rejectTimeStep<Inputs...>();
        return predict(execution::sequential, inputs...);
    }

    template <ExecutionPolicy P, typename... Inputs>
        requires(isFixedStep<Model_> && !internal::leadsWithSteps<Inputs...>)
    [[nodiscard]] auto
    predict(P && policy, Inputs const &... inputs)
    {
        internal::rejectTimeStep<Inputs...>();
        assertRunsUnder<P>();
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the mixture must have a component");
        return predictComponents(std::forward<P>(policy), [&](ComponentType & component) {
            return m_recursion.predict(component, inputs...);
        });
    }

    template <typename... Inputs>
        requires isFixedStep<Model_>
    [[nodiscard]] auto
    predict(Steps steps, Inputs const &... inputs)
    {
        return predict(execution::sequential, steps, inputs...);
    }

    template <ExecutionPolicy P, typename... Inputs>
        requires isFixedStep<Model_>
    [[nodiscard]] auto
    predict(P && policy, Steps steps, Inputs const &... inputs)
    {
        EMRYS_FILTER_ASSERT(steps.count >= 0 && "the number of steps must be nonnegative");
        using Result = decltype(predict(policy, inputs...));
        if constexpr (std::is_void_v<Result>)
        {
            for (Eigen::Index i = 0; i < steps.count; ++i)
            {
                predict(policy, inputs...);
            }
        }
        else
        {
            for (Eigen::Index i = 0; i < steps.count; ++i)
            {
                if (Result result = predict(policy, inputs...); !result)
                {
                    return result;
                }
            }
            return Result{};
        }
    }

    template <typename... Inputs>
        requires(!isFixedStep<Model_> && isGaussian<ComponentType>)
    [[nodiscard]] JointType
    predictJoint(TimeStep<Scalar> auto const & dt, Inputs const &... inputs)
    {
        auto joint = emptyJoint();
        predictJoint(execution::sequential, dt, joint, inputs...);
        return joint;
    }

    template <typename JointAllocator, typename... Inputs>
        requires(!isFixedStep<Model_> && isGaussian<ComponentType>)
    void
    predictJoint(
        TimeStep<Scalar> auto const & dt,
        GaussianMixturePredictedJoint<SpaceType, Scalar, JointAllocator> & joint,
        Inputs const &... inputs
    )
    {
        predictJoint(execution::sequential, dt, joint, inputs...);
    }

    template <ExecutionPolicy P, typename JointAllocator, typename... Inputs>
        requires(!isFixedStep<Model_> && isGaussian<ComponentType>)
    void
    predictJoint(
        P && policy,
        TimeStep<Scalar> auto const & dt,
        GaussianMixturePredictedJoint<SpaceType, Scalar, JointAllocator> & joint,
        Inputs const &... inputs
    )
    {
        Scalar const step = internal::toScalarTime<Scalar>(dt);
        recordJoint(std::forward<P>(policy), joint, [&](ComponentType & component) {
            return m_recursion.predictJoint(component, step, inputs...);
        });
    }

    template <typename... Inputs>
        requires(
            isFixedStep<Model_> && isGaussian<ComponentType> && !internal::leadsWithPolicy<Inputs...>
            && !internal::leadsWithSteps<Inputs...>
        )
    [[nodiscard]] JointType
    predictJoint(Inputs const &... inputs)
    {
        internal::rejectTimeStep<Inputs...>();
        auto joint = emptyJoint();
        predictJoint(execution::sequential, joint, inputs...);
        return joint;
    }

    template <typename JointAllocator, typename... Inputs>
        requires(isFixedStep<Model_> && isGaussian<ComponentType>)
    void
    predictJoint(GaussianMixturePredictedJoint<SpaceType, Scalar, JointAllocator> & joint, Inputs const &... inputs)
    {
        predictJoint(execution::sequential, joint, inputs...);
    }

    template <ExecutionPolicy P, typename JointAllocator, typename... Inputs>
        requires(isFixedStep<Model_> && isGaussian<ComponentType>)
    void
    predictJoint(
        P && policy,
        GaussianMixturePredictedJoint<SpaceType, Scalar, JointAllocator> & joint,
        Inputs const &... inputs
    )
    {
        recordJoint(std::forward<P>(policy), joint, [&](ComponentType & component) {
            return m_recursion.predictJoint(component, inputs...);
        });
    }

    template <Measurement M, typename Policy = reductions::Collapse>
        requires(
            !isGaussian<ComponentType>
            || MixtureReductionPolicy<Policy, internal::MeasuredUpdate<ComponentType, M>, Component_>
        )
    [[nodiscard]] UpdateType<M>
    update(M const & measurement, Policy const & policy = {})
    {
        return update(execution::sequential, measurement, policy);
    }

    template <ExecutionPolicy P, Measurement M, typename Policy = reductions::Collapse>
        requires(
            !isGaussian<ComponentType>
            || MixtureReductionPolicy<Policy, internal::MeasuredUpdate<ComponentType, M>, Component_>
        )
    [[nodiscard]] UpdateType<M>
    update(P && executionPolicy, M const & measurement, Policy const & policy = {})
    {
        auto result = emptyResult<ResultType<M>>();
        auto const updated = update(std::forward<P>(executionPolicy), measurement, policy, result);
        if (!updated)
        {
            return std::unexpected{updated.error()};
        }
        return result;
    }

    template <Measurement M, typename Policy, typename ResultAllocator>
        requires(
            !isGaussian<ComponentType>
            || MixtureReductionPolicy<Policy, internal::MeasuredUpdate<ComponentType, M>, Component_>
        )
    [[nodiscard]] std::expected<void, NumericError>
    update(
        M const & measurement,
        Policy const & policy,
        MixtureFilterUpdate<typename ResultType<M>::ComponentUpdateType, Scalar, ResultAllocator> & result
    )
    {
        return update(execution::sequential, measurement, policy, result);
    }

    template <ExecutionPolicy P, Measurement M, typename Policy, typename ResultAllocator>
        requires(
            !isGaussian<ComponentType>
            || MixtureReductionPolicy<Policy, internal::MeasuredUpdate<ComponentType, M>, Component_>
        )
    [[nodiscard]] std::expected<void, NumericError>
    update(
        P && executionPolicy,
        M const & measurement,
        Policy const & policy,
        MixtureFilterUpdate<typename ResultType<M>::ComponentUpdateType, Scalar, ResultAllocator> & result
    )
    {
        if constexpr (keepsBranches<M, Policy>())
        {
            return updateBranches(std::forward<P>(executionPolicy), measurement, policy, result);
        }
        else
        {
            // The reweighting reads each component's likelihood, so the recursion computes it
            // from the gain's factorization (decisions.md [innovation-statistics]).
            auto const reading = internal::WithStatistics<Policy>{policy};
            return updateComponents(
                std::forward<P>(executionPolicy),
                result,
                [&]<typename Update>(ComponentType & component, Update & slot) {
                    return m_recursion.update(component, measurement, reading, slot);
                },
                [&](ComponentType & component) {
                    return m_recursion.update(component, measurement, reading);
                }
            );
        }
    }
 
private:

    void
    checkState()
    {
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the initial mixture must have a component");
        internal::assertPartitions(m_recursion.model(), m_state.space());
        if (auto const capacity = m_state.capacity())
        {
            m_posterior.setCapacity(*capacity);
        }
        else
        {
            m_posterior.reserve(m_state.size());
        }
    }

    template <typename P>
    static constexpr void
    assertRunsUnder()
    {
        static_assert(
            internal::isSequential<P> || RecursionType::Reentrant,
            "a recursion with scratch, such as the IMM's, runs the components in order: pass Sequential"
        );
    }

    static RecursionType
    makeRecursion(Model_ model, MethodFor<Model_, Method_> method, StateType const & state)
    {
        if constexpr (std::is_constructible_v<
                          RecursionType,
                          Model_,
                          MethodFor<Model_, Method_>,
                          decltype(componentAllocator(state))>)
        {
            return RecursionType{std::move(model), std::move(method), componentAllocator(state)};
        }
        else
        {
            return RecursionType{std::move(model), std::move(method)};
        }
    }

    [[nodiscard]] static auto
    componentAllocator(StateType const & state)
    {
        if constexpr (requires(ComponentType const & c) { c.get_allocator(); })
        {
            EMRYS_FILTER_ASSERT(!state.empty() && "the initial mixture must have a component");
            return state.begin()->component.get_allocator();
        }
        else
        {
            return state.get_allocator();
        }
    }

    [[nodiscard]] auto
    componentAllocator() const
    {
        return componentAllocator(m_state);
    }

    template <Measurement M, typename Policy>
    [[nodiscard]] static consteval bool
    keepsBranches()
    {
        if constexpr (isGaussian<ComponentType>)
        {
            using Measured = typename decltype(internal::updateBranches<true>(
                std::declval<ComponentType const &>(),
                std::declval<M const &>()
            ))::value_type;
            return !internal::keepsOne<Policy, typename Measured::Container, ComponentType>;
        }
        else
        {
            return false;
        }
    }

    template <ExecutionPolicy P, typename PredictOne>
    [[nodiscard]] auto
    predictComponents(P && policy, PredictOne const & predictOne)
    {
        using Moved = std::invoke_result_t<PredictOne const &, ComponentType &>;
        if constexpr (std::is_void_v<Moved>)
        {
            m_state.forEach(std::forward<P>(policy), [&](auto & entry) {
                entry.updateComponent(predictOne);
            });
        }
        else
        {
            // The first error any component reports, written once under any policy.
            auto failed = std::atomic_flag{};
            auto error = typename Moved::error_type{};
            m_state.forEach(std::forward<P>(policy), [&](auto & entry) {
                Moved const moved = entry.updateComponent(predictOne);
                if (!moved && !failed.test_and_set())
                {
                    error = moved.error();
                }
            });
            return failed.test() ? Moved{std::unexpected{error}} : Moved{};
        }
    }

    [[nodiscard]] JointType
    emptyJoint() const
    {
        return {
            internal::bufferOf<Scalar>(get_allocator()),
            internal::bufferOf<PredictedJoint<SpaceType, Scalar>>(get_allocator()),
            typename JointType::LineageType(get_allocator())
        };
    }

    template <ExecutionPolicy P, typename JointAllocator, typename PredictOne>
    void
    recordJoint(
        P && policy,
        GaussianMixturePredictedJoint<SpaceType, Scalar, JointAllocator> & joint,
        PredictOne const & predictOne
    )
    {
        assertRunsUnder<P>();
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the mixture must have a component");
        std::size_t const count = std::size_t(m_state.size());
        internal::MixtureLineageAccess::assign(joint.lineage, m_state.lineage());
        joint.logWeights.clear();
        for (auto const & entry : m_state)
        {
            joint.logWeights.push_back(entry.logWeight);
        }
        joint.joints.resize(count);
        m_state.forEach(std::forward<P>(policy), [&](auto & entry) {
            joint.joints[std::size_t(entry.index())] = entry.updateComponent([&](ComponentType & component) {
                return predictOne(component);
            });
        });
        m_state.startLineage();
    }

    template <typename Result>
    [[nodiscard]] Result
    emptyResult() const
    {
        auto result = Result{get_allocator()};
        result.components.reserve(std::size_t(m_state.size()));
        internal::UpdateAccess::posteriorEnds(result).reserve(std::size_t(m_state.size()));
        return result;
    }

    template <typename Result>
    [[nodiscard]] auto &
    posteriorEndsFor(Result & result) const
    {
        auto & ends = internal::UpdateAccess::posteriorEnds(result);
        if (auto const capacity = m_state.capacity())
        {
            ends.reserve(std::size_t(*capacity));
        }
        ends.resize(std::size_t(m_state.size()));
        return ends;
    }

    template <typename Result>
    void
    keepPriorIndices(Result & result) const
    {
        auto & ends = internal::UpdateAccess::posteriorEnds(result);
        std::iota(ends.begin(), ends.end(), Eigen::Index{1});
    }

    template <ExecutionPolicy P, typename Result, typename UpdateInto, typename UpdateNew>
    [[nodiscard]] std::expected<void, NumericError>
    updateComponents(P && executionPolicy, Result & result, UpdateInto const & updateInto, UpdateNew const & updateNew)
    {
        assertRunsUnder<P>();
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the mixture must have a component");
        // Weights the caller left undefined fail the update before anything is written.
        auto const logTotal = m_state.logTotalWeight();
        if (!logTotal)
        {
            return std::unexpected{logTotal.error()};
        }
        Eigen::Index const count = m_state.size();
        using Update = typename Result::ComponentUpdateType;
 
        if constexpr (internal::isSequential<P>)
        {
            // The slots a result holds from the last update are reused, with the storage a
            // recursion's result keeps; the components beyond them get new ones, from the filter's
            // allocator.
            Eigen::Index const held = std::min(Eigen::Index(result.components.size()), count);
            result.components.erase(result.components.begin() + held, result.components.end());
            auto failure = std::optional<GaussianError>{};
            m_state.forEach(execution::sequential, [&](auto & entry) {
                if (failure)
                {
                    return;
                }
                auto const updated = [&]() -> std::expected<void, NumericError> {
                    if (entry.index() < held)
                    {
                        return entry.updateComponent([&](ComponentType & component) {
                            return updateInto(component, result.components[std::size_t(entry.index())]);
                        });
                    }
                    else if constexpr (internal::EmptyFromAllocator<Update, Allocator_>)
                    {
                        result.components.push_back(internal::emptyUpdate<Update>(get_allocator()));
                        return entry.updateComponent([&](ComponentType & component) {
                            return updateInto(component, result.components.back());
                        });
                    }
                    else
                    {
                        auto update = entry.updateComponent([&](ComponentType & component) {
                            return updateNew(component);
                        });
                        if (!update)
                        {
                            return std::unexpected{update.error()};
                        }
                        result.components.push_back(std::move(*update));
                        return {};
                    }
                }();
                if (!updated)
                {
                    failure = updated.error();
                }
            });
            if (failure)
            {
                return std::unexpected{*failure};
            }
        }
        else
        {
            using Updated = std::remove_cvref_t<decltype(updateNew(std::declval<ComponentType &>()))>;
            auto slots = internal::bufferOf<std::optional<Updated>>(get_allocator());
            slots.resize(std::size_t(count));
            m_state.forEach(std::forward<P>(executionPolicy), [&](auto & entry) {
                slots[std::size_t(entry.index())].emplace(entry.updateComponent([&](ComponentType & component) {
                    return updateNew(component);
                }));
            });
            result.components.clear();
            for (auto & slot : slots)
            {
                if (!*slot)
                {
                    return std::unexpected{slot->error()};
                }
                result.components.push_back(std::move(**slot));
            }
        }
 
        // Component i's weight times its likelihood, and that of the branches its policy kept.
        auto const logWeightOf = [&](Eigen::Index i) {
            return m_state[i].logWeight - *logTotal + result.components[std::size_t(i)].logLikelihood();
        };
        auto const logKeptOf = [&](Eigen::Index i) {
            return logWeightOf(i) + logKeptWeightOf(result.components[std::size_t(i)]);
        };
        // Whether the measurement is accepted is known before a weight is written: some component
        // accepted it and keeps a nonzero weight (decisions.md [update-outcomes]).
        auto evidence = LogSumExp<Scalar>{};
        bool accepted = false;
        for (Eigen::Index i = 0; i < count; ++i)
        {
            evidence.add(logWeightOf(i));
            accepted =
                accepted
                || (result.components[std::size_t(i)].accepted && logKeptOf(i) > -Eigen::NumTraits<Scalar>::infinity());
        }
        internal::UpdateAccess::logLikelihood(result) = evidence.value();
        result.accepted = accepted;
        auto & ends = posteriorEndsFor(result);
        if (!accepted)
        {
            // Only components of weight zero accepted it, if any, so the mixture is as it was.
            for (auto & component : result.components)
            {
                component.accepted = false;
            }
            keepPriorIndices(result);
            return {};
        }
        // The accepted components keep their order, so prior i ends where the count of those up
        // to it does. The scan starts from an Index: without a start, libc++ sums in the
        // components' own type.
        std::transform_inclusive_scan(
            result.components.begin(),
            result.components.end(),
            ends.begin(),
            std::plus<>{},
            [](auto const & component) {
                return Eigen::Index(component.accepted);
            },
            Eigen::Index{0}
        );
        for (Eigen::Index i = 0; i < count; ++i)
        {
            if (result.components[std::size_t(i)].accepted)
            {
                m_state.setLogWeight(i, logKeptOf(i));
            }
        }
        m_state.removeIf([&](Eigen::Index i, auto const &) {
            return !result.components[std::size_t(i)].accepted;
        });
        [[maybe_unused]] auto const normalized = m_state.normalize();
        EMRYS_FILTER_ASSERT(normalized && "an accepted component keeps a nonzero weight");
        return {};
    }

    template <ExecutionPolicy P, Measurement M, typename Policy, typename Result>
    [[nodiscard]] std::expected<void, NumericError>
    updateBranches(P && executionPolicy, M const & measurement, Policy const & policy, Result & result)
    {
        EMRYS_FILTER_ASSERT(!m_state.empty() && "the mixture must have a component");
        auto const logTotal = m_state.logTotalWeight();
        if (!logTotal)
        {
            return std::unexpected{logTotal.error()};
        }
        Eigen::Index const count = m_state.size();
        using Update = typename Result::ComponentUpdateType;
 
        m_posterior.clear();
        internal::MixtureLineageAccess::startFrom(m_posterior, m_state);
        auto & ends = posteriorEndsFor(result);
        auto evidence = LogSumExp<Scalar>{};
        // Reduces component i's measured branches into the posterior, each kept one weighted by the
        // likelihood of the branches it stands for; the evidence counts them all.
        // \returns whether any was kept
        auto const reduce = [&](Eigen::Index i, auto & measured) {
            Scalar const logWeight = m_state[i].logWeight - *logTotal + measured.logLikelihood();
            evidence.add(logWeight);
            auto out = PosteriorOutput{*this, i, logWeight};
            internal::reduceInto(policy, std::move(measured.branches), out);
            // The components are reduced in order, each into the end of the posterior.
            ends[std::size_t(i)] = m_posterior.size();
            return out.added;
        };
        if constexpr (internal::isSequential<P>)
        {
            Eigen::Index const held = std::min(Eigen::Index(result.components.size()), count);
            result.components.erase(result.components.begin() + held, result.components.end());
            for (Eigen::Index i = 0; i < count; ++i)
            {
                // The reweighting reads each component's likelihood, so the update computes it with its gain.
                auto const updated = internal::updateSlot<true>(
                    result.components,
                    std::size_t(i),
                    get_allocator(),
                    m_state[i].component,
                    measurement,
                    [&](auto & measured) {
                        return reduce(i, measured);
                    }
                );
                if (!updated)
                {
                    return updated;
                }
            }
        }
        else
        {
            using Measured =
                decltype(internal::updateBranches<true>(std::declval<ComponentType const &>(), measurement));
            auto measuredAll = internal::bufferOf<std::optional<Measured>>(get_allocator());
            measuredAll.resize(std::size_t(count));
            internal::forEachIndex(std::forward<P>(executionPolicy), count, [&](Eigen::Index i) {
                measuredAll[std::size_t(i)].emplace(internal::updateBranches<true>(m_state[i].component, measurement));
            });
            for (auto const & measured : measuredAll)
            {
                if (!*measured)
                {
                    return std::unexpected{measured->error()};
                }
            }
            result.components.clear();
            for (Eigen::Index i = 0; i < count; ++i)
            {
                auto & measured = **measuredAll[std::size_t(i)];
                bool const accepted = reduce(i, measured);
                if constexpr (internal::EmptyFromAllocator<Update, Allocator_>)
                {
                    result.components.push_back(internal::emptyUpdate<Update>(get_allocator()));
                    internal::filterUpdateInto(measured, accepted, result.components.back());
                }
                else
                {
                    result.components.push_back(internal::filterUpdate(std::move(measured), accepted));
                }
            }
        }
        internal::UpdateAccess::logLikelihood(result) = evidence.value();
        // The posterior is built apart, so a measurement that only components of weight zero accept
        // leaves the state as it was (decisions.md [update-outcomes]).
        result.accepted = !m_posterior.empty() && m_posterior.normalize().has_value();
        if (!result.accepted)
        {
            for (auto & component : result.components)
            {
                component.accepted = false;
            }
            keepPriorIndices(result);
            return {};
        }
        m_state.swap(m_posterior);
        return {};
    }

    struct PosteriorOutput
    {
        using ComponentType = Params<Component_>;   
        static constexpr bool HoldsSeveral = true;  
 
        MixtureFilter & filter;  
        Eigen::Index prior;      
        Scalar logWeight;        
        bool added = false;      

        void
        add(Scalar weight, ComponentType && component)
        {
            filter.m_posterior.add(logWeight + weight, std::move(component));
            internal::MixtureLineageAccess::inherit(filter.m_posterior, filter.m_state, prior);
            added = true;
        }

        void
        add(Scalar weight, ComponentType const & component)
        {
            add(weight, ComponentType(component));
        }
    };

    template <typename Update>
    [[nodiscard]] static Scalar
    logKeptWeightOf(Update const & update)
    {
        if constexpr (internal::isFilterUpdate<Update>)
        {
            return Scalar(internal::MixtureEvidenceAccess::logKeptWeight(update));
        }
        else
        {
            return Scalar(0);
        }
    }
 
    StateType m_state;          
    RecursionType m_recursion;  
    StateType m_posterior;      
};

template <Component C, typename Model, typename Method = DefaultMethodFor<Model>, typename Allocator>
    requires(!Recursion<Model> && !internal::UpkeepArgument<Method>)
MixtureFilter(Model, MixtureParams<C, Allocator>, Method = {}) -> MixtureFilter<
    C,
    internal::BoundStep<Model, typename Params<C>::SpaceType, typename Params<C>::Scalar>,
    Method,
    Allocator>;
 
template <Component C, RecursionOver<C> R, typename Allocator>
MixtureFilter(R, MixtureParams<C, Allocator>)
    -> MixtureFilter<C, typename R::ModelType, typename R::MethodType, Allocator, R>;
 
template <Component C, typename Model, typename Allocator, typename U>
    requires(!Recursion<Model> && internal::UpkeepArgument<U>)
MixtureFilter(Model, MixtureParams<C, Allocator>, U) -> MixtureFilter<
    C,
    internal::BoundStep<Model, typename Params<C>::SpaceType, typename Params<C>::Scalar>,
    DefaultMethodFor<Model>,
    Allocator,
    internal::RecursionForType<
        C,
        internal::BoundStep<Model, typename Params<C>::SpaceType, typename Params<C>::Scalar>,
        DefaultMethodFor<Model>,
        Allocator>,
    U>;
 
template <Component C, typename Model, typename Method, typename Allocator, typename U>
    requires(!Recursion<Model>)
MixtureFilter(Model, MixtureParams<C, Allocator>, Method, U) -> MixtureFilter<
    C,
    internal::BoundStep<Model, typename Params<C>::SpaceType, typename Params<C>::Scalar>,
    Method,
    Allocator,
    internal::RecursionForType<
        C,
        internal::BoundStep<Model, typename Params<C>::SpaceType, typename Params<C>::Scalar>,
        Method,
        Allocator>,
    U>;
 
template <Component C, RecursionOver<C> R, typename Allocator, typename U>
    requires internal::UpkeepArgument<U>
MixtureFilter(R, MixtureParams<C, Allocator>, U)
    -> MixtureFilter<C, typename R::ModelType, typename R::MethodType, Allocator, R, U>;


template <
    Space Space_,
    typename Model_,
    typename Method_ = DefaultMethodFor<Model_>,
    RealScalar Scalar_ = double,
    typename Allocator_ = DefaultAllocator,
    typename Upkeep_ = upkeep::None>
class GaussianMixtureFilter
    : public MixtureFilter<
          Gaussian<Space_, Scalar_>,
          Model_,
          Method_,
          Allocator_,
          internal::RecursionForType<Gaussian<Space_, Scalar_>, Model_, Method_, Allocator_>,
          Upkeep_>
{
public:

    using Base = MixtureFilter<
        Gaussian<Space_, Scalar_>,
        Model_,
        Method_,
        Allocator_,
        internal::RecursionForType<Gaussian<Space_, Scalar_>, Model_, Method_, Allocator_>,
        Upkeep_>;

    template <typename... Args>
        requires(
            std::constructible_from<Base, Args...>
            && !(sizeof...(Args) == 1 && (std::is_same_v<std::remove_cvref_t<Args>, GaussianMixtureFilter> && ...))
        )
    GaussianMixtureFilter(Args &&... args)
        : Base(std::forward<Args>(args)...)
    {}
};

// A guide with a defaulted method, its parameters in the class's order, and two with an upkeep.
template <
    typename Space,
    typename Model,
    typename Method = DefaultMethodFor<Model>,
    typename Scalar,
    typename Allocator>
    requires(!internal::UpkeepArgument<Method>)
GaussianMixtureFilter(Model, MixtureParams<Gaussian<Space, Scalar>, Allocator>, Method = {})
    -> GaussianMixtureFilter<Space, internal::BoundStep<Model, Space, Scalar>, Method, Scalar, Allocator>;
 
template <typename Space, typename Model, typename Scalar, typename Allocator, typename U>
    requires internal::UpkeepArgument<U>
GaussianMixtureFilter(Model, MixtureParams<Gaussian<Space, Scalar>, Allocator>, U) -> GaussianMixtureFilter<
    Space,
    internal::BoundStep<Model, Space, Scalar>,
    DefaultMethodFor<Model>,
    Scalar,
    Allocator,
    U>;
 
template <typename Space, typename Model, typename Method, typename Scalar, typename Allocator, typename U>
GaussianMixtureFilter(Model, MixtureParams<Gaussian<Space, Scalar>, Allocator>, Method, U)
    -> GaussianMixtureFilter<Space, internal::BoundStep<Model, Space, Scalar>, Method, Scalar, Allocator, U>;
 
}  // namespace Emrys
 
#ifndef EMRYS_FILTER_NO_IO
 
namespace Emrys
{
 
namespace internal
{

template <
    typename Component_,
    typename Model_,
    typename Method_,
    typename Allocator_,
    typename Recursion_,
    typename Upkeep_>
void
describe(
    std::ostream & stream,
    MixtureFilter<Component_, Model_, Method_, Allocator_, Recursion_, Upkeep_> const & filter,
    Indent indent
)
{
    describeFilter(stream, "mixture filter", filter.state(), indent);
}

template <typename Space_, typename Model_, typename Method_, typename Scalar_, typename Allocator_, typename Upkeep_>
void
describe(
    std::ostream & stream,
    GaussianMixtureFilter<Space_, Model_, Method_, Scalar_, Allocator_, Upkeep_> const & filter,
    Indent indent
)
{
    describeFilter(stream, "Gaussian mixture filter", filter.state(), indent);
}
 
}  // namespace internal
 
}  // namespace Emrys
 
#endif  // EMRYS_FILTER_NO_IO