Product.h

Go to the documentation of this file.

#pragma once

 
#include <algorithm>
#include <array>
#include <concepts>
#include <cstddef>
#include <expected>
#include <optional>
#include <random>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/NumericError.h>
#include <EmrysFilter/Distributions/Moments.h>
#include <EmrysFilter/Distributions/SampleSet.h>
#include <EmrysFilter/Distributions/Sampling.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/Subspace.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{

template <typename T>
concept HasLogDensity = requires(
    T const & density,
    Eigen::Matrix<
        typename T::Scalar,
        T::ValueDimensionAtCompileTime,
        1,
        Eigen::ColMajor,
        T::MaxValueDimensionAtCompileTime,
        1> const & x
) {
    typename T::SpaceType;
    { Emrys::valueDimension(density.space()) } -> std::convertible_to<Eigen::Index>;
    { density.logDensity(x) } -> std::convertible_to<typename T::Scalar>;
};
 
namespace internal
{

template <typename T>
concept CheckedParams = requires(T const & params) {
    { params.toDensity().has_value() } -> std::convertible_to<bool>;
    requires HasLogDensity<typename decltype(params.toDensity())::value_type>;
};

template <typename T>
concept HasMoments = requires(T const & distribution) { distribution.moments(); };
 
}  // namespace internal

template <typename T>
concept ProductFactor = (HasLogDensity<T> && internal::HasMoments<T>)
                     || (internal::CheckedParams<T>
                         && internal::HasMoments<typename decltype(std::declval<T const &>().toDensity())::value_type>);
 
namespace internal
{

template <typename Own, typename Block>
inline constexpr bool factorFits = std::is_same_v<Own, Block> || std::derived_from<Block, Own>
                                || std::derived_from<Own, Block> || (isLinear<Own> && isLinear<Block>);

template <typename Placed>
struct PlacedFactor
{};

template <ProductFactor F, typename Path>
struct PlacedFactor<OnPath<F, Path>>
{
    template <typename Space>
    using TypeIn = F;

    template <typename Space>
    [[nodiscard]] static F
    factorIn(OnPath<F, Path> && placed, Space const &)
    {
        return std::move(placed.value);
    }

    template <typename Space>
    [[nodiscard]] static auto
    indicesIn(OnPath<F, Path> const &, Space const & space)
    {
        return asIndexSet(idx<Path>(space));
    }

    template <typename Space>
    static constexpr bool fitsIn = factorFits<typename F::SpaceType, PathSubspaceType<Space, Path>>;
};

template <ProductFactor F, typename Indices>
struct FactorAt
{
    F factor;         
    Indices indices;  
};

template <ProductFactor F, typename Indices>
struct PlacedFactor<FactorAt<F, Indices>>
{
    template <typename Space>
    using TypeIn = F;

    template <typename Space>
    [[nodiscard]] static F
    factorIn(FactorAt<F, Indices> && placed, Space const &)
    {
        return std::move(placed.factor);
    }

    template <typename Space>
    [[nodiscard]] static auto
    indicesIn(FactorAt<F, Indices> const & placed, Space const &)
    {
        return asIndexSet(placed.indices);
    }

    template <typename Space>
    static constexpr bool fitsIn = true;
};

template <typename Indices, ProductFactor F>
[[nodiscard]] FactorAt<F, Indices>
placedAt(Indices indices, F factor)
{
    return {std::move(factor), std::move(indices)};
}

template <typename F>
[[nodiscard]] decltype(auto)
densityOf(F const & factor)
{
    if constexpr (HasLogDensity<F>)
    {
        return (factor);
    }
    else
    {
        auto const density = factor.toDensity();
        auto view = *density;
        return view;
    }
}

template <typename F>
inline constexpr bool momentsMayFail =
    isExpected<std::remove_cvref_t<decltype(densityOf(std::declval<F const &>()).moments())>>;

template <typename... T>
struct TypeList
{};

template <std::size_t N, typename... T>
struct DropFrontImpl
{
    using type = TypeList<T...>;
};

template <std::size_t N, typename First, typename... Rest>
    requires(N > 0)
struct DropFrontImpl<N, First, Rest...> : DropFrontImpl<N - 1, Rest...>
{};

template <std::size_t N, typename... T>
using DropFront = typename DropFrontImpl<N, T...>::type;

template <std::size_t N, typename Done, typename... T>
struct TakeFrontImpl
{
    using type = Done;
};

template <std::size_t N, typename... Done, typename First, typename... Rest>
    requires(N > 0)
struct TakeFrontImpl<N, TypeList<Done...>, First, Rest...> : TakeFrontImpl<N - 1, TypeList<Done..., First>, Rest...>
{};

template <std::size_t N, typename... T>
using TakeFront = typename TakeFrontImpl<N, TypeList<>, T...>::type;

template <int Size, int MaxSize, typename Indices, typename Target>
[[nodiscard]] IndexVector<Size, MaxSize>
positionsIn(Indices const & indices, Target const & target)
{
    auto result = IndexVector<Size, MaxSize>{};
    result.resize(Eigen::Index(indices.size()));
    for (Eigen::Index i = 0; i < result.size(); ++i)
    {
        auto const found = std::ranges::find(target, indices[i]);
        EMRYS_FILTER_ASSERT(found != std::ranges::end(target) && "a factor's coordinate is missing from the result");
        result(i) = Eigen::Index(found - std::ranges::begin(target));
    }
    return result;
}

template <typename R>
struct ConditionalValueOf
{
    using type = R;
};

template <typename T, typename E>
struct ConditionalValueOf<std::expected<T, E>>
{
    using type = T;
};

template <typename R>
using ConditionalValue = typename ConditionalValueOf<R>::type;

template <typename R>
struct ConditionalErrorOf
{
    using type = void;
};

template <typename T, typename E>
struct ConditionalErrorOf<std::expected<T, E>>
{
    using type = E;
};

template <CheckedParams R>
struct ConditionalErrorOf<R>
{
    using type = typename decltype(std::declval<R const &>().toDensity())::error_type;
};

template <typename R>
using ConditionalError = typename ConditionalErrorOf<R>::type;
 
}  // namespace internal
 
