
/// =====================================================================
/// !! DO NOT EDIT THIS FILE !!
/// This file was automatically generated on Mon Oct  7 16:20:55 CEST 2019
/// =====================================================================

/// MIT License
///
/// Copyright (c) 2015 Sebastian Dalleiger
///
/// Permission is hereby granted, free of charge, to any person obtaining a copy
/// of this software and associated documentation files (the "Software"), to deal
/// in the Software without restriction, including without limitation the rights
/// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
/// copies of the Software, and to permit persons to whom the Software is
/// furnished to do so, subject to the following conditions:
///
/// The above copyright notice and this permission notice shall be included in all
/// copies or substantial portions of the Software.
///
/// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
/// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
/// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
/// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
/// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
/// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
/// SOFTWARE.

/// =====================================================================
/// README
/// =====================================================================

/// # ndarray
///
/// ## Abstract: A library for multi-dimensional strided views (slices) and arrays for C++14/17
/// under MIT licence _without_ linear algebra.
///
/// ```c++
/// std::vector<double> vec(n * m * o);
/// sd::slice xs (vec);
/// assert(std::equals(xs.begin(), xs.end(), x.begin(), x.end()));
///
/// sd::slice sup = xs.cut(a, b).stride(2);
/// for (auto t : sub.iter()) {}
///
/// for (auto i : sub.subscripts()) {
///     auto t = sup[i];
/// }
///
/// sd::slice ts (vec, n, m, o); // sd::slice<double, 3>(vec, {n, m, o});
/// for (auto t : ts) {}
/// for (auto t : ts.iter()) {}
/// for (auto i : ts.subscripts()) {
///     // auto [s, t, u] = i;
///     // ts(s, t, u);
///     auto t = sup[i];
/// }
/// ```
///
///
/// # This library...
/// - provides multi-dimensional strided views (slice) and arrays (ndarray, fixed_ndarray)
/// - is similar to boost.multi-array
/// - is constexpr-aware
/// - uses C++14 and is aware of C++17 deduction guides and structured binding
/// - is lightweight, well-tested and used in practice by myself :)
/// - is easy on the compiler, which are able to generate efficient code. (see godbolt.org)
/// - is _header-only_ and has a _single-header_ variant, which is automatically generated using
/// an ugly bash script.
/// - has cmake support
///
/// # This library does __not__ provide
/// - math, linear algebra à la BLAS, Eigen, numpy, std::valarray and expression templates.
/// - built-in distributed computing facilities like MPI <!-- dash-project.org or
/// stellar-group.org/libraries/hpx  -->
/// - runtime bounds checking _by default_: It can be enabled by defining the macro
/// ```sd_expect``` before including the header for the first time.
///     - e.g.: ```sd_expect(BOOLEAN_EXPRESSION) if(!(BOOLEAN_EXPRESSION)) throw
///     std::runtime_error{"unexpected behaviour in" __FILE__ ":" __LINE__" (multi_array)"}```
///
/// # Dependencies and Installation
/// - C++14/17 compatible compiler (MSVC2017 works)
/// - Optionally: cmake and bash
///
/// ## Single Header
/// The generation depends on a hacky shell script: ```./scripts/meld.sh``` or just use the
/// custom target
/// ```sh
/// cmake --build ${build_dir} --target marray-single-header
/// ```
///
/// # Usage
///
/// ## Construction
///
/// ## Iterators
/// ### Index Ranges
/// ```c++
/// s.subscripts()
/// s.subscripts<1,0>() ->
/// s.indices() -> {1, 2, 3, ...} or strided {1, 3, 5, ...}
///
///
/// sd::layout<3> geo = t.shape();
/// for (auto idx : t.subscripts()) {
///     auto [i, j, k] = idx;
///     t(idx) == t[idx] == t[{i, j, k}] == t(i, j, k) == x[geo.linearize(idx)] ==
///     s(geo.linearize(idx) == s[geo.unlinearize(geo.linearize(idx))];
/// }
///
/// sd::ndarray y (0.0, n, m, o); // sd::ndarray<double, 3>(0.0, {n, m, o});
/// for (auto idx : y.subscripts<1, 0, 2>()) {
///     // 0 0 0
///     // 0 1 0
///     // 0 2 0
///     // ...
///     // 0 m 0
///     // 1 0 0
///     // 1 1 0
///     // ... etc.
/// }
///
/// ```
/// ### Strided Iterators
///
/// ## Slicing and Striding
/// ```c++
/// auto t0 = s[i];
/// auto t3 = s.cut({});
/// auto t2 = s.cut({}, {});
/// auto t1 = s.cut({}, {}, {});
/// auto s1 = s.flatten();
/// auto s2 = s.as({n1, m2, o3});
/// ```
///
/// ## ```const_propagating_slice``` (ndarray, fixed_ndarray)
/// ```c++
/// static_assert(!std::is_const_v<decltype(viewable)>);
///
/// const_propagating_slice<T> xs({viewable});
/// std::as_const(xs)[0] = T{}; // won't compile (good)
/// slice<T> ys(viewable);
/// std::as_const(ys)[0] = T{}; // will compile (bad!)
/// ```
/// (a hypothetical operator. would simplify the implementation a lot, btw)
///
/// <!-- Objects of this type will pass calls to member functions to the const-correct view-type
/// using ```::get()``` (or ```operator*, operator->```)
/// ```c++
///     auto function_name() const { return get().function_name(); }
///     auto function_name()       { return get().function_name(); }
/// ```
///
/// ```c++
/// const_propagating_slice<T, N> cps;
///
/// static_assert(std::is_same_v<decltype(              cps .get()), slice<      T, N>>);
/// static_assert(std::is_same_v<decltype(std::as_const(cps).get()), slice<const T, N>>);
/// static_assert(std::is_same_v<decltype(std::as_const(cps).get()), slice<const T, N>>);
/// ``` -->
/// ## Interface of ```slice```
/// ## Interface of ```layout```
///

#pragma once
#ifndef MULTI_ARRAY_GENERATED_SINGLE_HEADER_HXX
#define MULTI_ARRAY_GENERATED_SINGLE_HEADER_HXX

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <iterator>
#include <memory>
#include <numeric>
#include <tuple>
#include <type_traits>
#include <vector>
/// See "LICENSE" for the full license governing this code.
#ifndef sd_expect
#define sd_expect(_Expression) assert(_Expression)
#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_ALGORITHM_HXX
#define MARRAY_ALGORITHM_HXX

namespace sd
{
namespace stdext
{
//
// this exists because I happen to need constexpr algorithms.
//

template <class InputIt, class OutputIt, class UnaryOperation>
constexpr OutputIt
transform(InputIt first1, InputIt last1, OutputIt d_first, UnaryOperation unary_op)
{
    while (first1 != last1) { *d_first++ = unary_op(*first1++); }
    return d_first;
}

template <class InputIt1, class InputIt2, class OutputIt, class BinaryOperation>
constexpr OutputIt transform(InputIt1        first1,
                             InputIt1        last1,
                             InputIt2        first2,
                             OutputIt        d_first,
                             BinaryOperation binary_op)
{
    while (first1 != last1) { *d_first++ = binary_op(*first1++, *first2++); }
    return d_first;
}

template <class InputIt1, class InputIt2, class BinaryPredicate>
constexpr bool equal(InputIt1 first1, InputIt1 last1, InputIt2 first2, BinaryPredicate p)
{
    for (; first1 != last1; ++first1, ++first2)
    {
        if (!p(*first1, *first2)) { return false; }
    }
    return true;
}

template <class InputIt, class OutputIt>
constexpr OutputIt copy(InputIt first, InputIt last, OutputIt d_first)
{
    while (first != last) { *d_first++ = *first++; }
    return d_first;
}

template <class InputIt, class T>
constexpr void fill(InputIt first, InputIt last, T value)
{
    while (first != last) { *first++ = value; }
}
} // namespace stdext
} // namespace sd

#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_ARRAY_HXX
#define MARRAY_ARRAY_HXX

namespace sd
{

template <typename S, typename... Ts>
class are_convertible : public std::integral_constant<bool, true>
{
};

template <typename S, typename T, typename... Ts>
class are_convertible<S, T, Ts...>
    : public std::integral_constant<bool,
                                    std::is_convertible<S, T>::value &&
                                        are_convertible<S, Ts...>::value>

{
};

// template <typename S, typename... Args>
// class are_convertible
//     : public std::integral_constant<bool,
//                                     std::conjunction_v<std::is_convertible_v<S, Args>...>> {
// };

template <typename T, std::size_t rank, typename... Ts>
class is_valid_array
    : public std::integral_constant<bool,
                                    sizeof...(Ts) == rank && are_convertible<T, Ts...>::value>
{
};

// array is like a std::array but has
//    (1) constructors (which is not really important)
//    (2) for fixed dimension: array is a (i) partially ordered, (ii) Ring/Field-like
//    algebraic structure (I need coordinate-wise aritmetics on my arrays)
//    (3) additional optional tests
template <typename T, size_t Rank>
class array
{
public:
    static constexpr const size_t rank = Rank;

    using value_type      = T;
    using size_type       = size_t;
    using reference       = std::add_lvalue_reference_t<value_type>;
    using const_reference = std::add_lvalue_reference_t<std::add_const_t<value_type>>;

    constexpr array() noexcept                   = default;
    constexpr array(const array& other) noexcept = default;
    constexpr array(array&& other) noexcept      = default;
    constexpr array& operator=(const array& rhs) noexcept = default;
    constexpr array& operator=(array&& rhs) noexcept = default;

    template <typename... Ts,
              typename = std::enable_if_t<is_valid_array<value_type, rank, Ts...>::value>>
    constexpr array(const Ts&... ds) noexcept : elems{static_cast<value_type>(ds)...}
    {
    }

    constexpr array(const value_type (&values)[rank]) noexcept
    {
        stdext::copy(values, values + rank, elems);
    }

    constexpr array(const std::array<value_type, rank>& values) noexcept
    {
        stdext::copy(std::begin(values), std::end(values), elems);
    }

    constexpr array(const value_type* first, const value_type* last) noexcept
    {
        sd_expect(last - first == rank);
        stdext::copy(first, last, elems);
    }

    constexpr reference operator[](size_type component_idx) noexcept
    {
        sd_expect(component_idx < rank);
        return elems[component_idx];
    }

    constexpr const_reference operator[](size_type component_idx) const noexcept
    {
        sd_expect(component_idx < rank);
        return elems[component_idx];
    }

    constexpr reference operator()(size_type component_idx) noexcept
    {
        sd_expect(component_idx < rank);
        return elems[component_idx];
    }

    constexpr const_reference operator()(size_type component_idx) const noexcept
    {
        sd_expect(component_idx < rank);
        return elems[component_idx];
    }

    template <typename Predicate>
    constexpr bool compare(const array& rhs, Predicate&& p) const noexcept
    {
        return stdext::equal(begin(), end(), rhs.begin(), p);
    }

    constexpr bool operator==(const array& rhs) const noexcept
    {
        return compare(rhs, std::equal_to<value_type>{});
    }

    constexpr bool operator!=(const array& rhs) const noexcept { return !(*this == rhs); }

    constexpr bool operator<=(const array& rhs) const noexcept
    {
        return compare(rhs, std::less_equal<value_type>{});
    }

    constexpr bool operator<(const array& rhs) const noexcept
    {
        return compare(rhs, std::less<value_type>{});
    }

    constexpr bool operator>=(const array& rhs) const noexcept
    {
        return compare(rhs, std::greater_equal<value_type>{});
    }

    constexpr bool operator>(const array& rhs) const noexcept
    {
        return compare(rhs, std::greater<value_type>{});
    }

    constexpr array operator+() const noexcept { return *this; }

    constexpr array operator-() const noexcept
    {
        array ret = *this;
        stdext::transform(ret.begin(), ret.end(), ret.begin(), std::negate<value_type>{});
        return ret;
    }

    constexpr array operator+(const array& rhs) const noexcept
    {
        array ret = *this;
        ret += rhs;
        return ret;
    }

    constexpr array operator-(const array& rhs) const noexcept
    {
        array ret = *this;
        ret -= rhs;
        return ret;
    }

    constexpr array& operator+=(const array& rhs) noexcept
    {
        stdext::transform(begin(), end(), rhs.begin(), begin(), std::plus<value_type>{});
        return *this;
    }

    constexpr array& operator-=(const array& rhs) noexcept
    {
        stdext::transform(begin(), end(), rhs.begin(), begin(), std::minus<value_type>{});
        return *this;
    }

    constexpr array operator*(const array& rhs) const noexcept
    {
        array ret = *this;
        ret *= rhs;
        return ret;
    }

    constexpr array operator/(const array& rhs) const noexcept
    {
        array ret = *this;
        ret /= rhs;
        return ret;
    }

    constexpr array& operator*=(const array& rhs) noexcept
    {
        stdext::transform(begin(), end(), rhs.begin(), begin(), std::multiplies<value_type>{});
        return *this;
    }

    constexpr array& operator/=(const array& rhs) noexcept
    {
        stdext::transform(begin(), end(), rhs.begin(), begin(), std::divides<value_type>{});
        return *this;
    }

