Placement.h

Go to the documentation of this file.

#pragma once

 
#include <array>
#include <cstddef>
#include <optional>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Distributions/Moments.h>
#include <EmrysFilter/Models/MeasurementModel.h>
#include <EmrysFilter/Models/Methods.h>
#include <EmrysFilter/Models/Policies.h>
#include <EmrysFilter/Models/ProcessModel.h>
#include <EmrysFilter/Models/ResolvedReads.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Prefix, typename Name>
struct PlacedName : Name
{
    using Name::Name;
};

template <typename Prefix, typename Path>
struct JoinPathsOf;

template <typename... Prefix, typename... Tags>
struct JoinPathsOf<TagPath<Prefix...>, TagPath<Tags...>>
{
    using type = TagPath<Prefix..., Tags...>;
};


template <typename Prefix, typename Path>
using JoinPaths = typename JoinPathsOf<Prefix, Path>::type;

template <typename Path>
struct FirstTagOf;

template <typename First, typename... Rest>
struct FirstTagOf<TagPath<First, Rest...>>
{
    using type = First;
};


template <typename Path, typename S>
inline constexpr bool startsIn = countTag<typename FirstTagOf<Path>::type, S>() > 0;

template <typename T>
inline constexpr bool isRename = false;

template <PathValue X, typename Name>
inline constexpr bool isRename<OnPath<X, Name>> = isTagPath<PathOf<X>> || PathValue<X>;


template <typename Name, typename... Renames>
struct RenamedOf
{
    using type = Name;  
};

template <typename Name, typename X, typename New, typename... Rest>
struct RenamedOf<Name, OnPath<X, New>, Rest...>
{
    using type =
        std::conditional_t<std::is_same_v<PathOf<X>, TagPath<Name>>, New, typename RenamedOf<Name, Rest...>::type>;
};


template <typename Name, typename... Renames>
using Renamed = typename RenamedOf<Name, Renames...>::type;

template <typename Rename, typename S, typename Items>
struct RenamesAnInput
{
    static constexpr bool value = false;  
};

template <typename X, typename New, typename S, typename... Items>
struct RenamesAnInput<OnPath<X, New>, S, std::tuple<Items...>>
{
    static constexpr bool value =
        ((std::is_same_v<PathOf<X>, TagPath<typename Items::NameType>> && !startsIn<typename Items::PathType, S>) || ...
         || false);
};


template <typename A>
struct AffineNamesOf
{
    using type = std::tuple<>;  
};

template <typename... Names>
struct AffineNamesOf<AffineIn<Names...>>
{
    using type = std::tuple<Names...>;
};


template <typename Model>
using AffineNamesOfModel = typename AffineNamesOf<typename Model::template Policy<PolicyKind::affineIn>>::type;

template <typename Model, typename Name>
inline constexpr bool declaresAffine = []<typename... Names>(std::tuple<Names...> *) {
    return (std::is_same_v<Names, Name> || ...);
}(static_cast<AffineNamesOfModel<Model> *>(nullptr));

template <typename Name, typename... Names>
[[nodiscard]] consteval std::size_t
nameIndex()
{
    constexpr std::array<bool, sizeof...(Names) + 1> same{std::is_same_v<Names, Name>..., true};
    std::size_t i = 0;
    while (!same[i])
    {
        ++i;
    }
    return i;
}

template <typename Prefix, typename Model, typename... Renames>
struct Placement
{
    using PrefixSpace = PathTag<Prefix>;                                    
    using InnerReads = typename Model::template Policy<PolicyKind::reads>;  
    using InnerItems = typename ItemsOf<InnerReads>::type;                  
 
    static_assert(
        (isRename<Renames> && ...),
        "on<Path>(model, renames...) renames an input of the model with on<Name>(path<Tag>)"
    );
    static_assert(
        (RenamesAnInput<Renames, PrefixSpace, InnerItems>::value && ... && true),
        "a rename in on<Path>(model, on<Name>(path<Tag>)) names an input the model reads, a tag that is not a block of "
        "the block Path names"
    );

    template <typename Item>
    struct OuterItemOf
    {
        static constexpr bool State = startsIn<typename Item::PathType, PrefixSpace>;  
 
        using New = Renamed<typename Item::NameType, Renames...>;  
 
        static_assert(
            State || sameLeaf<New, typename Item::NameType>,
            "a renamed input has the same leaf beneath it as the name the model reads it under"
        );

