TransformArgs.h

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <array>
#include <cstddef>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <unsupported/Eigen/AutoDiff>
 
#include <EmrysFilter/Core/Forwarding.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Multiply.h>
#include <EmrysFilter/Distributions/Gaussian.h>
#include <EmrysFilter/Models/Linearization.h>
#include <EmrysFilter/Models/Methods.h>
#include <EmrysFilter/Spaces/Coordinates.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 T>
inline constexpr bool isWithJacobian = false;

template <typename Value, typename Jacobian>
inline constexpr bool isWithJacobian<WithJacobian<Value, Jacobian>>{true};


template <typename Returned>
[[nodiscard]] decltype(auto)
valueOf(Returned && returned)
{
    if constexpr (isWithJacobian<std::remove_cvref_t<Returned>>)
    {
        return (std::forward<Returned>(returned).value);
    }
    else
    {
        return std::forward<Returned>(returned);
    }
}

template <typename Blocks>
struct ArgumentJacobians
{
    Blocks blocks;  
};

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

template <typename Blocks>
inline constexpr bool isArgumentJacobians<ArgumentJacobians<Blocks>>{true};


template <std::size_t K, typename Jacobian>
[[nodiscard]] auto const &
jacobianBlock(Jacobian const & jacobian)
{
    if constexpr (requires { std::tuple_size<Jacobian>::value; })
    {
        return std::get<K>(jacobian);
    }
    else
    {
        static_assert(K == 0, "a function with more than one uncertain input must return a tuple of Jacobian blocks");
        return jacobian;
    }
}

template <int AffineSize, int AffineMaxSize, typename Points>
struct FactoredAt
{
    Points points;                                  
    IndexVector<AffineSize, AffineMaxSize> affine;  
};

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

template <int AffineSize, int AffineMaxSize, typename Points>
inline constexpr bool isFactoredAt<FactoredAt<AffineSize, AffineMaxSize, Points>>{true};

template <int... Sizes>
inline constexpr int sizeSum = (false || ... || (Sizes == Eigen::Dynamic)) ? Eigen::Dynamic : (0 + ... + Sizes);

template <int A, int B>
inline constexpr int sizeProduct = A == Eigen::Dynamic || B == Eigen::Dynamic ? Eigen::Dynamic : A * B;
 
}  // namespace internal

template <typename Input>
struct CrossMarked
{
    Input const & input;  
};

template <typename Reference>
struct LinearizedAbout
{
    using ReferenceType = Reference;

    static constexpr bool IsDistribution = UncertainInput<std::remove_cvref_t<Reference>>;
 
    Reference reference;  

    [[nodiscard]] decltype(auto)
    mean() const
    {
        if constexpr (IsDistribution)
        {
            return reference.mean();
        }
        else
        {
            return (reference);
        }
    }
};
 
namespace internal
{

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

template <typename Reference>
inline constexpr bool isLinearizedAbout<LinearizedAbout<Reference>>{true};


template <typename... Args>
inline constexpr bool leadsWithReference = false;

template <typename First, typename... Rest>
inline constexpr bool leadsWithReference<First, Rest...> = isLinearizedAbout<First>;


struct NoReference
{};

template <typename T>
concept IsSpace = Space<T>;

template <typename Arg>
struct ArgTraits
{
    using Input = Arg;

    using Scalar = void;
 
    static constexpr bool isUncertain = false;   
    static constexpr bool hasCross = false;      
    static constexpr int size = 0;               
    static constexpr int maxSize = 0;            
    static constexpr int differenceSize = 0;     
    static constexpr int differenceMaxSize = 0;  
};

template <UncertainInput Arg>
struct ArgTraits<Arg>
{
    using Input = Arg;

    using Scalar = typename Arg::Scalar;
 
