ProcessModel.h

Go to the documentation of this file.

#pragma once

 
#include <array>
#include <cstddef>
#include <optional>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <Eigen/LU>
 
#include <EmrysFilter/Core/Forwarding.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Multiply.h>
#include <EmrysFilter/Core/Solve.h>
#include <EmrysFilter/Core/Time.h>
#include <EmrysFilter/Distributions/Gaussian.h>
#include <EmrysFilter/Distributions/Moments.h>
#include <EmrysFilter/Distributions/Sampling.h>
#include <EmrysFilter/Models/Batched.h>
#include <EmrysFilter/Models/MeasurementModel.h>
#include <EmrysFilter/Models/Methods.h>
#include <EmrysFilter/Models/Policies.h>
#include <EmrysFilter/Models/ResolvedReads.h>
#include <EmrysFilter/Models/Transform.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Space>
class RowCoverage
{
public:

    explicit RowCoverage(Space const & space)
        : m_counts{CountsType::Zero(valueDimension(space))}
    {}

    void
    push_back(BlockRows const & rows)
    {
        if (rows.start < 0 || rows.size < 0 || rows.start + rows.size > m_counts.size())
        {
            m_inRange = false;
        }
        else
        {
            m_counts.segment(rows.start, rows.size) += 1;
        }
    }

    [[nodiscard]] bool
    partitioned() const
    {
        return m_inRange && (m_counts == 1).all();
    }
 
private:

    using CountsType = Eigen::
        Array<int, valueDimensionAtCompileTime<Space>, 1, Eigen::ColMajor, maxValueDimensionAtCompileTime<Space>, 1>;
 
    CountsType m_counts;    
    bool m_inRange = true;  
};

enum class StaticPartition
{
    unknown,  
    holds,    
    fails     
};

template <typename Space, typename... Lists>
consteval bool
pathsPartition()
{
    auto counts = std::array<int, std::size_t(Space::FixedValueDimension)>{};
    [[maybe_unused]] auto const count = [&]<typename... Paths>(TagPaths<Paths...> *) {  // unused for no Lists
        (
            [&] {
                constexpr ResolvedIndex rows = ResolvePath<Space, Paths>::value;
                for (int i = rows.offset; i < rows.offset + rows.size; ++i)
                {
                    ++counts[std::size_t(i)];
                }
            }(),
            ...
        );
    };
    (count(static_cast<Lists *>(nullptr)), ...);
    for (int const owners : counts)
    {
        if (owners != 1)
        {
            return false;
        }
    }
    return true;
}

template <typename Model, typename Space>
consteval StaticPartition
staticPartition()
{
    if constexpr (requires { Model::template staticPartition<Space>(); })
    {
        return Model::template staticPartition<Space>();
    }
    else
    {
        return StaticPartition::unknown;
    }
}

template <typename Model, typename Space>
void
assertPartitions(
    [[maybe_unused]] Model const & model,
    [[maybe_unused]] Space const & space,
    [[maybe_unused]] char const * message = "the process model must write every coordinate of the state space once"
)
{
    constexpr StaticPartition known = staticPartition<Model, Space>();
    static_assert(
        known != StaticPartition::fails,
        "the process model must write every coordinate of the state space exactly once: a block it writes twice, or "
        "one no model writes, is not a partition"
    );
    if constexpr (known == StaticPartition::unknown)
    {
        EMRYS_FILTER_ASSERT(model.partitions(space) && message);
    }
}

template <typename Space, typename Paths>
consteval StaticPartition
staticPartitionOfWrites()
{
    if constexpr (Space::ValueDimensionStep == 0)
    {
        return pathsPartition<Space, Paths>() ? StaticPartition::holds : StaticPartition::fails;
    }
    else
    {
        return StaticPartition::unknown;
    }
}

template <typename Space, typename... Paths>
[[nodiscard]] bool
writesPartition(Space const & space, TagPaths<Paths...> *)
{
    auto coverage = RowCoverage<Space>{space};
    (coverage.push_back(blockRows<Paths>(space)), ...);
    return coverage.partitioned();
}

template <typename P>
struct PathsOf;

template <typename... Paths>
struct PathsOf<TagPaths<Paths...>>
{
    using type = std::tuple<Paths...>;
};


template <typename Model>
inline constexpr bool writesWholeState = !Model::template has<PolicyKind::writes>;

template <typename Model>
using WritePathsOf = typename Model::template Policy<PolicyKind::writes>::PathList;

template <typename Model>
using OutputSpaceOf = typename Model::template Policy<PolicyKind::writes>::OutputSpace;

template <typename Model, std::size_t K>
using WritePathAt = std::tuple_element_t<K, typename PathsOf<WritePathsOf<Model>>::type>;

template <typename Model, std::size_t K>
using WriteNameAt = std::tuple_element_t<K, CompositeBlocks<OutputSpaceOf<Model>>>;

template <typename Model>
inline constexpr std::size_t writeCount = std::tuple_size_v<typename PathsOf<WritePathsOf<Model>>::type>;

template <typename Model, typename Space>
using OutputSpaceIn = typename decltype([] {
    if constexpr (writesWholeState<Model>)
    {
        return std::type_identity<Space>{};
    }
    else
    {
        return std::type_identity<OutputSpaceOf<Model>>{};
    }
}())::type;

template <typename Model, typename Space>
[[nodiscard]] OutputSpaceIn<Model, Space>
outputOf(Space const & space)
{
    if constexpr (writesWholeState<Model>)
    {
        return space;
    }
    else
    {
        return [&]<std::size_t... K>(std::index_sequence<K...>) {
            return OutputSpaceOf<Model>{nameOf<WriteNameAt<Model, K>>(pathBlock<WritePathAt<Model, K>>(space))...};
        }(std::make_index_sequence<writeCount<Model>>{});
    }
}

template <typename Target, typename Source>
[[nodiscard]] Target
renamedInstance(Source const & source)
{
    if constexpr (std::is_same_v<Target, Source>)
    {
        return source;
    }
    else
    {
        return []<typename... Names, typename... Blocks>(
                   Source const & from,
                   std::tuple<Names...> *,
                   std::tuple<Blocks...> *
               ) {
            return Target{nameOf<Names>(pathBlock<Blocks>(from))...};
        }(source, static_cast<CompositeBlocks<Target> *>(nullptr), static_cast<CompositeBlocks<Source> *>(nullptr));
    }
}

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

template <typename X, typename Path>
inline constexpr bool isOnPath<OnPath<X, Path>>{true};


template <typename Matrix, typename Value, typename Scalar>
void
addSpectralDensity(Matrix & q, BlockRows const & rows, Value const & value, Scalar dt)
{
    if constexpr (MatrixExpression<Value>)
    {
        EMRYS_FILTER_ASSERT(
            value.rows() == rows.size && value.cols() == rows.size
            && "a spectral density matrix is over the block's differences"
        );
        q.block(rows.start, rows.start, rows.size, rows.size) += dt * value.template cast<Scalar>();
    }
    else
    {
        q.block(rows.start, rows.start, rows.size, rows.size).diagonal().array() += dt * Scalar(value);
    }
}

template <typename Out, typename Matrix, typename Part, typename Scalar>
void
addContinuousPart(Out const & out, Matrix & q, Part const & part, Scalar dt)
{
    if constexpr (isOnPath<Part>)
    {
        [&]<typename X, typename Path>(OnPath<X, Path> const *) {
            addSpectralDensity(q, differenceBlockRows<Path>(out), part.value, dt);
        }(static_cast<Part const *>(nullptr));
    }
    else if constexpr (isPlacedByPath<Part>)
    {
        static_assert(
            std::is_void_v<Part>,
            "continuousNoise() takes no onEach<Path>(q); place the model on each repetition with onEach<Path>(model), "
            "and "
            "give it continuousNoise(q) for one repetition"
        );
    }
    else
    {
        addSpectralDensity(q, BlockRows{0, differenceDimension(out)}, part, dt);
    }
}

template <typename Set>
[[nodiscard]] consteval bool
checkProcessPolicies()
{
    static_assert(
        Set::template has<PolicyKind::reads> || !Set::template has<PolicyKind::writes>,
        "a process model that writes blocks with writes() says what it reads with reads(); give neither to read and "
        "write the whole state"
    );
    static_assert(
        Set::template has<PolicyKind::writes> || !Set::template has<PolicyKind::reads>,
        "a process model that reads blocks with reads() says which blocks of the state it writes with writes(); give "
        "neither to read and write the whole state"
    );
    static_assert(
        Set::template has<PolicyKind::reads> || !Set::template has<PolicyKind::affineIn>,
        "affineIn() names items of reads(), and a model without reads() reads the whole state under no names"
    );
    static_assert(
        !Set::template has<PolicyKind::logLikelihood>,
        "logLikelihood(g) is a measurement model's; a process model has no observation to weight points by"
    );
    static_assert(
        !Set::template has<PolicyKind::sample>,
        "sample(g) is a measurement model's; a process model's draw is its step and its noise"
    );
    static_assert(
        !Set::template has<PolicyKind::outputSpace>,
        "a process model's output is the blocks of the state it writes, whose spaces the state has, so it takes no "
        "outputSpace()"
    );
    static_assert(
        !(Set::template has<PolicyKind::continuousNoise> && Set::template has<PolicyKind::additiveNoise>),
        "a process model's noise is continuousNoise(), discretized over the step, or additiveNoise(), a fixed Q, not "
        "both"
    );
    return true;
}

template <typename Matrix>
class LinearStep
{
public:
 
    static constexpr bool IsAffine = true;  

    explicit LinearStep(Matrix matrix)
        : m_matrix{std::move(matrix)}
    {}

    [[nodiscard]] Matrix const &
    matrix() const
    {
        return m_matrix;
    }

    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] auto
    operator()([[maybe_unused]] In const & in, Out const &, X const & x, TimeScalar) const
    {
        static_assert(
            blocksAddValues<In>,
            "a matrix process model moves a state whose coordinates add their values; a state with a block such as a "
            "unit quaternion takes a function"
        );
        static_assert(
            Matrix::ColsAtCompileTime == Eigen::Dynamic || valueDimensionAtCompileTime<In> == Eigen::Dynamic
                || (Matrix::ColsAtCompileTime == valueDimensionAtCompileTime<In>
                    && Matrix::RowsAtCompileTime == valueDimensionAtCompileTime<In>),
            "a matrix process model's matrix is square, with a row and a column per coordinate of the state"
        );
        EMRYS_FILTER_ASSERT(
            m_matrix.rows() == x.size() && m_matrix.cols() == x.size()
            && "a matrix process model's matrix is square, with a row and a column per coordinate of the state"
        );
        using T = typename X::Scalar;
        using Value = Eigen::Matrix<T, Matrix::RowsAtCompileTime, 1, Eigen::ColMajor, Matrix::MaxRowsAtCompileTime, 1>;
        Value value = m_matrix.template cast<T>() * x;
        return WithJacobian<Value, Matrix const &>{std::move(value), m_matrix};
    }
 
private:
 
    Matrix m_matrix;  
};

template <bool Whole, typename Writes>
struct ProcessWrites
{
    using Paths = typename Writes::PathList;
    using OutputSpace = typename Writes::OutputSpace;
};

template <typename Writes>
struct ProcessWrites<true, Writes>
{
    using Paths = TagPaths<WholeSpace>;
    using OutputSpace = WholeSpace;
};
 
}  // namespace internal

template <typename F, typename... Policies>
class ProcessFunction
{
    using Set = internal::PolicySet<Policies...>;
 
public:

    using FunctionType = F;
 
    static constexpr bool IsModel = true;      
    static constexpr bool HasFunction = true;  

    template <internal::PolicyKind K>
    static constexpr bool has = Set::template has<K>;

    template <internal::PolicyKind K>
    using Policy = typename Set::template Type<K>;
 
    // processModel() and the placements give both reads() and writes(), or neither.
    static_assert(has<internal::PolicyKind::reads> == has<internal::PolicyKind::writes>);

    static constexpr bool Whole = !has<internal::PolicyKind::writes>;

    using WritePaths = typename internal::ProcessWrites<Whole, Policy<internal::PolicyKind::writes>>::Paths;

    using OutputSpace = typename internal::ProcessWrites<Whole, Policy<internal::PolicyKind::writes>>::OutputSpace;

    explicit ProcessFunction(F function, Policies... policies)
        : m_function{std::move(function)}
        , m_policies{std::move(policies)...}
    {}