        using type = std::conditional_t<
            State,
            ReadItem<JoinPaths<Prefix, typename Item::PathType>, PlacedName<Prefix, typename Item::NameType>>,
            ReadItem<TagPath<New>, New>>;
    };

    using ReadsPolicy = typename decltype([]<typename... Items>(std::tuple<Items...> *) {
        return std::type_identity<Reads<typename OuterItemOf<Items>::type...>>{};
    }(static_cast<InnerItems *>(nullptr)))::type;

    template <typename Name>
    using OuterName = typename decltype([]<typename... Items>(std::tuple<Items...> *) {
        return std::type_identity<std::tuple_element_t<
            nameIndex<Name, typename Items::NameType...>(),
            std::tuple<typename OuterItemOf<Items>::type::NameType...>>>{};
    }(static_cast<InnerItems *>(nullptr)))::type;

    using AffinePolicy = typename decltype([]<typename... Names>(std::tuple<Names...> *) {
        return std::type_identity<AffineIn<OuterName<Names>...>>{};
    }(static_cast<AffineNamesOfModel<Model> *>(nullptr)))::type;
};

template <typename Model>
class PlacedStep
{
public:
 
    using InnerIn = typename Model::template Policy<PolicyKind::reads>::InputSpace;  
    using InnerOut = OutputSpaceOf<Model>;                                           

    explicit PlacedStep(Model model)
        : m_model{std::move(model)}
    {}

    [[nodiscard]] Model const &
    model() const
    {
        return m_model;
    }

    [[nodiscard]] Model &
    model()
    {
        return m_model;
    }

    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] decltype(auto)
    operator()(In const & in, Out const & out, X const & x, TimeScalar dt) const
        requires MatrixExpression<X>
    {
        return callForm(m_model.function(), renamedInstance<InnerIn>(in), renamedInstance<InnerOut>(out), x, dt);
    }

    using Inner = std::remove_cvref_t<decltype(std::declval<Model const &>().function())>;

    static constexpr bool TakesBatches = takesBatches<Inner>;

    using MarkType = BatchMarkOrNone<Inner>;

    template <typename In, typename Out, typename X, typename TimeScalar, typename Y>
    void
    operator()(In const & in, Out const & out, X const & x, TimeScalar dt, Y & y) const
        requires TakesBatches
    {
        auto const innerIn = renamedInstance<InnerIn>(in);
        auto const innerOut = renamedInstance<InnerOut>(out);
        auto inner = y.template renamed<InnerOut>(innerIn);
        callForm(batchBody(m_model.function()), innerIn, innerOut, x, dt, inner);
    }

    template <typename Scalar, typename Out>
    [[nodiscard]] auto
    noiseOver(Out const & out, Scalar dt) const
    {
        return m_model.template noiseOver<Scalar>(renamedInstance<InnerOut>(out), dt);
    }
 
private:
 
    Model m_model;  
};

template <typename Prefix, typename Model>
using PlacedWrites = typename decltype([]<std::size_t... K>(std::index_sequence<K...>) {
    static_assert(
        (startsIn<WritePathAt<Model, K>, PathTag<Prefix>> && ...),
        "on<Path>(model) places a model whose writes() are blocks of the block Path names"
    );
    return std::type_identity<WritesAs<
        TagPaths<JoinPaths<Prefix, WritePathAt<Model, K>>...>,
        Composite<PlacedName<Prefix, WriteNameAt<Model, K>>...>>>{};
}(std::make_index_sequence<writeCount<Model>>{}))::type;

enum class EachKind
{
    state,     
    shared,    
    repeated,  
};

template <typename Prefix, typename Model, typename... Renames>
struct EachPlacement
{
    using RepeatedTag = PathTag<Prefix>;  
 
    static_assert(IsRepeated<RepeatedTag>, "onEach<Path>(model) places a model on a repeated block of the state");
 
    using Block = RepeatedBlock<RepeatedTag>;                               
    using InnerReads = typename Model::template Policy<PolicyKind::reads>;  
    using InnerItems = typename ItemsOf<InnerReads>::type;                  
    using InnerIn = typename InnerReads::InputSpace;                        
    using InnerOut = OutputSpaceOf<Model>;                                  
 
    static constexpr std::size_t ItemCount = std::tuple_size_v<InnerItems>;  
 