    static constexpr bool isUncertain = true;                                           
    static constexpr bool hasCross = false;                                             
    static constexpr int size = Arg::ValueDimensionAtCompileTime;                       
    static constexpr int maxSize = Arg::MaxValueDimensionAtCompileTime;                 
    static constexpr int differenceSize = Arg::DifferenceDimensionAtCompileTime;        
    static constexpr int differenceMaxSize = Arg::MaxDifferenceDimensionAtCompileTime;  
};

template <typename In>
struct ArgTraits<CrossMarked<In>> : ArgTraits<In>
{
    static constexpr bool hasCross = true;  
};

template <typename... Args>
inline constexpr std::size_t crossCount = (std::size_t{ArgTraits<Args>::hasCross} + ... + 0);

template <typename Arg>
[[nodiscard]] auto const &
argInput(Arg const & arg)
{
    if constexpr (ArgTraits<Arg>::hasCross)
    {
        return arg.input;
    }
    else
    {
        return arg;
    }
}

template <typename First, typename... Rest>
struct FirstNonVoidOf
{
    using type = std::conditional_t<std::is_void_v<First>, typename FirstNonVoidOf<Rest...>::type, First>;
};

template <typename First>
struct FirstNonVoidOf<First>
{
    using type = First;
};

template <typename... T>
using FirstNonVoid = typename FirstNonVoidOf<T...>::type;

template <typename... Args>
struct TransformArgs
{
    static constexpr std::size_t count = sizeof...(Args);                                 
    static constexpr std::array<bool, count> uncertain{ArgTraits<Args>::isUncertain...};  

    static constexpr std::size_t uncertainCount = (std::size_t{ArgTraits<Args>::isUncertain} + ... + 0);
 
    static_assert(uncertainCount > 0, "transform() needs at least one uncertain input");
    static_assert(crossCount<Args...> <= 1, "withCross() marks at most one input");

    using Scalar = FirstNonVoid<typename ArgTraits<Args>::Scalar...>;
 
    static_assert(
        ((std::is_void_v<typename ArgTraits<Args>::Scalar> || std::is_same_v<typename ArgTraits<Args>::Scalar, Scalar>)
         && ...),
        "every uncertain input must have the same scalar type"
    );

    static constexpr int size = sizeSum<ArgTraits<Args>::size...>;
 
    static constexpr int maxSize = sizeSum<ArgTraits<Args>::maxSize...>;  

    static constexpr int differenceSize = sizeSum<ArgTraits<Args>::differenceSize...>;

    static constexpr int differenceMaxSize = sizeSum<ArgTraits<Args>::differenceMaxSize...>;

    static constexpr std::array<std::size_t, count> ordinal{[] {
        auto result = std::array<std::size_t, count>{};
        std::size_t next = 0;
        for (std::size_t i = 0; i < count; ++i)
        {
            result[i] = uncertain[i] ? next++ : 0;
        }
        return result;
    }()};

    static constexpr std::size_t crossPosition = [] {
        constexpr auto marked = std::array<bool, count>{ArgTraits<Args>::hasCross...};
        for (std::size_t i = 0; i < count; ++i)
        {
            if (marked[i])
            {
                return i;
            }
        }
        return count;
    }();
};

template <typename F, typename OutputSpace, typename... Args>
class Transform
{
public:

    using ArgInfo = TransformArgs<Args...>;

    using Scalar = typename ArgInfo::Scalar;

    using Index = Eigen::Index;

    using OutputSpaceType = OutputSpace;
 
    static constexpr int InputSize = ArgInfo::size;                      
    static constexpr int InputMaxSize = ArgInfo::maxSize;                
    static constexpr int InputDifferenceSize = ArgInfo::differenceSize;  
    static constexpr int InputDifferenceMaxSize = ArgInfo::differenceMaxSize;          
    static constexpr int OutputSize = valueDimensionAtCompileTime<OutputSpace>;        
    static constexpr int OutputMaxSize = maxValueDimensionAtCompileTime<OutputSpace>;  

    static constexpr int OutputDifferenceSize = valueDimensionAtCompileTime<DifferenceSpace<OutputSpace>>;
 
    static constexpr int OutputDifferenceMaxSize =
        maxValueDimensionAtCompileTime<DifferenceSpace<OutputSpace>>;  

    template <typename T>
    using InputVector = Eigen::Matrix<T, InputSize, 1, Eigen::ColMajor, InputMaxSize, 1>;

    template <typename T>
    using InputDifferenceVector = Eigen::Matrix<T, InputDifferenceSize, 1, Eigen::ColMajor, InputDifferenceMaxSize, 1>;

    template <typename T>
    using OutputVector = Eigen::Matrix<T, OutputSize, 1, Eigen::ColMajor, OutputMaxSize, 1>;

    template <typename T>
    using OutputDifferenceVector =
        Eigen::Matrix<T, OutputDifferenceSize, 1, Eigen::ColMajor, OutputDifferenceMaxSize, 1>;

