ResolvedReads.h

Go to the documentation of this file.

#pragma once

 
#include <array>
#include <cstddef>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <unsupported/Eigen/AutoDiff>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Models/MeasurementModel.h>
#include <EmrysFilter/Models/Policies.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Space, typename Reads>
class ReadBlocks;

template <typename Space, typename... Paths>
class ReadBlocks<Space, TagPaths<Paths...>>
{
    static_assert(sizeof...(Paths) > 0, "a model reads at least one block of the state");

    using First = std::tuple_element_t<0, std::tuple<Paths...>>;
 
public:
 
    static constexpr std::size_t Count = sizeof...(Paths);  

    using JointSpace = std::conditional_t<Count == 1, PathTag<First>, Composite<PathTag<Paths>...>>;
 
    static constexpr int Size = valueDimensionAtCompileTime<JointSpace>;        
    static constexpr int MaxSize = maxValueDimensionAtCompileTime<JointSpace>;  
    static constexpr int DifferenceSize =
        valueDimensionAtCompileTime<DifferenceSpace<JointSpace>>;  
    static constexpr int DifferenceMaxSize =
        maxValueDimensionAtCompileTime<DifferenceSpace<JointSpace>>;           
    using Indices = IndexVector<Size, MaxSize>;                                
    using DifferenceIndices = IndexVector<DifferenceSize, DifferenceMaxSize>;  
    using Selection = Coordinates<Indices, DifferenceIndices>;                 

    explicit ReadBlocks(Space const & space)
        : m_rows{blockRows<Paths>(space)...}
        , m_differenceRows{differenceBlockRows<Paths>(space)...}
        , m_space{makeSpace(space)}
    {
        Eigen::Index offset = 0;
        Eigen::Index differenceOffset = 0;
        for (std::size_t i = 0; i < Count; ++i)
        {
            m_jointStarts[i] = offset;
            m_jointDifferenceStarts[i] = differenceOffset;
            offset += m_rows[i].size;
            differenceOffset += m_differenceRows[i].size;
        }
        m_indices = concatenate<Indices>(m_rows);
        m_differenceIndices = concatenate<DifferenceIndices>(m_differenceRows);
        EMRYS_FILTER_ASSERT(distinctIndices(m_indices) && "a model must not read a coordinate twice");
    }

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

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

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

    [[nodiscard]] Selection
    selection() const
    {
        return {m_indices, m_differenceIndices};
    }

    [[nodiscard]] BlockRows const &
    rows(std::size_t k) const
    {
        return m_rows[k];
    }

    [[nodiscard]] BlockRows const &
    differenceRows(std::size_t k) const
    {
        return m_differenceRows[k];
    }

    [[nodiscard]] Eigen::Index
    jointStart(std::size_t k) const
    {
        return m_jointStarts[k];
    }

    [[nodiscard]] Eigen::Index
    jointDifferenceStart(std::size_t k) const
    {
        return m_jointDifferenceStarts[k];
    }
 
private:

    [[nodiscard]] static JointSpace
    makeSpace(Space const & space)
    {
        if constexpr (Count == 1)
        {
            return pathBlock<First>(space);
        }
        else
        {
            return JointSpace{pathBlock<Paths>(space)...};
        }
    }

    template <typename Result>
    [[nodiscard]] static Result
    concatenate(std::array<BlockRows, Count> const & blocks)
    {
        Eigen::Index size = 0;
        for (BlockRows const & rows : blocks)
        {
            size += rows.size;
        }
        auto result = Result{};
        result.resize(size);
        Eigen::Index k = 0;
        for (BlockRows const & rows : blocks)
        {
            for (Eigen::Index j = 0; j < rows.size; ++j)
            {
                result(k++) = rows.start + j;
            }
        }
        return result;
    }
 
    std::array<BlockRows, Count> m_rows;                        
    std::array<BlockRows, Count> m_differenceRows;              
    std::array<Eigen::Index, Count> m_jointStarts{};            
    std::array<Eigen::Index, Count> m_jointDifferenceStarts{};  
    Indices m_indices;                                          
    DifferenceIndices m_differenceIndices;                      
    JointSpace m_space;                                         
};

template <typename Space>
class ReadBlocks<Space, TagPaths<WholeSpace>>
{
public:
 
    static constexpr std::size_t Count = 1;  

    using JointSpace = Space;
 