    static_assert(
        (isRename<Renames> && ...),
        "onEach<Path>(model, renames...) renames an input of the model with on<Name>(path<Tag>)"
    );
    static_assert(
        (RenamesAnInput<Renames, Block, InnerItems>::value && ... && true),
        "a rename in onEach<Path>(model, on<Name>(path<Tag>)) names an input the model reads, a tag that is not a "
        "block of one repetition of the block Path names"
    );
    static_assert(
        FixedSizeSpace<InnerIn> && FixedSizeSpace<InnerOut>,
        "onEach<Path>(model) repeats a model whose input and output are fixed-size"
    );
    static_assert(
        pathsPartition<Block, WritePathsOf<Model>>(),
        "onEach<Path>(model) places a model that writes every coordinate of one repetition of the block exactly once"
    );

    template <std::size_t K>
    using InnerItem = std::tuple_element_t<K, InnerItems>;

    template <std::size_t K>
    using NewName = Renamed<typename InnerItem<K>::NameType, Renames...>;

    template <std::size_t K>
    [[nodiscard]] static consteval EachKind
    kind()
    {
        if constexpr (startsIn<typename InnerItem<K>::PathType, Block>)
        {
            return EachKind::state;
        }
        else if constexpr (IsRepeated<NewName<K>>)
        {
            static_assert(
                sameLeaf<RepeatedBlock<NewName<K>>, typename InnerItem<K>::NameType>,
                "an input renamed to a Repeated tag repeats the leaf the model reads it under"
            );
            return EachKind::repeated;
        }
        else
        {
            static_assert(
                sameLeaf<NewName<K>, typename InnerItem<K>::NameType>,
                "a renamed input has the same leaf beneath it as the name the model reads it under"
            );
            return EachKind::shared;
        }
    }

    using ReadsPolicy = typename decltype([]<std::size_t... K>(std::index_sequence<K...>) {
        return std::type_identity<decltype(std::tuple_cat(
            std::tuple<ReadItem<Prefix, RepeatedTag>>{},
            std::conditional_t<
                kind<K>() == EachKind::state,
                std::tuple<>,
                std::tuple<ReadItem<TagPath<NewName<K>>, NewName<K>>>>{}...
        ))>{};
    }(std::make_index_sequence<ItemCount>{}))::type;

    using ReadsOfEach = typename decltype([]<typename... Items>(std::tuple<Items...> *) {
        return std::type_identity<Reads<Items...>>{};
    }(static_cast<ReadsPolicy *>(nullptr)))::type;

    using AffinePolicy = typename decltype([]<std::size_t... K>(std::index_sequence<K...>) {
        constexpr bool blockAffine =
            ((kind<K>() != EachKind::state || declaresAffine<Model, typename InnerItem<K>::NameType>) && ...);
        return std::type_identity<decltype(std::tuple_cat(
            std::conditional_t<blockAffine, std::tuple<RepeatedTag>, std::tuple<>>{},
            std::conditional_t < kind<K>() != EachKind::state && declaresAffine<Model, typename InnerItem<K>::NameType>,
            std::tuple<NewName<K>>,
            std::tuple < >> {}...
        ))>{};
    }(std::make_index_sequence<ItemCount>{}))::type;

    using AffineOfEach = typename decltype([]<typename... Names>(std::tuple<Names...> *) {
        return std::type_identity<AffineIn<Names...>>{};
    }(static_cast<AffinePolicy *>(nullptr)))::type;
};

template <typename Prefix, typename Model, typename... Renames>
class EachStep
{
    using P = EachPlacement<Prefix, Model, Renames...>;
    using RepeatedTag = typename P::RepeatedTag;
    using Block = typename P::Block;
    using InnerIn = typename P::InnerIn;
    using InnerOut = typename P::InnerOut;
 
    static constexpr int BlockSize = valueDimensionAtCompileTime<Block>;  
    static constexpr int BlockDifferenceSize = valueDimensionAtCompileTime<DifferenceSpace<Block>>;  
 
public:

    explicit EachStep(Model model)
        : m_model{std::move(model)}
    {}

    [[nodiscard]] Model const &
    model() const
    {
        return m_model;
    }