    using InputMatrix = Eigen::Matrix<
        Scalar,
        InputDifferenceSize,
        InputDifferenceSize,
        Eigen::ColMajor,
        InputDifferenceMaxSize,
        InputDifferenceMaxSize>;

    using OutputMatrix = Eigen::Matrix<
        Scalar,
        OutputDifferenceSize,
        OutputDifferenceSize,
        Eigen::ColMajor,
        OutputDifferenceMaxSize,
        OutputDifferenceMaxSize>;

    using JacobianMatrix =
        MatrixOf<Scalar, OutputDifferenceSize, InputDifferenceSize, OutputDifferenceMaxSize, InputDifferenceMaxSize>;

    using CrossMatrix =
        MatrixOf<Scalar, InputDifferenceSize, OutputDifferenceSize, InputDifferenceMaxSize, OutputDifferenceMaxSize>;

    template <int Points = InputDifferenceSize, int MaxPoints = InputDifferenceMaxSize>
    using OutputSamples = MatrixOf<Scalar, OutputSize, Points, OutputMaxSize, MaxPoints>;

    template <int Points = InputDifferenceSize, int MaxPoints = InputDifferenceMaxSize>
    using OutputDeviations = MatrixOf<Scalar, OutputDifferenceSize, Points, OutputDifferenceMaxSize, MaxPoints>;

    struct RootResult
    {
        OutputVector<Scalar> mean;  
        InputMatrix root;           
        CrossMatrix whitened;       
        OutputMatrix unexplained;   
    };

    Transform(F & f, OutputSpace const & outputSpace, Args const &... args)
        : m_f{f}
        , m_outputSpace{outputSpace}
        , m_args{args...}
    {
        forEachUncertain([this]<std::size_t I>(auto const & input) {
            m_offsets[I] = m_inputSize;
            m_differenceOffsets[I] = m_inputDifferenceSize;
            m_inputSize += input.valueDimension();
            m_inputDifferenceSize += input.differenceDimension();
        });
        // A resizable output space with one resizable block takes the size f returns.
        if constexpr (OutputSize == Eigen::Dynamic && requires(Index n) { m_outputSpace.resize(n); })
        {
            Index const outputSize = OutputVector<Scalar>{valueOf(call<Scalar>(m_f, mean()))}.size();
            if (valueDimension(m_outputSpace) != outputSize)
            {
                m_outputSpace.resize(outputSize);
            }
        }
        else if constexpr (OutputSize == Eigen::Dynamic)
        {
            // Several resizable blocks: f's output can't say how to split, so the space must fit it.
            EMRYS_FILTER_ASSERT(
                OutputVector<Scalar>{valueOf(call<Scalar>(m_f, mean()))}.size() == valueDimension(m_outputSpace)
                && "the function's output does not fit the output space"
            );
        }
    }

    [[nodiscard]] OutputSpace const &
    outputSpace() const
    {
        return m_outputSpace;
    }

    [[nodiscard]] Index
    inputSize() const
    {
        return m_inputSize;
    }

    [[nodiscard]] Index
    inputDifferenceSize() const
    {
        return m_inputDifferenceSize;
    }

    [[nodiscard]] Index
    outputDifferenceSize() const
    {
        return differenceDimension(m_outputSpace);
    }

    [[nodiscard]] InputVector<Scalar>
    mean() const
    {
        InputVector<Scalar> result = InputVector<Scalar>::Zero(m_inputSize);
        forEachUncertain([&]<std::size_t I>(auto const & input) {
            result.segment(m_offsets[I], input.valueDimension()) = input.mean();
        });
        return result;
    }

    [[nodiscard]] InputMatrix
    covarianceLower() const
    {
        InputMatrix result = InputMatrix::Zero(m_inputDifferenceSize, m_inputDifferenceSize);
        forEachUncertain([&]<std::size_t I>(auto const & input) {
            Index const d = input.differenceDimension();
            result.block(m_differenceOffsets[I], m_differenceOffsets[I], d, d) = symmetricCov(input);
        });
        return result;
    }

