Making a copy when casting a numpy array with negative strides to Eigen.

`EigenConformable::stride_compatible` returns false if the strides are
negative. In this case, do not use `EigenConformable::stride`, as it
is {0,0}. We cannot write negative strides in this element, as Eigen
will throw an assertion if we do.

The `type_caster` specialization for regular, dense Eigen matrices now
does a second `array_t::ensure` to copy data in case of negative strides.
I'm not sure that this is the best way to implement this.

I have added "TODO" tags linking these changes to Eigen bug #747, which,
when fixed, will allow Eigen to accept negative strides.

Author: Cris Luengo

include/pybind11/eigen.h

@@ -68,16 +68,22 @@ template <typename T> using is_eigen_other = all_of<
 template <bool EigenRowMajor> struct EigenConformable {
     bool conformable = false;
     EigenIndex rows = 0, cols = 0;
-    EigenDStride stride{0, 0};
+    EigenDStride stride{0, 0};      // Only valid if negativestridees is false!
+    bool negativestrides = false;   // If true, do not use stride!
 
     EigenConformable(bool fits = false) : conformable{fits} {}
     // Matrix type:
     EigenConformable(EigenIndex r, EigenIndex c,
             EigenIndex rstride, EigenIndex cstride) :
-        conformable{true}, rows{r}, cols{c},
-        stride(EigenRowMajor ? rstride : cstride /* outer stride */,
-               EigenRowMajor ? cstride : rstride /* inner stride */)
-        {}
+        conformable{true}, rows{r}, cols{c} {
+        // TODO: when Eigen bug #747 is fixed, remove the tests for non-negativity. http://eigen.tuxfamily.org/bz/show_bug.cgi?id=747
+        if (rstride < 0 || cstride < 0) {
+            negativestrides = true;
+        } else {
+            stride = {EigenRowMajor ? rstride : cstride /* outer stride */,
+                      EigenRowMajor ? cstride : rstride /* inner stride */ };
+        }
+    }
     // Vector type:
     EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride)
         : EigenConformable(r, c, r == 1 ? c*stride : stride, c == 1 ? r : r*stride) {}
@@ -86,6 +92,7 @@ template <bool EigenRowMajor> struct EigenConformable {
         // To have compatible strides, we need (on both dimensions) one of fully dynamic strides,
         // matching strides, or a dimension size of 1 (in which case the stride value is irrelevant)
         return
+            !negativestrides &&
             (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner() ||
                 (EigenRowMajor ? cols : rows) == 1) &&
             (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer() ||
@@ -138,8 +145,8 @@ template <typename Type_> struct EigenProps {
             EigenIndex
                 np_rows = a.shape(0),
                 np_cols = a.shape(1),
-                np_rstride = a.strides(0) / sizeof(Scalar),
-                np_cstride = a.strides(1) / sizeof(Scalar);
+                np_rstride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar)),
+                np_cstride = a.strides(1) / static_cast<ssize_t>(sizeof(Scalar));
             if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols))
                 return false;
 