    [[nodiscard]] Model &
    model()
    {
        return m_model;
    }

    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] auto
    operator()(In const & in, Out const & out, X const & x, TimeScalar dt) const
        requires MatrixExpression<X>
    {
        using T = ScalarOf<X>;
        using InnerVector = typename Moments<InnerIn, T>::VectorType;
        constexpr bool ReturnsJacobian = returnsJacobianWith<Model, TimeScalar, TimeScalar>;
        Eigen::Index const count = pathBlock<RepeatedTag>(in).count();
        BlockRows const inRows = blockRows<RepeatedTag>(in);
        BlockRows const outRows = blockRows<RepeatedTag>(out);
        auto y = zeros<T>(out);
        auto jacobian = [&] {
            if constexpr (ReturnsJacobian)
            {
                using Matrix = MatrixOf<
                    TimeScalar,
                    valueDimensionAtCompileTime<DifferenceSpace<Out>>,
                    valueDimensionAtCompileTime<DifferenceSpace<In>>,
                    maxValueDimensionAtCompileTime<DifferenceSpace<Out>>,
                    maxValueDimensionAtCompileTime<DifferenceSpace<In>>>;
                return Matrix{Matrix::Zero(differenceDimension(out), differenceDimension(in))};
            }
            else
            {
                return NoPolicy{};
            }
        }();
        InnerIn const innerIn{};
        InnerOut const innerOut{};
        for (Eigen::Index i = 0; i < count; ++i)
        {
            InnerVector xi{InnerVector::Zero(valueDimension(innerIn))};
            [&]<std::size_t... K>(std::index_sequence<K...>) {
                (gather<K>(in, x, i, inRows, innerIn, xi), ...);
            }(std::make_index_sequence<P::ItemCount>{});
            auto result = vectorReturned(callForm(m_model.function(), innerIn, innerOut, xi, dt));
            auto const & value = valueOf(result);
            [&]<std::size_t... W>(std::index_sequence<W...>) {
                ((y.segment(
                      outRows.start + i * BlockSize + blockRows<WritePathAt<Model, W>>(Block{}).start,
                      blockRows<WritePathAt<Model, W>>(Block{}).size
                  ) =
                      value.segment(
                          blockRows<WriteNameAt<Model, W>>(innerOut).start,
                          blockRows<WriteNameAt<Model, W>>(innerOut).size
                      )),
                 ...);
            }(std::make_index_sequence<writeCount<Model>>{});
            if constexpr (ReturnsJacobian)
            {
                [&]<std::size_t... W>(std::index_sequence<W...>) {
                    (jacobianRows<W>(in, out, i, innerIn, innerOut, result.jacobian, jacobian), ...);
                }(std::make_index_sequence<writeCount<Model>>{});
            }
        }
        if constexpr (ReturnsJacobian)
        {
            return withJacobian(std::move(y), std::move(jacobian));
        }
        else
        {
            return y;
        }
    }

    using Inner = std::remove_cvref_t<decltype(std::declval<Model const &>().function())>;

    static constexpr bool TakesBatches = takesBatches<Inner>;

    using MarkType = NoJacobians;

    template <typename In, typename Out, typename X, typename TimeScalar, typename Y>
    void
    operator()(In const & in, Out const & out, X const & x, TimeScalar dt, Y & y) const
        requires TakesBatches
    {
        using T = typename X::Scalar;
        using InnerVector = typename Moments<InnerIn, T>::VectorType;
        using XI = MatrixOf<
            T,
            valueDimensionAtCompileTime<InnerIn>,
            X::ColsAtCompileTime,
            maxValueDimensionAtCompileTime<InnerIn>,
            X::MaxColsAtCompileTime>;
        using YI = MatrixOf<
            T,
            valueDimensionAtCompileTime<InnerOut>,
            X::ColsAtCompileTime,
            maxValueDimensionAtCompileTime<InnerOut>,
            X::MaxColsAtCompileTime>;
        Eigen::Index const repetitions = pathBlock<RepeatedTag>(in).count();
        Eigen::Index const points = x.cols();
        BlockRows const inRows = blockRows<RepeatedTag>(in);
        BlockRows const outRows = blockRows<RepeatedTag>(out);
        auto & values = y.values();
        values.setZero();
        InnerIn const innerIn{};
        InnerOut const innerOut{};
        auto xi = XI{};
        xi.resize(valueDimension(innerIn), points);
        auto yi = YI{};
        yi.resize(valueDimension(innerOut), points);
        auto column = InnerVector{};
        column.resize(valueDimension(innerIn));
        for (Eigen::Index i = 0; i < repetitions; ++i)
        {
            for (Eigen::Index j = 0; j < points; ++j)
            {
                column.setZero();
                [&]<std::size_t... K>(std::index_sequence<K...>) {
                    (gather<K>(in, x.col(j), i, inRows, innerIn, column), ...);
                }(std::make_index_sequence<P::ItemCount>{});
                xi.col(j) = column;
            }
            callBatched(m_model.function(), innerIn, innerOut, xi, yi, dt);
            [&]<std::size_t... W>(std::index_sequence<W...>) {
                ((values.middleRows(
                      outRows.start + i * BlockSize + blockRows<WritePathAt<Model, W>>(Block{}).start,
                      blockRows<WritePathAt<Model, W>>(Block{}).size
                  ) =
                      yi.middleRows(
                          blockRows<WriteNameAt<Model, W>>(innerOut).start,
                          blockRows<WriteNameAt<Model, W>>(innerOut).size
                      )),
                 ...);
            }(std::make_index_sequence<writeCount<Model>>{});
        }
    }

    template <typename Scalar, typename Out>
    [[nodiscard]] typename Moments<Out, Scalar>::MatrixType
    noiseOver(Out const & out, Scalar dt) const
    {
        using Matrix = typename Moments<Out, Scalar>::MatrixType;
        Eigen::Index const n = differenceDimension(out);
        Matrix q = Matrix::Zero(n, n);
        InnerOut const innerOut{};
        auto const inner = m_model.template noiseOver<Scalar>(innerOut, dt);
        Eigen::Index const count = pathBlock<RepeatedTag>(out).count();
        Eigen::Index const start = differenceBlockRows<RepeatedTag>(out).start;
        for (Eigen::Index i = 0; i < count; ++i)
        {
            [&]<std::size_t... W>(std::index_sequence<W...>) {
                (noiseRows<W>(q, start + i * BlockDifferenceSize, innerOut, inner), ...);
            }(std::make_index_sequence<writeCount<Model>>{});
        }
        return q;
    }
 
