Fixed strides adaptor leading to a broadcasting issue

Author: JohanMabille

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))
  {

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,49 +115,55 @@ namespace xt
  inline self_type& operator++()
  {
  ++p_current;
+ ++p_shape;
  return *this;
  }
 
  inline self_type& operator--()
  {
  --p_current;
+ --p_shape;
  return *this;
  }
 
  inline self_type operator++(int)
  {
  self_type tmp(*this);
  ++p_current;
+ ++p_shape;
  return tmp;
  }
 
  inline self_type operator--(int)
  {
  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_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); }));

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):