Composite.h

Go to the documentation of this file.

#pragma once

 
#include <array>
#include <cstddef>
#include <limits>
#include <numeric>
#include <optional>
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Spaces/Coordinates.h>
#include <EmrysFilter/Spaces/Space.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Part>
using StoredPart = std::conditional_t<std::is_empty_v<Part>, std::tuple<>, std::tuple<Part>>;

template <typename Stored>
struct StoredParts
{
    Stored parts;  

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

template <>
struct StoredParts<std::tuple<>>
{
    StoredParts() = default;

    constexpr explicit StoredParts(std::tuple<>) {}

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

template <typename... Parts>
class PartStorage : private StoredParts<decltype(std::tuple_cat(std::declval<StoredPart<Parts>>()...))>
{
    using Stored = decltype(std::tuple_cat(std::declval<StoredPart<Parts>>()...));

    using Base = StoredParts<Stored>;

    template <std::size_t I>
    using Part = std::tuple_element_t<I, std::tuple<Parts...>>;

    static consteval std::size_t
    storedIndex(std::size_t i)
    {
        constexpr auto empty = std::array<bool, sizeof...(Parts)>{std::is_empty_v<Parts>...};
        std::size_t result = 0;
        for (std::size_t j = 0; j < i; ++j)
        {
            result += empty[j] ? 0 : 1;
        }
        return result;
    }

    template <typename P>
    static constexpr StoredPart<P>
    keep(P & part)
    {
        if constexpr (std::is_empty_v<P>)
        {
            return {};
        }
        else
        {
            return StoredPart<P>{std::move(part)};
        }
    }

    template <std::size_t I>
    static inline Part<I> instance{};
 
public:
 
    PartStorage() = default;

    constexpr explicit PartStorage(Parts... parts)
        : Base{std::tuple_cat(keep(parts)...)}
    {}

    template <std::size_t I>
    [[nodiscard]] constexpr Part<I> const &
    part() const
    {
        if constexpr (std::is_empty_v<Part<I>>)
        {
            return instance<I>;
        }
        else
        {
            return std::get<storedIndex(I)>(Base::parts);
        }
    }

    template <std::size_t I>
    [[nodiscard]] Part<I> &
    part()
    {
        if constexpr (std::is_empty_v<Part<I>>)
        {
            return instance<I>;
        }
        else
        {
            return std::get<storedIndex(I)>(Base::parts);
        }
    }

    friend bool
    operator==(PartStorage const &, PartStorage const &) = default;
};
 
}  // namespace internal

template <Space... Blocks>
class Composite : private internal::PartStorage<Blocks...>
{
    using Storage = internal::PartStorage<Blocks...>;
 
    static_assert(sizeof...(Blocks) > 0, "a composite needs at least one block");
    static constexpr std::size_t BlockCount = sizeof...(Blocks);  
    static constexpr std::array<bool, BlockCount> Resizable{
        (Blocks::ValueDimensionStep != 0)...
    };  

    static consteval int
    strides()
    {
        int result = 0;
        ((result = std::gcd(result, Blocks::ValueDimensionStep)), ...);
        return result;
    }

    static consteval std::size_t
    firstResizable()
    {
        for (std::size_t i = 0; i < BlockCount; ++i)
        {
            if (Resizable[i])
            {
                return i;
            }
        }
        return BlockCount;
    }
 
public:
 
    static constexpr int FixedValueDimension =
        (Blocks::FixedValueDimension + ...);  

    static constexpr int ValueDimensionStep = strides();

    static constexpr int MaxValueDimension =
        ((Blocks::MaxValueDimension == Eigen::Dynamic) || ...) ? Eigen::Dynamic : (Blocks::MaxValueDimension + ...);
 
    Composite() = default;

    explicit Composite(Blocks... blocks)
        : Storage{std::move(blocks)...}
    {}

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

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

    [[nodiscard]] Eigen::Index
    valueDimension() const
    {
        return offset<BlockCount>();
    }

    template <std::size_t I>
    [[nodiscard]] Eigen::Index
    offset() const
    {
        return [&]<std::size_t... J>(std::index_sequence<J...>) {
            return (Eigen::Index{0} + ... + block<J>().valueDimension());
        }(std::make_index_sequence<I>{});
    }