    [[nodiscard]] F const &
    function() const
    {
        return m_function;
    }

    template <internal::PolicyKind K>
    [[nodiscard]] Policy<K> const &
    policy() const
    {
        return Set::template get<K>(m_policies);
    }

    [[nodiscard]] std::tuple<Policies...> const &
    policies() const
    {
        return m_policies;
    }

    template <typename Scalar, typename Out>
    [[nodiscard]] typename Moments<Out, Scalar>::MatrixType
    noiseOver(Out const & out, [[maybe_unused]] Scalar dt) const
    {
        using Matrix = typename Moments<Out, Scalar>::MatrixType;
        Eigen::Index const n = differenceDimension(out);
        Matrix q = Matrix::Zero(n, n);
        if constexpr (requires { m_function.template noiseOver<Scalar>(out, dt); })
        {
            q = m_function.template noiseOver<Scalar>(out, dt);
        }
        else if constexpr (has<internal::PolicyKind::discreteNoise>)
        {
            q = policy<internal::PolicyKind::discreteNoise>().function(out, dt);
        }
        else if constexpr (has<internal::PolicyKind::continuousNoise>)
        {
            std::apply(
                [&](auto const &... parts) {
                    (internal::addContinuousPart(out, q, parts, dt), ...);
                },
                policy<internal::PolicyKind::continuousNoise>().parts()
            );
        }
        else if constexpr (has<internal::PolicyKind::additiveNoise>)
        {
            auto const placed = internal::placeNoise(out, policy<internal::PolicyKind::additiveNoise>());
            using Placed = std::remove_cvref_t<decltype(placed)>;
            static_assert(
                !internal::isGaussianMixture<typename Placed::Distribution>,
                "a process model's additive noise is one Gaussian; a mixture over the step is a mixture argument of "
                "predict() (roadmap 17z)"
            );
            auto const & distribution = placed.distribution();
            q = internal::symmetricFromLower<Matrix>(distribution.covarianceLower().template cast<Scalar>());
        }
        return q;
    }

    template <typename Space>
    [[nodiscard]] static consteval internal::StaticPartition
    staticPartition()
    {
        if constexpr (Whole)
        {
            return internal::StaticPartition::holds;
        }
        else
        {
            return internal::staticPartitionOfWrites<Space, WritePaths>();
        }
    }

    template <typename Space>
    [[nodiscard]] bool
    partitions([[maybe_unused]] Space const & space) const
    {
        if constexpr (Whole)
        {
            return true;
        }
        else
        {
            return internal::writesPartition(space, static_cast<WritePaths *>(nullptr));
        }
    }
 
private:
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS F m_function;  
    std::tuple<Policies...> m_policies;           
};

template <typename F, typename... Policies>
    requires(!MatrixExpression<F> && !internal::ModelPolicy<F>)
[[nodiscard]] auto
processModel(F && function, Policies &&... policies)
{
    using Set = internal::PolicySet<std::remove_cvref_t<Policies>...>;
    static_assert(internal::checkProcessPolicies<Set>());
    return ProcessFunction<std::decay_t<F>, std::remove_cvref_t<Policies>...>{
        std::forward<F>(function),
        std::forward<Policies>(policies)...
    };
}

template <MatrixExpression Matrix, typename... Policies>
[[nodiscard]] auto
processModel(Matrix && matrix, Policies &&... policies)
{
    using Plain = internal::PlainResult<std::remove_cvref_t<Matrix>>;
    using Set = internal::PolicySet<std::remove_cvref_t<Policies>...>;
    static_assert(
        !Set::template has<internal::PolicyKind::reads> && !Set::template has<internal::PolicyKind::writes>,
        "a matrix process model is over the whole state, so it takes no reads() or writes(); give a function for a "
        "model of some blocks"
    );
    static_assert(internal::checkProcessPolicies<Set>());
    return ProcessFunction<internal::LinearStep<Plain>, std::remove_cvref_t<Policies>...>{
        internal::LinearStep<Plain>{Plain{std::forward<Matrix>(matrix)}},
        std::forward<Policies>(policies)...
    };
}
 
namespace internal
{

template <typename F, typename... Policies>
std::true_type
processFunctionBase(ProcessFunction<F, Policies...> const *);

std::false_type
processFunctionBase(void const *);

template <typename Path>
using ReadsOne = Reads<ReadItem<AsTagPath<Path>, PathTag<Path>>>;

template <typename... Paths>
using ReadsOf = Reads<ReadItem<AsTagPath<Paths>, PathTag<Paths>>...>;

template <typename... Paths>
using WritesOf = Writes<AsTagPath<Paths>...>;

struct IdentityStep
{
    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] auto
    operator()(In const &, Out const & out, X const & x, TimeScalar) const
    {
        using T = ScalarOf<X>;
        using JacobianMatrix = typename Moments<Out, TimeScalar>::MatrixType;
        Eigen::Index const n = differenceDimension(out);
        return withJacobian(typename Moments<Out, T>::VectorType{x}, JacobianMatrix{JacobianMatrix::Identity(n, n)});
    }
};

template <int Dimension, typename Scalar>
struct RandomWalkNoise
{
    Eigen::Matrix<Scalar, Dimension, Dimension> spectralDensity;  

    template <typename Out, typename TimeScalar>
    [[nodiscard]] typename Moments<Out, TimeScalar>::MatrixType
    operator()(Out const &, TimeScalar dt) const
    {
        return typename Moments<Out, TimeScalar>::MatrixType{dt * spectralDensity.template cast<TimeScalar>()};
    }
};

[[nodiscard]] constexpr int
factorial(int n)
{
    return n <= 1 ? 1 : n * factorial(n - 1);
}

template <int N, typename TimeScalar>
[[nodiscard]] std::array<TimeScalar, N>
powersOf(TimeScalar dt)
{
    auto result = std::array<TimeScalar, N>{};
    result[0] = TimeScalar(1);
    for (int k = 1; k < N; ++k)
    {
        result[k] = result[k - 1] * dt;
    }
    return result;
}

template <int Dimension, typename... Tags>
struct KinematicStep
{
    static constexpr int Order = int(sizeof...(Tags));  

    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] auto
    operator()(In const & in, Out const & out, X const & x, TimeScalar dt) const
    {
        using T = ScalarOf<X>;
        using Vector = typename Moments<Out, T>::VectorType;
        using JacobianMatrix = typename Moments<Out, TimeScalar>::MatrixType;
        auto const starts = std::array<Eigen::Index, Order>{blockRows<Tags>(in).start...};
        auto const differenceStarts = std::array<Eigen::Index, Order>{differenceBlockRows<Tags>(out).start...};
        auto const dtPowers = powersOf<Order>(dt);
        Eigen::Index const n = differenceDimension(out);
        Vector y{x};
        JacobianMatrix jacobian = JacobianMatrix::Identity(n, n);
        for (int i = 0; i < Order; ++i)
        {
            for (int j = i + 1; j < Order; ++j)
            {
                TimeScalar const scale = dtPowers[j - i] / TimeScalar(factorial(j - i));
                y.template segment<Dimension>(starts[i]) += T(scale) * x.template segment<Dimension>(starts[j]);
                jacobian.template block<Dimension, Dimension>(differenceStarts[i], differenceStarts[j])
                    .diagonal()
                    .setConstant(scale);
            }
        }
        return withJacobian(std::move(y), std::move(jacobian));
    }
};

template <int Dimension, typename Scalar, typename... Tags>
struct KinematicNoise
{
    static constexpr int Order = int(sizeof...(Tags));  
 
    Eigen::Matrix<Scalar, Dimension, Dimension> spectralDensity;  

    [[nodiscard]] static constexpr int
    denominator(int i, int j)
    {
        return (2 * Order - 1 - i - j) * factorial(Order - 1 - i) * factorial(Order - 1 - j);
    }

    template <typename Out, typename TimeScalar>
    [[nodiscard]] typename Moments<Out, TimeScalar>::MatrixType
    operator()(Out const & out, TimeScalar dt) const
    {
        using Matrix = typename Moments<Out, TimeScalar>::MatrixType;
        auto const starts = std::array<Eigen::Index, Order>{differenceBlockRows<Tags>(out).start...};
        auto const dtPowers = powersOf<2 * Order>(dt);
        Eigen::Index const n = differenceDimension(out);
        Matrix q = Matrix::Zero(n, n);
        for (int i = 0; i < Order; ++i)
        {
            for (int j = 0; j < Order; ++j)
            {
                TimeScalar const scale = dtPowers[2 * Order - 1 - i - j] / TimeScalar(denominator(i, j));
                q.template block<Dimension, Dimension>(starts[i], starts[j]) =
                    scale * spectralDensity.template cast<TimeScalar>();
            }
        }
        return q;
    }
};

template <typename Scalar>
struct OrnsteinUhlenbeckStep
{
    Scalar timeConstant;  

    template <typename TimeScalar>
    [[nodiscard]] TimeScalar
    decay(TimeScalar dt) const
    {
        return Eigen::numext::exp(TimeScalar(-dt / TimeScalar(timeConstant)));
    }

    template <typename In, typename Out, typename X, typename TimeScalar>
    [[nodiscard]] auto
    operator()(In const &, Out const & out, X const & x, TimeScalar dt) const
    {
        using T = ScalarOf<X>;
        using JacobianMatrix = typename Moments<Out, TimeScalar>::MatrixType;
        TimeScalar const phi = decay(dt);
        Eigen::Index const n = differenceDimension(out);
        return withJacobian(
            typename Moments<Out, T>::VectorType{T(phi) * x},
            JacobianMatrix{phi * JacobianMatrix::Identity(n, n)}
        );
    }
};

template <typename Scalar>
struct OrnsteinUhlenbeckNoise
{
    Scalar timeConstant;  
    Scalar variance;      

    template <typename Out, typename TimeScalar>
    [[nodiscard]] typename Moments<Out, TimeScalar>::MatrixType
    operator()(Out const & out, TimeScalar dt) const
    {
        using Matrix = typename Moments<Out, TimeScalar>::MatrixType;
        TimeScalar const phi = Eigen::numext::exp(TimeScalar(-dt / TimeScalar(timeConstant)));
        Eigen::Index const n = differenceDimension(out);
        return Matrix{TimeScalar(variance) * (1 - phi * phi) * Matrix::Identity(n, n)};
    }
};
 
}  // namespace internal

template <typename Path, RealScalar Scalar_ = double>
class RandomWalk
    : public ProcessFunction<
          internal::IdentityStep,
          internal::ReadsOne<Path>,
          internal::WritesOf<Path>,
          AffineIn<internal::PathTag<Path>>,
          internal::DiscreteNoise<internal::RandomWalkNoise<
              internal::DifferenceSpaceOf<internal::PathTag<Path>>::type::FixedValueDimension,
              Scalar_>>>
{
    using Differences = typename internal::DifferenceSpaceOf<internal::PathTag<Path>>::type;
 
    static_assert(Differences::ValueDimensionStep == 0, "the block of a random walk must have a fixed size");

    using Noise = internal::RandomWalkNoise<Differences::FixedValueDimension, Scalar_>;

    using Base = ProcessFunction<
        internal::IdentityStep,
        internal::ReadsOne<Path>,
        internal::WritesOf<Path>,
        AffineIn<internal::PathTag<Path>>,
        internal::DiscreteNoise<Noise>>;
 
public:

    static constexpr int Dimension = Differences::FixedValueDimension;

    using Scalar = Scalar_;

    using SpectralDensityMatrix = Eigen::Matrix<Scalar_, Dimension, Dimension>;

    explicit RandomWalk(Scalar_ spectralDensity)
        : RandomWalk{SpectralDensityMatrix{spectralDensity * SpectralDensityMatrix::Identity()}}
    {}

    explicit RandomWalk(SpectralDensityMatrix const & spectralDensity)
        : Base{
              internal::IdentityStep{},
              internal::ReadsOne<Path>{},
              internal::WritesOf<Path>{},
              AffineIn<internal::PathTag<Path>>{},
              internal::DiscreteNoise<Noise>{Noise{spectralDensity}}
          }
    {}

    [[nodiscard]] SpectralDensityMatrix const &
    spectralDensity() const
    {
        return this->template policy<internal::PolicyKind::discreteNoise>().function.spectralDensity;
    }
};