    constexpr array operator*(value_type v) const noexcept
    {
        array ret = *this;
        ret *= v;
        return ret;
    }

    constexpr array operator/(value_type v) const noexcept
    {
        array ret = *this;
        ret /= v;
        return ret;
    }

    friend constexpr array operator*(value_type v, const array& rhs) noexcept
    {
        return rhs * v;
    }

    constexpr array& operator*=(value_type v) noexcept
    {
        stdext::transform(begin(),
                          end(),
                          begin(),
                          [v](value_type x) { return std::multiplies<value_type>{}(x, v); });
        return *this;
    }

    constexpr array& operator/=(value_type v) noexcept
    {
        stdext::transform(begin(),
                          end(),
                          begin(),
                          [v](value_type x) { return std::divides<value_type>{}(x, v); });
        return *this;
    }

    constexpr array& fill(value_type v) noexcept
    {
        stdext::fill(begin(), end(), v);
        return *this;
    }

    constexpr size_type size() const noexcept { return rank; }

    constexpr value_type* begin() noexcept { return elems; }

    constexpr value_type* end() noexcept { return elems + size(); }

    constexpr const value_type* begin() const noexcept { return elems; }

    constexpr const value_type* end() const noexcept { return elems + size(); }

    constexpr const value_type* data() const { return begin(); }

    constexpr value_type* data() { return begin(); }

private:
    value_type elems[rank] = {};

    static_assert(Rank > 0, "Rank must be greater than 0!");
    template <typename U, size_t OtherRank>
    friend class array;
};

// template <size_t Rank>
// using subscript = array<size_t, Rank>;

#ifdef __cpp_deduction_guides
template <typename... Ts>
array(const Ts&... ds) -> array<std::common_type_t<Ts...>, sizeof...(Ts)>;
#endif

#ifdef __cpp_structured_bindings
template <size_t Idx, class T, size_t Rank>
inline constexpr T& get(array<T, Rank>& arr) noexcept
{
    static_assert(Idx < Rank, "Index out of bounds in std::get<> (sd::array)");
    return arr.data()[Idx];
}

template <size_t Idx, class T, size_t Rank>
inline constexpr const T& get(const array<T, Rank>& arr) noexcept
{
    static_assert(Idx < Rank, "Index out of bounds in std::get<> (const sd::array)");
    return arr.data()[Idx];
}

template <size_t Idx, class T, size_t Rank>
inline constexpr T&& get(array<T, Rank>&& arr) noexcept
{
    static_assert(Idx < Rank, "Index out of bounds in std::get<> (const sd::array)");
    return std::move(arr.data()[Idx]);
}

template <size_t Idx, class T, size_t Rank>
inline constexpr const T&& get(const array<T, Rank>&& arr) noexcept
{
    static_assert(Idx < Rank, "Index out of bounds in std::get<> (const sd::array)");
    return std::move(arr.data()[Idx]);
}
} // namespace sd

namespace std
{
template <class T, size_t Rank>
struct tuple_size<sd::array<T, Rank>> : public integral_constant<size_t, Rank>
{
};

template <size_t Idx, class T, size_t Rank>
struct tuple_element<Idx, sd::array<T, Rank>>
{
    static_assert(Idx < Rank, "Index out of bounds in std::tuple_element<> (sd::array)");

public:
    using type = T;
};
#endif
} // namespace std

#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_LAYOUT_HXX
#define MARRAY_LAYOUT_HXX

namespace sd
{

template <size_t Rank>
class layout
{
public:
    constexpr static size_t rank = Rank;
    using index_type             = array<size_t, rank>;
    using size_type              = typename index_type::value_type;
    using sindex_type            = array<std::ptrdiff_t, rank>;

    enum class initial_ordering : uint8_t
    {
        RowMajor,
        ColMajor
    };

    constexpr layout() noexcept              = default;
    constexpr layout(const layout&) noexcept = default;
    constexpr layout& operator=(const layout&) noexcept = default;
    constexpr layout(layout&&) noexcept                 = default;
    constexpr layout& operator=(layout&&) noexcept = default;

    constexpr layout(const index_type& bounds, const index_type& strides)
        : m_extents(bounds), m_strides(strides)
    {
    }

    constexpr layout(const index_type& bounds, initial_ordering t = initial_ordering::RowMajor)
        : layout(bounds,
                 t == initial_ordering::RowMajor ? make_last_major_strides(bounds)
                                                 : make_first_major_strides(bounds))
    {
    }

    constexpr bool contains(const index_type& i) const noexcept { return i < m_extents; }

    template <size_t index, typename = std::enable_if_t<(index < Rank)>>
    constexpr auto extent() const noexcept
    {
        return m_extents[index];
    }

    template <size_t index, typename = std::enable_if_t<(index < Rank)>>
    constexpr auto strideof() const noexcept
    {
        return m_strides[index];
    }

    constexpr auto        extent(size_t index) const noexcept { return m_extents[index]; }
    constexpr auto        strideof(size_t index) const noexcept { return m_strides[index]; }
    constexpr const auto& extents() const noexcept { return m_extents; }
    constexpr const auto& strides() const noexcept { return m_strides; }

    constexpr auto& extent(size_t index) noexcept { return m_extents[index]; }
    constexpr auto& strideof(size_t index) noexcept { return m_strides[index]; }
    constexpr auto& extents() noexcept { return m_extents; }
    constexpr auto& strides() noexcept { return m_strides; }

    constexpr size_t size() const noexcept
    {
        size_type n = 0;
        for (size_t i = 0; i < rank; ++i)
        {
            if (m_extents[i] == 0) return 0;
            n += (m_extents[i] - 1) * m_strides[i];
        }
        return n + 1;
    }

    constexpr size_t total_size() const noexcept
    {
        size_type n = 1;
        for (size_t i = 0; i < rank; ++i) { n *= m_extents[i]; }
        return n;
    }

    constexpr size_t linearize(const index_type& ci) const noexcept
    {
        size_t offset = 0;
        for (size_t i = 0; i < rank; ++i) offset += m_strides[i] * ci[i];
        return offset;
    }

    constexpr ptrdiff_t linearize(const sindex_type& ci) const noexcept
    {
        ptrdiff_t offset = 0;
        for (size_t i = 0; i < rank; ++i)
            offset += m_strides[i] * (ci[i] < 0 ? m_extents[i] + ci[i] : ci[i]);
        return offset;
    }

    template <typename... Is>
    constexpr auto linearize(const Is&... is) const noexcept
    {
        return linearize(array<std::common_type_t<Is...>, rank>{is...});
    }

    constexpr index_type unlinearize(size_t index) const noexcept
    {
        if (index >= total_size()) return extents();

        index_type id{};
        for (size_t i = 0; i < rank; ++i)
        {
            id[i] = index / m_strides[i];
            index = index % m_strides[i];
        }
        return id;
    }

    constexpr void stride(const index_type& steps) noexcept
    {
        for (size_t i = 0; i < rank; ++i)
        {
            m_strides[i] *= steps[i];
            m_extents[i] =
                steps[i] == 0 ? 0 : std::ceil(static_cast<double>(m_extents[i]) / steps[i]);
        }
    }

    constexpr auto stride(const index_type& steps) const noexcept
    {
        auto copy = *this;
        copy.stride(steps);
        return copy;
    }

    static constexpr index_type make_last_major_strides(const index_type& ci) noexcept
    {
        index_type id{};
        id[Rank - 1] = 1;
        for (size_t i = rank - 1; i > 0; --i) id[i - 1] = id[i] * ci[i];
        return id;
    }

    static constexpr index_type make_first_major_strides(const index_type& ci) noexcept
    {
        index_type id{};
        id[0] = 1;
        for (size_t i = 1; i < rank; ++i) id[i] = id[i - 1] * ci[i - 1];
        return id;
    }

    constexpr bool operator==(const layout<rank>& r) const
    {
        return extents() == r.extents() && strides() == r.strides();
    }

    constexpr bool operator!=(const layout<rank>& r) const { return !(*this == r); }

    void transpose()
    {
        for (std::size_t j = 0; j < rank / 2; ++j)
        {
            std::size_t k = rank - j - 1;
            // transpose shape
            std::swap(extent(j), extent(k));
            // transpose strides
            std::swap(strideof(j), strideof(k));
        }
    }

    auto transposed() const
    {
        auto geo = *this;
        geo.transpose();
        return geo;
    }

private:
    index_type m_extents, m_strides;
};

#ifdef __cpp_deduction_guides

template <typename T, size_t R>
layout(const array<T, R>& a, ...) -> layout<R>;

template <typename T, size_t R>
layout(const std::array<T, R>& a, ...) -> layout<R>;

// template <typename T0, typename... T>
// layout(T0&&, T&&..., ...)->layout<1 + sizeof...(T)>;

// template <typename T>
// layout(std::initializer_list<T> a, ...)->layout<a.size()>;

#endif

} // namespace sd

#endif

#ifndef MARRAY_INDEX_ITERATOR_HXX
#define MARRAY_INDEX_ITERATOR_HXX

namespace sd
{

template <size_t Rank, size_t... is>
class index_iterator
{
public:
    static const size_t rank = Rank;
    using value_type         = array<size_t, Rank>;
    using reference          = value_type&;
    using pointer            = value_type*;
    using difference_type    = std::ptrdiff_t;
    using iterator_category  = std::random_access_iterator_tag;
    using const_pointer      = value_type const*;
    using const_reference    = value_type const&;
    using index_type         = value_type;
    using index_size_type    = typename index_type::value_type;

    constexpr index_iterator(const value_type& extents, value_type pos) noexcept
        : ext(extents), pos(std::move(pos))
    {
    }

    constexpr index_iterator(const value_type& extents, bool first) noexcept
        : index_iterator(extents, first ? index_type{} : extents)
    {
    }

    constexpr reference operator*() noexcept { return pos; }

    constexpr pointer operator->() noexcept { return &pos; }

    constexpr const_reference operator*() const noexcept { return pos; }

    constexpr const_pointer operator->() const noexcept { return &pos; }

    constexpr index_iterator& operator++() noexcept
    {
        for (std::size_t i = rank; i-- > 0;)
        {
            if (pos[i] + 1 < ext[i])
            {
                pos[i]++;
                return *this;
            }
            pos[i] = 0;
        }
        // If we're here we've wrapped over - set to past-the-end.
        pos = ext;
        return *this;

        // for (const auto i : order) {
        //     if (pos[i] + 1 < extent(i)) {
        //         ++pos[i];
        //         return *this;
        //     }
        //     pos[i] = 0;
        // }
        // pos = extents();
        // return *this;
    }

    constexpr index_iterator operator++(int) noexcept
    {
        auto ret = *this;
        ++(*this);
        return ret;
    }

    constexpr index_iterator& operator--() noexcept
    {

        if (pos >= ext)
        {
            // if at the past-the-end, set to last element
            for (std::size_t i = 0; i < rank; ++i) { pos[i] = ext[i] - 1; }
            return *this;
        }
        for (std::size_t i = rank; i-- > 0;)
        {
            if (pos[i] >= 1)
            {
                pos[i]--;
                return *this;
            }
            pos[i] = ext[i] - 1;
        }
        sd_expect(false);
        return *this;

        // if (pos >= extents()) {
        //     pos = extents() - index_type().fill(1);
        //     return *this;
        // }
        // sd_expect(pos < extents());
        // for (const auto i : order) {
        //     // for (size_t i = rank; i-- > 0;) {
        //     if (pos[i] >= 1) {
        //         --pos[i];
        //         return *this;
        //     }
        //     pos[i] = extent(i) - 1;
        // }
        // sd_expect(false);
        // return *this;
    }

    constexpr index_iterator operator--(int) noexcept
    {
        auto ret = *this;
        --(*this);
        return ret;
    }

    constexpr index_iterator operator+(difference_type n) const noexcept
    {
        index_iterator ret{*this};
        return ret += n;
    }

    constexpr index_iterator& operator+=(difference_type n) noexcept
    {
        if (n < 0)
            pos -= layout<rank>{extents()}.unlinearize(-n);
        else if (n > 0)
            pos += layout<rank>{extents()}.unlinearize(n);
        return *this;
    }

    constexpr index_iterator operator-(difference_type n) const noexcept
    {
        index_iterator ret{*this};
        return ret -= n;
    }

    constexpr index_iterator& operator-=(difference_type n) noexcept { return *this += -n; }

    constexpr difference_type operator-(const index_iterator& rhs) const noexcept
    {
        sd_expect(extents() == rhs.extents());
        const auto p = pos - rhs.pos;
        size_t     n = 1;
        for (size_t i = 0; i < Rank; ++i) n *= p[i];
        return n;
    }

    constexpr value_type operator[](difference_type n) const noexcept { return *(*this + n); }

    constexpr bool operator==(const index_iterator& rhs) const noexcept
    {
        return pos == rhs.pos;
    }

    constexpr bool operator!=(const index_iterator& rhs) const noexcept
    {
        return !(*this == rhs);
    }

    constexpr bool operator<(const index_iterator& rhs) const noexcept { return pos < rhs.pos; }