template <Space Space_, RealScalar Scalar_, ProductFactor... Factors_>
class Product;
 
namespace internal
{

template <
    typename Scalar,
    typename Space,
    std::size_t K,
    typename Replacement,
    typename Factors_,
    typename Sequence = std::make_index_sequence<std::tuple_size_v<Factors_>>>
struct ReplacedProductOf;

template <typename Scalar, typename Space, std::size_t K, typename Replacement, typename... Factors_, std::size_t... J>
struct ReplacedProductOf<Scalar, Space, K, Replacement, std::tuple<Factors_...>, std::index_sequence<J...>>
{
    using type = Product<Space, Scalar, std::conditional_t<J == K, Replacement, Factors_>...>;
};

template <typename Scalar, typename Space, std::size_t K, typename Replacement, typename Factors_>
using ReplacedProduct = typename ReplacedProductOf<Scalar, Space, K, Replacement, Factors_>::type;
 
}  // namespace internal

template <ProductFactor F, Space BlockSpace>
class EachRepetition
{
    static_assert(
        requires { typename BlockSpace::BlockType; } && BlockSpace::FixedValueDimension == 0,
        "onEach<Path>() places a factor on a repeated block"
    );

    static constexpr int Step = BlockSpace::ValueDimensionStep;

    static constexpr bool MomentsMayFail = internal::momentsMayFail<F>;
 
public:
 
    using Scalar = typename F::Scalar;                    
    using Index = Eigen::Index;                           
    using SpaceType = BlockSpace;                         
    using MomentsType = Moments<BlockSpace, Scalar>;      
    using VectorType = typename MomentsType::VectorType;  

    static constexpr int ValueDimensionAtCompileTime = MomentsType::ValueDimensionAtCompileTime;

    static constexpr int MaxValueDimensionAtCompileTime = MomentsType::MaxValueDimensionAtCompileTime;

    EachRepetition(F factor, BlockSpace space)
        : m_factor{std::move(factor)}
        , m_space{std::move(space)}
    {
        EMRYS_FILTER_ASSERT(
            Emrys::valueDimension(m_factor.space()) == Step && "the factor must cover one repetition of the block"
        );
        if constexpr (!HasLogDensity<F>)
        {
            [[maybe_unused]] bool const checked = m_factor.toDensity().has_value();
            EMRYS_FILTER_ASSERT(checked && "a factor failed its check; call its toDensity() first to handle this");
        }
    }

    [[nodiscard]] F const &
    factor() const
    {
        return m_factor;
    }

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

    [[nodiscard]] Index
    count() const
    {
        return m_space.count();
    }

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

    template <typename Derived>
    [[nodiscard]] Scalar
    logDensity(Eigen::MatrixBase<Derived> const & x) const
    {
        static_assert(Derived::IsVectorAtCompileTime, "logDensity() expects a vector");
        EMRYS_FILTER_ASSERT(x.size() == valueDimension());
        decltype(auto) density = internal::densityOf(m_factor);
        Scalar result = 0;
        for (Index k = 0; k < count(); ++k)
        {
            result += Scalar(density.logDensity(x.derived().template segment<Step>(k * Step)));
        }
        return result;
    }

    [[nodiscard]] auto
    moments() const -> std::conditional_t<MomentsMayFail, std::expected<MomentsType, MomentsError>, MomentsType>
    {
        decltype(auto) density = internal::densityOf(m_factor);
        auto const & one = density.moments();
        auto const place = [&](auto const & moments) {
            auto result = MomentsType{valueDimension(), valueDimension()};
            result.mean().setZero();
            result.covarianceLower().setZero();
            for (Index k = 0; k < count(); ++k)
            {
                result.mean().template segment<Step>(k * Step) = moments.mean();
                result.covarianceLower().template block<Step, Step>(k * Step, k * Step) =
                    moments.covarianceLower().template selfadjointView<Eigen::Lower>();
            }
            return result;
        };
        if constexpr (MomentsMayFail)
        {
            if (!one)
            {
                return std::unexpected{one.error()};
            }
            return place(*one);
        }
        else
        {
            return place(one);
        }
    }
 
private:
 