template <RealScalar Scalar_, typename... Paths>
class Kinematic
    : public ProcessFunction<
          internal::KinematicStep<
              internal::PathTag<std::tuple_element_t<0, std::tuple<Paths...>>>::FixedValueDimension,
              internal::PathTag<Paths>...>,
          internal::ReadsOf<Paths...>,
          internal::WritesOf<Paths...>,
          AffineIn<internal::PathTag<Paths>...>,
          internal::DiscreteNoise<internal::KinematicNoise<
              internal::PathTag<std::tuple_element_t<0, std::tuple<Paths...>>>::FixedValueDimension,
              Scalar_,
              internal::PathTag<Paths>...>>>
{
    static_assert(sizeof...(Paths) > 1, "a kinematic chain needs at least two blocks; one block is a RandomWalk");

    using First = internal::PathTag<std::tuple_element_t<0, std::tuple<Paths...>>>;
 
public:

    static constexpr int Order = int(sizeof...(Paths));

    static constexpr int Dimension = First::FixedValueDimension;
 
    static_assert(
        ((internal::PathTag<Paths>::ValueDimensionStep == 0
          && internal::PathTag<Paths>::FixedValueDimension == Dimension)
         && ...),
        "the blocks of a kinematic chain must have the same fixed size"
    );
 
    static_assert(
        (internal::hasEqualSizes<internal::PathTag<Paths>> && ...),
        "the blocks of a kinematic chain must have differences of the same size as their values"
    );

    using Scalar = Scalar_;

    using SpectralDensityMatrix = Eigen::Matrix<Scalar_, Dimension, Dimension>;
 
private:

    using Step = internal::KinematicStep<Dimension, internal::PathTag<Paths>...>;

    using Noise = internal::KinematicNoise<Dimension, Scalar_, internal::PathTag<Paths>...>;

    using Base = ProcessFunction<
        Step,
        internal::ReadsOf<Paths...>,
        internal::WritesOf<Paths...>,
        AffineIn<internal::PathTag<Paths>...>,
        internal::DiscreteNoise<Noise>>;
 
public:

    explicit Kinematic(Scalar_ spectralDensity)
        : Kinematic{SpectralDensityMatrix{spectralDensity * SpectralDensityMatrix::Identity()}}
    {}

    explicit Kinematic(SpectralDensityMatrix const & spectralDensity)
        : Base{
              Step{},
              internal::ReadsOf<Paths...>{},
              internal::WritesOf<Paths...>{},
              AffineIn<internal::PathTag<Paths>...>{},
              internal::DiscreteNoise<Noise>{Noise{spectralDensity}}
          }
    {}

    [[nodiscard]] SpectralDensityMatrix const &
    spectralDensity() const
    {
        return this->template policy<internal::PolicyKind::discreteNoise>().function.spectralDensity;
    }
};

template <typename Position, typename Velocity, RealScalar Scalar = double>
using ConstantVelocity = Kinematic<Scalar, Position, Velocity>;

template <typename Position, typename Velocity, typename Acceleration, RealScalar Scalar = double>
using ConstantAcceleration = Kinematic<Scalar, Position, Velocity, Acceleration>;

template <RealScalar Scalar = double>
struct TimeConstant
{
    Scalar value;  
};

template <RealScalar Scalar = double>
struct StationaryVariance
{
    Scalar value;  
};

template <typename Value>
    requires std::is_arithmetic_v<Value>
TimeConstant(Value) -> TimeConstant<internal::DeducedScalar<Value>>;
 
template <typename Value>
    requires std::is_arithmetic_v<Value>
StationaryVariance(Value) -> StationaryVariance<internal::DeducedScalar<Value>>;


template <typename Path, RealScalar Scalar_ = double>
class OrnsteinUhlenbeck
    : public ProcessFunction<
          internal::OrnsteinUhlenbeckStep<Scalar_>,
          internal::ReadsOne<Path>,
          internal::WritesOf<Path>,
          AffineIn<internal::PathTag<Path>>,
          internal::DiscreteNoise<internal::OrnsteinUhlenbeckNoise<Scalar_>>>
{
    using Base = ProcessFunction<
        internal::OrnsteinUhlenbeckStep<Scalar_>,
        internal::ReadsOne<Path>,
        internal::WritesOf<Path>,
        AffineIn<internal::PathTag<Path>>,
        internal::DiscreteNoise<internal::OrnsteinUhlenbeckNoise<Scalar_>>>;
 
public:

    using Scalar = Scalar_;

    OrnsteinUhlenbeck(TimeConstant<Scalar_> timeConstant, StationaryVariance<Scalar_> variance)
        : Base{
              internal::OrnsteinUhlenbeckStep<Scalar_>{timeConstant.value},
              internal::ReadsOne<Path>{},
              internal::WritesOf<Path>{},
              AffineIn<internal::PathTag<Path>>{},
              internal::DiscreteNoise<internal::OrnsteinUhlenbeckNoise<Scalar_>>{
                  internal::OrnsteinUhlenbeckNoise<Scalar_>{timeConstant.value, variance.value}
              }
          }
    {
        EMRYS_FILTER_ASSERT(timeConstant.value > 0 && "the time constant must be positive");
        EMRYS_FILTER_ASSERT(variance.value >= 0 && "the variance must be nonnegative");
    }

    [[nodiscard]] Scalar_
    timeConstant() const
    {
        return this->function().timeConstant;
    }

    [[nodiscard]] Scalar_
    variance() const
    {
        return this->template policy<internal::PolicyKind::discreteNoise>().function.variance;
    }
};

template <typename Path>
class Constant
    : public ProcessFunction<
          internal::IdentityStep,
          internal::ReadsOne<Path>,
          internal::WritesOf<Path>,
          AffineIn<internal::PathTag<Path>>>
{
    using Base = ProcessFunction<
        internal::IdentityStep,
        internal::ReadsOne<Path>,
        internal::WritesOf<Path>,
        AffineIn<internal::PathTag<Path>>>;
 
public:

    Constant()
        : Base{
              internal::IdentityStep{},
              internal::ReadsOne<Path>{},
              internal::WritesOf<Path>{},
              AffineIn<internal::PathTag<Path>>{}
          }
    {}
};
 
namespace internal
{

template <typename R>
struct ItemsOf;

template <typename... Items>
struct ItemsOf<Reads<Items...>>
{
    using type = std::tuple<Items...>;
};


template <typename Acc, typename... Items>
struct ReadUnion
{
    using type = Acc;  
};

template <typename... Kept, typename First, typename... Rest>
struct ReadUnion<Reads<Kept...>, First, Rest...>
{
    static constexpr bool present = (std::is_same_v<typename Kept::NameType, typename First::NameType> || ...);
 
    static_assert(
        ((!std::is_same_v<typename Kept::NameType, typename First::NameType>
          || std::is_same_v<typename Kept::PathType, typename First::PathType>)
         && ...),
        "two parts of a composition read different blocks under one name; read one under another name with "
        "on<Name>(path) in its reads()"
    );
 
    using type = typename ReadUnion<std::conditional_t<present, Reads<Kept...>, Reads<Kept..., First>>, Rest...>::type;
};


template <typename Tuple>
struct ReadUnionOf;

template <typename... Items>
struct ReadUnionOf<std::tuple<Items...>>
{
    using type = typename ReadUnion<Reads<>, Items...>::type;
};


template <typename Paths, typename Names>
struct WritesOfTuple;

template <typename... Paths, typename... Names>
struct WritesOfTuple<std::tuple<Paths...>, std::tuple<Names...>>
{
    static_assert(distinctTypes<Paths...>(), "two parts of a composition write the same block");
    static_assert(
        distinctTypes<Names...>(),
        "two parts of a composition write blocks under one name, their paths' last tag; place each model on its "
        "block with on<Path>(model)"
    );
    using type = WritesAs<TagPaths<Paths...>, Composite<Names...>>;
};


template <typename Items, typename Name>
inline constexpr bool itemsName = false;

// A partial specialization rather than a lambda invoked in the initializer, which MSVC 19.50
// fails on with an internal compiler error.
template <typename... Items, typename Name>
inline constexpr bool itemsName<std::tuple<Items...>, Name> = (std::is_same_v<typename Items::NameType, Name> || ...);


template <typename Model, typename Name>
inline constexpr bool modelReads =
    itemsName<typename ItemsOf<typename Model::template Policy<PolicyKind::reads>>::type, Name>;

template <typename Model, typename Name>
inline constexpr bool modelAffineIn =
    !modelReads<Model, Name> || affineNames<typename Model::template Policy<PolicyKind::affineIn>, Name>;

template <typename Acc, typename Models, typename... Names>
struct AffineUnion
{
    using type = Acc;  
};

template <typename... Kept, typename... Models, typename First, typename... Rest>
struct AffineUnion<AffineIn<Kept...>, std::tuple<Models...>, First, Rest...>
{
    static constexpr bool affine = (modelAffineIn<Models, First> && ...);
 