    constexpr bool operator<=(const index_iterator& rhs) const noexcept
    {
        return !(rhs < *this);
    }

    constexpr bool operator>(const index_iterator& rhs) const noexcept { return rhs < *this; }

    constexpr bool operator>=(const index_iterator& rhs) const noexcept
    {
        return !(rhs > *this);
    }

    void swap(index_iterator& rhs) noexcept
    {
        sd_expect(rhs.ext == ext);
        std::swap(pos, rhs.pos);
    }

private:
    const value_type& extents() const { return ext; }
    index_size_type   extent(index_size_type index) const { return (ext)[index]; }

    static_assert(sizeof...(is) == 0 || sizeof...(is) == rank,
                  "invalid rank of custom iteration order (index_iteratior)");

    template <size_t R>
    static constexpr array<size_t, R> make_default_order()
    {
        array<size_t, R> arr;
        for (size_t i = 0; i < R; ++i) arr[i] = (R - 1) - i;
        return arr;
    }

    static constexpr const auto order =
        sizeof...(is) == 0 ? make_default_order<rank>() : array<size_t, rank>{is...};

    value_type ext;
    value_type pos;
};

template <size_t R, size_t... is>
struct index_range
{
    array<size_t, R>         m_extents;
    index_iterator<R, is...> begin() const { return {m_extents, total_size() > 0}; }
    index_iterator<R, is...> end() const { return {m_extents, false}; }
    constexpr size_t         total_size() const noexcept
    {
        size_t n = 1;
        for (auto const& i : m_extents) n *= i;
        return n;
    }
};

namespace internal
{
template <size_t i, size_t... is>
struct forward_iteration
{
    using type = typename forward_iteration<i - 1, i - 1, is...>::type;
};

template <size_t... is>
struct forward_iteration<0, is...>
{
    using type = index_range<sizeof...(is), is...>;
};

template <size_t i, size_t... is>
struct backward_iteration
{
    using type = typename backward_iteration<i - 1, is..., i - 1>::type;
};

template <size_t... is>
struct backward_iteration<0, is...>
{
    using type = index_range<sizeof...(is), is...>;
};
} // namespace internal

template <size_t r>
using last_major_range = typename internal::forward_iteration<r>::type;

template <size_t r>
using first_major_range = typename internal::backward_iteration<r>::type;
} // namespace sd

#endif
/// See LICENSE" for the full license governing this code.

#ifndef MARRAY_SLICE_HXX
#define MARRAY_SLICE_HXX

namespace sd
{

template <typename V,
          typename = decltype(std::declval<V>().data()),
          typename = decltype(std::declval<V>().size())>
struct is_viewable
{
    using type              = V;
    static const bool value = true;
    constexpr         operator bool() const { return value; }
    constexpr bool    operator()() const { return value; }
};

// template <typename U, typename T>
// struct is_safely_convertable_value
//     : public std::integral_constant<
//           bool,
//           std::is_convertible<std::add_pointer_t<U>, std::add_pointer_t<T>>::value>
// {
// };

template <typename... Ts>
class are_integral : public std::integral_constant<bool, true>
{
};

template <typename T, typename... Ts>
class are_integral<T, Ts...>
    : public std::integral_constant<bool,
                                    std::is_integral<T>::value && are_integral<Ts...>::value>
{
};

template <typename T>
struct pointer_extractor_type
{
    constexpr pointer_extractor_type(std::shared_ptr<T>& p) noexcept : raw{p.get()} {}
    constexpr pointer_extractor_type(std::shared_ptr<T[]>& p) noexcept : raw{p.get()} {}
    constexpr pointer_extractor_type(std::unique_ptr<T>& p) noexcept : raw{p.get()} {}
    constexpr pointer_extractor_type(std::unique_ptr<T[]>& p) noexcept : raw{p.get()} {}
    constexpr pointer_extractor_type(std::nullptr_t) noexcept : raw{nullptr} {}
    constexpr pointer_extractor_type(T* p) noexcept : raw{p} {}
    template <size_t S>
    constexpr pointer_extractor_type(T p[S]) noexcept : raw{p}
    {
    }
    template <typename X, typename = decltype(std::declval<X>().data())>
    constexpr pointer_extractor_type(X& x) noexcept : raw{x.data()}
    {
        static_assert(std::is_convertible<decltype(x.data()), std::add_pointer_t<T>>::value,
                      "illegal viewable. (slice.pointer_extractor_type) ");
    }

    constexpr T* get() { return raw; }
    constexpr    operator T*() { return raw; }

private:
    T* raw{nullptr};
};

template <class T, size_t Rank>
class slice;

template <typename T, size_t Rank>
class cpslice;

template <typename T, size_t... is>
class basic_slice_iterator;

template <typename T>
struct data_span
{
    using pointer = T*;

    pointer begin() const { return f; }
    pointer end() const { return l; }
    size_t  size() const { return l - f; }
    pointer data() const { return f; }

    pointer f, l;
};

template <typename T, size_t rank, size_t... is>
struct ordered_range
{
    auto begin() const { return f; }
    auto end() const { return l; }

    basic_slice_iterator<slice<T, rank> const&, is...> f, l;
};

template <class T, size_t Rank = 1>
class slice
{
public:
    static constexpr size_t rank = Rank;

    using layout_type = layout<rank>;
    using index_type  = typename layout_type::index_type;
    using size_type   = typename index_type::value_type;
    using sindex_type = typename layout_type::sindex_type;
    using value_type  = T;
    using pointer     = T*;
    using reference   = T&;
    using iterator    = basic_slice_iterator<slice<T, rank> const&>;
    using extract_ptr = pointer_extractor_type<T>;

    constexpr slice() noexcept             = default;
    constexpr slice(const slice&) noexcept = default;
    constexpr slice(slice&&) noexcept      = default;
    constexpr slice& operator=(const slice&) noexcept = default;
    constexpr slice& operator=(slice&&) noexcept = default;
    ~slice() noexcept                            = default;

    constexpr slice(extract_ptr ptr, layout_type geo) noexcept
        : elems(ptr), geometry(std::move(geo))
    {
    }
    constexpr slice(extract_ptr ptr, index_type bounds) noexcept
        : elems(ptr), geometry(std::move(bounds))
    {
    }
    constexpr slice(extract_ptr ptr, index_type bounds, index_type strides) noexcept
        : elems(ptr), geometry(std::move(bounds), std::move(strides))
    {
    }
    // template <typename U>
    // constexpr slice(slice<U, rank> rhs) noexcept : slice(rhs, rhs.shape())
    // {
    // }
    template <typename... Is,
              typename = std::enable_if_t<sizeof...(Is) == rank && are_integral<Is...>::value>>
    constexpr slice(extract_ptr v, Is&&... is) noexcept
        : elems(v), geometry(index_type{std::forward<Is>(is)...})
    {
    }
    template <bool cond = (rank == 1), typename = std::enable_if_t<cond>>
    constexpr slice(extract_ptr a, extract_ptr b)
        : elems(a), geometry(index_type{static_cast<size_t>(b.get() - a.get())})
    {
    }
    template <typename Viewable,
              bool cond = (rank == 1),
              typename  = std::enable_if_t<cond>,
              typename  = is_viewable<Viewable>>
    constexpr slice(Viewable& vw) noexcept : slice(vw, layout_type{{vw.size()}})
    {
    }

    //  constexpr slice(pointer_extractor_type<T> ptr, index_type bounds) noexcept
    //      : elems(ptr), geometry(std::move(bounds))
    //  {
    //  }

    // template<typename... Ts>
    // constexpr slice(pointer_extractor_type<T> ptr, Ts&&... ts) noexcept
    //     : elems(ptr), geometry(std::forward<Ts>(ts)...)
    // {
    // }

    // constexpr slice(pointer_extractor_type<T> ptr, layout_type geo) noexcept
    //     : elems(ptr), geometry(std::move(geo))
    // {
    // }
    // constexpr slice(pointer ptr, index_type bounds) noexcept
    //     : elems(ptr), geometry(std::move(bounds))
    // {
    // }
    // constexpr slice(pointer ptr, index_type bounds, index_type strides) noexcept
    //     : elems(ptr), geometry(std::move(bounds), std::move(strides))
    // {
    // }
    // constexpr slice(pointer ptr, layout_type geo) noexcept
    //     : elems(ptr), geometry(std::move(geo))
    // {
    // }
    // constexpr slice(pointer_extractor_type<T> ptr, layout_type geo) noexcept
    //     : elems(ptr), geometry(std::move(geo))
    // {
    // }
    // template <bool cond = (rank == 1), typename = std::enable_if_t<cond>>
    // constexpr slice(pointer a, pointer b)
    //     : elems(b), geometry(index_type{static_cast<size_t>(b - a)})
    // {
    // }
    // template <size_t Extent, typename = std::enable_if_t<(Extent > 0 && rank == 1)>>
    // constexpr slice(value_type arr[Extent]) noexcept : elems(arr), geometry({Extent}, {1})
    // {
    // }
    // template <class U,
    //           typename = std::enable_if_t<is_safely_convertable_value<U, T>::value &&
    //                                       !std::is_same<U, T>::value>>
    // constexpr slice(const slice<U, rank>& rhs) noexcept
    //     : elems(rhs.data()), geometry(rhs.shape())
    // {
    // }
    // template <
    //     class U,
    //     size_t S,
    //     typename = std::enable_if_t<(rank == 1) && is_safely_convertable_value<U, T>::value>>
    // constexpr slice(const slice<U, S>& rhs) noexcept : slice(rhs.data(), rhs.size())
    // {
    // }
    // template <class Viewable, typename = is_viewable_t<Viewable>>
    // constexpr slice(Viewable& vw, layout_type bounds) noexcept
    //     : elems(vw.data()), geometry(std::move(bounds))
    // {
    //     sd_expect(geometry.total_size() <= vw.size());
    //     static_assert(std::is_convertible<decltype(vw.size()), size_type>::value,
    //                   "illegal viewable (slice)");
    //     using U = std::remove_pointer_t<decltype(vw.data())>;
    //     static_assert(is_safely_convertable_value<U, T>::value, "illegal viewable (slice)");
    // }

    // template <class Viewable, typename = is_viewable_t<Viewable>>
    // constexpr slice(Viewable& vw, index_type bounds) noexcept
    //     : slice(vw, layout_type{std::move(bounds)})
    // {
    // }
    // template <class Viewable, typename = is_viewable_t<Viewable>>
    // constexpr slice(Viewable& vw, index_type bounds, index_type strides) noexcept
    //     : slice(vw, layout_type{std::move(bounds), std::move(strides)})
    // {
    // }
    // template <typename... Is,
    //           typename = std::enable_if_t<sizeof...(Is) == rank &&
    //           are_integral<Is...>::value>>
    // constexpr slice(pointer ptr, const Is&... is) noexcept : slice(ptr, index_type{is...})
    // {
    // }

    // template <typename Viewable,
    //           typename... Is,
    //           typename = is_viewable_t<Viewable>,
    //           typename = std::enable_if_t<sizeof...(Is) == rank &&
    //           are_integral<Is...>::value>>
    // constexpr slice(Viewable& v, const Is&... is) noexcept
    //     : slice(v, index_type{std::forward<Is...>(is)...})
    // {
    // }

    template <size_t... is>
    constexpr index_range<rank, is...> subscripts() const noexcept
    {
        return {extents()};
    }
    constexpr last_major_range<rank> last_major_subscripts() const noexcept
    {
        return {extents()};
    }
    constexpr index_range<1> indices() const noexcept { return {size()}; }
    template <size_t... is>
    constexpr ordered_range<T, rank, is...> iter() const noexcept
    {
        return {sbegin(), send()};
    }
    constexpr data_span<T> span() const noexcept { return {data(), data() + size()}; }

    constexpr iterator sbegin() const noexcept { return {this, size() > 0}; }
    constexpr iterator send() const noexcept { return {this, false}; }
    constexpr pointer  begin() const noexcept { return data(); }
    constexpr pointer  end() const noexcept { return data() + size(); }

    constexpr std::reverse_iterator<iterator> rsbegin() const noexcept
    {
        return std::make_reverse_iterator(send());
    }
    constexpr std::reverse_iterator<iterator> rsend() const noexcept
    {
        return std::make_reverse_iterator(sbegin());
    }
    constexpr std::reverse_iterator<pointer> rbegin() const noexcept
    {
        return std::make_reverse_iterator(end());
    }
    constexpr std::reverse_iterator<pointer> rend() const noexcept
    {
        return std::make_reverse_iterator(begin());
    }

    explicit constexpr  operator bool() const noexcept { return !empty(); }
    constexpr bool      empty() const noexcept { return elems == nullptr || size() == 0; }
    constexpr size_type size() const noexcept { return geometry.total_size(); }
    constexpr size_type length() const noexcept { return geometry.extent(0); }
    constexpr pointer   data() const noexcept { return elems; }
    constexpr const index_type&  strides() const noexcept { return geometry.strides(); }
    constexpr const index_type&  extents() const noexcept { return geometry.extents(); }
    constexpr const layout_type& shape() const noexcept { return geometry; }