@@ -149,7 +156,7 @@ template <typename Type_> struct EigenProps {
         // Otherwise we're storing an n-vector.  Only one of the strides will be used, but whichever
         // is used, we want the (single) numpy stride value.
         const EigenIndex n = a.shape(0),
-              stride = a.strides(0) / sizeof(Scalar);
+              stride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar));
 
         if (vector) { // Eigen type is a compile-time vector
             if (fixed && size != n)
@@ -255,7 +262,23 @@ struct type_caster<Type, enable_if_t<is_eigen_dense_plain<Type>::value>> {
         if (!fits)
             return false; // Non-comformable vector/matrix types
 
-        value = Eigen::Map<const Type, 0, EigenDStride>(buf.data(), fits.rows, fits.cols, fits.stride);
+        if (fits.negativestrides) {
+
+            // Eigen does not support negative strides, so we need to make a copy here with normal strides.
+            // TODO: when Eigen bug #747 is fixed, remove this if case, always execute the else part.
+            // http://eigen.tuxfamily.org/bz/show_bug.cgi?id=747
+            auto buf2 = array_t<Scalar,array::forcecast || array::f_style>::ensure(src);
+            if (!buf2)
+                return false;
+            // not checking sizes, we already did that
+            fits = props::conformable(buf2);
+            value = Eigen::Map<const Type, 0, EigenDStride>(buf2.data(), fits.rows, fits.cols, fits.stride);
+
+        } else {
+
+           value = Eigen::Map<const Type, 0, EigenDStride>(buf.data(), fits.rows, fits.cols, fits.stride);
+
+        }
 
         return true;
     }

include/pybind11/numpy.h

@@ -248,10 +248,10 @@ template <typename T> using is_pod_struct = all_of<
     satisfies_none_of<T, std::is_reference, std::is_array, is_std_array, std::is_arithmetic, is_complex, std::is_enum>
 >;
 
-template <size_t Dim = 0, typename Strides> size_t byte_offset_unsafe(const Strides &) { return 0; }
+template <size_t Dim = 0, typename Strides> ssize_t byte_offset_unsafe(const Strides &) { return 0; }
 template <size_t Dim = 0, typename Strides, typename... Ix>
-size_t byte_offset_unsafe(const Strides &strides, size_t i, Ix... index) {
-    return i * strides[Dim] + byte_offset_unsafe<Dim + 1>(strides, index...);
+ssize_t byte_offset_unsafe(const Strides &strides, size_t i, Ix... index) {
+    return static_cast<ssize_t>(i) * strides[Dim] + byte_offset_unsafe<Dim + 1>(strides, index...);
 }
 
 /** Proxy class providing unsafe, unchecked const access to array data.  This is constructed through
@@ -615,18 +615,18 @@ class array : public buffer {
 
     /// Byte offset from beginning of the array to a given index (full or partial).
     /// May throw if the index would lead to out of bounds access.
-    template<typename... Ix> size_t offset_at(Ix... index) const {
+    template<typename... Ix> ssize_t offset_at(Ix... index) const {
         if (sizeof...(index) > ndim())
             fail_dim_check(sizeof...(index), "too many indices for an array");
         return byte_offset(size_t(index)...);
     }
 
-    size_t offset_at() const { return 0; }
+    ssize_t offset_at() const { return 0; }
 
     /// Item count from beginning of the array to a given index (full or partial).
     /// May throw if the index would lead to out of bounds access.
-    template<typename... Ix> size_t index_at(Ix... index) const {
-        return offset_at(index...) / itemsize();
+    template<typename... Ix> ssize_t index_at(Ix... index) const {
+        return offset_at(index...) / static_cast<ssize_t>(itemsize());
     }
 
     /** Returns a proxy object that provides access to the array's data without bounds or
@@ -692,7 +692,7 @@ class array : public buffer {
                           " (ndim = " + std::to_string(ndim()) + ")");
     }
 
-    template<typename... Ix> size_t byte_offset(Ix... index) const {
+    template<typename... Ix> ssize_t byte_offset(Ix... index) const {
         check_dimensions(index...);
         return detail::byte_offset_unsafe(strides(), size_t(index)...);
     }
@@ -773,8 +773,8 @@ template <typename T, int ExtraFlags = array::forcecast> class array_t : public
         return sizeof(T);
     }
 
-    template<typename... Ix> size_t index_at(Ix... index) const {
-        return offset_at(index...) / itemsize();
+    template<typename... Ix> ssize_t index_at(Ix... index) const {
+        return offset_at(index...) / static_cast<ssize_t>(itemsize());
     }
 
     template<typename... Ix> const T* data(Ix... index) const {
@@ -789,14 +789,14 @@ template <typename T, int ExtraFlags = array::forcecast> class array_t : public
     template<typename... Ix> const T& at(Ix... index) const {
         if (sizeof...(index) != ndim())
             fail_dim_check(sizeof...(index), "index dimension mismatch");
-        return *(static_cast<const T*>(array::data()) + byte_offset(size_t(index)...) / itemsize());
+        return *(static_cast<const T*>(array::data()) + byte_offset(size_t(index)...) / static_cast<ssize_t>(itemsize()));
     }
 
     // Mutable reference to element at a given index
     template<typename... Ix> T& mutable_at(Ix... index) {
         if (sizeof...(index) != ndim())
             fail_dim_check(sizeof...(index), "index dimension mismatch");
-        return *(static_cast<T*>(array::mutable_data()) + byte_offset(size_t(index)...) / itemsize());
+        return *(static_cast<T*>(array::mutable_data()) + byte_offset(size_t(index)...) / static_cast<ssize_t>(itemsize()));
     }
 
     /** Returns a proxy object that provides access to the array's data without bounds or

include/pybind11/stl_bind.h

@@ -350,14 +350,14 @@ vector_buffer(Class_& cl) {
 
     cl.def("__init__", [](Vector& vec, buffer buf) {
         auto info = buf.request();
-        if (info.ndim != 1 || info.strides[0] <= 0 || info.strides[0] % sizeof(T))
+        if (info.ndim != 1 || info.strides[0] <= 0 || info.strides[0] % static_cast<ssize_t>(sizeof(T)))
             throw type_error("Only valid 1D buffers can be copied to a vector");
         if (!detail::compare_buffer_info<T>::compare(info) || sizeof(T) != info.itemsize)
             throw type_error("Format mismatch (Python: " + info.format + " C++: " + format_descriptor<T>::format() + ")");
         new (&vec) Vector();
         vec.reserve(info.shape[0]);
         T *p = static_cast<T*>(info.ptr);
-        auto step = info.strides[0] / sizeof(T);
+        auto step = info.strides[0] / static_cast<ssize_t>(sizeof(T));
         T *end = p + info.shape[0] * step;
         for (; p < end; p += step)
             vec.push_back(*p);

tests/test_eigen.py

@@ -154,6 +154,53 @@ def test_nonunit_stride_from_python():
     np.testing.assert_array_equal(counting_mat, [[0., 2, 2], [6, 16, 10], [6, 14, 8]])
 
 
+def test_negative_stride_from_python(msg):
+    from pybind11_tests import (
+        double_row, double_col, double_complex, double_mat_cm, double_mat_rm,
+        double_threec, double_threer)
+
+    # Eigen doesn't support (as of yet) negative strides. When a function takes an Eigen
+    # matrix by copy or const reference, we can pass a numpy array that has negative strides.
+    # Otherwise, an exception will be thrown as Eigen will not be able to map the numpy array.
+
+    counting_mat = np.arange(9.0, dtype=np.float32).reshape((3, 3))
+    counting_mat = counting_mat[::-1, ::-1]
+    second_row = counting_mat[1, :]
+    second_col = counting_mat[:, 1]
+    np.testing.assert_array_equal(double_row(second_row), 2.0 * second_row)
+    np.testing.assert_array_equal(double_col(second_row), 2.0 * second_row)
+    np.testing.assert_array_equal(double_complex(second_row), 2.0 * second_row)
+    np.testing.assert_array_equal(double_row(second_col), 2.0 * second_col)
+    np.testing.assert_array_equal(double_col(second_col), 2.0 * second_col)
+    np.testing.assert_array_equal(double_complex(second_col), 2.0 * second_col)
+
+    counting_3d = np.arange(27.0, dtype=np.float32).reshape((3, 3, 3))
+    counting_3d = counting_3d[::-1, ::-1, ::-1]
+    slices = [counting_3d[0, :, :], counting_3d[:, 0, :], counting_3d[:, :, 0]]
+    for slice_idx, ref_mat in enumerate(slices):
+        np.testing.assert_array_equal(double_mat_cm(ref_mat), 2.0 * ref_mat)
+        np.testing.assert_array_equal(double_mat_rm(ref_mat), 2.0 * ref_mat)
+
+    # Mutator:
+    with pytest.raises(TypeError) as excinfo:
+        double_threer(second_row)
+    assert msg(excinfo.value) == """
+    double_threer(): incompatible function arguments. The following argument types are supported:
+        1. (numpy.ndarray[float32[1, 3], flags.writeable]) -> arg0: None
+
+    Invoked with: array([ 5.,  4.,  3.], dtype=float32)
+"""
+
+    with pytest.raises(TypeError) as excinfo:
+        double_threec(second_col)
+    assert msg(excinfo.value) == """
+    double_threec(): incompatible function arguments. The following argument types are supported:
+        1. (numpy.ndarray[float32[3, 1], flags.writeable]) -> arg0: None
+
+    Invoked with: array([ 7.,  4.,  1.], dtype=float32)
+"""
+
+
 def test_nonunit_stride_to_python():
     from pybind11_tests import diagonal, diagonal_1, diagonal_n, block