    static constexpr int Size = valueDimensionAtCompileTime<Space>;        
    static constexpr int MaxSize = maxValueDimensionAtCompileTime<Space>;  
    static constexpr int DifferenceSize =
        valueDimensionAtCompileTime<DifferenceSpace<Space>>;  
    static constexpr int DifferenceMaxSize = maxValueDimensionAtCompileTime<DifferenceSpace<Space>>;  
    using Indices = IndexVector<Size, MaxSize>;                                
    using DifferenceIndices = IndexVector<DifferenceSize, DifferenceMaxSize>;  
    using Selection = Coordinates<Indices, DifferenceIndices>;                 

    explicit ReadBlocks(Space const & space)
        : m_rows{0, valueDimension(space)}
        , m_differenceRows{0, differenceDimension(space)}
        , m_indices{Indices::LinSpaced(m_rows.size, 0, m_rows.size - 1)}
        , m_differenceIndices{DifferenceIndices::LinSpaced(m_differenceRows.size, 0, m_differenceRows.size - 1)}
        , m_space{space}
    {}

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

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

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

    [[nodiscard]] Selection
    selection() const
    {
        return {m_indices, m_differenceIndices};
    }

    [[nodiscard]] BlockRows const &
    rows(std::size_t) const
    {
        return m_rows;
    }

    [[nodiscard]] BlockRows const &
    differenceRows(std::size_t) const
    {
        return m_differenceRows;
    }

    [[nodiscard]] Eigen::Index
    jointStart(std::size_t) const
    {
        return 0;
    }

    [[nodiscard]] Eigen::Index
    jointDifferenceStart(std::size_t) const
    {
        return 0;
    }
 
private:
 
    BlockRows m_rows;                       
    BlockRows m_differenceRows;             
    Indices m_indices;                      
    DifferenceIndices m_differenceIndices;  
    Space m_space;                          
};

template <typename Space>
struct ReadsWhole
{
    static constexpr std::size_t Count = 1;  

    using InputSpace = Space;

    template <std::size_t I>
    using Item = ReadItem<WholeSpace, Space>;
};

template <typename Model>
inline constexpr bool readsWholeState = !Model::template has<PolicyKind::reads>;

template <typename Model>
struct WholeStateRefusal
{
    static_assert(
        !readsWholeState<Model>,
        "a model without reads() and writes() reads and writes the whole state, so it is not placed by on() or "
        "onEach() or composed by compose(): give it reads() and writes() to name its blocks"
    );
 
    static constexpr bool value = true;  
};

template <typename Model>
inline constexpr bool wholeStateRefused = WholeStateRefusal<Model>::value;

template <typename Model, typename Space>
using ReadsIn =
    std::conditional_t<readsWholeState<Model>, ReadsWhole<Space>, typename Model::template Policy<PolicyKind::reads>>;

template <typename Tuple>
struct TagPathsOfTuple;

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


template <typename Policy, typename Name>
inline constexpr bool affineNames = false;

template <typename... Names, typename Name>
inline constexpr bool affineNames<AffineIn<Names...>, Name> = (std::is_same_v<Names, Name> || ...);


template <typename Name, typename Block>
[[nodiscard]] Name
nameOf(Block const & block)
{
    if constexpr (std::is_same_v<Name, Block>)
    {
        return block;
    }
    else if constexpr (FixedSizeSpace<Name>)
    {
        return Name{};
    }
    else
    {
        auto name = Name{};
        fitValueDimension(name, valueDimension(block));
        return name;
    }
}

template <typename T>
struct PlainScalarOf
{
    using type = T;  
};

template <typename Derivatives>
struct PlainScalarOf<Eigen::AutoDiffScalar<Derivatives>>
{
    using type = Eigen::AutoDiffScalar<typename std::remove_cvref_t<Derivatives>::PlainObject>;
};


template <typename R>
[[nodiscard]] auto
vectorReturned(R && returned)
{
    using Plain = std::remove_cvref_t<R>;
    if constexpr (!MatrixExpression<Plain> && !isWithJacobian<Plain>)
    {
        static_assert(
            std::is_same_v<Plain, typename PlainScalarOf<Plain>::type>,
            "a function that returns a scalar returns it as its input's scalar type, T{...}: an expression of "
            "automatic-differentiation scalars refers to temporaries of the function; or return a vector of one "
            "coordinate"
        );
        return Eigen::Matrix<Plain, 1, 1>{std::forward<R>(returned)};
    }
    else
    {
        return std::forward<R>(returned);
    }
}

