Discrete.h

Go to the documentation of this file.

#pragma once

 
#include <array>
#include <concepts>
#include <cstddef>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Space.h>
#include <EmrysFilter/Spaces/TagPath.h>
 
namespace Emrys
{


template <typename S>
concept DiscreteSpace =
    requires {
        typename std::integral_constant<int, S::Cardinality>;
        typename std::integral_constant<int, S::MaxCardinality>;
    } && (S::Cardinality == Eigen::Dynamic || (S::Cardinality >= 1 && S::MaxCardinality == S::Cardinality))
    && std::regular<S> && requires(S const & space) {
           { space.cardinality() } -> std::convertible_to<Eigen::Index>;
       };
 
namespace internal
{

template <int Cardinality, int MaxCardinality>
class FiniteExtent
{
public:

    [[nodiscard]] static constexpr Eigen::Index
    cardinality()
    {
        return Cardinality;
    }

    friend bool
    operator==(FiniteExtent const &, FiniteExtent const &) = default;
};

template <int MaxCardinality>
class FiniteExtent<Eigen::Dynamic, MaxCardinality>
{
public:

    FiniteExtent() = default;

    constexpr explicit FiniteExtent(Eigen::Index n)
        : m_cardinality{n}
    {
        EMRYS_FILTER_ASSERT(
            n >= 1 && (MaxCardinality == Eigen::Dynamic || n <= MaxCardinality)
            && "the cardinality does not fit the space"
        );
    }

    [[nodiscard]] constexpr Eigen::Index
    cardinality() const
    {
        return m_cardinality;
    }

    friend bool
    operator==(FiniteExtent const &, FiniteExtent const &) = default;
 
private:
 
    Eigen::Index m_cardinality = 1;  
};

template <int... Cardinalities>
consteval int
productCardinality()
{
    int result = 1;
    ((result = result == Eigen::Dynamic || Cardinalities == Eigen::Dynamic ? Eigen::Dynamic : result * Cardinalities),
     ...);
    return result;
}
 
}  // namespace internal

template <int Cardinality_, int MaxCardinality_ = Cardinality_>
struct Finite : internal::FiniteExtent<Cardinality_, MaxCardinality_>
{
    static_assert(Cardinality_ == Eigen::Dynamic || Cardinality_ >= 1, "a discrete space has at least one value");
 
    static_assert(
        Cardinality_ == Eigen::Dynamic || MaxCardinality_ == Cardinality_,
        "a fixed cardinality is its own maximum"
    );
 
    static_assert(MaxCardinality_ == Eigen::Dynamic || MaxCardinality_ >= 1);
    static constexpr int Cardinality = Cardinality_;        
    static constexpr int MaxCardinality = MaxCardinality_;  
 
    Finite() = default;

    constexpr explicit Finite(Eigen::Index n)
        requires(Cardinality_ == Eigen::Dynamic)
        : internal::FiniteExtent<Cardinality_, MaxCardinality_>(n)
    {}

    friend bool
    operator==(Finite const &, Finite const &) = default;
};

template <typename... Values>
struct Named
{
    static_assert(sizeof...(Values) >= 1, "a discrete space has at least one value");
    static constexpr int Cardinality = int(sizeof...(Values));  
    static constexpr int MaxCardinality = Cardinality;          

    template <typename Value>
    static constexpr bool contains = (std::is_same_v<Value, Values> || ...);
 
private:

    template <typename Value>
    static constexpr std::size_t occurrences = (std::size_t{0} + ... + std::size_t{std::is_same_v<Value, Values>});
 
    static_assert(((occurrences<Values> == 1) && ...), "the values of a Named space must be distinct");
 
public:

    [[nodiscard]] static constexpr Eigen::Index
    cardinality()
    {
        return Cardinality;
    }

    template <typename Value>
    [[nodiscard]] static constexpr Eigen::Index
    index()
    {
        static_assert(contains<Value>, "the space has no value with this name");
        Eigen::Index result = 0;
        bool found = false;
        ((found = found || std::is_same_v<Value, Values>, result += found ? 0 : 1), ...);
        return result;
    }

