Repository: xtensor-stack/xtensor-python (language: python)

Issue:
Broadcasting problem on basic operations?
Hi there,

It seems that there may be some broadcasting problems with `xt::pyarray` on basic operations with `operator+`, `operator-`, `operator*`, and `operator/`. 

For example, with `operator-`:
```cpp
#define FORCE_IMPORT_ARRAY
#include "xtensor-python/pyarray.hpp"

xt::pyarray<double> arr1
        {{ 4., 5., 6., 7. }};
xt::pyarray<double> arr2
        {{ 4., 3., 2., 1. },
         { 4., 3., 2., 1. },
         { 4., 3., 2., 1. }};

std::cout << arr1 - arr2 << std::endl;
```
returns

```text
{{ 0.000000e+000,  2.000000e+000,  4.000000e+000,  6.000000e+000},
 { 2.889806e+204, -3.000000e+000, -2.000000e+000,  1.000000e+000},
 { 0.000000e+000,  3.000000e+000,  1.722192e+097, -1.000000e+000}}
```

For simplicity, I have added some tests to your existing test suite in a fork of mine, so that you can try to reproduce the problem yourself. See:
https://github.com/ThibHlln/xtensor-python/blob/broadcasting-problem/test_python/main.cpp#L114-L132
https://github.com/ThibHlln/xtensor-python/blob/broadcasting-problem/test_python/test_pyarray.py#L66-L108 

All 4 tests fail for me.

Could you tell me if these tests pass for you please? And if so, any idea what might be wrong here?

Thank you for your help,
Thibault

*Note: I am building/running this on Windows 10, version 1803 with `pybind11==2.9.2`, `numpy==1.22.3`, `xtensor==0.24.1`, `xtl==0.7.4` and using `xtensor-python` from the master branch.*

**EDIT**
If this helps, it works fine if I use `xt::pytensor<double, 2>` in place of `xt::pyarray<double>`.


Discussion / hints:
That sounds like a bug. Thanks for reporting.
Hi,

Does the following work?

```cpp
py::array<double> res = arr1 - arr2;
std::cout << res << std::endl;
```
Hi,

Sorry, my example was misleading: this is when the Python interface is used that the problem occurs as I tried to demonstrate in your existing testing suite with:
https://github.com/ThibHlln/xtensor-python/blob/broadcasting-problem/test_python/main.cpp#L114-L132
https://github.com/ThibHlln/xtensor-python/blob/broadcasting-problem/test_python/test_pyarray.py#L66-L108

If the pyarrays are defined directly in C++, the broadcasting works. This is when working with pyarrays defined via the Python bindings that the problem occurs (observable both in the returned numpy array and if I do a `std::cout << arr1 - arr2 << std::endl;` in C++).

I can reproduce the problem, I am investigating.

