From 46da5e5f524f3efe584a9ee0c9465461754e5efe Mon Sep 17 00:00:00 2001
From: serge-sans-paille <serge.guelton@telecom-bretagne.eu>
Date: Wed, 29 Mar 2023 22:26:48 +0200
Subject: [PATCH] Detect potential copyto usage and improve numpy.copyto
 support

Detect a[:] = expr pattern and use np.copyto(a, expr) instead.
Use memcpy when possible (through broadcast_copy improvement).
Accept moved reference.
---
 pythran/optimizations/__init__.py             |   1 +
 pythran/optimizations/copyto.py               |  91 ++++++++++++
 pythran/pythonic/include/numpy/copyto.hpp     |  18 ++-
 pythran/pythonic/include/types/list.hpp       |   9 +-
 pythran/pythonic/include/types/ndarray.hpp    |   3 +
 .../include/types/numpy_broadcast.hpp         |   2 +
 pythran/pythonic/include/types/numpy_expr.hpp |   1 +
 .../pythonic/include/types/numpy_gexpr.hpp    |   6 +
 .../pythonic/include/types/numpy_iexpr.hpp    |   5 +
 .../pythonic/include/types/numpy_texpr.hpp    |   1 +
 .../pythonic/include/types/numpy_vexpr.hpp    |   1 +
 pythran/pythonic/include/types/tuple.hpp      |   1 +
 pythran/pythonic/numpy/copyto.hpp             |  56 +++++++-
 pythran/pythonic/types/list.hpp               | 136 +++++++++---------
 pythran/pythonic/types/ndarray.hpp            |   2 +-
 pythran/pythonic/utils/broadcast_copy.hpp     |  44 +++++-
 pythran/pythran.cfg                           |   1 +
 pythran/tests/test_ndarray.py                 |   5 +
 pythran/tests/test_numpy_func3.py             |   5 +
 19 files changed, 306 insertions(+), 82 deletions(-)
 create mode 100644 pythran/optimizations/copyto.py