    template <size_t index, typename = std::enable_if_t<(index < Rank)>>
    constexpr size_type extent() const noexcept
    {
        return geometry.template extent<index>();
    }
    template <size_t index, typename = std::enable_if_t<(index < Rank)>>
    constexpr size_type strideof() const noexcept
    {
        return geometry.template strideof<index>();
    }
    constexpr size_type extent(size_type index) const noexcept
    {
        return geometry.extent(index);
    }
    constexpr size_type strideof(size_type index) const noexcept
    {
        return geometry.strideof(index);
    }

    constexpr reference operator[](const index_type& idx) const noexcept
    {
        sd_expect(geometry.contains(idx));
        return elems[geometry.linearize(idx)];
    }

    template <bool cond = (rank == 1), typename = std::enable_if_t<cond>>
    constexpr reference operator[](const size_type idx) const noexcept
    {
        sd_expect(idx < size());
        return operator[](index_type{idx});
    }

    template <bool cond = (rank > 1), typename = std::enable_if_t<cond>>
    constexpr slice<T, rank - 1> operator[](const size_type i) const noexcept
    {
        sd_expect(i < extent(0));

        const index_type& e = geometry.extents();
        const index_type& s = geometry.strides();

        return {elems + s[0] * i, {e.begin() + 1, e.end()}, {s.begin() + 1, s.end()}};
    }

    constexpr reference operator()(const size_type idx) const noexcept
    {
        // !!direct access, not strided!!
        sd_expect(idx < size());
        return elems[idx];
    }

    constexpr reference operator()(const index_type& idx) const noexcept
    {
        return operator[](idx);
    }

    template <typename... Is,
              typename = std::enable_if_t<are_integral<Is...>::value && sizeof...(Is) == rank>>
    constexpr reference operator()(const Is&... is) const noexcept
    {
        return operator[](index_type{is...});
    }

    constexpr slice<T, rank> cut_extent(const index_type& origin,
                                        const index_type& slice_extents,
                                        const index_type& steps) const noexcept
    {
        return cut_extent(origin, slice_extents).stride(steps);
    }

    constexpr slice<T, rank> cut_extent(const index_type& origin,
                                        const index_type& slice_extents) const noexcept
    {
        sd_expect(extents() - origin >= slice_extents);
        return {elems + geometry.linearize(origin), slice_extents, strides()};
    }

    constexpr slice<T, rank> cut(const index_type& origin,
                                 const index_type& bound,
                                 const index_type& steps) const noexcept
    {
        return cut(origin, bound).stride(steps);
    }

    constexpr slice<T, rank> cut(const index_type& origin,
                                 const index_type& bound) const noexcept
    {
        sd_expect(bound <= extents());
        sd_expect(origin <= bound);
        return cut_extent(origin, bound - origin);
    }

    constexpr slice<T, rank> cut(const index_type& origin) const noexcept
    {
        return cut_extent(origin, extents() - origin);
    }

    constexpr slice<T, rank> subspan(const sindex_type& origin,
                                     const sindex_type& bound,
                                     const index_type&  steps) const noexcept
    {
        return cut(subspan_index_to_cut_index(origin), subspan_index_to_cut_index(bound))
            .stride(steps);
    }

    constexpr slice<T, rank> subspan(const sindex_type& origin,
                                     const sindex_type& bound) const noexcept
    {
        // return {elems + geometry.linearize(origin), subspan_index_to_cut_index(bound)};
        return cut(subspan_index_to_cut_index(origin), subspan_index_to_cut_index(bound));
        //- subspan_index_to_cut_index(origin)
    }

    constexpr slice<T, rank> subspan(const sindex_type& origin) const noexcept
    {
        return cut(subspan_index_to_cut_index(origin));
    }

    constexpr slice<T, rank> stride(const index_type& steps) const noexcept
    {
        return {elems, geometry.stride(steps)};
    }

    constexpr slice<T, 1> flatten() const noexcept { return {data(), size()}; }

    template <typename... Is, typename = std::enable_if_t<are_integral<Is...>::value>>
    constexpr auto as(Is const&... s) const noexcept
    {
        return as(layout<sizeof...(Is)>{{s...}});
    }
    template <size_t S>
    constexpr auto as(const array<size_t, S>& other_extents) const noexcept
    {
        return as(layout<S>{other_extents});
    }
    template <size_t S>
    constexpr auto as(const array<size_t, S>& other_extents,
                      const array<size_t, S>& other_strides) const noexcept
    {
        return as(layout<S>{other_extents, other_strides});
    }
    template <size_t S>
    constexpr auto as(const layout<S>& l) const noexcept
    {
        sd_expect(l.total_size() <= geometry.total_size());
        return slice<T, S>(data(), l);
    }

    constexpr slice<T, rank> transpose() const noexcept
    {
        return {elems, geometry.transposed()};
    }
    constexpr cpslice<T, rank> propagate_const() const noexcept { return {*this}; }

    constexpr decltype(auto) front() const noexcept { return operator()(0); }
    constexpr decltype(auto) back() const noexcept { return operator()(size() - 1); }
    constexpr reference      front_element() const noexcept { return operator[](index_type{}); }
    constexpr reference      back_element() const noexcept
    {
        return operator[](extents() - index_type{}.fill(1));
    }

private:
    index_type subspan_index_to_cut_index(const sindex_type& s) const noexcept
    {
        index_type t;
        for (size_t i = 0; i < rank; ++i) { t[i] = s[i] < 0 ? extent(i) + s[i] : s[i]; }
        return t;
    }

    pointer     elems{nullptr};
    layout_type geometry;
    template <typename T2, size_t R2>
    friend class slice;
};

template <typename T, size_t r, typename S, size_t s>
constexpr bool operator==(const slice<T, r>& a, const slice<S, s>& b) noexcept
{
    return a.size() == b.size() && std::equal(a.begin(), a.end(), b.begin(), b.end());
    // return a.extents() == b.extents() && a.strides() == b.strides() && a.data() == b.data();
}

template <typename T, size_t r, typename S, size_t s>
constexpr bool operator!=(const slice<T, r>& a, const slice<S, s>& b) noexcept
{
    return !(a == b);
}

template <typename T, size_t r>
constexpr decltype(auto) begin(const sd::slice<T, r>& s)
{
    return s.begin();
}

template <typename T, size_t r>
constexpr decltype(auto) end(const sd::slice<T, r>& s)
{
    return s.end();
}

template <typename T, size_t r>
constexpr decltype(auto) data(const sd::slice<T, r>& s)
{
    return s.data();
}

template <typename T, size_t r>
constexpr decltype(auto) size(const sd::slice<T, r>& s)
{
    return s.size();
}

template <typename T, size_t r>
constexpr decltype(auto) empty(const sd::slice<T, r>& s)
{
    return s.empty();
}

template <typename V, size_t... is>
class basic_slice_iterator
{
public:
    using slice_type = std::decay_t<V>;

    using iterator_category = std::random_access_iterator_tag;
    using value_type        = typename slice_type::value_type;
    using reference         = value_type&;
    using pointer           = value_type*;
    using difference_type   = std::ptrdiff_t;
    using const_pointer     = value_type const*;
    using const_reference   = value_type const&;

    using index_type      = typename slice_type::index_type;
    using index_size_type = typename index_type::value_type;
    using iterator        = index_iterator<slice_type::rank, is...>;

    constexpr basic_slice_iterator(const slice_type* view, bool first)
        : m_view(view)
        , m_iter(first ? m_view->template subscripts<is...>().begin()
                       : m_view->template subscripts<is...>().end())
    {
    }

    constexpr basic_slice_iterator(const slice_type*                       view,
                                   index_iterator<slice_type::rank, is...> iter)
        : m_view(view), m_iter(iter)
    {
    }

    template <typename = decltype(std::declval<slice_type&>().operator[](index_type{}))>
    constexpr decltype(auto) operator*() noexcept
    {
        return (*m_view)[*m_iter];
    }
    constexpr decltype(auto) operator*() const noexcept { return (*m_view)[*m_iter]; }
    constexpr auto           operator->() noexcept { return &(*m_view)[*m_iter]; }

    constexpr basic_slice_iterator& operator++() noexcept
    {
        ++m_iter;
        return *this;
    }

    constexpr basic_slice_iterator operator++(int) noexcept
    {
        auto ret = *this;
        ++(*this);
        return ret;
    }

    constexpr basic_slice_iterator& operator--() noexcept
    {
        --m_iter;
        return *this;
    }

    constexpr basic_slice_iterator operator--(int) noexcept
    {
        auto ret = *this;
        --(*this);
        return ret;
    }

    constexpr basic_slice_iterator operator+(difference_type n) const noexcept
    {
        basic_slice_iterator ret{*this};
        return ret += n;
    }

    constexpr basic_slice_iterator& operator+=(difference_type n) noexcept
    {
        m_iter += n;
        return *this;
    }

    constexpr basic_slice_iterator operator-(difference_type n) const noexcept
    {
        basic_slice_iterator ret{*this};
        return ret -= n;
    }

    constexpr basic_slice_iterator& operator-=(difference_type n) noexcept
    {
        return *this += -n;
    }

    constexpr difference_type operator-(const basic_slice_iterator& rhs) const noexcept
    {
        sd_expect(m_view == rhs.m_view);
        return m_iter - rhs.m_iter;
    }

    constexpr value_type operator[](difference_type n) const noexcept
    {
        return (*m_view)[m_iter[n]];
    }

    constexpr bool operator==(const basic_slice_iterator& rhs) const noexcept
    {
        sd_expect(m_view == rhs.m_view);
        return m_iter == rhs.m_iter;
    }

    constexpr bool operator!=(const basic_slice_iterator& rhs) const noexcept
    {
        return !(*this == rhs);
    }

    constexpr bool operator<(const basic_slice_iterator& rhs) const noexcept
    {
        sd_expect(m_view == rhs.m_view);
        return m_iter < rhs.m_iter;
    }

    constexpr bool operator<=(const basic_slice_iterator& rhs) const noexcept
    {
        return !(rhs < *this);
    }

    constexpr bool operator>(const basic_slice_iterator& rhs) const noexcept
    {
        return rhs < *this;
    }

    constexpr bool operator>=(const basic_slice_iterator& rhs) const noexcept
    {
        return !(rhs > *this);
    }

    void swap(basic_slice_iterator& rhs) noexcept
    {
        std::swap(m_iter, rhs.m_iter);
        std::swap(m_view, rhs.m_view);
    }

    auto index() const { return *m_iter; }

    constexpr friend basic_slice_iterator operator+(difference_type             n,
                                                    const basic_slice_iterator& rhs)
    {
        return rhs + n;
    }

private:
    template <typename U, size_t Rank>
    friend class slice;

    slice_type const* m_view;
    iterator          m_iter;
};
} // namespace sd

#ifdef __cpp_deduction_guides
namespace sd
{
template <typename V, typename = std::enable_if_t<is_viewable<V>::value>>
slice(V&& v) -> slice<std::remove_reference_t<decltype(*v.data())>, 1>;

template <typename V, size_t R, typename = std::enable_if_t<is_viewable<V>::value>>
slice(V&& v, array<size_t, R> i, ...) -> slice<std::remove_reference_t<decltype(*v.data())>, R>;

template <typename V, size_t R>
slice(V* v, array<size_t, R> i, ...) -> slice<std::remove_reference_t<V>, R>;

template <typename V, size_t R, typename = std::enable_if_t<is_viewable<V>::value>>
slice(V&& v, layout<R> geo) -> slice<std::remove_reference_t<decltype(*v.data())>, R>;

template <typename V, size_t R>
slice(V* v, layout<R> geo) -> slice<std::remove_reference_t<V>, R>;

template <typename V,
          typename... Is,
          typename = std::enable_if_t<are_integral<Is...>::value && sizeof...(Is) != 0>>
slice(V* ptr, Is&&... is) -> slice<std::remove_reference_t<V>, sizeof...(Is)>;

template <typename V,
          typename... Is,
          typename = std::enable_if_t<are_integral<Is...>::value && sizeof...(Is) != 0>>
slice(V&& v, Is&&... is) -> slice<std::remove_reference_t<decltype(*v.data())>, sizeof...(Is)>;
} // namespace sd
#endif
#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_CONST_PROPAGATING_SLICE_HXX
#define MARRAY_CONST_PROPAGATING_SLICE_HXX

namespace sd
{
template <typename T, size_t Rank = 1>
class cpslice
{
public:
    constexpr static const size_t rank = Rank;

    using view_type       = slice<T, rank>;
    using const_view_type = slice<std::add_const_t<T>, rank>;
    using index_type      = typename view_type::index_type;
    using size_type       = typename view_type::size_type;
    using value_type      = typename view_type::value_type;
    using layout_type     = typename view_type::layout_type;
    using iterator        = typename view_type::iterator;
    using pointer         = typename view_type::pointer;
    using sindex_type     = typename view_type::sindex_type;

    union storage_type
    {
        view_type       view{};
        const_view_type view_const;
    } storage{};

    constexpr cpslice(view_type x) noexcept { storage.view = x; }

