Transform.h

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <expected>
#include <limits>
#include <optional>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <Eigen/Cholesky>
#include <Eigen/Eigenvalues>
#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/Core/Solve.h>
#include <EmrysFilter/Distributions/Gaussian.h>
#include <EmrysFilter/Distributions/SampleSet.h>
#include <EmrysFilter/Distributions/Sampling.h>
#include <EmrysFilter/Models/Batched.h>
#include <EmrysFilter/Models/Condition.h>
#include <EmrysFilter/Models/Linearization.h>
#include <EmrysFilter/Models/MeasurementModel.h>
#include <EmrysFilter/Models/Methods.h>
#include <EmrysFilter/Models/ResolvedReads.h>
#include <EmrysFilter/Models/TransformArgs.h>
#include <EmrysFilter/Spaces/Coordinates.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Reset.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{

 
namespace internal
{

template <typename Problem>
struct LinearizationAt
{
    typename Problem::template OutputVector<typename Problem::Scalar> mean;  
    typename Problem::JacobianMatrix jacobian;  
};

template <typename Problem>
void
normalizeLinearization(Problem const & problem, LinearizationAt<Problem> & result)
{
    normalizeValues(problem.outputSpace(), result.mean);
}

template <typename Problem>
inline constexpr bool returnsJacobian = isWithJacobian<
    std::remove_cvref_t<decltype(std::declval<Problem const &>().template returned<typename Problem::Scalar>(
        std::declval<typename Problem::template InputVector<typename Problem::Scalar> const &>()
    ))>>;

template <typename Problem>
using DualOf = Eigen::AutoDiffScalar<Eigen::Matrix<
    typename Problem::Scalar,
    Problem::InputDifferenceSize,
    1,
    Eigen::ColMajor,
    Problem::InputDifferenceMaxSize,
    1>>;

template <typename Problem, typename Dual = DualOf<Problem>>
inline constexpr bool returnsDuals = requires {
    requires std::same_as<typename std::remove_cvref_t<typename Problem::template Returned<Dual>>::Scalar, Dual>;
};

template <typename Problem>
[[nodiscard]] LinearizationAt<Problem>
linearizeReturned(Problem const & problem)
{
    using Scalar = typename Problem::Scalar;
    using JacobianMatrix = typename Problem::JacobianMatrix;
    auto returned = problem.template returned<Scalar>(problem.mean());
    using Blocks = decltype(returned.jacobian);
    auto result = LinearizationAt<Problem>{
        std::move(returned.value),
        JacobianMatrix(problem.outputDifferenceSize(), problem.inputDifferenceSize())
    };
    EMRYS_FILTER_ASSERT(
        result.mean.size() == valueDimension(problem.outputSpace()) && "the output does not fit the output space"
    );
    problem.forEachUncertain([&]<std::size_t I>(auto const & input) {
        auto const & block = [&]() -> auto const & {
            if constexpr (isArgumentJacobians<Blocks>)
            {
                return std::get<I>(returned.jacobian.blocks);
            }
            else
            {
                return jacobianBlock<Problem::ArgInfo::ordinal[I]>(returned.jacobian);
            }
        }();
        EMRYS_FILTER_ASSERT(block.rows() == problem.outputDifferenceSize() && block.cols() == input.differenceDimension() &&
                        "a Jacobian block must have a row per output difference coordinate and a column per input "
                        "difference coordinate");
        result.jacobian.middleCols(problem.differenceOffset(I), input.differenceDimension()) = block;
    });
    normalizeLinearization(problem, result);
    return result;
}

template <typename Problem>
[[nodiscard]] LinearizationAt<Problem>
linearizeByAutoDiff(Problem const & problem)
{
    using Scalar = typename Problem::Scalar;
    using Dual = DualOf<Problem>;
    using OutputSpace = typename Problem::OutputSpaceType;
    Eigen::Index const n = problem.inputDifferenceSize();
    typename Problem::template InputVector<Scalar> const mu{problem.mean()};
    typename Problem::template InputDifferenceVector<Dual> delta;
    delta.resize(n);
    for (Eigen::Index i = 0; i < n; ++i)
    {
        delta(i) = Dual{Scalar(0), int(n), int(i)};
    }
    typename Problem::template InputVector<Dual> const x{problem.stepInputs(mu.template cast<Dual>(), delta)};
    typename Problem::template OutputVector<Dual> const y{problem.template evaluate<Dual>(x)};
 
    auto result = LinearizationAt<Problem>{
        Problem::template OutputVector<Scalar>::Zero(y.size()),
        Problem::JacobianMatrix::Zero(problem.outputDifferenceSize(), n)
    };
    for (Eigen::Index r = 0; r < y.size(); ++r)
    {
        result.mean(r) = y(r).value();
    }
    auto const deviation = [&] {
        if constexpr (HasMinus<OutputSpace>)
        {
            return typename Problem::template OutputDifferenceVector<Dual>{
                Emrys::minus(problem.outputSpace(), y, result.mean.template cast<Dual>())
            };
        }
        else
        {
            return y;
        }
    }();
    EMRYS_FILTER_ASSERT(deviation.size() == result.jacobian.rows());
    for (Eigen::Index r = 0; r < deviation.size(); ++r)
    {
        // An output that does not depend on the inputs may carry no derivatives at all.
        if (deviation(r).derivatives().size() > 0)
        {
            result.jacobian.row(r) = deviation(r).derivatives().transpose();
        }
    }
    normalizeLinearization(problem, result);
    return result;
}

template <typename Problem>
[[nodiscard]] LinearizationAt<Problem>
linearizeAt(methods::Linearized const &, Problem const & problem)
{
    if constexpr (returnsJacobian<Problem>)
    {
        return linearizeReturned(problem);
    }
    else
    {
        static_assert(
            returnsDuals<Problem>,
            "methods::linearized() differentiates f automatically, so f must take its uncertain inputs as const auto& "
            "and build its result with their scalar, which it may take as its template parameter, "
            "[]<typename T>(auto const & x), or return its Jacobian with withJacobian()"
        );
        if constexpr (returnsDuals<Problem>)
        {
            return linearizeByAutoDiff(problem);
        }
        else
        {
            std::unreachable();
        }
    }
}

template <typename Method, typename Problem>
    requires isLinearMethod<Method>
[[nodiscard]] LinearizationAt<Problem>
transformResult(Method const & method, Problem const & problem)
{
    return linearizeAt(method, problem);
}

template <typename Method, typename Scalar, typename Problem>
[[nodiscard]] typename Problem::RootResult
transformResult(methods::Iterated<Method, Scalar> const & /*method*/, Problem const & /*problem*/)
{
    static_assert(
        sizeof(Method) == 0,
        "Iterated is a measurement update's method: a transform or a prediction has no posterior to relinearize at, "
        "and an update given a block does not support it yet; use its inner method there"
    );
    return {};
}

template <typename Scalar>
struct SigmaWeights
{
    Scalar scale;      
    Scalar center;     
    Scalar centerCov;  
    Scalar point;      
};

template <typename Scalar, UnscentedRule Rule, int N>
[[nodiscard]] consteval SigmaWeights<Scalar>
unscentedWeights()
{
    if constexpr (N == 0)
    {
        return SigmaWeights<Scalar>{Scalar(0), Scalar(1), Scalar(1), Scalar(0)};
    }
    else
    {
        constexpr double n = N;
        constexpr UnscentedParams<double> params = Rule::template params<double>(N);
        constexpr double alphaSquared = params.alpha * params.alpha;
        constexpr double lambda = alphaSquared * (n + params.kappa) - n;
        static_assert(n + lambda > 0, "the sigma-point spread must be positive");
        constexpr double center = lambda / (n + lambda);
        return SigmaWeights<Scalar>{
            Scalar(compileTimeSqrt(n + lambda)),
            Scalar(center),
            Scalar(center + 1 - alphaSquared + params.beta),
            Scalar(1 / (2 * (n + lambda)))
        };
    }
}

template <typename Scalar, int N, UnscentedRule Rule>
[[nodiscard]] SigmaWeights<Scalar>
sigmaWeights(Rule const & /*method*/, Eigen::Index n)
{
    if constexpr (N != Eigen::Dynamic && takesDoubleConstants<Scalar>)
    {
        EMRYS_FILTER_ASSERT(n == N);
        return unscentedWeights<Scalar, Rule, N>();
    }
    else
    {
        if (n == 0)
        {
            return {Scalar(0), Scalar(1), Scalar(1), Scalar(0)};
        }
        UnscentedParams<Scalar> const params = Rule::template params<Scalar>(n);
        Scalar const alphaSquared = params.alpha * params.alpha;
        Scalar const lambda = alphaSquared * (Scalar(n) + params.kappa) - Scalar(n);
        EMRYS_FILTER_ASSERT(Scalar(n) + lambda > Scalar(0) && "the sigma-point spread must be positive");
        Scalar const center = lambda / (Scalar(n) + lambda);
        return {
            Eigen::numext::sqrt(Scalar(n) + lambda),
            center,
            center + Scalar(1) - alphaSquared + params.beta,
            Scalar(1) / (Scalar(2) * (Scalar(n) + lambda))
        };
    }
}

template <typename Problem, int Rows, int MaxRows>
using WhitenedOf =
    MatrixOf<typename Problem::Scalar, Rows, Problem::OutputDifferenceSize, MaxRows, Problem::OutputDifferenceMaxSize>;

template <typename Problem, typename Deviations>
struct SigmaMoments
{
    typename Problem::template OutputVector<typename Problem::Scalar> mean;  