    void
    resize(Eigen::Index n)
        requires((internal::resizableBlockCount<Blocks> + ...) == 1)
    {
        auto & resizable = block<firstResizable()>();
        Eigen::Index const rest = valueDimension() - resizable.valueDimension();
        EMRYS_FILTER_ASSERT(n >= rest && "the dimension does not fit the space");
        resizable.resize(n - rest);
    }

    template <typename Derived>
    void
    normalize(Eigen::MatrixBase<Derived> const & rows) const
        requires((internal::HasNormalize<Blocks> || ...))
    {
        forEachBlock(rows, [](auto const & block, auto && sub) {
            Emrys::normalizeValues(block, sub);
        });
    }

    template <typename X, typename D, typename Out>
    void
    plus(Eigen::MatrixBase<X> const & x, Eigen::MatrixBase<D> const & d, Eigen::MatrixBase<Out> const & out) const
        requires((internal::HasPlus<Blocks> || ...))
    {
        forEachBlockOf<false, true, false>(
            x,
            d,
            out,
            [](auto const & block, auto const & xs, auto const & ds, auto const & os) {
                internal::plusInto(block, xs, ds, os);
            }
        );
    }

    template <typename A, typename B, typename Out>
    void
    minus(Eigen::MatrixBase<A> const & a, Eigen::MatrixBase<B> const & b, Eigen::MatrixBase<Out> const & out) const
        requires((internal::HasMinus<Blocks> || ...))
    {
        forEachBlockOf<false, false, true>(
            a,
            b,
            out,
            [](auto const & block, auto const & as, auto const & bs, auto const & os) {
                internal::minusInto(block, as, bs, os);
            }
        );
    }

    [[nodiscard]] Composite<DifferenceSpace<Blocks>...>
    difference() const
    {
        return [&]<std::size_t... I>(std::index_sequence<I...>) {
            return Composite<DifferenceSpace<Blocks>...>{differenceSpace(block<I>())...};
        }(std::make_index_sequence<BlockCount>{});
    }

    friend bool
    operator==(Composite const &, Composite const &) = default;
 
private:

    template <typename Derived, typename F>
    void
    forEachBlock(Eigen::MatrixBase<Derived> const & rows, F f) const
    {
        auto & m = writable(rows);
        EMRYS_FILTER_ASSERT(m.rows() == valueDimension() && "the rows do not match the space");
        Eigen::Index offset = 0;
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((f(block<I>(),
                internal::middleRowsOf<valueDimensionAtCompileTime<std::tuple_element_t<I, std::tuple<Blocks...>>>>(
                    m,
                    offset,
                    block<I>().valueDimension()
                )),
              offset += block<I>().valueDimension()),
             ...);
        }(std::make_index_sequence<BlockCount>{});
    }

    template <
        bool FirstDifference,
        bool SecondDifference,
        bool OutDifference,
        typename First,
        typename Second,
        typename Out,
        typename F>
    void
    forEachBlockOf(
        Eigen::MatrixBase<First> const & first,
        Eigen::MatrixBase<Second> const & second,
        Eigen::MatrixBase<Out> const & out,
        F f
    ) const
    {
        auto & o = writable(out);
        Eigen::Index firstOffset = 0;
        Eigen::Index secondOffset = 0;
        Eigen::Index outOffset = 0;
        auto const step = [&]<typename B>(B const & b) {
            Eigen::Index const firstRows = internal::rowsOf<FirstDifference>(b);
            Eigen::Index const secondRows = internal::rowsOf<SecondDifference>(b);
            Eigen::Index const outRows = internal::rowsOf<OutDifference>(b);
            f(b,
              internal::middleRowsOf<internal::rowsAtCompileTime<FirstDifference, B>>(first, firstOffset, firstRows),
              internal::middleRowsOf<internal::rowsAtCompileTime<SecondDifference, B>>(
                  second,
                  secondOffset,
                  secondRows
              ),
              internal::middleRowsOf<internal::rowsAtCompileTime<OutDifference, B>>(o, outOffset, outRows));
            firstOffset += firstRows;
            secondOffset += secondRows;
            outOffset += outRows;
        };
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            (step(block<I>()), ...);
        }(std::make_index_sequence<BlockCount>{});
        EMRYS_FILTER_ASSERT(
            firstOffset == first.rows() && secondOffset == second.rows() && outOffset == o.rows()
            && "the rows do not match the space"
        );
    }
};