    template <typename Arg0, typename... Args>
    constexpr cpslice(Arg0&& arg0, Args&&... args) noexcept
    {
        static_assert(std::is_constructible<view_type, Arg0, Args...>::value,
                      "cannot construct cpslice");
        storage.view = view_type{std::forward<Arg0>(arg0), std::forward<Args>(args)...};
    }

    constexpr cpslice() noexcept               = default;
    ~cpslice() noexcept                        = default;
    constexpr cpslice(cpslice&&) noexcept      = default;
    constexpr cpslice(const cpslice&) noexcept = default;
    constexpr cpslice& operator=(cpslice&&) noexcept = default;
    constexpr cpslice& operator=(const cpslice&) noexcept = default;

    //////////////////////////////////////////////////
    // propagate_const, I want to have operator.
    //////////////////////////////////////////////////
    // clang-format off
    constexpr const const_view_type& get() const { return storage.view_const; }
    constexpr             view_type& get()       { return storage.view; }

    constexpr auto operator*() const  -> decltype(auto) { return get(); }
    constexpr auto operator*()        -> decltype(auto) { return get(); }
    constexpr auto operator->() const -> decltype(auto) { return &get(); }
    constexpr auto operator->()       -> decltype(auto) { return &get(); }

    constexpr operator view_type()             { return {get()}; }
    constexpr operator const_view_type() const { return {get()}; }
    //////////////////////////////////////////////////
    // slice interface                              
    //////////////////////////////////////////////////
                            constexpr auto last_major_subscripts() const noexcept -> decltype(auto) { return get().last_major_subscripts(); }     
    template <size_t... is> constexpr auto subscripts()            const noexcept -> decltype(auto) { return get().template subscripts<is...>(); }
    template <size_t... is> constexpr auto iter()                  const noexcept -> decltype(auto) { return get().template iter<is...>(); }    
    template <size_t index> constexpr auto extent()                const noexcept -> decltype(auto) { return get().template extent<index>(); }
    template <size_t index> constexpr auto strideof()              const noexcept -> decltype(auto) { return get().template strideof<index>(); }

    constexpr explicit operator bool()       const noexcept                   { return get().operator bool(); }
    constexpr auto indices()                 const noexcept -> decltype(auto) { return get().indices(); }                       
    constexpr auto span()                    const noexcept -> decltype(auto) { return get().span(); }
    constexpr auto data()                    const noexcept -> decltype(auto) { return get().data(); }
    constexpr auto data()                          noexcept -> decltype(auto) { return get().data(); }
    constexpr auto size()                    const noexcept -> decltype(auto) { return get().size(); }
    constexpr auto length()                  const noexcept -> decltype(auto) { return get().length(); }
    constexpr auto extents()                 const noexcept -> decltype(auto) { return get().extents(); }
    constexpr auto strides()                 const noexcept -> decltype(auto) { return get().strides(); }    
    constexpr auto shape()                   const noexcept -> decltype(auto) { return get().shape(); }
    constexpr auto extent(size_type index)   const noexcept -> decltype(auto) { return get().extent(index); }
    constexpr auto strideof(size_type index) const noexcept -> decltype(auto) { return get().strideof(index); }
    constexpr auto empty()                   const noexcept -> decltype(auto) { return get().empty(); }
    constexpr auto begin()                   const noexcept -> decltype(auto) { return get().begin(); }
    constexpr auto begin()                         noexcept -> decltype(auto) { return get().begin(); }
    constexpr auto end()                     const noexcept -> decltype(auto) { return get().end(); }    
    constexpr auto end()                           noexcept -> decltype(auto) { return get().end(); }
    constexpr auto sbegin()                  const noexcept -> decltype(auto) { return get().sbegin(); }
    constexpr auto sbegin()                        noexcept -> decltype(auto) { return get().sbegin(); }
    constexpr auto send()                    const noexcept -> decltype(auto) { return get().send(); }    
    constexpr auto send()                          noexcept -> decltype(auto) { return get().send(); }
    constexpr auto rbegin()                  const noexcept -> decltype(auto) { return get().rbegin(); }
    constexpr auto rbegin()                        noexcept -> decltype(auto) { return get().rbegin(); }
    constexpr auto rend()                    const noexcept -> decltype(auto) { return get().rend(); }    
    constexpr auto rend()                          noexcept -> decltype(auto) { return get().rend(); }
    constexpr auto rsbegin()                 const noexcept -> decltype(auto) { return get().rsbegin(); }
    constexpr auto rsbegin()                       noexcept -> decltype(auto) { return get().rsbegin(); }
    constexpr auto rsend()                   const noexcept -> decltype(auto) { return get().rsend(); }    
    constexpr auto rsend()                         noexcept -> decltype(auto) { return get().rsend(); }
    constexpr auto flatten()                 const noexcept -> decltype(auto) { return get().flatten(); }
    constexpr auto flatten()                       noexcept -> decltype(auto) { return get().flatten(); }
    constexpr auto transpose()               const noexcept -> decltype(auto) { return get().transpose(); }    
    constexpr auto transpose()                     noexcept -> decltype(auto) { return get().transpose(); }
    constexpr auto propagate_const()         const noexcept -> decltype(auto) { return *this; }

                              constexpr auto operator() (const index_type& idx) const noexcept -> decltype(auto) { return get()(idx); }
                              constexpr auto operator() (const index_type& idx)       noexcept -> decltype(auto) { return get()(idx); }
    template <typename... Is> constexpr auto operator() (const Is&... is)       const noexcept -> decltype(auto) { return get()(is...); }
    template <typename... Is> constexpr auto operator() (const Is&... is)             noexcept -> decltype(auto) { return get()(is...); }

    constexpr auto operator[](const size_type& idx)  const noexcept -> decltype(auto) { return get()[idx]; }
    constexpr auto operator[](const size_type& idx)        noexcept -> decltype(auto) { return get()[idx]; }
    constexpr auto operator[](const index_type& idx) const noexcept -> decltype(auto) { return get()[idx]; }
    constexpr auto operator[](const index_type& idx)       noexcept -> decltype(auto) { return get()[idx]; }
    
    constexpr decltype(auto) front()         const noexcept { return get().front();         }
    constexpr decltype(auto) back()          const noexcept { return get().back();          }
    constexpr decltype(auto) front_element() const noexcept { return get().front_element(); }
    constexpr decltype(auto) back_element()  const noexcept { return get().back_element();  }
    constexpr decltype(auto) front()               noexcept { return get().front();         }
    constexpr decltype(auto) back()                noexcept { return get().back();          }
    constexpr decltype(auto) front_element()       noexcept { return get().front_element(); }
    constexpr decltype(auto) back_element()        noexcept { return get().back_element();  }

    constexpr decltype(auto) cut        (const index_type& origin, const index_type& slice_extents, const index_type& strides)   const noexcept { return get().cut(origin, slice_extents, strides); }
    constexpr decltype(auto) cut        (const index_type& origin, const index_type& slice_extents, const index_type& strides)         noexcept { return get().cut(origin, slice_extents, strides); }
    constexpr decltype(auto) cut        (const index_type& origin, const index_type& slice_extents)                              const noexcept { return get().cut(origin, slice_extents); }
    constexpr decltype(auto) cut        (const index_type& origin, const index_type& slice_extents)                                    noexcept { return get().cut(origin, slice_extents); }
    constexpr decltype(auto) cut        (const index_type& origin)                                                               const noexcept { return get().cut(origin); }
    constexpr decltype(auto) cut        (const index_type& origin)                                                                     noexcept { return get().cut(origin); }
    constexpr decltype(auto) cut_extent (const index_type& origin, const index_type& slice_extents, const index_type& strides)   const noexcept { return get().cut_extent(origin, slice_extents, strides); }
    constexpr decltype(auto) cut_extent (const index_type& origin, const index_type& slice_extents, const index_type& strides)         noexcept { return get().cut_extent(origin, slice_extents, strides); }
    constexpr decltype(auto) cut_extent (const index_type& origin, const index_type& slice_extents)                              const noexcept { return get().cut_extent(origin, slice_extents); }
    constexpr decltype(auto) cut_extent (const index_type& origin, const index_type& slice_extents)                                    noexcept { return get().cut_extent(origin, slice_extents); }
    constexpr decltype(auto) cut_extent (const index_type& origin)                                                               const noexcept { return get().cut_extent(origin); }
    constexpr decltype(auto) cut_extent (const index_type& origin)                                                                     noexcept { return get().cut_extent(origin); }
    constexpr decltype(auto) subspan    (const sindex_type& origin, const sindex_type& slice_extents, const index_type& strides) const noexcept { return get().subspan(origin, slice_extents, strides); }
    constexpr decltype(auto) subspan    (const sindex_type& origin, const sindex_type& slice_extents, const index_type& strides)       noexcept { return get().subspan(origin, slice_extents, strides); }
    constexpr decltype(auto) subspan    (const sindex_type& origin, const sindex_type& slice_extents)                            const noexcept { return get().subspan(origin, slice_extents); }
    constexpr decltype(auto) subspan    (const sindex_type& origin, const sindex_type& slice_extents)                                  noexcept { return get().subspan(origin, slice_extents); }
    constexpr decltype(auto) subspan    (const sindex_type& origin)                                                              const noexcept { return get().subspan(origin); }
    constexpr decltype(auto) subspan    (const sindex_type& origin)                                                                    noexcept { return get().subspan(origin); }

    // template<typename... IndexTypes>
    // constexpr decltype(auto) cut_extent(IndexTypes&&... is) const noexcept { return get().cut_extent(std::forward<IndexTypes>(is)...); }
    // template<typename... IndexTypes>
    // constexpr decltype(auto) cut_extent(IndexTypes&&... is)       noexcept { return get().cut_extent(std::forward<IndexTypes>(is)...); }

    // template<typename... IndexTypes>
    // constexpr decltype(auto) cut(IndexTypes&&... is) const noexcept { return get().cut(std::forward<IndexTypes>(is)...); }
    // template<typename... IndexTypes>
    // constexpr decltype(auto) cut(IndexTypes&&... is)       noexcept { return get().cut(std::forward<IndexTypes>(is)...); }

    // template<typename... IndexTypes>
    // constexpr decltype(auto) subspan(IndexTypes&&... is) const noexcept { return get().subspan(std::forward<IndexTypes>(is)...); }
    // template<typename... IndexTypes>
    // constexpr decltype(auto) subspan(IndexTypes&&... is)       noexcept { return get().subspan(std::forward<IndexTypes>(is)...); }

                                    constexpr auto stride    (const index_type& steps)  const noexcept -> decltype(auto) { return get().stride(steps); }
                                    constexpr auto stride    (const index_type& steps)        noexcept -> decltype(auto) { return get().stride(steps); }
    template <typename... Is>       constexpr auto as        (const Is&... s)           const noexcept -> decltype(auto) { return get().as(s...); }
    template <typename... Is>       constexpr auto as        (const Is&... s)                 noexcept -> decltype(auto) { return get().as(s...); }
    template <size_t r, typename S> constexpr auto operator!=(const slice<S, r>& b)     const noexcept -> bool           { return this->operator*() != b; }
    template <size_t r, typename S> constexpr auto operator==(const slice<S, r>& b)     const noexcept -> bool           { return this->operator*() == b; }
    // clang-format on
};

template <typename T, size_t r, typename S>
constexpr bool operator==(const cpslice<T, r>& a, const cpslice<S, r>& b) noexcept
{
    return a.extents() == b.extents() && a.strides() == b.strides() && a.data() == b.data();
}

template <typename T, size_t r, typename S, size_t s>
constexpr bool operator!=(const cpslice<T, r>& a, const cpslice<S, s>& b) noexcept
{
    return !(a == b);
}

template <typename T, size_t r>
constexpr decltype(auto) begin(const sd::cpslice<T, r>& s)
{
    return s.begin();
}

template <typename T, size_t r>
constexpr decltype(auto) end(const sd::cpslice<T, r>& s)
{
    return s.end();
}

template <typename T, size_t r>
constexpr decltype(auto) data(const sd::cpslice<T, r>& s)
{
    return s.data();
}

template <typename T, size_t r>
constexpr decltype(auto) size(const sd::cpslice<T, r>& s)
{
    return s.size();
}

template <typename T, size_t r>
constexpr decltype(auto) empty(const sd::cpslice<T, r>& s)
{
    return s.empty();
}

} // namespace sd