    F m_factor;                                         
    EMRYS_FILTER_NO_UNIQUE_ADDRESS BlockSpace m_space;  
};

template <typename F, typename BlockSpace, std::uniform_random_bit_generator G>
[[nodiscard]] SampleSet<BlockSpace, typename F::Scalar>
drawSamples(EachRepetition<F, BlockSpace> const & each, Eigen::Index n, G & rng)
{
    EMRYS_FILTER_ASSERT(n >= 0);
    constexpr int Step = BlockSpace::ValueDimensionStep;
    auto samples = SampleSet<BlockSpace, typename F::Scalar>{each.valueDimension(), n};
    decltype(auto) density = internal::densityOf(each.factor());
    for (Eigen::Index k = 0; k < each.count(); ++k)
    {
        samples.template middleRows<Step>(k * Step) = drawSamples(density, n, rng);
    }
    return samples;
}

template <typename F, typename BlockSpace, std::uniform_random_bit_generator G>
[[nodiscard]] typename Moments<BlockSpace, typename F::Scalar>::VectorType
drawSample(EachRepetition<F, BlockSpace> const & each, G & rng)
{
    constexpr int Step = BlockSpace::ValueDimensionStep;
    auto sample = typename Moments<BlockSpace, typename F::Scalar>::VectorType{};
    sample.resize(each.valueDimension());
    decltype(auto) density = internal::densityOf(each.factor());
    for (Eigen::Index k = 0; k < each.count(); ++k)
    {
        sample.template segment<Step>(k * Step) = drawSample(density, rng);
    }
    return sample;
}

template <typename F, typename BlockSpace, std::uniform_random_bit_generator G>
[[nodiscard]] SamplesWithLogDensity<BlockSpace, typename F::Scalar>
drawSamplesWithLogDensity(EachRepetition<F, BlockSpace> const & each, Eigen::Index n, G & rng)
{
    EMRYS_FILTER_ASSERT(n >= 0);
    using Scalar = typename F::Scalar;
    constexpr int Step = BlockSpace::ValueDimensionStep;
    auto result = SamplesWithLogDensity<BlockSpace, Scalar>{
        SampleSet<BlockSpace, Scalar>{each.valueDimension(), n},
        Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n)
    };
    decltype(auto) density = internal::densityOf(each.factor());
    for (Eigen::Index k = 0; k < each.count(); ++k)
    {
        auto draws = drawSamplesWithLogDensity(density, n, rng);
        result.samples.template middleRows<Step>(k * Step) = draws.samples;
        result.logDensity += draws.logDensity;
    }
    return result;
}

template <typename F, typename BlockSpace, std::uniform_random_bit_generator G>
[[nodiscard]] SampleWithLogDensity<BlockSpace, typename F::Scalar>
drawSampleWithLogDensity(EachRepetition<F, BlockSpace> const & each, G & rng)
{
    constexpr int Step = BlockSpace::ValueDimensionStep;
    auto result = SampleWithLogDensity<BlockSpace, typename F::Scalar>{};
    result.sample.resize(each.valueDimension());
    decltype(auto) density = internal::densityOf(each.factor());
    for (Eigen::Index k = 0; k < each.count(); ++k)
    {
        auto const draw = drawSampleWithLogDensity(density, rng);
        result.sample.template segment<Step>(k * Step) = draw.sample;
        result.logDensity += draw.logDensity;
    }
    return result;
}
 
namespace internal
{

template <typename Space, typename Path>
struct EachBlockOf
{
    using type = PathSubspaceType<Space, Path>;

    [[nodiscard]] static type
    in(Space const & space)
    {
        return pathSubspace<Path>(space);
    }

    [[nodiscard]] static auto
    indicesIn(Space const & space)
    {
        return asIndexSet(idx<Path>(space));
    }
};

template <typename Space>
struct EachBlockOf<Space, WholeSpace>
{
    using type = Space;

    [[nodiscard]] static Space const &
    in(Space const & space)
    {
        return space;
    }

    [[nodiscard]] static auto
    indicesIn(Space const & space)
    {
        return asIndexSet(Eigen::seqN(Eigen::Index{0}, valueDimension(space)));
    }
};

template <ProductFactor F, typename Path>
struct PlacedFactor<OnEach<F, Path>>
{
    template <typename Space>
    using TypeIn = EachRepetition<F, typename EachBlockOf<Space, Path>::type>;

    template <typename Space>
    [[nodiscard]] static TypeIn<Space>
    factorIn(OnEach<F, Path> && placed, Space const & space)
    {
        return {std::move(placed.value), EachBlockOf<Space, Path>::in(space)};
    }

