Roots.h

Go to the documentation of this file.

#pragma once

 
#include <concepts>
#include <cstddef>
#include <expected>
#include <optional>
#include <tuple>
#include <type_traits>
#include <utility>
#include <variant>
 
#include <Eigen/Cholesky>
#include <Eigen/Eigenvalues>
#include <Eigen/SVD>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/NumericError.h>
#include <EmrysFilter/Core/Triangularization.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Derived>
using SquareMatrixOf = Eigen::Matrix<
    typename Derived::Scalar,
    Derived::RowsAtCompileTime,
    Derived::RowsAtCompileTime,
    Eigen::ColMajor,
    Derived::MaxRowsAtCompileTime,
    Derived::MaxRowsAtCompileTime>;

template <typename Matrix>
using RootColumnOf = Eigen::
    Matrix<typename Matrix::Scalar, Matrix::RowsAtCompileTime, 1, Eigen::ColMajor, Matrix::MaxRowsAtCompileTime, 1>;

template <bool CanFail, typename Matrix>
using RootResult = std::conditional_t<CanFail, std::expected<Matrix, GaussianError>, Matrix>;

struct DiscardReport
{
    template <typename Report>
    constexpr void
    operator()(Report const & /*report*/) const noexcept
    {}
};

template <typename Scalar, double Value>
[[nodiscard]] Scalar
scalarFromConstant()
{
    using Literal = typename Eigen::NumTraits<Scalar>::Literal;
    return Scalar(static_cast<Literal>(Value));
}

template <typename Derived>
[[nodiscard]] bool
lowerFinite(Eigen::MatrixBase<Derived> const & a)
{
    for (Eigen::Index j = 0; j < a.cols(); ++j)
    {
        if (!allFinite(a.col(j).tail(a.rows() - j)))
        {
            return false;
        }
    }
    return true;
}

template <double Relative = 1.0, double Absolute = 0.0, typename Derived>
[[nodiscard]] std::expected<typename Derived::PlainObject, GaussianError>
ridgedCovariance(Eigen::MatrixBase<Derived> const & cov, typename Derived::Scalar * delta = nullptr)
{
    using Scalar = typename Derived::Scalar;
    EMRYS_FILTER_ASSERT(cov.rows() == cov.cols() && "a covariance must be square");
    Eigen::Index const n = cov.rows();
 
    // Row i's off-diagonal elements are row i left of the diagonal and, by symmetry, column i
    // below it. NaN or infinity anywhere in the lower triangle reaches some row's bound.
    auto const diagonal = cov.diagonal();
    Scalar bound = Eigen::NumTraits<Scalar>::infinity();
    for (Eigen::Index i = 0; i < n; ++i)
    {
        Scalar const radius = cov.row(i).head(i).cwiseAbs().sum() + cov.col(i).tail(n - i - 1).cwiseAbs().sum();
        Scalar const rowBound = diagonal(i) - radius;
        if (!isFinite(rowBound))
        {
            return std::unexpected{GaussianError::NonFinite};
        }
        bound = Eigen::numext::mini(bound, rowBound);
    }
 
    Scalar const shift = Eigen::numext::maxi(-bound, Scalar(0));
    Scalar const scale = Eigen::numext::maxi(n > 0 ? diagonal.cwiseAbs().maxCoeff() : Scalar(0), shift);
    Scalar const margin =
        scalarFromConstant<Scalar, Relative>() * Eigen::numext::sqrt(Eigen::NumTraits<Scalar>::epsilon()) * scale
        + scalarFromConstant<Scalar, Absolute>();
    typename Derived::PlainObject ridged = cov;
    ridged.diagonal().array() += shift + margin;
    if (delta != nullptr)
    {
        *delta = shift + margin;
    }
    return ridged;
}

template <typename Matrix>
void
fixColumnSigns(Matrix & root)
{
    for (Eigen::Index j = 0; j < root.cols(); ++j)
    {
        Eigen::Index largest = 0;
        root.col(j).cwiseAbs().maxCoeff(&largest);
        if (root(largest, j) < typename Matrix::Scalar(0))
        {
            root.col(j) *= typename Matrix::Scalar(-1);
        }
    }
}