Fix patch:
diff --git a/include/xtensor-python/pyarray.hpp b/include/xtensor-python/pyarray.hpp
index d80b19b..3e6a6ba 100644
--- a/include/xtensor-python/pyarray.hpp
+++ b/include/xtensor-python/pyarray.hpp
@@ -519,7 +519,8 @@ namespace xt
         m_shape = inner_shape_type(reinterpret_cast<size_type*>(PyArray_SHAPE(this->python_array())),
                                    static_cast<size_type>(PyArray_NDIM(this->python_array())));
         m_strides = inner_strides_type(reinterpret_cast<difference_type*>(PyArray_STRIDES(this->python_array())),
-                                       static_cast<size_type>(PyArray_NDIM(this->python_array())));
+                                       static_cast<size_type>(PyArray_NDIM(this->python_array())),
+                                       reinterpret_cast<size_type*>(PyArray_SHAPE(this->python_array())));
 
         if (L != layout_type::dynamic && !do_strides_match(m_shape, m_strides, L, 1))
         {
diff --git a/include/xtensor-python/pystrides_adaptor.hpp b/include/xtensor-python/pystrides_adaptor.hpp
index 161255b..fde929a 100644
--- a/include/xtensor-python/pystrides_adaptor.hpp
+++ b/include/xtensor-python/pystrides_adaptor.hpp
@@ -41,8 +41,10 @@ namespace xt
         using const_reverse_iterator = std::reverse_iterator<const_iterator>;
         using reverse_iterator = const_reverse_iterator;
 
+        using shape_type = size_t*;
+
         pystrides_adaptor() = default;
-        pystrides_adaptor(const_pointer data, size_type size);
+        pystrides_adaptor(const_pointer data, size_type size, shape_type shape);
 
         bool empty() const noexcept;
         size_type size() const noexcept;
@@ -66,6 +68,7 @@ namespace xt
 
         const_pointer p_data;
         size_type m_size;
+        shape_type p_shape;
     };
 
     /**********************************
@@ -84,21 +87,23 @@ namespace xt
         using reference = typename pystrides_adaptor<N>::const_reference;
         using difference_type = typename pystrides_adaptor<N>::difference_type;
         using iterator_category = std::random_access_iterator_tag;
+        using shape_pointer = typename pystrides_adaptor<N>::shape_type;
 
-        inline pystrides_iterator(pointer current)
+        inline pystrides_iterator(pointer current, shape_pointer shape)
             : p_current(current)
+            , p_shape(shape)
         {
         }
 
         inline reference operator*() const
         {
-            return *p_current / N;
+            return *p_shape == size_t(1) ? 0 : *p_current / N;
         }
 
         inline pointer operator->() const
         {
             // Returning the address of a temporary
-            value_type res = *p_current / N;
+            value_type res = this->operator*();
             return &res;
         }
 
@@ -110,12 +115,14 @@ namespace xt
         inline self_type& operator++()
         {
             ++p_current;
+            ++p_shape;
             return *this;
         }
 
         inline self_type& operator--()
         {
             --p_current;
+            --p_shape;
             return *this;
         }
 
@@ -123,6 +130,7 @@ namespace xt
         {
             self_type tmp(*this);
             ++p_current;
+            ++p_shape;
             return tmp;
         }
 
@@ -130,29 +138,32 @@ namespace xt
         {
             self_type tmp(*this);
             --p_current;
+            --p_shape;
             return tmp;
         }
 
         inline self_type& operator+=(difference_type n)
         {
             p_current += n;
+            p_shape += n;
             return *this;
         }
 
         inline self_type& operator-=(difference_type n)
         {
             p_current -= n;
+            p_shape -= n;
             return *this;
         }
 
         inline self_type operator+(difference_type n) const
         {
-            return self_type(p_current + n);
+            return self_type(p_current + n, p_shape + n);
         }
 
         inline self_type operator-(difference_type n) const
         {
-            return self_type(p_current - n);
+            return self_type(p_current - n, p_shape - n);
         }
 
         inline difference_type operator-(const self_type& rhs) const
@@ -166,6 +177,7 @@ namespace xt
     private:
 
         pointer p_current;
+        shape_pointer p_shape;
     };
 
     template <std::size_t N>
@@ -215,8 +227,8 @@ namespace xt
      ************************************/
 
     template <std::size_t N>
-    inline pystrides_adaptor<N>::pystrides_adaptor(const_pointer data, size_type size)
-        : p_data(data), m_size(size)
+    inline pystrides_adaptor<N>::pystrides_adaptor(const_pointer data, size_type size, shape_type shape)
+        : p_data(data), m_size(size), p_shape(shape)
     {
     }
 
@@ -235,19 +247,19 @@ namespace xt
     template <std::size_t N>
     inline auto pystrides_adaptor<N>::operator[](size_type i) const -> const_reference
     {
-        return p_data[i] / N;
+        return p_shape[i] == size_t(1) ? 0 : p_data[i] / N;
     }
 
     template <std::size_t N>
     inline auto pystrides_adaptor<N>::front() const -> const_reference
     {
-        return p_data[0] / N;
+        return this->operator[](0);
     }
 
     template <std::size_t N>
     inline auto pystrides_adaptor<N>::back() const -> const_reference
     {
-        return p_data[m_size - 1] / N;
+        return this->operator[](m_size - 1);
     }
 
     template <std::size_t N>
@@ -265,13 +277,13 @@ namespace xt
     template <std::size_t N>
     inline auto pystrides_adaptor<N>::cbegin() const -> const_iterator
     {
-        return const_iterator(p_data);
+        return const_iterator(p_data, p_shape);
     }
 
     template <std::size_t N>
     inline auto pystrides_adaptor<N>::cend() const -> const_iterator
     {
-        return const_iterator(p_data + m_size);
+        return const_iterator(p_data + m_size, p_shape + m_size);
     }
 
     template <std::size_t N>