    template <typename Space>
    [[nodiscard]] static auto
    indicesIn(OnEach<F, Path> const &, Space const & space)
    {
        return EachBlockOf<Space, Path>::indicesIn(space);
    }

    template <typename Space>
    static constexpr bool fitsIn = [] {
        using Block = typename EachBlockOf<Space, Path>::type;
        if constexpr (requires { typename Block::BlockType; })
        {
            return factorFits<typename F::SpaceType, typename Block::BlockType>;
        }
        else
        {
            return false;
        }
    }();
};

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

template <ProductFactor F, typename Path>
inline constexpr bool isPlacedFactor<OnPath<F, Path>>{true};
 
template <ProductFactor F, typename Path>
inline constexpr bool isPlacedFactor<OnEach<F, Path>>{true};
 
}  // namespace internal

template <typename... Placed>
class Factors
{
    static_assert(sizeof...(Placed) > 0, "a product needs at least one factor");
    static_assert(
        (internal::isPlacedFactor<Placed> && ...),
        "each factor is placed by tag path, with on<Path>() or onEach<Path>()"
    );
 
public:

    explicit Factors(Placed... placed)
        : m_placed{std::move(placed)...}
    {}

    [[nodiscard]] std::tuple<Placed...> const &
    placed() const
    {
        return m_placed;
    }
 
private:
 
    std::tuple<Placed...> m_placed;  
};
 
namespace internal
{

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

template <typename... Placed>
inline constexpr bool isFactors<Factors<Placed...>>{true};
 
}  // namespace internal

template <Space Space_, RealScalar Scalar_, ProductFactor... Factors_>
class Product
{
    static_assert(sizeof...(Factors_) > 0, "a product needs at least one factor");
    static_assert(
        (std::is_same_v<typename Factors_::Scalar, Scalar_> && ...),
        "every factor must have the same scalar"
    );

    static constexpr bool MomentsMayFail = (internal::momentsMayFail<Factors_> || ...);
 
public:
 
    using Scalar = Scalar_;                               
    using Index = Eigen::Index;                           
    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;
 
    static constexpr std::size_t FactorCount = sizeof...(Factors_);  
 
    static_assert(
        internal::hasEqualSizes<Space_>,
        "a product places its factors by coordinate, so its space's values and differences must have the "
        "same coordinates"
    );

    template <std::size_t K>
    using FactorType = std::tuple_element_t<K, std::tuple<Factors_...>>;

    template <typename... Placed>
        requires(sizeof...(Placed) == sizeof...(Factors_)
                 && (std::is_same_v<typename internal::PlacedFactor<Placed>::template TypeIn<Space_>, Factors_> && ...))
    explicit Product(Space_ space, Placed... placed)
        : m_space{std::move(space)}
        , m_indices{resolve<Factors_>(placed)...}
        , m_factors{internal::PlacedFactor<Placed>::factorIn(std::move(placed), m_space)...}
    {
        static_assert(
            (internal::PlacedFactor<Placed>::template fitsIn<Space_> && ...),
            "a factor's space must be the space of the block it covers"
        );
        EMRYS_FILTER_ASSERT(sizesMatch() && "a factor's dimension must match the number of coordinates it covers");
        EMRYS_FILTER_ASSERT(partitionsSpace() && "the factors' indices must partition the space");
        std::apply(
            [](auto const &... factors) {
                (checkParams(factors), ...);
            },
            m_factors
        );
    }

    template <typename... Placed>
    explicit Product(Space_ space, Factors<Placed...> const & factors)
        : Product{std::move(space), factors, std::index_sequence_for<Placed...>{}}
    {}

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

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

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

    template <std::size_t K>
    [[nodiscard]] FactorType<K> const &
    factor() const
    {
        return std::get<K>(m_factors);
    }

    template <std::size_t K>
    [[nodiscard]] auto const &
    indices() const
    {
        return std::get<K>(m_indices);
    }

    template <typename Fn>
    void
    forEachFactor(Fn && fn) const
    {
        [&]<std::size_t... K>(std::index_sequence<K...>) {
            (fn(internal::densityOf(std::get<K>(m_factors)), std::get<K>(m_indices)), ...);
        }(std::index_sequence_for<Factors_...>{});
    }

    template <typename Derived>
    [[nodiscard]] Scalar
    logDensity(Eigen::MatrixBase<Derived> const & x) const
    {
        static_assert(Derived::IsVectorAtCompileTime, "logDensity() expects a vector");
        EMRYS_FILTER_ASSERT(x.size() == valueDimension());
        Scalar result = 0;
        forEachFactor([&](auto const & density, auto const & indices) {
            result += Scalar(density.logDensity(x.derived()(indices)));
        });
        return result;
    }

    template <typename Derived>
    [[nodiscard]] Scalar
    density(Eigen::MatrixBase<Derived> const & x) const
    {
        return Eigen::numext::exp(logDensity(x));
    }