template <typename Derived>
[[nodiscard]] bool
isTriangularRoot(Eigen::MatrixBase<Derived> const & a)
{
    using Scalar = typename Derived::Scalar;
    if (a.rows() != a.cols())
    {
        return false;
    }
    for (Eigen::Index j = 0; j < a.cols(); ++j)
    {
        if (a(j, j) < Scalar(0))
        {
            return false;
        }
        for (Eigen::Index i = 0; i < j; ++i)
        {
            if (a(i, j) != Scalar(0))
            {
                return false;
            }
        }
    }
    return true;
}

template <typename Derived>
[[nodiscard]] SquareMatrixOf<Derived>
triangularRootOf(Eigen::MatrixBase<Derived> const & factor)
{
    if constexpr (Derived::ColsAtCompileTime == Eigen::Dynamic
                  || Derived::ColsAtCompileTime == Derived::RowsAtCompileTime)
    {
        if (isTriangularRoot(factor))
        {
            return SquareMatrixOf<Derived>{factor};
        }
    }
    return lowerTriangularRoot(factor);
}

template <typename Derived>
struct PrincipalAxes
{
    SquareMatrixOf<Derived> vectors;              
    RootColumnOf<SquareMatrixOf<Derived>> roots;  
    Eigen::Index clamped = 0;                     
};

template <typename Derived>
[[nodiscard]] PrincipalAxes<Derived>
principalAxesOfCovariance(Eigen::MatrixBase<Derived> const & cov)
{
    using Matrix = SquareMatrixOf<Derived>;
    using Scalar = typename Matrix::Scalar;
    auto const solver = Eigen::SelfAdjointEigenSolver<Matrix>{Matrix{cov}, Eigen::ComputeEigenvectors};
    Eigen::Index const n = cov.rows();
    auto axes = PrincipalAxes<Derived>{};
    axes.vectors.resize(n, n);
    axes.roots.resize(n);
    for (Eigen::Index j = 0; j < n; ++j)
    {
        // The solver sorts its eigenvalues smallest first.
        Scalar const value = solver.eigenvalues()(n - 1 - j);
        if (value < Scalar(0))
        {
            ++axes.clamped;
        }
        axes.vectors.col(j) = solver.eigenvectors().col(n - 1 - j);
        axes.roots(j) = Eigen::numext::sqrt(Eigen::numext::maxi(value, Scalar(0)));
    }
    return axes;
}

template <typename Derived>
[[nodiscard]] PrincipalAxes<Derived>
principalAxesOfFactor(Eigen::MatrixBase<Derived> const & factor)
{
    using Matrix = SquareMatrixOf<Derived>;
    Matrix const lower = triangularRootOf(factor);
    auto const svd = Eigen::JacobiSVD<Matrix, Eigen::ComputeFullU>{lower};
    return PrincipalAxes<Derived>{svd.matrixU(), svd.singularValues(), 0};
}

template <typename T>
inline constexpr bool isFallbackList = false;
 
}  // namespace internal

template <typename T>
concept Axes = requires {
    { std::bool_constant<T::canFail>{} } -> std::convertible_to<bool>;
} && requires(T const & axes, Eigen::Matrix3d const & cov) {
    { axes.root(cov) } -> std::same_as<internal::RootResult<T::canFail, Eigen::Matrix3d>>;
};

template <typename T>
concept Repair = requires(T const & repair, Eigen::Matrix3d const & cov) {
    { repair.repaired(cov) } -> std::same_as<std::expected<Eigen::Matrix3d, GaussianError>>;
};
 
namespace internal
{

struct NoRepair
{};
 
}  // namespace internal

namespace roots
{

template <typename Derived>
using RootMatrix = internal::SquareMatrixOf<Derived>;

struct NoReport
{};

struct ClampReport
{
    Eigen::Index clamped = 0;  
};

template <typename Scalar>
struct RidgeReport
{
    Scalar ridge = 0;  
};

struct Cholesky
{
    static constexpr bool canFail = true;  

