diff --git a/docs/advanced/cast/eigen.rst b/docs/advanced/cast/eigen.rst index b83ca9af9..5b0b08ca6 100644 --- a/docs/advanced/cast/eigen.rst +++ b/docs/advanced/cast/eigen.rst @@ -1,48 +1,308 @@ Eigen -===== +##### `Eigen `_ is C++ header-based library for dense and sparse linear algebra. Due to its popularity and widespread adoption, pybind11 -provides transparent conversion support between Eigen and Scientific Python linear -algebra data types. +provides transparent conversion and limited mapping support between Eigen and +Scientific Python linear algebra data types. -Specifically, when including the optional header file :file:`pybind11/eigen.h`, -pybind11 will automatically and transparently convert +To enable the built-in Eigen support you must include the optional header file +:file:`pybind11/eigen.h`. -1. Static and dynamic Eigen dense vectors and matrices to instances of - ``numpy.ndarray`` (and vice versa). +Pass-by-value +============= -2. Returned matrix expressions such as blocks (including columns or rows) and - diagonals will be converted to ``numpy.ndarray`` of the expression - values. +When binding a function with ordinary Eigen dense object arguments (for +example, ``Eigen::MatrixXd``), pybind11 will accept any input value that is +already (or convertible to) a ``numpy.ndarray`` with dimensions compatible with +the Eigen type, copy its values into a temporary Eigen variable of the +appropriate type, then call the function with this temporary variable. -3. Returned matrix-like objects such as Eigen::DiagonalMatrix or - Eigen::SelfAdjointView will be converted to ``numpy.ndarray`` containing the - expressed value. +Sparse matrices are similarly copied to or from +``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` objects. -4. Eigen sparse vectors and matrices to instances of - ``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` (and vice versa). +Pass-by-reference +================= -This makes it possible to bind most kinds of functions that rely on these types. -One major caveat are functions that take Eigen matrices *by reference* and modify -them somehow, in which case the information won't be propagated to the caller. +One major limitation of the above is that every data conversion implicitly +involves a copy, which can be both expensive (for large matrices) and disallows +binding functions that change their (Matrix) arguments. Pybind11 allows you to +work around this by using Eigen's ``Eigen::Ref`` class much as you +would when writing a function taking a generic type in Eigen itself (subject to +some limitations discussed below). + +When calling a bound function accepting a ``Eigen::Ref`` +type, pybind11 will attempt to avoid copying by using an ``Eigen::Map`` object +that maps into the source ``numpy.ndarray`` data: this requires both that the +data types are the same (e.g. ``dtype='float64'`` and ``MatrixType::Scalar`` is +``double``); and that the storage is layout compatible. The latter limitation +is discussed in detail in the section below, and requires careful +consideration: by default, numpy matrices and eigen matrices are *not* storage +compatible. + +If the numpy matrix cannot be used as is (either because its types differ, e.g. +passing an array of integers to an Eigen paramater requiring doubles, or +because the storage is incompatible), pybind11 makes a temporary copy and +passes the copy instead. + +When a bound function parameter is instead ``Eigen::Ref`` (note the +lack of ``const``), pybind11 will only allow the function to be called if it +can be mapped *and* if the numpy array is writeable (that is +``a.flags.writeable`` is true). Any access (including modification) made to +the passed variable will be transparently carried out directly on the +``numpy.ndarray``. + +This means you can can write code such as the following and have it work as +expected: .. code-block:: cpp - /* The Python bindings of these functions won't replicate - the intended effect of modifying the function arguments */ - void scale_by_2(Eigen::Vector3f &v) { - v *= 2; - } - void scale_by_2(Eigen::Ref &v) { + void scale_by_2(Eigen::Ref m) { v *= 2; } -To see why this is, refer to the section on :ref:`opaque` (although that -section specifically covers STL data types, the underlying issue is the same). -The :ref:`numpy` sections discuss an efficient alternative for exposing the -underlying native Eigen types as opaque objects in a way that still integrates -with NumPy and SciPy. +Note, however, that you will likely run into limitations due to numpy and +Eigen's difference default storage order for data; see the below section on +:ref:`storage_orders` for details on how to bind code that won't run into such +limitations. + +.. note:: + + Passing by reference is not supported for sparse types. + +Returning values to Python +========================== + +When returning an ordinary dense Eigen matrix type to numpy (e.g. +``Eigen::MatrixXd`` or ``Eigen::RowVectorXf``) pybind11 keeps the matrix and +returns a numpy array that directly references the Eigen matrix: no copy of the +data is performed. The numpy array will have ``array.flags.owndata`` set to +``False`` to indicate that it does not own the data, and the lifetime of the +stored Eigen matrix will be tied to the returned ``array``. + +If you bind a function with a non-reference, ``const`` return type (e.g. +``const Eigen::MatrixXd``), the same thing happens except that pybind11 also +sets the numpy array's ``writeable`` flag to false. + +If you return an lvalue reference or pointer, the usual pybind11 rules apply, +as dictated by the binding function's return value policy (see the +documentation on :ref:`return_value_policies` for full details). That means, +without an explicit return value policy, lvalue references will be copied and +pointers will be managed by pybind11. In order to avoid copying, you should +explictly specify an appropriate return value policy, as in the following +example: + +.. code-block:: cpp + + class MyClass { + Eigen::MatrixXd big_mat = Eigen::MatrixXd::Zero(10000, 10000); + public: + Eigen::MatrixXd &getMatrix() { return big_mat; } + const Eigen::MatrixXd &viewMatrix() { return big_mat; } + }; + + // Later, in binding code: + py::class_(m, "MyClass") + .def(py::init<>()) + .def("copy_matrix", &MyClass::getMatrix) // Makes a copy! + .def("get_matrix", &MyClass::getMatrix, py::return_value_policy::reference_internal) + .def("view_matrix", &MyClass::viewMatrix, py::return_value_policy::reference_internal) + ; + +.. code-block:: python + + a = MyClass() + m = a.get_matrix() # flags.writeable = True, flags.owndata = False + v = a.view_matrix() # flags.writeable = False, flags.owndata = False + c = a.copy_matrix() # flags.writeable = True, flags.owndata = True + # m[5,6] and v[5,6] refer to the same element, c[5,6] does not. + +Note in this example that ``py::return_value_policy::reference_internal`` is +used to tie the life of the MyClass object to the life of the returned arrays. + +You may also return an ``Eigen::Ref``, ``Eigen::Map`` or other map-like Eigen +object (for example, the return value of ``matrix.block()`` and related +methods) that map into a dense Eigen type. When doing so, the default +behaviour of pybind11 is to simply reference the returned data: you must take +care to ensure that this data remains valid! You may ask pybind11 to +explicitly *copy* such a return value by using the +``py::return_value_policy::copy`` policy when binding the function. You may +also use ``py::return_value_policy::reference_internal`` or a +``py::keep_alive`` to ensure the data stays valid as long as the returned numpy +array does. + +When returning such a reference of map, pybind11 additionally respects the +readonly-status of the returned value, marking the numpy array as non-writeable +if the reference or map was itself read-only. + +.. note:: + + Sparse types are always copied when returned. + +.. _storage_orders: + +Storage orders +============== + +Passing arguments via ``Eigen::Ref`` has some limitations that you must be +aware of in order to effectively pass matrices by reference. First and +foremost is that the default ``Eigen::Ref`` class requires +contiguous storage along columns (for column-major types, the default in Eigen) +or rows if ``MatrixType`` is specifically an ``Eigen::RowMajor`` storage type. +The former, Eigen's default, is incompatible with ``numpy``'s default row-major +storage, and so you will not be able to pass numpy arrays to Eigen by reference +without making one of two changes. + +(Note that this does not apply to vectors (or column or row matrices): for such +types the "row-major" and "column-major" distinction is meaningless). + +The first approach is to change the use of ``Eigen::Ref`` to the +more general ``Eigen::Ref>`` (or similar type with a fully dynamic stride type in the +third template argument). Since this is a rather cumbersome type, pybind11 +provides a ``py::EigenDRef`` type alias for your convenience (along +with EigenDMap for the equivalent Map, and EigenDStride for just the stride +type). + +This type allows Eigen to map into any arbitrary storage order. This is not +the default in Eigen for performance reasons: contiguous storage allows +vectorization that cannot be done when storage is not known to be contiguous at +compile time. The default ``Eigen::Ref`` stride type allows non-contiguous +storage along the outer dimension (that is, the rows of a column-major matrix +or columns of a row-major matrix), but not along the inner dimension. + +This type, however, has the added benefit of also being able to map numpy array +slices. For example, the following (contrived) example uses Eigen with a numpy +slice to multiply by 2 all coefficients that are both on even rows (0, 2, 4, +...) and in columns 2, 5, or 8: + +.. code-block:: cpp + + m.def("scale", [](py::EigenDRef m, double c) { m *= c; }); + +.. code-block:: python + + # a = np.array(...) + scale_by_2(myarray[0::2, 2:9:3]) + +The second approach to avoid copying is more intrusive: rearranging the +underlying data types to not run into the non-contiguous storage problem in the +first place. In particular, that means using matrices with ``Eigen::RowMajor`` +storage, where appropriate, such as: + +.. code-block:: cpp + + using RowMatrixXd = Eigen::Matrix; + // Use RowMatrixXd instead of MatrixXd + +Now bound functions accepting ``Eigen::Ref`` arguments will be +callable with numpy's (default) arrays without involving a copying. + +You can, alternatively, change the storage order that numpy arrays use by +adding the ``order='F'`` option when creating an array: + +.. code-block:: python + + myarray = np.array(source, order='F') + +Such an object will be passable to a bound function accepting an +``Eigen::Ref`` (or similar column-major Eigen type). + +One major caveat with this approach, however, is that it is not entirely as +easy as simply flipping all Eigen or numpy usage from one to the other: some +operations may alter the storage order of a numpy array. For example, ``a2 = +array.transpose()`` results in ``a2`` being a view of ``array`` that references +the same data, but in the opposite storage order! + +While this approach allows fully optimized vectorized calculations in Eigen, it +cannot be used with array slices, unlike the first approach. + +When *returning* a matrix to Python (either a regular matrix, a reference via +``Eigen::Ref<>``, or a map/block into a matrix), no special storage +consideration is required: the created numpy array will have the required +stride that allows numpy to properly interpret the array, whatever its storage +order. + +Failing rather than copying +=========================== + +The default behaviour when binding ``Eigen::Ref`` eigen +references is to copy matrix values when passed a numpy array that does not +conform to the element type of ``MatrixType`` or does not have a compatible +stride layout. If you want to explicitly avoid copying in such a case, you +should bind arguments using the ``py::arg().noconvert()`` annotation (as +described in the :ref:`nonconverting_arguments` documentation). + +The following example shows an example of arguments that don't allow data +copying to take place: + +.. code-block:: cpp + + // The method and function to be bound: + class MyClass { + // ... + double some_method(const Eigen::Ref &matrix) { /* ... */ } + }; + float some_function(const Eigen::Ref &big, + const Eigen::Ref &small) { + // ... + } + + // The associated binding code: + using namespace pybind11::literals; // for "arg"_a + py::class_(m, "MyClass") + // ... other class definitions + .def("some_method", &MyClass::some_method, py::arg().nocopy()); + + m.def("some_function", &some_function, + "big"_a.nocopy(), // <- Don't allow copying for this arg + "small"_a // <- This one can be copied if needed + ); + +With the above binding code, attempting to call the the ``some_method(m)`` +method on a ``MyClass`` object, or attempting to call ``some_function(m, m2)`` +will raise a ``RuntimeError`` rather than making a temporary copy of the array. +It will, however, allow the ``m2`` argument to be copied into a temporary if +necessary. + +Note that explicitly specifying ``.noconvert()`` is not required for *mutable* +Eigen references (e.g. ``Eigen::Ref`` without ``const`` on the +``MatrixXd``): mutable references will never be called with a temporary copy. + +Vectors versus column/row matrices +================================== + +Eigen and numpy have fundamentally different notions of a vector. In Eigen, a +vector is simply a matrix with the number of columns or rows set to 1 at +compile time (for a column vector or row vector, respectively). Numpy, in +contast, has comparable 2-dimensional 1xN and Nx1 arrays, but *also* has +1-dimensional arrays of size N. + +When passing a 2-dimensional 1xN or Nx1 array to Eigen, the Eigen type must +have matching dimensions: That is, you cannot pass a 2-dimensional Nx1 numpy +array to an Eigen value expecting a row vector, or a 1xN numpy array as a +column vector argument. + +On the other hand, pybind11 allows you to pass 1-dimensional arrays of length N +as Eigen parameters. If the Eigen type can hold a column vector of length N it +will be passed as such a column vector. If not, but the Eigen type constraints +will accept a row vector, it will be passed as a row vector. (The column +vector takes precendence when both are supported, for example, when passing a +1D numpy array to a MatrixXd argument). Note that the type need not be +expicitly a vector: it is permitted to pass a 1D numpy array of size 5 to an +Eigen ``Matrix``: you would end up with a 1x5 Eigen matrix. +Passing the same to an ``Eigen::MatrixXd`` would result in a 5x1 Eigen matrix. + +When returning an eigen vector to numpy, the conversion is ambiguous: a row +vector of length 4 could be returned as either a 1D array of length 4, or as a +2D array of size 1x4. When encoutering such a situation, pybind11 compromises +by considering the returned Eigen type: if it is a compile-time vector--that +is, the type has either the number of rows or columns set to 1 at compile +time--pybind11 converts to a 1D numpy array when returning the value. For +instances that are a vector only at run-time (e.g. ``MatrixXd``, +``Matrix``), pybind11 returns the vector as a 2D array to +numpy. If this isn't want you want, you can use ``array.reshape(...)`` to get +a view of the same data in the desired dimensions. .. seealso:: diff --git a/docs/advanced/functions.rst b/docs/advanced/functions.rst index 6c52ab1fa..e0b0fe095 100644 --- a/docs/advanced/functions.rst +++ b/docs/advanced/functions.rst @@ -6,6 +6,8 @@ with the basics of binding functions and classes, as explained in :doc:`/basics` and :doc:`/classes`. The following guide is applicable to both free and member functions, i.e. *methods* in Python. +.. _return_value_policies: + Return value policies ===================== @@ -320,6 +322,8 @@ like so: py::class_("MyClass") .def("myFunction", py::arg("arg") = (SomeType *) nullptr); +.. _nonconverting_arguments: + Non-converting arguments ======================== diff --git a/include/pybind11/eigen.h b/include/pybind11/eigen.h index a0d1951e8..5ffbc7e09 100644 --- a/include/pybind11/eigen.h +++ b/include/pybind11/eigen.h @@ -30,156 +30,487 @@ # pragma warning(disable: 4127) // warning C4127: Conditional expression is constant #endif +// Eigen prior to 3.2.7 doesn't have proper move constructors--but worse, some classes get implicit +// move constructors that break things. We could detect this an explicitly copy, but an extra copy +// of matrices seems highly undesirable. +static_assert(EIGEN_VERSION_AT_LEAST(3,2,7), "Eigen support in pybind11 requires Eigen >= 3.2.7"); + NAMESPACE_BEGIN(pybind11) + +// Provide a convenience alias for easier pass-by-ref usage with fully dynamic strides: +using EigenDStride = Eigen::Stride; +template using EigenDRef = Eigen::Ref; +template using EigenDMap = Eigen::Map; + NAMESPACE_BEGIN(detail) -template using is_eigen_dense = is_template_base_of; -template using is_eigen_sparse = is_template_base_of; -template using is_eigen_ref = is_template_base_of; +#if EIGEN_VERSION_AT_LEAST(3,3,0) +using EigenIndex = Eigen::Index; +#else +using EigenIndex = EIGEN_DEFAULT_DENSE_INDEX_TYPE; +#endif +// Matches Eigen::Map, Eigen::Ref, blocks, etc: +template using is_eigen_dense_map = all_of, std::is_base_of, T>>; +template using is_eigen_mutable_map = std::is_base_of, T>; +template using is_eigen_dense_plain = all_of>, is_template_base_of>; +template using is_eigen_sparse = is_template_base_of; // Test for objects inheriting from EigenBase that aren't captured by the above. This // basically covers anything that can be assigned to a dense matrix but that don't have a typical // matrix data layout that can be copied from their .data(). For example, DiagonalMatrix and // SelfAdjointView fall into this category. -template using is_eigen_base = all_of< +template using is_eigen_other = all_of< is_template_base_of, - negation>, - negation> + negation, is_eigen_dense_plain, is_eigen_sparse>> >; +// Captures numpy/eigen conformability status (returned by EigenProps::conformable()): +template struct EigenConformable { + bool conformable = false; + EigenIndex rows = 0, cols = 0; + EigenDStride stride{0, 0}; + + 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 */) + {} + // Vector type: + EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride) : EigenConformable(r, c, r == 1 ? c*stride : stride, c == 1 ? r : r*stride) {} + template bool stride_compatible() const { + return + (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner()) && + (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer()); + } + operator bool() const { return conformable; } +}; + +template struct eigen_extract_stride { using type = Type; }; +template +struct eigen_extract_stride> { using type = StrideType; }; +template +struct eigen_extract_stride> { using type = StrideType; }; + +// Helper struct for extracting information from an Eigen type +template struct EigenProps { + using Type = Type_; + using Scalar = typename Type::Scalar; + using StrideType = typename eigen_extract_stride::type; + static constexpr EigenIndex + rows = Type::RowsAtCompileTime, + cols = Type::ColsAtCompileTime, + size = Type::SizeAtCompileTime; + static constexpr bool + row_major = Type::IsRowMajor, + vector = Type::IsVectorAtCompileTime, // At least one dimension has fixed size 1 + fixed_rows = rows != Eigen::Dynamic, + fixed_cols = cols != Eigen::Dynamic, + fixed = size != Eigen::Dynamic, // Fully-fixed size + dynamic = !fixed_rows && !fixed_cols; // Fully-dynamic size + + template using if_zero = std::integral_constant; + static constexpr EigenIndex inner_stride = if_zero::value, + outer_stride = if_zero::value; + static constexpr bool dynamic_stride = inner_stride == Eigen::Dynamic && outer_stride == Eigen::Dynamic; + static constexpr bool requires_row_major = !dynamic_stride && !vector && (row_major ? inner_stride : outer_stride) == 1; + static constexpr bool requires_col_major = !dynamic_stride && !vector && (row_major ? outer_stride : inner_stride) == 1; + + // Takes an input array and determines whether we can make it fit into the Eigen type. If + // the array is a vector, we attempt to fit it into either an Eigen 1xN or Nx1 vector + // (preferring the latter if it will fit in either, i.e. for a fully dynamic matrix type). + static EigenConformable conformable(const array &a) { + const auto dims = a.ndim(); + if (dims < 1 || dims > 2) + return false; + + if (dims == 2) { // Matrix type: require exact match (or dynamic) + + 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); + if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols)) + return false; + + return {np_rows, np_cols, np_rstride, np_cstride}; + } + + // 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); + + if (vector) { // Eigen type is a compile-time vector + if (fixed && size != n) + return false; // Vector size mismatch + return {rows == 1 ? 1 : n, cols == 1 ? 1 : n, stride}; + } + else if (fixed) { + // The type has a fixed size, but is not a vector: abort + return false; + } + else if (fixed_cols) { + // Since this isn't a vector, cols must be != 1. We allow this only if it exactly + // equals the number of elements (rows is Dynamic, and so 1 row is allowed). + if (cols != n) return false; + return {1, n, stride}; + } + else { + // Otherwise it's either fully dynamic, or column dynamic; both become a column vector + if (fixed_rows && rows != n) return false; + return {n, 1, stride}; + } + } + + static PYBIND11_DESCR descriptor() { + constexpr bool show_writeable = is_eigen_dense_map::value && is_eigen_mutable_map::value; + constexpr bool show_order = is_eigen_dense_map::value; + constexpr bool show_c_contiguous = show_order && requires_row_major; + constexpr bool show_f_contiguous = !show_c_contiguous && show_order && requires_col_major; + + return _("numpy.ndarray[") + npy_format_descriptor::name() + + _("[") + _(_<(size_t) rows>(), _("m")) + + _(", ") + _(_<(size_t) cols>(), _("n")) + + _("]") + + // For a reference type (e.g. Ref) we have other constraints that might need to be + // satisfied: writeable=True (for a mutable reference), and, depending on the map's stride + // options, possibly f_contiguous or c_contiguous. We include them in the descriptor output + // to provide some hint as to why a TypeError is occurring (otherwise it can be confusing to + // see that a function accepts a 'numpy.ndarray[float64[3,2]]' and an error message that you + // *gave* a numpy.ndarray of the right type and dimensions. + _(", flags.writeable", "") + + _(", flags.c_contiguous", "") + + _(", flags.f_contiguous", "") + + _("]"); + } +}; + +// Casts an Eigen type to numpy array. If given a base, the numpy array references the src data, +// otherwise it'll make a copy. writeable lets you turn off the writeable flag for the array. +template handle eigen_array_cast(typename props::Type const &src, handle base = handle(), bool writeable = true) { + constexpr size_t elem_size = sizeof(typename props::Scalar); + std::vector shape, strides; + if (props::vector) { + shape.push_back(src.size()); + strides.push_back(elem_size * src.innerStride()); + } + else { + shape.push_back(src.rows()); + shape.push_back(src.cols()); + strides.push_back(elem_size * src.rowStride()); + strides.push_back(elem_size * src.colStride()); + } + array a(std::move(shape), std::move(strides), src.data(), base); + if (!writeable) + array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; + + return a.release(); +} + +// Takes an lvalue ref to some Eigen type and a (python) base object, creating a numpy array that +// reference the Eigen object's data with `base` as the python-registered base class (if omitted, +// the base will be set to None, and lifetime management is up to the caller). The numpy array is +// non-writeable if the given type is const. +template +handle eigen_ref_array(Type &src, handle parent = none()) { + // none here is to get past array's should-we-copy detection, which currently always + // copies when there is no base. Setting the base to None should be harmless. + return eigen_array_cast(src, parent, !std::is_const::value); +} + +// Takes a pointer to some dense, plain Eigen type, builds a capsule around it, then returns a numpy +// array that references the encapsulated data with a python-side reference to the capsule to tie +// its destruction to that of any dependent python objects. Const-ness is determined by whether or +// not the Type of the pointer given is const. +template ::value>> +handle eigen_encapsulate(Type *src) { + capsule base(src, [](PyObject *o) { delete reinterpret_steal(o).operator Type*(); }); + return eigen_ref_array(*src, base); +} + +// Type caster for regular, dense matrix types (e.g. MatrixXd), but not maps/refs/etc. of dense +// types. template -struct type_caster::value && !is_eigen_ref::value>> { - typedef typename Type::Scalar Scalar; - static constexpr bool rowMajor = Type::IsRowMajor; - static constexpr bool isVector = Type::IsVectorAtCompileTime; - static constexpr auto nRows = Type::RowsAtCompileTime, nCols = Type::ColsAtCompileTime; +struct type_caster::value>> { + using Scalar = typename Type::Scalar; + using props = EigenProps; bool load(handle src, bool) { auto buf = array_t::ensure(src); if (!buf) return false; - using namespace Eigen; - - using Strides = Stride; - using AMap = Map; - - if (buf.ndim() == 1) { // A one-dimensional array - Index n_elts = buf.shape(0); - size_t str = buf.strides(0) / sizeof(Scalar); - Strides stride(str, str); // Whether we map to inner or outer is irrelevant - if (isVector) { - if (Type::SizeAtCompileTime != Dynamic && Type::SizeAtCompileTime != n_elts) - return false; // Vector size mismatch - value = AMap(buf.mutable_data(), nRows == 1 ? 1 : n_elts, nCols == 1 ? 1 : n_elts, stride); - } - else if (Type::SizeAtCompileTime != Dynamic) { - // The type has a fixed size, but is not a vector: abort - return false; - } - else if (nRows == Dynamic && nCols == Dynamic) { - // Fully dynamic size. numpy doesn't distinguish between a row vector and column - // vector, so we'll (arbitrarily) choose a column vector. - value = AMap(buf.mutable_data(), n_elts, 1, stride); - } - else if (nRows != Dynamic) { - // Since this isn't a vector, nRows must be != 1. We allow this only if it exactly - // equals the number of elements (nCols is Dynamic, and so 1 is allowed). - if (nRows != n_elts) return false; - value = AMap(buf.mutable_data(), n_elts, 1, stride); - } - else { // nCols != Dynamic; same as above, but for fixed columns - if (nCols != n_elts) return false; - value = AMap(buf.mutable_data(), 1, n_elts, stride); - } - } else if (buf.ndim() == 2) { - - if ((Type::RowsAtCompileTime != Dynamic && buf.shape(0) != (size_t) Type::RowsAtCompileTime) || - (Type::ColsAtCompileTime != Dynamic && buf.shape(1) != (size_t) Type::ColsAtCompileTime)) - return false; - - value = AMap( - buf.mutable_data(), buf.shape(0), buf.shape(1), - Strides(buf.strides(rowMajor ? 0 : 1) / sizeof(Scalar), - buf.strides(rowMajor ? 1 : 0) / sizeof(Scalar)) - ); - } else { + auto dims = buf.ndim(); + if (dims < 1 || dims > 2) return false; - } + + auto fits = props::conformable(buf); + if (!fits) + return false; // Non-comformable vector/matrix types + + value = Eigen::Map(buf.data(), fits.rows, fits.cols, fits.stride); + return true; } - static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) { - if (isVector) { - return array( - { (size_t) src.size() }, // shape - { sizeof(Scalar) * static_cast(src.innerStride()) }, // strides - src.data() // data - ).release(); - } else { - return array( - { (size_t) src.rows(), // shape - (size_t) src.cols() }, - { sizeof(Scalar) * static_cast(src.rowStride()), // strides - sizeof(Scalar) * static_cast(src.colStride()) }, - src.data() // data - ).release(); +private: + + // Cast implementation + template + static handle cast_impl(CType *src, return_value_policy policy, handle parent) { + switch (policy) { + case return_value_policy::take_ownership: + case return_value_policy::automatic: + return eigen_encapsulate(src); + case return_value_policy::move: + return eigen_encapsulate(new CType(std::move(*src))); + case return_value_policy::copy: + return eigen_array_cast(*src); + case return_value_policy::reference: + case return_value_policy::automatic_reference: + return eigen_ref_array(*src); + case return_value_policy::reference_internal: + return eigen_ref_array(*src, parent); + default: + throw cast_error("unhandled return_value_policy: should not happen!"); + }; + } + +public: + + // Normal returned non-reference, non-const value: + static handle cast(Type &&src, return_value_policy /* policy */, handle parent) { + return cast_impl(&src, return_value_policy::move, parent); + } + // If you return a non-reference const, we mark the numpy array readonly: + static handle cast(const Type &&src, return_value_policy /* policy */, handle parent) { + return cast_impl(&src, return_value_policy::move, parent); + } + // lvalue reference return; default (automatic) becomes copy + static handle cast(Type &src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference) + policy = return_value_policy::copy; + return cast_impl(&src, policy, parent); + } + // const lvalue reference return; default (automatic) becomes copy + static handle cast(const Type &src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference) + policy = return_value_policy::copy; + return cast(&src, policy, parent); + } + // non-const pointer return + static handle cast(Type *src, return_value_policy policy, handle parent) { + return cast_impl(src, policy, parent); + } + // const pointer return + static handle cast(const Type *src, return_value_policy policy, handle parent) { + return cast_impl(src, policy, parent); + } + + static PYBIND11_DESCR name() { return type_descr(props::descriptor()); } + + operator Type*() { return &value; } + operator Type&() { return value; } + template using cast_op_type = cast_op_type; + +private: + Type value; +}; + +// Eigen Ref/Map classes have slightly different policy requirements, meaning we don't want to force +// `move` when a Ref/Map rvalue is returned; we treat Ref<> sort of like a pointer (we care about +// the underlying data, not the outer shell). +template +struct return_value_policy_override::value>> { + static return_value_policy policy(return_value_policy p) { return p; } +}; + +// Base class for casting reference/map/block/etc. objects back to python. +template struct eigen_map_caster { +private: + using props = EigenProps; + +public: + + // Directly referencing a ref/map's data is a bit dangerous (whatever the map/ref points to has + // to stay around), but we'll allow it under the assumption that you know what you're doing (and + // have an appropriate keep_alive in place). We return a numpy array pointing directly at the + // ref's data (The numpy array ends up read-only if the ref was to a const matrix type.) Note + // that this means you need to ensure you don't destroy the object in some other way (e.g. with + // an appropriate keep_alive, or with a reference to a statically allocated matrix). + static handle cast(const MapType &src, return_value_policy policy, handle parent) { + switch (policy) { + case return_value_policy::copy: + return eigen_array_cast(src); + case return_value_policy::reference_internal: + return eigen_array_cast(src, parent, is_eigen_mutable_map::value); + case return_value_policy::reference: + case return_value_policy::automatic: + case return_value_policy::automatic_reference: + return eigen_array_cast(src, none(), is_eigen_mutable_map::value); + default: + // move, take_ownership don't make any sense for a ref/map: + pybind11_fail("Invalid return_value_policy for Eigen Map/Ref/Block type"); } } - PYBIND11_TYPE_CASTER(Type, _("numpy.ndarray[") + npy_format_descriptor::name() + - _("[") + rows() + _(", ") + cols() + _("]]")); + static PYBIND11_DESCR name() { return props::descriptor(); } -protected: - template = 0> - static PYBIND11_DESCR rows() { return _("m"); } - template = 0> - static PYBIND11_DESCR rows() { return _(); } - template = 0> - static PYBIND11_DESCR cols() { return _("n"); } - template = 0> - static PYBIND11_DESCR cols() { return _(); } + // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return + // types but not bound arguments). We still provide them (with an explicitly delete) so that + // you end up here if you try anyway. + bool load(handle, bool) = delete; + operator MapType() = delete; + template using cast_op_type = MapType; }; -// Eigen::Ref satisfies is_eigen_dense, but isn't constructable, so it needs a special -// type_caster to handle argument copying/forwarding. -template -struct type_caster> { -protected: - using Type = Eigen::Ref; - using Derived = typename std::remove_const::type; - using DerivedCaster = make_caster; - DerivedCaster derived_caster; - std::unique_ptr value; +// We can return any map-like object (but can only load Refs, specialized next): +template struct type_caster::value>> + : eigen_map_caster {}; + +// Loader for Ref<...> arguments. See the documentation for info on how to make this work without +// copying (it requires some extra effort in many cases). +template +struct type_caster< + Eigen::Ref, + enable_if_t>::value> +> : public eigen_map_caster> { +private: + using Type = Eigen::Ref; + using props = EigenProps; + using Scalar = typename props::Scalar; + using MapType = Eigen::Map; + using Array = array_t; + static constexpr bool need_writeable = is_eigen_mutable_map::value; + // Delay construction (these have no default constructor) + std::unique_ptr map; + std::unique_ptr ref; + // Our array. When possible, this is just a numpy array pointing to the source data, but + // sometimes we can't avoid copying (e.g. input is not a numpy array at all, has an incompatible + // layout, or is an array of a type that needs to be converted). Using a numpy temporary + // (rather than an Eigen temporary) saves an extra copy when we need both type conversion and + // storage order conversion. (Note that we refuse to use this temporary copy when loading an + // argument for a Ref with M non-const, i.e. a read-write reference). + Array copy_or_ref; public: - bool load(handle src, bool convert) { if (derived_caster.load(src, convert)) { value.reset(new Type(derived_caster.operator Derived&())); return true; } return false; } - static handle cast(const Type &src, return_value_policy policy, handle parent) { return DerivedCaster::cast(src, policy, parent); } - static handle cast(const Type *src, return_value_policy policy, handle parent) { return DerivedCaster::cast(*src, policy, parent); } + bool load(handle src, bool convert) { + // First check whether what we have is already an array of the right type. If not, we can't + // avoid a copy (because the copy is also going to do type conversion). + bool need_copy = !isinstance(src); - static PYBIND11_DESCR name() { return DerivedCaster::name(); } + EigenConformable fits; + if (!need_copy) { + // We don't need a converting copy, but we also need to check whether the strides are + // compatible with the Ref's stride requirements + Array aref = reinterpret_borrow(src); - operator Type*() { return value.get(); } - operator Type&() { if (!value) pybind11_fail("Eigen::Ref<...> value not loaded"); return *value; } + if (aref && (!need_writeable || aref.writeable())) { + fits = props::conformable(aref); + if (!fits) return false; // Incompatible dimensions + if (!fits.template stride_compatible()) + need_copy = true; + else + copy_or_ref = std::move(aref); + } + else { + need_copy = true; + } + } + + if (need_copy) { + // We need to copy: If we need a mutable reference, or we're not supposed to convert + // (either because we're in the no-convert overload pass, or because we're explicitly + // instructed not to copy (via `py::arg().noconvert()`) we have to fail loading. + if (!convert || need_writeable) return false; + + Array copy = Array::ensure(src); + if (!copy) return false; + fits = props::conformable(copy); + if (!fits || !fits.template stride_compatible()) + return false; + copy_or_ref = copy; + } + + ref.reset(); + map.reset(new MapType(data(copy_or_ref), fits.rows, fits.cols, make_stride(fits.stride.outer(), fits.stride.inner()))); + ref.reset(new Type(*map)); + + return true; + } + + operator Type*() { return ref.get(); } + operator Type&() { return *ref; } template using cast_op_type = pybind11::detail::cast_op_type<_T>; + +private: + template ::value, int> = 0> + Scalar *data(Array &a) { return a.mutable_data(); } + + template ::value, int> = 0> + const Scalar *data(Array &a) { return a.data(); } + + // Attempt to figure out a constructor of `Stride` that will work. + // If both strides are fixed, use a default constructor: + template using stride_ctor_default = bool_constant< + S::InnerStrideAtCompileTime != Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic && + std::is_default_constructible::value>; + // Otherwise, if there is a two-index constructor, assume it is (outer,inner) like + // Eigen::Stride, and use it: + template using stride_ctor_dual = bool_constant< + !stride_ctor_default::value && std::is_constructible::value>; + // Otherwise, if there is a one-index constructor, and just one of the strides is dynamic, use + // it (passing whichever stride is dynamic). + template using stride_ctor_outer = bool_constant< + !any_of, stride_ctor_dual>::value && + S::OuterStrideAtCompileTime == Eigen::Dynamic && S::InnerStrideAtCompileTime != Eigen::Dynamic && + std::is_constructible::value>; + template using stride_ctor_inner = bool_constant< + !any_of, stride_ctor_dual>::value && + S::InnerStrideAtCompileTime == Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic && + std::is_constructible::value>; + + template ::value, int> = 0> + static S make_stride(EigenIndex, EigenIndex) { return S(); } + template ::value, int> = 0> + static S make_stride(EigenIndex outer, EigenIndex inner) { return S(outer, inner); } + template ::value, int> = 0> + static S make_stride(EigenIndex outer, EigenIndex) { return S(outer); } + template ::value, int> = 0> + static S make_stride(EigenIndex, EigenIndex inner) { return S(inner); } + }; -// type_caster for special matrix types (e.g. DiagonalMatrix): load() is not supported, but we can -// cast them into the python domain by first copying to a regular Eigen::Matrix, then casting that. +// type_caster for special matrix types (e.g. DiagonalMatrix), which are EigenBase, but not +// EigenDense (i.e. they don't have a data(), at least not with the usual matrix layout). +// load() is not supported, but we can cast them into the python domain by first copying to a +// regular Eigen::Matrix, then casting that. template -struct type_caster::value && !is_eigen_ref::value>> { +struct type_caster::value>> { protected: - using Matrix = Eigen::Matrix; - using MatrixCaster = make_caster; + using Matrix = Eigen::Matrix; + using props = EigenProps; public: - [[noreturn]] bool load(handle, bool) { pybind11_fail("Unable to load() into specialized EigenBase object"); } - static handle cast(const Type &src, return_value_policy policy, handle parent) { return MatrixCaster::cast(Matrix(src), policy, parent); } - static handle cast(const Type *src, return_value_policy policy, handle parent) { return MatrixCaster::cast(Matrix(*src), policy, parent); } + static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) { + handle h = eigen_encapsulate(new Matrix(src)); + return h; + } + static handle cast(const Type *src, return_value_policy policy, handle parent) { return cast(*src, policy, parent); } - static PYBIND11_DESCR name() { return MatrixCaster::name(); } + static PYBIND11_DESCR name() { return props::descriptor(); } - [[noreturn]] operator Type*() { pybind11_fail("Loading not supported for specialized EigenBase object"); } - [[noreturn]] operator Type&() { pybind11_fail("Loading not supported for specialized EigenBase object"); } - template using cast_op_type = pybind11::detail::cast_op_type<_T>; + // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return + // types but not bound arguments). We still provide them (with an explicitly delete) so that + // you end up here if you try anyway. + bool load(handle, bool) = delete; + operator Type() = delete; + template using cast_op_type = Type; }; template diff --git a/tests/test_eigen.cpp b/tests/test_eigen.cpp index b10d025d3..329054338 100644 --- a/tests/test_eigen.cpp +++ b/tests/test_eigen.cpp @@ -8,29 +8,47 @@ */ #include "pybind11_tests.h" +#include "constructor_stats.h" #include #include -Eigen::VectorXf double_col(const Eigen::VectorXf& x) -{ return 2.0f * x; } +using MatrixXdR = Eigen::Matrix; -Eigen::RowVectorXf double_row(const Eigen::RowVectorXf& x) -{ return 2.0f * x; } -Eigen::MatrixXf double_mat_cm(const Eigen::MatrixXf& x) -{ return 2.0f * x; } -// Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended -Eigen::MatrixXd cholesky1(Eigen::Ref &x) { return x.llt().matrixL(); } -Eigen::MatrixXd cholesky2(const Eigen::Ref &x) { return x.llt().matrixL(); } -Eigen::MatrixXd cholesky3(const Eigen::Ref &x) { return x.llt().matrixL(); } -Eigen::MatrixXd cholesky4(Eigen::Ref &x) { return x.llt().matrixL(); } -Eigen::MatrixXd cholesky5(Eigen::Ref x) { return x.llt().matrixL(); } -Eigen::MatrixXd cholesky6(Eigen::Ref x) { return x.llt().matrixL(); } +// Sets/resets a testing reference matrix to have values of 10*r + c, where r and c are the +// (1-based) row/column number. +template void reset_ref(M &x) { + for (int i = 0; i < x.rows(); i++) for (int j = 0; j < x.cols(); j++) + x(i, j) = 11 + 10*i + j; +} -typedef Eigen::Matrix MatrixXfRowMajor; -MatrixXfRowMajor double_mat_rm(const MatrixXfRowMajor& x) -{ return 2.0f * x; } +// Returns a static, column-major matrix +Eigen::MatrixXd &get_cm() { + static Eigen::MatrixXd *x; + if (!x) { + x = new Eigen::MatrixXd(3, 3); + reset_ref(*x); + } + return *x; +} +// Likewise, but row-major +MatrixXdR &get_rm() { + static MatrixXdR *x; + if (!x) { + x = new MatrixXdR(3, 3); + reset_ref(*x); + } + return *x; +} +// Resets the values of the static matrices returned by get_cm()/get_rm() +void reset_refs() { + reset_ref(get_cm()); + reset_ref(get_rm()); +} + +// Returns element 2,1 from a matrix (used to test copy/nocopy) +double get_elem(Eigen::Ref m) { return m(2, 1); }; test_initializer eigen([](py::module &m) { typedef Eigen::Matrix FixedMatrixR; @@ -46,21 +64,84 @@ test_initializer eigen([](py::module &m) { m.attr("have_eigen") = true; - // Non-symmetric matrix with zero elements - Eigen::MatrixXf mat(5, 6); - mat << 0, 3, 0, 0, 0, 11, 22, 0, 0, 0, 17, 11, 7, 5, 0, 1, 0, 11, 0, - 0, 0, 0, 0, 11, 0, 0, 14, 0, 8, 11; + m.def("double_col", [](const Eigen::VectorXf &x) -> Eigen::VectorXf { return 2.0f * x; }); + m.def("double_row", [](const Eigen::RowVectorXf &x) -> Eigen::RowVectorXf { return 2.0f * x; }); + m.def("double_threec", [](py::EigenDRef x) { x *= 2; }); + m.def("double_threer", [](py::EigenDRef x) { x *= 2; }); + m.def("double_mat_cm", [](Eigen::MatrixXf x) -> Eigen::MatrixXf { return 2.0f * x; }); + m.def("double_mat_rm", [](DenseMatrixR x) -> DenseMatrixR { return 2.0f * x; }); - m.def("double_col", &double_col); - m.def("double_row", &double_row); - m.def("double_mat_cm", &double_mat_cm); - m.def("double_mat_rm", &double_mat_rm); - m.def("cholesky1", &cholesky1); - m.def("cholesky2", &cholesky2); - m.def("cholesky3", &cholesky3); - m.def("cholesky4", &cholesky4); - m.def("cholesky5", &cholesky5); - m.def("cholesky6", &cholesky6); + // Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended + m.def("cholesky1", [](Eigen::Ref x) -> Eigen::MatrixXd { return x.llt().matrixL(); }); + m.def("cholesky2", [](const Eigen::Ref &x) -> Eigen::MatrixXd { return x.llt().matrixL(); }); + m.def("cholesky3", [](const Eigen::Ref &x) -> Eigen::MatrixXd { return x.llt().matrixL(); }); + m.def("cholesky4", [](Eigen::Ref x) -> Eigen::MatrixXd { return x.llt().matrixL(); }); + + // Mutators: these add some value to the given element using Eigen, but Eigen should be mapping into + // the numpy array data and so the result should show up there. There are three versions: one that + // works on a contiguous-row matrix (numpy's default), one for a contiguous-column matrix, and one + // for any matrix. + auto add_rm = [](Eigen::Ref x, int r, int c, double v) { x(r,c) += v; }; + auto add_cm = [](Eigen::Ref x, int r, int c, double v) { x(r,c) += v; }; + + // Mutators (Eigen maps into numpy variables): + m.def("add_rm", add_rm); // Only takes row-contiguous + m.def("add_cm", add_cm); // Only takes column-contiguous + // Overloaded versions that will accept either row or column contiguous: + m.def("add1", add_rm); + m.def("add1", add_cm); + m.def("add2", add_cm); + m.def("add2", add_rm); + // This one accepts a matrix of any stride: + m.def("add_any", [](py::EigenDRef x, int r, int c, double v) { x(r,c) += v; }); + + // Return mutable references (numpy maps into eigen varibles) + m.def("get_cm_ref", []() { return Eigen::Ref(get_cm()); }); + m.def("get_rm_ref", []() { return Eigen::Ref(get_rm()); }); + // The same references, but non-mutable (numpy maps into eigen variables, but is !writeable) + m.def("get_cm_const_ref", []() { return Eigen::Ref(get_cm()); }); + m.def("get_rm_const_ref", []() { return Eigen::Ref(get_rm()); }); + // Just the corners (via a Map instead of a Ref): + m.def("get_cm_corners", []() { + auto &x = get_cm(); + return py::EigenDMap( + x.data(), + py::EigenDStride(x.outerStride() * (x.rows() - 1), x.innerStride() * (x.cols() - 1))); + }); + m.def("get_cm_corners_const", []() { + const auto &x = get_cm(); + return py::EigenDMap( + x.data(), + py::EigenDStride(x.outerStride() * (x.rows() - 1), x.innerStride() * (x.cols() - 1))); + }); + + m.def("reset_refs", reset_refs); // Restores get_{cm,rm}_ref to original values + + // Increments and returns ref to (same) matrix + m.def("incr_matrix", [](Eigen::Ref m, double v) { + m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v); + return m; + }, py::return_value_policy::reference); + + // Same, but accepts a matrix of any strides + m.def("incr_matrix_any", [](py::EigenDRef m, double v) { + m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v); + return m; + }, py::return_value_policy::reference); + + // Returns an eigen slice of even rows + m.def("even_rows", [](py::EigenDRef m) { + return py::EigenDMap( + m.data(), (m.rows() + 1) / 2, m.cols(), + py::EigenDStride(m.outerStride(), 2 * m.innerStride())); + }, py::return_value_policy::reference); + + // Returns an eigen slice of even columns + m.def("even_cols", [](py::EigenDRef m) { + return py::EigenDMap( + m.data(), m.rows(), (m.cols() + 1) / 2, + py::EigenDStride(2 * m.outerStride(), m.innerStride())); + }, py::return_value_policy::reference); // Returns diagonals: a vector-like object with an inner stride != 1 m.def("diagonal", [](const Eigen::Ref &x) { return x.diagonal(); }); @@ -72,6 +153,52 @@ test_initializer eigen([](py::module &m) { return x.block(start_row, start_col, block_rows, block_cols); }); + // return value referencing/copying tests: + class ReturnTester { + Eigen::MatrixXd mat = create(); + public: + ReturnTester() { print_created(this); } + ~ReturnTester() { print_destroyed(this); } + static Eigen::MatrixXd create() { return Eigen::MatrixXd::Ones(10, 10); } + static const Eigen::MatrixXd createConst() { return Eigen::MatrixXd::Ones(10, 10); } + Eigen::MatrixXd &get() { return mat; } + Eigen::MatrixXd *getPtr() { return &mat; } + const Eigen::MatrixXd &view() { return mat; } + const Eigen::MatrixXd *viewPtr() { return &mat; } + Eigen::Ref ref() { return mat; } + Eigen::Ref refConst() { return mat; } + Eigen::Block block(int r, int c, int nrow, int ncol) { return mat.block(r, c, nrow, ncol); } + Eigen::Block blockConst(int r, int c, int nrow, int ncol) const { return mat.block(r, c, nrow, ncol); } + py::EigenDMap corners() { return py::EigenDMap(mat.data(), + py::EigenDStride(mat.outerStride() * (mat.outerSize()-1), mat.innerStride() * (mat.innerSize()-1))); } + py::EigenDMap cornersConst() const { return py::EigenDMap(mat.data(), + py::EigenDStride(mat.outerStride() * (mat.outerSize()-1), mat.innerStride() * (mat.innerSize()-1))); } + }; + using rvp = py::return_value_policy; + py::class_(m, "ReturnTester") + .def(py::init<>()) + .def_static("create", &ReturnTester::create) + .def_static("create_const", &ReturnTester::createConst) + .def("get", &ReturnTester::get, rvp::reference_internal) + .def("get_ptr", &ReturnTester::getPtr, rvp::reference_internal) + .def("view", &ReturnTester::view, rvp::reference_internal) + .def("view_ptr", &ReturnTester::view, rvp::reference_internal) + .def("copy_get", &ReturnTester::get) // Default rvp: copy + .def("copy_view", &ReturnTester::view) // " + .def("ref", &ReturnTester::ref) // Default for Ref is to reference + .def("ref_const", &ReturnTester::refConst) // Likewise, but const + .def("ref_safe", &ReturnTester::ref, rvp::reference_internal) + .def("ref_const_safe", &ReturnTester::refConst, rvp::reference_internal) + .def("copy_ref", &ReturnTester::ref, rvp::copy) + .def("copy_ref_const", &ReturnTester::refConst, rvp::copy) + .def("block", &ReturnTester::block) + .def("block_safe", &ReturnTester::block, rvp::reference_internal) + .def("block_const", &ReturnTester::blockConst, rvp::reference_internal) + .def("copy_block", &ReturnTester::block, rvp::copy) + .def("corners", &ReturnTester::corners, rvp::reference_internal) + .def("corners_const", &ReturnTester::cornersConst, rvp::reference_internal) + ; + // Returns a DiagonalMatrix with diagonal (1,2,3,...) m.def("incr_diag", [](int k) { Eigen::DiagonalMatrix m(k); @@ -88,64 +215,49 @@ test_initializer eigen([](py::module &m) { return m.selfadjointView(); }); - m.def("fixed_r", [mat]() -> FixedMatrixR { - return FixedMatrixR(mat); - }); + // Test matrix for various functions below. + Eigen::MatrixXf mat(5, 6); + mat << 0, 3, 0, 0, 0, 11, + 22, 0, 0, 0, 17, 11, + 7, 5, 0, 1, 0, 11, + 0, 0, 0, 0, 0, 11, + 0, 0, 14, 0, 8, 11; - m.def("fixed_c", [mat]() -> FixedMatrixC { - return FixedMatrixC(mat); - }); + m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); }); + m.def("fixed_r_const", [mat]() -> const FixedMatrixR { return FixedMatrixR(mat); }); + m.def("fixed_c", [mat]() -> FixedMatrixC { return FixedMatrixC(mat); }); + m.def("fixed_copy_r", [](const FixedMatrixR &m) -> FixedMatrixR { return m; }); + m.def("fixed_copy_c", [](const FixedMatrixC &m) -> FixedMatrixC { return m; }); + m.def("fixed_mutator_r", [](Eigen::Ref) {}); + m.def("fixed_mutator_c", [](Eigen::Ref) {}); + m.def("fixed_mutator_a", [](py::EigenDRef) {}); + m.def("dense_r", [mat]() -> DenseMatrixR { return DenseMatrixR(mat); }); + m.def("dense_c", [mat]() -> DenseMatrixC { return DenseMatrixC(mat); }); + m.def("dense_copy_r", [](const DenseMatrixR &m) -> DenseMatrixR { return m; }); + m.def("dense_copy_c", [](const DenseMatrixC &m) -> DenseMatrixC { return m; }); + m.def("sparse_r", [mat]() -> SparseMatrixR { return Eigen::SparseView(mat); }); + m.def("sparse_c", [mat]() -> SparseMatrixC { return Eigen::SparseView(mat); }); + m.def("sparse_copy_r", [](const SparseMatrixR &m) -> SparseMatrixR { return m; }); + m.def("sparse_copy_c", [](const SparseMatrixC &m) -> SparseMatrixC { return m; }); + m.def("partial_copy_four_rm_r", [](const FourRowMatrixR &m) -> FourRowMatrixR { return m; }); + m.def("partial_copy_four_rm_c", [](const FourColMatrixR &m) -> FourColMatrixR { return m; }); + m.def("partial_copy_four_cm_r", [](const FourRowMatrixC &m) -> FourRowMatrixC { return m; }); + m.def("partial_copy_four_cm_c", [](const FourColMatrixC &m) -> FourColMatrixC { return m; }); - m.def("fixed_passthrough_r", [](const FixedMatrixR &m) -> FixedMatrixR { - return m; - }); + // Test that we can cast a numpy object to a Eigen::MatrixXd explicitly + m.def("cpp_copy", [](py::handle m) { return m.cast()(1, 0); }); + m.def("cpp_ref_c", [](py::handle m) { return m.cast>()(1, 0); }); + m.def("cpp_ref_r", [](py::handle m) { return m.cast>()(1, 0); }); + m.def("cpp_ref_any", [](py::handle m) { return m.cast>()(1, 0); }); - m.def("fixed_passthrough_c", [](const FixedMatrixC &m) -> FixedMatrixC { - return m; - }); - m.def("dense_r", [mat]() -> DenseMatrixR { - return DenseMatrixR(mat); - }); - - m.def("dense_c", [mat]() -> DenseMatrixC { - return DenseMatrixC(mat); - }); - - m.def("dense_passthrough_r", [](const DenseMatrixR &m) -> DenseMatrixR { - return m; - }); - - m.def("dense_passthrough_c", [](const DenseMatrixC &m) -> DenseMatrixC { - return m; - }); - - m.def("sparse_r", [mat]() -> SparseMatrixR { - return Eigen::SparseView(mat); - }); - - m.def("sparse_c", [mat]() -> SparseMatrixC { - return Eigen::SparseView(mat); - }); - - m.def("sparse_passthrough_r", [](const SparseMatrixR &m) -> SparseMatrixR { - return m; - }); - - m.def("sparse_passthrough_c", [](const SparseMatrixC &m) -> SparseMatrixC { - return m; - }); - - m.def("partial_passthrough_four_rm_r", [](const FourRowMatrixR &m) -> FourRowMatrixR { - return m; - }); - m.def("partial_passthrough_four_rm_c", [](const FourColMatrixR &m) -> FourColMatrixR { - return m; - }); - m.def("partial_passthrough_four_cm_r", [](const FourRowMatrixC &m) -> FourRowMatrixC { - return m; - }); - m.def("partial_passthrough_four_cm_c", [](const FourColMatrixC &m) -> FourColMatrixC { - return m; - }); + // Test that we can prevent copying into an argument that would normally copy: First a version + // that would allow copying (if types or strides don't match) for comparison: + m.def("get_elem", &get_elem); + // Now this alternative that calls the tells pybind to fail rather than copy: + m.def("get_elem_nocopy", [](Eigen::Ref m) -> double { return get_elem(m); }, + py::arg().noconvert()); + // Also test a row-major-only no-copy const ref: + m.def("get_elem_rm_nocopy", [](Eigen::Ref> &m) -> long { return m(2, 1); }, + py::arg().noconvert()); }); diff --git a/tests/test_eigen.py b/tests/test_eigen.py index 5d4d94a52..f3ac9f02b 100644 --- a/tests/test_eigen.py +++ b/tests/test_eigen.py @@ -3,7 +3,7 @@ import pytest with pytest.suppress(ImportError): import numpy as np - ref = np.array([[ 0, 3, 0, 0, 0, 11], + ref = np.array([[ 0., 3, 0, 0, 0, 11], [22, 0, 0, 0, 17, 11], [ 7, 5, 0, 1, 0, 11], [ 0, 0, 0, 0, 0, 11], @@ -20,57 +20,130 @@ def assert_sparse_equal_ref(sparse_mat): @pytest.requires_eigen_and_numpy def test_fixed(): - from pybind11_tests import fixed_r, fixed_c, fixed_passthrough_r, fixed_passthrough_c + from pybind11_tests import fixed_r, fixed_c, fixed_copy_r, fixed_copy_c assert_equal_ref(fixed_c()) assert_equal_ref(fixed_r()) - assert_equal_ref(fixed_passthrough_r(fixed_r())) - assert_equal_ref(fixed_passthrough_c(fixed_c())) - assert_equal_ref(fixed_passthrough_r(fixed_c())) - assert_equal_ref(fixed_passthrough_c(fixed_r())) + assert_equal_ref(fixed_copy_r(fixed_r())) + assert_equal_ref(fixed_copy_c(fixed_c())) + assert_equal_ref(fixed_copy_r(fixed_c())) + assert_equal_ref(fixed_copy_c(fixed_r())) @pytest.requires_eigen_and_numpy def test_dense(): - from pybind11_tests import dense_r, dense_c, dense_passthrough_r, dense_passthrough_c + from pybind11_tests import dense_r, dense_c, dense_copy_r, dense_copy_c assert_equal_ref(dense_r()) assert_equal_ref(dense_c()) - assert_equal_ref(dense_passthrough_r(dense_r())) - assert_equal_ref(dense_passthrough_c(dense_c())) - assert_equal_ref(dense_passthrough_r(dense_c())) - assert_equal_ref(dense_passthrough_c(dense_r())) + assert_equal_ref(dense_copy_r(dense_r())) + assert_equal_ref(dense_copy_c(dense_c())) + assert_equal_ref(dense_copy_r(dense_c())) + assert_equal_ref(dense_copy_c(dense_r())) + @pytest.requires_eigen_and_numpy def test_partially_fixed(): - from pybind11_tests import partial_passthrough_four_rm_r, partial_passthrough_four_rm_c, partial_passthrough_four_cm_r, partial_passthrough_four_cm_c + from pybind11_tests import (partial_copy_four_rm_r, partial_copy_four_rm_c, + partial_copy_four_cm_r, partial_copy_four_cm_c) - ref2 = np.array([[0,1,2,3], [4,5,6,7], [8,9,10,11], [12,13,14,15]]) - np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2), ref2) - np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2), ref2) - np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2[:, 1]), ref2[:, [1]]) - np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2[0, :]), ref2[[0], :]) - np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2[:, (0, 2)]), ref2[:, (0,2)]) - np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2[(3,1,2), :]), ref2[(3,1,2), :]) + ref2 = np.array([[0., 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]]) + np.testing.assert_array_equal(partial_copy_four_rm_r(ref2), ref2) + np.testing.assert_array_equal(partial_copy_four_rm_c(ref2), ref2) + np.testing.assert_array_equal(partial_copy_four_rm_r(ref2[:, 1]), ref2[:, [1]]) + np.testing.assert_array_equal(partial_copy_four_rm_c(ref2[0, :]), ref2[[0], :]) + np.testing.assert_array_equal(partial_copy_four_rm_r(ref2[:, (0, 2)]), ref2[:, (0, 2)]) + np.testing.assert_array_equal( + partial_copy_four_rm_c(ref2[(3, 1, 2), :]), ref2[(3, 1, 2), :]) + + np.testing.assert_array_equal(partial_copy_four_cm_r(ref2), ref2) + np.testing.assert_array_equal(partial_copy_four_cm_c(ref2), ref2) + np.testing.assert_array_equal(partial_copy_four_cm_r(ref2[:, 1]), ref2[:, [1]]) + np.testing.assert_array_equal(partial_copy_four_cm_c(ref2[0, :]), ref2[[0], :]) + np.testing.assert_array_equal(partial_copy_four_cm_r(ref2[:, (0, 2)]), ref2[:, (0, 2)]) + np.testing.assert_array_equal( + partial_copy_four_cm_c(ref2[(3, 1, 2), :]), ref2[(3, 1, 2), :]) + + +@pytest.requires_eigen_and_numpy +def test_mutator_descriptors(): + from pybind11_tests import fixed_mutator_r, fixed_mutator_c, fixed_mutator_a + zr = np.arange(30, dtype='float32').reshape(5, 6) # row-major + zc = zr.reshape(6, 5).transpose() # column-major + + fixed_mutator_r(zr) + fixed_mutator_c(zc) + fixed_mutator_a(zr) + fixed_mutator_a(zc) + with pytest.raises(TypeError) as excinfo: + fixed_mutator_r(zc) + assert ('(numpy.ndarray[float32[5, 6], flags.writeable, flags.c_contiguous]) -> arg0: None' + in str(excinfo.value)) + with pytest.raises(TypeError) as excinfo: + fixed_mutator_c(zr) + assert ('(numpy.ndarray[float32[5, 6], flags.writeable, flags.f_contiguous]) -> arg0: None' + in str(excinfo.value)) + with pytest.raises(TypeError) as excinfo: + fixed_mutator_a(np.array([[1, 2], [3, 4]], dtype='float32')) + assert ('(numpy.ndarray[float32[5, 6], flags.writeable]) -> arg0: None' + in str(excinfo.value)) + zr.flags.writeable = False + with pytest.raises(TypeError): + fixed_mutator_r(zr) + with pytest.raises(TypeError): + fixed_mutator_a(zr) + + +@pytest.requires_eigen_and_numpy +def test_cpp_casting(): + from pybind11_tests import (cpp_copy, cpp_ref_c, cpp_ref_r, cpp_ref_any, + fixed_r, fixed_c, get_cm_ref, get_rm_ref, ReturnTester) + assert cpp_copy(fixed_r()) == 22. + assert cpp_copy(fixed_c()) == 22. + z = np.array([[5., 6], [7, 8]]) + assert cpp_copy(z) == 7. + assert cpp_copy(get_cm_ref()) == 21. + assert cpp_copy(get_rm_ref()) == 21. + assert cpp_ref_c(get_cm_ref()) == 21. + assert cpp_ref_r(get_rm_ref()) == 21. + with pytest.raises(RuntimeError) as excinfo: + # Can't reference fixed_c: it contains floats, cpp_ref_any wants doubles + cpp_ref_any(fixed_c()) + assert 'Unable to cast Python instance' in str(excinfo.value) + with pytest.raises(RuntimeError) as excinfo: + # Can't reference fixed_r: it contains floats, cpp_ref_any wants doubles + cpp_ref_any(fixed_r()) + assert 'Unable to cast Python instance' in str(excinfo.value) + assert cpp_ref_any(ReturnTester.create()) == 1. + + assert cpp_ref_any(get_cm_ref()) == 21. + assert cpp_ref_any(get_cm_ref()) == 21. + + +@pytest.requires_eigen_and_numpy +def test_pass_readonly_array(): + from pybind11_tests import fixed_copy_r, fixed_r, fixed_r_const + z = np.full((5, 6), 42.0) + z.flags.writeable = False + np.testing.assert_array_equal(z, fixed_copy_r(z)) + np.testing.assert_array_equal(fixed_r_const(), fixed_r()) + assert not fixed_r_const().flags.writeable + np.testing.assert_array_equal(fixed_copy_r(fixed_r_const()), fixed_r_const()) - np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2), ref2) - np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2), ref2) - np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2[:, 1]), ref2[:, [1]]) - np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2[0, :]), ref2[[0], :]) - np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2[:, (0, 2)]), ref2[:, (0,2)]) - np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2[(3,1,2), :]), ref2[(3,1,2), :]) @pytest.requires_eigen_and_numpy def test_nonunit_stride_from_python(): - from pybind11_tests import double_row, double_col, double_mat_cm, double_mat_rm + from pybind11_tests import ( + double_row, double_col, double_mat_cm, double_mat_rm, + double_threec, double_threer) counting_mat = np.arange(9.0, dtype=np.float32).reshape((3, 3)) - first_row = counting_mat[0, :] - first_col = counting_mat[:, 0] - np.testing.assert_array_equal(double_row(first_row), 2.0 * first_row) - np.testing.assert_array_equal(double_col(first_row), 2.0 * first_row) - np.testing.assert_array_equal(double_row(first_col), 2.0 * first_col) - np.testing.assert_array_equal(double_col(first_col), 2.0 * first_col) + 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_row(second_col), 2.0 * second_col) + np.testing.assert_array_equal(double_col(second_col), 2.0 * second_col) counting_3d = np.arange(27.0, dtype=np.float32).reshape((3, 3, 3)) slices = [counting_3d[0, :, :], counting_3d[:, 0, :], counting_3d[:, :, 0]] @@ -78,6 +151,11 @@ def test_nonunit_stride_from_python(): 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: + double_threer(second_row) + double_threec(second_col) + np.testing.assert_array_equal(counting_mat, [[0., 2, 2], [6, 16, 10], [6, 14, 8]]) + @pytest.requires_eigen_and_numpy def test_nonunit_stride_to_python(): @@ -95,21 +173,396 @@ def test_nonunit_stride_to_python(): @pytest.requires_eigen_and_numpy def test_eigen_ref_to_python(): - from pybind11_tests import cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6 + from pybind11_tests import cholesky1, cholesky2, cholesky3, cholesky4 - chols = [cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6] + chols = [cholesky1, cholesky2, cholesky3, cholesky4] for i, chol in enumerate(chols, start=1): - mymat = chol(np.array([[1, 2, 4], [2, 13, 23], [4, 23, 77]])) + mymat = chol(np.array([[1., 2, 4], [2, 13, 23], [4, 23, 77]])) assert np.all(mymat == np.array([[1, 0, 0], [2, 3, 0], [4, 5, 6]])), "cholesky{}".format(i) +def assign_both(a1, a2, r, c, v): + a1[r, c] = v + a2[r, c] = v + + +def array_copy_but_one(a, r, c, v): + z = np.array(a, copy=True) + z[r, c] = v + return z + + +@pytest.requires_eigen_and_numpy +def test_eigen_return_references(): + """Tests various ways of returning references and non-referencing copies""" + from pybind11_tests import ReturnTester + master = np.ones((10, 10)) + a = ReturnTester() + a_get1 = a.get() + assert not a_get1.flags.owndata and a_get1.flags.writeable + assign_both(a_get1, master, 3, 3, 5) + a_get2 = a.get_ptr() + assert not a_get2.flags.owndata and a_get2.flags.writeable + assign_both(a_get1, master, 2, 3, 6) + + a_view1 = a.view() + assert not a_view1.flags.owndata and not a_view1.flags.writeable + with pytest.raises(ValueError): + a_view1[2, 3] = 4 + a_view2 = a.view_ptr() + assert not a_view2.flags.owndata and not a_view2.flags.writeable + with pytest.raises(ValueError): + a_view2[2, 3] = 4 + + a_copy1 = a.copy_get() + assert a_copy1.flags.owndata and a_copy1.flags.writeable + np.testing.assert_array_equal(a_copy1, master) + a_copy1[7, 7] = -44 # Shouldn't affect anything else + c1want = array_copy_but_one(master, 7, 7, -44) + a_copy2 = a.copy_view() + assert a_copy2.flags.owndata and a_copy2.flags.writeable + np.testing.assert_array_equal(a_copy2, master) + a_copy2[4, 4] = -22 # Shouldn't affect anything else + c2want = array_copy_but_one(master, 4, 4, -22) + + a_ref1 = a.ref() + assert not a_ref1.flags.owndata and a_ref1.flags.writeable + assign_both(a_ref1, master, 1, 1, 15) + a_ref2 = a.ref_const() + assert not a_ref2.flags.owndata and not a_ref2.flags.writeable + with pytest.raises(ValueError): + a_ref2[5, 5] = 33 + a_ref3 = a.ref_safe() + assert not a_ref3.flags.owndata and a_ref3.flags.writeable + assign_both(a_ref3, master, 0, 7, 99) + a_ref4 = a.ref_const_safe() + assert not a_ref4.flags.owndata and not a_ref4.flags.writeable + with pytest.raises(ValueError): + a_ref4[7, 0] = 987654321 + + a_copy3 = a.copy_ref() + assert a_copy3.flags.owndata and a_copy3.flags.writeable + np.testing.assert_array_equal(a_copy3, master) + a_copy3[8, 1] = 11 + c3want = array_copy_but_one(master, 8, 1, 11) + a_copy4 = a.copy_ref_const() + assert a_copy4.flags.owndata and a_copy4.flags.writeable + np.testing.assert_array_equal(a_copy4, master) + a_copy4[8, 4] = 88 + c4want = array_copy_but_one(master, 8, 4, 88) + + a_block1 = a.block(3, 3, 2, 2) + assert not a_block1.flags.owndata and a_block1.flags.writeable + a_block1[0, 0] = 55 + master[3, 3] = 55 + a_block2 = a.block_safe(2, 2, 3, 2) + assert not a_block2.flags.owndata and a_block2.flags.writeable + a_block2[2, 1] = -123 + master[4, 3] = -123 + a_block3 = a.block_const(6, 7, 4, 3) + assert not a_block3.flags.owndata and not a_block3.flags.writeable + with pytest.raises(ValueError): + a_block3[2, 2] = -44444 + + a_copy5 = a.copy_block(2, 2, 2, 3) + assert a_copy5.flags.owndata and a_copy5.flags.writeable + np.testing.assert_array_equal(a_copy5, master[2:4, 2:5]) + a_copy5[1, 1] = 777 + c5want = array_copy_but_one(master[2:4, 2:5], 1, 1, 777) + + a_corn1 = a.corners() + assert not a_corn1.flags.owndata and a_corn1.flags.writeable + a_corn1 *= 50 + a_corn1[1, 1] = 999 + master[0, 0] = 50 + master[0, 9] = 50 + master[9, 0] = 50 + master[9, 9] = 999 + a_corn2 = a.corners_const() + assert not a_corn2.flags.owndata and not a_corn2.flags.writeable + with pytest.raises(ValueError): + a_corn2[1, 0] = 51 + + # All of the changes made all the way along should be visible everywhere + # now (except for the copies, of course) + np.testing.assert_array_equal(a_get1, master) + np.testing.assert_array_equal(a_get2, master) + np.testing.assert_array_equal(a_view1, master) + np.testing.assert_array_equal(a_view2, master) + np.testing.assert_array_equal(a_ref1, master) + np.testing.assert_array_equal(a_ref2, master) + np.testing.assert_array_equal(a_ref3, master) + np.testing.assert_array_equal(a_ref4, master) + np.testing.assert_array_equal(a_block1, master[3:5, 3:5]) + np.testing.assert_array_equal(a_block2, master[2:5, 2:4]) + np.testing.assert_array_equal(a_block3, master[6:10, 7:10]) + np.testing.assert_array_equal(a_corn1, master[0::master.shape[0] - 1, 0::master.shape[1] - 1]) + np.testing.assert_array_equal(a_corn2, master[0::master.shape[0] - 1, 0::master.shape[1] - 1]) + + np.testing.assert_array_equal(a_copy1, c1want) + np.testing.assert_array_equal(a_copy2, c2want) + np.testing.assert_array_equal(a_copy3, c3want) + np.testing.assert_array_equal(a_copy4, c4want) + np.testing.assert_array_equal(a_copy5, c5want) + + +def assert_keeps_alive(cl, method, *args): + from pybind11_tests import ConstructorStats + cstats = ConstructorStats.get(cl) + start_with = cstats.alive() + a = cl() + assert cstats.alive() == start_with + 1 + z = method(a, *args) + assert cstats.alive() == start_with + 1 + del a + # Here's the keep alive in action: + assert cstats.alive() == start_with + 1 + del z + # Keep alive should have expired: + assert cstats.alive() == start_with + + +@pytest.requires_eigen_and_numpy +def test_eigen_keepalive(): + from pybind11_tests import ReturnTester, ConstructorStats + a = ReturnTester() + + cstats = ConstructorStats.get(ReturnTester) + assert cstats.alive() == 1 + unsafe = [a.ref(), a.ref_const(), a.block(1, 2, 3, 4)] + copies = [a.copy_get(), a.copy_view(), a.copy_ref(), a.copy_ref_const(), + a.copy_block(4, 3, 2, 1)] + del a + assert cstats.alive() == 0 + del unsafe + del copies + + for meth in [ReturnTester.get, ReturnTester.get_ptr, ReturnTester.view, + ReturnTester.view_ptr, ReturnTester.ref_safe, ReturnTester.ref_const_safe, + ReturnTester.corners, ReturnTester.corners_const]: + assert_keeps_alive(ReturnTester, meth) + + for meth in [ReturnTester.block_safe, ReturnTester.block_const]: + assert_keeps_alive(ReturnTester, meth, 4, 3, 2, 1) + + +@pytest.requires_eigen_and_numpy +def test_eigen_ref_mutators(): + """Tests whether Eigen can mutate numpy values""" + from pybind11_tests import add_rm, add_cm, add_any, add1, add2 + orig = np.array([[1., 2, 3], [4, 5, 6], [7, 8, 9]]) + zr = np.array(orig) + zc = np.array(orig, order='F') + add_rm(zr, 1, 0, 100) + assert np.all(zr == np.array([[1., 2, 3], [104, 5, 6], [7, 8, 9]])) + add_cm(zc, 1, 0, 200) + assert np.all(zc == np.array([[1., 2, 3], [204, 5, 6], [7, 8, 9]])) + + add_any(zr, 1, 0, 20) + assert np.all(zr == np.array([[1., 2, 3], [124, 5, 6], [7, 8, 9]])) + add_any(zc, 1, 0, 10) + assert np.all(zc == np.array([[1., 2, 3], [214, 5, 6], [7, 8, 9]])) + + # Can't reference a col-major array with a row-major Ref, and vice versa: + with pytest.raises(TypeError): + add_rm(zc, 1, 0, 1) + with pytest.raises(TypeError): + add_cm(zr, 1, 0, 1) + + # Overloads: + add1(zr, 1, 0, -100) + add2(zr, 1, 0, -20) + assert np.all(zr == orig) + add1(zc, 1, 0, -200) + add2(zc, 1, 0, -10) + assert np.all(zc == orig) + + # a non-contiguous slice (this won't work on either the row- or + # column-contiguous refs, but should work for the any) + cornersr = zr[0::2, 0::2] + cornersc = zc[0::2, 0::2] + + assert np.all(cornersr == np.array([[1., 3], [7, 9]])) + assert np.all(cornersc == np.array([[1., 3], [7, 9]])) + + with pytest.raises(TypeError): + add_rm(cornersr, 0, 1, 25) + with pytest.raises(TypeError): + add_cm(cornersr, 0, 1, 25) + with pytest.raises(TypeError): + add_rm(cornersc, 0, 1, 25) + with pytest.raises(TypeError): + add_cm(cornersc, 0, 1, 25) + add_any(cornersr, 0, 1, 25) + add_any(cornersc, 0, 1, 44) + assert np.all(zr == np.array([[1., 2, 28], [4, 5, 6], [7, 8, 9]])) + assert np.all(zc == np.array([[1., 2, 47], [4, 5, 6], [7, 8, 9]])) + + # You shouldn't be allowed to pass a non-writeable array to a mutating Eigen method: + zro = zr[0:4, 0:4] + zro.flags.writeable = False + with pytest.raises(TypeError): + add_rm(zro, 0, 0, 0) + with pytest.raises(TypeError): + add_any(zro, 0, 0, 0) + with pytest.raises(TypeError): + add1(zro, 0, 0, 0) + with pytest.raises(TypeError): + add2(zro, 0, 0, 0) + + # integer array shouldn't be passable to a double-matrix-accepting mutating func: + zi = np.array([[1, 2], [3, 4]]) + with pytest.raises(TypeError): + add_rm(zi) + + +@pytest.requires_eigen_and_numpy +def test_numpy_ref_mutators(): + """Tests numpy mutating Eigen matrices (for returned Eigen::Ref<...>s)""" + from pybind11_tests import ( + get_cm_ref, get_cm_const_ref, get_rm_ref, get_rm_const_ref, reset_refs) + reset_refs() # In case another test already changed it + + zc = get_cm_ref() + zcro = get_cm_const_ref() + zr = get_rm_ref() + zrro = get_rm_const_ref() + + assert [zc[1, 2], zcro[1, 2], zr[1, 2], zrro[1, 2]] == [23] * 4 + + assert not zc.flags.owndata and zc.flags.writeable + assert not zr.flags.owndata and zr.flags.writeable + assert not zcro.flags.owndata and not zcro.flags.writeable + assert not zrro.flags.owndata and not zrro.flags.writeable + + zc[1, 2] = 99 + expect = np.array([[11., 12, 13], [21, 22, 99], [31, 32, 33]]) + # We should have just changed zc, of course, but also zcro and the original eigen matrix + assert np.all(zc == expect) + assert np.all(zcro == expect) + assert np.all(get_cm_ref() == expect) + + zr[1, 2] = 99 + assert np.all(zr == expect) + assert np.all(zrro == expect) + assert np.all(get_rm_ref() == expect) + + # Make sure the readonly ones are numpy-readonly: + with pytest.raises(ValueError): + zcro[1, 2] = 6 + with pytest.raises(ValueError): + zrro[1, 2] = 6 + + # We should be able to explicitly copy like this (and since we're copying, + # the const should drop away) + y1 = np.array(get_cm_const_ref()) + + assert y1.flags.owndata and y1.flags.writeable + # We should get copies of the eigen data, which was modified above: + assert y1[1, 2] == 99 + y1[1, 2] += 12 + assert y1[1, 2] == 111 + assert zc[1, 2] == 99 # Make sure we aren't referencing the original + + +@pytest.requires_eigen_and_numpy +def test_both_ref_mutators(): + """Tests a complex chain of nested eigen/numpy references""" + from pybind11_tests import ( + incr_matrix, get_cm_ref, incr_matrix_any, even_cols, even_rows, reset_refs) + reset_refs() # In case another test already changed it + + z = get_cm_ref() # numpy -> eigen + z[0, 2] -= 3 + z2 = incr_matrix(z, 1) # numpy -> eigen -> numpy -> eigen + z2[1, 1] += 6 + z3 = incr_matrix(z, 2) # (numpy -> eigen)^3 + z3[2, 2] += -5 + z4 = incr_matrix(z, 3) # (numpy -> eigen)^4 + z4[1, 1] -= 1 + z5 = incr_matrix(z, 4) # (numpy -> eigen)^5 + z5[0, 0] = 0 + assert np.all(z == z2) + assert np.all(z == z3) + assert np.all(z == z4) + assert np.all(z == z5) + expect = np.array([[0., 22, 20], [31, 37, 33], [41, 42, 38]]) + assert np.all(z == expect) + + y = np.array(range(100), dtype='float64').reshape(10, 10) + y2 = incr_matrix_any(y, 10) # np -> eigen -> np + y3 = incr_matrix_any(y2[0::2, 0::2], -33) # np -> eigen -> np slice -> np -> eigen -> np + y4 = even_rows(y3) # numpy -> eigen slice -> (... y3) + y5 = even_cols(y4) # numpy -> eigen slice -> (... y4) + y6 = incr_matrix_any(y5, 1000) # numpy -> eigen -> (... y5) + + # Apply same mutations using just numpy: + yexpect = np.array(range(100), dtype='float64').reshape(10, 10) + yexpect += 10 + yexpect[0::2, 0::2] -= 33 + yexpect[0::4, 0::4] += 1000 + assert np.all(y6 == yexpect[0::4, 0::4]) + assert np.all(y5 == yexpect[0::4, 0::4]) + assert np.all(y4 == yexpect[0::4, 0::2]) + assert np.all(y3 == yexpect[0::2, 0::2]) + assert np.all(y2 == yexpect) + assert np.all(y == yexpect) + + +@pytest.requires_eigen_and_numpy +def test_nocopy_wrapper(): + from pybind11_tests import get_elem, get_elem_nocopy, get_elem_rm_nocopy + # get_elem requires a column-contiguous matrix reference, but should be + # callable with other types of matrix (via copying): + int_matrix_colmajor = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], order='F') + dbl_matrix_colmajor = np.array(int_matrix_colmajor, dtype='double', order='F', copy=True) + int_matrix_rowmajor = np.array(int_matrix_colmajor, order='C', copy=True) + dbl_matrix_rowmajor = np.array(int_matrix_rowmajor, dtype='double', order='C', copy=True) + + # All should be callable via get_elem: + assert get_elem(int_matrix_colmajor) == 8 + assert get_elem(dbl_matrix_colmajor) == 8 + assert get_elem(int_matrix_rowmajor) == 8 + assert get_elem(dbl_matrix_rowmajor) == 8 + + # All but the second should fail with get_elem_nocopy: + with pytest.raises(TypeError) as excinfo: + get_elem_nocopy(int_matrix_colmajor) + assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.f_contiguous' in str(excinfo.value)) + assert get_elem_nocopy(dbl_matrix_colmajor) == 8 + with pytest.raises(TypeError) as excinfo: + get_elem_nocopy(int_matrix_rowmajor) + assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.f_contiguous' in str(excinfo.value)) + with pytest.raises(TypeError) as excinfo: + get_elem_nocopy(dbl_matrix_rowmajor) + assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.f_contiguous' in str(excinfo.value)) + + # For the row-major test, we take a long matrix in row-major, so only the third is allowed: + with pytest.raises(TypeError) as excinfo: + get_elem_rm_nocopy(int_matrix_colmajor) + assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.c_contiguous' in str(excinfo.value)) + with pytest.raises(TypeError) as excinfo: + get_elem_rm_nocopy(dbl_matrix_colmajor) + assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.c_contiguous' in str(excinfo.value)) + assert get_elem_rm_nocopy(int_matrix_rowmajor) == 8 + with pytest.raises(TypeError) as excinfo: + get_elem_rm_nocopy(dbl_matrix_rowmajor) + assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and + ', flags.c_contiguous' in str(excinfo.value)) + + @pytest.requires_eigen_and_numpy def test_special_matrix_objects(): from pybind11_tests import incr_diag, symmetric_upper, symmetric_lower - assert np.all(incr_diag(7) == np.diag([1, 2, 3, 4, 5, 6, 7])) + assert np.all(incr_diag(7) == np.diag([1., 2, 3, 4, 5, 6, 7])) - asymm = np.array([[ 1, 2, 3, 4], + asymm = np.array([[ 1., 2, 3, 4], [ 5, 6, 7, 8], [ 9, 10, 11, 12], [13, 14, 15, 16]]) @@ -141,23 +594,23 @@ def test_dense_signature(doc): @pytest.requires_eigen_and_scipy def test_sparse(): - from pybind11_tests import sparse_r, sparse_c, sparse_passthrough_r, sparse_passthrough_c + from pybind11_tests import sparse_r, sparse_c, sparse_copy_r, sparse_copy_c assert_sparse_equal_ref(sparse_r()) assert_sparse_equal_ref(sparse_c()) - assert_sparse_equal_ref(sparse_passthrough_r(sparse_r())) - assert_sparse_equal_ref(sparse_passthrough_c(sparse_c())) - assert_sparse_equal_ref(sparse_passthrough_r(sparse_c())) - assert_sparse_equal_ref(sparse_passthrough_c(sparse_r())) + assert_sparse_equal_ref(sparse_copy_r(sparse_r())) + assert_sparse_equal_ref(sparse_copy_c(sparse_c())) + assert_sparse_equal_ref(sparse_copy_r(sparse_c())) + assert_sparse_equal_ref(sparse_copy_c(sparse_r())) @pytest.requires_eigen_and_scipy def test_sparse_signature(doc): - from pybind11_tests import sparse_passthrough_r, sparse_passthrough_c + from pybind11_tests import sparse_copy_r, sparse_copy_c - assert doc(sparse_passthrough_r) == """ - sparse_passthrough_r(arg0: scipy.sparse.csr_matrix[float32]) -> scipy.sparse.csr_matrix[float32] + assert doc(sparse_copy_r) == """ + sparse_copy_r(arg0: scipy.sparse.csr_matrix[float32]) -> scipy.sparse.csr_matrix[float32] """ # noqa: E501 line too long - assert doc(sparse_passthrough_c) == """ - sparse_passthrough_c(arg0: scipy.sparse.csc_matrix[float32]) -> scipy.sparse.csc_matrix[float32] + assert doc(sparse_copy_c) == """ + sparse_copy_c(arg0: scipy.sparse.csc_matrix[float32]) -> scipy.sparse.csc_matrix[float32] """ # noqa: E501 line too long