diff --git a/pythran/optimizations/__init__.py b/pythran/optimizations/__init__.py
index a7a8aec9a9..899e6ce876 100644
--- a/pythran/optimizations/__init__.py
+++ b/pythran/optimizations/__init__.py
@@ -11,6 +11,7 @@
 """
 
 from .constant_folding import ConstantFolding, PartialConstantFolding
+from .copyto import CopyTo
 from .dead_code_elimination import DeadCodeElimination
 from .forward_substitution import ForwardSubstitution, PreInliningForwardSubstitution
 from .iter_transformation import IterTransformation
diff --git a/pythran/optimizations/copyto.py b/pythran/optimizations/copyto.py
new file mode 100644
index 0000000000..33aa79fcc5
--- /dev/null
+++ b/pythran/optimizations/copyto.py
@@ -0,0 +1,91 @@
+""" Replaces a[:] = b by a call to numpy.copyto. """
+
+from pythran.passmanager import Transformation
+from pythran.analyses.ast_matcher import ASTMatcher, AST_any
+from pythran.conversion import mangle
+from pythran.utils import isnum
+
+import gast as ast
+import copy
+
+
+class CopyTo(Transformation):
+
+    """
+    Replaces a[:] = b by a call to numpy.copyto.
+
+    This is a slight extension to numpy.copyto as it assumes it also supports
+    string and list as first argument.
+
+    >>> import gast as ast
+    >>> from pythran import passmanager, backend
+    >>> node = ast.parse('a[:] = b')
+    >>> pm = passmanager.PassManager("test")
+    >>> _, node = pm.apply(CopyTo, node)
+    >>> print(pm.dump(backend.Python, node))
+    import numpy as __pythran_import_numpy
+    __pythran_import_numpy.copyto(a, b)
+    >>> node = ast.parse('a[:] = b[:]')
+    >>> pm = passmanager.PassManager("test")
+    >>> _, node = pm.apply(CopyTo, node)
+    >>> print(pm.dump(backend.Python, node))
+    import numpy as __pythran_import_numpy
+    __pythran_import_numpy.copyto(a, b)
+    """
+
+    def isNone(self, node):
+        if node is None:
+            return True
+        return isinstance(node, ast.Constant) and node.value is None
+
+    def is_full_slice(self, node):
+        # FIXME: could accept a call to len for node.upper
+        return (
+                isinstance(node, ast.Slice) and
+                (node.lower == 0 or self.isNone(node.lower)) and
+                (self.isNone(node.upper)) and
+                (self.isNone(node.step) or node.step == 1)
+                )
+
+    def is_fully_sliced(self, node):
+        if not isinstance(node, ast.Subscript):
+            return False
+        if not isinstance(node.value, ast.Name):
+            return False
+        if self.is_full_slice(node.slice):
+            return True
+        elif isinstance(node.slice, ast.Tuple):
+            return all(self.is_full_slice(elt) for elt in node.slice.elts)
+        else:
+            return False
+
+    def visit_Module(self, node):
+        self.generic_visit(node)
+        if self.update:
+            import_alias = ast.alias(name='numpy', asname=mangle('numpy'))
+            importIt = ast.Import(names=[import_alias])
+            node.body.insert(0, importIt)
+        return node
+
+    def visit_Assign(self, node):
+        if len(node.targets) != 1:
+            return node
+        target, = node.targets
+        if not self.is_fully_sliced(target):
+            return node
+        if self.is_fully_sliced(node.value):
+            value = node.value.value
+        else:
+            value = node.value
+
+        self.update = True
+        return ast.Expr(
+                ast.Call(
+                    ast.Attribute(ast.Name(mangle('numpy'), ast.Load(), None, None),
+                                  'copyto',
+                                  ast.Load()),
+                    [target.value, value],
+                    [])
+                )
+
+
diff --git a/pythran/pythonic/include/numpy/copyto.hpp b/pythran/pythonic/include/numpy/copyto.hpp
index e0b03ff824..2fefe9470c 100644
--- a/pythran/pythonic/include/numpy/copyto.hpp
+++ b/pythran/pythonic/include/numpy/copyto.hpp
@@ -8,7 +8,23 @@ PYTHONIC_NS_BEGIN
 namespace numpy
 {
   template <class T, class pS, class E>
-  types::ndarray<T, pS> copyto(types::ndarray<T, pS> &out, E const &expr);
+  types::none_type copyto(types::ndarray<T, pS> &out, E const &expr);
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::ndarray<T, pS> &&out, E const &expr);
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::numpy_texpr<types::ndarray<T, pS>> &out, E const &expr);
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::numpy_texpr<types::ndarray<T, pS>> &&out, E const &expr);
+
+  // pythran extensions
+  template <class E, class F>
+  types::none_type copyto(E &out, F const &expr) {
+    out[types::fast_contiguous_slice(0, types::none_type{})] = expr;
+    return {};
+  }
 
   DEFINE_FUNCTOR(pythonic::numpy, copyto);
 }
diff --git a/pythran/pythonic/include/types/list.hpp b/pythran/pythonic/include/types/list.hpp
index 1a8080f5b1..a4204c7f0f 100644
--- a/pythran/pythonic/include/types/list.hpp
+++ b/pythran/pythonic/include/types/list.hpp
@@ -68,7 +68,7 @@ namespace types
         typename std::remove_cv<typename std::remove_reference<T>::type>::type
             _type;
     typedef container<_type> container_type;
-    utils::shared_ref<container_type> data;
+    utils::shared_ref<container_type> _data;
 
     template <class U>
     friend class list;
@@ -100,6 +100,7 @@ namespace types
     static const bool is_vectorizable =
         types::is_vectorizable_dtype<dtype>::value &&
         !std::is_same<S, slice>::value;
+    static const bool is_flat = std::is_same<slice, S>::value;
     static const bool is_strided = std::is_same<slice, S>::value;
 
     using shape_t = types::array<long, value>;
@@ -192,7 +193,7 @@ namespace types
         typename std::remove_cv<typename std::remove_reference<T>::type>::type
             _type;
     typedef container<_type> container_type;
-    utils::shared_ref<container_type> data;
+    utils::shared_ref<container_type> _data;
 
     template <class U, class S>
     friend class sliced_list;
@@ -220,6 +221,7 @@ namespace types
     typedef typename utils::nested_container_value_type<list>::type dtype;
     static const size_t value = utils::nested_container_depth<list>::value;
     static const bool is_vectorizable = types::is_vectorizable<dtype>::value;
+    static const bool is_flat = true;
     static const bool is_strided = false;
 
     // constructors
@@ -325,6 +327,9 @@ namespace types
       return fast(index);
     }
 
+    dtype* data() { return _data->data();}
+    const dtype* data() const { return _data->data();}
+
     // modifiers
     template <class Tp>
     void push_back(Tp &&x);
diff --git a/pythran/pythonic/include/types/ndarray.hpp b/pythran/pythonic/include/types/ndarray.hpp
index 337c7892ff..f6779f7831 100644
--- a/pythran/pythonic/include/types/ndarray.hpp
+++ b/pythran/pythonic/include/types/ndarray.hpp
@@ -216,6 +216,7 @@ namespace types
   template <class T, class pS>
   struct ndarray {
     static const bool is_vectorizable = types::is_vectorizable<T>::value;
+    static const bool is_flat = true;
     static const bool is_strided = false;
 
     /* types */
@@ -646,6 +647,8 @@ namespace types
     flat_iterator fend();
 
     /* member functions */
+    T* data() { return buffer;}
+    T const* data() const { return buffer;}
     long flat_size() const;
     bool may_overlap(ndarray const &) const;
 
diff --git a/pythran/pythonic/include/types/numpy_broadcast.hpp b/pythran/pythonic/include/types/numpy_broadcast.hpp
index 0138f197e2..620c5d8e7f 100644
--- a/pythran/pythonic/include/types/numpy_broadcast.hpp
+++ b/pythran/pythonic/include/types/numpy_broadcast.hpp
@@ -78,6 +78,7 @@ namespace types
     using broadcast_or_broadcasted_t =
         typename broadcast_or_broadcasted<Tp, is_dtype<Tp>::value>::type;
     static const bool is_vectorizable = true;
+    static const bool is_flat = false;
     static const bool is_strided = false;
     using dtype = typename std::remove_reference<T>::type::dtype;
     using value_type = typename std::remove_reference<T>::type::value_type;
@@ -279,6 +280,7 @@ namespace types
     // always)
     using dtype = typename broadcast_dtype<T, B>::type;
     static const bool is_vectorizable = types::is_vectorizable<dtype>::value;
+    static const bool is_flat = false;
     static const bool is_strided = false;
     using value_type = dtype;
     using const_iterator = const_broadcast_iterator<dtype>;
diff --git a/pythran/pythonic/include/types/numpy_expr.hpp b/pythran/pythonic/include/types/numpy_expr.hpp
index 4303a2559c..c6ae90eb17 100644
--- a/pythran/pythonic/include/types/numpy_expr.hpp
+++ b/pythran/pythonic/include/types/numpy_expr.hpp
@@ -604,6 +604,7 @@ namespace types
                              Args>::type>::type::dtype>::value...>::value &&
         types::is_vector_op<
             Op, typename std::remove_reference<Args>::type::dtype...>::value;
+    static const bool is_flat = false;
     static const bool is_strided =
         utils::any_of<std::remove_reference<Args>::type::is_strided...>::value;
 
diff --git a/pythran/pythonic/include/types/numpy_gexpr.hpp b/pythran/pythonic/include/types/numpy_gexpr.hpp
index 14e8af9c98..6dc7c6327a 100644
--- a/pythran/pythonic/include/types/numpy_gexpr.hpp
+++ b/pythran/pythonic/include/types/numpy_gexpr.hpp
@@ -583,6 +583,9 @@ namespace types
          std::is_same<contiguous_normalized_slice,
                       typename std::tuple_element<
                           sizeof...(S) - 1, std::tuple<S...>>::type>::value);
+    static const bool is_flat =
+        std::remove_reference<Arg>::type::is_flat && value == 1 &&
+         utils::all_of<std::is_same<contiguous_normalized_slice, S>::value...>::value;
     static const bool is_strided =
         std::remove_reference<Arg>::type::is_strided ||
         (((sizeof...(S) - count_long<S...>::value) == value) &&
@@ -874,6 +877,9 @@ namespace types
 
     explicit operator bool() const;
 
+
+    dtype* data() { return buffer;}
+    const dtype* data() const { return buffer;}
     long flat_size() const;
     long size() const;
     ndarray<dtype, shape_t> copy() const
diff --git a/pythran/pythonic/include/types/numpy_iexpr.hpp b/pythran/pythonic/include/types/numpy_iexpr.hpp
index eae50de39c..e2004abf99 100644
--- a/pythran/pythonic/include/types/numpy_iexpr.hpp
+++ b/pythran/pythonic/include/types/numpy_iexpr.hpp
@@ -36,6 +36,8 @@ namespace types
     static constexpr size_t value = std::remove_reference<Arg>::type::value - 1;
     static const bool is_vectorizable =
         std::remove_reference<Arg>::type::is_vectorizable;
+    static const bool is_flat =
+        std::remove_reference<Arg>::type::is_flat;
     using dtype = typename std::remove_reference<Arg>::type::dtype;
     using value_type =
         typename std::remove_reference<decltype(numpy_iexpr_helper<value>::get(
@@ -334,6 +336,9 @@ namespace types
       return (*this)[std::get<0>(index)];
     }
 
+    dtype* data() { return buffer;}
+    const dtype* data() const { return buffer;}
+
   private:
     /* compute the buffer offset, returning the offset between the
      * first element of the iexpr and the start of the buffer.
diff --git a/pythran/pythonic/include/types/numpy_texpr.hpp b/pythran/pythonic/include/types/numpy_texpr.hpp
index a9cbe9e8ba..d076a0a0e6 100644
--- a/pythran/pythonic/include/types/numpy_texpr.hpp
+++ b/pythran/pythonic/include/types/numpy_texpr.hpp
@@ -28,6 +28,7 @@ namespace types
   struct numpy_texpr_2 {
     static_assert(E::value == 2, "texpr only implemented for matrices");
     static const bool is_vectorizable = false;
+    static const bool is_flat = false;
     static const bool is_strided = true;
     using Arg = E;
 
diff --git a/pythran/pythonic/include/types/numpy_vexpr.hpp b/pythran/pythonic/include/types/numpy_vexpr.hpp
index ef7e9fa528..02c5f1e5cf 100644
--- a/pythran/pythonic/include/types/numpy_vexpr.hpp
+++ b/pythran/pythonic/include/types/numpy_vexpr.hpp
@@ -12,6 +12,7 @@ namespace types
 
     static constexpr size_t value = T::value;
     static const bool is_vectorizable = false;
+    static const bool is_flat = false;
     using dtype = typename dtype_of<T>::type;
     using value_type = T;
     static constexpr bool is_strided = T::is_strided;
diff --git a/pythran/pythonic/include/types/tuple.hpp b/pythran/pythonic/include/types/tuple.hpp
index 2a50c1100f..72a7c42e62 100644
--- a/pythran/pythonic/include/types/tuple.hpp
+++ b/pythran/pythonic/include/types/tuple.hpp
@@ -327,6 +327,7 @@ namespace types
     static const size_t value =
         utils::nested_container_depth<array_base>::value;
     static const bool is_vectorizable = true;
+    static const bool is_flat = true;
     static const bool is_strided = false;
 
     // flat_size implementation
diff --git a/pythran/pythonic/numpy/copyto.hpp b/pythran/pythonic/numpy/copyto.hpp
index 11647ad7a6..1a455451a0 100644
--- a/pythran/pythonic/numpy/copyto.hpp
+++ b/pythran/pythonic/numpy/copyto.hpp
@@ -2,6 +2,7 @@
 #define PYTHONIC_NUMPY_COPYTO_HPP
 
 #include "pythonic/include/numpy/copyto.hpp"
+#include "pythonic//numpy/asarray.hpp"
 
 #include "pythonic/utils/functor.hpp"
 #include "pythonic/types/ndarray.hpp"
@@ -10,10 +11,59 @@ PYTHONIC_NS_BEGIN
 namespace numpy
 {
   template <class T, class pS, class E>
-  types::ndarray<T, pS> copyto(types::ndarray<T, pS> &out, E const &expr)
+  types::none_type copyto(types::ndarray<T, pS> &out, E const &expr)
   {
-    out[types::contiguous_slice(0, types::none_type{})] = expr;
-    return out;
+    using out_type = types::ndarray<T, pS>;
+    if (may_overlap(out, expr)) {
+      auto aexpr = asarray(expr);
+      utils::broadcast_copy < out_type &, decltype(aexpr), out_type::value,
+             (int)out_type::value - (int)utils::dim_of<E>::value,
+             out_type::is_vectorizable &&
+                 std::is_same<typename out_type::dtype, typename types::dtype_of<E>::type>::value &&
+                 types::is_vectorizable<E>::value > (out, aexpr);
+    }
+    else {
+      utils::broadcast_copy < out_type &, E, out_type::value,
+             (int)out_type::value - (int)utils::dim_of<E>::value,
+             out_type::is_vectorizable &&
+                 std::is_same<typename out_type::dtype, typename types::dtype_of<E>::type>::value &&
+                 types::is_vectorizable<E>::value > (out, expr);
+    }
+    return {};
+  }
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::ndarray<T, pS> &&out, E const &expr)
+  {
+    return copyto(out, expr);
+  }
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::numpy_texpr<types::ndarray<T, pS>> &out, E const &expr)
+  {
+    using out_type = types::numpy_texpr<types::ndarray<T, pS>>;
+    if (may_overlap(out, expr)) {
+      auto aexpr = asarray(expr);
+      utils::broadcast_copy < out_type &, decltype(aexpr), out_type::value,
+             (int)out_type::value - (int)utils::dim_of<E>::value,
+             out_type::is_vectorizable &&
+                 std::is_same<typename out_type::dtype, typename types::dtype_of<E>::type>::value &&
+                 types::is_vectorizable<E>::value > (out, aexpr);
+    }
+    else {
+      utils::broadcast_copy < out_type &, E, out_type::value,
+             (int)out_type::value - (int)utils::dim_of<E>::value,
+             out_type::is_vectorizable &&
+                 std::is_same<typename out_type::dtype, typename types::dtype_of<E>::type>::value &&
+                 types::is_vectorizable<E>::value > (out, expr);
+    }
+    return {};
+  }
+
+  template <class T, class pS, class E>
+  types::none_type copyto(types::numpy_texpr<types::ndarray<T, pS>> &&out, E const &expr)
+  {
+    return copyto(out, expr);
   }
 }
 PYTHONIC_NS_END
diff --git a/pythran/pythonic/types/list.hpp b/pythran/pythonic/types/list.hpp
index 678737f896..98beb401fb 100644
--- a/pythran/pythonic/types/list.hpp
+++ b/pythran/pythonic/types/list.hpp
@@ -23,24 +23,24 @@ namespace types
 
   // Constructors
   template <class T, class S>
-  sliced_list<T, S>::sliced_list() : data(utils::no_memory())
+  sliced_list<T, S>::sliced_list() : _data(utils::no_memory())
   {
   }
   template <class T, class S>
   sliced_list<T, S>::sliced_list(sliced_list<T, S> const &s)
-      : data(s.data), slicing(s.slicing)
+      : _data(s._data), slicing(s.slicing)
   {
   }
   template <class T, class S>
   sliced_list<T, S>::sliced_list(list<T> const &other, S const &s)
-      : data(other.data), slicing(s.normalize(other.size()))
+      : _data(other._data), slicing(s.normalize(other.size()))
   {
   }
   template <class T, class S>
   template <class Sn>
   sliced_list<T, S>::sliced_list(utils::shared_ref<container_type> const &other,
                                  Sn const &s)
-      : data(other), slicing(s)
+      : _data(other), slicing(s)
   {
   }
 
@@ -85,8 +85,8 @@ namespace types
   {
     assert(0 <= i && i < size());
     auto const index = slicing.get(i);
-    assert(0 <= index && index < (long)data->size());
-    return (*data)[index];
+    assert(0 <= index && index < (long)_data->size());
+    return (*_data)[index];
   }
   template <class T, class S>
   typename sliced_list<T, S>::const_reference
@@ -94,16 +94,16 @@ namespace types
   {
     assert(i < size());
     auto const index = slicing.get(i);
-    assert(0 <= index && index < (long)data->size());
-    return (*data)[index];
+    assert(0 <= index && index < (long)_data->size());
+    return (*_data)[index];
   }
   template <class T, class S>
   typename sliced_list<T, S>::reference sliced_list<T, S>::operator[](long i)
   {
     assert(i < size());
     auto const index = slicing.get(i);
-    assert(0 <= index && index < (long)data->size());
-    return (*data)[index];
+    assert(0 <= index && index < (long)_data->size());
+    return (*_data)[index];
   }
 
   template <class T, class S>
@@ -113,7 +113,7 @@ namespace types
       sliced_list<T, decltype(std::declval<S>() * std::declval<Sp>())>>::type
   sliced_list<T, S>::operator[](Sp s) const
   {
-    return {data, slicing * s.normalize(this->size())};
+    return {_data, slicing * s.normalize(this->size())};
   }
 
   // io
@@ -152,10 +152,10 @@ namespace types
   {
     if (slicing.step == 1) {
       // inserting before erasing in case of self-copy
-      auto insert_pt = data->begin() + slicing.lower;
-      data->insert(insert_pt, s.begin(), s.end());
-      auto erase_pt = data->begin() + s.size();
-      data->erase(erase_pt + slicing.lower, erase_pt + slicing.upper);
+      auto insert_pt = _data->begin() + slicing.lower;
+      _data->insert(insert_pt, s.begin(), s.end());
+      auto erase_pt = _data->begin() + s.size();
+      _data->erase(erase_pt + slicing.lower, erase_pt + slicing.upper);
     } else
       assert(!"not implemented yet");
     return *this;
@@ -165,10 +165,10 @@ namespace types
   {
     if (slicing.step == 1) {
       // inserting before erasing in case of self-copy
-      auto insert_pt = data->begin() + slicing.lower;
-      data->insert(insert_pt, seq.begin(), seq.end());
-      auto erase_pt = data->begin() + seq.size();
-      data->erase(erase_pt + slicing.lower, erase_pt + slicing.upper);
+      auto insert_pt = _data->begin() + slicing.lower;
+      _data->insert(insert_pt, seq.begin(), seq.end());
+      auto erase_pt = _data->begin() + seq.size();
+      _data->erase(erase_pt + slicing.lower, erase_pt + slicing.upper);
     } else
       assert(!"not implemented yet");
     return *this;
@@ -208,7 +208,7 @@ namespace types
   typename sliced_list<T, S>::simd_iterator
   sliced_list<T, S>::vbegin(vectorizer) const
   {
-    return {data->data() + slicing.lower};
+    return {_data->data() + slicing.lower};
   }
 
   template <class T, class S>
@@ -218,7 +218,7 @@ namespace types
   {
     using vector_type = typename xsimd::batch<dtype>;
     static const std::size_t vector_size = vector_type::size;
-    return {data->data() + slicing.lower +
+    return {_data->data() + slicing.lower +
             long(size() / vector_size * vector_size)};
   }
 
@@ -229,7 +229,7 @@ namespace types
   template <class V>
   bool sliced_list<T, S>::contains(V const &v) const
   {
-    return std::find(data->begin(), data->end(), v) != data->end();
+    return std::find(_data->begin(), _data->end(), v) != _data->end();
   }
   template <class T, class S>
   intptr_t sliced_list<T, S>::id() const
@@ -248,55 +248,55 @@ namespace types
 
   // constructors
   template <class T>
-  list<T>::list() : data(utils::no_memory())
+  list<T>::list() : _data(utils::no_memory())
   {
   }
   template <class T>
   template <class InputIterator>
-  list<T>::list(InputIterator start, InputIterator stop) : data()
+  list<T>::list(InputIterator start, InputIterator stop) : _data()
   {
     if (std::is_same<
             typename std::iterator_traits<InputIterator>::iterator_category,
             std::random_access_iterator_tag>::value)
-      data->reserve(std::distance(start, stop));
+      _data->reserve(std::distance(start, stop));
     else
-      data->reserve(DEFAULT_LIST_CAPACITY);
-    std::copy(start, stop, std::back_inserter(*data));
+      _data->reserve(DEFAULT_LIST_CAPACITY);
+    std::copy(start, stop, std::back_inserter(*_data));
   }
   template <class T>
-  list<T>::list(empty_list const &) : data(0)
+  list<T>::list(empty_list const &) : _data(0)
   {
   }
   template <class T>
-  list<T>::list(size_type sz) : data(sz)
+  list<T>::list(size_type sz) : _data(sz)
   {
   }
   template <class T>
-  list<T>::list(T const &value, single_value, size_type sz) : data(sz, value)
+  list<T>::list(T const &value, single_value, size_type sz) : _data(sz, value)
   {
   }
   template <class T>
-  list<T>::list(std::initializer_list<T> l) : data(std::move(l))
+  list<T>::list(std::initializer_list<T> l) : _data(std::move(l))
   {
   }
   template <class T>
-  list<T>::list(list<T> &&other) : data(std::move(other.data))
+  list<T>::list(list<T> &&other) : _data(std::move(other._data))
   {
   }
   template <class T>
-  list<T>::list(list<T> const &other) : data(other.data)
+  list<T>::list(list<T> const &other) : _data(other._data)
   {
   }
   template <class T>
   template <class F>
-  list<T>::list(list<F> const &other) : data(other.size())
+  list<T>::list(list<F> const &other) : _data(other.size())
   {
     std::copy(other.begin(), other.end(), begin());
   }
   template <class T>
   template <class Tp, class S>
   list<T>::list(sliced_list<Tp, S> const &other)
-      : data(other.begin(), other.end())
+      : _data(other.begin(), other.end())
   {
   }
 
@@ -304,45 +304,45 @@ namespace types
   template <class T>
   list<T> &list<T>::operator=(list<T> &&other)
   {
-    data = std::move(other.data);
+    _data = std::move(other._data);
     return *this;
   }
   template <class T>
   template <class F>
   list<T> &list<T>::operator=(list<F> const &other)
   {
-    data = utils::shared_ref<container_type>{other.size()};
+    _data = utils::shared_ref<container_type>{other.size()};
     std::copy(other.begin(), other.end(), begin());
     return *this;
   }
   template <class T>
   list<T> &list<T>::operator=(list<T> const &other)
   {
-    data = other.data;
+    _data = other._data;
     return *this;
   }
   template <class T>
   list<T> &list<T>::operator=(empty_list const &)
   {
-    data = utils::shared_ref<container_type>();
+    _data = utils::shared_ref<container_type>();
     return *this;
   }
   template <class T>
   template <class Tp, size_t N, class V>
   list<T> &list<T>::operator=(array_base<Tp, N, V> const &other)
   {
-    data = utils::shared_ref<container_type>(other.begin(), other.end());
+    _data = utils::shared_ref<container_type>(other.begin(), other.end());
     return *this;
   }
   template <class T>
   template <class Tp, class S>
   list<T> &list<T>::operator=(sliced_list<Tp, S> const &other)
   {
-    if (other.data == data) {
-      auto it = std::copy(other.begin(), other.end(), data->begin());
-      data->resize(it - data->begin());
+    if (other._data == _data) {
+      auto it = std::copy(other.begin(), other.end(), _data->begin());
+      _data->resize(it - _data->begin());
     } else
-      data = utils::shared_ref<container_type>(other.begin(), other.end());
+      _data = utils::shared_ref<container_type>(other.begin(), other.end());
     return *this;
   }
 
@@ -350,8 +350,8 @@ namespace types
   template <class S>
   list<T> &list<T>::operator+=(sliced_list<T, S> const &other)
   {
-    data->resize(size() + other.size());
-    std::copy(other.begin(), other.end(), data->begin());
+    _data->resize(size() + other.size());
+    std::copy(other.begin(), other.end(), _data->begin());
     return *this;
   }
 
@@ -419,42 +419,42 @@ namespace types
   template <class T>
   typename list<T>::iterator list<T>::begin()
   {
-    return data->begin();
+    return _data->begin();
   }
   template <class T>
   typename list<T>::const_iterator list<T>::begin() const
   {
-    return data->begin();
+    return _data->begin();
   }
   template <class T>
   typename list<T>::iterator list<T>::end()
   {
-    return data->end();
+    return _data->end();
   }
   template <class T>
   typename list<T>::const_iterator list<T>::end() const
   {
-    return data->end();
+    return _data->end();
   }
   template <class T>
   typename list<T>::reverse_iterator list<T>::rbegin()
   {
-    return data->rbegin();
+    return _data->rbegin();
   }
   template <class T>
   typename list<T>::const_reverse_iterator list<T>::rbegin() const
   {
-    return data->rbegin();
+    return _data->rbegin();
   }
   template <class T>
   typename list<T>::reverse_iterator list<T>::rend()
   {
-    return data->rend();
+    return _data->rend();
   }
   template <class T>
   typename list<T>::const_reverse_iterator list<T>::rend() const
   {
-    return data->rend();
+    return _data->rend();
   }
 
   // comparison
@@ -487,7 +487,7 @@ namespace types
   template <class vectorizer>
   typename list<T>::simd_iterator list<T>::vbegin(vectorizer) const
   {
-    return {data->data()};
+    return {_data->data()};
   }
 
   template <class T>
@@ -496,14 +496,14 @@ namespace types
   {
     using vector_type = typename xsimd::batch<dtype>;
     static const std::size_t vector_size = vector_type::size;
-    return {data->data() + long(size() / vector_size * vector_size)};
+    return {_data->data() + long(size() / vector_size * vector_size)};
   }
 
 #endif
   template <class T>
   typename list<T>::reference list<T>::fast(long n)
   {
-    return (*data)[n];
+    return (*_data)[n];
   }
   template <class T>
   typename list<T>::reference list<T>::operator[](long n)
@@ -517,7 +517,7 @@ namespace types
   typename list<T>::const_reference list<T>::fast(long n) const
   {
     assert(n < size());
-    return (*data)[n];
+    return (*_data)[n];
   }
   template <class T>
   typename list<T>::const_reference list<T>::operator[](long n) const
@@ -542,31 +542,31 @@ namespace types
   void list<T>::push_back(Tp &&x)
   {
     // FIXME: clang-3.4 doesn't support emplace_back for vector of bool
-    data->push_back(std::forward<Tp>(x));
+    _data->push_back(std::forward<Tp>(x));
   }
   template <class T>
   template <class Tp>
   void list<T>::insert(long i, Tp &&x)
   {
     if (i == size())
-      data->emplace_back(std::forward<Tp>(x));
+      _data->emplace_back(std::forward<Tp>(x));
     else
-      data->insert(data->begin() + i, std::forward<Tp>(x));
+      _data->insert(_data->begin() + i, std::forward<Tp>(x));
   }
   template <class T>
   void list<T>::reserve(size_t n)
   {
-    data->reserve(n);
+    _data->reserve(n);
   }
   template <class T>
   void list<T>::resize(size_t n)
   {
-    data->resize(n);
+    _data->resize(n);
   }
   template <class T>
   typename list<T>::iterator list<T>::erase(size_t n)
   {
-    return data->erase(data->begin() + n);
+    return _data->erase(_data->begin() + n);
   }
   template <class T>
   T list<T>::pop(long x)
@@ -598,7 +598,7 @@ namespace types
   template <class T>
   list<T>::operator bool() const
   {
-    return !data->empty();
+    return !_data->empty();
   }
 
   template <class T>
@@ -672,7 +672,7 @@ namespace types
   template <class T>
   long list<T>::size() const
   {
-    return data->size();
+    return _data->size();
   }
   template <class T>
   template <class E>
@@ -697,12 +697,12 @@ namespace types
   template <class V>
   bool list<T>::contains(V const &v) const
   {
-    return std::find(data->begin(), data->end(), v) != data->end();
+    return std::find(_data->begin(), _data->end(), v) != _data->end();
   }
   template <class T>
   intptr_t list<T>::id() const
   {
-    return reinterpret_cast<intptr_t>(&(*data));
+    return reinterpret_cast<intptr_t>(&(*_data));
   }
 
   template <class T>
diff --git a/pythran/pythonic/types/ndarray.hpp b/pythran/pythonic/types/ndarray.hpp
index 55032719b5..fb661da723 100644
--- a/pythran/pythonic/types/ndarray.hpp
+++ b/pythran/pythonic/types/ndarray.hpp
@@ -1055,7 +1055,7 @@ namespace types
   template <class pS>
   list<T> &list<T>::operator=(ndarray<T, pshape<pS>> const &other)
   {
-    data = utils::shared_ref<T>(other.begin(), other.end());
+    _data = utils::shared_ref<T>(other.begin(), other.end());
     return *this;
   }
 } // namespace types
diff --git a/pythran/pythonic/utils/broadcast_copy.hpp b/pythran/pythonic/utils/broadcast_copy.hpp
index 60a0abe4e7..d44ef84bd0 100644
--- a/pythran/pythonic/utils/broadcast_copy.hpp
+++ b/pythran/pythonic/utils/broadcast_copy.hpp
@@ -244,32 +244,62 @@ namespace utils
 
   template <class E, class F, size_t N, int D, bool vector_form>
   E &broadcast_copy_helper(E &self, F const &other,
-                           std::integral_constant<bool, true>)
+                           std::integral_constant<bool, true>,
+                           std::integral_constant<bool, false>)
   {
     static_assert(D >= 0, "downcasting already happened");
-    if (self.size())
+    if (self.size()) {
 #ifdef USE_XSIMD
-      broadcast_copy_dispatcher<E, F, N, (size_t)D, vector_form>{}(self, other);
+      constexpr bool vectorize = vector_form;
 #else
-      broadcast_copy_dispatcher<E, F, N, (size_t)D, false>{}(self, other);
+      constexpr bool vectorize = false;;
 #endif
+      broadcast_copy_dispatcher<E, F, N, (size_t)D, vectorize>{}(self, other);
+    }
     return self;
   }
 
   template <class E, class F, size_t N, int D, bool vector_form>
   E &broadcast_copy_helper(E &self, F const &other,
-                           std::integral_constant<bool, false>)
+                           std::integral_constant<bool, true>,
+                           std::integral_constant<bool, true>)
+  {
+    if(D==0) {
+      std::copy(other.data(), other.data() + other.flat_size(), self.data());
+      return self;
+    }
+    else {
+      return broadcast_copy_helper<E, F, N, D, vector_form>(
+          self, other, std::integral_constant<bool, true>(),
+          std::integral_constant<bool, false>{});
+    }
+  }
+
+  template <class E, class F, size_t N, int D, bool vector_form, bool plain>
+  E &broadcast_copy_helper(E &self, F const &other,
+                           std::integral_constant<bool, false>,
+                           std::integral_constant<bool, plain> is_plain)
   {
     auto reshaped = other.reshape(sutils::getshape(self));
     return broadcast_copy_helper<E, decltype(reshaped), N, 0, vector_form>(
-        self, reshaped, std::true_type());
+        self, reshaped, std::true_type(), is_plain);
   }
 
+  template<class T, bool = types::is_dtype<T>::value>
+  struct is_flat {
+    static const bool value = T::is_flat;
+  };
+  template<class T>
+  struct is_flat<T, true> {
+    static const bool value = false;
+  };
+
   template <class E, class F, size_t N, int D, bool vector_form>
   E &broadcast_copy(E &self, F const &other)
   {
     return broadcast_copy_helper<E, F, N, D, vector_form>(
-        self, other, std::integral_constant<bool, (D >= 0)>());
+        self, other, std::integral_constant<bool, (D >= 0)>(),
+        std::integral_constant<bool, std::decay<E>::type::is_flat && is_flat<typename std::decay<F>::type>::value>{});
   }
 
   /* update
diff --git a/pythran/pythran.cfg b/pythran/pythran.cfg
index dc4ae0e999..6858afce68 100644
--- a/pythran/pythran.cfg
+++ b/pythran/pythran.cfg
@@ -17,6 +17,7 @@ optimizations = pythran.optimizations.InlineBuiltins
                 pythran.transformations.FalsePolymorphism
                 pythran.optimizations.PatternTransform
                 pythran.optimizations.Square
+                pythran.optimizations.CopyTo
                 pythran.optimizations.RangeLoopUnfolding
                 pythran.optimizations.RangeBasedSimplify
                 pythran.optimizations.ListToTuple
diff --git a/pythran/tests/test_ndarray.py b/pythran/tests/test_ndarray.py
index a819b57dc7..d43580291e 100644
--- a/pythran/tests/test_ndarray.py
+++ b/pythran/tests/test_ndarray.py
@@ -722,6 +722,11 @@ def test_gexpr_copy1(self):
                       numpy.arange(16),
                       gexpr_copy1=[NDArray[int, :]])
 
+    def test_gexpr_copy2(self):
+        self.run_test("def gexpr_copy2(a,b): a[:] = b[:]; return a",
+                      numpy.arange(5), [1] * 5,
+                      gexpr_copy2=[NDArray[int, :], List[int]])
+
     def test_ndarray_iter0(self):
         self.run_test("def ndarray_iter0(a): return list(map(str, a))",
                       numpy.arange(16),
diff --git a/pythran/tests/test_numpy_func3.py b/pythran/tests/test_numpy_func3.py
index 61e81c573b..59dde50093 100644
--- a/pythran/tests/test_numpy_func3.py
+++ b/pythran/tests/test_numpy_func3.py
@@ -527,6 +527,11 @@ def test_copyto_1(self):
                       numpy.array([[1,2], [7, 8]]), numpy.array([3,4]),
                       np_copyto1=[NDArray[int, :, :], NDArray[int, :]])
 
+    def test_copyto_2(self):
+        self.run_test("def np_copyto2(x, y): from numpy import copyto; return copyto(x.T, y), x",
+                      numpy.array([[1,2], [7, 8]]), numpy.array([3,4]),
+                      np_copyto2=[NDArray[int, :, :], NDArray[int, :]])
+
     def test_numpy_pow0(self):
         self.run_test('def numpy_pow0(a): return a ** 2',
                       numpy.arange(100).reshape((10, 10)),