private:

    template <std::size_t K, typename In, typename X, typename XI>
    static void
    gather(In const & in, X const & x, Eigen::Index i, BlockRows const & inRows, InnerIn const & innerIn, XI & xi)
    {
        using Item = typename P::template InnerItem<K>;
        BlockRows const rows = blockRows<typename Item::NameType>(innerIn);
        if constexpr (P::template kind<K>() == EachKind::state)
        {
            BlockRows const within = blockRows<typename Item::PathType>(Block{});
            xi.segment(rows.start, rows.size) = x.segment(inRows.start + i * BlockSize + within.start, within.size);
        }
        else if constexpr (P::template kind<K>() == EachKind::repeated)
        {
            xi.segment(rows.start, rows.size) =
                block<typename P::template NewName<K>>(in, x).segment(i * rows.size, rows.size);
        }
        else
        {
            xi.segment(rows.start, rows.size) = block<typename P::template NewName<K>>(in, x);
        }
    }

    template <std::size_t W, typename In, typename Out, typename J, typename Whole>
    static void
    jacobianRows(
        In const & in,
        Out const & out,
        Eigen::Index i,
        InnerIn const & innerIn,
        InnerOut const & innerOut,
        J const & inner,
        Whole & jacobian
    )
    {
        BlockRows const innerRows = differenceBlockRows<WriteNameAt<Model, W>>(innerOut);
        Eigen::Index const row = differenceBlockRows<RepeatedTag>(out).start + i * BlockDifferenceSize
                               + differenceBlockRows<WritePathAt<Model, W>>(Block{}).start;
        [&]<std::size_t... K>(std::index_sequence<K...>) {
            (
                [&] {
                    using Item = typename P::template InnerItem<K>;
                    BlockRows const innerColumns = differenceBlockRows<typename Item::NameType>(innerIn);
                    Eigen::Index column = 0;
                    if constexpr (P::template kind<K>() == EachKind::state)
                    {
                        column = differenceBlockRows<RepeatedTag>(in).start + i * BlockDifferenceSize
                               + differenceBlockRows<typename Item::PathType>(Block{}).start;
                    }
                    else if constexpr (P::template kind<K>() == EachKind::repeated)
                    {
                        column = differenceBlockRows<typename P::template NewName<K>>(in).start + i * innerColumns.size;
                    }
                    else
                    {
                        column = differenceBlockRows<typename P::template NewName<K>>(in).start;
                    }
                    jacobian.block(row, column, innerRows.size, innerColumns.size) =
                        inner.block(innerRows.start, innerColumns.start, innerRows.size, innerColumns.size);
                }(),
                ...
            );
        }(std::make_index_sequence<P::ItemCount>{});
    }

    template <std::size_t W, typename Q, typename Inner>
    static void
    noiseRows(Q & q, Eigen::Index start, InnerOut const & innerOut, Inner const & inner)
    {
        BlockRows const innerRows = differenceBlockRows<WriteNameAt<Model, W>>(innerOut);
        Eigen::Index const row = start + differenceBlockRows<WritePathAt<Model, W>>(Block{}).start;
        [&]<std::size_t... V>(std::index_sequence<V...>) {
            (
                [&] {
                    BlockRows const innerColumns = differenceBlockRows<WriteNameAt<Model, V>>(innerOut);
                    Eigen::Index const column = start + differenceBlockRows<WritePathAt<Model, V>>(Block{}).start;
                    q.block(row, column, innerRows.size, innerColumns.size) =
                        inner.block(innerRows.start, innerColumns.start, innerRows.size, innerColumns.size);
                }(),
                ...
            );
        }(std::make_index_sequence<writeCount<Model>>{});
    }
 
    Model m_model;  
};