    template <typename Derived>
    [[nodiscard]] std::expected<internal::SquareMatrixOf<Derived>, GaussianError>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        using Matrix = internal::SquareMatrixOf<Derived>;
        if (!internal::lowerFinite(cov))
        {
            return std::unexpected{GaussianError::NonFinite};
        }
        auto const llt = Eigen::LLT<Matrix>{Matrix{cov}};
        if (llt.info() != Eigen::Success)
        {
            return std::unexpected{GaussianError::NotPositiveDefinite};
        }
        return Matrix{llt.matrixL()};
    }

    template <typename Derived>
    [[nodiscard]] std::expected<internal::SquareMatrixOf<Derived>, GaussianError>
    rootFromFactor(Eigen::MatrixBase<Derived> const & factor) const
    {
        if (!allFinite(factor))
        {
            return std::unexpected{GaussianError::NonFinite};
        }
        return internal::triangularRootOf(factor);
    }
};

struct Pivoted
{
    static constexpr bool canFail = false;  

    template <typename Scalar>
    using Report = ClampReport;

    template <typename Derived>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        return root(cov, internal::DiscardReport{});
    }

    template <typename Derived, typename Sink>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov, Sink && sink) const
    {
        using Matrix = internal::SquareMatrixOf<Derived>;
        using Scalar = typename Matrix::Scalar;
        auto const ldlt = Eigen::LDLT<Matrix>{Matrix{cov}};
        auto const pivots = ldlt.vectorD();
        sink(ClampReport{(pivots.array() < Scalar(0)).count()});
        auto const d = pivots.cwiseMax(Scalar(0)).cwiseSqrt().eval();
        Matrix lowerD = ldlt.matrixL();
        lowerD *= d.asDiagonal();
        return ldlt.transpositionsP().transpose() * lowerD;
    }
};

struct Principal
{
    static constexpr bool canFail = false;  

    template <typename Scalar>
    using Report = ClampReport;

    template <typename Derived>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        return root(cov, internal::DiscardReport{});
    }

    template <typename Derived, typename Sink>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov, Sink && sink) const
    {
        auto const axes = internal::principalAxesOfCovariance(cov);
        sink(ClampReport{axes.clamped});
        return fromAxes<Derived>(axes);
    }

    template <typename Derived>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    rootFromFactor(Eigen::MatrixBase<Derived> const & factor) const
    {
        return fromAxes<Derived>(internal::principalAxesOfFactor(factor));
    }
 
private:

    template <typename Derived>
    [[nodiscard]] static internal::SquareMatrixOf<Derived>
    fromAxes(internal::PrincipalAxes<Derived> const & axes)
    {
        internal::SquareMatrixOf<Derived> root = axes.vectors * axes.roots.asDiagonal();
        internal::fixColumnSigns(root);
        return root;
    }
};

struct Symmetric
{
    static constexpr bool canFail = false;  

    template <typename Scalar>
    using Report = ClampReport;

    template <typename Derived>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        return root(cov, internal::DiscardReport{});
    }

    template <typename Derived, typename Sink>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    root(Eigen::MatrixBase<Derived> const & cov, Sink && sink) const
    {
        auto const axes = internal::principalAxesOfCovariance(cov);
        sink(ClampReport{axes.clamped});
        return fromAxes<Derived>(axes);
    }

    template <typename Derived>
    [[nodiscard]] internal::SquareMatrixOf<Derived>
    rootFromFactor(Eigen::MatrixBase<Derived> const & factor) const
    {
        return fromAxes<Derived>(internal::principalAxesOfFactor(factor));
    }
 
private:

    template <typename Derived>
    [[nodiscard]] static internal::SquareMatrixOf<Derived>
    fromAxes(internal::PrincipalAxes<Derived> const & axes)
    {
        return axes.vectors * axes.roots.asDiagonal() * axes.vectors.transpose();
    }
};

template <typename Vector>
struct Along
{
    static constexpr bool canFail = true;  
 
    Vector direction;  

    template <typename Derived>
    [[nodiscard]] std::expected<internal::SquareMatrixOf<Derived>, GaussianError>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        auto const lower = Cholesky{}.root(cov);
        if (!lower)
        {
            return std::unexpected{lower.error()};
        }
        return rotated(*lower);
    }

    template <typename Derived>
    [[nodiscard]] std::expected<internal::SquareMatrixOf<Derived>, GaussianError>
    rootFromFactor(Eigen::MatrixBase<Derived> const & factor) const
    {
        auto const lower = Cholesky{}.rootFromFactor(factor);
        if (!lower)
        {
            return std::unexpected{lower.error()};
        }
        return rotated(*lower);
    }
 