    [[nodiscard]] auto
    moments() const -> std::conditional_t<MomentsMayFail, std::expected<MomentsType, MomentsError>, MomentsType>
    {
        auto result = MomentsType{valueDimension(), differenceDimension()};
        result.mean().setZero();
        result.covarianceLower().setZero();
        auto error = std::optional<MomentsError>{};
        auto const place = [&](auto const & moments, auto const & indices) {
            using FactorMatrix = std::remove_cvref_t<decltype(moments.covarianceLower())>;
            result.mean()(indices) = moments.mean();
            result.covarianceLower()(indices, indices) =
                FactorMatrix{moments.covarianceLower().template selfadjointView<Eigen::Lower>()};
        };
        forEachFactor([&](auto const & density, auto const & indices) {
            if (error)
            {
                return;
            }
            auto const & factorMoments = density.moments();
            if constexpr (internal::isExpected<std::remove_cvref_t<decltype(factorMoments)>>)
            {
                if (!factorMoments)
                {
                    error = factorMoments.error();
                }
                else
                {
                    place(*factorMoments, indices);
                }
            }
            else
            {
                place(factorMoments, indices);
            }
        });
        if constexpr (MomentsMayFail)
        {
            if (error)
            {
                return std::unexpected{*error};
            }
        }
        return result;
    }

    template <typename Path>
    [[nodiscard]] decltype(auto)
    marginal() const &
    {
        return internal::withPathTags<Path>([&]<typename... Tags>() -> decltype(auto) {
            return marginalAlong<Tags...>();
        });
    }

    template <typename Path>
    void
    marginal() const && = delete;

    template <typename GivenPath, typename Derived>
    [[nodiscard]] auto
    conditional(Eigen::MatrixBase<Derived> const & values) const
    {
        return internal::withPathTags<GivenPath>([&]<typename... Tags>() {
            return conditionalAlong<Tags...>(values);
        });
    }
 
private:

    template <typename... Placed, std::size_t... K>
    Product(Space_ space, Factors<Placed...> const & factors, std::index_sequence<K...>)
        : Product{std::move(space), std::get<K>(factors.placed())...}
    {}

    template <typename... Path>
    [[nodiscard]] decltype(auto)
    marginalAlong() const
    {
        using Match = MatchType<Path...>;
        if constexpr (Match::value.factor != None)
        {
            assertOnOwnBlock<Match::value.factor>(OwnBlockPath<Match::value, Path...>{});
            if constexpr (Match::value.whole)
            {
                return std::get<Match::value.factor>(m_factors);
            }
            else
            {
                return factorMarginal(
                    std::get<Match::value.factor>(m_factors),
                    internal::DropFront<Match::value.start, Path...>{}
                );
            }
        }
        else
        {
            static_assert(
                std::ranges::count(withinOf<Path...>(), true) > 0,
                "the path neither reaches into a factor nor contains one"
            );
            auto const target = internal::resolveIndices<ValueDimensionAtCompileTime, MaxValueDimensionAtCompileTime>(
                valueDimension(),
                internal::asIndexSet(idx<TagPath<Path...>>(m_space))
            );
            EMRYS_FILTER_ASSERT(
                coversExactly(withinOf<Path...>(), target) && "the path's block must be made of whole factors"
            );
            return keepOnly<withinOf<Path...>()>(target, internal::subspace<Space_, Path...>(m_space));
        }
    }

    template <typename... GivenPath, typename Derived>
    [[nodiscard]] auto
    conditionalAlong(Eigen::MatrixBase<Derived> const & values) const
    {
        return internal::tagConditional<TagPath<GivenPath...>>(
            m_space,
            [&](auto const & targetCoordinates, auto const & givenCoordinates, auto targetSpace, auto const &) {
                // The values and the differences have the same coordinates.
                auto const & target = targetCoordinates.values;
                auto const & given = givenCoordinates.values;
                EMRYS_FILTER_ASSERT(values.size() == given.size());
                using Match = MatchType<GivenPath...>;
                if constexpr (Match::value.factor != None && !Match::value.whole)
                {
                    assertOnOwnBlock<Match::value.factor>(OwnBlockPath<Match::value, GivenPath...>{});
                    return replaceFactor<Match::value.factor>(
                        target,
                        given,
                        std::move(targetSpace),
                        internal::DropFront<Match::value.start, GivenPath...>{},
                        values
                    );
                }
                else
                {
                    static_assert(
                        std::ranges::count(droppedOf<GivenPath...>(), true) > 0,
                        "the path neither reaches into a factor nor contains one"
                    );
                    static_assert(
                        std::ranges::count(droppedOf<GivenPath...>(), true) < std::ptrdiff_t(FactorCount),
                        "conditioning on every factor leaves nothing"
                    );
                    EMRYS_FILTER_ASSERT(
                        coversExactly(droppedOf<GivenPath...>(), given)
                        && "the given block must be made of whole factors"
                    );
                    return keepOnly<complement(droppedOf<GivenPath...>())>(target, std::move(targetSpace));
                }
            }
        );
    }
 
    static constexpr std::size_t None = FactorCount;  