#ifdef __cpp_deduction_guides
namespace sd
{
template <typename T, size_t R>
cpslice(sd::slice<T, R>& x) -> cpslice<T, R>;

template <typename T, size_t R>
cpslice(const sd::slice<T, R>& x) -> cpslice<const T, R>;

template <typename V, typename = std::enable_if_t<is_viewable<V>::value>>
cpslice(V&& v) -> cpslice<std::remove_reference_t<decltype(*v.data())>, 1>;

template <typename V, size_t R, typename = std::enable_if_t<is_viewable<V>::value>>
cpslice(V&& v, array<size_t, R> i, ...)
    -> cpslice<std::remove_reference_t<decltype(*v.data())>, R>;

template <typename V, size_t R>
cpslice(V* v, array<size_t, R> i, ...) -> cpslice<std::remove_reference_t<V>, R>;

template <typename V, size_t R, typename = std::enable_if_t<is_viewable<V>::value>>
cpslice(V&& v, layout<R> geo) -> cpslice<std::remove_reference_t<decltype(*v.data())>, R>;

template <typename V, size_t R>
cpslice(V* v, layout<R> geo) -> cpslice<std::remove_reference_t<V>, R>;

template <typename V,
          typename... Is,
          typename = std::enable_if_t<are_integral<Is...>::value && sizeof...(Is) != 0>>
cpslice(V* ptr, Is&&... is) -> cpslice<std::remove_reference_t<V>, sizeof...(Is)>;

template <typename V,
          typename... Is,
          typename = std::enable_if_t<are_integral<Is...>::value && sizeof...(Is) != 0>>
cpslice(V&& v, Is&&... is)
    -> cpslice<std::remove_reference_t<decltype(*v.data())>, sizeof...(Is)>;
} // namespace sd
#endif

#endif

#ifndef MARRAY_COMPONENTWISE_MATH_HXX

#if __cplusplus >= 201703L

#define MARRAY_COMPONENTWISE_MATH_HXX
#define MARRAY_WITH_COMPONENTWISE_MATH 1

namespace sd
{

template <typename T, typename F, std::size_t... I>
void for_each_impl(F&& f, T&& tuple, std::index_sequence<I...>)
{
    (void)std::initializer_list<int>{
        (f(std::get<I>(std::forward<T>(tuple))), void(), int{})...};
}

template <typename T, typename F>
void for_each_in_tuple(F&& f, T&& tuple)
{
    constexpr std::size_t N = std::tuple_size<std::decay_t<T>>::value;
    for_each_impl(std::forward<F>(f), std::forward<T>(tuple), std::make_index_sequence<N>{});
}

template <typename F, typename... X>
class lazy_function
{
public:
    constexpr static size_t rank = std::max<size_t>({std::decay_t<X>::rank...});
    static_assert(rank > 0, "rank == 0");
    using iterator_type = basic_slice_iterator<lazy_function<F, X...>>;
    // using result_type = std::invoke_result_t<F, decltype(std::declval<X>()[{}])...>;
    using index_type = sd::array<size_t, rank>;
    using value_type = std::common_type_t<typename std::decay_t<X>::value_type...>;

    template <typename... Args>
    constexpr lazy_function(F f, Args&&... args) : m_f(f), m_args(std::forward<Args>(args)...)
    {
        std::fill(this->m_extents.begin(), this->m_extents.end(), 1);
        auto broadcast_extents = [this](const auto& v) constexpr
        {
            bool        valid_offset = this->extents().size() >= v.extents().size();
            std::size_t offset       = this->extents().size() - v.extents().size();

            for (std::size_t i = 0; i < v.extents().size(); ++i)
            {
                if (valid_offset && this->extents()[offset + i] == 1)
                {
                    this->m_extents[offset + i] = v.extents()[i];
                }
                else if (v.extents()[i] != this->extents()[offset + i] && v.extents()[i] != 1)
                {
                    throw std::runtime_error("Broadcast error.");
                }
            }
            return true;
        };

        for_each_in_tuple(broadcast_extents, m_args);
    }

    template <std::size_t... I, typename... Args>
    constexpr decltype(auto) access_impl(std::index_sequence<I...>, Args... args) const
    {
        return m_f(std::get<I>(m_args)(args...)...);
    }

    template <typename... Args>
    constexpr decltype(auto) operator()(Args... args) const
    {
        return access_impl(std::make_index_sequence<sizeof...(X)>(), args...);
    }

    template <typename Arg>
    constexpr decltype(auto) operator[](const Arg& index) const
    {
        return access_impl(std::make_index_sequence<sizeof...(X)>(), index);
    }
    constexpr auto begin() const
    {
        return iterator_type(this, index_iterator<rank>(extents(), size() > 0));
    }
    constexpr auto end() const
    {
        return iterator_type(this, index_iterator<rank>(extents(), false));
    }

    constexpr size_t size() const
    {
        return std::accumulate(
            m_extents.begin(), m_extents.end(), size_t(1), std::multiplies<>{});
    }

    constexpr const auto& extents() const { return m_extents; }

    F                m_f;
    index_type       m_extents;
    std::tuple<X...> m_args;
};

template <typename S>
struct closure_type
{
    // using underlying_type =
    // std::conditional_t<std::is_const<std::remove_reference_t<S>>::value,
    //                                            const std::decay_t<S>,
    //                                            std::decay_t<S>>;
    // using type            = typename std::conditional<std::is_lvalue_reference<S>::value,
    //                                        underlying_type&,
    //                                        underlying_type>::type;

    using type =
        typename std::conditional<std::is_lvalue_reference<S>::value, S, std::decay_t<S>>::type;
};

template <typename T>
using closure_type_t = typename closure_type<T>::type;

template <typename T, size_t R>
struct constant_function
{
    constexpr static size_t rank = R;
    using value_type             = T;

    template <typename... Args>
    auto operator()(Args&&...) const
    {
        return value;
    }
    template <typename Arg>
    auto operator[](Arg&&) const
    {
        return value;
    }

    sd::array<size_t, R> extents() const { return sd::array<size_t, R>{}.fill(1); }

    T value;
};

template <size_t rank, typename T>
decltype(auto) wrap_if_scalar(T&& arg)
{
    using value_type = std::decay_t<decltype(arg)>;
    if constexpr (std::is_scalar<value_type>::value)
        return constant_function<value_type, rank>{arg};
    else
        return std::forward<T>(arg);
}

template <typename... Ts>
class are_viewable : public std::integral_constant<bool, true>
{
};

template <typename T, typename... Ts>
class are_viewable<T, Ts...>
    : public std::integral_constant<bool,
                                    sd::is_viewable<T>::value && are_viewable<Ts...>::value>
{
};

template <typename T, typename = void>
struct is_operand_valid : std::false_type
{
};

template <typename T>
struct is_operand_valid<
    T,
    std::void_t<decltype(std::decay_t<T>::rank), decltype(std::declval<T>().extents())>>
    : std::true_type
{
};

template <typename A, typename B>
constexpr bool are_operands_valid()
{
    return (is_operand_valid<A>::value || is_operand_valid<B>::value);
    // && !(std::is_scalar<A>::value && std::is_scalar<B>::value);
}

template <typename A, typename B>
constexpr size_t get_rank()
{
    size_t rank = 1;

    if constexpr (is_operand_valid<B>::value) { rank = std::max(rank, std::decay_t<B>::rank); }
    if constexpr (is_operand_valid<A>::value) { rank = std::max(rank, std::decay_t<A>::rank); }

    return rank;
}

template <typename L, typename... Args>
constexpr auto make_lazy_function(L func, Args&&... args)
{
    return lazy_function<L, closure_type_t<Args>...>(func, std::forward<Args>(args)...);
}

#define BINARY_OP_SFINAE(OP)                                                                   \
    template <typename A, typename B, typename = std::enable_if_t<are_operands_valid<A, B>()>> \
    constexpr auto operator OP(A&& a, B&& b)                                                   \
    {                                                                                          \
        return make_lazy_function([](auto&& x, auto y) { return x OP y; },                     \
                                  wrap_if_scalar<get_rank<A, B>()>(std::forward<A>(a)),        \
                                  wrap_if_scalar<get_rank<A, B>()>(std::forward<B>(b)));       \
    }

#define BINARY_OP(OP)                                                                          \
    template <typename A, typename B>                                                          \
    constexpr auto operator OP(A&& a, B&& b)                                                   \
    {                                                                                          \
        return make_lazy_function([](auto&& x, auto&& y) { return x OP y; },                   \
                                  wrap_if_scalar<get_rank<A, B>()>(std::forward<A>(a)),        \
                                  wrap_if_scalar<get_rank<A, B>()>(std::forward<B>(b)));       \
    }

#define UNARY_OP_SFINAE(OP)                                                                    \
    template <typename A, typename = std::enable_if_t<is_operand_valid<A>::value>>             \
    constexpr auto operator OP(A&& a)                                                          \
    {                                                                                          \
        return make_lazy_function([](auto&& x) { return OP x; },                               \
                                  wrap_if_scalar<std::decay_t<A>::rank>(std::forward<A>(a)));  \
    }

#define UNARY_OP(OP)                                                                           \
    template <typename A>                                                                      \
    constexpr auto operator OP(A&& a)                                                          \
    {                                                                                          \
        return make_lazy_function([](auto&& x) { return OP x; },                               \
                                  wrap_if_scalar<std::decay_t<A>::rank>(std::forward<A>(a)));  \
    }

#define UNARY_FUNC_SFINAE(FUNC, CALL)                                                          \
    template <typename A,                                                                      \
              typename = std::enable_if_t<!std::is_scalar<std::decay_t<A>>::value>,            \
              typename = std::enable_if_t<is_operand_valid<A>::value>>                         \
    constexpr auto FUNC(A&& a)                                                                 \
    {                                                                                          \
        return make_lazy_function(                                                             \
            [](auto&& x)                                                                       \
            {                                                                                  \
                using namespace std;                                                           \
                return CALL(x);                                                                \
            },                                                                                 \
            wrap_if_scalar<std::decay_t<A>::rank>(std::forward<A>(a)));                        \
    }

#define UNARY_FUNC(FUNC, CALL)                                                                 \
    template <typename A>                                                                      \
    constexpr auto FUNC(A&& a)                                                                 \
    {                                                                                          \
        return make_lazy_function(                                                             \
            [](auto&& x)                                                                       \
            {                                                                                  \
                using namespace std;                                                           \
                return CALL(x);                                                                \
            },                                                                                 \
            wrap_if_scalar<std::decay_t<A>::rank>(std::forward<A>(a)));                        \
    }

#define ACCUMULATE_FUNC_SFINAE(FUNC, INITIAL, OP)                                              \
    template <typename A, typename = std::enable_if_t<is_operand_valid<A>::value>>             \
    constexpr auto FUNC(A&& a)                                                                 \
    {                                                                                          \
        using value_type = typename std::decay_t<A>::value_type;                               \
        using namespace std;                                                                   \
        return std::accumulate(cbegin(a), cend(a), static_cast<value_type>(INITIAL), OP);      \
    }

#define ACCUMULATE_FUNC(FUNC, INITIAL, OP)                                                     \
    template <typename A>                                                                      \
    constexpr auto FUNC(A&& a)                                                                 \
    {                                                                                          \
        using value_type = typename std::decay_t<A>::value_type;                               \
        using namespace std;                                                                   \
        return std::accumulate(cbegin(a), cend(a), static_cast<value_type>(INITIAL), OP);      \
    }

namespace math
{

BINARY_OP(+)
BINARY_OP(-)
BINARY_OP(*)
BINARY_OP(/)
// BINARY_OP(%)
// BINARY_OP(&)
// BINARY_OP(|)
// BINARY_OP(^)
// BINARY_OP(==)
// BINARY_OP(!=)
BINARY_OP(<=)
BINARY_OP(<)
BINARY_OP(>=)
BINARY_OP(>)
UNARY_OP(+)
UNARY_OP(-)

UNARY_FUNC(fabs, ::fabs)
UNARY_FUNC(fmod, ::fmod)
UNARY_FUNC(remainder, ::remainder)
UNARY_FUNC(remquo, ::remquo)
UNARY_FUNC(fma, ::fma)
UNARY_FUNC(fmax, ::fmax)
UNARY_FUNC(fmin, ::fmin)
UNARY_FUNC(fdim, ::fdim)
UNARY_FUNC(exp, ::exp)
UNARY_FUNC(exp2, ::exp2)
UNARY_FUNC(expm1, ::expm1)
UNARY_FUNC(log, ::log)
UNARY_FUNC(log10, ::log10)
UNARY_FUNC(log2, ::log2)
UNARY_FUNC(log1p, ::log1p)
UNARY_FUNC(pow, ::pow)
UNARY_FUNC(sqrt, ::sqrt)
UNARY_FUNC(cbrt, ::cbrt)
UNARY_FUNC(hypot, ::hypot)
UNARY_FUNC(sin, ::sin)
UNARY_FUNC(cos, ::cos)
UNARY_FUNC(tan, ::tan)
UNARY_FUNC(asin, ::asin)
UNARY_FUNC(acos, ::acos)
UNARY_FUNC(atan, ::atan)
UNARY_FUNC(atan2, ::atan2)
UNARY_FUNC(sinh, ::sinh)
UNARY_FUNC(cosh, ::cosh)
UNARY_FUNC(tanh, ::tanh)
UNARY_FUNC(asinh, ::asinh)
UNARY_FUNC(acosh, ::acosh)
UNARY_FUNC(atanh, ::atanh)
UNARY_FUNC(erf, ::erf)
UNARY_FUNC(erfc, ::erfc)
UNARY_FUNC(tgamma, ::tgamma)
UNARY_FUNC(lgamma, ::lgamma)
UNARY_FUNC(ceil, ::ceil)
UNARY_FUNC(floor, ::floor)
UNARY_FUNC(trunc, ::trunc)
UNARY_FUNC(round, ::round)
UNARY_FUNC(nearbyint, ::nearbyint)
UNARY_FUNC(rint, ::rint)
UNARY_FUNC(frexp, ::frexp)
UNARY_FUNC(ldexp, ::ldexp)
UNARY_FUNC(modf, ::modf)
UNARY_FUNC(scalbn, ::scalbn)
UNARY_FUNC(ilogb, ::ilogb)
UNARY_FUNC(logb, ::logb)
UNARY_FUNC(nextafter, ::nextafter)
UNARY_FUNC(copysign, ::copysign)
UNARY_FUNC(fpclassify, std::fpclassify)
UNARY_FUNC(isfinite, std::isfinite)
UNARY_FUNC(isinf, std::isinf)
UNARY_FUNC(isnan, std::isnan)
UNARY_FUNC(isnormal, std::isnormal)
UNARY_FUNC(signbit, std::signbit)
UNARY_FUNC(iszero, [](const auto& x) { return x == 0; })

template <typename A,
          typename Fn,
          typename = std::enable_if_t<!std::is_scalar<std::decay_t<A>>::value>,
          typename = std::enable_if_t<is_operand_valid<A>::value>>
constexpr auto map(A&& a, Fn&& fn)
{
    return make_lazy_function([f = std::forward<Fn>(fn)](const auto& x) { return f(x); },
                              wrap_if_scalar<std::decay_t<A>::rank>(std::forward<A>(a)));
}

ACCUMULATE_FUNC(sum, 0, std::plus<>())
ACCUMULATE_FUNC(mul, 1, std::multiplies<>())

} // namespace math

#undef BINARY_OP
#undef UNARY_OP
#undef UNARY_FUNC
#undef ACCUMULATE_FUNC
#undef BINARY_OP_SFINAE
#undef UNARY_OP_SFINAE
#undef UNARY_FUNC_SFINAE
#undef ACCUMULATE_FUNC_SFINAE

template <typename T, typename Mask>
class masking_function
{
public:
    using underlying_type = std::decay_t<T>;