    using type = typename AffineUnion<
        std::conditional_t<affine, AffineIn<Kept..., First>, AffineIn<Kept...>>,
        std::tuple<Models...>,
        Rest...>::type;
};


template <typename Union, typename Models>
struct AffineOf;

template <typename... Items, typename Models>
struct AffineOf<Reads<Items...>, Models>
{
    using type = typename AffineUnion<AffineIn<>, Models, typename Items::NameType...>::type;
};


template <typename Model, typename Scalar, typename... Extra>
inline constexpr bool returnsJacobianWith = isWithJacobian<std::remove_cvref_t<decltype(callForm(
    std::declval<typename Model::FunctionType const &>(),
    std::declval<typename Model::template Policy<PolicyKind::reads>::InputSpace const &>(),
    std::declval<OutputSpaceOf<Model> const &>(),
    std::declval<
        typename Moments<typename Model::template Policy<PolicyKind::reads>::InputSpace, Scalar>::VectorType const &>(),
    std::declval<Extra>()...
))>>;

template <typename Model, typename Scalar>
inline constexpr bool returnsStepJacobian = returnsJacobianWith<Model, Scalar, Scalar>;

template <typename Model, typename Space, typename Scalar>
inline constexpr bool returnsStepJacobianIn = [] {
    if constexpr (readsWholeState<Model>)
    {
        return isWithJacobian<std::remove_cvref_t<decltype(callForm(
            std::declval<typename Model::FunctionType const &>(),
            std::declval<Space const &>(),
            std::declval<Space const &>(),
            std::declval<typename Moments<Space, Scalar>::VectorType const &>(),
            std::declval<Scalar>()
        ))>>;
    }
    else
    {
        return returnsStepJacobian<Model, Scalar>;
    }
}();

template <typename Model>
using ReadItemsOf = typename ItemsOf<typename Model::template Policy<PolicyKind::reads>>::type;

template <typename Part, typename In>
[[nodiscard]] auto
partIn(In const & in)
{
    using PartIn = typename Part::template Policy<PolicyKind::reads>::InputSpace;
    return [&]<typename... Items>(std::tuple<Items...> *) {
        return PartIn{pathBlock<typename Items::NameType>(in)...};
    }(static_cast<ReadItemsOf<Part> *>(nullptr));
}

template <typename Part, typename In, typename PartIn, typename X>
[[nodiscard]] auto
partVector(In const & in, PartIn const & pin, X const & x)
{
    auto px = zeros<typename X::Scalar>(pin);
    [&]<typename... Items>(std::tuple<Items...> *) {
        ((block<typename Items::NameType>(pin, px) = block<typename Items::NameType>(in, x)), ...);
    }(static_cast<ReadItemsOf<Part> *>(nullptr));
    return px;
}

template <typename Part, typename Out>
[[nodiscard]] OutputSpaceOf<Part>
partOut(Out const & out)
{
    return [&]<typename... Names>(std::tuple<Names...> *) {
        return OutputSpaceOf<Part>{pathBlock<Names>(out)...};
    }(static_cast<CompositeBlocks<OutputSpaceOf<Part>> *>(nullptr));
}

template <
    typename Row,
    typename Out,
    typename In,
    typename J,
    typename POut,
    typename PIn,
    typename PJ,
    typename... Items>
void
partJacobianRow(
    Out const & out,
    In const & in,
    J & jacobian,
    POut const & pout,
    PIn const & pin,
    PJ const & partJacobian,
    std::tuple<Items...> *
)
{
    ((diffBlock<Row, typename Items::NameType>(out, in, jacobian) =
          diffBlock<Row, typename Items::NameType>(pout, pin, partJacobian)),
     ...);
}

template <
    std::size_t I,
    typename Parts,
    typename In,
    typename Out,
    typename X,
    typename Y,
    typename J,
    typename... Extra>
void
stepPart(
    Parts const & parts,
    In const & in,
    Out const & out,
    X const & x,
    Y & y,
    [[maybe_unused]] J & jacobian,
    Extra const &... extra
)
{
    using Part = std::tuple_element_t<I, Parts>;
    auto const pin = partIn<Part>(in);
    auto const pout = partOut<Part>(out);
    auto const px = partVector<Part>(in, pin, x);
    auto result = vectorReturned(callForm(std::get<I>(parts).function(), pin, pout, px, extra...));
    auto const & value = valueOf(result);
    using Names = CompositeBlocks<OutputSpaceOf<Part>>;
    [&]<typename... Rows>(std::tuple<Rows...> *) {
        ((block<Rows>(out, y) = block<Rows>(pout, value)), ...);
        if constexpr (!std::is_same_v<J, NoPolicy>)
        {
            (partJacobianRow<
                 Rows>(out, in, jacobian, pout, pin, result.jacobian, static_cast<ReadItemsOf<Part> *>(nullptr)),
             ...);
        }
    }(static_cast<Names *>(nullptr));
}

template <typename JacobianScalar, typename... Models, typename In, typename Out, typename X, typename... Extra>
[[nodiscard]] auto
stepParts(std::tuple<Models...> const & parts, In const & in, Out const & out, X const & x, Extra const &... extra)
{
    constexpr bool AllJacobians = (returnsJacobianWith<Models, JacobianScalar, Extra...> && ...);
    auto y = zeros<typename X::Scalar>(out);
    auto jacobian = [&] {
        if constexpr (AllJacobians)
        {
            using JacobianMatrix = MatrixOf<
                JacobianScalar,
                valueDimensionAtCompileTime<DifferenceSpace<Out>>,
                valueDimensionAtCompileTime<DifferenceSpace<In>>,
                maxValueDimensionAtCompileTime<DifferenceSpace<Out>>,
                maxValueDimensionAtCompileTime<DifferenceSpace<In>>>;
            return JacobianMatrix{JacobianMatrix::Zero(differenceDimension(out), differenceDimension(in))};
        }
        else
        {
            return NoPolicy{};
        }
    }();
    [&]<std::size_t... I>(std::index_sequence<I...>) {
        (stepPart<I>(parts, in, out, x, y, jacobian, extra...), ...);
    }(std::index_sequence_for<Models...>{});
    if constexpr (AllJacobians)
    {
        return withJacobian(std::move(y), std::move(jacobian));
    }
    else
    {
        return y;
    }
}

template <typename... Models>
inline constexpr bool anyTakesBatches =
    (takesBatches<std::remove_cvref_t<decltype(std::declval<Models const &>().function())>> || ...);

template <typename... Models, typename In, typename Out, typename X, typename Y, typename... Extra>
void
stepPartsBatch(
    std::tuple<Models...> const & parts,
    In const & in,
    Out const & out,
    X const & x,
    Y & y,
    Extra const &... extra
)
{
    using T = typename X::Scalar;
    // Used in the parts' types only, which MSVC 19.51 takes for unused (C4189).
    [[maybe_unused]] constexpr int Cols = X::ColsAtCompileTime;
    [[maybe_unused]] constexpr int MaxCols = X::MaxColsAtCompileTime;
    Eigen::Index const count = x.cols();
    auto & values = y.values();
    values.setZero();
    [&]<std::size_t... I>(std::index_sequence<I...>) {
        (
            [&] {
                using Part = std::tuple_element_t<I, std::tuple<Models...>>;
                auto const pin = partIn<Part>(in);
                auto const pout = partOut<Part>(out);
                using PIn = std::remove_cvref_t<decltype(pin)>;
                using POut = OutputSpaceOf<Part>;
                using PX =
                    MatrixOf<T, valueDimensionAtCompileTime<PIn>, Cols, maxValueDimensionAtCompileTime<PIn>, MaxCols>;
                using PY =
                    MatrixOf<T, valueDimensionAtCompileTime<POut>, Cols, maxValueDimensionAtCompileTime<POut>, MaxCols>;
                auto px = PX{};
                px.resize(valueDimension(pin), count);
                for (Eigen::Index j = 0; j < count; ++j)
                {
                    px.col(j) = partVector<Part>(in, pin, x.col(j));
                }
                auto py = PY{};
                py.resize(valueDimension(pout), count);
                auto const & function = std::get<I>(parts).function();
                if constexpr (takesBatches<std::remove_cvref_t<decltype(function)>>)
                {
                    callBatched(function, pin, pout, px, py, extra...);
                }
                else
                {
                    for (Eigen::Index j = 0; j < count; ++j)
                    {
                        typename Moments<PIn, T>::VectorType const column = px.col(j);
                        py.col(j) = valueOf(vectorReturned(callForm(function, pin, pout, column, extra...)));
                    }
                }
                [&]<typename... Rows>(std::tuple<Rows...> *) {
                    ((block<Rows>(out, values) = block<Rows>(pout, py)), ...);
                }(static_cast<CompositeBlocks<POut> *>(nullptr));
            }(),
            ...
        );
    }(std::index_sequence_for<Models...>{});
}
 
}  // namespace internal
 
template <typename... Models>
class Composition;
 
namespace internal
{

template <typename Model>
inline constexpr bool isComposition = false;


template <typename... Models>
inline constexpr bool isComposition<Composition<Models...>> = true;


template <typename Model>
struct DeclaredScalars
{
    using type = std::tuple<>;  
};


template <typename Model>
    requires(!isComposition<Model> && requires { typename Model::Scalar; })
struct DeclaredScalars<Model>
{
    using type = std::tuple<typename Model::Scalar>;  
};

template <typename... Models>
struct DeclaredScalars<Composition<Models...>>
{
    using type = decltype(std::tuple_cat(std::declval<typename DeclaredScalars<Models>::type>()...));  
};


template <typename... Models>
using DeclaredScalarsOf = decltype(std::tuple_cat(std::declval<typename DeclaredScalars<Models>::type>()...));

template <typename... Models>
[[nodiscard]] consteval bool
oneDeclaredScalar()
{
    return []<typename... Scalars>(std::tuple<Scalars...> *) {
        if constexpr (sizeof...(Scalars) == 0)
        {
            return true;
        }
        else
        {
            using First = std::tuple_element_t<0, std::tuple<Scalars...>>;
            return (std::is_same_v<First, Scalars> && ...);
        }
    }(static_cast<DeclaredScalarsOf<Models...> *>(nullptr));
}

template <typename Scalars>
struct FirstScalarOr
{
    using type = double;  
};


template <typename First, typename... Rest>
struct FirstScalarOr<std::tuple<First, Rest...>>
{
    using type = First;  
};

 
}  // namespace internal

template <typename... Models>
class Composition
{
    static_assert(sizeof...(Models) > 0, "a composition needs a process model");
    static_assert(
        internal::oneDeclaredScalar<Models...>(),
        "the parts of a composition must have one scalar: give every built-in model the same Scalar"
    );

    using AllItems = decltype(std::tuple_cat(
        std::declval<
            typename internal::ItemsOf<typename Models::template Policy<internal::PolicyKind::reads>>::type>()...
    ));

    using AllPaths =
        decltype(std::tuple_cat(std::declval<typename internal::PathsOf<typename Models::WritePaths>::type>()...));

    using AllNames =
        decltype(std::tuple_cat(std::declval<internal::CompositeBlocks<typename Models::OutputSpace>>()...));
 
public:

    using Scalar = typename internal::FirstScalarOr<internal::DeclaredScalarsOf<Models...>>::type;
 
    using FunctionType = Composition;  
 
    static constexpr bool IsModel = true;      
    static constexpr bool HasFunction = true;  

    using ReadsPolicy = typename internal::ReadUnionOf<AllItems>::type;

    using WritesPolicy = typename internal::WritesOfTuple<AllPaths, AllNames>::type;

    using AffinePolicy = typename internal::AffineOf<ReadsPolicy, std::tuple<Models...>>::type;

    template <internal::PolicyKind K>
    static constexpr bool has =
        K == internal::PolicyKind::reads || K == internal::PolicyKind::writes || K == internal::PolicyKind::affineIn;

    template <internal::PolicyKind K>
    using Policy = std::conditional_t<
        K == internal::PolicyKind::reads,
        ReadsPolicy,
        std::conditional_t<
            K == internal::PolicyKind::writes,
            WritesPolicy,
            std::conditional_t<K == internal::PolicyKind::affineIn, AffinePolicy, internal::NoPolicy>>>;

    using WritePaths = typename WritesPolicy::PathList;

    using OutputSpace = typename WritesPolicy::OutputSpace;

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

    [[nodiscard]] Composition const &
    function() const
    {
        return *this;
    }

    template <internal::PolicyKind K>
    [[nodiscard]] Policy<K> const &
    policy() const
    {
        static_assert(has<K>, "the composition was not given this policy");
        static constexpr Policy<K> value{};
        return value;
    }

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

    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>
    {
        return internal::stepParts<TimeScalar>(m_parts, in, out, x, dt);
    }

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

    using MarkType = internal::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
    {
        internal::stepPartsBatch(m_parts, in, out, x, y, dt);
    }

    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);
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (noisePart<I, Scalar>(out, dt, q), ...);
        }(std::index_sequence_for<Models...>{});
        return q;
    }

    template <typename Space>
    [[nodiscard]] static consteval internal::StaticPartition
    staticPartition()
    {
        return internal::staticPartitionOfWrites<Space, WritePaths>();
    }

    template <typename Space>
    [[nodiscard]] bool
    partitions(Space const & space) const
    {
        return internal::writesPartition(space, static_cast<WritePaths *>(nullptr));
    }
 
private:

    template <std::size_t I>
    using Part = std::tuple_element_t<I, std::tuple<Models...>>;

    template <std::size_t I, typename Scalar, typename Out, typename Q>
    void
    noisePart(Out const & out, Scalar dt, Q & q) const
    {
        auto const pout = internal::partOut<Part<I>>(out);
        auto const partNoise = std::get<I>(m_parts).template noiseOver<Scalar>(pout, dt);
        using Names = internal::CompositeBlocks<typename Part<I>::OutputSpace>;
        [&]<typename... Rows>(std::tuple<Rows...> *) {
            (noiseRow<Rows>(out, q, pout, partNoise, static_cast<Names *>(nullptr)), ...);
        }(static_cast<Names *>(nullptr));
    }

    template <typename Row, typename Out, typename Q, typename POut, typename PQ, typename... Columns>
    static void
    noiseRow(Out const & out, Q & q, POut const & pout, PQ const & partNoise, std::tuple<Columns...> *)
    {
        ((diffBlock<Row, Columns>(out, out, q) = diffBlock<Row, Columns>(pout, pout, partNoise)), ...);
    }
 
    std::tuple<Models...> m_parts;  
};
 
namespace internal
{

template <typename... Models>
std::true_type
compositionBase(Composition<Models...> const *);

std::false_type
compositionBase(void const *);
 
}  // namespace internal

template <typename M>
concept ProcessModel = decltype(internal::processFunctionBase(std::declval<std::remove_cvref_t<M> const *>()))::value
                    || decltype(internal::compositionBase(std::declval<std::remove_cvref_t<M> const *>()))::value;

struct NoProcessModel
{
    template <typename Space>
    [[nodiscard]] static consteval internal::StaticPartition
    staticPartition()
    {
        if constexpr (Space::ValueDimensionStep == 0)
        {
            return Space::FixedValueDimension == 0 ? internal::StaticPartition::holds
                                                   : internal::StaticPartition::fails;
        }
        else
        {
            return internal::StaticPartition::unknown;
        }
    }

    template <typename Space>
    [[nodiscard]] bool
    partitions(Space const & space) const
    {
        return valueDimension(space) == 0;
    }
};

template <ProcessModel... Models>
[[nodiscard]] Composition<std::remove_cvref_t<Models>...>
compose(Models &&... models)
{
    static_assert((internal::wholeStateRefused<std::remove_cvref_t<Models>> && ...));
    return Composition<std::remove_cvref_t<Models>...>{std::forward<Models>(models)...};
}
 