template <typename Model>
[[nodiscard]] consteval bool
affineInEverything()
{
    if constexpr (readsWholeState<Model>)
    {
        return requires { requires Model::FunctionType::IsAffine; };
    }
    else
    {
        using Affine = typename Model::template Policy<PolicyKind::affineIn>;
        return []<typename... Items>(Reads<Items...> const *) {
            return (affineNames<Affine, typename Items::NameType> && ... && true);
        }(static_cast<typename Model::template Policy<PolicyKind::reads> const *>(nullptr));
    }
}

template <typename Model, typename Arguments, typename Space>
class ResolvedReads
{
public:
 
    using ModelType = Model;  
    using SpaceType = Space;  

    static constexpr bool Whole = readsWholeState<Model>;

    using ReadsPolicy = ReadsIn<Model, Space>;
 
    static constexpr std::size_t ItemCount = ReadsPolicy::Count;  

    template <std::size_t I>
    using Item = typename ReadsPolicy::template Item<I>;

    using In = typename ReadsPolicy::InputSpace;

    using ArgumentsType = Arguments;
 
    static constexpr std::size_t ArgumentCount = std::tuple_size_v<Arguments>;  

    template <std::size_t J>
    using Argument = std::tuple_element_t<J, Arguments>;
 
protected:

    template <std::size_t I>
    [[nodiscard]] static consteval bool
    itemInState()
    {
        if constexpr (Whole)
        {
            return true;
        }
        else
        {
            return []<typename First, typename... Rest>(TagPath<First, Rest...> *) {
                return countTag<First, Space>() > 0;
            }(static_cast<typename Item<I>::PathType *>(nullptr));
        }
    }

    template <std::size_t I>
    [[nodiscard]] static consteval std::size_t
    argumentOfItem()
    {
        if constexpr (itemInState<I>())
        {
            return ArgumentCount;
        }
        else
        {
            static_assert(
                Item<I>::SingleTag,
                "an item of reads() that names no block of the state names an input of the call by one tag, so a "
                "path of several tags must be in the state"
            );
            using Tag = PathTag<typename Item<I>::PathType>;
            constexpr auto matches = []<std::size_t... J>(std::index_sequence<J...>) {
                return std::array<bool, ArgumentCount + 1>{std::is_same_v<InputTag<Argument<J>>, Tag>..., false};
            }(std::make_index_sequence<ArgumentCount>{});
            for (std::size_t j = 0; j < ArgumentCount; ++j)
            {
                if (matches[j])
                {
                    return j;
                }
            }
            return ArgumentCount;
        }
    }
 
    // These are consteval functions of an index rather than static constexpr arrays: MSVC reads a
    // class template's constexpr array members in a class-scope static_assert before it has
    // initialized them ("read of an uninitialized symbol"), and a function has nothing to read.

    [[nodiscard]] static consteval bool
    inState(std::size_t i)
    {
        return []<std::size_t... I>(std::size_t item, std::index_sequence<I...>) {
            return ((I == item && itemInState<I>()) || ... || false);
        }(i, std::make_index_sequence<ItemCount>{});
    }

    [[nodiscard]] static consteval std::size_t
    argumentOf(std::size_t i)
    {
        return []<std::size_t... I>(std::size_t item, std::index_sequence<I...>) {
            std::size_t found = ArgumentCount;
            ((I == item ? (found = argumentOfItem<I>(), 0) : 0), ...);
            return found;
        }(i, std::make_index_sequence<ItemCount>{});
    }
 
    static_assert(
        []<std::size_t... I>(std::index_sequence<I...>) {
            return ((itemInState<I>() || argumentOfItem<I>() < ArgumentCount) && ...);
        }(std::make_index_sequence<ItemCount>{}),
        "the model reads an input that the state lacks, which the call must give as an uncertain input over the "
        "tag's space or as sure<Tag>(value)"
    );

    template <std::size_t J>
    [[nodiscard]] static consteval bool
    argumentRead()
    {
        if constexpr (std::is_void_v<InputTag<Argument<J>>>)
        {
            return true;
        }
        else
        {
            return []<std::size_t... I>(std::index_sequence<I...>) {
                return ((argumentOfItem<I>() == J) || ... || false);
            }(std::make_index_sequence<ItemCount>{});
        }
    }
 