    template <typename A = roots::DefaultAxes>
    [[nodiscard]] InputMatrix
    root(A const & axes = {}) const
    {
        internal::refuseFailingAxes<A>();
        InputMatrix result = InputMatrix::Zero(m_inputDifferenceSize, m_inputDifferenceSize);
        forEachUncertain([&]<std::size_t I>(auto const & input) {
            Index const d = input.differenceDimension();
            auto const take = [&](auto const & entry) {
                if constexpr (requires { input.root(); })
                {
                    return internal::rootValue(internal::rootFromFactor(entry, input.root()));
                }
                else
                {
                    return internal::rootValue(entry.root(covarianceLowerOf(input)));
                }
            };
            constexpr bool spread =
                ArgInfo::crossPosition < ArgInfo::count ? I == ArgInfo::crossPosition : ArgInfo::ordinal[I] == 0;
            if constexpr (spread)
            {
                result.block(m_differenceOffsets[I], m_differenceOffsets[I], d, d) = take(axes);
            }
            else
            {
                result.block(m_differenceOffsets[I], m_differenceOffsets[I], d, d) = take(roots::DefaultAxes{});
            }
        });
        return result;
    }

    template <typename X, typename D>
    [[nodiscard]] InputVector<typename X::Scalar>
    plusInputs(Eigen::MatrixBase<X> const & x, Eigen::MatrixBase<D> const & d) const
    {
        InputVector<typename X::Scalar> result = x;
        forEachUncertain([&]<std::size_t I>(auto const & input) {
            using InputSpace = std::remove_cvref_t<decltype(input.space())>;
            constexpr int Values = internal::rowsAtCompileTime<false, InputSpace>;
            constexpr int Differences = internal::rowsAtCompileTime<true, InputSpace>;
            internal::middleRowsOf<Values>(result, m_offsets[I], input.valueDimension()) = Emrys::plus(
                input.space(),
                internal::middleRowsOf<Values>(x, m_offsets[I], input.valueDimension()),
                internal::middleRowsOf<Differences>(d, m_differenceOffsets[I], input.differenceDimension())
            );
        });
        return result;
    }

    template <typename X, typename D>
    [[nodiscard]] InputVector<typename X::Scalar>
    stepInputs(Eigen::MatrixBase<X> const & x, Eigen::MatrixBase<D> const & d) const
    {
        InputVector<typename X::Scalar> result = x;
        forEachUncertain([&]<std::size_t I>(auto const & input) {
            using InputSpace = std::remove_cvref_t<decltype(input.space())>;
            constexpr int Values = internal::rowsAtCompileTime<false, InputSpace>;
            constexpr int Differences = internal::rowsAtCompileTime<true, InputSpace>;
            auto values = internal::middleRowsOf<Values>(result, m_offsets[I], input.valueDimension());
            auto const step =
                internal::middleRowsOf<Differences>(d, m_differenceOffsets[I], input.differenceDimension());
            if constexpr (HasPlus<InputSpace>)
            {
                values = Emrys::plus(
                    input.space(),
                    internal::middleRowsOf<Values>(x, m_offsets[I], input.valueDimension()),
                    step
                );
            }
            else
            {
                values += step;
            }
        });
        return result;
    }

    template <typename T>
    [[nodiscard]] OutputVector<T>
    evaluate(InputVector<T> const & x) const
    {
        OutputVector<T> y = valueOf(call<T>(m_f, x));
        EMRYS_FILTER_ASSERT(y.size() == valueDimension(m_outputSpace) && "the output does not fit the output space");
        return y;
    }

    static constexpr bool EvaluatesBatches = requires { requires std::remove_cvref_t<F>::EvaluatesBatches; };

    template <typename T, int Points, int MaxPoints>
    using InputSamples = MatrixOf<T, InputSize, Points, InputMaxSize, MaxPoints>;

    template <typename Outputs, typename PointAt>
    void
    evaluateEach(Index count, PointAt const & pointAt, Outputs & outputs) const
    {
        using T = typename Outputs::Scalar;
        outputs.resize(valueDimension(m_outputSpace), count);
        if constexpr (EvaluatesBatches)
        {
            auto x = InputSamples<T, Outputs::ColsAtCompileTime, Outputs::MaxColsAtCompileTime>{};
            x.resize(m_inputSize, count);
            for (Index j = 0; j < count; ++j)
            {
                x.col(j) = pointAt(j);
            }
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                m_f.evaluateBatch(outputs, argumentColumns<I>(x)...);
            }(std::index_sequence_for<Args...>{});
        }
        else
        {
            for (Index j = 0; j < count; ++j)
            {
                outputs.col(j) = evaluate<T>(pointAt(j));
            }
        }
    }