template <Space Block, int MaxCount_ = Eigen::Dynamic>
class Repeated
{
    static_assert(
        Block::ValueDimensionStep == 0 && Block::FixedValueDimension > 0,
        "a repeated block must have a fixed, nonzero size"
    );
    static_assert(MaxCount_ == Eigen::Dynamic || MaxCount_ >= 0);
 
public:

    using BlockType = Block;

    static constexpr int MaxCount = MaxCount_;
 
    static constexpr int FixedValueDimension = 0;  
    static constexpr int ValueDimensionStep = Block::FixedValueDimension;  

    static constexpr int MaxValueDimension =
        MaxCount == Eigen::Dynamic ? Eigen::Dynamic : MaxCount * Block::FixedValueDimension;
 
    Repeated() = default;

    explicit Repeated(Block block, Eigen::Index count = 0)
        : m_block{std::move(block)}
    {
        setCount(count);
    }

    [[nodiscard]] Block const &
    block() const
    {
        return m_block;
    }

    [[nodiscard]] Eigen::Index
    count() const
    {
        return m_count;
    }

    [[nodiscard]] Eigen::Index
    valueDimension() const
    {
        return m_count * ValueDimensionStep;
    }

    void
    setCount(Eigen::Index count)
    {
        EMRYS_FILTER_ASSERT(
            count >= 0 && (MaxCount == Eigen::Dynamic || count <= MaxCount) && "the dimension does not fit the space"
        );
        m_count = count;
    }

    void
    resize(Eigen::Index n)
    {
        EMRYS_FILTER_ASSERT(n >= 0 && n % ValueDimensionStep == 0 && "the dimension does not fit the space");
        setCount(n / ValueDimensionStep);
    }

    template <typename Derived>
    void
    normalize(Eigen::MatrixBase<Derived> const & rows) const
        requires(internal::HasNormalize<Block>)
    {
        forEachRepetition(rows, [this](auto && sub) {
            Emrys::normalizeValues(m_block, sub);
        });
    }

    template <typename X, typename D, typename Out>
    void
    plus(Eigen::MatrixBase<X> const & x, Eigen::MatrixBase<D> const & d, Eigen::MatrixBase<Out> const & out) const
        requires(internal::HasPlus<Block>)
    {
        forEachRepetitionOf<false, true, false>(x, d, out, [this](auto const & xs, auto const & ds, auto const & os) {
            internal::plusInto(m_block, xs, ds, os);
        });
    }

    template <typename A, typename B, typename Out>
    void
    minus(Eigen::MatrixBase<A> const & a, Eigen::MatrixBase<B> const & b, Eigen::MatrixBase<Out> const & out) const
        requires(internal::HasMinus<Block>)
    {
        forEachRepetitionOf<false, false, true>(a, b, out, [this](auto const & as, auto const & bs, auto const & os) {
            internal::minusInto(m_block, as, bs, os);
        });
    }

    [[nodiscard]] Repeated<DifferenceSpace<Block>, MaxCount_>
    difference() const
    {
        return Repeated<DifferenceSpace<Block>, MaxCount_>{differenceSpace(m_block), m_count};
    }

    friend bool
    operator==(Repeated const &, Repeated const &) = default;
 
private:

    template <typename Derived, typename F>
    void
    forEachRepetition(Eigen::MatrixBase<Derived> const & rows, F f) const
    {
        auto & m = writable(rows);
        EMRYS_FILTER_ASSERT(m.rows() == valueDimension() && "the rows do not match the space");
        for (Eigen::Index k = 0; k < m_count; ++k)
        {
            f(internal::middleRowsOf<ValueDimensionStep>(m, k * ValueDimensionStep, ValueDimensionStep));
        }
    }