    static_assert(
        []<std::size_t... J>(std::index_sequence<J...>) {
            return (argumentRead<J>() && ... && true);
        }(std::make_index_sequence<ArgumentCount>{}),
        "the call was given an input the model does not read"
    );

    [[nodiscard]] static consteval bool
    uncertainArgument(std::size_t j)
    {
        return []<std::size_t... J>(std::size_t argument, std::index_sequence<J...>) {
            return ((J == argument && argumentKind<Argument<J>> == ArgumentKind::uncertain) || ... || false);
        }(j, std::make_index_sequence<ArgumentCount>{});
    }

    [[nodiscard]] static consteval bool
    argumentIsUncertain(std::size_t i)
    {
        return !inState(i) && uncertainArgument(argumentOf(i));
    }

    [[nodiscard]] static consteval std::size_t
    countItems(bool state)
    {
        std::size_t count = 0;
        for (std::size_t i = 0; i < ItemCount; ++i)
        {
            if (state ? inState(i) : argumentIsUncertain(i))
            {
                ++count;
            }
        }
        return count;
    }
 
public:
 
    static constexpr std::size_t StateCount = countItems(true);   
    static constexpr std::size_t InputCount = countItems(false);  
 
protected:

    [[nodiscard]] static consteval std::size_t
    itemWhere(bool state, std::size_t n)
    {
        std::size_t seen = 0;
        for (std::size_t i = 0; i < ItemCount; ++i)
        {
            if (state ? inState(i) : argumentIsUncertain(i))
            {
                if (seen == n)
                {
                    return i;
                }
                ++seen;
            }
        }
        return ItemCount;
    }

    [[nodiscard]] static consteval std::size_t
    stateItems(std::size_t k)
    {
        return itemWhere(true, k);
    }

    [[nodiscard]] static consteval std::size_t
    inputItems(std::size_t p)
    {
        return itemWhere(false, p);
    }

    [[nodiscard]] static consteval std::size_t
    stateIndexOf(std::size_t i)
    {
        for (std::size_t k = 0; k < StateCount; ++k)
        {
            if (stateItems(k) == i)
            {
                return k;
            }
        }
        return StateCount;
    }

    [[nodiscard]] static consteval std::size_t
    inputIndexOf(std::size_t i)
    {
        for (std::size_t p = 0; p < InputCount; ++p)
        {
            if (inputItems(p) == i)
            {
                return p;
            }
        }
        return InputCount;
    }
 
public:

    using StatePaths = typename TagPathsOfTuple<decltype([]<std::size_t... K>(std::index_sequence<K...>) {
        return std::tuple<typename Item<stateItems(K)>::PathType...>{};
    }(std::make_index_sequence<StateCount>{}))>::type;

    using Reads = ReadBlocks<Space, StatePaths>;

    template <typename T>
    using InVector = typename Moments<In, T>::VectorType;