    template <typename Outputs, typename Jacobians, typename PointAt>
    void
    evaluateEachWithAffineJacobians(
        Index count,
        PointAt const & pointAt,
        Outputs & outputs,
        Jacobians & jacobians
    ) const
        requires EvaluatesBatches
    {
        outputs.resize(valueDimension(m_outputSpace), count);
        auto x = InputSamples<Scalar, Outputs::ColsAtCompileTime, Outputs::MaxColsAtCompileTime>{};
        x.resize(m_inputSize, count);
        for (Index j = 0; j < count; ++j)
        {
            x.col(j) = pointAt(j);
        }
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            m_f.evaluateBatchWithAffineJacobians(outputs, jacobians, argumentColumns<I>(x)...);
        }(std::index_sequence_for<Args...>{});
    }

    template <typename T>
    [[nodiscard]] decltype(auto)
    returned(InputVector<T> const & x) const
    {
        return call<T>(m_f, x);
    }

    template <typename T, typename G>
    [[nodiscard]] decltype(auto)
    call(G & g, InputVector<T> const & x) const
    {
        return [&]<std::size_t... I>(std::index_sequence<I...>) -> decltype(auto) {
            return callWithScalar<T>(g, argument<T, I>(x)...);
        }(std::index_sequence_for<Args...>{});
    }

    template <typename Fn>
    void
    forEachUncertain(Fn && fn) const
    {
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (visitIfUncertain<I>(fn), ...);
        }(std::index_sequence_for<Args...>{});
    }

    [[nodiscard]] Index
    offset(std::size_t i) const
    {
        return m_offsets[i];
    }

    [[nodiscard]] Index
    differenceOffset(std::size_t i) const
    {
        return m_differenceOffsets[i];
    }

    template <typename MethodResult>
    [[nodiscard]] auto
    assemble(MethodResult const & result) const
    {
        constexpr bool Root = std::is_same_v<MethodResult, RootResult>;
        Index const m = outputDifferenceSize();
        if constexpr (ArgInfo::crossPosition == ArgInfo::count)
        {
            OutputMatrix cov = [&] {
                if constexpr (Root)
                {
                    OutputMatrix parts = result.unexplained;
                    parts.noalias() += multiply(result.whitened.transpose(), result.whitened);
                    return parts;
                }
                else
                {
                    CrossMatrix const cross = multiply(covarianceLower(), result.jacobian.transpose());
                    return OutputMatrix{multiply(result.jacobian, cross)};
                }
            }();
            cov = (Scalar(0.5) * (cov + cov.transpose())).eval();
            return GaussianParams<OutputSpace, Scalar>{m_outputSpace, result.mean, std::move(cov)};
        }
        else
        {
            constexpr std::size_t position = ArgInfo::crossPosition;
            using InputSpace = typename ArgTraits<ArgAt<position>>::Input::SpaceType;
            using Form = std::conditional_t<Root, RootForm, JacobianForm>;
            using WithCrossResult = OutputWithCross<InputSpace, OutputSpace, Form, Scalar>;
            Index const offset = m_differenceOffsets[position];
            Index const rows = argInput(std::get<position>(m_args)).differenceDimension();
            auto const symmetric = [](OutputMatrix & a) {
                a = (Scalar(0.5) * (a + a.transpose())).eval();
            };
            // N: each other input's share, from its own columns of J or rows of Psi.
            OutputMatrix inputs = OutputMatrix::Zero(m, m);
            forEachUncertain([&]<std::size_t I>(auto const & input) {
                if constexpr (I != position)
                {
                    using Input = std::remove_cvref_t<decltype(input)>;
                    using Shared = MatrixOf<
                        Scalar,
                        Input::DifferenceDimensionAtCompileTime,
                        OutputDifferenceSize,
                        Input::MaxDifferenceDimensionAtCompileTime,
                        OutputDifferenceMaxSize>;
                    Index const at = m_differenceOffsets[I];
                    Index const d = input.differenceDimension();
                    if constexpr (Root)
                    {
                        Shared const block = result.whitened.middleRows(at, d);
                        inputs.noalias() += multiply(block.transpose(), block);
                    }
                    else
                    {
                        Shared const shared =
                            multiply(symmetricCov(input), result.jacobian.middleCols(at, d).transpose());
                        inputs.noalias() += multiply(result.jacobian.middleCols(at, d), shared);
                    }
                }
            });
            symmetric(inputs);
            if constexpr (Root)
            {
                typename WithCrossResult::CrossMatrix whitened = result.whitened.middleRows(offset, rows);
                OutputMatrix explained = multiply(whitened.transpose(), whitened);
                symmetric(explained);
                return WithCrossResult{
                    m_outputSpace,
                    result.mean,
                    typename WithCrossResult::RootMatrix{result.root.block(offset, offset, rows, rows)},
                    std::move(whitened),
                    std::move(explained),
                    std::move(inputs),
                    result.unexplained
                };
            }
            else
            {
                typename WithCrossResult::JacobianMatrix jacobian = result.jacobian.middleCols(offset, rows);
                typename WithCrossResult::CrossMatrix cross =
                    multiply(symmetricCov(argInput(std::get<position>(m_args))), jacobian.transpose());
                OutputMatrix explained = multiply(jacobian, cross);
                symmetric(explained);
                return WithCrossResult{
                    m_outputSpace,
                    result.mean,
                    std::move(jacobian),
                    std::move(cross),
                    std::move(explained),
                    std::move(inputs),
                    OutputMatrix::Zero(m, m)
                };
            }
        }
    }
 