    static constexpr auto rank = underlying_type::rank;
    using value_type           = typename underlying_type::value_type;
    using index_type           = typename underlying_type::index_type;
    using iterator             = basic_slice_iterator<masking_function<T, Mask>>;

    template <typename S>
    constexpr decltype(auto) operator[](S&& s) const
    {
        auto pos = what.shape().unlinearize(mask[std::forward<S>(s)]);
        return what[pos];
        // return what[what[std::forward<S>(s)]];
    }

    template <typename... S>
    constexpr decltype(auto) operator()(S&&... s) const
    {
        auto pos = what.shape().unlinearize(mask(std::forward<S>(s)...));
        return what[pos];
        // return what(mask(std::forward<S>(s)...));
    }

    constexpr decltype(auto) extents() const { return mask.extents(); }
    constexpr decltype(auto) size() const { return mask.size(); }

    constexpr auto begin() const
    {
        return iterator(this, index_iterator<rank>(extents(), size() > 0));
    }
    constexpr auto end() const
    {
        return iterator(this, index_iterator<rank>(extents(), false));
    }

    T    what;
    Mask mask;
};

template <typename A,
          typename MASK,
          typename = std::enable_if_t<!std::is_scalar<std::decay_t<A>>::value>,
          typename = std::enable_if_t<is_operand_valid<A>::value>>
constexpr auto mask(A&& a, MASK&& b)
{
    return masking_function<closure_type_t<A>, closure_type_t<MASK>>({a, b});
}

// struct masked_iterator_sentinel {};

template <typename data_iterator, typename mask_iterator>
struct masked_iterator
{
    data_iterator iter_data, data_end;
    mask_iterator iter_mask, mask_end;

    masked_iterator(data_iterator iter_data,
                    data_iterator data_end,
                    mask_iterator iter_mask,
                    mask_iterator mask_end)
        : iter_data(iter_data), data_end(data_end), iter_mask(iter_mask), mask_end(mask_end)
    {
        if (iter_mask != mask_end && !*iter_mask) operator++();
    }

    void operator++()
    {
        while (iter_mask != mask_end && iter_data != data_end)
        {
            ++iter_data;
            ++iter_mask;

            if (iter_mask == mask_end || *iter_mask) break;
        }

        if (iter_mask == mask_end) { iter_data = data_end; }
        if (iter_data == data_end) { iter_mask = mask_end; }
    }
    decltype(auto) operator*() { return *iter_data; }
    bool           operator!=(const masked_iterator& other)
    {
        return iter_mask != mask_end && iter_data != data_end && iter_mask != other.iter_mask &&
               iter_data != other.iter_data;
    }
};

template <typename A, typename MASK>
/*constexpr*/ auto lazy_mask(A& a, MASK&& b)
{
    assert(a.extents() == b.extents());

    using data_iter_t = std::remove_reference_t<decltype(a.begin())>;
    using mask_iter_t = std::remove_reference_t<decltype(b.begin())>;
    using iter_t      = masked_iterator<data_iter_t, mask_iter_t>;
    struct range
    {
        A&     a;
        MASK&& b;
        iter_t begin() { return {a.begin(), a.end(), b.begin(), b.end()}; }
        iter_t end() { return {a.end(), a.end(), b.end(), b.end()}; }
        // masked_iterator_sentinel end() { return {}; }
    };

    return range{a, std::forward<MASK>(b)};
}

template <typename A,
          typename MASK,
          typename = std::enable_if_t<!std::is_scalar<std::decay_t<A>>::value>,
          typename = std::enable_if_t<is_operand_valid<A>::value>>
/*constexpr*/ auto operator&(A& a, MASK&& b)
{
    std::vector<size_t> mask;
    mask.reserve(b.size());

    auto flatten_out = a.flatten();

    for (auto coord : a.subscripts())
    {
        const bool use_coord = static_cast<bool>(b[coord]);
        if (use_coord)
        {
            auto linearized_index = a.shape().linearize(coord);
            mask.push_back(linearized_index);
        }
    }
    using masking_function_type =
        masking_function<closure_type_t<decltype(std::move(flatten_out))>, std::vector<size_t>>;

    return masking_function_type({std::move(flatten_out), std::move(mask)});

    // return masking_function<closure_type_t<A>, closure_type_t<MASK>>({a, b});
}

template <typename A, size_t R, typename B>
auto mask_predicate(sd::cpslice<A, R> a, const B& b)
{
    return masking_function<sd::cpslice<A>, const B&>({a, b});
}

// // template<typename F, typename ... Args>
// // struct generic_lazy_function
// // {
// //     using Lambda = F;
// //     using iterator_type = basic_slice_iterator<generic_lazy_function<F, Args...>>;

// //     Lambda              lambda;
// //     std::tuple<Args...> args;
// //     sd::layout<Rank> geo;

// //     constexpr auto extents() const noexcept { return geo.extents(); }
// //     constexpr auto strides() const noexcept { return geo.strides(); }

// //     constexpr iterator begin() const noexcept { return {this, true}; }
// //     constexpr iterator end() const noexcept { return {this, false}; }

// //     template <typename... Pos>
// //     decltype(auto) operator()(Pos&&... pos)
// //     {
// //         return lambda(args, pos...);
// //     }
// //     template <typename... Pos>
// //     decltype(auto) operator()(index_type i)
// //     {
// //         return lambda(args, i);
// //     }
// // };

// // template<typename F>
// // struct generic_lazy_function
// // {
// //     using Lambda = F;
// //     using iterator_type = basic_slice_iterator<generic_lazy_function<F>>;

// //     Lambda              lambda;
// //     sd::layout<Rank> geo;

// //     constexpr auto extents() const noexcept { return geo.extents(); }
// //     constexpr auto strides() const noexcept { return geo.strides(); }

// //     constexpr iterator begin() const noexcept { return {this, true}; }
// //     constexpr iterator end() const noexcept { return {this, false}; }

// //     template <typename... Pos>
// //     decltype(auto) operator()(Pos&&... pos)
// //     {
// //         return lambda(pos...);
// //     }
// //     template <typename... Pos>
// //     decltype(auto) operator()(index_type i)
// //     {
// //         return lambda(i);
// //     }
// // };

// template <class A, class B, typename = std::enable_if_t<are_operands_valid<A, B>()>>
// auto operator select(A&& a, B&& b)
// {

//     template <typename Lambda>
//     struct lambda_wrapper
//     {
//         Lambda lambda;
//         template <typename... Pos>
//         decltype(auto) operator()(index_type i)
//         {
//             return lambda(std::get<I>(args)(pos...)...);
//         }
//     }

//     auto
//     make_lazy_function()
//     {
// generic_lazy_function(lambda_wrapper([](){});;;
//     }

//     template<typename S>
//     auto select()
//     {
//         return generic_lazy_function(
//             B.shape(),
//             [](auto& A, const auto& F, Pos... pos) { return A(F(pos...)); },
//             std::make_tuple(A, B));
//     }
// }

} // namespace sd

#endif
#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_NDARRAY_HXX
#define MARRAY_NDARRAY_HXX

#if __has_include(<marray/componentwise_math.hxx>)
#endif

namespace sd
{

template <typename V,
          typename = std::enable_if_t<sd::is_viewable<V>::value>,
          typename = decltype(V::rank),
          typename = decltype(std::declval<V>().shape()),
          typename = decltype(std::declval<V>().data()),
          typename = decltype(std::declval<V>().begin()),
          typename = decltype(std::declval<V>().end())>
struct ndarray_constructable_from
{
};

template <typename T, size_t R = 1, typename A = std::allocator<T>>
class ndarray : public cpslice<T, R>
{
    using base = cpslice<T, R>;
    using base::get;

public:
    using allocator = A;
    using base::rank;
    using typename base::index_type;
    using typename base::iterator;
    using typename base::layout_type;
    using typename base::pointer;
    using typename base::size_type;
    using typename base::value_type;

    ndarray()          = default;
    ndarray(ndarray&&) = default;
    ndarray& operator=(ndarray&&) = default;
    ndarray(const ndarray& s) : base{}, elems(s.elems) { get() = {elems, s.shape()}; }
    ndarray& operator=(const ndarray& s)
    {
        elems = s.elems;
        get() = {elems, s.shape()};
        return *this;
    }
    ndarray(const T& init, layout<R> geo) : elems(geo.total_size(), init)
    {
        get() = slice<T, R>{elems.data(), std::move(geo)};
    }
    ndarray(const T& init, const array<size_t, R>& extents, const array<size_t, R>& strides)
        : ndarray(init, {extents, strides})
    {
    }
    ndarray(const T& init, const array<size_t, R>& extents) : ndarray(init, layout<R>{extents})
    {
    }
    template <typename... Is,
              typename = std::enable_if_t<sizeof...(Is) == R && are_integral<Is...>::value>>
    ndarray(const T& init, const Is&... is) : ndarray(init, index_type{is...})
    {
    }
    template <typename... Is,
              typename = std::enable_if_t<sizeof...(Is) == R && are_integral<Is...>::value>>
    ndarray(const Is&... is) : ndarray(T{}, index_type{is...})
    {
    }
    template <typename S, typename = sd::ndarray_constructable_from<S>>
    ndarray(const S& s) : elems(s.begin(), s.end())
    {
        get() = {elems.data(), s.shape()};
    }
    template <typename S, typename = sd::ndarray_constructable_from<S>>
    ndarray& operator=(const S& s)
    {
        elems.assign(s.begin(), s.end());
        get() = {elems.data(), s.shape()};
        return *this;
    }
    template <typename Alloc, bool cond = R == 1, typename = std::enable_if_t<cond>>
    ndarray(std::vector<T, Alloc> s) : elems(std::move(s))
    {
        get() = elems;
    }

#if defined(MARRAY_WITH_COMPONENTWISE_MATH)
    template <typename... Ts>
    ndarray(const lazy_function<Ts...>& f)
        : ndarray(layout<lazy_function<Ts...>::rank>(f.extents()))
    {
        for (auto i : this->subscripts()) { this->operator[](i) = f(i); }
    }
    template <typename... Ts>
    ndarray& operator=(const lazy_function<Ts...>& f)
    {
        resize(layout<lazy_function<Ts...>::rank>(f.extents()));
        for (auto i : this->subscripts()) { this->operator[](i) = f(i); }
        return *this;
    }
#endif

    // template<typename ... Ts,
    //     typename = std::enable_if_t<std::is_constructible<slice<const T, rank>,
    //     Ts...>::value>>
    // ndarray(Ts&&... ts) : ndarray(slice<const T, rank>(std::forward<Ts>(ts)...)) {
    // }

    // ndarray(slice<T, R> a) : elems(a.begin(), a.end())
    // {
    //     get() = base{elems.data(), a.shape()};
    // }
    // ndarray(slice<const T, R> a) : elems(a.begin(), a.end())
    // {
    //     get() = base{elems.data(), a.shape()};
    // }
    // ndarray(cpslice<T, R> a) : elems(a.begin(), a.end())
    // {
    //     get() = base{elems.data(), a.shape()};
    // }