    ResolvedReads(Arguments const & arguments, Space const & space)
        : m_arguments{arguments}
        , m_reads{space}
        , m_in{makeIn(arguments, space)}
        , m_stateDifferenceSize{differenceDimension(space)}
    {
        if constexpr (Whole)
        {
            m_valueRows[0] = BlockRows{0, valueDimension(m_in)};
            m_differenceRows[0] = BlockRows{0, differenceDimension(m_in)};
        }
        else
        {
            static_assert(checkReadNames<ReadsPolicy, Space>());
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                ((m_valueRows[I] = blockRows<typename Item<I>::NameType>(m_in),
                  m_differenceRows[I] = differenceBlockRows<typename Item<I>::NameType>(m_in)),
                 ...);
            }(std::make_index_sequence<ItemCount>{});
        }
    }

    [[nodiscard]] Arguments const &
    arguments() const
    {
        return m_arguments;
    }

    [[nodiscard]] Reads const &
    reads() const
    {
        return m_reads;
    }

    [[nodiscard]] In const &
    in() const
    {
        return m_in;
    }

    [[nodiscard]] Eigen::Index
    stateDifferenceSize() const
    {
        return m_stateDifferenceSize;
    }

    [[nodiscard]] BlockRows const &
    valueRows(std::size_t i) const
    {
        return m_valueRows[i];
    }

    [[nodiscard]] BlockRows const &
    differenceRows(std::size_t i) const
    {
        return m_differenceRows[i];
    }

    [[nodiscard]] auto
    inputs() const
    {
        return [&]<std::size_t... P>(std::index_sequence<P...>) {
            return std::tuple<Argument<argumentOf(inputItems(P))> const &...>{
                std::get<argumentOf(inputItems(P))>(m_arguments)...
            };
        }(std::make_index_sequence<InputCount>{});
    }

    [[nodiscard]] static consteval bool
    inputsAffine()
    {
        using Affine = typename Model::template Policy<PolicyKind::affineIn>;
        return []<std::size_t... P>(std::index_sequence<P...>) {
            return (affineNames<Affine, typename Item<inputItems(P)>::NameType> && ... && true);
        }(std::make_index_sequence<InputCount>{});
    }

    template <std::size_t I>
    [[nodiscard]] static consteval bool
    itemFactored()
    {
        using Affine = typename Model::template Policy<PolicyKind::affineIn>;
        using Name = typename Item<I>::NameType;
        return affineNames<Affine, Name> && isLinear<Name> && (inState(I) || argumentIsUncertain(I));
    }

    [[nodiscard]] static consteval int
    affineSize()
    {
        return []<std::size_t... I>(std::index_sequence<I...>) {
            return sizeSum<
                (itemFactored<I>() ? valueDimensionAtCompileTime<DifferenceSpace<typename Item<I>::NameType>> : 0)...>;
        }(std::make_index_sequence<ItemCount>{});
    }

    [[nodiscard]] static consteval int
    affineMaxSize()
    {
        return []<std::size_t... I>(std::index_sequence<I...>) {
            return sizeSum<(
                itemFactored<I>() ? maxValueDimensionAtCompileTime<DifferenceSpace<typename Item<I>::NameType>> : 0
            )...>;
        }(std::make_index_sequence<ItemCount>{});
    }

    [[nodiscard]] static consteval bool
    factorsAffine()
    {
        return []<std::size_t... I>(std::index_sequence<I...>) {
            return (itemFactored<I>() || ... || false);
        }(std::make_index_sequence<ItemCount>{});
    }

    template <bool WholeState>
    [[nodiscard]] IndexVector<affineSize(), affineMaxSize()>
    affineCoordinates() const
    {
        std::array<Eigen::Index, InputCount + 1> inputStarts{};
        inputStarts[0] = WholeState ? m_stateDifferenceSize : m_reads.differenceIndices().size();
        std::apply(
            [&](auto const &... input) {
                std::size_t p = 0;
                ((inputStarts[p + 1] = inputStarts[p] + inputOf(input).differenceDimension(), ++p), ...);
            },
            inputs()
        );
        Eigen::Index count = 0;
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((count += itemFactored<I>() ? m_differenceRows[I].size : 0), ...);
        }(std::make_index_sequence<ItemCount>{});
        auto result = IndexVector<affineSize(), affineMaxSize()>{};
        result.resize(count);
        Eigen::Index next = 0;
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (
                [&] {
                    if constexpr (itemFactored<I>())
                    {
                        Eigen::Index start = 0;
                        if constexpr (inState(I))
                        {
                            start = WholeState ? m_reads.differenceRows(stateIndexOf(I)).start
                                               : m_reads.jointDifferenceStart(stateIndexOf(I));
                        }
                        else
                        {
                            start = inputStarts[inputIndexOf(I)];
                        }
                        for (Eigen::Index r = 0; r < m_differenceRows[I].size; ++r)
                        {
                            result(next++) = start + r;
                        }
                    }
                }(),
                ...
            );
        }(std::make_index_sequence<ItemCount>{});
        return result;
    }

    template <typename T, typename Joint, typename Inputs>
    [[nodiscard]] InVector<T>
    gather(Joint const & joint, Inputs const & inputs) const
    {
        auto x = InVector<T>{};
        x.resize(valueDimension(m_in));
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (gatherItem<I, T>(x, joint, inputs), ...);
        }(std::make_index_sequence<ItemCount>{});
        return x;
    }

    template <typename T, typename Rows, typename Inputs>
    [[nodiscard]] InVector<T>
    gatherFromRows(Rows const & rows, Eigen::Index start, Inputs const & inputs) const
    {
        auto x = InVector<T>{};
        x.resize(valueDimension(m_in));
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (gatherItemFromRows<I, T>(x, rows, start, inputs), ...);
        }(std::make_index_sequence<ItemCount>{});
        return x;
    }

    template <typename T, typename X, typename RowsOf>
    void
    gatherColumnsFromRows(X & x, RowsOf const & rowsOf, Eigen::Index start) const
    {
        EMRYS_FILTER_ASSERT(
            x.rows() == valueDimension(m_in) && "the batch has a row per value coordinate of the input"
        );
        for (Eigen::Index j = 0; j < x.cols(); ++j)
        {
            auto column = x.col(j);
            auto const & rows = rowsOf(j);
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                (gatherItemFromRows<I, T>(column, rows, start, std::tuple<>{}), ...);
            }(std::make_index_sequence<ItemCount>{});
        }
    }
 