private:

    template <typename Matrix>
    [[nodiscard]] std::expected<Matrix, GaussianError>
    rotated(Matrix const & lower) const
    {
        using Scalar = typename Matrix::Scalar;
        using ColumnType = internal::RootColumnOf<Matrix>;
        EMRYS_FILTER_ASSERT(
            direction.size() == lower.rows() && "a direction has a coordinate per coordinate of the covariance"
        );
        ColumnType const whitened =
            lower.template triangularView<Eigen::Lower>().solve(direction.template cast<Scalar>());
        Scalar const length = whitened.norm();
        if (!isFinite(length) || !(length > Scalar(0)))
        {
            return std::unexpected{GaussianError::NotPositiveDefinite};
        }
        ColumnType const unit = whitened / length;
        // The reflection H = I - 2 v v^T / v^T v with v = e1 - s u takes e1 to s u. The sign s is
        // the one that makes v^T v at least 2, so that no cancellation reaches it, and the first
        // column is then flipped by s so that it points along d.
        Scalar const sign = unit(0) < Scalar(0) ? Scalar(1) : Scalar(-1);
        ColumnType v = -sign * unit;
        v(0) += Scalar(1);
        Matrix root = lower - (Scalar(2) / v.squaredNorm()) * (lower * v) * v.transpose();
        root.col(0) *= sign;
        return root;
    }
};

template <double Relative = 1.0, double Absolute = 0.0>
struct Ridge
{
    static_assert(Relative >= 0 && Absolute >= 0, "a ridge's margins are at least zero");
 
    static constexpr double relative = Relative;  
    static constexpr double absolute = Absolute;  

    template <typename Scalar>
    using Report = RidgeReport<Scalar>;

    template <typename Derived>
    [[nodiscard]] std::expected<typename Derived::PlainObject, GaussianError>
    repaired(Eigen::MatrixBase<Derived> const & cov) const
    {
        return internal::ridgedCovariance<Relative, Absolute>(cov);
    }

    template <typename Derived, typename Sink>
    [[nodiscard]] std::expected<typename Derived::PlainObject, GaussianError>
    repaired(Eigen::MatrixBase<Derived> const & cov, Sink && sink) const
    {
        auto delta = typename Derived::Scalar(0);
        auto result = internal::ridgedCovariance<Relative, Absolute>(cov, &delta);
        sink(RidgeReport<typename Derived::Scalar>{delta});
        return result;
    }
};
 
}  // namespace roots
 
