#ifndef PYTHONIC_INCLUDE_TYPES_NUMPY_GEXPR_HPP
#define PYTHONIC_INCLUDE_TYPES_NUMPY_GEXPR_HPP

#include "pythonic/include/utils/meta.hpp"
#include "pythonic/include/utils/array_helper.hpp"

PYTHONIC_NS_BEGIN

namespace types
{

  /* helper to count new axis
   */
  template <class... T>
  struct count_new_axis;
  template <>
  struct count_new_axis<> {
    static constexpr size_t value = 0;
  };

  template <>
  struct count_new_axis<types::none_type> {
    static constexpr size_t value = 1;
  };

  template <class T>
  struct count_new_axis<T> {
    static constexpr size_t value = 0;
  };

  template <class T0, class... T>
  struct count_new_axis<T0, T...> {
    static constexpr size_t value =
        count_new_axis<T0>::value + count_new_axis<T...>::value;
  };

  /* helper to turn a new axis into a slice
   */
  template <class T>
  struct to_slice {
    using type = T;
    static constexpr bool is_new_axis = false;
    T operator()(T value);
  };

  template <>
  struct to_slice<none_type> {
    using type = fast_contiguous_slice;
    static constexpr bool is_new_axis = true;
    fast_contiguous_slice operator()(none_type);
  };

  template <class T>
  struct to_normalized_slice {
    using type = T;
    T operator()(T value);
  };

  template <>
  struct to_normalized_slice<none_type> {
    using type = contiguous_normalized_slice;
    contiguous_normalized_slice operator()(none_type);
  };

  /* helper to build a new shape out of a shape and a slice with new axis
   */
  template <size_t N, class pS, class IsNewAxis>
  auto make_reshape(pS const &shape, IsNewAxis is_new_axis)
      -> decltype(sutils::copy_new_axis<pS::value + N>(shape, is_new_axis));

  /* helper to build an extended slice aka numpy_gexpr out of a subscript
   */
  template <size_t C>
  struct extended_slice {
    template <class E, class... S>
    auto operator()(E &&expr, S const &... s)
        -> decltype(std::forward<E>(expr).reshape(make_reshape<C>(
            expr,
            std::tuple<
                std::integral_constant<bool, to_slice<S>::is_new_axis>...>()))(
            to_slice<S>{}(s)...))

    {
      return std::forward<E>(expr).reshape(make_reshape<C>(
          expr,
          std::tuple<
              std::integral_constant<bool, to_slice<S>::is_new_axis>...>()))(
          to_slice<S>{}(s)...);
    }
  };

  template <>
  struct extended_slice<0> {
    template <class E, class... S>
    auto operator()(E &&expr, long const &s0, S const &... s) ->
        typename std::enable_if<
            utils::all_of<std::is_integral<S>::value...>::value,
            decltype(std::forward<E>(expr)[types::make_tuple(s0, s...)])>::type
    {
      return std::forward<E>(expr)[types::make_tuple(s0, s...)];
    }
    template <class E, class... S>
    auto operator()(E &&expr, long const &s0, S const &... s) ->
        typename std::enable_if<
            !utils::all_of<std::is_integral<S>::value...>::value,
            decltype(std::forward<E>(expr)[s0](s...))>::type
    {
      return std::forward<E>(expr)[s0](s...);
    }

    template <class E, class... S, size_t... Is>
    numpy_gexpr<typename std::decay<E>::type, normalize_t<S>...>
    fwd(E &&expr, std::tuple<S...> const &s, utils::index_sequence<Is...>)
    {
      return {std::forward<E>(expr),
              std::get<Is>(s).normalize(expr.template shape<Is>())...};
    }

    template <class E, class Sp, class... S>
    typename std::enable_if<
        is_slice<Sp>::value,
        numpy_gexpr<E, normalize_t<Sp>, normalize_t<S>...>>::type
    operator()(E &&expr, Sp const &s0, S const &... s)
    {
      return make_gexpr(std::forward<E>(expr), s0, s...);
    }