private:

    template <std::size_t I>
    using ArgAt = std::tuple_element_t<I, std::tuple<Args...>>;

    template <typename Input>
    [[nodiscard]] static auto
    symmetricCov(Input const & input)
    {
        using Block = Eigen::Matrix<
            Scalar,
            Input::DifferenceDimensionAtCompileTime,
            Input::DifferenceDimensionAtCompileTime,
            Eigen::ColMajor,
            Input::MaxDifferenceDimensionAtCompileTime,
            Input::MaxDifferenceDimensionAtCompileTime>;
        Block cov = covarianceLowerOf(input);
        cov.template triangularView<Eigen::StrictlyUpper>() = cov.transpose();
        return cov;
    }

    template <std::size_t I, typename Fn>
    void
    visitIfUncertain(Fn & fn) const
    {
        if constexpr (ArgInfo::uncertain[I])
        {
            fn.template operator()<I>(argInput(std::get<I>(m_args)));
        }
    }

    template <typename T, std::size_t I>
    [[nodiscard]] decltype(auto)
    argument(InputVector<T> const & x) const
    {
        auto const & arg = std::get<I>(m_args);
        if constexpr (ArgInfo::uncertain[I])
        {
            using Input = typename ArgTraits<ArgAt<I>>::Input;
            using Block = Eigen::Matrix<
                T,
                Input::ValueDimensionAtCompileTime,
                1,
                Eigen::ColMajor,
                Input::MaxValueDimensionAtCompileTime,
                1>;
            return Block{x.segment(m_offsets[I], argInput(arg).valueDimension())};
        }
        else
        {
            return arg;
        }
    }

    template <std::size_t I, typename X>
    [[nodiscard]] decltype(auto)
    argumentColumns(X const & x) const
    {
        auto const & arg = std::get<I>(m_args);
        if constexpr (ArgInfo::uncertain[I])
        {
            using Input = typename ArgTraits<ArgAt<I>>::Input;
            return internal::middleRowsOf<Input::ValueDimensionAtCompileTime>(
                x,
                m_offsets[I],
                argInput(arg).valueDimension()
            );
        }
        else
        {
            return (arg);
        }
    }

    template <typename T, std::size_t... I>
    static auto returnedAt(std::index_sequence<I...>) -> decltype(std::declval<WithScalar<T, F> const &>()(
        std::declval<Transform const &>().template argument<T, I>(std::declval<InputVector<T> const &>())...
    ));
 
public:

    template <typename T>
    using Returned = decltype(returnedAt<T>(std::index_sequence_for<Args...>{}));
 
private:
 
    F & m_f;                                                   
    OutputSpace m_outputSpace;                                 
    std::tuple<Args const &...> m_args;                        
    std::array<Index, sizeof...(Args)> m_offsets{};            
    std::array<Index, sizeof...(Args)> m_differenceOffsets{};  
    Index m_inputSize = 0;                                     
    Index m_inputDifferenceSize = 0;                           
};
 
}  // namespace internal
 
}  // namespace Emrys