    friend bool
    operator==(Named const &, Named const &) = default;
};

template <DiscreteSpace... Axes>
class DiscreteProduct : private internal::PartStorage<Axes...>
{
    static_assert(sizeof...(Axes) > 0, "a product needs at least one axis");

    using Storage = internal::PartStorage<Axes...>;
 
public:
 
    static constexpr std::size_t AxisCount = sizeof...(Axes);  

    static constexpr int Cardinality = internal::productCardinality<Axes::Cardinality...>();

    static constexpr int MaxCardinality = internal::productCardinality<Axes::MaxCardinality...>();
 
    DiscreteProduct() = default;

    constexpr explicit DiscreteProduct(Axes... axes)
        : Storage{std::move(axes)...}
    {}

    template <std::size_t I>
    [[nodiscard]] constexpr auto const &
    axis() const
    {
        return Storage::template part<I>();
    }

    [[nodiscard]] constexpr Eigen::Index
    cardinality() const
    {
        return [&]<std::size_t... I>(std::index_sequence<I...>) {
            return (Eigen::Index{1} * ... * axis<I>().cardinality());
        }(std::make_index_sequence<AxisCount>{});
    }

    template <std::size_t I>
    [[nodiscard]] constexpr Eigen::Index
    stride() const
    {
        return [&]<std::size_t... J>(std::index_sequence<J...>) {
            return (Eigen::Index{1} * ... * axis<I + 1 + J>().cardinality());
        }(std::make_index_sequence<AxisCount - 1 - I>{});
    }

    [[nodiscard]] constexpr Eigen::Index
    value(std::array<Eigen::Index, AxisCount> const & coordinates) const
    {
        auto const cardinalities = [&]<std::size_t... I>(std::index_sequence<I...>) {
            return std::array<Eigen::Index, AxisCount>{axis<I>().cardinality()...};
        }(std::make_index_sequence<AxisCount>{});
        Eigen::Index result = 0;
        for (std::size_t i = 0; i < AxisCount; ++i)
        {
            EMRYS_FILTER_ASSERT(
                coordinates[i] >= 0 && coordinates[i] < cardinalities[i] && "a coordinate is out of range"
            );
            result = result * cardinalities[i] + coordinates[i];
        }
        return result;
    }

    [[nodiscard]] constexpr std::array<Eigen::Index, AxisCount>
    coordinates(Eigen::Index value) const
    {
        EMRYS_FILTER_ASSERT(value >= 0 && value < cardinality() && "the value is out of range");
        auto result = std::array<Eigen::Index, AxisCount>{};
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((result[I] = (value / stride<I>()) % axis<I>().cardinality()), ...);
        }(std::make_index_sequence<AxisCount>{});
        return result;
    }

    friend bool
    operator==(DiscreteProduct const &, DiscreteProduct const &) = default;
};
 