    template <class E, class F, class... S>
    typename std::enable_if<
        !is_slice<F>::value,
        numpy_gexpr<ndarray<typename std::decay<E>::type::dtype,
                            array<long, std::decay<E>::type::value>>,
                    contiguous_normalized_slice, normalize_t<S>...>>::type
    operator()(E &&expr, F const &s0, S const &... s)
    {
      return numpy_vexpr<ndarray<typename std::decay<E>::type::dtype,
                                 array<long, std::decay<E>::type::value>>,
                         F>{std::forward<E>(expr), s0}(
          fast_contiguous_slice(none_type{}, none_type{}), s...);
    }
  };

  /* Meta-Function to count the number of slices in a type list
   */
  template <class... Types>
  struct count_long;

  template <>
  struct count_long<long> {
    static constexpr size_t value = 1;
  };

  template <>
  struct count_long<normalized_slice> {
    static constexpr size_t value = 0;
  };

  template <>
  struct count_long<contiguous_normalized_slice> {
    static constexpr size_t value = 0;
  };

  template <class T, class... Types>
  struct count_long<T, Types...> {
    static constexpr size_t value =
        count_long<T>::value + count_long<Types...>::value;
  };

  template <>
  struct count_long<> {
    static constexpr size_t value = 0;
  };

  /* helper to get the type of the nth element of an array
   */
  template <class T, size_t N>
  struct nth_value_type {
    using type = typename nth_value_type<typename T::value_type, N - 1>::type;
  };

  template <class T>
  struct nth_value_type<T, 0> {
    using type = T;
  };

  /* helper that yields true if the first slice of a pack is a contiguous
   * slice
   */
  template <class... S>
  struct is_contiguous {
    static const bool value = false;
  };

  template <class... S>
  struct is_contiguous<contiguous_normalized_slice, S...> {
    static const bool value = true;
  };

  /* numpy_gexpr factory
   *
   * replaces the constructor, in order to properly merge gexpr composition
   *into a single gexpr
   */
  namespace details
  {

    template <class T, class Ts, size_t... Is>
    std::tuple<T, typename std::tuple_element<Is, Ts>::type...>
    tuple_push_head(T const &val, Ts const &vals, utils::index_sequence<Is...>)
    {
      return std::tuple<T, typename std::tuple_element<Is, Ts>::type...>{
          val, std::get<Is>(vals)...};
    }

    template <class T, class Ts>
    auto tuple_push_head(T const &val, Ts const &vals)
        -> decltype(tuple_push_head(
            val, vals,
            utils::make_index_sequence<std::tuple_size<Ts>::value>()))
    {
      return tuple_push_head(
          val, vals, utils::make_index_sequence<std::tuple_size<Ts>::value>());
    }

    // this struct is specialized for every type combination && takes care of
    // the slice merge
    template <class T, class Tp>
    struct merge_gexpr;

    template <>
    struct merge_gexpr<std::tuple<>, std::tuple<>> {
      template <size_t I, class S>
      std::tuple<> run(S const &, std::tuple<> const &t0, std::tuple<> const &);
    };

    template <class... T0>
    struct merge_gexpr<std::tuple<T0...>, std::tuple<>> {
      template <size_t I, class S>
      std::tuple<T0...> run(S const &, std::tuple<T0...> const &t0,
                            std::tuple<>);
      static_assert(
          utils::all_of<std::is_same<T0, normalize_t<T0>>::value...>::value,
          "all slices are normalized");
    };

    template <class... T1>
    struct merge_gexpr<std::tuple<>, std::tuple<T1...>> {
      template <size_t I, class S>
      std::tuple<normalize_t<T1>...> run(S const &, std::tuple<>,
                                         std::tuple<T1...> const &t1);
    };