    struct PathMatch
    {
        std::size_t factor = None;  
        std::size_t start = 0;      
        bool whole = false;         
        bool named = false;         
    };

    template <PathMatch Match, typename... Path>
    using OwnBlockPath = std::conditional_t<
        Match.named,
        internal::TakeFront<Match.start + 1, Path...>,
        internal::TypeList<typename FactorType<Match.factor>::SpaceType>>;

    template <typename T>
    static consteval std::size_t
    factorNamed()
    {
        constexpr auto same = std::array<bool, FactorCount>{std::is_same_v<T, typename Factors_::SpaceType>...};
        static_assert(std::ranges::count(same, true) <= 1, "two factors share this space type, so it names neither");
        auto const found = std::ranges::find(same, true);
        return found == same.end() ? None : std::size_t(found - same.begin());
    }

    template <typename T>
    static consteval std::size_t
    factorContaining()
    {
        constexpr auto inside =
            std::array<bool, FactorCount>{(internal::countTag<T, typename Factors_::SpaceType>() > 0)...};
        static_assert(std::ranges::count(inside, true) <= 1, "this tag lies in more than one factor's space");
        auto const found = std::ranges::find(inside, true);
        return found == inside.end() ? None : std::size_t(found - inside.begin());
    }

    template <typename T>
    static consteval std::array<bool, FactorCount>
    factorsWithin()
    {
        return {(internal::countTag<typename Factors_::SpaceType, T>() > 0)...};
    }

    template <typename... Path>
    static consteval PathMatch
    matchPath()
    {
        constexpr auto named = std::array<std::size_t, sizeof...(Path)>{factorNamed<Path>()...};
        for (std::size_t j = 0; j < named.size(); ++j)
        {
            if (named[j] != None)
            {
                return {named[j], j, j + 1 == named.size(), true};
            }
        }
        return {factorContaining<std::tuple_element_t<0, std::tuple<Path...>>>(), 0, false};
    }

    template <PathMatch Match, typename Last>
    static consteval std::array<bool, FactorCount>
    droppedBy()
    {
        if constexpr (Match.factor != None)
        {
            auto result = std::array<bool, FactorCount>{};
            result[Match.factor] = true;
            return result;
        }
        else
        {
            return factorsWithin<Last>();
        }
    }
 
    // The path analyses are types and consteval functions, not local constexpr variables: GCC can
    // conflate a local used as a template argument across instantiations of a generic lambda.

    template <typename... Path>
    using MatchType = std::integral_constant<PathMatch, matchPath<Path...>()>;

    template <typename... Path>
    using LastTag = std::tuple_element_t<sizeof...(Path) - 1, std::tuple<Path...>>;

    template <typename... Path>
    static consteval std::array<bool, FactorCount>
    withinOf()
    {
        return factorsWithin<LastTag<Path...>>();
    }

    template <typename... Path>
    static consteval std::array<bool, FactorCount>
    droppedOf()
    {
        return droppedBy<MatchType<Path...>::value, LastTag<Path...>>();
    }

    static consteval std::array<bool, FactorCount>
    complement(std::array<bool, FactorCount> mask)
    {
        for (bool & b : mask)
        {
            b = !b;
        }
        return mask;
    }

    template <std::array<bool, FactorCount> Mask>
    static consteval auto
    maskedFactors()
    {
        auto result = std::array<std::size_t, std::size_t(std::ranges::count(Mask, true))>{};
        std::size_t out = 0;
        for (std::size_t k = 0; k < FactorCount; ++k)
        {
            if (Mask[k])
            {
                result[out++] = k;
            }
        }
        return result;
    }

    template <typename Factor>
    static void
    checkParams(Factor const & factor)
    {
        if constexpr (!HasLogDensity<Factor>)
        {
            [[maybe_unused]] bool const checked = factor.toDensity().has_value();
            EMRYS_FILTER_ASSERT(checked && "a factor failed its check; call its toDensity() first to handle this");
        }
    }

    template <typename Factor, typename Placed>
    [[nodiscard]] internal::IndexVector<Factor::ValueDimensionAtCompileTime, Factor::MaxValueDimensionAtCompileTime>
    resolve(Placed const & placed) const
    {
        return internal::IndexVector<Factor::ValueDimensionAtCompileTime, Factor::MaxValueDimensionAtCompileTime>{
            internal::resolveIndices<ValueDimensionAtCompileTime, MaxValueDimensionAtCompileTime>(
                valueDimension(),
                internal::PlacedFactor<Placed>::indicesIn(placed, m_space)
            )
        };
    }

    template <typename Fn>
    void
    forEachStored(Fn && fn) const
    {
        [&]<std::size_t... K>(std::index_sequence<K...>) {
            (fn(std::get<K>(m_factors), std::get<K>(m_indices)), ...);
        }(std::index_sequence_for<Factors_...>{});
    }

    [[nodiscard]] bool
    sizesMatch() const
    {
        bool result = true;
        forEachStored([&](auto const & factor, auto const & indices) {
            result = result && Emrys::valueDimension(factor.space()) == indices.size();
        });
        return result;
    }