namespace internal
{

template <typename Entry, typename Scalar>
struct ReportOf
{
    using type = roots::NoReport;  
};

template <typename Entry, typename Scalar>
    requires requires { typename Entry::template Report<Scalar>; }
struct ReportOf<Entry, Scalar>
{
    using type = typename Entry::template Report<Scalar>;
};


template <typename Entry, typename Scalar>
using ReportOfT = typename ReportOf<Entry, Scalar>::type;

template <typename A, typename Scalar>
struct FailureReportOf
{
    using type = ReportOfT<A, Scalar>;  
};

template <typename A, typename Scalar>
    requires requires { typename A::template FailureReport<Scalar>; }
struct FailureReportOf<A, Scalar>
{
    using type = typename A::template FailureReport<Scalar>;
};


template <typename A, typename Scalar>
using FailureReportOfT = typename FailureReportOf<A, Scalar>::type;

template <typename Entry>
[[nodiscard]] consteval bool
entryCanFail()
{
    if constexpr (requires { Entry::canFail; })
    {
        return Entry::canFail;
    }
    else
    {
        return true;
    }
}

template <typename... Fallbacks>
[[nodiscard]] consteval bool
fallbacksReached()
{
    constexpr bool canFail[] = {entryCanFail<Fallbacks>()..., true};
    for (std::size_t i = 0; i + 1 < sizeof...(Fallbacks); ++i)
    {
        if (!canFail[i])
        {
            return false;
        }
    }
    return true;
}

template <typename A, typename Derived, typename Sink>
[[nodiscard]] auto
rootWithSink(A const & axes, Eigen::MatrixBase<Derived> const & cov, Sink && sink)
{
    if constexpr (requires { axes.root(cov, sink); })
    {
        return axes.root(cov, sink);
    }
    else
    {
        sink(roots::NoReport{});
        return axes.root(cov);
    }
}

template <typename R, typename Derived, typename Sink>
[[nodiscard]] auto
repairWithSink(R const & repair, Eigen::MatrixBase<Derived> const & cov, Sink && sink)
{
    if constexpr (requires { repair.repaired(cov, sink); })
    {
        return repair.repaired(cov, sink);
    }
    else
    {
        sink(roots::NoReport{});
        return repair.repaired(cov);
    }
}

template <typename A, typename Derived>
[[nodiscard]] auto
rootFromFactor(A const & axes, Eigen::MatrixBase<Derived> const & factor)
{
    if constexpr (requires { axes.rootFromFactor(factor); })
    {
        return axes.rootFromFactor(factor);
    }
    else
    {
        return axes.root(SquareMatrixOf<Derived>{factor * factor.transpose()});
    }
}

template <typename Matrix, typename Result>
[[nodiscard]] std::expected<Matrix, GaussianError>
asExpected(Result && result)
{
    if constexpr (requires { result.has_value(); })
    {
        if (!result)
        {
            return std::unexpected{result.error()};
        }
        return Matrix{std::move(*result)};
    }
    else
    {
        return Matrix{std::forward<Result>(result)};
    }
}
 
}  // namespace internal
 
namespace roots
{

template <typename Head, typename... Fallbacks>
struct WithFallbacks
{
    static_assert(Axes<Head>, "a fallback list starts with axes");
    static_assert(internal::entryCanFail<Head>(), "axes that cannot fail take no fallback");
    static_assert(((Axes<Fallbacks> || Repair<Fallbacks>) && ...), "a fallback is axes or a repair");
    static_assert(
        (!internal::isFallbackList<Fallbacks> && ...),
        "a fallback list is flat: give a fallback's fallbacks to the list itself"
    );
    static_assert(
        internal::fallbacksReached<Fallbacks...>(),
        "a fallback after axes that cannot fail is never reached"
    );

    static constexpr bool canFail =
        internal::entryCanFail<std::tuple_element_t<sizeof...(Fallbacks) - 1, std::tuple<Fallbacks...>>>();

    template <typename Scalar>
    using Report = std::variant<internal::ReportOfT<Head, Scalar>, internal::ReportOfT<Fallbacks, Scalar>...>;

    template <typename Scalar>
    using FailureReport = std::tuple<
        std::optional<internal::ReportOfT<Head, Scalar>>,
        std::optional<internal::ReportOfT<Fallbacks, Scalar>>...>;
 
    Head head;                           
    std::tuple<Fallbacks...> fallbacks;  

    template <typename Derived>
    [[nodiscard]] internal::RootResult<canFail, internal::SquareMatrixOf<Derived>>
    root(Eigen::MatrixBase<Derived> const & cov) const
    {
        return root(cov, internal::DiscardReport{});
    }