    template <class S0, class... T0, class S1, class... T1>
    struct merge_gexpr<std::tuple<S0, T0...>, std::tuple<S1, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<S0, T0...> const &t0,
               std::tuple<S1, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t0) * std::get<0>(t1),
              merge_gexpr<std::tuple<T0...>, std::tuple<T1...>>{}
                  .template run<I + 1>(s, tuple_tail(t0), tuple_tail(t1))))
      {
        return tuple_push_head(
            std::get<0>(t0) * std::get<0>(t1),
            merge_gexpr<std::tuple<T0...>, std::tuple<T1...>>{}
                .template run<I + 1>(s, tuple_tail(t0), tuple_tail(t1)));
      }
      static_assert(
          std::is_same<decltype(std::declval<S0>() * std::declval<S1>()),
                       normalize_t<decltype(std::declval<S0>() *
                                            std::declval<S1>())>>::value,
          "all slices are normalized");
    };

    template <class S0, class... T0, class... T1>
    struct merge_gexpr<std::tuple<S0, T0...>, std::tuple<none_type, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<S0, T0...> const &t0,
               std::tuple<none_type, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t1),
              merge_gexpr<std::tuple<S0, T0...>, std::tuple<T1...>>{}
                  .template run<I + 1>(s, t0, tuple_tail(t1))))
      {
        return tuple_push_head(
            std::get<0>(t1),
            merge_gexpr<std::tuple<S0, T0...>, std::tuple<T1...>>{}
                .template run<I + 1>(s, t0, tuple_tail(t1)));
      }
    };

    template <class... T0, class S1, class... T1>
    struct merge_gexpr<std::tuple<long, T0...>, std::tuple<S1, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<long, T0...> const &t0,
               std::tuple<S1, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t0),
              merge_gexpr<std::tuple<T0...>, std::tuple<S1, T1...>>{}
                  .template run<I>(s, tuple_tail(t0), t1)))
      {
        return tuple_push_head(
            std::get<0>(t0),
            merge_gexpr<std::tuple<T0...>, std::tuple<S1, T1...>>{}
                .template run<I>(s, tuple_tail(t0), t1));
      }
    };
    template <class... T0, class... T1>
    struct merge_gexpr<std::tuple<long, T0...>, std::tuple<none_type, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<long, T0...> const &t0,
               std::tuple<none_type, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t0),
              merge_gexpr<std::tuple<T0...>, std::tuple<none_type, T1...>>{}
                  .template run<I>(s, tuple_tail(t0), t1)))
      {
        return tuple_push_head(
            std::get<0>(t0),
            merge_gexpr<std::tuple<T0...>, std::tuple<none_type, T1...>>{}
                .template run<I>(s, tuple_tail(t0), t1));
      }
    };

    template <class S0, class... T0, class... T1>
    struct merge_gexpr<std::tuple<S0, T0...>, std::tuple<long, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<S0, T0...> const &t0,
               std::tuple<long, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t1) * std::get<0>(t0).step + std::get<0>(t0).lower,
              merge_gexpr<std::tuple<T0...>, std::tuple<T1...>>{}
                  .template run<I + 1>(s, tuple_tail(t0), tuple_tail(t1))))
      {
        return tuple_push_head(
            std::get<0>(t1) * std::get<0>(t0).step + std::get<0>(t0).lower,
            merge_gexpr<std::tuple<T0...>, std::tuple<T1...>>{}
                .template run<I + 1>(s, tuple_tail(t0), tuple_tail(t1)));
      }
    };

    template <class... T0, class... T1>
    struct merge_gexpr<std::tuple<long, T0...>, std::tuple<long, T1...>> {
      template <size_t I, class S>
      auto run(S const &s, std::tuple<long, T0...> const &t0,
               std::tuple<long, T1...> const &t1)
          -> decltype(tuple_push_head(
              std::get<0>(t0),
              merge_gexpr<std::tuple<T0...>, std::tuple<long, T1...>>{}
                  .template run<I>(s, tuple_tail(t0), t1)))
      {
        return tuple_push_head(
            std::get<0>(t0),
            merge_gexpr<std::tuple<T0...>, std::tuple<long, T1...>>{}
                .template run<I>(s, tuple_tail(t0), t1));
      }
    };

    template <class Arg, class... Sp>
    typename std::enable_if<count_new_axis<Sp...>::value == 0,
                            numpy_gexpr<Arg, Sp...>>::type
    _make_gexpr(Arg arg, std::tuple<Sp...> const &t);

    template <class Arg, class S, size_t... Is>
    numpy_gexpr<Arg, typename to_normalized_slice<
                         typename std::tuple_element<Is, S>::type>::type...>
    _make_gexpr_helper(Arg arg, S const &s, utils::index_sequence<Is...>);

    template <class Arg, class... Sp>
    auto _make_gexpr(Arg arg, std::tuple<Sp...> const &s) ->
        typename std::enable_if<
            count_new_axis<Sp...>::value != 0,
            decltype(_make_gexpr_helper(
                arg.reshape(make_reshape<count_new_axis<Sp...>::value>(
                    arg, std::tuple<std::integral_constant<
                             bool, to_slice<Sp>::is_new_axis>...>())),
                s, utils::make_index_sequence<sizeof...(Sp)>()))>::type;

    template <class Arg, class... S>
    struct make_gexpr {
      template <size_t... Is>
      numpy_gexpr<Arg, normalize_t<S>...>
      operator()(Arg arg, std::tuple<S...>, utils::index_sequence<Is...>);
      numpy_gexpr<Arg, normalize_t<S>...> operator()(Arg arg, S const &... s);
    };

    // this specialization is in charge of merging gexpr
    template <class Arg, class... S, class... Sp>
    struct make_gexpr<numpy_gexpr<Arg, S...> const &, Sp...> {
      auto operator()(numpy_gexpr<Arg, S...> const &arg, Sp const &... s)
          -> decltype(
              _make_gexpr(std::declval<Arg>(),
                          merge_gexpr<std::tuple<S...>, std::tuple<Sp...>>{}
                              .template run<0>(arg, std::tuple<S...>(),
                                               std::tuple<Sp...>())))
      {
        return _make_gexpr(
            arg.arg,
            merge_gexpr<std::tuple<S...>, std::tuple<Sp...>>{}.template run<0>(
                arg, arg.slices, std::make_tuple(s...)));
      }
    };
  }

  template <class Arg, class... S>
  auto make_gexpr(Arg &&arg, S const &... s)
      -> decltype(details::make_gexpr<Arg, S...>{}(std::forward<Arg>(arg),
                                                   s...));

  /* type-based compile time overlapping detection: detect if a type may
   *overlap with another
   * the goal is to detect whether the following operation
   *
   * a[...] = b
   *
   * requires a copy.
   *
   * It requires a copy if b = a[...], as in
   *
   * a[1:] = a[:-1]
   *
   * because this is *!* equivalent to for i in range(0, n-1): a[i+1] = a[i]
   *
   * to avoid the copy, we rely on the lhs type
   */

  template <class E>
  struct may_overlap_gexpr : std::integral_constant<bool, !is_dtype<E>::value> {
  };

  template <class T0, class T1>
  struct may_overlap_gexpr<broadcast<T0, T1>> : std::false_type {
  };

  template <class E>
  struct may_overlap_gexpr<broadcasted<E>> : std::false_type {
  };

  template <class E>
  struct may_overlap_gexpr<E &> : may_overlap_gexpr<E> {
  };

  template <class E>
  struct may_overlap_gexpr<E const &> : may_overlap_gexpr<E> {
  };

  template <class T, class pS>
  struct may_overlap_gexpr<ndarray<T, pS>> : std::false_type {
  };

  template <class E>
  struct may_overlap_gexpr<numpy_iexpr<E>> : may_overlap_gexpr<E> {
  };

  template <class E>
  struct may_overlap_gexpr<numpy_texpr<E>> : may_overlap_gexpr<E> {
  };

  template <class E>
  struct may_overlap_gexpr<list<E>>
      : std::integral_constant<bool, !is_dtype<E>::value> {
  };

  template <class E, size_t N, class V>
  struct may_overlap_gexpr<array_base<E, N, V>> : may_overlap_gexpr<E> {
  };

  template <class Op, class... Args>
  struct may_overlap_gexpr<numpy_expr<Op, Args...>>
      : utils::any_of<may_overlap_gexpr<Args>::value...> {
  };

  template <class OpS, class pS, class... S>
  struct gexpr_shape;

  template <class... Tys, class... oTys>
  struct gexpr_shape<pshape<Tys...>, pshape<oTys...>> {
    using type = pshape<Tys..., oTys...>;
  };

  template <class... Tys>
  struct gexpr_shape<pshape<Tys...>, array<long, 0>> {
    using type = pshape<Tys...>;
  };

  template <class... Tys, size_t N>
  struct gexpr_shape<pshape<Tys...>, array<long, N>>
      : gexpr_shape<pshape<Tys..., long>, array<long, N - 1>> {
  };

  template <class... Tys, class... oTys, class... S>
  struct gexpr_shape<pshape<Tys...>,
                     pshape<std::integral_constant<long, 1>, oTys...>,
                     contiguous_normalized_slice, S...>
      : gexpr_shape<pshape<Tys..., std::integral_constant<long, 1>>,
                    pshape<oTys...>, S...> {
  };
  template <class... Tys, class... oTys, class... S>
  struct gexpr_shape<pshape<Tys...>,
                     pshape<std::integral_constant<long, 1>, oTys...>,
                     normalized_slice, S...>
      : gexpr_shape<pshape<Tys..., std::integral_constant<long, 1>>,
                    pshape<oTys...>, S...> {
  };

  template <class... Tys, class oT, class... oTys, class... S>
  struct gexpr_shape<pshape<Tys...>, pshape<oT, oTys...>, long, S...>
      : gexpr_shape<pshape<Tys...>, pshape<oTys...>, S...> {
  };
  template <class... Tys, class oT, class... oTys, class cS, class... S>
  struct gexpr_shape<pshape<Tys...>, pshape<oT, oTys...>, cS, S...>
      : gexpr_shape<pshape<Tys..., long>, pshape<oTys...>, S...> {
  };
  template <class... Tys, size_t N, class... S>
  struct gexpr_shape<pshape<Tys...>, array<long, N>, long, S...>
      : gexpr_shape<pshape<Tys...>, array<long, N - 1>, S...> {
  };
  template <class... Tys, size_t N, class cS, class... S>
  struct gexpr_shape<pshape<Tys...>, array<long, N>, cS, S...>
      : gexpr_shape<pshape<Tys..., long>, array<long, N - 1>, S...> {
  };

  template <class pS, class... S>
  using gexpr_shape_t = typename gexpr_shape<pshape<>, pS, S...>::type;

  /* Expression template for numpy expressions - extended slicing operators
   */
  template <class Arg, class... S>
  struct numpy_gexpr {
    static_assert(
        utils::all_of<std::is_same<S, normalize_t<S>>::value...>::value,
        "all slices are normalized");
    static_assert(
        utils::all_of<(std::is_same<S, long>::value ||
                       std::is_same<S, contiguous_normalized_slice>::value ||
                       std::is_same<S, normalized_slice>::value)...>::value,
        "all slices are valid");
    static_assert(std::decay<Arg>::type::value >= sizeof...(S),
                  "slicing respects array shape");

    // numpy_gexpr is a wrapper for extended sliced array around a numpy
    // expression.
    // It contains compacted sorted slices value in lower, step && upper is
    // the same as shape.
    // indices for long index are store in the indices array.
    // position for slice and long value in the extended slice can be found
    // through the S... template
    // && compacted values as we know that first S is a slice.

    static_assert(
        utils::all_of<
            std::is_same<S, typename std::decay<S>::type>::value...>::value,
        "no modifiers on slices");

    using dtype = typename std::remove_reference<Arg>::type::dtype;
    static constexpr size_t value =
        std::remove_reference<Arg>::type::value - count_long<S...>::value;

    // It is not possible to vectorize everything. We only vectorize if the
    // last dimension is contiguous, which happens if
    // 1. Arg is an ndarray (this is too strict)
    // 2. the size of the gexpr is lower than the dim of arg, || it's the
    // same, but the last slice is contiguous
    static const bool is_vectorizable =
        std::remove_reference<Arg>::type::is_vectorizable &&
        (sizeof...(S) < std::remove_reference<Arg>::type::value ||
         std::is_same<contiguous_normalized_slice,
                      typename std::tuple_element<
                          sizeof...(S)-1, std::tuple<S...>>::type>::value);
    static const bool is_strided =
        std::remove_reference<Arg>::type::is_strided ||
        (((sizeof...(S)-count_long<S...>::value) == value) &&
         !std::is_same<contiguous_normalized_slice,
                       typename std::tuple_element<
                           sizeof...(S)-1, std::tuple<S...>>::type>::value);

    using value_type = typename std::decay<decltype(
        numpy_iexpr_helper<value>::get(std::declval<numpy_gexpr>(), 1))>::type;

    using iterator =
        typename std::conditional<is_strided || value != 1,
                                  nditerator<numpy_gexpr>, dtype *>::type;
    using const_iterator =
        typename std::conditional<is_strided || value != 1,
                                  const_nditerator<numpy_gexpr>,
                                  dtype const *>::type;

    typename std::remove_cv<Arg>::type arg;

    std::tuple<S...> slices;

    using shape_t =
        gexpr_shape_t<typename std::remove_reference<Arg>::type::shape_t, S...>;

    shape_t _shape;
    dtype *buffer;
    array<long, value> _strides;
    template <size_t I>
    auto shape() const -> decltype(std::get<I>(_shape))
    {
      return std::get<I>(_shape);
    }
    template <size_t I>
    auto strides() const -> decltype(std::get<I>(_strides))
    {
      return std::get<I>(_strides);
    }

    numpy_gexpr();
    numpy_gexpr(numpy_gexpr const &) = default;
    numpy_gexpr(numpy_gexpr &&) = default;

    template <class Argp> // ! using the default one, to make it possible to
    // accept reference && non reference version of
    // Argp
    numpy_gexpr(numpy_gexpr<Argp, S...> const &other);

    template <size_t J, class Slice>
    typename std::enable_if<
        std::is_same<Slice, normalized_slice>::value ||
            std::is_same<Slice, contiguous_normalized_slice>::value,
        void>::type
    init_shape(Slice const &s, utils::int_<1>, utils::int_<J>);

    template <size_t I, size_t J, class Slice>
    typename std::enable_if<
        std::is_same<Slice, normalized_slice>::value ||
            std::is_same<Slice, contiguous_normalized_slice>::value,
        void>::type
    init_shape(Slice const &s, utils::int_<I>, utils::int_<J>);

    template <size_t J>
    void init_shape(long cs, utils::int_<1>, utils::int_<J>);

    template <size_t I, size_t J>
    void init_shape(long cs, utils::int_<I>, utils::int_<J>);

    // private because we must use the make_gexpr factory to create a gexpr
  private:
    template <class _Arg, class... _other_classes>
    friend struct details::make_gexpr;

    friend struct array_base_slicer;
    template <class _Arg, class... _other_classes>
    friend
        typename std::enable_if<count_new_axis<_other_classes...>::value == 0,
                                numpy_gexpr<_Arg, _other_classes...>>::type
        details::_make_gexpr(_Arg arg, std::tuple<_other_classes...> const &t);
    template <class _Arg, class _other_classes, size_t... Is>
    friend numpy_gexpr<_Arg,
                       typename to_normalized_slice<typename std::tuple_element<
                           Is, _other_classes>::type>::type...>
    details::_make_gexpr_helper(_Arg arg, _other_classes const &s,
                                utils::index_sequence<Is...>);

    template <size_t C>
    friend struct extended_slice;

#ifdef ENABLE_PYTHON_MODULE
    template <typename T>
    friend struct pythonic::from_python;
#endif

    // When we create a new numpy_gexpr, we deduce step, lower && shape from
    // slices
    // && indices from long value.
    // Also, last shape information are set from origin array like in :
    // >>> a = numpy.arange(2*3*4).reshape(2,3,4)
    // >>> a[:, 1]
    // the last dimension (4) is missing from slice information
    // Finally, if origin expression was already sliced, lower bound && step
    // have to
    // be increased
    numpy_gexpr(Arg const &arg, std::tuple<S const &...> const &values);
    numpy_gexpr(Arg const &arg, S const &... s);

  public:
    template <class Argp, class... Sp>
    numpy_gexpr(numpy_gexpr<Argp, Sp...> const &expr, Arg arg);

    template <class G>
    numpy_gexpr(G const &expr, Arg &&arg);
    template <class pS>
    ndarray<dtype, pS> reshape(pS const &shape) const
    {
      return copy().reshape(shape);
    }

    template <class E>
    typename std::enable_if<may_overlap_gexpr<E>::value, numpy_gexpr &>::type
    _copy(E const &expr);

    template <class E>
    typename std::enable_if<!may_overlap_gexpr<E>::value, numpy_gexpr &>::type
    _copy(E const &expr);

    template <class E>
    numpy_gexpr &operator=(E const &expr);

    numpy_gexpr &operator=(numpy_gexpr const &expr);

    template <class Argp>
    numpy_gexpr &operator=(numpy_gexpr<Argp, S...> const &expr);

    template <class Op, class E>
    typename std::enable_if<may_overlap_gexpr<E>::value, numpy_gexpr &>::type
    update_(E const &expr);

    template <class Op, class E>
    typename std::enable_if<!may_overlap_gexpr<E>::value, numpy_gexpr &>::type
    update_(E const &expr);

    template <class E>
    numpy_gexpr &operator+=(E const &expr);

    numpy_gexpr &operator+=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator-=(E const &expr);

    numpy_gexpr &operator-=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator*=(E const &expr);

    numpy_gexpr &operator*=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator/=(E const &expr);

    numpy_gexpr &operator/=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator|=(E const &expr);

    numpy_gexpr &operator|=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator&=(E const &expr);

    numpy_gexpr &operator&=(numpy_gexpr const &expr);

    template <class E>
    numpy_gexpr &operator^=(E const &expr);

    numpy_gexpr &operator^=(numpy_gexpr const &expr);

    const_iterator begin() const;
    const_iterator end() const;

    iterator begin();
    iterator end();

    auto fast(long i) const
        & -> decltype(numpy_iexpr_helper<value>::get(*this, i))
    {
      return numpy_iexpr_helper<value>::get(*this, i);
    }

    auto fast(long i) & -> decltype(numpy_iexpr_helper<value>::get(*this, i))
    {
      return numpy_iexpr_helper<value>::get(*this, i);
    }

    template <class E, class... Indices>
    void store(E elt, Indices... indices)
    {
      static_assert(is_dtype<E>::value, "valid store");
      *(buffer + noffset<value>{}(*this, array<long, value>{{indices...}})) =
          static_cast<E>(elt);
    }
    template <class... Indices>
    dtype load(Indices... indices) const
    {
      return *(buffer +
               noffset<value>{}(*this, array<long, value>{{indices...}}));
    }
    template <class Op, class E, class... Indices>
    void update(E elt, Indices... indices) const
    {
      static_assert(is_dtype<E>::value, "valid store");
      Op{}(
          *(buffer + noffset<value>{}(*this, array<long, value>{{indices...}})),
          static_cast<E>(elt));
    }

#ifdef USE_XSIMD
    using simd_iterator = const_simd_nditerator<numpy_gexpr>;
    using simd_iterator_nobroadcast = simd_iterator;
    template <class vectorizer>
    simd_iterator vbegin(vectorizer) const;
    template <class vectorizer>
    simd_iterator vend(vectorizer) const;
#endif

    template <class... Sp>
    auto operator()(Sp const &... s) const -> decltype(make_gexpr(*this, s...));

    template <class Sp>
    auto operator[](Sp const &s) const -> typename std::enable_if<
        is_slice<Sp>::value, decltype(make_gexpr(*this, (s.lower, s)))>::type;

    template <size_t M>
    auto fast(array<long, M> const &indices) const
        & -> decltype(nget<M - 1>().fast(*this, indices));

    template <size_t M>
        auto fast(array<long, M> const &indices) &&
        -> decltype(nget<M - 1>().fast(std::move(*this), indices));

    template <size_t M>
    auto operator[](array<long, M> const &indices) const
        & -> decltype(nget<M - 1>()(*this, indices));

    template <size_t M>
        auto operator[](array<long, M> const &indices) &&
        -> decltype(nget<M - 1>()(std::move(*this), indices));

    template <class F> // indexing through an array of indices -- a view
    typename std::enable_if<is_numexpr_arg<F>::value &&
                                !is_array_index<F>::value &&
                                !std::is_same<bool, typename F::dtype>::value,
                            numpy_vexpr<numpy_gexpr, F>>::type
    operator[](F const &filter) const
    {
      return {*this, filter};
    }
    template <class F> // indexing through an array of indices -- a view
    typename std::enable_if<is_numexpr_arg<F>::value &&
                                !is_array_index<F>::value &&
                                !std::is_same<bool, typename F::dtype>::value,
                            numpy_vexpr<numpy_gexpr, F>>::type
    fast(F const &filter) const
    {
      return {*this, filter};
    }

    template <class F>
    typename std::enable_if<
        is_numexpr_arg<F>::value &&
            std::is_same<bool, typename F::dtype>::value,
        numpy_vexpr<numpy_gexpr, ndarray<long, pshape<long>>>>::type
    fast(F const &filter) const;

    template <class F>
    typename std::enable_if<
        is_numexpr_arg<F>::value &&
            std::is_same<bool, typename F::dtype>::value,
        numpy_vexpr<numpy_gexpr, ndarray<long, pshape<long>>>>::type
    operator[](F const &filter) const;
    auto operator[](long i) const -> decltype(this->fast(i));
    auto operator[](long i) -> decltype(this->fast(i));

    // template <class... Sp>
    // auto operator()(long i, Sp const &... s) const
    //    -> decltype((*this)[i](s...));

    explicit operator bool() const;

    long flat_size() const;
    long size() const;
    ndarray<dtype, shape_t> copy() const
    {
      return {*this};
    }
  };
}