    WhitenedOf<Problem, Deviations::ColsAtCompileTime, Deviations::MaxColsAtCompileTime> whitened;
 
    typename Problem::OutputMatrix unexplained;  
    Deviations plus;                             
    Deviations minus;                            
    typename Problem::template OutputDifferenceVector<typename Problem::Scalar> center;  
};

template <typename Problem, typename Samples>
[[nodiscard]] auto
sigmaMoments(
    Problem const & problem,
    SigmaWeights<typename Problem::Scalar> const & weights,
    typename Problem::template OutputVector<typename Problem::Scalar> const & center,
    Samples const & plusPoints,
    Samples const & minusPoints
)
{
    using Scalar = typename Problem::Scalar;
    using OutputVector = typename Problem::template OutputVector<Scalar>;
    using OutputDeviations =
        typename Problem::template OutputDeviations<Samples::ColsAtCompileTime, Samples::MaxColsAtCompileTime>;
    auto const & space = problem.outputSpace();
    OutputDeviations plusDeviations = Emrys::minus(space, plusPoints, center);
    OutputDeviations minusDeviations = Emrys::minus(space, minusPoints, center);
    OutputVector mean =
        Emrys::plus(space, center, weights.point * (plusDeviations.rowwise().sum() + minusDeviations.rowwise().sum()));
 
    typename Problem::template OutputDifferenceVector<Scalar> const centerDeviation = Emrys::minus(space, center, mean);
    plusDeviations = Emrys::minus(space, plusPoints, mean);
    minusDeviations = Emrys::minus(space, minusPoints, mean);
 
    using Result = SigmaMoments<Problem, OutputDeviations>;
    OutputDeviations const residual = Scalar(0.5) * (plusDeviations + minusDeviations);
    decltype(Result::whitened) whitened =
        (weights.point * weights.scale) * (plusDeviations - minusDeviations).transpose();
    typename Problem::OutputMatrix unexplained = weights.centerCov * centerDeviation * centerDeviation.transpose();
    if constexpr (Samples::ColsAtCompileTime != 0)
    {
        // A product over no points has an inner size of zero, which MSVC's Eigen refuses to block.
        unexplained.noalias() += (Scalar(2) * weights.point) * multiply(residual, residual.transpose());
    }
    unexplained = (Scalar(0.5) * (unexplained + unexplained.transpose())).eval();
    return Result{
        std::move(mean),
        std::move(whitened),
        std::move(unexplained),
        std::move(plusDeviations),
        std::move(minusDeviations),
        centerDeviation
    };
}

template <
    typename Problem,
    int Points_ =
        Problem::InputDifferenceSize == Eigen::Dynamic ? Eigen::Dynamic : 2 * Problem::InputDifferenceSize + 1,
    int MaxPoints_ =
        Problem::InputDifferenceMaxSize == Eigen::Dynamic ? Eigen::Dynamic : 2 * Problem::InputDifferenceMaxSize + 1>
struct SigmaRoot
{
    static constexpr int Points = Points_;

    static constexpr int MaxPoints = MaxPoints_;

    using RootMatrix = typename Problem::template OutputDeviations<Points, MaxPoints>;
 
    typename Problem::template OutputVector<typename Problem::Scalar> mean;  
    RootMatrix root;  
    bool centerNegative = false;  
};

template <typename Problem>
struct SigmaOutputs
{
    typename Problem::template OutputVector<typename Problem::Scalar> center;  
    typename Problem::template OutputSamples<> plus;   
    typename Problem::template OutputSamples<> minus;  
};

template <typename Problem, typename Spread>
[[nodiscard]] SigmaOutputs<Problem>
evaluateSigmaPoints(
    Problem const & problem,
    typename Problem::template InputVector<typename Problem::Scalar> const & mu,
    Spread const & spread
)
{
    using Scalar = typename Problem::Scalar;
    using InputVector = typename Problem::template InputVector<Scalar>;
    constexpr int N = Problem::InputDifferenceSize;
    constexpr int NMax = Problem::InputDifferenceMaxSize;
    using Outputs = typename Problem::template OutputSamples<
        N == Eigen::Dynamic ? Eigen::Dynamic : 2 * N + 1,
        NMax == Eigen::Dynamic ? Eigen::Dynamic : 2 * NMax + 1>;
    Eigen::Index const n = spread.cols();
    auto outputs = Outputs{};
    problem.evaluateEach(
        2 * n + 1,
        [&](Eigen::Index k) -> InputVector {
            if (k == 0)
            {
                return mu;
            }
            Eigen::Index const j = (k - 1) / 2;
            return k % 2 == 1 ? problem.plusInputs(mu, spread.col(j)) : problem.plusInputs(mu, -spread.col(j));
        },
        outputs
    );
    auto result = SigmaOutputs<Problem>{};
    result.center = outputs.col(0);
    // Sized by resize(), since braces would set the coefficients of a fixed-size vector.
    result.plus.resize(outputs.rows(), n);
    result.minus.resize(outputs.rows(), n);
    // Over no coordinates there are no plus and minus points, at compile time.
    if constexpr (N != 0)
    {
        for (Eigen::Index j = 0; j < n; ++j)
        {
            result.plus.col(j) = outputs.col(1 + 2 * j);
            result.minus.col(j) = outputs.col(2 + 2 * j);
        }
    }
    return result;
}

template <UnscentedRule Rule, typename Problem, typename A = roots::DefaultAxes>
[[nodiscard]] SigmaRoot<Problem>
transformRoot(Rule const & method, Problem const & problem, A const & axes = {})
{
    using Scalar = typename Problem::Scalar;
    using Result = SigmaRoot<Problem>;
    Eigen::Index const n = problem.inputDifferenceSize();
    SigmaWeights<Scalar> const weights = sigmaWeights<Scalar, Problem::InputDifferenceSize>(method, n);
 
    typename Problem::template InputVector<Scalar> const mu{problem.mean()};
    typename Problem::InputMatrix const spread = weights.scale * problem.root(axes);
    auto const outputs = evaluateSigmaPoints(problem, mu, spread);
 
    auto moments = sigmaMoments(problem, weights, outputs.center, outputs.plus, outputs.minus);
    Eigen::Index const m = problem.outputDifferenceSize();
    Scalar const pointRoot = Eigen::numext::sqrt(weights.point);
    Scalar const centerRoot = Eigen::numext::sqrt(Eigen::numext::abs(weights.centerCov));
    typename Result::RootMatrix root(m, 2 * n + 1);
    root.leftCols(n) = pointRoot * moments.plus;
    root.middleCols(n, n) = pointRoot * moments.minus;
    root.col(2 * n) = centerRoot * moments.center;
    return Result{std::move(moments.mean), std::move(root), weights.centerCov < Scalar(0)};
}

template <UnscentedRule Rule, typename Problem, typename A = roots::DefaultAxes>
[[nodiscard]] typename Problem::RootResult
transformResult(Rule const & method, Problem const & problem, A const & axes = {})
{
    using Scalar = typename Problem::Scalar;
    Eigen::Index const n = problem.inputDifferenceSize();
    SigmaWeights<Scalar> const weights = sigmaWeights<Scalar, Problem::InputDifferenceSize>(method, n);
 
    typename Problem::template InputVector<Scalar> const mu{problem.mean()};
    typename Problem::InputMatrix root = problem.root(axes);
    typename Problem::InputMatrix const spread = weights.scale * root;
    auto const outputs = evaluateSigmaPoints(problem, mu, spread);
 
    auto moments = sigmaMoments(problem, weights, outputs.center, outputs.plus, outputs.minus);
    return {
        std::move(moments.mean),
        std::move(root),
        typename Problem::CrossMatrix{std::move(moments.whitened)},
        std::move(moments.unexplained)
    };
}

template <PointSet Points>
[[nodiscard]] consteval int
pointCount(int n)
{
    if (n == Eigen::Dynamic)
    {
        return Eigen::Dynamic;
    }
    if (n == 0)
    {
        return 1;
    }
    return int(std::min<Eigen::Index>(Points::pointCount(n), MaxPointCount + 1));
}

template <PointSet Points>
[[nodiscard]] Eigen::Index
pointCountAt(Eigen::Index n)
{
    return n == 0 ? Eigen::Index(1) : Points::pointCount(n);
}

template <typename Scalar, typename Points, int N, int NMax>
struct UnitPoints
{
    static constexpr int Count = pointCount<Points>(N);        
    static constexpr int MaxCount = pointCount<Points>(NMax);  
 
    static_assert(
        NMax == Eigen::Dynamic || MaxCount <= std::max(2 * NMax + 2, MaxPointCount),
        "a point set over this many coordinates has more than MaxPointCount (4096) points; declare the "
        "coordinates the model is affine in with affineIn(), or use a set with fewer points, such as a Gauss-Hermite "
        "rule of a lower order"
    );
 