Test patch:
diff --git a/test_python/main.cpp b/test_python/main.cpp
index 5a5edfb..7a0c524 100644
--- a/test_python/main.cpp
+++ b/test_python/main.cpp
@@ -108,6 +108,28 @@ auto no_complex_overload_reg(const double& a)
 {
     return a;
 }
+//
+// Operator examples
+//
+xt::pyarray<double> array_addition(const xt::pyarray<double>& m, const xt::pyarray<double>& n)
+{
+    return m + n;
+}
+
+xt::pyarray<double> array_subtraction(xt::pyarray<double>& m, xt::pyarray<double>& n)
+{
+    return m - n;
+}
+
+xt::pyarray<double> array_multiplication(xt::pyarray<double>& m, xt::pyarray<double>& n)
+{
+    return m * n;
+}
+
+xt::pyarray<double> array_division(xt::pyarray<double>& m, xt::pyarray<double>& n)
+{
+    return m / n;
+}
 
 // Vectorize Examples
 
@@ -310,6 +332,11 @@ PYBIND11_MODULE(xtensor_python_test, m)
     m.def("readme_example1", readme_example1);
     m.def("readme_example2", xt::pyvectorize(readme_example2));
 
+    m.def("array_addition", array_addition);
+    m.def("array_subtraction", array_subtraction);
+    m.def("array_multiplication", array_multiplication);
+    m.def("array_division", array_division);
+    
     m.def("vectorize_example1", xt::pyvectorize(add));
 
     m.def("rect_to_polar", xt::pyvectorize([](complex_t x) { return std::abs(x); }));
diff --git a/test_python/test_pyarray.py b/test_python/test_pyarray.py
index c413002..9955ee7 100644
--- a/test_python/test_pyarray.py
+++ b/test_python/test_pyarray.py
@@ -23,7 +23,7 @@
 import numpy as np
 
 class XtensorTest(TestCase):
-
+    """
     def test_rm(self):
         xt.test_rm(np.array([10], dtype=int))
 
@@ -62,6 +62,50 @@ def test_example3(self):
         with self.assertRaises(TypeError):
             x = np.arange(3*2).reshape(3, 2)
             xt.example3_xfixed2(x)
+    """
+    def test_broadcast_addition(self):
+        x = np.array([[2., 3., 4., 5.]])
+        y = np.array([[1., 2., 3., 4.],
+                      [1., 2., 3., 4.],
+                      [1., 2., 3., 4.]])
+        res = np.array([[3., 5., 7., 9.],
+                        [3., 5., 7., 9.],
+                        [3., 5., 7., 9.]])
+        z = xt.array_addition(x, y)
+        np.testing.assert_allclose(z, res, 1e-12)
+    """
+    def test_broadcast_subtraction(self):
+        x = np.array([[4., 5., 6., 7.]])
+        y = np.array([[4., 3., 2., 1.],
+                      [4., 3., 2., 1.],
+                      [4., 3., 2., 1.]])
+        res = np.array([[0., 2., 4., 6.],
+                        [0., 2., 4., 6.],
+                        [0., 2., 4., 6.]])
+        z = xt.array_subtraction(x, y)
+        np.testing.assert_allclose(z, res, 1e-12)
+
+    def test_broadcast_multiplication(self):
+        x = np.array([[1., 2., 3., 4.]])
+        y = np.array([[3., 2., 3., 2.],
+                      [3., 2., 3., 2.],
+                      [3., 2., 3., 2.]])
+        res = np.array([[3., 4., 9., 8.],
+                        [3., 4., 9., 8.],
+                        [3., 4., 9., 8.]])
+        z = xt.array_multiplication(x, y)
+        np.testing.assert_allclose(z, res, 1e-12)
+
+    def test_broadcast_division(self):
+        x = np.array([[8., 6., 4., 2.]])
+        y = np.array([[2., 2., 2., 2.],
+                      [2., 2., 2., 2.],
+                      [2., 2., 2., 2.]])
+        res = np.array([[4., 3., 2., 1.],
+                        [4., 3., 2., 1.],
+                        [4., 3., 2., 1.]])
+        z = xt.array_division(x, y)
+        np.testing.assert_allclose(z, res, 1e-12)
 
     def test_vectorize(self):
         x1 = np.array([[0, 1], [2, 3]])
@@ -263,7 +307,7 @@ def test_native_casters(self):
         self.assertEqual(adapter.shape, (2, 2))
         adapter[1, 1] = -3
         self.assertEqual(arr[0, 5], -3)
-
+     """
 
 class AttributeTest(TestCase):