template <class E, class... S>
struct assignable_noescape<types::numpy_gexpr<E, S...>> {
  using type = types::numpy_gexpr<E, S...>;
};

template <class T, class pS, class... S>
struct assignable<types::numpy_gexpr<types::ndarray<T, pS> const &, S...>> {
  using type = types::numpy_gexpr<types::ndarray<T, pS>, S...>;
};

template <class T, class pS, class... S>
struct assignable<types::numpy_gexpr<types::ndarray<T, pS> &, S...>> {
  using type = types::numpy_gexpr<types::ndarray<T, pS>, S...>;
};

template <class Arg, class... S>
struct assignable<types::numpy_gexpr<Arg, S...>> {
  using type = types::numpy_gexpr<typename assignable<Arg>::type, S...>;
};

template <class Arg, class... S>
struct lazy<types::numpy_gexpr<Arg, S...>>
    : assignable<types::numpy_gexpr<Arg, S...>> {
};

PYTHONIC_NS_END
/* type inference stuff  {*/
#include "pythonic/include/types/combined.hpp"

template <class Arg, class... S>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>,
                  pythonic::types::numpy_gexpr<Arg, S...>> {
  using type = pythonic::types::numpy_gexpr<Arg, S...>;
};

template <class Arg, class... S, class Argp, class... Sp>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>,
                  pythonic::types::numpy_gexpr<Argp, Sp...>> {
  using t0 = pythonic::types::numpy_gexpr<Arg, S...>;
  using t1 = pythonic::types::numpy_gexpr<Argp, Sp...>;
  using type =
      pythonic::types::ndarray <
      typename __combined<typename t0::dtype, typename t1::dtype>::type,
        pythonic::types::array < long,
        t0::value < t1::value ? t1::value : t0::value >>
      ;
};