    template <typename Derived, typename Sink>
    [[nodiscard]] internal::RootResult<canFail, internal::SquareMatrixOf<Derived>>
    root(Eigen::MatrixBase<Derived> const & cov, Sink && sink) const
    {
        using Matrix = internal::SquareMatrixOf<Derived>;
        using Scalar = typename Derived::Scalar;
        auto tried = FailureReport<Scalar>{};
        auto const recordInto = [&tried]<std::size_t I>(std::integral_constant<std::size_t, I>) {
            return [&tried](auto const & report) {
                std::get<I>(tried) = report;
            };
        };
        std::optional<Matrix> found;
        std::size_t used = 0;
        auto error = GaussianError::NotPositiveDefinite;
        auto const keep = [&](std::size_t index, auto && result) {
            auto const root = internal::asExpected<Matrix>(std::forward<decltype(result)>(result));
            if (root)
            {
                found = *root;
                used = index;
            }
            else
            {
                error = root.error();
            }
        };
        keep(0, internal::rootWithSink(head, cov, recordInto(std::integral_constant<std::size_t, 0>{})));
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (
                [&] {
                    if (found)
                    {
                        return;
                    }
                    auto const & entry = std::get<I>(fallbacks);
                    auto const record = recordInto(std::integral_constant<std::size_t, I + 1>{});
                    using Entry = std::tuple_element_t<I, std::tuple<Fallbacks...>>;
                    if constexpr (Repair<Entry>)
                    {
                        auto const fixed = internal::repairWithSink(entry, cov, record);
                        if (!fixed)
                        {
                            error = fixed.error();
                            return;
                        }
                        keep(I + 1, head.root(*fixed));
                    }
                    else
                    {
                        keep(I + 1, internal::rootWithSink(entry, cov, record));
                    }
                }(),
                ...
            );
        }(std::index_sequence_for<Fallbacks...>{});
        if (found)
        {
            sink(successReport<Scalar>(tried, used));
            return std::move(*found);
        }
        if constexpr (canFail)
        {
            sink(tried);
            return std::unexpected{error};
        }
        else
        {
            EMRYS_FILTER_ASSERT(false && "the last fallback cannot fail");
            return Matrix{};
        }
    }

    template <typename Derived>
    [[nodiscard]] internal::RootResult<canFail, internal::SquareMatrixOf<Derived>>
    rootFromFactor(Eigen::MatrixBase<Derived> const & factor) const
    {
        using Matrix = internal::SquareMatrixOf<Derived>;
        auto first = internal::asExpected<Matrix>(internal::rootFromFactor(head, factor));
        if (first)
        {
            return std::move(*first);
        }
        std::optional<Matrix> found;
        auto error = first.error();
        std::optional<Matrix> covariance;
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (
                [&] {
                    if (found)
                    {
                        return;
                    }
                    auto const & entry = std::get<I>(fallbacks);
                    using Entry = std::tuple_element_t<I, std::tuple<Fallbacks...>>;
                    std::expected<Matrix, GaussianError> root = std::unexpected{error};
                    if constexpr (Repair<Entry>)
                    {
                        if (!covariance)
                        {
                            covariance = Matrix{factor * factor.transpose()};
                        }
                        auto const fixed = entry.repaired(*covariance);
                        if (!fixed)
                        {
                            error = fixed.error();
                            return;
                        }
                        root = internal::asExpected<Matrix>(head.root(*fixed));
                    }
                    else
                    {
                        root = internal::asExpected<Matrix>(internal::rootFromFactor(entry, factor));
                    }
                    if (root)
                    {
                        found = std::move(*root);
                    }
                    else
                    {
                        error = root.error();
                    }
                }(),
                ...
            );
        }(std::index_sequence_for<Fallbacks...>{});
        if (found)
        {
            return std::move(*found);
        }
        if constexpr (canFail)
        {
            return std::unexpected{error};
        }
        else
        {
            EMRYS_FILTER_ASSERT(false && "the last fallback cannot fail");
            return Matrix{};
        }
    }
 
private:

    template <typename Scalar, std::size_t I = 0>
    [[nodiscard]] static Report<Scalar>
    successReport(FailureReport<Scalar> const & tried, std::size_t used)
    {
        if constexpr (I + 1 < std::tuple_size_v<FailureReport<Scalar>>)
        {
            if (used != I)
            {
                return successReport<Scalar, I + 1>(tried, used);
            }
        }
        EMRYS_FILTER_ASSERT(std::get<I>(tried).has_value());
        return Report<Scalar>{std::in_place_index<I>, *std::get<I>(tried)};
    }
};
 
}  // namespace roots
 