    MatrixOf<Scalar, N, Count, NMax, MaxCount> points;  
    MatrixOf<Scalar, Count, 1, MaxCount, 1> weights;    
    Eigen::Index reference;  
};

template <typename Scalar, int N, int NMax, PointSet Points>
[[nodiscard]] UnitPoints<Scalar, Points, N, NMax>
unitPoints(Points const & /*set*/, Eigen::Index n)
{
    using Result = UnitPoints<Scalar, Points, N, NMax>;
    Eigen::Index const count = pointCountAt<Points>(n);
    auto result = Result{};
    result.points.resize(n, count);
    result.weights.resize(count);
    result.reference = 0;
    if (n == 0)
    {
        result.weights(0) = Scalar(1);
        return result;
    }
    result.reference = Points::unitPoints(n, result.points, result.weights);
    EMRYS_FILTER_ASSERT(
        result.reference >= 0 && result.reference < count
        && "the point the mean is taken about must be one of the points"
    );
    return result;
}

template <typename Problem, typename Deviations, typename Whitened>
struct PointMoments
{
    typename Problem::template OutputVector<typename Problem::Scalar> mean;  
    Whitened whitened;                           
    typename Problem::OutputMatrix unexplained;  
    Deviations deviations;                       
};

template <typename Problem, typename Samples, typename Points, typename Weights, typename Root>
[[nodiscard]] auto
pointMoments(
    Problem const & problem,
    Points const & points,
    Weights const & weights,
    Eigen::Index reference,
    Samples const & outputs,
    Root const & root
)
{
    using Scalar = typename Problem::Scalar;
    using OutputVector = typename Problem::template OutputVector<Scalar>;
    using OutputDeviations =
        typename Problem::template OutputDeviations<Samples::ColsAtCompileTime, Samples::MaxColsAtCompileTime>;
    auto const & space = problem.outputSpace();
    OutputVector const center = outputs.col(reference);
    OutputDeviations deviations = Emrys::minus(space, outputs, center);
    typename Problem::template OutputDifferenceVector<Scalar> step;
    step.noalias() = multiply(deviations, weights);
    OutputVector mean = Emrys::plus(space, center, step);
    deviations = Emrys::minus(space, outputs, mean);
    using Whitened = WhitenedOf<Problem, Points::RowsAtCompileTime, Points::MaxRowsAtCompileTime>;
    Whitened whitened = multiply(points * weights.asDiagonal(), deviations.transpose());
    EMRYS_FILTER_ASSERT(root.cols() == whitened.rows() && "the root must have a column per coordinate of the points");
    for (Eigen::Index j = 0; j < root.cols(); ++j)
    {
        if ((root.col(j).array() == Scalar(0)).all())
        {
            whitened.row(j).setZero();
        }
    }
    OutputDeviations residuals = deviations;
    residuals.noalias() -= multiply(whitened.transpose(), points);
    typename Problem::OutputMatrix unexplained = multiply(residuals * weights.asDiagonal(), residuals.transpose());
    unexplained = (Scalar(0.5) * (unexplained + unexplained.transpose())).eval();
    return PointMoments<Problem, OutputDeviations, Whitened>{
        std::move(mean),
        std::move(whitened),
        std::move(unexplained),
        std::move(deviations)
    };
}

template <typename Problem, typename Unit, typename Spread>
[[nodiscard]] auto
evaluatePoints(Problem const & problem, Unit const & unit, Spread const & spread)
{
    using Scalar = typename Problem::Scalar;
    using OutputSamples = typename Problem::template OutputSamples<Unit::Count, Unit::MaxCount>;
    typename Problem::template InputVector<Scalar> const mu{problem.mean()};
    auto outputs = OutputSamples{};
    problem.evaluateEach(
        unit.weights.size(),
        [&](Eigen::Index j) {
            return problem.plusInputs(mu, spread.col(j));
        },
        outputs
    );
    return outputs;
}

template <typename Points, typename Problem, typename A = roots::DefaultAxes>
    requires isPointSet<Points>
[[nodiscard]] auto
transformRoot(Points const & set, Problem const & problem, A const & axes = {})
{
    using Scalar = typename Problem::Scalar;
    using Unit = UnitPoints<Scalar, Points, Problem::InputDifferenceSize, Problem::InputDifferenceMaxSize>;
    using Result = SigmaRoot<Problem, Unit::Count, Unit::MaxCount>;
    using Spread =
        MatrixOf<Scalar, Problem::InputDifferenceSize, Unit::Count, Problem::InputDifferenceMaxSize, Unit::MaxCount>;
    Unit const unit = unitPoints<Scalar, Problem::InputDifferenceSize, Problem::InputDifferenceMaxSize>(
        set,
        problem.inputDifferenceSize()
    );
    typename Problem::InputMatrix const inputRoot = problem.root(axes);
    auto spread = Spread{};
    spread.resize(unit.points.rows(), unit.points.cols());
    spread.noalias() = multiply(inputRoot, unit.points);
    auto moments = pointMoments(
        problem,
        unit.points,
        unit.weights,
        unit.reference,
        evaluatePoints(problem, unit, spread),
        inputRoot
    );
    typename Result::RootMatrix root = moments.deviations * unit.weights.cwiseSqrt().asDiagonal();
    return Result{std::move(moments.mean), std::move(root), false};
}

template <typename Points, typename Problem, typename A = roots::DefaultAxes>
    requires isPointSet<Points>
[[nodiscard]] typename Problem::RootResult
transformResult(Points const & set, Problem const & problem, A const & axes = {})
{
    using Scalar = typename Problem::Scalar;
    using Unit = UnitPoints<Scalar, Points, Problem::InputDifferenceSize, Problem::InputDifferenceMaxSize>;
    using Spread =
        MatrixOf<Scalar, Problem::InputDifferenceSize, Unit::Count, Problem::InputDifferenceMaxSize, Unit::MaxCount>;
    Unit const unit = unitPoints<Scalar, Problem::InputDifferenceSize, Problem::InputDifferenceMaxSize>(
        set,
        problem.inputDifferenceSize()
    );
    typename Problem::InputMatrix root = problem.root(axes);
    auto spread = Spread{};
    spread.resize(unit.points.rows(), unit.points.cols());
    spread.noalias() = multiply(root, unit.points);
    auto moments =
        pointMoments(problem, unit.points, unit.weights, unit.reference, evaluatePoints(problem, unit, spread), root);
    return {
        std::move(moments.mean),
        std::move(root),
        typename Problem::CrossMatrix{std::move(moments.whitened)},
        std::move(moments.unexplained)
    };
}

template <typename Rule, typename A, typename Problem>
[[nodiscard]] typename Problem::RootResult
transformResult(methods::Spread<Rule, A> const & method, Problem const & problem)
{
    return transformResult(method.rule, problem, method.axes);
}

template <typename Rule, typename A, typename Problem>
[[nodiscard]] auto
transformRoot(methods::Spread<Rule, A> const & method, Problem const & problem)
{
    return transformRoot(method.rule, problem, method.axes);
}

template <typename Problem, typename X, typename Affine, typename Columns>
[[nodiscard]] typename Problem::template OutputVector<typename Problem::Scalar>
valueWithAffineJacobian(Problem const & problem, X const & x, Affine const & affine, Columns && columns)
{
    using Scalar = typename Problem::Scalar;
    using OutputVector = typename Problem::template OutputVector<Scalar>;
    using Index = Eigen::Index;
    Index const l = affine.size();
    if constexpr (returnsJacobian<Problem>)
    {
        auto returned = problem.template returned<Scalar>(x);
        using Blocks = decltype(returned.jacobian);
        problem.forEachUncertain([&]<std::size_t I>(auto const & input) {
            auto const & block = [&]() -> auto const & {
                if constexpr (isArgumentJacobians<Blocks>)
                {
                    return std::get<I>(returned.jacobian.blocks);
                }
                else
                {
                    return jacobianBlock<Problem::ArgInfo::ordinal[I]>(returned.jacobian);
                }
            }();
            Index const start = problem.differenceOffset(I);
            Index const size = input.differenceDimension();
            for (Index j = 0; j < l; ++j)
            {
                if (affine(j) >= start && affine(j) < start + size)
                {
                    columns.col(j) = block.col(affine(j) - start);
                }
            }
        });
        OutputVector y = std::move(returned.value);
        EMRYS_FILTER_ASSERT(
            y.size() == valueDimension(problem.outputSpace()) && "the output does not fit the output space"
        );
        return y;
    }
    else
    {
        using Derivatives =
            Eigen::Matrix<Scalar, Affine::SizeAtCompileTime, 1, Eigen::ColMajor, Affine::MaxSizeAtCompileTime, 1>;
        using Dual = Eigen::AutoDiffScalar<Derivatives>;
        static_assert(
            returnsDuals<Problem, Dual>,
            "a function that a sigma-point method factors affine coordinates out of returns its Jacobian, or is "
            "generic "
            "in its scalar, to be differentiated in those coordinates"
        );
        using OutputSpace = typename Problem::OutputSpaceType;
        Index const total = problem.inputDifferenceSize();
        typename Problem::template InputDifferenceVector<Dual> delta;
        delta.resize(total);
        for (Index i = 0; i < total; ++i)
        {
            delta(i) = Dual{Scalar(0), Derivatives::Zero(l)};
        }
        for (Index j = 0; j < l; ++j)
        {
            delta(affine(j)) = Dual{Scalar(0), int(l), int(j)};
        }
        typename Problem::template InputVector<Dual> const dual{problem.stepInputs(x.template cast<Dual>(), delta)};
        typename Problem::template OutputVector<Dual> const y{problem.template evaluate<Dual>(dual)};
        OutputVector value = OutputVector::Zero(y.size());
        for (Index r = 0; r < y.size(); ++r)
        {
            value(r) = y(r).value();
        }
        auto const deviation = [&] {
            if constexpr (HasMinus<OutputSpace>)
            {
                return typename Problem::template OutputDifferenceVector<Dual>{
                    Emrys::minus(problem.outputSpace(), y, value.template cast<Dual>())
                };
            }
            else
            {
                return y;
            }
        }();
        for (Index r = 0; r < deviation.size(); ++r)
        {
            // An output that does not depend on the affine coordinates may carry no derivatives.
            if (deviation(r).derivatives().size() > 0)
            {
                columns.row(r) = deviation(r).derivatives().transpose();
            }
            else
            {
                columns.row(r).setZero();
            }
        }
        return value;
    }
}

template <typename Problem, typename Affine, typename ShiftAt, typename Outputs, typename Jacobians>
void
evaluateWithAffineJacobians(
    Problem const & problem,
    typename Problem::template InputVector<typename Problem::Scalar> const & mu,
    Affine const & affine,
    Eigen::Index count,
    ShiftAt const & shiftAt,
    Outputs & outputs,
    Jacobians & jacobians
)
{
    using Index = Eigen::Index;
    Index const l = affine.size();
    if constexpr (Problem::EvaluatesBatches)
    {
        EMRYS_FILTER_ASSERT(jacobians.cols() == l * count && "a Jacobian over the affine coordinates per point");
        problem.evaluateEachWithAffineJacobians(
            count,
            [&](Index k) {
                return problem.plusInputs(mu, shiftAt(k));
            },
            outputs,
            jacobians
        );
    }
    else
    {
        for (Index k = 0; k < count; ++k)
        {
            auto const y = valueWithAffineJacobian(
                problem,
                problem.plusInputs(mu, shiftAt(k)),
                affine,
                jacobians.middleCols(k * l, l)
            );
            if (k == 0)
            {
                outputs.resize(y.size(), count);
            }
            outputs.col(k) = y;
        }
    }
}

template <typename Problem, int AffineSize, int AffineMaxSize, typename Points>
[[nodiscard]] typename Problem::RootResult
transformResult(FactoredAt<AffineSize, AffineMaxSize, Points> const & method, Problem const & problem)
{
    using Scalar = typename Problem::Scalar;
    using Index = Eigen::Index;
    using InputVector = typename Problem::template InputVector<Scalar>;
    using InputDifferenceVector = typename Problem::template InputDifferenceVector<Scalar>;
    using OutputVector = typename Problem::template OutputVector<Scalar>;
    using OutputMatrix = typename Problem::OutputMatrix;
    constexpr int OutputSize = Problem::OutputDifferenceSize;
    constexpr int OutputMaxSize = Problem::OutputDifferenceMaxSize;
    // The rest of the joint input's differences: the coordinates that get sigma points.
    constexpr int RestSize = Problem::InputDifferenceSize == Eigen::Dynamic || AffineSize == Eigen::Dynamic
                               ? Eigen::Dynamic
                               : Problem::InputDifferenceSize - AffineSize;
    constexpr int RestMaxSize = RestSize == Eigen::Dynamic ? Problem::InputDifferenceMaxSize : RestSize;
    using RestIndices = IndexVector<RestSize, RestMaxSize>;
    using OutputSamples = typename Problem::template OutputSamples<RestSize, RestMaxSize>;
    using RestMatrix = MatrixOf<Scalar, RestSize, RestSize, RestMaxSize, RestMaxSize>;
    using RestAffine = MatrixOf<Scalar, RestSize, AffineSize, RestMaxSize, AffineMaxSize>;
    using RestCross = MatrixOf<Scalar, RestSize, OutputSize, RestMaxSize, OutputMaxSize>;
    using AffineMatrix = MatrixOf<Scalar, AffineSize, AffineSize, AffineMaxSize, AffineMaxSize>;
    using AffineRest = MatrixOf<Scalar, AffineSize, RestSize, AffineMaxSize, RestMaxSize>;
    using AffineCross = MatrixOf<Scalar, AffineSize, OutputSize, AffineMaxSize, OutputMaxSize>;
    using AffineJacobian = MatrixOf<Scalar, OutputSize, AffineSize, OutputMaxSize, AffineMaxSize>;
 
    auto const & affine = method.affine;
    Index const total = problem.inputDifferenceSize();
    Index const l = affine.size();
    Index const n = total - l;
    EMRYS_FILTER_ASSERT(
        distinctIndices(affine) && (l == 0 || (affine.minCoeff() >= 0 && affine.maxCoeff() < total))
        && "affine coordinates must be distinct coordinates of the joint input"
    );
    auto rest = RestIndices{};
    rest.resize(n);
    for (Index i = 0, k = 0; i < total; ++i)
    {
        if ((affine.array() != i).all())
        {
            rest(k++) = i;
        }
    }
 
    // x_a given x_n: mean mu_a + gain (x_n - mu_n), covariance conditional; points from a root of P_nn.
    InputVector const mu = problem.mean();
    typename Problem::InputMatrix const cov = problem.covarianceLower();
    AffineMatrix conditional = lowerSymmetricBlock(cov, affine, affine);
    RestMatrix restRoot = RestMatrix::Zero(n, n);
    AffineRest gain = AffineRest::Zero(l, n);
    if constexpr (RestSize != 0)
    {
        if (n > 0)
        {
            auto const ldlt = Eigen::LDLT<RestMatrix>{lowerSymmetricBlock(cov, rest, rest)};
            RestAffine const restAffine = lowerSymmetricBlock(cov, rest, affine);
            gain = transposedSolve<AffineRest>(ldlt, restAffine);
            conditional.noalias() -= multiply(gain, restAffine);
            refuseFailingAxes<AxesOf<Points>>();
            restRoot = rootValue(axesOf(method.points).root(lowerSymmetricBlock(cov, rest, rest)));
        }
    }
    conditional = (Scalar(0.5) * (conditional + conditional.transpose())).eval();
    // L_c, the root of P_a|n, which spreads x_a given x_n as the rest's root spreads x_n.
    AffineMatrix conditionalRoot = AffineMatrix::Zero(l, l);
    if constexpr (AffineSize != 0)
    {
        if (l > 0)
        {
            conditionalRoot = defaultRoot(conditional);
        }
    }
 
    // f at each point, and its Jacobian in x_a there kept as column block k of jacobians, in one
    // evaluation of the points (evaluateWithAffineJacobians()).
    Index const m = problem.outputDifferenceSize();
    // The affine coordinates' parts (decisions.md [linearization-forms]): with the mean Jacobian H_a = sum_i w_i h_i,
    // Psi_a = L_c^T H_a^T, and Omega_a = sum_i w_i (h_i - H_a) P_a|n (h_i - H_a)^T in a second pass
    // over the Jacobians kept, so that neither is recovered by subtraction.
    auto const affineParts = [&](auto const & jacobians, auto const & pointWeights) {
        Index const count = pointWeights.size();
        AffineJacobian meanJacobian = AffineJacobian::Zero(m, l);
        for (Index k = 0; k < count; ++k)
        {
            meanJacobian += pointWeights(k) * jacobians.middleCols(k * l, l);
        }
        OutputMatrix omega = OutputMatrix::Zero(m, m);
        for (Index k = 0; k < count; ++k)
        {
            AffineJacobian const deviation = jacobians.middleCols(k * l, l) - meanJacobian;
            AffineCross const spread = multiply(conditional, deviation.transpose());
            omega.noalias() += pointWeights(k) * multiply(deviation, spread);
        }
        AffineCross whitened = multiply(conditionalRoot.transpose(), meanJacobian.transpose());
        return std::pair{std::move(whitened), std::move(omega)};
    };
 
    // The output mean, the rows of x_n in Psi and the points' Omega, and x_a's parts.
    // With every coordinate affine there is nothing to spread: the center alone, by any rule.
    auto [mean, restWhitened, pointOmega, affineWhitened, affineOmega] = [&] {
        if constexpr (isUnscented<RuleOf<Points>> || RestSize == 0)
        {
            constexpr int Count = RestSize == Eigen::Dynamic ? Eigen::Dynamic : 2 * RestSize + 1;
            constexpr int MaxCount = RestMaxSize == Eigen::Dynamic ? Eigen::Dynamic : 2 * RestMaxSize + 1;
            using Jacobians = MatrixOf<
                Scalar,
                OutputSize,
                sizeProduct<AffineSize, Count>,
                OutputMaxSize,
                sizeProduct<AffineMaxSize, MaxCount>>;
            using PointWeights = MatrixOf<Scalar, Count, 1, MaxCount, 1>;
            SigmaWeights<Scalar> const weights = [&] {
                if constexpr (RestSize == 0)
                {
                    return SigmaWeights<Scalar>{Scalar(0), Scalar(1), Scalar(1), Scalar(0)};
                }
                else
                {
                    return sigmaWeights<Scalar, RestSize>(ruleOf(method.points), n);
                }
            }();
            RestMatrix const spread = weights.scale * restRoot;
            AffineRest const affineSpread = [&] {
                if constexpr (RestSize != 0)
                {
                    return AffineRest{multiply(gain, spread)};
                }
                else
                {
                    return AffineRest{};
                }
            }();
            // The center, then the plus and minus points of each column of the spread in turn.
            auto const shiftAt = [&]([[maybe_unused]] Index k) {
                InputDifferenceVector shift = InputDifferenceVector::Zero(total);
                if constexpr (RestSize != 0)
                {
                    if (k > 0)
                    {
                        Index const j = (k - 1) / 2;
                        shift(rest) = spread.col(j);
                        shift(affine) = affineSpread.col(j);
                        if (k % 2 == 0)
                        {
                            shift = -shift;
                        }
                    }
                }
                return shift;
            };
            auto jacobians = Jacobians{};
            jacobians.resize(m, l * (2 * n + 1));
            auto outputs = typename Problem::template OutputSamples<Count, MaxCount>{};
            evaluateWithAffineJacobians(problem, mu, affine, 2 * n + 1, shiftAt, outputs, jacobians);
            auto pointWeights = PointWeights{};
            pointWeights.resize(2 * n + 1);
            pointWeights(0) = weights.center;
            OutputVector const center = outputs.col(0);
            // Sized by resize(), since braces would set the coefficients of a fixed-size vector.
            auto plus = OutputSamples{};
            auto minus = OutputSamples{};
            plus.resize(center.size(), n);
            minus.resize(center.size(), n);
            // With every coordinate affine there are no plus and minus points, at compile time.
            if constexpr (RestSize != 0)
            {
                for (Index j = 0; j < n; ++j)
                {
                    plus.col(j) = outputs.col(1 + 2 * j);
                    minus.col(j) = outputs.col(2 + 2 * j);
                    pointWeights(1 + 2 * j) = weights.point;
                    pointWeights(2 + 2 * j) = weights.point;
                }
            }
            auto sigma = sigmaMoments(problem, weights, center, plus, minus);
            auto [whitened, omega] = affineParts(jacobians, pointWeights);
            return std::tuple{
                std::move(sigma.mean),
                RestCross{std::move(sigma.whitened)},
                std::move(sigma.unexplained),
                std::move(whitened),
                std::move(omega)
            };
        }
        else
        {
            using Unit = UnitPoints<Scalar, RuleOf<Points>, RestSize, RestMaxSize>;
            using RestSpread = MatrixOf<Scalar, RestSize, Unit::Count, RestMaxSize, Unit::MaxCount>;
            using AffineSpread = MatrixOf<Scalar, AffineSize, Unit::Count, AffineMaxSize, Unit::MaxCount>;
            using Jacobians = MatrixOf<
                Scalar,
                OutputSize,
                sizeProduct<AffineSize, Unit::Count>,
                OutputMaxSize,
                sizeProduct<AffineMaxSize, Unit::MaxCount>>;
            Unit const unit = unitPoints<Scalar, RestSize, RestMaxSize>(ruleOf(method.points), n);
            Index const count = unit.weights.size();
            auto spread = RestSpread{};
            spread.resize(n, count);
            spread.noalias() = multiply(restRoot, unit.points);
            auto affineSpread = AffineSpread{};
            affineSpread.resize(l, count);
            affineSpread.noalias() = multiply(gain, spread);
            auto jacobians = Jacobians{};
            jacobians.resize(m, l * count);
            auto outputs = typename Problem::template OutputSamples<Unit::Count, Unit::MaxCount>{};
            auto const shiftAt = [&](Index j) {
                InputDifferenceVector shift = InputDifferenceVector::Zero(total);
                shift(rest) = spread.col(j);
                shift(affine) = affineSpread.col(j);
                return shift;
            };
            evaluateWithAffineJacobians(problem, mu, affine, count, shiftAt, outputs, jacobians);
            auto points = pointMoments(problem, unit.points, unit.weights, unit.reference, outputs, restRoot);
            auto [whitened, omega] = affineParts(jacobians, unit.weights);
            return std::tuple{
                std::move(points.mean),
                RestCross{std::move(points.whitened)},
                std::move(points.unexplained),
                std::move(whitened),
                std::move(omega)
            };
        }
    }();
 
    // The joint root, with x_n = mu_n + R_n xi_n and x_a = mu_a + G R_n xi_n + L_c xi_a: block-diagonal
    // by input, since G and L_c are, and Psi's rows at x_n from the points and at x_a from L_c.
    typename Problem::InputMatrix root = Problem::InputMatrix::Zero(total, total);
    typename Problem::CrossMatrix whitened = Problem::CrossMatrix::Zero(total, m);
    if constexpr (RestSize != 0)
    {
        AffineRest const affineRoot = multiply(gain, restRoot);
        for (Index k = 0; k < n; ++k)
        {
            whitened.row(rest(k)) = restWhitened.row(k);
            for (Index j = 0; j < n; ++j)
            {
                root(rest(k), rest(j)) = restRoot(k, j);
            }
            for (Index j = 0; j < l; ++j)
            {
                root(affine(j), rest(k)) = affineRoot(j, k);
            }
        }
    }
    for (Index k = 0; k < l; ++k)
    {
        whitened.row(affine(k)) = affineWhitened.row(k);
        for (Index j = 0; j < l; ++j)
        {
            root(affine(k), affine(j)) = conditionalRoot(k, j);
        }
    }
    OutputMatrix unexplained = pointOmega + affineOmega;
    unexplained = (Scalar(0.5) * (unexplained + unexplained.transpose())).eval();
    return {std::move(mean), std::move(root), std::move(whitened), std::move(unexplained)};
}

template <Space Space_, RealScalar Scalar_>
class PointInput
{
public:
 