template <typename F, typename InnerIn>
class PlacedMeasure
{
public:

    explicit PlacedMeasure(F function)
        : m_function{std::move(function)}
    {}

    template <typename In, typename Out, typename X>
    [[nodiscard]] decltype(auto)
    operator()(In const & in, Out const & out, X const & x) const
        requires MatrixExpression<X>
    {
        return callForm(m_function, renamedInstance<InnerIn>(in), out, x);
    }

    static constexpr bool TakesBatches = takesBatches<F>;

    using MarkType = BatchMarkOrNone<F>;

    template <typename In, typename Out, typename X, typename Y>
    void
    operator()(In const & in, Out const & out, X const & x, Y & y) const
        requires TakesBatches
    {
        auto const innerIn = renamedInstance<InnerIn>(in);
        auto inner = y.template renamed<Out>(innerIn);
        callForm(batchBody(m_function), innerIn, out, x, inner);
    }

    [[nodiscard]] Eigen::Index
    outputRows() const
        requires requires(F const & f) { f.outputRows(); }
    {
        return m_function.outputRows();
    }
 
private:
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS F m_function;  
};

template <typename G, typename InnerIn>
struct PlacedLikelihood
{
    G function;  

    template <typename In, typename X, typename Z>
    [[nodiscard]] auto
    operator()(In const & in, X const & x, Z const & z) const
    {
        return function(renamedInstance<InnerIn>(in), x, z);
    }
};

template <typename G, typename InnerIn>
struct PlacedSample
{
    G function;  

    template <typename In, typename X, typename R>
    [[nodiscard]] auto
    operator()(In const & in, X const & x, R & rng) const
    {
        return function(renamedInstance<InnerIn>(in), x, rng);
    }
};

template <typename Prefix, typename G>
struct PlacedOutputSpace
{
    G function;  

    template <typename StateSpace>
    [[nodiscard]] auto
    operator()(StateSpace const & space) const
    {
        return function(pathSubspace<Prefix>(space));
    }
};

template <typename Prefix, typename InnerIn, typename ReadsPolicy, typename AffinePolicy, typename Policy>
[[nodiscard]] auto
placedPolicy(Policy const & policy)
{
    constexpr std::optional<PolicyKind> kind = policyKind<Policy>;
    if constexpr (kind == PolicyKind::reads)
    {
        return ReadsPolicy{};
    }
    else if constexpr (kind == PolicyKind::affineIn)
    {
        return AffinePolicy{};
    }
    else if constexpr (kind == PolicyKind::logLikelihood)
    {
        using G = std::remove_cvref_t<decltype(policy.function)>;
        return LogLikelihood<PlacedLikelihood<G, InnerIn>>{PlacedLikelihood<G, InnerIn>{policy.function}};
    }
    else if constexpr (kind == PolicyKind::sample)
    {
        using G = std::remove_cvref_t<decltype(policy.function)>;
        return Sample<PlacedSample<G, InnerIn>>{PlacedSample<G, InnerIn>{policy.function}};
    }
    else if constexpr (kind == PolicyKind::outputSpace)
    {
        if constexpr (Policy::Computed)
        {
            using G = std::remove_cvref_t<decltype(policy.value)>;
            return OutputSpace<PlacedOutputSpace<Prefix, G>>{PlacedOutputSpace<Prefix, G>{policy.value}};
        }
        else
        {
            return policy;
        }
    }
    else
    {
        return policy;
    }
}

template <typename... Models>
class MeasurementParts
{
public:

    explicit MeasurementParts(Models... models)
        : m_parts{std::move(models)...}
    {}

    [[nodiscard]] std::tuple<Models...> const &
    parts() const
    {
        return m_parts;
    }

    template <typename In, typename Out, typename X>
    [[nodiscard]] auto
    operator()(In const & in, Out const & out, X const & x) const
        requires MatrixExpression<X>
    {
        return stepParts<ScalarOf<X>>(m_parts, in, out, x);
    }

    static constexpr bool TakesBatches = anyTakesBatches<Models...>;

    using MarkType = NoJacobians;

    template <typename In, typename Out, typename X, typename Y>
    void
    operator()(In const & in, Out const & out, X const & x, Y & y) const
        requires TakesBatches
    {
        stepPartsBatch(m_parts, in, out, x, y);
    }
 
private:
 
    std::tuple<Models...> m_parts;  
};

template <typename Out, typename... Models>
struct SumOfLikelihoods
{
    std::tuple<Models...> parts;  

    template <typename In, typename X, typename Z>
    [[nodiscard]] auto
    operator()(In const & in, X const & x, Z const & z) const
    {
        return [&]<std::size_t... I>(std::index_sequence<I...>) {
            return (likelihoodOf<I>(in, x, z) + ...);
        }(std::index_sequence_for<Models...>{});
    }
 
private:

    template <std::size_t I, typename In, typename X, typename Z>
    [[nodiscard]] auto
    likelihoodOf(In const & in, X const & x, Z const & z) const
    {
        using Part = std::tuple_element_t<I, std::tuple<Models...>>;
        using PartOut = OutputSpaceOf<Part>;
        auto const pin = partIn<Part>(in);
        auto const px = partVector<Part>(in, pin, x);
        Out const out{};
        PartOut const pout{};
        auto pz = zeros<typename Z::Scalar>(pout);
        [&]<typename... Names>(std::tuple<Names...> *) {
            ((block<Names>(pout, pz) = block<Names>(out, z)), ...);
        }(static_cast<CompositeBlocks<PartOut> *>(nullptr));
        return callWithScalar<typename X::Scalar>(
            std::get<I>(parts).template policy<PolicyKind::logLikelihood>().function,
            pin,
            px,
            pz
        );
    }
};

template <typename Model>
[[nodiscard]] auto
placedNoiseParts(Model const & model)
{
    if constexpr (!Model::template has<PolicyKind::additiveNoise>)
    {
        return std::tuple<>{};
    }
    else
    {
        using Noise = typename Model::template Policy<PolicyKind::additiveNoise>;
        if constexpr (Noise::Placed)
        {
            return model.template policy<PolicyKind::additiveNoise>().parts();
        }
        else
        {
            static_assert(
                writeCount<Model> == 1,
                "a part of a composition whose additiveNoise() covers its whole output writes one block; place each "
                "part of its noise with on<Path>() to write several"
            );
            return std::tuple{on<WritePathAt<Model, 0>>(
                typename Noise::Distribution{model.template policy<PolicyKind::additiveNoise>().distribution()}
            )};
        }
    }
}
 
}  // namespace internal

template <typename Path, ProcessModel Model, typename... Renames>
[[nodiscard]] auto
on(Model const & model, Renames const &...)
{
    static_assert(internal::wholeStateRefused<Model>);
    using Prefix = internal::AsTagPath<Path>;
    using Placement = internal::Placement<Prefix, Model, Renames...>;
    using Step = internal::PlacedStep<Model>;
    return ProcessFunction<
        Step,
        typename Placement::ReadsPolicy,
        internal::PlacedWrites<Prefix, Model>,
        typename Placement::AffinePolicy>{
        Step{model},
        typename Placement::ReadsPolicy{},
        internal::PlacedWrites<Prefix, Model>{},
        typename Placement::AffinePolicy{}
    };
}

template <typename Path, MeasurementModel Model, typename... Renames>
[[nodiscard]] auto
on(Model const & model, Renames const &...)
{
    static_assert(internal::wholeStateRefused<Model>);
    using Prefix = internal::AsTagPath<Path>;
    using Placement = internal::Placement<Prefix, Model, Renames...>;
    using InnerIn = typename Model::template Policy<internal::PolicyKind::reads>::InputSpace;
    using ReadsPolicy = typename Placement::ReadsPolicy;
    using AffinePolicy = typename Placement::AffinePolicy;
    auto const function = [&] {
        if constexpr (Model::HasFunction)
        {
            return internal::PlacedMeasure<typename Model::FunctionType, InnerIn>{model.function()};
        }
        else
        {
            return internal::NoFunction{};
        }
    }();
    return std::apply(
        [&](auto const &... policies) {
            return MeasurementFunction<
                std::remove_cvref_t<decltype(function)>,
                decltype(internal::placedPolicy<Prefix, InnerIn, ReadsPolicy, AffinePolicy>(policies))...>{
                function,
                internal::placedPolicy<Prefix, InnerIn, ReadsPolicy, AffinePolicy>(policies)...
            };
        },
        model.policies()
    );
}