    template <
        bool FirstDifference,
        bool SecondDifference,
        bool OutDifference,
        typename First,
        typename Second,
        typename Out,
        typename F>
    void
    forEachRepetitionOf(
        Eigen::MatrixBase<First> const & first,
        Eigen::MatrixBase<Second> const & second,
        Eigen::MatrixBase<Out> const & out,
        F f
    ) const
    {
        auto & o = writable(out);
        Eigen::Index const firstRows = internal::rowsOf<FirstDifference>(m_block);
        Eigen::Index const secondRows = internal::rowsOf<SecondDifference>(m_block);
        Eigen::Index const outRows = internal::rowsOf<OutDifference>(m_block);
        EMRYS_FILTER_ASSERT(
            first.rows() == m_count * firstRows && second.rows() == m_count * secondRows
            && o.rows() == m_count * outRows && "the rows do not match the space"
        );
        for (Eigen::Index k = 0; k < m_count; ++k)
        {
            f(internal::middleRowsOf<internal::rowsAtCompileTime<FirstDifference, Block>>(
                  first,
                  k * firstRows,
                  firstRows
              ),
              internal::middleRowsOf<internal::rowsAtCompileTime<SecondDifference, Block>>(
                  second,
                  k * secondRows,
                  secondRows
              ),
              internal::middleRowsOf<internal::rowsAtCompileTime<OutDifference, Block>>(o, k * outRows, outRows));
        }
    }
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Block m_block;  
    Eigen::Index m_count = 0;                      
};
 
namespace internal
{

template <typename... Blocks>
std::tuple<Blocks...> *
compositeBlocks(Composite<Blocks...> const *);

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

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

template <typename Blocks>
inline constexpr int resizableBlockCountIn = 0;

template <typename... Blocks>
inline constexpr int resizableBlockCountIn<std::tuple<Blocks...>> = (0 + ... + resizableBlockCount<Blocks>);

template <IsComposite S>
inline constexpr int resizableBlockCount<S> = resizableBlockCountIn<CompositeBlocks<S>>;

template <typename S>
concept SizedByLength = resizableBlockCount<S> <= 1;

template <typename Block, int MaxCount>
Block *
repeatedBlock(Repeated<Block, MaxCount> const *);

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

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

template <int Dimension, int MaxDimension>
std::true_type
leafExtent(LeafExtent<Dimension, MaxDimension> const *);

std::false_type
leafExtent(void const *);

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

template <PeriodicInterval Interval, int Dimension, int MaxDimension>
std::true_type
periodicLeaf(Periodic<Interval, Dimension, MaxDimension> const *);

std::false_type
periodicLeaf(void const *);

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

struct LeafSupport
{
    double lower = -std::numeric_limits<double>::infinity();  
    double upper = std::numeric_limits<double>::infinity();   
    bool periodic = false;                                    
};

// constexpr rather than consteval: MSVC 14.50 can't apply -> to a consteval call's optional result.
template <typename Leaf>
[[nodiscard]] constexpr std::optional<LeafSupport>
leafSupport()
{
    if constexpr (IsPeriodic<Leaf>)
    {
        return LeafSupport{.lower = Leaf::lower(), .upper = Leaf::upper(), .periodic = true};
    }
    else if constexpr (isLinear<Leaf>)
    {
        return LeafSupport{};
    }
    else
    {
        return std::nullopt;
    }
}

template <typename Space, typename F>
void
forEachLeaf(Space const & space, F && f, Eigen::Index offset = 0)
{
    if constexpr (IsComposite<Space>)
    {
        [&]<std::size_t... I>(std::index_sequence<I...>) {
            ((forEachLeaf(space.template block<I>(), f, offset), offset += space.template block<I>().valueDimension()),
             ...);
        }(std::make_index_sequence<std::tuple_size_v<CompositeBlocks<Space>>>{});
    }
    else if constexpr (IsRepeated<Space>)
    {
        for (Eigen::Index k = 0; k < space.count(); ++k)
        {
            forEachLeaf(space.block(), f, offset + k * Space::ValueDimensionStep);
        }
    }
    else
    {
        f(offset, space);
    }
}
 
}  // namespace internal
 
}  // namespace Emrys