    [[nodiscard]] bool
    partitionsSpace() const
    {
        Index total = 0;
        forEachStored([&](auto const &, auto const & indices) {
            total += indices.size();
        });
        bool const distinct = std::apply(
            [](auto const &... indices) {
                return internal::distinctIndices(indices...);
            },
            m_indices
        );
        return total == valueDimension() && distinct;
    }

    template <std::size_t K, typename... Path>
    void
    assertOnOwnBlock(internal::TypeList<Path...>) const
    {
        [[maybe_unused]] auto const block =
            internal::resolveIndices<ValueDimensionAtCompileTime, MaxValueDimensionAtCompileTime>(
                valueDimension(),
                internal::asIndexSet(idx<TagPath<Path...>>(m_space))
            );
        EMRYS_FILTER_ASSERT(
            std::equal(block.begin(), block.end(), std::get<K>(m_indices).begin(), std::get<K>(m_indices).end())
            && "a factor reached by a tag path must be placed on that tag's block"
        );
    }

    template <typename Block>
    [[nodiscard]] bool
    coversExactly(std::array<bool, FactorCount> const & mask, Block const & block) const
    {
        std::size_t k = 0;
        std::size_t covered = 0;
        bool inBlock = true;
        forEachStored([&](auto const &, auto const & indices) {
            if (mask[k++])
            {
                covered += std::size_t(indices.size());
                for (Index const i : indices)
                {
                    inBlock = inBlock && std::ranges::find(block, i) != std::ranges::end(block);
                }
            }
        });
        return inBlock && covered == std::size_t(std::ranges::size(block));
    }

    template <typename Factor, typename... Path>
    [[nodiscard]] static auto
    factorMarginal(Factor const & factor, internal::TypeList<Path...>)
    {
        static_assert(
            requires { factor.template marginal<TagPath<Path...>>(); },
            "this factor has no marginal<Path>()"
        );
        return factor.template marginal<TagPath<Path...>>();
    }

    template <std::array<bool, FactorCount> Keep, typename Target, typename SubSpace>
    [[nodiscard]] decltype(auto)
    keepOnly(Target const & target, SubSpace space) const
    {
        constexpr auto kept = maskedFactors<Keep>();
        if constexpr (kept.size() == 1)
        {
            return std::get<kept[0]>(m_factors);
        }
        else
        {
            return [&]<std::size_t... J>(std::index_sequence<J...>) {
                return Product<SubSpace, Scalar, FactorType<kept[J]>...>{
                    std::move(space),
                    internal::placedAt(positionsOf<kept[J]>(target), std::get<kept[J]>(m_factors))...
                };
            }(std::make_index_sequence<kept.size()>{});
        }
    }

    template <std::size_t K, typename Target>
    [[nodiscard]] auto
    positionsOf(Target const & target) const
    {
        using Factor = FactorType<K>;
        return internal::positionsIn<Factor::ValueDimensionAtCompileTime, Factor::MaxValueDimensionAtCompileTime>(
            std::get<K>(m_indices),
            target
        );
    }

    template <std::size_t K, typename Target, typename Given, typename SubSpace, typename... Path, typename Derived>
    [[nodiscard]] auto
    replaceFactor(
        Target const & target,
        Given const & given,
        SubSpace space,
        internal::TypeList<Path...>,
        Eigen::MatrixBase<Derived> const & values
    ) const
    {
        auto const & factor = std::get<K>(m_factors);
        static_assert(
            requires { factor.template conditional<TagPath<Path...>>(values); },
            "this factor has no conditional<Path>()"
        );
        if constexpr (FactorCount == 1)
        {
            return factor.template conditional<TagPath<Path...>>(values);
        }
        else
        {
            using Raw = decltype(factor.template conditional<TagPath<Path...>>(values));
            using Replacement = internal::ConditionalValue<Raw>;
            using Error = internal::ConditionalError<Raw>;
            static_assert(ProductFactor<Replacement>, "this factor's conditional cannot be a factor of a product");
            using Result = internal::ReplacedProduct<Scalar, SubSpace, K, Replacement, std::tuple<Factors_...>>;
 
            auto const build = [&](Replacement replacement) {
                // The replacement covers factor K's coordinates that were not given, in order.
                auto const & own = std::get<K>(m_indices);
                auto const isRemaining = [&](Index i) {
                    return std::ranges::find(given, i) == std::ranges::end(given);
                };
                auto remaining = internal::IndexVector<
                    Replacement::ValueDimensionAtCompileTime,
                    Replacement::MaxValueDimensionAtCompileTime>{};
                remaining.resize(Index(std::ranges::count_if(own, isRemaining)));
                Index r = 0;
                for (Index const i : own)
                {
                    if (isRemaining(i))
                    {
                        remaining(r++) = i;
                    }
                }
                auto positions = internal::positionsIn<
                    Replacement::ValueDimensionAtCompileTime,
                    Replacement::MaxValueDimensionAtCompileTime>(remaining, target);
                return [&]<std::size_t... J>(std::index_sequence<J...>) {
                    return Result{std::move(space), placeFactor<J, K>(target, positions, replacement)...};
                }(std::index_sequence_for<Factors_...>{});
            };
 
            auto raw = factor.template conditional<TagPath<Path...>>(values);
            if constexpr (std::is_void_v<Error>)
            {
                return build(std::move(raw));
            }
            else
            {
                using Return = std::expected<Result, Error>;
                if constexpr (internal::isExpected<Raw>)
                {
                    if (!raw)
                    {
                        return Return{std::unexpect, raw.error()};
                    }
                    return Return{build(std::move(*raw))};
                }
                else
                {
                    if (auto checked = raw.toDensity(); !checked)
                    {
                        return Return{std::unexpect, checked.error()};
                    }
                    return Return{build(std::move(raw))};
                }
            }
        }
    }