namespace internal
{

template <typename... Axes>
std::tuple<Axes...> *
discreteProductAxes(DiscreteProduct<Axes...> const *);

template <typename T>
concept IsDiscreteProduct = requires(T const * t) { discreteProductAxes(t); };

template <typename T>
using DiscreteProductAxes = std::remove_pointer_t<decltype(discreteProductAxes(std::declval<T const *>()))>;

template <typename... Values>
std::true_type
namedValues(Named<Values...> const *);

std::false_type
namedValues(void const *);

template <typename T>
concept IsNamed = decltype(namedValues(std::declval<T const *>()))::value;

template <typename T>
inline constexpr bool IsUntaggedAxis = false;
template <int Cardinality, int MaxCardinality>
inline constexpr bool IsUntaggedAxis<Finite<Cardinality, MaxCardinality>>{true};
 
template <typename... Values>
inline constexpr bool IsUntaggedAxis<Named<Values...>>{true};
 
template <typename... Axes>
inline constexpr bool IsUntaggedAxis<DiscreteProduct<Axes...>>{true};


template <typename Tag, typename S>
consteval int
countAxis()
{
    if constexpr (std::is_same_v<Tag, S>)
    {
        return 1;
    }
    else if constexpr (IsDiscreteProduct<S>)
    {
        return []<typename... Axes>(std::tuple<Axes...> *) {
            return (countAxis<Tag, Axes>() + ... + 0);
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        return 0;
    }
}

template <typename Value, typename S>
consteval int
countValue()
{
    if constexpr (IsNamed<S>)
    {
        return S::template contains<Value> ? 1 : 0;
    }
    else if constexpr (IsDiscreteProduct<S>)
    {
        return []<typename... Axes>(std::tuple<Axes...> *) {
            return (countValue<Value, Axes>() + ... + 0);
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        return 0;
    }
}

template <typename S, typename F>
constexpr void
forEachAxis(S const & space, F && f)
{
    [&]<std::size_t... I>(std::index_sequence<I...>) {
        (f(space.template axis<I>(), space.template stride<I>()), ...);
    }(std::make_index_sequence<std::tuple_size_v<DiscreteProductAxes<S>>>{});
}

template <typename Tag, typename S>
[[nodiscard]] constexpr std::pair<Eigen::Index, Tag const *>
locateAxis(S const & space)
{
    if constexpr (std::is_same_v<Tag, S>)
    {
        return {1, &space};
    }
    else
    {
        static_assert(
            !IsUntaggedAxis<Tag>,
            "a path names axes by tag; derive one, such as struct Maneuver : Named<Cruise, Turn> {}"
        );
        static_assert(countAxis<Tag, S>() > 0, "the space does not have this axis");
        static_assert(countAxis<Tag, S>() < 2, "the axis appears more than once; qualify it with an enclosing tag");
        auto found = std::pair<Eigen::Index, Tag const *>{0, nullptr};
        forEachAxis(space, [&]<typename Axis>(Axis const & axis, Eigen::Index stride) {
            if constexpr (countAxis<Tag, Axis>() > 0)
            {
                auto const inner = locateAxis<Tag>(axis);
                found = {inner.first * stride, inner.second};
            }
        });
        return found;
    }
}

template <typename Tag, typename... Rest, typename S>
[[nodiscard]] constexpr auto
locateAxisPath(S const & space)
{
    auto const [stride, tag] = locateAxis<Tag>(space);
    if constexpr (sizeof...(Rest) == 0)
    {
        return std::pair{stride, tag};
    }
    else
    {
        auto const [innerStride, inner] = locateAxisPath<Rest...>(*tag);
        return std::pair{stride * innerStride, inner};
    }
}

struct ValueLocation
{
    Eigen::Index stride;       
    Eigen::Index cardinality;  
    Eigen::Index index;        

    [[nodiscard]] constexpr bool
    matches(Eigen::Index value) const
    {
        return (value / stride) % cardinality == index;
    }
};

template <typename Value, typename S>
[[nodiscard]] constexpr ValueLocation
locateValue(S const & space)
{
    static_assert(countValue<Value, S>() > 0, "the space has no value with this name");
    static_assert(
        countValue<Value, S>() < 2,
        "the value names a value of more than one axis; qualify it with a TagPath to its axis"
    );
    if constexpr (IsNamed<S>)
    {
        return {1, space.cardinality(), S::template index<Value>()};
    }
    else
    {
        auto found = ValueLocation{};
        forEachAxis(space, [&]<typename Axis>(Axis const & axis, Eigen::Index stride) {
            if constexpr (countValue<Value, Axis>() > 0)
            {
                ValueLocation const inner = locateValue<Value>(axis);
                found = {inner.stride * stride, inner.cardinality, inner.index};
            }
        });
        return found;
    }
}

template <typename P, typename S>
[[nodiscard]] constexpr ValueLocation
locateValuePath(S const & space)
{
    return withPathTags<P>([&]<typename... Tags>() {
        constexpr std::size_t count = sizeof...(Tags);
        using Value = std::tuple_element_t<count - 1, std::tuple<Tags...>>;
        if constexpr (count == 1)
        {
            return locateValue<Value>(space);
        }
        else
        {
            return [&]<std::size_t... I>(std::index_sequence<I...>) {
                auto const [stride, axis] = locateAxisPath<std::tuple_element_t<I, std::tuple<Tags...>>...>(space);
                ValueLocation inner = locateValue<Value>(*axis);
                inner.stride *= stride;
                return inner;
            }(std::make_index_sequence<count - 1>{});
        }
    });
}

template <typename P, typename S>
[[nodiscard]] constexpr auto
locateAxisByPath(S const & space)
{
    return withPathTags<P>([&]<typename... Tags>() {
        return locateAxisPath<Tags...>(space);
    });
}

template <typename... Axes>
[[nodiscard]] auto
joinAxes(std::tuple<Axes...> axes)
{
    static_assert(sizeof...(Axes) > 0, "nothing is left of the space");
    if constexpr (sizeof...(Axes) == 1)
    {
        return std::get<0>(std::move(axes));
    }
    else
    {
        return std::apply(
            [](Axes &... d) {
                return DiscreteProduct<Axes...>{std::move(d)...};
            },
            axes
        );
    }
}
 
template <typename Tag, typename... Rest, typename S>
[[nodiscard]] auto
axesWithout(S const & space);  // defined below; axisWithout() recurses into it

template <typename Tag, typename... Rest, typename Axis>
[[nodiscard]] auto
axisWithout(Axis const & axis)
{
    if constexpr (countAxis<Tag, Axis>() > 0)
    {
        return axesWithout<Tag, Rest...>(axis);
    }
    else
    {
        return std::tuple<Axis>{axis};
    }
}

template <typename Tag, typename... Rest, typename S>
[[nodiscard]] auto
axesWithout(S const & space)
{
    if constexpr (std::is_same_v<Tag, S>)
    {
        if constexpr (sizeof...(Rest) == 0)
        {
            return std::tuple<>{};
        }
        else
        {
            return axesWithout<Rest...>(space);
        }
    }
    else
    {
        static_assert(
            !IsUntaggedAxis<Tag>,
            "a path names axes by tag; derive one, such as struct Maneuver : Named<Cruise, Turn> {}"
        );
        static_assert(IsDiscreteProduct<S>, "the space does not have this axis");
        return [&]<std::size_t... I>(std::index_sequence<I...>) {
            return std::tuple_cat(axisWithout<Tag, Rest...>(space.template axis<I>())...);
        }(std::make_index_sequence<std::tuple_size_v<DiscreteProductAxes<S>>>{});
    }
}

template <typename P, typename S>
[[nodiscard]] auto
complementAxes(S const & space)
{
    return withPathTags<P>([&]<typename... Tags>() {
        return joinAxes(axesWithout<Tags...>(space));
    });
}
 
template <typename Value, typename S>
struct ValueAxisOf;

template <typename Value, typename Axes>
struct ValueAxisIn
{
    using type = void;
};

template <typename Value, typename Axis, typename... Rest>
struct ValueAxisIn<Value, std::tuple<Axis, Rest...>>
{
    using type = typename std::conditional_t<
        (countValue<Value, Axis>() > 0),
        ValueAxisOf<Value, Axis>,
        ValueAxisIn<Value, std::tuple<Rest...>>>::type;
};

template <typename Value, typename S>
struct ValueAxisOf
{
    using type = typename ValueAxisIn<Value, DiscreteProductAxes<S>>::type;
};

template <typename Value, typename S>
    requires IsNamed<S>
struct ValueAxisOf<Value, S>
{
    using type = S;
};

template <typename Value, typename S>
using ValueAxis = typename ValueAxisOf<Value, S>::type;

template <typename P, typename S>
struct ValueAxisPathOf
{
    using type = TagPath<ValueAxis<P, S>>;
};

template <typename... Tags, typename S>
struct ValueAxisPathOf<TagPath<Tags...>, S>
{
    static constexpr std::size_t Count = sizeof...(Tags);  

    using Value = std::tuple_element_t<Count - 1, std::tuple<Tags...>>;

    template <std::size_t... I>
    static consteval auto
    make(std::index_sequence<I...>)
    {
        if constexpr (Count == 1)
        {
            return TagPath<ValueAxis<Value, S>>{};
        }
        else
        {
            using Last = std::tuple_element_t<Count - 2, std::tuple<Tags...>>;
            return TagPath<std::tuple_element_t<I, std::tuple<Tags...>>..., ValueAxis<Value, Last>>{};
        }
    }

    using type = decltype(make(std::make_index_sequence<Count - 1>{}));
};

template <typename P, typename S>
using ValueAxisPath = typename ValueAxisPathOf<P, S>::type;

template <typename S>
consteval int
leafAxisCount()
{
    if constexpr (IsDiscreteProduct<S>)
    {
        return []<typename... Axes>(std::tuple<Axes...> *) {
            return (leafAxisCount<Axes>() + ... + 0);
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        return 1;
    }
}

struct LeafRange
{
    int first = 0;  
    int count = 0;  
};

template <typename Tag, typename S>
consteval LeafRange
leafRangeOfAxis()
{
    if constexpr (std::is_same_v<Tag, S>)
    {
        return {0, leafAxisCount<S>()};
    }
    else if constexpr (IsDiscreteProduct<S>)
    {
        return []<typename... Axes>(std::tuple<Axes...> *) {
            LeafRange result{};
            int offset = 0;
            auto const visit = [&]<typename Axis>() {
                if constexpr (countAxis<Tag, Axis>() > 0)
                {
                    LeafRange const inner = leafRangeOfAxis<Tag, Axis>();
                    result = {offset + inner.first, inner.count};
                }
                offset += leafAxisCount<Axis>();
            };
            (visit.template operator()<Axes>(), ...);
            return result;
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        return {};
    }
}

template <typename S, typename Tag, typename... Rest>
consteval LeafRange
leafRangeOfAxisPath()
{
    LeafRange const outer = leafRangeOfAxis<Tag, S>();
    if constexpr (sizeof...(Rest) == 0)
    {
        return outer;
    }
    else
    {
        LeafRange const inner = leafRangeOfAxisPath<Tag, Rest...>();
        return {outer.first + inner.first, inner.count};
    }
}

template <typename P, typename S>
consteval LeafRange
leafRangeOfPath()
{
    return []<typename... Tags>(TagPath<Tags...> *) {
        return leafRangeOfAxisPath<S, Tags...>();
    }(static_cast<AsTagPath<P> *>(nullptr));
}

template <typename Value, typename S>
consteval int
leafOfValue()
{
    if constexpr (IsDiscreteProduct<S>)
    {
        return []<typename... Axes>(std::tuple<Axes...> *) {
            int result = 0;
            int offset = 0;
            auto const visit = [&]<typename Axis>() {
                if constexpr (countValue<Value, Axis>() > 0)
                {
                    result = offset + leafOfValue<Value, Axis>();
                }
                offset += leafAxisCount<Axis>();
            };
            (visit.template operator()<Axes>(), ...);
            return result;
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        return 0;
    }
}

template <typename P, typename S>
consteval int
leafOfValuePath()
{
    return []<typename... Tags>(TagPath<Tags...> *) {
        constexpr std::size_t count = sizeof...(Tags);
        using Value = std::tuple_element_t<count - 1, std::tuple<Tags...>>;
        if constexpr (count == 1)
        {
            return leafOfValue<Value, S>();
        }
        else
        {
            return []<std::size_t... I>(std::index_sequence<I...>) {
                using Axis = std::tuple_element_t<count - 2, std::tuple<Tags...>>;
                return leafRangeOfAxisPath<S, std::tuple_element_t<I, std::tuple<Tags...>>...>().first
                     + leafOfValue<Value, Axis>();
            }(std::make_index_sequence<count - 1>{});
        }
    }(static_cast<AsTagPath<P> *>(nullptr));
}

template <typename S, typename Tag, typename... Rest>
consteval bool
axisPathResolves()
{
    if constexpr (!std::is_same_v<Tag, S> && (IsUntaggedAxis<Tag> || countAxis<Tag, S>() != 1))
    {
        return false;
    }
    else if constexpr (sizeof...(Rest) == 0)
    {
        return true;
    }
    else
    {
        return axisPathResolves<Tag, Rest...>();
    }
}

template <typename P, typename S>
consteval bool
pathResolves()
{
    return []<typename... Tags>(TagPath<Tags...> *) {
        return axisPathResolves<S, Tags...>();
    }(static_cast<AsTagPath<P> *>(nullptr));
}

template <typename P, typename S>
consteval bool
valuePathResolves()
{
    return []<typename... Tags>(TagPath<Tags...> *) {
        constexpr std::size_t count = sizeof...(Tags);
        using Value = std::tuple_element_t<count - 1, std::tuple<Tags...>>;
        if constexpr (count == 1)
        {
            return countValue<Value, S>() == 1;
        }
        else
        {
            return []<std::size_t... I>(std::index_sequence<I...>) {
                using Axis = std::tuple_element_t<count - 2, std::tuple<Tags...>>;
                if constexpr (!axisPathResolves<S, std::tuple_element_t<I, std::tuple<Tags...>>...>())
                {
                    return false;
                }
                else
                {
                    return countValue<Value, Axis>() == 1;
                }
            }(std::make_index_sequence<count - 1>{});
        }
    }(static_cast<AsTagPath<P> *>(nullptr));
}

template <typename S>
consteval std::array<bool, std::size_t(leafAxisCount<S>())>
namedLeaves()
{
    auto result = std::array<bool, std::size_t(leafAxisCount<S>())>{};
    if constexpr (IsDiscreteProduct<S>)
    {
        [&]<typename... Axes>(std::tuple<Axes...> *) {
            std::size_t offset = 0;
            auto const append = [&]<typename Axis>() {
                for (bool const named : namedLeaves<Axis>())
                {
                    result[offset++] = named;
                }
            };
            (append.template operator()<Axes>(), ...);
        }(static_cast<DiscreteProductAxes<S> *>(nullptr));
    }
    else
    {
        result[0] = IsNamed<S>;
    }
    return result;
}

enum class Naming
{
    complete,        
    unresolved,      
    repeated,        
    missingNamed,    
    missingUnnamed,  
};

template <std::size_t N>
consteval void
countLeaves(std::array<int, N> & times, LeafRange range)
{
    for (int k = range.first; k < range.first + range.count; ++k)
    {
        ++times[std::size_t(k)];
    }
}

template <typename S, typename... Values, typename... Paths>
consteval Naming
namingOf(std::tuple<Values...> *, std::tuple<Paths...> *)
{
    if constexpr (!((valuePathResolves<Values, S>() && ...) && (pathResolves<Paths, S>() && ...)))
    {
        return Naming::unresolved;
    }
    else
    {
        auto times = std::array<int, std::size_t(leafAxisCount<S>())>{};
        (countLeaves(times, LeafRange{leafOfValuePath<Values, S>(), 1}), ...);
        (countLeaves(times, leafRangeOfPath<Paths, S>()), ...);
        constexpr auto named = namedLeaves<S>();
        for (int const n : times)
        {
            if (n > 1)
            {
                return Naming::repeated;
            }
        }
        for (std::size_t k = 0; k < times.size(); ++k)
        {
            if (times[k] == 0 && named[k])
            {
                return Naming::missingNamed;
            }
        }
        for (int const n : times)
        {
            if (n == 0)
            {
                return Naming::missingUnnamed;
            }
        }
        return Naming::complete;
    }
}

template <typename P>
struct LastTagOf;

template <typename... Tags>
struct LastTagOf<TagPath<Tags...>>
{
    using type = std::tuple_element_t<sizeof...(Tags) - 1, std::tuple<Tags...>>;
};

template <typename P>
using LastTag = typename LastTagOf<AsTagPath<P>>::type;

template <typename S>
[[nodiscard]] constexpr std::array<Eigen::Index, std::size_t(leafAxisCount<S>())>
leafCardinalities(S const & space)
{
    auto result = std::array<Eigen::Index, std::size_t(leafAxisCount<S>())>{};
    if constexpr (IsDiscreteProduct<S>)
    {
        std::size_t offset = 0;
        forEachAxis(space, [&]<typename Axis>(Axis const & axis, Eigen::Index) {
            for (Eigen::Index const cardinality : leafCardinalities(axis))
            {
                result[offset++] = cardinality;
            }
        });
    }
    else
    {
        result[0] = space.cardinality();
    }
    return result;
}

template <typename P, typename S>
[[nodiscard]] constexpr Eigen::Index
namedOffset(S const & space)
{
    ValueLocation const location = locateValuePath<P>(space);
    return location.index * location.stride;
}
 
}  // namespace internal

template <typename Path, std::size_t N>
struct AxisIndex
{
    std::array<Eigen::Index, N> indices;  
};

template <typename Path, std::integral... I>
[[nodiscard]] constexpr AxisIndex<Path, sizeof...(I)>
at(I... indices)
{
    using Axis = internal::LastTag<Path>;
    static_assert(DiscreteSpace<Axis>, "at<Axis>() names an axis of a discrete space");
    static_assert(
        sizeof...(I) == std::size_t(internal::leafAxisCount<Axis>()),
        "give one index for each leaf axis of the axis named"
    );
    return {{Eigen::Index(indices)...}};
}
 
namespace internal
{

template <typename Path, std::size_t N, typename S>
[[nodiscard]] constexpr Eigen::Index
indexedOffset(S const & space, AxisIndex<Path, N> const & index)
{
    auto const [stride, axis] = locateAxisByPath<Path>(space);
    auto const cardinalities = leafCardinalities(*axis);
    Eigen::Index flat = 0;
    for (std::size_t k = 0; k < N; ++k)
    {
        EMRYS_FILTER_ASSERT(index.indices[k] >= 0 && index.indices[k] < cardinalities[k] && "an index is out of range");
        flat = flat * cardinalities[k] + index.indices[k];
    }
    return flat * stride;
}
 
}  // namespace internal

template <typename... Values, DiscreteSpace S, typename... Paths, std::size_t... N>
[[nodiscard]] constexpr Eigen::Index
valueOf(S const & space, AxisIndex<Paths, N> const &... indices)
{
    constexpr bool named = sizeof...(Values) + sizeof...(Paths) > 0;
    static_assert(named, "name a value");
    constexpr internal::Naming naming = !named ? internal::Naming::unresolved
                                               : internal::namingOf<S>(
                                                     static_cast<std::tuple<Values...> *>(nullptr),
                                                     static_cast<std::tuple<Paths...> *>(nullptr)
                                                 );
    static_assert(naming != internal::Naming::repeated, "an axis is named more than once");
    static_assert(naming != internal::Naming::missingNamed, "name one value in each axis");
    static_assert(naming != internal::Naming::missingUnnamed, "an axis without names takes its index from at<Axis>(i)");
    return (Eigen::Index{0} + ... + internal::namedOffset<Values>(space))
         + (Eigen::Index{0} + ... + internal::indexedOffset(space, indices));
}

template <DiscreteSpace S, typename... Values, typename... Paths, std::size_t... N>
    requires FixedSizeSpace<S>
[[nodiscard]] constexpr Eigen::Index
valueOf(AxisIndex<Paths, N> const &... indices)
{
    return valueOf<Values...>(S{}, indices...);
}
 
}  // namespace Emrys