    using Scalar = Scalar_;                               
    using SpaceType = Space_;                             
    using MomentsType = Moments<Space_, Scalar_>;         
    using VectorType = typename MomentsType::VectorType;  
    using MatrixType = typename MomentsType::MatrixType;  

    static constexpr int ValueDimensionAtCompileTime = MomentsType::ValueDimensionAtCompileTime;

    static constexpr int MaxValueDimensionAtCompileTime = MomentsType::MaxValueDimensionAtCompileTime;

    static constexpr int DifferenceDimensionAtCompileTime = MomentsType::DifferenceDimensionAtCompileTime;

    static constexpr int MaxDifferenceDimensionAtCompileTime = MomentsType::MaxDifferenceDimensionAtCompileTime;

    PointInput(Space_ space, VectorType point)
        : m_space{std::move(space)}
        , m_point{std::move(point)}
    {}

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

    [[nodiscard]] VectorType const &
    mean() const
    {
        return m_point;
    }

    [[nodiscard]] MatrixType
    covarianceLower() const
    {
        return MatrixType::Zero(differenceDimension(), differenceDimension());
    }

    [[nodiscard]] Eigen::Index
    valueDimension() const
    {
        return Emrys::valueDimension(m_space);
    }

    [[nodiscard]] Eigen::Index
    differenceDimension() const
    {
        return Emrys::differenceDimension(m_space);
    }
 
private:
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Space_ m_space;  
    VectorType m_point;                             
};

template <typename InputSpace, typename OutputSpace, typename Scalar>
struct LinearizationAbout
{
    using Linearization = StatisticalLinearization<InputSpace, OutputSpace, Scalar>;
 