template <typename Path, ProcessModel Model, typename... Renames>
[[nodiscard]] auto
onEach(Model const & model, Renames const &...)
{
    static_assert(internal::wholeStateRefused<Model>);
    using Prefix = internal::AsTagPath<Path>;
    using Each = internal::EachPlacement<Prefix, Model, Renames...>;
    using Step = internal::EachStep<Prefix, Model, Renames...>;
    return ProcessFunction<Step, typename Each::ReadsOfEach, Writes<Prefix>, typename Each::AffineOfEach>{
        Step{model},
        typename Each::ReadsOfEach{},
        Writes<Prefix>{},
        typename Each::AffineOfEach{}
    };
}
 
namespace internal
{

template <typename... Models>
[[nodiscard]] auto
composeMeasurements(Models... models)
{
    static_assert((wholeStateRefused<Models> && ...));
    // A whole-state model has no writes() either, which its own refusal says.
    static_assert(
        ((Models::template has<PolicyKind::writes> || readsWholeState<Models>) && ...),
        "compose() takes measurement models that name their output with writes()"
    );
    static_assert(
        (Models::HasFunction && ...),
        "compose() takes measurement models with a function; one known by its likelihood alone has no output"
    );
    static_assert(
        (!Models::template has<PolicyKind::outputSpace> && ...),
        "compose() takes measurement models whose output writes() fixes, without outputSpace()"
    );
    constexpr std::size_t likelihoods = (std::size_t{Models::template has<PolicyKind::logLikelihood>} + ...);
    static_assert(
        likelihoods == 0 || likelihoods == sizeof...(Models),
        "in a composition of measurement models, every model gives logLikelihood() or none does"
    );
    using AllItems = decltype(std::tuple_cat(std::declval<ReadItemsOf<Models>>()...));
    using ReadsPolicy = typename ReadUnionOf<AllItems>::type;
    using WritesPolicy = typename WritesOfTuple<
        decltype(std::tuple_cat(std::declval<typename PathsOf<WritePathsOf<Models>>::type>()...)),
        decltype(std::tuple_cat(std::declval<CompositeBlocks<OutputSpaceOf<Models>>>()...))>::type;
    using AffinePolicy = typename AffineOf<ReadsPolicy, std::tuple<Models...>>::type;
    using Function = MeasurementParts<Models...>;
    auto const noise = std::tuple_cat(placedNoiseParts(models)...);
    auto const base = std::tuple{ReadsPolicy{}, WritesPolicy{}, AffinePolicy{}};
    auto const withNoise = [&] {
        if constexpr (std::tuple_size_v<std::remove_cvref_t<decltype(noise)>> > 0)
        {
            return std::tuple_cat(
                base,
                std::tuple{std::apply(
                    [](auto const &... parts) {
                        return AdditiveNoise<std::remove_cvref_t<decltype(parts)>...>{parts...};
                    },
                    noise
                )}
            );
        }
        else
        {
            return base;
        }
    }();
    auto const policies = [&] {
        if constexpr (likelihoods > 0)
        {
            using Sum = SumOfLikelihoods<typename WritesPolicy::OutputSpace, Models...>;
            return std::tuple_cat(withNoise, std::tuple{LogLikelihood<Sum>{Sum{std::tuple{models...}}}});
        }
        else
        {
            return withNoise;
        }
    }();
    return std::apply(
        [&](auto const &... each) {
            return MeasurementFunction<Function, std::remove_cvref_t<decltype(each)>...>{
                Function{std::move(models)...},
                each...
            };
        },
        policies
    );
}
 
}  // namespace internal

template <MeasurementModel... Models>
    requires(sizeof...(Models) > 0)
[[nodiscard]] auto
compose(Models &&... models)
{
    return internal::composeMeasurements(std::remove_cvref_t<Models>{std::forward<Models>(models)}...);
}
 
}  // namespace Emrys