namespace internal
{

template <typename Model, typename Space>
[[nodiscard]] consteval bool
readsOnlyState()
{
    if constexpr (readsWholeState<Model>)
    {
        return true;
    }
    else
    {
        return []<typename... Items>(std::tuple<Items...> *) {
            return (
                []<typename First, typename... Rest>(TagPath<First, Rest...> *) {
                    return countTag<First, Space>() > 0;
                }(static_cast<typename Items::PathType *>(nullptr))
                && ... && true
            );
        }(static_cast<typename ItemsOf<typename Model::template Policy<PolicyKind::reads>>::type *>(nullptr));
    }
}

template <typename Scalar, ProcessModel Model, Space S>
[[nodiscard]] typename Moments<S, Scalar>::MatrixType
processNoiseOfBlocks(Model const & model, S const & space, Scalar dt)
{
    using Matrix = typename Moments<S, Scalar>::MatrixType;
    auto const out = outputOf<Model>(space);
    auto const partial = model.template noiseOver<Scalar>(out, dt);
    Eigen::Index const n = differenceDimension(space);
    Matrix q = Matrix::Zero(n, n);
    [&]<std::size_t... K>(std::index_sequence<K...>) {
        std::array<BlockRows, sizeof...(K)> const rows{differenceBlockRows<WritePathAt<Model, K>>(space)...};
        std::array<BlockRows, sizeof...(K)> const outRows{differenceBlockRows<WriteNameAt<Model, K>>(out)...};
        for (std::size_t k = 0; k < sizeof...(K); ++k)
        {
            for (std::size_t l = 0; l < sizeof...(K); ++l)
            {
                q.block(rows[k].start, rows[l].start, rows[k].size, rows[l].size) =
                    partial.block(outRows[k].start, outRows[l].start, outRows[k].size, outRows[l].size);
            }
        }
    }(std::make_index_sequence<writeCount<Model>>{});
    return q;
}

template <typename Scalar, ProcessModel Model, Space S>
[[nodiscard]] typename Moments<S, Scalar>::MatrixType
processNoise(Model const & model, S const & space, Scalar dt)
{
    using Matrix = typename Moments<S, Scalar>::MatrixType;
    if constexpr (writesWholeState<Model>)
    {
        return Matrix{model.template noiseOver<Scalar>(space, dt)};
    }
    else
    {
        return processNoiseOfBlocks<Scalar>(model, space, dt);
    }
}

template <typename Call, typename Space>
class ResolvedProcess : public ResolvedReads<typename Call::ModelType, typename Call::ArgumentsType, Space>
{
    using Base = ResolvedReads<typename Call::ModelType, typename Call::ArgumentsType, Space>;
 
public:
 
    using CallType = Call;                          
    using Model = typename Call::ModelType;         
    using Function = typename Model::FunctionType;  
 
    using Base::InputCount;
    using Base::ItemCount;
    using Base::StateCount;
    using typename Base::In;
    using typename Base::Reads;

    template <std::size_t J>
    using Argument = typename Base::template Argument<J>;

    template <typename T>
    using InVector = typename Base::template InVector<T>;

    using WritePaths = typename Model::WritePaths;
 
    static constexpr std::size_t WriteCount = std::tuple_size_v<typename PathsOf<WritePaths>::type>;  

    static constexpr bool Whole = writesWholeState<Model>;

    using Out = OutputSpaceIn<Model, Space>;

    static constexpr bool Affine = affineInEverything<Model>();

    static constexpr bool ReturnsJacobian = returnsStepJacobianIn<Model, Space, double>;
 
private:
 
    using Base::argumentOf;
    using Base::inputItems;
    using Base::m_differenceRows;
    using Base::m_in;
    using Base::m_reads;
    using Base::m_stateDifferenceSize;
    using Base::stateIndexOf;
    using Base::stateItems;
 
public:

    ResolvedProcess(Call const & call, Space const & space)
        : Base{call.arguments(), space}
        , m_call{call}
        , m_out{outputOf<Model>(space)}
        , m_stateSize{valueDimension(space)}
    {
        if constexpr (Whole)
        {
            m_writeRows[0] = m_outRows[0] = BlockRows{0, m_stateSize};
            m_writeDifferenceRows[0] = m_outDifferenceRows[0] = BlockRows{0, differenceDimension(space)};
        }
        else
        {
            [&]<std::size_t... K>(std::index_sequence<K...>) {
                ((m_writeRows[K] = blockRows<WritePathAt<Model, K>>(space),
                  m_writeDifferenceRows[K] = differenceBlockRows<WritePathAt<Model, K>>(space),
                  m_outRows[K] = blockRows<WriteNameAt<Model, K>>(m_out),
                  m_outDifferenceRows[K] = differenceBlockRows<WriteNameAt<Model, K>>(m_out)),
                 ...);
            }(std::make_index_sequence<WriteCount>{});
        }
    }

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

    [[nodiscard]] Out const &
    out() const
    {
        return m_out;
    }

    template <typename Scalar>
    [[nodiscard]] typename Moments<Space, Scalar>::MatrixType
    noise(Scalar dt) const
    {
        using Matrix = typename Moments<Space, Scalar>::MatrixType;
        auto const partial = model().template noiseOver<Scalar>(m_out, dt);
        Matrix q = Matrix::Zero(m_stateDifferenceSize, m_stateDifferenceSize);
        for (std::size_t k = 0; k < WriteCount; ++k)
        {
            for (std::size_t l = 0; l < WriteCount; ++l)
            {
                q.block(
                    m_writeDifferenceRows[k].start,
                    m_writeDifferenceRows[l].start,
                    m_writeDifferenceRows[k].size,
                    m_writeDifferenceRows[l].size
                ) =
                    partial.block(
                        m_outDifferenceRows[k].start,
                        m_outDifferenceRows[l].start,
                        m_outDifferenceRows[k].size,
                        m_outDifferenceRows[l].size
                    );
            }
        }
        return q;
    }

    template <typename State, typename TimeScalar, typename... Inputs>
    [[nodiscard]] auto
    propagate(State const & state, TimeScalar dt, Inputs const &... inputs) const
    {
        static_assert(
            sizeof...(Inputs) == InputCount || sizeof...(Inputs) == 0,
            "the step takes one vector per uncertain input, or none to read each at its mean"
        );
        using T = ScalarOf<State>;
        using Vector = typename Moments<Space, T>::VectorType;
        EMRYS_FILTER_ASSERT(state.size() == m_stateSize && "the state does not fit the space");
        auto const x =
            this->template gatherFromRows<T>(state, Eigen::Index{0}, std::tuple<Inputs const &...>{inputs...});
        auto result = vectorReturned(callForm(model().function(), m_in, m_out, x, dt));
        using Result = std::remove_cvref_t<decltype(result)>;
        Vector y{state};
        scatter(y, valueOf(result));
        if constexpr (isWithJacobian<Result>)
        {
            return withJacobian(std::move(y), expand(result.jacobian));
        }
        else
        {
            return y;
        }
    }

    static constexpr bool BatchedFunction = takesBatches<Function>;

    template <typename S, typename T>
    using Columns =
        MatrixOf<T, valueDimensionAtCompileTime<S>, Eigen::Dynamic, maxValueDimensionAtCompileTime<S>, Eigen::Dynamic>;

    template <typename T, typename Column, typename State, typename... Inputs>
    void
    gatherColumn(Column && column, State const & state, Inputs const &... inputs) const
    {
        static_assert(sizeof...(Inputs) == InputCount, "the step takes one vector per uncertain input");
        column = this->template gatherFromRows<T>(state, Eigen::Index{0}, std::tuple<Inputs const &...>{inputs...});
    }

    template <typename X, typename TimeScalar, typename Y>
    void
    stepBatch(X const & x, TimeScalar dt, Y && y) const
        requires BatchedFunction
    {
        callBatched(model().function(), m_in, m_out, x, y, dt);
    }

    template <typename Y, typename States, typename TimeScalar, typename... Inputs>
    void
    propagateColumns(Y & y, States const & states, TimeScalar dt, Inputs const &... inputs) const
        requires BatchedFunction
    {
        auto const x = gatherStates<Y>(states, inputs...);
        auto values = ValuesOf<Y, typename States::Scalar>{};
        values.resize(valueDimension(m_out), states.cols());
        stepBatch(x, dt, values);
        scatterColumns(y, states, values);
    }

    template <typename Y, typename J, typename States, typename TimeScalar, typename... Inputs>
    void
    propagateColumnsWithAffineJacobians(
        Y & y,
        J & jacobians,
        States const & states,
        TimeScalar dt,
        Inputs const &... inputs
    ) const
        requires BatchedFunction
    {
        using T = typename States::Scalar;
        auto const x = gatherStates<Y>(states, inputs...);
        Eigen::Index const count = states.cols();
        auto values = ValuesOf<Y, T>{};
        values.resize(valueDimension(m_out), count);
        auto wanted = std::array<bool, std::tuple_size_v<InputBlocks<In>>>{};
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((wanted[I] = Base::template itemFactored<I>()), ...);
        }(std::make_index_sequence<ItemCount>{});
        Eigen::Index const l = jacobians.cols() / std::max<Eigen::Index>(count, 1);
        using OutJacobians = MatrixOf<
            T,
            valueDimensionAtCompileTime<DifferenceSpace<Out>>,
            J::ColsAtCompileTime,
            maxValueDimensionAtCompileTime<DifferenceSpace<Out>>,
            J::MaxColsAtCompileTime>;
        auto outJacobians = OutJacobians{};
        outJacobians.resize(differenceDimension(m_out), l * count);
        callBatchedWithJacobians(model().function(), m_in, m_out, x, wanted, values, outJacobians, dt);
        scatterColumns(y, states, values);
        if constexpr (Whole)
        {
            jacobians = outJacobians;
        }
        else
        {
            jacobians.setZero();
            Eigen::Index offset = 0;
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                (
                    [&] {
                        if constexpr (Base::template itemFactored<I>())
                        {
                            if constexpr (Base::inState(I))
                            {
                                BlockRows const rows = this->m_reads.differenceRows(Base::stateIndexOf(I));
                                for (Eigen::Index j = 0; j < count; ++j)
                                {
                                    for (Eigen::Index r = 0; r < rows.size; ++r)
                                    {
                                        jacobians(rows.start + r, j * l + offset + r) = T(1);
                                    }
                                }
                            }
                            offset += this->m_differenceRows[I].size;
                        }
                    }(),
                    ...
                );
            }(std::make_index_sequence<ItemCount>{});
            for (std::size_t k = 0; k < WriteCount; ++k)
            {
                BlockRows const & written = m_writeDifferenceRows[k];
                BlockRows const & rows = m_outDifferenceRows[k];
                jacobians.middleRows(written.start, written.size) = outJacobians.middleRows(rows.start, rows.size);
            }
        }
    }

    template <typename Y, typename States, typename... Inputs>
    [[nodiscard]] auto
    gatherStates(States const & states, Inputs const &... inputs) const
    {
        static_assert(sizeof...(Inputs) == InputCount);
        using T = typename States::Scalar;
        using X = MatrixOf<
            T,
            valueDimensionAtCompileTime<In>,
            Y::ColsAtCompileTime,
            maxValueDimensionAtCompileTime<In>,
            Y::MaxColsAtCompileTime>;
        EMRYS_FILTER_ASSERT(states.rows() == m_stateSize && "the states do not fit the space");
        auto x = X{};
        x.resize(valueDimension(m_in), states.cols());
        for (Eigen::Index j = 0; j < states.cols(); ++j)
        {
            gatherColumn<T>(x.col(j), states.col(j), inputs.col(j)...);
        }
        return x;
    }

    template <typename Y, typename T>
    using ValuesOf = MatrixOf<
        T,
        valueDimensionAtCompileTime<Out>,
        Y::ColsAtCompileTime,
        maxValueDimensionAtCompileTime<Out>,
        Y::MaxColsAtCompileTime>;

    template <typename Y, typename States, typename Values>
    void
    scatterColumns(Y & y, States const & states, Values const & values) const
    {
        y = states;
        for (Eigen::Index j = 0; j < states.cols(); ++j)
        {
            auto column = y.col(j);
            scatter(column, values.col(j));
        }
    }

    template <typename State, typename Value>
    [[nodiscard]] auto
    scattered(State const & state, Value const & value) const
    {
        using Vector = typename Moments<Space, ScalarOf<State>>::VectorType;
        Vector y{state};
        scatter(y, value);
        return y;
    }

    template <typename Scalar>
    [[nodiscard]] typename Moments<Space, Scalar>::MatrixType
    transition(Scalar dt) const
        requires(Affine && ReturnsJacobian && InputCount == 0)
    {
        using Vector = typename Moments<Space, Scalar>::VectorType;
        auto const zero = Vector{Vector::Zero(m_stateSize)};
        return propagate(zero, dt).jacobian;
    }
 