    using JacobianMatrix = typename Linearization::JacobianMatrix;                  
    using ErrorMatrix = typename Linearization::ErrorMatrix;                        
    using OutputVector = typename GaussianParams<OutputSpace, Scalar>::VectorType;  
 
    OutputSpace space;        
    OutputVector value;       
    JacobianMatrix jacobian;  
    ErrorMatrix inputs;       

    ErrorMatrix unexplained;

    [[nodiscard]] ErrorMatrix
    error() const
    {
        return inputs + unexplained;
    }
};

template <typename InputSpace, typename Point, typename Base, typename OutputSpace, typename Scalar>
void
toChartOf(
    InputSpace const & space,
    Eigen::MatrixBase<Base> const & base,
    Point const & point,
    LinearizationAbout<InputSpace, OutputSpace, Scalar> & linearization
)
{
    using Linearization = LinearizationAbout<InputSpace, OutputSpace, Scalar>;
    using Difference = typename Moments<InputSpace, Scalar>::DifferenceVectorType;
    using OutputDifference = typename GaussianParams<OutputSpace, Scalar>::DifferenceVectorType;
    Difference const step = Emrys::minus(space, point, base);
    if constexpr (!hasIdentityReset<InputSpace>)
    {
        auto carried = typename Linearization::JacobianMatrix{};
        carried.resize(linearization.jacobian.rows(), linearization.jacobian.cols());
        carried.noalias() = multiply(linearization.jacobian, resetJacobian(space, base, step));
        linearization.jacobian = std::move(carried);
    }
    auto back = OutputDifference{};
    back.resize(linearization.jacobian.rows());
    back.noalias() = -multiply(linearization.jacobian, step);
    linearization.value = Emrys::plus(linearization.space, linearization.value, back);
}

template <typename Method, typename About, typename InputSpace>
[[nodiscard]] consteval bool
checkReference()
{
    using Reference = typename About::ReferenceType;
    static_assert(
        !isIterated<Method>,
        "an Iterated update chooses its own points to linearize at, so give linearizedAbout() to a step by its "
        "inner method instead"
    );
    if constexpr (About::IsDistribution)
    {
        static_assert(
            std::is_same_v<typename std::remove_cvref_t<Reference>::SpaceType, InputSpace>,
            "linearizedAbout() takes a point or a distribution in the space of the input it is for"
        );
    }
    static_assert(
        isLinearMethod<Method> || About::IsDistribution,
        "a sigma-point method linearizes over a distribution, so give linearizedAbout() one, such as another "
        "filter's state, rather than a point"
    );
    return true;
}

template <
    typename Method,
    typename F,
    typename OutputSpace,
    typename InputSpace,
    typename Base,
    typename Reference,
    typename... Args>
[[nodiscard]] auto
linearizeAbout(
    Method const & method,
    F & f,
    OutputSpace const & outputSpace,
    InputSpace const & space,
    Eigen::MatrixBase<Base> const & base,
    LinearizedAbout<Reference> const & about,
    Args const &... args
)
{
    using Scalar = typename Base::Scalar;
    using About = LinearizedAbout<Reference>;
    using Point = PointInput<InputSpace, Scalar>;
    using Result = LinearizationAbout<InputSpace, OutputSpace, Scalar>;
    static_assert(
        crossCount<Args...> == 0,
        "withCross() may mark only the input linearizedAbout() is for, which is the argument right after it"
    );
    static_assert(checkReference<Method, About, InputSpace>());
    typename Point::VectorType const point = about.mean();
    EMRYS_FILTER_ASSERT(
        point.size() == valueDimension(space) && "a point to linearize about must have the input's value dimension"
    );
    Eigen::Index const n = differenceDimension(space);
    auto result = [&] {
        if constexpr (isLinearMethod<Method>)
        {
            auto const input = Point{space, point};
            using Problem = Transform<F, OutputSpace, Point, Args...>;
            auto const problem = Problem{f, outputSpace, input, args...};
            auto const linearized = linearizeAt(method, problem);
            // The other inputs' share, J Sigma J^T over the joint input, whose block at the point is zero.
            typename Problem::CrossMatrix const cross =
                multiply(problem.covarianceLower(), linearized.jacobian.transpose());
            typename Result::ErrorMatrix inputs = multiply(linearized.jacobian, cross);
            inputs = (Scalar(0.5) * (inputs + inputs.transpose())).eval();
            Eigen::Index const m = inputs.rows();
            return Result{
                problem.outputSpace(),
                typename Result::OutputVector{linearized.mean},
                typename Result::JacobianMatrix{linearized.jacobian.middleCols(problem.differenceOffset(0), n)},
                std::move(inputs),
                Result::ErrorMatrix::Zero(m, m)
            };
        }
        else if constexpr (About::IsDistribution)
        {
            using Distribution = std::remove_cvref_t<Reference>;
            using Marked = CrossMarked<Distribution>;
            using Problem = Transform<F, OutputSpace, Marked, Args...>;
            auto const marked = Marked{about.reference};
            auto const problem = Problem{f, outputSpace, marked, args...};
            auto const predicted = problem.assemble(transformResult(method, problem));
            return Result{
                problem.outputSpace(),
                typename Result::OutputVector{predicted.output().mean()},
                predicted.jacobian(),
                predicted.inputsCovariance(),
                predicted.unexplainedCovariance()
            };
        }
    }();
    toChartOf(space, base, point, result);
    return result;
}

template <typename Method, typename F, typename OutputSpace, typename Reference, typename X, typename... Args>
[[nodiscard]] auto
transformAbout(
    Method const & method,
    F & f,
    OutputSpace const & outputSpace,
    LinearizedAbout<Reference> const & about,
    X const & x,
    Args const &... args
)
{
    static_assert(ArgTraits<X>::isUncertain, "linearizedAbout() is followed by the uncertain input it is for");
    using Input = typename ArgTraits<X>::Input;
    using InputSpace = typename Input::SpaceType;
    using Scalar = typename Input::Scalar;
    using Output = GaussianParams<OutputSpace, Scalar>;
    using Covariance = typename Moments<InputSpace, Scalar>::MatrixType;
    using Cross = CrossMatrixOf<InputSpace, OutputSpace, Scalar>;
    auto const & input = argInput(x);
    auto linearization = linearizeAbout(method, f, outputSpace, input.space(), input.mean(), about, args...);
    Covariance const cov = symmetricFromLower<Covariance>(covarianceLowerOf(input));
    Cross cross = multiply(cov, linearization.jacobian.transpose());
    typename Output::MatrixType explained = multiply(linearization.jacobian, cross);
    explained = (Scalar(0.5) * (explained + explained.transpose())).eval();
    if constexpr (ArgTraits<X>::hasCross)
    {
        return OutputWithCross<InputSpace, OutputSpace, JacobianForm, Scalar>{
            std::move(linearization.space),
            std::move(linearization.value),
            std::move(linearization.jacobian),
            std::move(cross),
            std::move(explained),
            std::move(linearization.inputs),
            std::move(linearization.unexplained)
        };
    }
    else
    {
        typename Output::MatrixType outputCov = explained + linearization.inputs + linearization.unexplained;
        return Output{linearization.space, linearization.value, std::move(outputCov)};
    }
}
 
}  // namespace internal