namespace internal
{

template <typename Head, typename... Fallbacks>
inline constexpr bool isFallbackList<roots::WithFallbacks<Head, Fallbacks...>>{true};


template <typename Derived>
[[nodiscard]] SquareMatrixOf<Derived>
defaultRoot(Eigen::MatrixBase<Derived> const & cov)
{
    return roots::WithFallbacks<roots::Cholesky, roots::Pivoted>{}.root(cov);
}

template <typename Head, typename... Fallbacks>
[[nodiscard]] constexpr auto
withFallbacks(Head head, Fallbacks... fallbacks)
{
    if constexpr (sizeof...(Fallbacks) == 0)
    {
        return head;
    }
    else
    {
        return roots::WithFallbacks<Head, Fallbacks...>{
            std::move(head),
            std::tuple<Fallbacks...>{std::move(fallbacks)...}
        };
    }
}
 
}  // namespace internal
 
namespace roots
{

template <typename... Fallbacks>
[[nodiscard]] constexpr auto
cholesky(Fallbacks... fallbacks)
{
    return internal::withFallbacks(Cholesky{}, std::move(fallbacks)...);
}

template <typename... Fallbacks>
[[nodiscard]] constexpr auto
pivoted(Fallbacks... fallbacks)
{
    return internal::withFallbacks(Pivoted{}, std::move(fallbacks)...);
}

template <typename... Fallbacks>
[[nodiscard]] constexpr auto
principal(Fallbacks... fallbacks)
{
    return internal::withFallbacks(Principal{}, std::move(fallbacks)...);
}

template <typename... Fallbacks>
[[nodiscard]] constexpr auto
symmetric(Fallbacks... fallbacks)
{
    return internal::withFallbacks(Symmetric{}, std::move(fallbacks)...);
}

template <typename Derived, typename... Fallbacks>
[[nodiscard]] auto
along(Eigen::MatrixBase<Derived> const & direction, Fallbacks... fallbacks)
{
    static_assert(Derived::ColsAtCompileTime == 1, "along() takes a direction as a column vector");
    return internal::withFallbacks(
        Along<typename Derived::PlainObject>{typename Derived::PlainObject{direction}},
        std::move(fallbacks)...
    );
}

template <double Relative = 1.0, double Absolute = 0.0>
[[nodiscard]] constexpr Ridge<Relative, Absolute>
ridge()
{
    return {};
}

template <typename Head, typename First, typename... Rest>
[[nodiscard]] constexpr WithFallbacks<Head, First, Rest...>
withFallbacks(Head axes, First first, Rest... rest)
{
    return {std::move(axes), std::tuple<First, Rest...>{std::move(first), std::move(rest)...}};
}

using DefaultAxes = WithFallbacks<Cholesky, Pivoted>;

template <Axes A, typename Derived>
[[nodiscard]] auto
squareRoot(A const & axes, Eigen::MatrixBase<Derived> const & cov)
{
    EMRYS_FILTER_ASSERT(cov.rows() == cov.cols() && "a covariance must be square");
    return axes.root(cov);
}

template <typename Matrix, typename Report>
struct RootWithReport
{
    Matrix root;    
    Report report;  
};

template <typename Report>
struct RootFailure
{
    GaussianError error;  
    Report report;        
};

template <Axes A, typename Derived>
[[nodiscard]] auto
squareRootWithReport(A const & axes, Eigen::MatrixBase<Derived> const & cov)
{
    using Scalar = typename Derived::Scalar;
    using Matrix = internal::SquareMatrixOf<Derived>;
    using Report = internal::ReportOfT<A, Scalar>;
    using Failure = internal::FailureReportOfT<A, Scalar>;
    std::optional<Report> report;
    std::optional<Failure> failure;
    auto result = internal::rootWithSink(axes, cov, [&](auto const & given) {
        using Given = std::remove_cvref_t<decltype(given)>;
        if constexpr (std::same_as<Given, Report>)
        {
            report = given;
        }
        if constexpr (std::same_as<Given, Failure>)
        {
            failure = given;
        }
    });
    if constexpr (A::canFail)
    {
        using Result = std::expected<RootWithReport<Matrix, Report>, RootFailure<Failure>>;
        if (!result)
        {
            return Result{std::unexpect, RootFailure<Failure>{result.error(), failure.value_or(Failure{})}};
        }
        return Result{RootWithReport<Matrix, Report>{std::move(*result), report.value_or(Report{})}};
    }
    else
    {
        return RootWithReport<Matrix, Report>{std::move(result), report.value_or(Report{})};
    }
}
 
}  // namespace roots
 
}  // namespace Emrys