private:

    template <typename Y, typename Value>
    void
    scatter(Y & y, Value const & value) const
    {
        EMRYS_FILTER_ASSERT(
            value.size() == valueDimension(m_out) && "a process model's function returns a vector in its output space"
        );
        if constexpr (Whole)
        {
            y = value;
        }
        else
        {
            scatterBlocks(y, value);
        }
    }

    template <typename Y, typename Value>
    void
    scatterBlocks(Y & y, Value const & value) const
    {
        [&]<std::size_t... K>(std::index_sequence<K...>) {
            ((middleRowsOf<
                  valueDimensionAtCompileTime<PathTag<std::tuple_element_t<K, typename PathsOf<WritePaths>::type>>>>(
                  y,
                  m_writeRows[K].start,
                  m_writeRows[K].size
              ) =
                  middleRowsOf<valueDimensionAtCompileTime<
                      PathTag<std::tuple_element_t<K, typename PathsOf<WritePaths>::type>>>>(
                      value,
                      m_outRows[K].start,
                      m_outRows[K].size
                  )),
             ...);
        }(std::make_index_sequence<WriteCount>{});
    }

    template <typename J>
    [[nodiscard]] auto
    expand(J const & jacobian) const
    {
        static_assert(
            MatrixExpression<J>,
            "a process model's function returns its Jacobian as one matrix over its input's differences, "
            "withJacobian(y, J)"
        );
        EMRYS_FILTER_ASSERT(
            jacobian.rows() == differenceDimension(m_out) && jacobian.cols() == differenceDimension(m_in)
            && "a returned Jacobian has a row per difference coordinate of the output and a column per one of the input"
        );
        using Scalar = typename J::Scalar;
        using StateMatrix = typename Moments<Space, Scalar>::MatrixType;
        StateMatrix f = StateMatrix::Identity(m_stateDifferenceSize, m_stateDifferenceSize);
        for (std::size_t k = 0; k < WriteCount; ++k)
        {
            f.middleRows(m_writeDifferenceRows[k].start, m_writeDifferenceRows[k].size).setZero();
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                ((f.block(
                      m_writeDifferenceRows[k].start,
                      m_reads.differenceRows(I).start,
                      m_writeDifferenceRows[k].size,
                      m_reads.differenceRows(I).size
                  ) =
                      jacobian.block(
                          m_outDifferenceRows[k].start,
                          m_differenceRows[stateItems(I)].start,
                          m_outDifferenceRows[k].size,
                          m_differenceRows[stateItems(I)].size
                      )),
                 ...);
            }(std::make_index_sequence<StateCount>{});
        }
        if constexpr (InputCount == 0)
        {
            return f;
        }
        else
        {
            return [&]<std::size_t... P>(std::index_sequence<P...>) {
                return std::tuple{std::move(f), inputBlock<P>(jacobian)...};
            }(std::make_index_sequence<InputCount>{});
        }
    }

    template <std::size_t P, typename J>
    [[nodiscard]] auto
    inputBlock(J const & jacobian) const
    {
        using Input =
            std::remove_cvref_t<decltype(inputOf(std::declval<Argument<argumentOf(inputItems(P))> const &>()))>;
        using Block = MatrixOf<
            typename J::Scalar,
            valueDimensionAtCompileTime<DifferenceSpace<Space>>,
            Input::DifferenceDimensionAtCompileTime,
            maxValueDimensionAtCompileTime<DifferenceSpace<Space>>,
            Input::MaxDifferenceDimensionAtCompileTime>;
        BlockRows const & columns = m_differenceRows[inputItems(P)];
        Block block = Block::Zero(m_stateDifferenceSize, columns.size);
        for (std::size_t k = 0; k < WriteCount; ++k)
        {
            block.middleRows(m_writeDifferenceRows[k].start, m_writeDifferenceRows[k].size) =
                jacobian.block(m_outDifferenceRows[k].start, columns.start, m_outDifferenceRows[k].size, columns.size);
        }
        return block;
    }
 
    Call const & m_call;                                        
    Out m_out;                                                  
    Eigen::Index m_stateSize;                                   
    std::array<BlockRows, WriteCount> m_writeRows{};            
    std::array<BlockRows, WriteCount> m_writeDifferenceRows{};  
    std::array<BlockRows, WriteCount> m_outRows{};              
    std::array<BlockRows, WriteCount> m_outDifferenceRows{};    
};

template <typename R, typename TimeScalar>
class ProcessStep
{
public:

    static constexpr bool EvaluatesBatches = R::BatchedFunction;

    ProcessStep(R const & resolved, TimeScalar dt)
        : m_resolved{resolved}
        , m_dt{dt}
    {}

    template <typename X, typename... Inputs>
    [[nodiscard]] auto
    operator()(X const & x, Inputs const &... inputs) const
    {
        return m_resolved.propagate(x, m_dt, inputs...);
    }

    template <typename Y, typename States, typename... Inputs>
    void
    evaluateBatch(Y & y, States const & states, Inputs const &... inputs) const
        requires EvaluatesBatches
    {
        m_resolved.propagateColumns(y, states, m_dt, inputs...);
    }

    template <typename Y, typename J, typename States, typename... Inputs>
    void
    evaluateBatchWithAffineJacobians(Y & y, J & jacobians, States const & states, Inputs const &... inputs) const
        requires EvaluatesBatches
    {
        m_resolved.propagateColumnsWithAffineJacobians(y, jacobians, states, m_dt, inputs...);
    }
 
private:
 
    R const & m_resolved;  
    TimeScalar m_dt;       
};

template <typename R, typename TimeScalar>
class ProcessStepValues : public ProcessStep<R, TimeScalar>
{
public:
 
    using ProcessStep<R, TimeScalar>::ProcessStep;

    template <typename X, typename... Inputs>
    [[nodiscard]] auto
    operator()(X const & x, Inputs const &... inputs) const
    {
        return valueOf(ProcessStep<R, TimeScalar>::operator()(x, inputs...));
    }
};

template <typename Space, typename Call>
[[nodiscard]] ResolvedProcess<Call, Space>
resolveProcess(Call const & call, Space const & space)
{
    return ResolvedProcess<Call, Space>{call, space};
}

template <typename... Args>
inline constexpr bool referenceAfterFirst = []<std::size_t... I>(std::index_sequence<I...>) {
    return ((I > 0 && isLinearizedAbout<std::tuple_element_t<I, std::tuple<Args...>>>) || ... || false);
}(std::index_sequence_for<Args...>{});

template <typename F, typename... Args>
[[nodiscard]] decltype(auto)
withReference(F && f, Args const &... args)
{
    static_assert(
        !referenceAfterFirst<Args...>,
        "a prediction's linearizedAbout() reference comes first among the arguments after the step, once"
    );
    if constexpr (leadsWithReference<Args...>)
    {
        return std::forward<F>(f)(args...);
    }
    else
    {
        return std::forward<F>(f)(NoReference{}, args...);
    }
}
 
}  // namespace internal

template <Space Space_, typename Derived, ProcessModel Model, typename... Args>
[[nodiscard]] typename Moments<Space_, typename Derived::Scalar>::VectorType
propagate(
    Space_ const & space,
    Eigen::MatrixBase<Derived> const & x,
    Model const & model,
    TimeStep<typename Derived::Scalar> auto const & dt,
    Args const &... args
)
{
    using Scalar = typename Derived::Scalar;
    EMRYS_FILTER_ASSERT(x.size() == valueDimension(space) && "the state does not fit the space");
    auto const call = internal::modelWith(std::cref(model), args...);
    auto const resolved = internal::resolveProcess(call, space);
    return internal::valueOf(resolved.propagate(x.derived(), internal::toScalarTime<Scalar>(dt)));
}
 
template <Space Space_, typename Model_, typename Method_, RealScalar Scalar_>
class FixedStepModel;

template <typename Model>
inline constexpr bool isFixedStep = false;

template <typename Space_, typename Model_, typename Method_, typename Scalar_>
inline constexpr bool isFixedStep<FixedStepModel<Space_, Model_, Method_, Scalar_>>{true};
 
template <ProcessModel Model_, typename Method_, typename Time_>
class UnboundFixedStep;

template <typename Model>
inline constexpr bool isUnboundStep = false;

template <typename Model_, typename Method_, typename Time_>
inline constexpr bool isUnboundStep<UnboundFixedStep<Model_, Method_, Time_>>{true};

template <Space Space_, RealScalar Scalar_ = double>
struct PredictedJoint
{
    using Scalar = Scalar_;                                
    using SpaceType = Space_;                              
    using GaussianType = GaussianParams<Space_, Scalar_>;  
    using MatrixType = typename GaussianType::MatrixType;  
 
    GaussianType filtered;   
    GaussianType predicted;  

    MatrixType cross;
};
 