    template <std::size_t J, std::size_t K, typename Target, typename Positions, typename Replacement>
    [[nodiscard]] auto
    placeFactor(Target const & target, Positions const & positions, Replacement & replacement) const
    {
        if constexpr (J == K)
        {
            return internal::placedAt(positions, std::move(replacement));
        }
        else
        {
            return internal::placedAt(positionsOf<J>(target), std::get<J>(m_factors));
        }
    }
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Space_ m_space;  

    std::tuple<
        internal::IndexVector<Factors_::ValueDimensionAtCompileTime, Factors_::MaxValueDimensionAtCompileTime>...>
        m_indices;
 
    std::tuple<Factors_...> m_factors;  
};

template <typename Space, typename First, typename... Rest>
Product(Space, First, Rest...) -> Product<
    Space,
    typename internal::PlacedFactor<First>::template TypeIn<Space>::Scalar,
    typename internal::PlacedFactor<First>::template TypeIn<Space>,
    typename internal::PlacedFactor<Rest>::template TypeIn<Space>...>;
 
template <typename Space, typename First, typename... Rest>
Product(Space, Factors<First, Rest...>) -> Product<
    Space,
    typename internal::PlacedFactor<First>::template TypeIn<Space>::Scalar,
    typename internal::PlacedFactor<First>::template TypeIn<Space>,
    typename internal::PlacedFactor<Rest>::template TypeIn<Space>...>;


template <typename Scalar, typename Space, typename... Factors_, std::uniform_random_bit_generator G>
[[nodiscard]] SampleSet<Space, Scalar>
drawSamples(Product<Space, Scalar, Factors_...> const & product, Eigen::Index n, G & rng)
{
    EMRYS_FILTER_ASSERT(n >= 0);
    auto samples = SampleSet<Space, Scalar>{product.valueDimension(), n};
    product.forEachFactor([&](auto const & density, auto const & indices) {
        samples(indices, Eigen::placeholders::all) = drawSamples(density, n, rng);
    });
    return samples;
}

template <typename Scalar, typename Space, typename... Factors_, std::uniform_random_bit_generator G>
[[nodiscard]] typename Moments<Space, Scalar>::VectorType
drawSample(Product<Space, Scalar, Factors_...> const & product, G & rng)
{
    auto sample = typename Moments<Space, Scalar>::VectorType{};
    sample.resize(product.valueDimension());
    product.forEachFactor([&](auto const & density, auto const & indices) {
        sample(indices) = drawSample(density, rng);
    });
    return sample;
}

template <typename Scalar, typename Space, typename... Factors_, std::uniform_random_bit_generator G>
[[nodiscard]] SamplesWithLogDensity<Space, Scalar>
drawSamplesWithLogDensity(Product<Space, Scalar, Factors_...> const & product, Eigen::Index n, G & rng)
{
    EMRYS_FILTER_ASSERT(n >= 0);
    auto result = SamplesWithLogDensity<Space, Scalar>{
        SampleSet<Space, Scalar>{product.valueDimension(), n},
        Eigen::Array<Scalar, Eigen::Dynamic, 1>::Zero(n)
    };
    product.forEachFactor([&](auto const & density, auto const & indices) {
        auto draws = drawSamplesWithLogDensity(density, n, rng);
        result.samples(indices, Eigen::placeholders::all) = draws.samples;
        result.logDensity += draws.logDensity;
    });
    return result;
}

template <typename Scalar, typename Space, typename... Factors_, std::uniform_random_bit_generator G>
[[nodiscard]] SampleWithLogDensity<Space, Scalar>
drawSampleWithLogDensity(Product<Space, Scalar, Factors_...> const & product, G & rng)
{
    auto result = SampleWithLogDensity<Space, Scalar>{};
    result.sample.resize(product.valueDimension());
    product.forEachFactor([&](auto const & density, auto const & indices) {
        auto const draw = drawSampleWithLogDensity(density, rng);
        result.sample(indices) = draw.sample;
        result.logDensity += draw.logDensity;
    });
    return result;
}
 
}  // namespace Emrys