#ifndef XTENSOR_FUNCTION_HPP

#define XTENSOR_FUNCTION_HPP

#include <algorithm>

#include <cstddef>

#include <iterator>

#include <numeric>

#include <tuple>

#include <type_traits>

#include <utility>

#include <xtl/xsequence.hpp>

#include <xtl/xtype_traits.hpp>

#include "xaccessible.hpp"

#include "xexpression_traits.hpp"

#include "xiterable.hpp"

#include "xlayout.hpp"

#include "xscalar.hpp"

#include "xshape.hpp"

#include "xstrides.hpp"

#include "xtensor_simd.hpp"

#include "xutils.hpp"

namespace xt

{

namespace detail

{

template <bool... B>

using conjunction_c = xtl::conjunction<std::integral_constant<bool, B>...>;

/************************

template<class S, class is_shape_trivial>

struct xfunction_cache_impl

};

template<std::size_t... N, class is_shape_trivial>

struct xfunction_cache_impl<fixed_shape<N...>, is_shape_trivial>

{

XTENSOR_CONSTEXPR_ENHANCED_STATIC fixed_shape<N...> shape = fixed_shape<N...>();

XTENSOR_CONSTEXPR_ENHANCED_STATIC bool is_trivial = is_shape_trivial::value;

XTENSOR_CONSTEXPR_ENHANCED_STATIC bool is_initialized = true;

};

#ifdef XTENSOR_HAS_CONSTEXPR_ENHANCED

// Out of line definitions to prevent linker errors prior to C++17

template <std::size_t... N, class is_shape_trivial>

constexpr fixed_shape<N...> xfunction_cache_impl<fixed_shape<N...>, is_shape_trivial>::shape;

template <std::size_t... N, class is_shape_trivial>

constexpr bool xfunction_cache_impl<fixed_shape<N...>, is_shape_trivial>::is_trivial;

template <std::size_t... N, class is_shape_trivial>

constexpr bool xfunction_cache_impl<fixed_shape<N...>, is_shape_trivial>::is_initialized;

#endif

template <class... CT>

struct xfunction_bool_load_type

{

using type = xtl::promote_type_t<typename std::decay_t<CT>::bool_load_type...>;

};

template <class CT>

struct xfunction_bool_load_type<CT>

{

using type = typename std::decay_t<CT>::bool_load_type;

};

template <class... CT>

using xfunction_bool_load_type_t = typename xfunction_bool_load_type<CT...>::type;

}

/************************

namespace extension

{

template <class Tag, class F, class... CT>

struct xfunction_base_impl;

template <class F, class... CT>

struct xfunction_base_impl<xtensor_expression_tag, F, CT...>

{

using type = xtensor_empty_base;

};

template <class F, class... CT>

struct xfunction_base

: xfunction_base_impl<xexpression_tag_t<CT...>, F, CT...>

{

};

template <class F, class... CT>

using xfunction_base_t = typename xfunction_base<F, CT...>::type;

}

template<class promote>

struct xfunction_cache : detail::xfunction_cache_impl<typename promote::type, promote> {};

template <class F, class... CT>

class xfunction_iterator;

template <class F, class... CT>

class xfunction_stepper;

template <class F, class... CT>

class xfunction;

template <class F, class... CT>

struct xiterable_inner_types<xfunction<F, CT...>>

{

using inner_shape_type = promote_shape_t<typename std::decay_t<CT>::shape_type...>;

using const_stepper = xfunction_stepper<F, CT...>;

using stepper = const_stepper;

};

template <class F, class... CT>

struct xcontainer_inner_types<xfunction<F, CT...>>

{

// Added indirection for MSVC 2017 bug with the operator value_type()

using value_type = typename meta_identity<decltype(std::declval<F>()(std::declval<xvalue_type_t<std::decay_t<CT>>>()...))>::type;

using reference = value_type;

using const_reference = value_type;

using size_type = common_size_type_t<std::decay_t<CT>...>;

};

template <class T, class F, class... CT>

struct has_simd_interface<xfunction<F, CT...>, T>

: xtl::conjunction<has_simd_type<T>,

has_simd_apply<F, xt_simd::simd_type<T>>,

has_simd_interface<std::decay_t<CT>, T>...>

{

};

/*************

/**

template <class F, class... CT>

class xfunction : private xconst_iterable<xfunction<F, CT...>>,

};

/**********************

template <class F, class... CT>

class xfunction_iterator : public xtl::xrandom_access_iterator_base<xfunction_iterator<F, CT...>,

};

template <class F, class... CT>

bool operator==(const xfunction_iterator<F, CT...>& it1,

const xfunction_iterator<F, CT...>& it2);

template <class F, class... CT>

bool operator<(const xfunction_iterator<F, CT...>& it1,

const xfunction_iterator<F, CT...>& it2);

/*********************

template <class F, class... CT>

class xfunction_stepper

};

/*********************************

/**

//@{

/**

template <class F, class... CT>

template <class Func, class... CTA, class U>

inline xfunction<F, CT...>::xfunction(Func&& f, CTA&&... e) noexcept

: m_e(std::forward<CTA>(e)...), m_f(std::forward<Func>(f)) {}

/**

template <class F, class... CT>

template <class FA, class... CTA>

inline xfunction<F, CT...>::xfunction(xfunction<FA, CTA...> xf) noexcept

: m_e(xf.arguments()), m_f(xf.functor()) {}

//@}

/**

//@{

/**

template <class F, class... CT>

inline auto xfunction<F, CT...>::dimension() const noexcept -> size_type

{

size_type dimension = m_cache.is_initialized ? m_cache.shape.size() : compute_dimension();

return dimension;

}

template <class F, class... CT>

inline void xfunction<F, CT...>::compute_cached_shape() const

{

static_assert(!detail::is_fixed<shape_type>::value, "Calling compute_cached_shape on fixed!");

m_cache.shape = uninitialized_shape<xindex_type_t<inner_shape_type>>(compute_dimension());

m_cache.is_trivial = broadcast_shape(m_cache.shape, false);

m_cache.is_initialized = true;

}

/**

template <class F, class... CT>

inline auto xfunction<F, CT...>::shape() const -> const inner_shape_type&

{

xtl::mpl::static_if<!detail::is_fixed<inner_shape_type>::value>([&](auto self)

{

if(!m_cache.is_initialized)

{

self(this)->compute_cached_shape();

}

},

[](auto /*self*/){});