namespace internal
{

template <typename Method, typename R>
inline constexpr bool predictsByJacobian =
    std::is_same_v<Method, methods::Linearized> && R::InputCount == 0 && R::ReturnsJacobian;

template <typename Method, typename SpaceType, typename Prior, typename R, typename TimeScalar>
[[nodiscard]] auto
transformStep(Method const & method, SpaceType const & space, Prior const & prior, R const & resolved, TimeScalar dt)
{
    auto const f = ProcessStep{resolved, dt};
    return std::apply(
        [&](auto const &... inputs) {
            if constexpr (isSigmaPointMethod<Method> && R::factorsAffine())
            {
                auto const factored = FactoredAt<R::affineSize(), R::affineMaxSize(), Method>{
                    method,
                    resolved.template affineCoordinates<true>()
                };
                return transformBy(factored, f, space, prior, inputs...);
            }
            else
            {
                return transformBy(method, f, space, prior, inputs...);
            }
        },
        resolved.inputs()
    );
}

template <UncertainInput Prior>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
gaussianOf(Prior const & prior)
{
    using Result = GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>;
    if constexpr (std::is_same_v<Prior, Result>)
    {
        return prior;
    }
    else
    {
        return {
            prior.space(),
            typename Result::VectorType{prior.mean()},
            typename Result::MatrixType{covarianceLowerOf(prior)}
        };
    }
}

template <UncertainInput Prior, typename Form, typename Noise>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
jointOf(
    Prior const & prior,
    OutputWithCross<typename Prior::SpaceType, typename Prior::SpaceType, Form, typename Prior::Scalar> step,
    Noise const & noise
)
{
    using Matrix = typename PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>::MatrixType;
    Matrix cross = step.crossCovariance();
    auto predicted = std::move(step).output();
    predicted.setCovariance(predicted.covarianceLower() + noise);
    return {gaussianOf(prior), std::move(predicted), std::move(cross)};
}

template <UncertainInput Prior, typename Form>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
jointOf(
    Prior const & prior,
    OutputWithCross<typename Prior::SpaceType, typename Prior::SpaceType, Form, typename Prior::Scalar> step
)
{
    using Matrix = typename PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>::MatrixType;
    Matrix cross = step.crossCovariance();
    return {gaussianOf(prior), std::move(step).output(), std::move(cross)};
}

template <typename R, typename Derived, typename TimeScalar>
[[nodiscard]] auto
stepWithJacobian(R const & resolved, Eigen::MatrixBase<Derived> const & x, TimeScalar dt)
{
    return resolved.propagate(x.derived(), dt);
}

template <typename Space, typename Scalar, typename R, typename AddNoise>
void
predictInPlaceByJacobianWithNoise(
    GaussianParams<Space, Scalar> & state,
    R const & resolved,
    Scalar dt,
    AddNoise const & addNoise
)
{
    using Matrix = typename GaussianParams<Space, Scalar>::MatrixType;
    auto const & space = state.space();
    auto step = stepWithJacobian(resolved, state.mean(), dt);
    auto & moments = GaussianAccess::moments(state);
    Matrix & cov = moments.covarianceLower();
    cov.template triangularView<Eigen::StrictlyUpper>() = cov.transpose();
    auto propagated = Matrix{};
    propagated.noalias() = multiply(step.jacobian, cov);
    cov.noalias() = multiply(propagated, step.jacobian.transpose());
    addNoise(cov);
    moments.mean() = std::move(step.value);
    normalizeValues(space, moments.mean());
}

template <typename Space, typename Scalar, typename R>
void
predictInPlaceByJacobian(GaussianParams<Space, Scalar> & state, R const & resolved, Scalar dt)
{
    auto const addNoise = [&](typename GaussianParams<Space, Scalar>::MatrixType & cov) {
        cov += resolved.template noise<Scalar>(dt);
    };
    predictInPlaceByJacobianWithNoise(state, resolved, dt, addNoise);
}

template <typename Method, typename R, typename Space, typename Base, typename Reference>
[[nodiscard]] LinearizationAbout<Space, Space, typename Base::Scalar>
linearizeStepAbout(
    Method const & method,
    R const & resolved,
    Space const & space,
    typename Base::Scalar dt,
    Eigen::MatrixBase<Base> const & base,
    LinearizedAbout<Reference> const & about
)
{
    using Scalar = typename Base::Scalar;
    using Result = LinearizationAbout<Space, Space, Scalar>;
    if constexpr (predictsByJacobian<Method, R>)
    {
        typename Moments<Space, Scalar>::VectorType const point = about.mean();
        auto step = stepWithJacobian(resolved, point, dt);
        Eigen::Index const n = differenceDimension(space);
        auto result = Result{
            space,
            typename Result::OutputVector{std::move(step.value)},
            typename Result::JacobianMatrix{std::move(step.jacobian)},
            Result::ErrorMatrix::Zero(n, n),
            Result::ErrorMatrix::Zero(n, n)
        };
        normalizeValues(space, result.value);
        toChartOf(space, base, point, result);
        return result;
    }
    else
    {
        auto const f = ProcessStep{resolved, dt};
        return std::apply(
            [&](auto const &... inputs) {
                if constexpr (isSigmaPointMethod<Method> && R::factorsAffine())
                {
                    auto const factored = FactoredAt<R::affineSize(), R::affineMaxSize(), Method>{
                        method,
                        resolved.template affineCoordinates<true>()
                    };
                    return linearizeAbout(factored, f, space, space, base, about, inputs...);
                }
                else
                {
                    return linearizeAbout(method, f, space, space, base, about, inputs...);
                }
            },
            resolved.inputs()
        );
    }
}

template <typename Method, UncertainInput Prior, typename R, typename Derived, typename Reference>
[[nodiscard]] OutputWithCross<
    typename Prior::SpaceType,
    typename Prior::SpaceType,
    JacobianForm,
    typename Prior::Scalar>
predictAbout(
    Method const & method,
    Prior const & prior,
    R const & resolved,
    typename Prior::Scalar dt,
    Eigen::MatrixBase<Derived> const & noise,
    LinearizedAbout<Reference> const & about
)
{
    using Space = typename Prior::SpaceType;
    using Scalar = typename Prior::Scalar;
    using Result = OutputWithCross<Space, Space, JacobianForm, Scalar>;
    using Matrix = typename GaussianParams<Space, Scalar>::MatrixType;
    auto const & space = prior.space();
    auto linearization = linearizeStepAbout(method, resolved, space, dt, prior.mean(), about);
    Matrix const cov = symmetricFromLower<Matrix>(covarianceLowerOf(prior));
    typename Result::CrossMatrix cross = multiply(cov, linearization.jacobian.transpose());
    Matrix explained = multiply(linearization.jacobian, cross);
    explained = (Scalar(0.5) * (explained + explained.transpose())).eval();
    // The process noise is one more input the step does not explain.
    Matrix others = linearization.inputs + symmetricFromLower<Matrix>(noise);
    return Result{
        space,
        std::move(linearization.value),
        std::move(linearization.jacobian),
        std::move(cross),
        std::move(explained),
        std::move(others),
        std::move(linearization.unexplained)
    };
}

template <typename Method, UncertainInput Prior, typename R, typename Derived, typename Reference = NoReference>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
predictWithNoise(
    Method const & method,
    Prior const & prior,
    R const & resolved,
    typename Prior::Scalar dt,
    Eigen::MatrixBase<Derived> const & noise,
    Reference const & reference = {}
)
{
    if constexpr (isLinearizedAbout<Reference>)
    {
        return predictAbout(method, prior, resolved, dt, noise, reference).output();
    }
    else
    {
        auto predicted = transformStep(method, prior.space(), prior, resolved, dt);
        predicted.setCovariance(predicted.covarianceLower() + noise);
        return predicted;
    }
}

template <typename Method, UncertainInput Prior, typename R, typename Derived, typename Reference = NoReference>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
predictJointWithNoise(
    Method const & method,
    Prior const & prior,
    R const & resolved,
    typename Prior::Scalar dt,
    Eigen::MatrixBase<Derived> const & noise,
    Reference const & reference = {}
)
{
    if constexpr (isLinearizedAbout<Reference>)
    {
        return jointOf(prior, predictAbout(method, prior, resolved, dt, noise, reference));
    }
    else
    {
        return jointOf(prior, transformStep(method, prior.space(), withCross(prior), resolved, dt), noise);
    }
}

template <typename Method, UncertainInput Prior, typename R, typename Derived, typename Reference>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
predictResolved(
    Method const & method,
    Prior const & prior,
    R const & resolved,
    typename Prior::Scalar dt,
    Eigen::MatrixBase<Derived> const & noise,
    Reference const & reference
)
{
    using Result = GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>;
    if constexpr (!isLinearizedAbout<Reference> && predictsByJacobian<Method, R>)
    {
        // The function gives its Jacobian at the mean, so the Linearized prediction needs no
        // differentiation: x' = f(x), P' = F P F^T + Q.
        auto result = Result{
            prior.space(),
            typename Result::VectorType{prior.mean()},
            typename Result::MatrixType{covarianceLowerOf(prior)}
        };
        predictInPlaceByJacobianWithNoise(result, resolved, dt, [&](typename Result::MatrixType & cov) {
            cov += noise;
        });
        return result;
    }
    else
    {
        return predictWithNoise(method, prior, resolved, dt, noise, reference);
    }
}

template <typename Model, typename Method>
inline constexpr bool predictsBy = true;

template <typename Space_, typename Model_, typename ModelMethod, typename Scalar_, typename Method>
inline constexpr bool predictsBy<FixedStepModel<Space_, Model_, ModelMethod, Scalar_>, Method>{
    std::is_same_v<ModelMethod, Method>
};

template <typename Method>
[[nodiscard]] consteval bool
checkPredictionMethod()
{
    static_assert(
        !isIterated<Method>,
        "Iterated is a measurement update's method: a prediction has no posterior to relinearize at, so give the "
        "filter or the fixed-step model its inner method"
    );
    return true;
}

template <typename Model>
struct DefaultMethodOf
{
    using type = methods::Linearized;  
};

template <typename Space_, typename Model_, typename Method_, typename Scalar_>
struct DefaultMethodOf<FixedStepModel<Space_, Model_, Method_, Scalar_>>
{
    using type = Method_;  
};
 
}  // namespace internal

template <typename Model>
using DefaultMethodFor = typename internal::DefaultMethodOf<Model>::type;

struct ModelsMethod
{};

template <typename Model, typename Method>
using MethodFor = std::conditional_t<isFixedStep<Model>, ModelsMethod, Method>;
 
namespace internal
{

template <typename Model, typename Method>
    requires(!std::is_same_v<Method, ModelsMethod>)
[[nodiscard]] Method
predictionMethod(Model const &, Method method)
{
    return method;
}

template <typename Space_, typename Model_, typename Method_, typename Scalar_>
[[nodiscard]] Method_
predictionMethod(FixedStepModel<Space_, Model_, Method_, Scalar_> const & model, ModelsMethod)
{
    return model.method();
}

template <typename Model>
inline constexpr bool cachesTransition = false;

template <typename Space_, typename Model_, typename Method_, typename Scalar_>
inline constexpr bool cachesTransition<FixedStepModel<Space_, Model_, Method_, Scalar_>> =
    FixedStepModel<Space_, Model_, Method_, Scalar_>::CachesTransition;

struct NoTransition
{};
 
}  // namespace internal

template <typename Method, UncertainInput Prior, ProcessModel Model, typename... Args>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
predict(
    Method const & method,
    Prior const & prior,
    Model const & model,
    TimeStep<typename Prior::Scalar> auto const & time,
    Args const &... args
)
{
    if constexpr (!requires { prior.covarianceLower(); })
    {
        return predict(method, internal::withCovariance(prior), model, time, args...);
    }
    else
    {
        using Scalar = typename Prior::Scalar;
        auto const & space = prior.space();
        Scalar const dt = internal::toScalarTime<Scalar>(time);
        internal::assertPartitions(model, space);
        return internal::withReference(
            [&](auto const & reference, auto const &... inputs) {
                auto const call = internal::modelWith(std::cref(model), inputs...);
                auto const resolved = internal::resolveProcess(call, space);
                return internal::predictResolved(
                    method,
                    prior,
                    resolved,
                    dt,
                    resolved.template noise<Scalar>(dt),
                    reference
                );
            },
            args...
        );
    }
}

template <typename Method, UncertainInput Prior, ProcessModel Model, typename... Args>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
predictJoint(
    Method const & method,
    Prior const & prior,
    Model const & model,
    TimeStep<typename Prior::Scalar> auto const & time,
    Args const &... args
)
{
    if constexpr (!requires { prior.covarianceLower(); })
    {
        return predictJoint(method, internal::withCovariance(prior), model, time, args...);
    }
    else
    {
        using Scalar = typename Prior::Scalar;
        auto const & space = prior.space();
        Scalar const dt = internal::toScalarTime<Scalar>(time);
        internal::assertPartitions(model, space);
        return internal::withReference(
            [&](auto const & reference, auto const &... inputs) {
                auto const call = internal::modelWith(std::cref(model), inputs...);
                auto const resolved = internal::resolveProcess(call, space);
                return internal::predictJointWithNoise(
                    method,
                    prior,
                    resolved,
                    dt,
                    resolved.template noise<Scalar>(dt),
                    reference
                );
            },
            args...
        );
    }
}

template <ProcessModel Model_, typename Method_, typename Time_>
class UnboundFixedStep
{
public:
 
    using ModelType = Model_;    
    using MethodType = Method_;  
    using TimeType = Time_;      

    UnboundFixedStep(Model_ model, Time_ dt, Method_ method)
        : m_model{std::move(model)}
        , m_method{std::move(method)}
        , m_timeStep{std::move(dt)}
    {}

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

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

    [[nodiscard]] Time_ const &
    timeStep() const
    {
        return m_timeStep;
    }

    template <typename Space>
    [[nodiscard]] bool
    partitions(Space const & space) const
    {
        return m_model.partitions(space);
    }

    template <typename Space>
    [[nodiscard]] static consteval internal::StaticPartition
    staticPartition()
    {
        return internal::staticPartition<Model_, Space>();
    }
 
private:
 
    Model_ m_model;                                   
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Method_ m_method;  
    Time_ m_timeStep;                                 
};
 
namespace internal
{

template <typename Space_>
[[nodiscard]] Space_
boundSpace()
{
    static_assert(
        FixedSizeSpace<Space_> && std::default_initializable<Space_>,
        "fixedStep(model, dt) binds to a fixed-size state space; over a space sized at run time, give "
        "FixedStepModel{model, space, dt} the space"
    );
    return Space_{};
}
 
}  // namespace internal

template <Space Space_, typename Model_, typename Method_ = methods::Linearized, RealScalar Scalar_ = double>
class FixedStepModel
{
    static_assert(ProcessModel<Model_>, "a fixed-step model is a process model taken over one step");
    static_assert(internal::checkPredictionMethod<Method_>());
 
public:
 
    using Scalar = Scalar_;                                            
    using SpaceType = Space_;                                          
    using ModelType = Model_;                                          
    using MethodType = Method_;                                        
    using MatrixType = typename Moments<Space_, Scalar_>::MatrixType;  

    static constexpr bool ReadsOnlyState = internal::readsOnlyState<Model_, Space_>();

    static constexpr bool CachesTransition = std::is_same_v<Method_, methods::Linearized> && ReadsOnlyState
                                          && internal::returnsStepJacobianIn<Model_, Space_, Scalar_>
                                          && internal::affineInEverything<Model_>();