template <class Arg, class... S, class O>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>, O> {
  using type = pythonic::types::numpy_gexpr<Arg, S...>;
};
template <class Arg, class... S, class T>
struct __combined<pythonic::types::list<T>,
                  pythonic::types::numpy_gexpr<Arg, S...>> {
  using type = pythonic::types::list<typename __combined<
      typename pythonic::types::numpy_gexpr<Arg, S...>::value_type, T>::type>;
};
template <class Arg, class... S, class T>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>,
                  pythonic::types::list<T>> {
  using type = pythonic::types::list<typename __combined<
      typename pythonic::types::numpy_gexpr<Arg, S...>::value_type, T>::type>;
};
template <class Arg, class... S>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>,
                  pythonic::types::none_type> {
  using type = pythonic::types::none<pythonic::types::numpy_gexpr<Arg, S...>>;
};

template <class Arg, class... S, class O>
struct __combined<O, pythonic::types::numpy_gexpr<Arg, S...>> {
  using type = pythonic::types::numpy_gexpr<Arg, S...>;
};
template <class Arg, class... S>
struct __combined<pythonic::types::none_type,
                  pythonic::types::numpy_gexpr<Arg, S...>> {
  using type = pythonic::types::none<pythonic::types::numpy_gexpr<Arg, S...>>;
};

/* combined are sorted such that the assigned type comes first */
template <class Arg, class... S, class T, class pS>
struct __combined<pythonic::types::numpy_gexpr<Arg, S...>,
                  pythonic::types::ndarray<T, pS>> {
  using type = pythonic::types::ndarray<T, pS>;
};

template <class Arg, class... S, class T, class pS>
struct __combined<pythonic::types::ndarray<T, pS>,
                  pythonic::types::numpy_gexpr<Arg, S...>> {
  using type = pythonic::types::ndarray<T, pS>;
};

#endif