template <UncertainInput Input>
[[nodiscard]] CrossMarked<Input>
withCross(Input const & input)
{
    return {input};
}

template <UncertainInput Input>
void
withCross(Input const && input) = delete;

template <UncertainInput Distribution>
[[nodiscard]] LinearizedAbout<Distribution const &>
linearizedAbout(Distribution const & distribution)
{
    return {distribution};
}

template <UncertainInput Distribution>
void
linearizedAbout(Distribution const && distribution) = delete;

template <typename Derived>
[[nodiscard]] LinearizedAbout<typename Derived::PlainObject>
linearizedAbout(Eigen::MatrixBase<Derived> const & point)
{
    static_assert(Derived::ColsAtCompileTime == 1, "a point to linearize at is a column vector");
    return {typename Derived::PlainObject{point}};
}
 
namespace internal
{

template <typename Method, typename F, typename OutputSpace, typename... Args>
    requires(IsSpace<OutputSpace> && crossCount<Args...> <= 1)
[[nodiscard]] auto
transformBy(Method const & method, F && f, OutputSpace const & outputSpace, Args const &... args)
{
    if constexpr (leadsWithReference<Args...>)
    {
        return transformAbout(method, f, outputSpace, args...);
    }
    else
    {
        using Problem = Transform<std::remove_reference_t<F>, OutputSpace, Args...>;
        Problem const problem{f, outputSpace, args...};
        return problem.assemble(transformResult(method, problem));
    }
}

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

template <typename T>
    requires(ArgTraits<T>::isUncertain)
struct TransformInputTagOf<T>
{
    using type = typename ArgTraits<T>::Input::SpaceType;
};
 
template <typename S, typename Vector>
struct TransformInputTagOf<Sure<S, Vector>>
{
    using type = S;
};


template <typename T>
using TransformInputTag = typename TransformInputTagOf<T>::type;

template <typename... Args>
inline constexpr std::size_t leadingPolicies = [] {
    constexpr std::array<bool, sizeof...(Args) + 1> policies{ModelPolicy<Args>..., false};
    std::size_t n = 0;
    while (policies[n])
    {
        ++n;
    }
    return n;
}();

template <typename... Args>
inline constexpr bool policyAfterInput = [] {
    constexpr std::array<bool, sizeof...(Args) + 1> policies{ModelPolicy<Args>..., false};
    for (std::size_t i = leadingPolicies<Args...>; i < sizeof...(Args); ++i)
    {
        if (policies[i])
        {
            return true;
        }
    }
    return false;
}();

template <typename... Ts>
struct FirstTypeOf
{
    using type = void;  
};

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


template <bool Uncertain, typename A>
struct UncertainScalarOf
{
    using type = void;  
};

template <typename A>
struct UncertainScalarOf<true, A>
{
    using type = typename A::Scalar;
};


template <typename Model, typename... Inputs>
[[nodiscard]] consteval auto
transformReadsPolicy()
{
    if constexpr (readsWholeState<Model>)
    {
        static_assert(
            sizeof...(Inputs) == 1 && (ArgTraits<Inputs>::isUncertain && ...),
            "transform() of a function without reads() takes one uncertain input, the vector the function reads; "
            "name several inputs with reads(), each over a tagged space"
        );
        return std::type_identity<ReadsWhole<TransformInputTag<std::tuple_element_t<0, std::tuple<Inputs...>>>>>{};
    }
    else
    {
        return std::type_identity<typename Model::template Policy<PolicyKind::reads>>{};
    }
}

template <typename Model>
[[nodiscard]] consteval bool
transformOutputFromNoise()
{
    if constexpr (Model::template has<PolicyKind::outputSpace> || Model::template has<PolicyKind::writes>)
    {
        return false;
    }
    else if constexpr (Model::template has<PolicyKind::additiveNoise>)
    {
        return Model::template Policy<PolicyKind::additiveNoise>::HasOwnSpace;
    }
    else
    {
        return false;
    }
}

template <typename Model, typename In, bool OutputFromNoise>
[[nodiscard]] consteval auto
transformOutputSpace()
{
    using Function = typename Model::FunctionType;
    if constexpr (Model::template has<PolicyKind::outputSpace>)
    {
        using Policy = typename Model::template Policy<PolicyKind::outputSpace>;
        return std::type_identity<decltype(std::declval<Policy const &>().in(std::declval<In const &>()))>{};
    }
    else if constexpr (Model::template has<PolicyKind::writes>)
    {
        return std::type_identity<typename Model::template Policy<PolicyKind::writes>::OutputSpace>{};
    }
    else if constexpr (OutputFromNoise)
    {
        return std::type_identity<
            typename Model::template Policy<PolicyKind::additiveNoise>::Distribution::SpaceType>{};
    }
    else if constexpr (requires { typename Function::OutputSpaceType; })
    {
        return std::type_identity<typename Function::OutputSpaceType>{};
    }
    else
    {
        using Returned = decltype(vectorReturned(callForm(
            std::declval<Function const &>(),
            std::declval<In const &>(),
            std::declval<NoOutputSpace const &>(),
            std::declval<typename Moments<In, double>::VectorType const &>()
        )));
        using Evaluated = std::remove_cvref_t<decltype(valueOf(std::declval<Returned>()).eval())>;
        static_assert(Evaluated::ColsAtCompileTime == 1, "a transformed function must return a column vector");
        return std::type_identity<Linear<Evaluated::RowsAtCompileTime, Evaluated::MaxRowsAtCompileTime>>{};
    }
}

template <typename Model, typename... Inputs>
class ModelTransform
{
public:

    static constexpr bool Whole = readsWholeState<Model>;
 
    static constexpr std::size_t InputCount = sizeof...(Inputs);  

    template <std::size_t J>
    using InputAt = std::tuple_element_t<J, std::tuple<Inputs...>>;

    template <std::size_t J>
    using TagAt = TransformInputTag<InputAt<J>>;
 
private:

    [[nodiscard]] static consteval bool
    uncertain(std::size_t j)
    {
        constexpr std::array<bool, InputCount + 1> flags{ArgTraits<Inputs>::isUncertain..., false};
        return flags[j];
    }
 
    static_assert(
        ((!std::is_void_v<TransformInputTag<Inputs>>) && ...),
        "an input of transform() is an uncertain input, one marked by withCross(), or sure<Tag>(value), after a "
        "linearizedAbout() reference if there is one; its output space is a policy, outputSpace(space), before the "
        "inputs"
    );
    static_assert(
        distinctTypes<TransformInputTag<Inputs>...>(),
        "two inputs of transform() have one tag; an input's tag names it once"
    );

    using ReadsPolicy = typename decltype(transformReadsPolicy<Model, Inputs...>())::type;
 
public:
 
    static constexpr std::size_t ItemCount = ReadsPolicy::Count;  

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

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

    using In = typename ReadsPolicy::InputSpace;

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

    template <std::size_t I>
    [[nodiscard]] static consteval std::size_t
    inputOfItem()
    {
        if constexpr (Whole)
        {
            return 0;
        }
        else
        {
            using Path = typename Item<I>::PathType;
            static_assert(
                std::is_same_v<Path, TagPath<PathTag<Path>>>,
                "transform() matches each item of reads() to an input by its one tag, path<Tag> or "
                "on<Name>(path<Tag>); "
                "it has no state for a path of several tags to lie in"
            );
            constexpr std::array<bool, InputCount + 1> matches{
                std::is_same_v<TransformInputTag<Inputs>, PathTag<Path>>...,
                false
            };
            std::size_t j = 0;
            while (j < InputCount && !matches[j])
            {
                ++j;
            }
            return j;
        }
    }

    [[nodiscard]] static consteval std::size_t
    inputOf(std::size_t i)
    {
        return []<std::size_t... I>(std::size_t item, std::index_sequence<I...>) {
            std::size_t found = InputCount;
            ((I == item ? (found = inputOfItem<I>(), 0) : 0), ...);
            return found;
        }(i, std::make_index_sequence<ItemCount>{});
    }

    template <std::size_t... I>
    [[nodiscard]] static consteval bool
    itemsGiven(std::index_sequence<I...>)
    {
        return ((inputOfItem<I>() < InputCount) && ...);
    }
 
public:

    static constexpr bool ItemsGiven = itemsGiven(std::make_index_sequence<ItemCount>{});
 
private:

    [[nodiscard]] static consteval bool
    inputRead(std::size_t j)
    {
        for (std::size_t i = 0; i < ItemCount; ++i)
        {
            if (inputOf(i) == j)
            {
                return true;
            }
        }
        return false;
    }

    template <std::size_t... J>
    [[nodiscard]] static consteval bool
    inputsRead(std::index_sequence<J...>)
    {
        return (inputRead(J) && ...);
    }
 
public:

    static constexpr bool InputsRead = inputsRead(std::make_index_sequence<InputCount>{});
 
private:

    [[nodiscard]] static consteval std::size_t
    uncertainIndex(std::size_t j)
    {
        std::size_t k = 0;
        for (std::size_t i = 0; i < j; ++i)
        {
            if (uncertain(i))
            {
                ++k;
            }
        }
        return k;
    }

    static constexpr std::size_t UncertainCount = (std::size_t{ArgTraits<Inputs>::isUncertain} + ... + 0);
 
public:

    using Function = typename Model::FunctionType;

    static constexpr bool OutputFromNoise = transformOutputFromNoise<Model>();

    using Out = typename decltype(transformOutputSpace<Model, In, OutputFromNoise>())::type;

    explicit ModelTransform(Model const & model, Inputs const &... inputs)
        : m_model{model}
        , m_inputs{inputs...}
        , m_in{makeIn()}
        , m_out{makeOut()}
    {
        if constexpr (Whole)
        {
            m_valueRows[0] = BlockRows{0, valueDimension(m_in)};
            m_differenceRows[0] = BlockRows{0, differenceDimension(m_in)};
        }
        else
        {
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                ((m_valueRows[I] = blockRows<Name<I>>(m_in), m_differenceRows[I] = differenceBlockRows<Name<I>>(m_in)),
                 ...);
            }(std::make_index_sequence<ItemCount>{});
        }
    }

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

    template <typename... A>
    [[nodiscard]] auto
    operator()(A const &... args) const
    {
        static_assert(sizeof...(A) == InputCount);
        using T = typename FirstTypeOf<typename UncertainScalarOf<ArgTraits<Inputs>::isUncertain, A>::type...>::type;
        auto const arguments = std::tuple<A const &...>{args...};
        auto x = InVector<T>{};
        x.resize(valueDimension(m_in));
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (gatherItem<I, T>(x, std::get<inputOfItem<I>()>(arguments)), ...);
        }(std::make_index_sequence<ItemCount>{});
        return split(vectorReturned(callForm(m_model.function(), m_in, m_out, x)));
    }

    static constexpr bool EvaluatesBatches = takesBatches<Function>;

    template <typename Y, typename... A>
    void
    evaluateBatch(Y & y, A const &... args) const
        requires EvaluatesBatches
    {
        callBatched(m_model.function(), m_in, m_out, gatherColumns<Y>(y.cols(), args...), y);
    }

    template <typename Y, typename J, typename... A>
    void
    evaluateBatchWithAffineJacobians(Y & y, J & jacobians, A const &... args) const
        requires EvaluatesBatches
    {
        auto wanted = std::array<bool, std::tuple_size_v<InputBlocks<In>>>{};
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((wanted[I] = itemFactored<I>()), ...);
        }(std::make_index_sequence<ItemCount>{});
        auto const x = gatherColumns<Y>(y.cols(), args...);
        callBatchedWithJacobians(m_model.function(), m_in, m_out, x, wanted, y, jacobians);
    }

    template <typename Y, typename... A>
    [[nodiscard]] auto
    gatherColumns(Eigen::Index count, A const &... args) const
    {
        static_assert(sizeof...(A) == InputCount);
        using T = typename Y::Scalar;
        using X = MatrixOf<
            T,
            valueDimensionAtCompileTime<In>,
            Y::ColsAtCompileTime,
            maxValueDimensionAtCompileTime<In>,
            Y::MaxColsAtCompileTime>;
        auto const arguments = std::tuple<A const &...>{args...};
        auto x = X{};
        x.resize(valueDimension(m_in), count);
        auto column = InVector<T>{};
        column.resize(valueDimension(m_in));
        for (Eigen::Index j = 0; j < count; ++j)
        {
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                (gatherItem<I, T>(column, columnOf<inputOfItem<I>()>(std::get<inputOfItem<I>()>(arguments), j)), ...);
            }(std::make_index_sequence<ItemCount>{});
            x.col(j) = column;
        }
        return x;
    }

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

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

    [[nodiscard]] static consteval int
    affineSize()
    {
        return []<std::size_t... I>(std::index_sequence<I...>) {
            return sizeSum<(itemFactored<I>() ? valueDimensionAtCompileTime<DifferenceSpace<Name<I>>> : 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<Name<I>>> : 0)...>;
        }(std::make_index_sequence<ItemCount>{});
    }

    [[nodiscard]] IndexVector<affineSize(), affineMaxSize()>
    affineCoordinates() const
    {
        std::array<Eigen::Index, InputCount + 1> starts{};
        [&]<std::size_t... J>(std::index_sequence<J...>) {
            ((starts[J + 1] = starts[J] + (uncertain(J) ? argInput(std::get<J>(m_inputs)).differenceDimension() : 0)),
             ...);
        }(std::make_index_sequence<InputCount>{});
        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>())
                    {
                        for (Eigen::Index r = 0; r < m_differenceRows[I].size; ++r)
                        {
                            result(next++) = starts[inputOfItem<I>()] + r;
                        }
                    }
                }(),
                ...
            );
        }(std::make_index_sequence<ItemCount>{});
        return result;
    }

    template <typename Result>
    [[nodiscard]] auto
    withNoise(Result && result) const
    {
        if constexpr (!Model::template has<PolicyKind::additiveNoise>)
        {
            return std::forward<Result>(result);
        }
        else
        {
            auto const placed = placeNoise(m_out, m_model.template policy<PolicyKind::additiveNoise>());
            auto const & noise = placed.distribution();
            static_assert(
                !isGaussianMixture<std::remove_cvref_t<decltype(noise)>>,
                "transform() adds one Gaussian noise; mixture noise gives a branch per component, which update() keeps"
            );
            return plusNoise(result, noise);
        }
    }
 