protected:

    [[nodiscard]] static In
    makeIn([[maybe_unused]] Arguments const & arguments, Space const & space)
    {
        if constexpr (Whole)
        {
            return space;
        }
        else
        {
            return [&]<std::size_t... I>(std::index_sequence<I...>) {
                return In{nameOf<typename Item<I>::NameType>(blockOf<I>(arguments, space))...};
            }(std::make_index_sequence<ItemCount>{});
        }
    }

    template <std::size_t I>
    [[nodiscard]] static auto
    blockOf([[maybe_unused]] Arguments const & arguments, [[maybe_unused]] Space const & space)
    {
        if constexpr (inState(I))
        {
            return pathSubspace<typename Item<I>::PathType>(space);
        }
        else
        {
            auto const & argument = std::get<argumentOf(I)>(arguments);
            if constexpr (isSure<std::remove_cvref_t<decltype(argument)>>)
            {
                return argument.space;
            }
            else
            {
                return inputOf(argument).space();
            }
        }
    }

    template <std::size_t I, typename T, typename Joint, typename Inputs>
    void
    gatherItem(InVector<T> & x, Joint const & joint, Inputs const & inputs) const
    {
        constexpr int Rows = valueDimensionAtCompileTime<typename Item<I>::NameType>;
        auto target = middleRowsOf<Rows>(x, m_valueRows[I].start, m_valueRows[I].size);
        if constexpr (inState(I))
        {
            target = middleRowsOf<Rows>(joint, m_reads.jointStart(stateIndexOf(I)), m_valueRows[I].size);
        }
        else
        {
            gatherArgument<I, T>(target, inputs);
        }
    }

    template <std::size_t I, typename T, typename X, typename Rows, typename Inputs>
    void
    gatherItemFromRows(X & x, Rows const & rows, Eigen::Index start, Inputs const & inputs) const
    {
        constexpr int Size = valueDimensionAtCompileTime<typename Item<I>::NameType>;
        auto target = middleRowsOf<Size>(x, m_valueRows[I].start, m_valueRows[I].size);
        if constexpr (inState(I))
        {
            BlockRows const & block = m_reads.rows(stateIndexOf(I));
            EMRYS_FILTER_ASSERT(
                block.start >= start && block.start + block.size <= start + rows.size()
                && "the rows must hold every block read"
            );
            target = middleRowsOf<Size>(rows, block.start - start, block.size).template cast<T>();
        }
        else
        {
            gatherArgument<I, T>(target, inputs);
        }
    }

    template <std::size_t I, typename T, typename Target, typename Inputs>
    void
    gatherArgument(Target & target, [[maybe_unused]] Inputs const & inputs) const
    {
        auto const & argument = std::get<argumentOf(I)>(m_arguments);
        if constexpr (!argumentIsUncertain(I))
        {
            target = argument.value.template cast<T>();
        }
        else if constexpr (std::tuple_size_v<Inputs> > 0)
        {
            target = std::get<inputIndexOf(I)>(inputs);
        }
        else
        {
            target = inputOf(argument).mean().template cast<T>();
        }
    }
 
    Arguments const & m_arguments;                        
    Reads m_reads;                                        
    In m_in;                                              
    Eigen::Index m_stateDifferenceSize;                   
    std::array<BlockRows, ItemCount> m_valueRows{};       
    std::array<BlockRows, ItemCount> m_differenceRows{};  
};
 
}  // namespace internal
 
}  // namespace Emrys