return m_cache.shape;

}

/**

template <class F, class... CT>

inline layout_type xfunction<F, CT...>::layout() const noexcept

{

return layout_impl(std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

inline bool xfunction<F, CT...>::is_contiguous() const noexcept

{

return layout() != layout_type::dynamic && accumulate([](bool r, const auto& exp) { return r && exp.is_contiguous(); }, true, m_e);

}

//@}

/**

/**

template <class F, class... CT>

template <class... Args>

inline auto xfunction<F, CT...>::operator()(Args... args) const -> const_reference

{

// The static cast prevents the compiler from instantiating the template methods with signed integers,

// leading to warning about signed/unsigned conversions in the deeper layers of the access methods

return access_impl(std::make_index_sequence<sizeof...(CT)>(), static_cast<size_type>(args)...);

}

/**

template <class F, class... CT>

template <class... Args>

inline auto xfunction<F, CT...>::unchecked(Args... args) const -> const_reference

{

// The static cast prevents the compiler from instantiating the template methods with signed integers,

// leading to warning about signed/unsigned conversions in the deeper layers of the access methods

return unchecked_impl(std::make_index_sequence<sizeof...(CT)>(), static_cast<size_type>(args)...);

}

/**

template <class F, class... CT>

template <class It>

inline auto xfunction<F, CT...>::element(It first, It last) const -> const_reference

{

return element_access_impl(std::make_index_sequence<sizeof...(CT)>(), first, last);

}

//@}

/**

//@{

/**

template <class F, class... CT>

template <class S>

inline bool xfunction<F, CT...>::broadcast_shape(S& shape, bool reuse_cache) const

{

if(m_cache.is_initialized && reuse_cache)

{

std::copy(m_cache.shape.cbegin(), m_cache.shape.cend(), shape.begin());

return m_cache.is_trivial;

}

else

{

// e.broadcast_shape must be evaluated even if b is false

auto func = [&shape](bool b, auto&& e) { return e.broadcast_shape(shape) && b; };

return accumulate(func, true, m_e);

}

}

/**

template <class F, class... CT>

template <class S>

inline bool xfunction<F, CT...>::has_linear_assign(const S& strides) const noexcept

{

auto func = [&strides](bool b, auto&& e) { return b && e.has_linear_assign(strides); };

return accumulate(func, true, m_e);

}

//@}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_begin() const noexcept -> const_storage_iterator

{

return storage_cbegin();

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_end() const noexcept -> const_storage_iterator

{

return storage_cend();

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_cbegin() const noexcept -> const_storage_iterator

{

auto f = [](const auto& e) noexcept { return linear_begin(e); };

return build_iterator(f, std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_cend() const noexcept -> const_storage_iterator

{

auto f = [](const auto& e) noexcept { return linear_end(e); };

return build_iterator(f, std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_rbegin() const noexcept -> const_reverse_storage_iterator

{

return storage_crbegin();

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_rend() const noexcept -> const_reverse_storage_iterator

{

return storage_crend();

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_crbegin() const noexcept -> const_reverse_storage_iterator

{

return const_reverse_storage_iterator(storage_cend());

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::storage_crend() const noexcept -> const_reverse_storage_iterator

{

return const_reverse_storage_iterator(storage_cbegin());

}

template <class F, class... CT>

template <class S>

inline auto xfunction<F, CT...>::stepper_begin(const S& shape) const noexcept -> const_stepper

{

auto f = [&shape](const auto& e) noexcept { return e.stepper_begin(shape); };

return build_stepper(f, std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

template <class S>

inline auto xfunction<F, CT...>::stepper_end(const S& shape, layout_type l) const noexcept -> const_stepper

{

auto f = [&shape, l](const auto& e) noexcept { return e.stepper_end(shape, l); };

return build_stepper(f, std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::data_element(size_type i) const -> const_reference

{

return data_element_impl(std::make_index_sequence<sizeof...(CT)>(), i);

}

template <class F, class... CT>

template <class UT, class>

inline xfunction<F, CT...>::operator value_type() const

{

return operator()();

}

template <class F, class... CT>

template <class align, class requested_type, std::size_t N>

inline auto xfunction<F, CT...>::load_simd(size_type i) const

-> simd_return_type<requested_type>

{

return load_simd_impl<align, requested_type, N>(std::make_index_sequence<sizeof...(CT)>(), i);

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::arguments() const noexcept -> const tuple_type&

{

return m_e;

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::functor() const noexcept -> const functor_type&

{

return m_f;

}

template <class F, class... CT>

template <std::size_t... I>

inline layout_type xfunction<F, CT...>::layout_impl(std::index_sequence<I...>) const noexcept

{

return compute_layout(std::get<I>(m_e).layout()...);

}

template <class F, class... CT>

template <std::size_t... I, class... Args>

inline auto xfunction<F, CT...>::access_impl(std::index_sequence<I...>, Args... args) const -> const_reference

{

XTENSOR_TRY(check_index(shape(), args...));

XTENSOR_CHECK_DIMENSION(shape(), args...);

return m_f(std::get<I>(m_e)(args...)...);

}

template <class F, class... CT>

template <std::size_t... I, class... Args>

inline auto xfunction<F, CT...>::unchecked_impl(std::index_sequence<I...>, Args... args) const -> const_reference

{

return m_f(std::get<I>(m_e).unchecked(args...)...);

}

template <class F, class... CT>

template <std::size_t... I, class It>

inline auto xfunction<F, CT...>::element_access_impl(std::index_sequence<I...>, It first, It last) const -> const_reference

{

XTENSOR_TRY(check_element_index(shape(), first, last));

return m_f((std::get<I>(m_e).element(first, last))...);

}

template <class F, class... CT>

template <std::size_t... I>

inline auto xfunction<F, CT...>::data_element_impl(std::index_sequence<I...>, size_type i) const -> const_reference

{

return m_f((std::get<I>(m_e).data_element(i))...);

}

template <class F, class... CT>

template <class align, class requested_type, std::size_t N, std::size_t... I>

inline auto xfunction<F, CT...>::load_simd_impl(std::index_sequence<I...>, size_type i) const

{

return m_f.simd_apply((std::get<I>(m_e)

.template load_simd<align, requested_type>(i))...);

}

template <class F, class... CT>

template <class Func, std::size_t... I>

inline auto xfunction<F, CT...>::build_stepper(Func&& f, std::index_sequence<I...>) const noexcept -> const_stepper

{

return const_stepper(this, f(std::get<I>(m_e))...);

}

template <class F, class... CT>

template <class Func, std::size_t... I>

inline auto xfunction<F, CT...>::build_iterator(Func&& f, std::index_sequence<I...>) const noexcept

{

return const_storage_iterator(this, f(std::get<I>(m_e))...);

}

template <class F, class... CT>

inline auto xfunction<F, CT...>::compute_dimension() const noexcept -> size_type

{

auto func = [](size_type d, auto&& e) noexcept { return (std::max)(d, e.dimension()); };

return accumulate(func, size_type(0), m_e);

}

/*************************************

template <class F, class... CT>

template <class... It>

inline xfunction_iterator<F, CT...>::xfunction_iterator(const xfunction_type* func, It&&... it) noexcept

: p_f(func), m_it(std::forward<It>(it)...)

{

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator++() -> self_type&

{

auto f = [](auto& it) { ++it; };

for_each(f, m_it);

return *this;

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator--() -> self_type&

{

auto f = [](auto& it) { return --it; };

for_each(f, m_it);

return *this;

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator+=(difference_type n) -> self_type&

{

auto f = [n](auto& it) { it += n; };

for_each(f, m_it);

return *this;

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator-=(difference_type n) -> self_type&

{

auto f = [n](auto& it) { it -= n; };

for_each(f, m_it);

return *this;

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator-(const self_type& rhs) const -> difference_type

{

return tuple_max_diff(std::make_index_sequence<sizeof...(CT)>(), m_it, rhs.m_it);

}

template <class F, class... CT>

inline auto xfunction_iterator<F, CT...>::operator*() const -> reference

{

return deref_impl(std::make_index_sequence<sizeof...(CT)>());

}

template <class F, class... CT>

inline bool xfunction_iterator<F, CT...>::equal(const self_type& rhs) const

{

// Optimization: no need to compare each subiterator since they all

// are incremented decremented together.

constexpr std::size_t temp = xtl::mpl::find_if<is_not_xdummy_iterator, data_type>::value;

constexpr std::size_t index = (temp == std::tuple_size<data_type>::value) ? 0 : temp;

return std::get<index>(m_it) == std::get<index>(rhs.m_it);

}

template <class F, class... CT>

inline bool xfunction_iterator<F, CT...>::less_than(const self_type& rhs) const

{

// Optimization: no need to compare each subiterator since they all

// are incremented decremented together.

constexpr std::size_t temp = xtl::mpl::find_if<is_not_xdummy_iterator, data_type>::value;

constexpr std::size_t index = (temp == std::tuple_size<data_type>::value) ? 0 : temp;

return std::get<index>(m_it) < std::get<index>(rhs.m_it);

}

template <class F, class... CT>

template <std::size_t... I>

inline auto xfunction_iterator<F, CT...>::deref_impl(std::index_sequence<I...>) const -> reference

{

return (p_f->m_f)(*std::get<I>(m_it)...);

}

template <class F, class... CT>

template <std::size_t... I>

inline auto xfunction_iterator<F, CT...>::tuple_max_diff(std::index_sequence<I...>,

const data_type& lhs,

const data_type& rhs) const -> difference_type

{

auto diff = std::make_tuple((std::get<I>(lhs) - std::get<I>(rhs))...);

auto func = [](difference_type n, auto&& v) { return (std::max)(n, v); };

return accumulate(func, difference_type(0), diff);

}

template <class F, class... CT>

inline bool operator==(const xfunction_iterator<F, CT...>& it1,

const xfunction_iterator<F, CT...>& it2)

{

return it1.equal(it2);

}

template <class F, class... CT>

inline bool operator<(const xfunction_iterator<F, CT...>& it1,

const xfunction_iterator<F, CT...>& it2)

{

return it1.less_than(it2);

}

/************************************

template <class F, class... CT>

template <class... St>

inline xfunction_stepper<F, CT...>::xfunction_stepper(const xfunction_type* func, St&&... st) noexcept

: p_f(func), m_st(std::forward<St>(st)...)

{

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::step(size_type dim)

{

auto f = [dim](auto& st) { st.step(dim); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::step_back(size_type dim)

{

auto f = [dim](auto& st) { st.step_back(dim); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::step(size_type dim, size_type n)

{

auto f = [dim, n](auto& st) { st.step(dim, n); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::step_back(size_type dim, size_type n)

{

auto f = [dim, n](auto& st) { st.step_back(dim, n); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::reset(size_type dim)

{

auto f = [dim](auto& st) { st.reset(dim); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::reset_back(size_type dim)

{

auto f = [dim](auto& st) { st.reset_back(dim); };

for_each(f, m_st);

}

template <class F, class... CT>

inline void xfunction_stepper<F, CT...>::to_begin()

{

auto f = [](auto& st) { st.to_begin(); };

for_each(f, m_st);