private:

    [[nodiscard]] In
    makeIn() const
    {
        if constexpr (Whole)
        {
            return argInput(std::get<0>(m_inputs)).space();
        }
        else
        {
            return [&]<std::size_t... I>(std::index_sequence<I...>) {
                return In{nameOf<Name<I>>(spaceOf<inputOfItem<I>()>())...};
            }(std::make_index_sequence<ItemCount>{});
        }
    }

    template <std::size_t J>
    [[nodiscard]] auto
    spaceOf() const
    {
        auto const & input = std::get<J>(m_inputs);
        if constexpr (isSure<std::remove_cvref_t<decltype(input)>>)
        {
            return input.space;
        }
        else
        {
            return argInput(input).space();
        }
    }

    [[nodiscard]] Out
    makeOut() const
    {
        if constexpr (Model::template has<PolicyKind::outputSpace>)
        {
            return m_model.template policy<PolicyKind::outputSpace>().in(m_in);
        }
        else if constexpr (OutputFromNoise)
        {
            return m_model.template policy<PolicyKind::additiveNoise>().distribution().space();
        }
        else if constexpr (!FixedSizeSpace<Out> && requires(Function const & f) { f.outputRows(); })
        {
            return Out{m_model.function().outputRows()};
        }
        else
        {
            return Out{};
        }
    }

    template <std::size_t I, typename T, typename Argument>
    void
    gatherItem(InVector<T> & x, Argument const & argument) const
    {
        constexpr int Rows = valueDimensionAtCompileTime<Name<I>>;
        auto target = middleRowsOf<Rows>(x, m_valueRows[I].start, m_valueRows[I].size);
        if constexpr (uncertain(inputOfItem<I>()))
        {
            target = argument;
        }
        else
        {
            target = argument.value.template cast<T>();
        }
    }

    template <std::size_t J, typename Argument>
    [[nodiscard]] static decltype(auto)
    columnOf(Argument const & argument, Eigen::Index j)
    {
        if constexpr (uncertain(J))
        {
            return argument.col(j);
        }
        else
        {
            return (argument);
        }
    }

    template <typename R>
    [[nodiscard]] auto
    split(R && returned) const
    {
        using Plain = std::remove_cvref_t<R>;
        if constexpr (!isWithJacobian<Plain> || Whole)
        {
            return std::forward<R>(returned);
        }
        else
        {
            using J = std::remove_cvref_t<decltype(returned.jacobian)>;
            static_assert(
                MatrixExpression<J>,
                "a transformed function returns its Jacobian as one matrix over its input's differences, "
                "withJacobian(y, J)"
            );
            auto const & jacobian = returned.jacobian;
            EMRYS_FILTER_ASSERT(
                jacobian.cols() == differenceDimension(m_in)
                && "a returned Jacobian has a column per difference coordinate of the function's input"
            );
            auto blocks = [&]<std::size_t... J_>(std::index_sequence<J_...>) {
                return std::tuple_cat(inputBlock<J_>(jacobian)...);
            }(std::make_index_sequence<InputCount>{});
            if constexpr (UncertainCount == 1)
            {
                return withJacobian(std::move(returned.value), std::move(std::get<0>(blocks)));
            }
            else
            {
                return withJacobian(std::move(returned.value), std::move(blocks));
            }
        }
    }

    template <std::size_t J, typename Jacobian>
    [[nodiscard]] auto
    inputBlock(Jacobian const & jacobian) const
    {
        if constexpr (!uncertain(J))
        {
            return std::tuple<>{};
        }
        else
        {
            using Traits = ArgTraits<InputAt<J>>;
            using Block = MatrixOf<
                typename Jacobian::Scalar,
                Jacobian::RowsAtCompileTime,
                Traits::differenceSize,
                Jacobian::MaxRowsAtCompileTime,
                Traits::differenceMaxSize>;
            Block block = Block::Zero(jacobian.rows(), argInput(std::get<J>(m_inputs)).differenceDimension());
            [&]<std::size_t... I>(std::index_sequence<I...>) {
                (
                    [&] {
                        if constexpr (inputOfItem<I>() == J)
                        {
                            block += jacobian.middleCols(m_differenceRows[I].start, m_differenceRows[I].size);
                        }
                    }(),
                    ...
                );
            }(std::make_index_sequence<ItemCount>{});
            return std::tuple{std::move(block)};
        }
    }
 
    Model const & m_model;                                
    std::tuple<Inputs const &...> m_inputs;               
    In m_in;                                              
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Out m_out;             
    std::array<BlockRows, ItemCount> m_valueRows{};       
    std::array<BlockRows, ItemCount> m_differenceRows{};  
};

template <typename Space, typename Scalar, typename Noise>
[[nodiscard]] GaussianParams<Space, Scalar>
plusNoise(GaussianParams<Space, Scalar> const & result, Noise const & noise)
{
    using Matrix = typename GaussianParams<Space, Scalar>::MatrixType;
    Matrix cov = symmetricFromLower<Matrix>(result.covarianceLower());
    cov += symmetricFromLower<Matrix>(noise.covarianceLower().template cast<Scalar>());
    return GaussianParams<Space, Scalar>{
        result.space(),
        Emrys::plus(result.space(), result.mean(), noise.mean().template cast<Scalar>()),
        std::move(cov)
    };
}

template <typename InputSpace, typename Space, typename Form, typename Scalar, typename Noise>
[[nodiscard]] OutputWithCross<InputSpace, Space, Form, Scalar>
plusNoise(OutputWithCross<InputSpace, Space, Form, Scalar> const & result, Noise const & noise)
{
    using Result = OutputWithCross<InputSpace, Space, Form, Scalar>;
    using Matrix = typename Result::OutputMatrix;
    Matrix inputs = result.inputsCovariance();
    inputs += symmetricFromLower<Matrix>(noise.covarianceLower().template cast<Scalar>());
    auto const & space = result.output().space();
    typename Result::OutputVector mean =
        Emrys::plus(space, result.output().mean(), noise.mean().template cast<Scalar>());
    if constexpr (std::is_same_v<Form, RootForm>)
    {
        return Result{
            space,
            std::move(mean),
            typename Result::RootMatrix{result.root()},
            typename Result::CrossMatrix{result.whitenedCross()},
            Matrix{result.explainedCovariance()},
            std::move(inputs),
            Matrix{result.unexplainedCovariance()}
        };
    }
    else
    {
        return Result{
            space,
            std::move(mean),
            typename Result::JacobianMatrix{result.jacobian()},
            typename Result::CrossMatrix{result.crossCovariance()},
            Matrix{result.explainedCovariance()},
            std::move(inputs),
            Matrix{result.unexplainedCovariance()}
        };
    }
}

template <typename G>
[[nodiscard]] consteval bool
inputsMatch()
{
    static_assert(
        G::ItemsGiven,
        "the function reads an input that transform() was not given: give each item of reads() as an uncertain "
        "input over its tag's space, or as sure<Tag>(value)"
    );
    static_assert(
        G::InputsRead,
        "transform() was given an input the function does not read: reads() names each input by the tag of its "
        "space, and an untagged input is read only by a function without reads()"
    );
    return G::ItemsGiven && G::InputsRead;
}

template <typename Method, typename Model, typename... Inputs>
[[nodiscard]] auto
transformModel(Method const & method, Model const & model, Inputs const &... inputs)
{
    static_assert(
        Model::HasFunction,
        "a model known by its likelihood alone has no function for transform() to push its inputs through"
    );
    if constexpr (leadsWithReference<Inputs...>)
    {
        return [&]<typename About, typename... Rest>(About const & about, Rest const &... rest) {
            using G = ModelTransform<Model, Rest...>;
            if constexpr (inputsMatch<G>())
            {
                auto const g = G{model, rest...};
                return g.withNoise(transformBy(method, g, g.out(), about, rest...));
            }
        }(inputs...);
    }
    else if constexpr (!inputsMatch<ModelTransform<Model, Inputs...>>())
    {
        return;
    }
    else
    {
        using G = ModelTransform<Model, Inputs...>;
        auto const g = G{model, inputs...};
        if constexpr (isSigmaPointMethod<Method> && G::factorsAffine())
        {
            auto const factored =
                FactoredAt<G::affineSize(), G::affineMaxSize(), Method>{method, g.affineCoordinates()};
            return g.withNoise(transformBy(factored, g, g.out(), inputs...));
        }
        else
        {
            return g.withNoise(transformBy(method, g, g.out(), inputs...));
        }
    }
}
 
}  // namespace internal

template <typename Method, typename F, typename... Args>
    requires(!MeasurementModel<F> && !internal::isModel<F> && internal::crossCount<Args...> <= 1)
[[nodiscard]] auto
transform(Method const & method, F const & f, Args const &... args)
{
    static_assert(
        !internal::policyAfterInput<Args...>,
        "transform() takes its policies before its inputs: transform(method, f, reads(...), additiveNoise(...), "
        "inputs...)"
    );
    constexpr std::size_t Policies = internal::leadingPolicies<Args...>;
    auto const all = std::tuple<Args const &...>{args...};
    return [&]<std::size_t... P, std::size_t... I>(std::index_sequence<P...>, std::index_sequence<I...>) {
        auto const model = measurementModel(f, std::get<P>(all)...);
        return internal::transformModel(method, model, std::get<Policies + I>(all)...);
    }(std::make_index_sequence<Policies>{}, std::make_index_sequence<sizeof...(Args) - Policies>{});
}

template <typename Method, typename Model, typename... Inputs>
    requires((MeasurementModel<Model> || internal::isModel<Model>) && internal::crossCount<Inputs...> <= 1)
[[nodiscard]] auto
transform(Method const & method, Model const & model, Inputs const &... inputs)
{
    static_assert(
        MeasurementModel<Model>,
        "transform() pushes its inputs through a measurement model or a function; a process model's step is "
        "predict()'s"
    );
    static_assert(
        internal::leadingPolicies<Inputs...> == 0,
        "transform() of a model takes the model's own policies, so give the inputs alone after it"
    );
    return internal::transformModel(method, model, inputs...);
}
 
}  // namespace Emrys