    // template<typename S, typename = std::enable_if_t<R != 1>>
    // ndarray(slice<const T> a) : elems(a.begin(), a.end())
    // {
    //     get() = base{elems.data(), a.shape()};
    // }
    // ndarray& operator=(const slice<const T, R> rhs)
    // {
    //     elems.assign(rhs.begin(), rhs.end());
    //     get() = base{elems.data(), rhs.shape()};
    // }

    // ndarray& operator=(const cpslice<T> rhs)
    // {
    //     elems.assign(rhs.begin(), rhs.end());
    //     get() = base{elems.data(), rhs.shape()};
    // }

    // ndarray(const layout<R>& geo) : elems(geo.total_size())
    // {
    //     get() = slice<T, R>{elems.data(), geo};
    // }
    //     ndarray(const array<size_t, R>& extents, const array<size_t, R>& strides)
    //         : ndarray({extents, strides})
    //     {
    //     }
    // template <
    //     typename U,
    //     size_t S>
    //     // ,
    //     // typename = std::enable_if_t</*S != R && R == 1 && */
    //     //                             is_safely_convertable_value<std::remove_cv_t<U>,
    //     T>::value>>
    // ndarray(slice<U, S> a) : elems(a.begin(), a.end())
    // {
    //     get() = {elems.data(), a.shape()};
    // }
    //     template <typename U, size_t S>
    //     ndarray(const cpslice<U, S>& s) : ndarray(s.get())
    //     {
    //     }

    //     template <typename Seq,
    //               typename = std::enable_if_t<!std::is_same<value_type,
    //               std::decay_t<Seq>>::value>, typename =
    //               std::enable_if_t<is_viewable<Seq>::value>>
    //     ndarray& operator=(const Seq& s)
    //     {
    //         return operator=(slice<const T>(s));
    //     }
    // template <typename U = T, size_t S = R>
    // ndarray& operator=(slice<U, S> rhs)
    // {
    //     elems = {rhs.begin(), rhs.end()};
    //     get() = {elems.data(), rhs.shape()};
    //     return *this;
    // }
    // template <typename U = T, size_t S = R>
    // ndarray& operator=(const cpslice<U, S> rhs)
    // {
    //     return operator=(rhs.get());
    // }

    //     ndarray& operator=(const ndarray& rhs) { return this->operator=(slice<const T,
    //     R>(rhs)); } ndarray(const ndarray& rhs) : ndarray(slice<const T, R>(rhs)) {}

    void reshape(const layout<R>& geo, const T& init = T())
    {
        auto next = *this;
        elems.clear();
        elems.resize(geo.total_size(), init);
        get() = {elems.data(), geo};

        for (auto s : next.subscripts()) { this->operator()(s) = next(s); }
    }

    void resize(const layout<R>& geo, const T& init = T())
    {
        elems.resize(geo.total_size(), init);
        get() = {elems.data(), geo};
    }
    template <typename... Is,
              typename = std::enable_if_t<sizeof...(Is) == R && are_integral<Is...>::value>>
    void resize(const Is&... is)
    {
        resize(array<size_t, R>{is...});
    }
    void clear()
    {
        elems.clear();
        get() = {};
    }

    void reserve(const layout<rank>& n) { elems.reserve(n.total_size()); }
    void reserve(const size_t n) { elems.reserve(n); }

    void shrink_to_fit()
    {
        elems.shrink_to_fit();
        get() = {elems.data(), get().shape()};
    }

    auto&       container() { return elems; }
    const auto& container() const { return elems; }

    /// insert a row with a length of extent(1) * ... * extent(rank) initialized with value of x
    /// Important Note: This function will reset the layout, i.e. *this will be in
    /// row-major-order after calling this function
    void push_back(value_type x = value_type{})
    {
        if (this->size() == 0)
        {
            array<size_t, rank> ext;
            ext.fill(1);
            --ext[0];
            get() = {elems.data(), {ext}};
        }

        auto next_layout = layout<rank>{this->extents() + first_component()};
        if /*constexpr*/ (rank == 1) { elems.emplace_back(std::move(x)); }
        else
        {
            auto increase_size = next_layout.total_size() - get().shape().total_size();
            if (increase_size == 0)
                throw std::runtime_error{"cannot increase size for this layout (slice)"};
            elems.insert(elems.end(), increase_size, x);
        }
        get() = {elems.data(), next_layout};
    }

    template <typename Seq,
              typename = std::enable_if_t<!std::is_same<value_type, std::decay_t<Seq>>::value>,
              typename = decltype(std::declval<Seq>().begin()),
              typename = decltype(std::declval<Seq>().end())>
    void push_back(const Seq& x)
    {
        push_back(x.begin(), x.end());
    }

    template <typename Iter>
    void push_back(Iter first, Iter last)
    {

        if (this->size() == 0)
        {
            array<size_t, rank> ext;
            ext.fill(1);
            --ext[0];
            ext[rank == 1 ? 0 : 1] = std::distance(first, last);
            get()                  = {elems.data(), {ext}};
        }

        auto next_layout   = layout<rank>{this->extents() + first_component()};
        auto increase_size = next_layout.total_size() - get().shape().total_size();
        if (increase_size == 0)
            throw std::runtime_error{"cannot increase size for this layout (slice)"};
        if (std::distance(first, last) != static_cast<std::ptrdiff_t>(increase_size))
            throw std::invalid_argument{"argument has wrong size (slice)"};
        elems.insert(elems.end(), first, last);
        get() = {elems.data(), next_layout};
    }

    /// Important Note: This function will reset the layout, i.e. *this will be in
    /// row-major-order after calling this function
    void pop_back()
    {
        auto next_layout = layout<rank>{this->extents() - first_component()};

        if /*constexpr*/ (rank == 1) { elems.pop_back(); }
        else
        {
            auto decrease_size = get().shape().total_size() - next_layout.total_size();
            if (decrease_size == 0)
                throw std::runtime_error{"cannot increase size for this layout (slice)"};
            elems.erase(elems.end() - decrease_size, elems.end());
        }
        get() = {elems.data(), next_layout};
    }

private:
    constexpr static array<size_t, rank> first_component()
    {
        sd::array<size_t, rank> arr;
        arr[0] = 1;
        return arr;
    }

    static_assert(!std::is_same<T, bool>::value, "T is of type bool (ndarray)");

    std::vector<T, allocator> elems;
};

#ifdef __cpp_deduction_guides
template <typename T,
          typename... Is,
          typename = std::enable_if_t<sizeof...(Is) != 0 && are_integral<Is...>::value>>
ndarray(const T& init, Is&&... is) -> ndarray<T, sizeof...(Is), std::allocator<T>>;

template <typename U, size_t R>
ndarray(const slice<U, R>& a)
    -> ndarray<std::remove_reference_t<std::remove_cv_t<U>>,
               R,
               std::allocator<std::remove_reference_t<std::remove_cv_t<U>>>>;

template <typename U, size_t R>
ndarray(const cpslice<U, R>& a)
    -> ndarray<std::remove_reference_t<std::remove_cv_t<U>>,
               R,
               std::allocator<std::remove_reference_t<std::remove_cv_t<U>>>>;

#if defined(MARRAY_WITH_COMPONENTWISE_MATH)
template <typename... Ts>
ndarray(const lazy_function<Ts...>&)
    -> ndarray<typename lazy_function<Ts...>::value_type, lazy_function<Ts...>::rank>;
#endif
#endif

} // namespace sd

#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_FIXED_NDARRAY_HXX
#define MARRAY_FIXED_NDARRAY_HXX

namespace sd
{

template <typename T, size_t Ext0, size_t... Ext>
class fixed_ndarray : public cpslice<T, sizeof...(Ext) + 1>
{
public:
    constexpr static size_t rank = sizeof...(Ext) + 1;

private:
    constexpr static layout<rank> geometry{{Ext0, Ext...}};
    constexpr static size_t       tot_size = geometry.total_size();

public:
    using container_type = sd::array<T, tot_size>;

    constexpr fixed_ndarray() { this->get() = {field.data(), geometry}; }
    constexpr fixed_ndarray(const fixed_ndarray& rhs)
    {
        stdext::copy(rhs.field.begin(), rhs.field.end(), field.begin());
        this->get() = {field.data(), geometry};
    }
    constexpr fixed_ndarray& operator=(const fixed_ndarray& rhs)
    {
        stdext::copy(rhs.field.begin(), rhs.field.end(), field.begin());
        this->get() = {field.data(), geometry};
        return *this;
    }
    constexpr fixed_ndarray(const container_type& rhs) : field(rhs)
    {
        this->get() = {field.data(), geometry};
    }
    constexpr fixed_ndarray& operator=(const container_type& rhs)
    {
        field       = rhs;
        this->get() = {field.data(), geometry};
        return *this;
    }

    const container_type& container() const { return field; }
    container_type&       container() { return field; }

private:
    container_type field{};
};

} // namespace sd

#endif
/// See "LICENSE" for the full license governing this code.

#ifndef MARRAY_ROW_ITERATOR_HXX
#define MARRAY_ROW_ITERATOR_HXX

namespace sd
{

template <typename T, size_t rank = 2>
class row_iterator
{
    slice<T, rank> data;
    const size_t   kfrom, kto;
    size_t         row;

    using Iterator = row_iterator<T, rank>;

    bool      is_reverse() const { return kfrom >= kto; }
    ptrdiff_t direction() const { return 2 * !is_reverse() - 1; }

    row_iterator(row_iterator other, size_t new_row)
        : data{other.data}, kfrom(other.kfrom), kto(other.kto), row(new_row)
    {
    }

public:
    using iterator_category = std::random_access_iterator_tag;
    using value_type        = slice<T>;
    using reference         = value_type&;
    using pointer           = value_type*;
    using difference_type   = std::ptrdiff_t;
    using const_pointer     = value_type const*;
    using const_reference   = value_type const&;

    row_iterator(slice<T, rank> data, bool last, bool reverse = false)
        : data(data)
        , kfrom(reverse ? data.extent(0) : 0)
        , kto(reverse ? 0 : data.extent(0))
        , row(!last ? kfrom : kto)
    {
    }

    operator difference_type() const { return row; }

    auto operator*() const { return data[row]; }
    auto operator->() const { return data[row].data(); }
    auto operator[](difference_type r) const { return data[r]; }

    auto operator++()
    {
        row += direction();
        return *this;
    }
    auto operator--()
    {
        row += -direction();
        return *this;
    }
    auto operator++(int)
    {
        auto retval = *this;
        ++(*this);
        return retval;
    }
    auto operator--(int)
    {
        auto retval = *this;
        --(*this);
        return retval;
    }

    Iterator& operator+=(difference_type rhs)
    {
        row += direction() * rhs;
        return *this;
    }
    Iterator& operator-=(difference_type rhs)
    {
        row -= direction() * rhs;
        return *this;
    }
    Iterator operator+(difference_type rhs) const
    {
        return Iterator(*this, row + direction() * rhs);
    }
    Iterator operator-(difference_type rhs) const
    {
        return Iterator(*this, row - direction() * rhs);
    }
    Iterator operator+(const Iterator& rhs) const
    {
        return {*this, row + direction() * rhs.row};
    }
    Iterator operator-(const Iterator& rhs) const
    {
        return {*this, row - direction() * rhs.row};
    }
    friend Iterator operator+(difference_type lhs, const Iterator& rhs)
    {
        return Iterator(rhs, rhs.direction() * lhs + rhs.row);
    }
    friend Iterator operator-(difference_type lhs, const Iterator& rhs)
    {
        return Iterator(rhs, rhs.direction() * lhs - rhs.row);
    }
    bool operator==(const Iterator& rhs) const { return row == rhs.row; }
    bool operator!=(const Iterator& rhs) const { return !(*this == rhs); }
    bool operator>(const Iterator& rhs) const { return row > rhs.row; }
    bool operator<(const Iterator& rhs) const { return row < rhs.row; }
    bool operator>=(const Iterator& rhs) const { return row >= rhs.row; }
    bool operator<=(const Iterator& rhs) const { return row <= rhs.row; }
};

template <typename T, size_t rank>
class row_range
{
public:
    using iterator       = row_iterator<T, rank>;
    using const_iterator = row_iterator<const T, rank>;

    iterator begin() { return {data, false}; }
    iterator end() { return {data, true}; }
    iterator rbegin() { return {data, false, true}; }
    iterator rend() { return {data, true, true}; }

    const_iterator begin() const { return {data, false}; }
    const_iterator end() const { return {data, true}; }
    const_iterator rbegin() const { return {data, false, true}; }
    const_iterator rend() const { return {data, true, true}; }

    slice<T, rank> data;
};

template <typename T, size_t rank>
row_range<T, rank> iterate_rows(slice<T, rank> a)
{
    return {a};
}

template <typename T, size_t rank>
row_range<T, rank> iterate_rows(cpslice<T, rank> a)
{
    return {a};
}

} // namespace sd

#endif

/// See "LICENSE" for the full license governing this code.
#ifndef MARRAY_MARRAY_HXX
#define MARRAY_MARRAY_HXX
#endif

#endif