    FixedStepModel(Model_ model, Space_ const & space, TimeStep<Scalar_> auto const & dt, Method_ method = {})
        : m_model{std::move(model)}
        , m_method{std::move(method)}
        , m_timeStep{internal::toScalarTime<Scalar_>(dt)}
        , m_noise{internal::processNoise<Scalar_>(m_model, space, m_timeStep)}
        , m_noiseRoot{internal::defaultRoot(m_noise)}
        , m_transition{transitionOf(m_model, space, m_timeStep)}
        , m_transitionInverse{inverseOf(m_transition)}
    {
        internal::assertPartitions(m_model, space);
    }

    template <TimeStep<Scalar_> Time>
    FixedStepModel(UnboundFixedStep<Model_, Method_, Time> const & step)
        : FixedStepModel{step.model(), internal::boundSpace<Space_>(), step.timeStep(), step.method()}
    {}

    template <typename Space>
    [[nodiscard]] static consteval internal::StaticPartition
    staticPartition()
    {
        return internal::staticPartition<Model_, Space>();
    }

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

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

    [[nodiscard]] Scalar_
    timeStep() const
    {
        return m_timeStep;
    }

    [[nodiscard]] MatrixType const &
    noise() const
    {
        return m_noise;
    }

    [[nodiscard]] MatrixType const &
    transition() const
        requires CachesTransition
    {
        return m_transition;
    }

    [[nodiscard]] MatrixType const &
    noiseRoot() const
    {
        return m_noiseRoot;
    }

    [[nodiscard]] std::optional<MatrixType> const &
    transitionInverse() const
        requires CachesTransition
    {
        return m_transitionInverse;
    }

    template <typename Derived>
    [[nodiscard]] typename Moments<Space_, typename Derived::Scalar>::VectorType
    propagate(Space_ const & space, Eigen::MatrixBase<Derived> const & x) const
        requires ReadsOnlyState
    {
        return Emrys::propagate(space, x, m_model, m_timeStep);
    }

    [[nodiscard]] bool
    partitions(Space_ const & space) const
    {
        return m_model.partitions(space) && m_noise.rows() == differenceDimension(space);
    }
 
private:

    [[nodiscard]] static auto
    transitionOf(
        [[maybe_unused]] Model_ const & model,
        [[maybe_unused]] Space_ const & space,
        [[maybe_unused]] Scalar_ dt
    )
    {
        if constexpr (CachesTransition)
        {
            auto const call = internal::modelWith(std::cref(model));
            auto const resolved = internal::resolveProcess(call, space);
            return MatrixType{resolved.template transition<Scalar_>(dt)};
        }
        else
        {
            return internal::NoTransition{};
        }
    }

    [[nodiscard]] static auto
    inverseOf([[maybe_unused]] auto const & transition)
    {
        if constexpr (CachesTransition)
        {
            auto const lu = Eigen::FullPivLU<MatrixType>{transition};
            return lu.isInvertible() ? std::optional<MatrixType>{internal::inverseFrom<MatrixType>(lu)} : std::nullopt;
        }
        else
        {
            return internal::NoTransition{};
        }
    }
 
    Model_ m_model;                                   
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Method_ m_method;  
    Scalar_ m_timeStep;                               
    MatrixType m_noise;                               
    MatrixType m_noiseRoot;                           

    EMRYS_FILTER_NO_UNIQUE_ADDRESS std::conditional_t<CachesTransition, MatrixType, internal::NoTransition>
        m_transition;

    EMRYS_FILTER_NO_UNIQUE_ADDRESS
    std::conditional_t<CachesTransition, std::optional<MatrixType>, internal::NoTransition>
        m_transitionInverse;
};

template <typename Model, typename Space, typename Time>
FixedStepModel(Model, Space, Time) -> FixedStepModel<Space, Model, methods::Linearized, double>;
 
template <typename Model, typename Space, typename Time, typename Method>
FixedStepModel(Model, Space, Time, Method) -> FixedStepModel<Space, Model, Method, double>;


template <ProcessModel Model, typename Time, typename Method = methods::Linearized>
    requires TimeStep<Time, double>
[[nodiscard]] UnboundFixedStep<Model, Method, Time>
fixedStep(Model model, Time dt, Method method = {})
{
    return UnboundFixedStep<Model, Method, Time>{std::move(model), std::move(dt), std::move(method)};
}
 
namespace internal
{

template <typename Model, typename Space, typename Scalar>
struct BoundStepOf
{
    using type = Model;  
};

template <typename Model_, typename Method_, typename Time_, typename Space, typename Scalar>
struct BoundStepOf<UnboundFixedStep<Model_, Method_, Time_>, Space, Scalar>
{
    using type = FixedStepModel<Space, Model_, Method_, Scalar>;  
};

template <typename Model, typename Space, typename Scalar>
using BoundStep = typename BoundStepOf<Model, Space, Scalar>::type;

template <typename Model_, typename Method_, typename Time_>
struct DefaultMethodOf<UnboundFixedStep<Model_, Method_, Time_>>
{
    using type = Method_;  
};
 
}  // namespace internal
 
namespace internal
{

template <UncertainInput Prior, typename Space_, typename Model_, typename Method_, typename Scalar_>
    requires FixedStepModel<Space_, Model_, Method_, Scalar_>::CachesTransition
[[nodiscard]] OutputWithCross<Space_, Space_, JacobianForm, Scalar_>
predictByTransition(Prior const & prior, FixedStepModel<Space_, Model_, Method_, Scalar_> const & model)
{
    using MatrixType = typename FixedStepModel<Space_, Model_, Method_, Scalar_>::MatrixType;
    auto const & space = prior.space();
    auto const & f = model.transition();
    typename Moments<Space_, Scalar_>::VectorType mean{model.propagate(space, prior.mean())};
    normalizeValues(space, mean);
    MatrixType cross = multiply(covarianceLowerOf(prior), f.transpose());
    MatrixType explained = multiply(f, cross);
    explained = (Scalar_(0.5) * (explained + explained.transpose())).eval();
    Eigen::Index const n = explained.rows();
    // The process noise is one more input the transition does not explain.
    return OutputWithCross<Space_, Space_, JacobianForm, Scalar_>{
        space,
        std::move(mean),
        MatrixType{f},
        std::move(cross),
        std::move(explained),
        MatrixType{model.noise()},
        MatrixType::Zero(n, n)
    };
}
 
}  // namespace internal

template <UncertainInput Prior, typename Space_, typename Model_, typename Method_, typename Scalar_, typename... Args>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
predict(Prior const & prior, FixedStepModel<Space_, Model_, Method_, Scalar_> const & model, Args const &... args)
{
    static_assert(std::is_same_v<typename Prior::SpaceType, Space_>, "the model must be over the prior's space");
    static_assert(std::is_same_v<typename Prior::Scalar, Scalar_>, "the model's scalar must be the prior's");
    using Model = FixedStepModel<Space_, Model_, Method_, Scalar_>;
    [[maybe_unused]] auto const & space = prior.space();
    internal::assertPartitions(model, space);
    if constexpr (Model::CachesTransition && sizeof...(Args) == 0)
    {
        return internal::predictByTransition(prior, model).output();
    }
    else
    {
        return internal::withReference(
            [&](auto const & reference, auto const &... inputs) {
                auto const call = internal::modelWith(std::cref(model.model()), inputs...);
                auto const resolved = internal::resolveProcess(call, space);
                return internal::predictResolved(
                    model.method(),
                    prior,
                    resolved,
                    model.timeStep(),
                    model.noise(),
                    reference
                );
            },
            args...
        );
    }
}

template <UncertainInput Prior, typename Space_, typename Model_, typename Method_, typename Scalar_, typename... Args>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
predictJoint(Prior const & prior, FixedStepModel<Space_, Model_, Method_, Scalar_> const & model, Args const &... args)
{
    static_assert(std::is_same_v<typename Prior::SpaceType, Space_>, "the model must be over the prior's space");
    static_assert(std::is_same_v<typename Prior::Scalar, Scalar_>, "the model's scalar must be the prior's");
    using Model = FixedStepModel<Space_, Model_, Method_, Scalar_>;
    [[maybe_unused]] auto const & space = prior.space();
    internal::assertPartitions(model, space);
    if constexpr (Model::CachesTransition && sizeof...(Args) == 0)
    {
        return internal::jointOf(prior, internal::predictByTransition(prior, model));
    }
    else
    {
        return internal::withReference(
            [&](auto const & reference, auto const &... inputs) {
                auto const call = internal::modelWith(std::cref(model.model()), inputs...);
                auto const resolved = internal::resolveProcess(call, space);
                return internal::predictJointWithNoise(
                    model.method(),
                    prior,
                    resolved,
                    model.timeStep(),
                    model.noise(),
                    reference
                );
            },
            args...
        );
    }
}
 
namespace internal
{

template <typename Prior, typename Model_, typename Method_, typename Time_>
[[nodiscard]] FixedStepModel<typename Prior::SpaceType, Model_, Method_, typename Prior::Scalar>
boundToPrior(Prior const & prior, UnboundFixedStep<Model_, Method_, Time_> const & step)
{
    return FixedStepModel<typename Prior::SpaceType, Model_, Method_, typename Prior::Scalar>{
        step.model(),
        prior.space(),
        step.timeStep(),
        step.method()
    };
}
 
}  // namespace internal

template <UncertainInput Prior, typename Model_, typename Method_, typename Time_, typename... Args>
[[nodiscard]] GaussianParams<typename Prior::SpaceType, typename Prior::Scalar>
predict(Prior const & prior, UnboundFixedStep<Model_, Method_, Time_> const & step, Args const &... args)
{
    return Emrys::predict(prior, internal::boundToPrior(prior, step), args...);
}

template <UncertainInput Prior, typename Model_, typename Method_, typename Time_, typename... Args>
[[nodiscard]] PredictedJoint<typename Prior::SpaceType, typename Prior::Scalar>
predictJoint(Prior const & prior, UnboundFixedStep<Model_, Method_, Time_> const & step, Args const &... args)
{
    return Emrys::predictJoint(prior, internal::boundToPrior(prior, step), args...);
}
 
namespace internal
{

template <typename Method, typename Space, typename Scalar, ProcessModel Model, typename... Args>
void
predictInPlace(
    Method const & method,
    GaussianParams<Space, Scalar> & state,
    Model const & model,
    Scalar dt,
    Args const &... args
)
{
    withReference(
        [&](auto const & reference, auto const &... inputs) {
            auto const call = modelWith(std::cref(model), inputs...);
            auto const resolved = resolveProcess(call, state.space());
            using R = std::remove_cvref_t<decltype(resolved)>;
            if constexpr (!isLinearizedAbout<std::remove_cvref_t<decltype(reference)>> && predictsByJacobian<Method, R>)
            {
                assertPartitions(model, state.space());
                predictInPlaceByJacobian(state, resolved, dt);
            }
            else
            {
                state = Emrys::predict(method, state, model, dt, args...);
            }
        },
        args...
    );
}

template <
    typename Space,
    typename Scalar,
    typename Space_,
    typename Model_,
    typename Method_,
    typename Scalar_,
    typename... Args>
void
predictInPlace(
    GaussianParams<Space, Scalar> & state,
    FixedStepModel<Space_, Model_, Method_, Scalar_> const & model,
    Args const &... args
)
{
    using Model = FixedStepModel<Space_, Model_, Method_, Scalar_>;
    static_assert(std::is_same_v<Space, Space_>);
    static_assert(std::is_same_v<Scalar, Scalar_>);
    if constexpr (Model::CachesTransition && sizeof...(Args) == 0)
    {
        state = Emrys::predict(state, model);
    }
    else
    {
        withReference(
            [&](auto const & reference, auto const &... inputs) {
                auto const call = modelWith(std::cref(model.model()), inputs...);
                auto const resolved = resolveProcess(call, state.space());
                using R = std::remove_cvref_t<decltype(resolved)>;
                if constexpr (!isLinearizedAbout<std::remove_cvref_t<decltype(reference)>>
                              && predictsByJacobian<Method_, R>)
                {
                    assertPartitions(model, state.space());
                    predictInPlaceByJacobianWithNoise(state, resolved, model.timeStep(), [&](auto & cov) {
                        cov += model.noise();
                    });
                }
                else
                {
                    state = Emrys::predict(state, model, args...);
                }
            },
            args...
        );
    }
}
 
}  // namespace internal
 
}  // namespace Emrys