From 9ffb3dda5fb9b97a8ca5dea288a7de2e532d408c Mon Sep 17 00:00:00 2001 From: Jason Rhinelander Date: Thu, 4 Aug 2016 15:24:41 -0400 Subject: [PATCH] Eigen support for special matrix objects Functions returning specialized Eigen matrices like Eigen::DiagonalMatrix and Eigen::SelfAdjointView--which inherit from EigenBase but not DenseBase--isn't currently allowed; such classes are explicitly copyable into a Matrix (by definition), and so we can support functions that return them by copying the value into a Matrix then casting that resulting dense Matrix into a numpy.ndarray. This commit does exactly that. --- docs/advanced.rst | 17 ++++++++++++++--- example/eigen.cpp | 16 ++++++++++++++++ example/eigen.py | 18 ++++++++++++++++++ example/eigen.ref | 3 +++ include/pybind11/eigen.h | 35 +++++++++++++++++++++++++++++++++-- 5 files changed, 84 insertions(+), 5 deletions(-) diff --git a/docs/advanced.rst b/docs/advanced.rst index 3ebf83e93..27df2d896 100644 --- a/docs/advanced.rst +++ b/docs/advanced.rst @@ -1098,6 +1098,14 @@ pybind11 will automatically and transparently convert 1. Static and dynamic Eigen dense vectors and matrices to instances of ``numpy.ndarray`` (and vice versa). +1. Returned matrix expressions such as blocks (including columns or rows) and + diagonals will be converted to ``numpy.ndarray`` of the expression + values. + +1. Returned matrix-like objects such as Eigen::DiagonalMatrix or + Eigen::SelfAdjointView will be converted to ``numpy.ndarray`` containing the + expressed value. + 1. Eigen sparse vectors and matrices to instances of ``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` (and vice versa). @@ -1107,10 +1115,13 @@ them somehow, in which case the information won't be propagated to the caller. .. code-block:: cpp - /* The Python bindings of this function won't replicate - the intended effect of modifying the function argument */ + /* 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; + v *= 2; + } + void scale_by_2(Eigen::Ref &v) { + v *= 2; } To see why this is, refer to the section on :ref:`opaque` (although that diff --git a/example/eigen.cpp b/example/eigen.cpp index 9c33c551e..c88bee2bb 100644 --- a/example/eigen.cpp +++ b/example/eigen.cpp @@ -66,6 +66,22 @@ void init_eigen(py::module &m) { return x.block(start_row, start_col, block_rows, block_cols); }); + // Returns a DiagonalMatrix with diagonal (1,2,3,...) + m.def("incr_diag", [](int k) { + Eigen::DiagonalMatrix m(k); + for (int i = 0; i < k; i++) m.diagonal()[i] = i+1; + return m; + }); + + // Returns a SelfAdjointView referencing the lower triangle of m + m.def("symmetric_lower", [](const Eigen::MatrixXi &m) { + return m.selfadjointView(); + }); + // Returns a SelfAdjointView referencing the lower triangle of m + m.def("symmetric_upper", [](const Eigen::MatrixXi &m) { + return m.selfadjointView(); + }); + m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); }); diff --git a/example/eigen.py b/example/eigen.py index 04078b9d9..5f7ec5141 100644 --- a/example/eigen.py +++ b/example/eigen.py @@ -14,6 +14,7 @@ from example import double_mat_cm, double_mat_rm from example import cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6 from example import diagonal, diagonal_1, diagonal_n from example import block +from example import incr_diag, symmetric_upper, symmetric_lower try: import numpy as np import scipy @@ -88,3 +89,20 @@ for i in range(-5, 7): print("block(2,1,3,3) %s" % ("OK" if (block(ref, 2, 1, 3, 3) == ref[2:5, 1:4]).all() else "FAILED")) print("block(1,4,4,2) %s" % ("OK" if (block(ref, 1, 4, 4, 2) == ref[1:, 4:]).all() else "FAILED")) print("block(1,4,3,2) %s" % ("OK" if (block(ref, 1, 4, 3, 2) == ref[1:4, 4:]).all() else "FAILED")) + +print("incr_diag %s" % ("OK" if (incr_diag(7) == np.diag([1,2,3,4,5,6,7])).all() else "FAILED")) + +asymm = np.array([ + [1, 2, 3, 4], + [5, 6, 7, 8], + [9, 10,11,12], + [13,14,15,16]]) +symm_lower = np.array(asymm) +symm_upper = np.array(asymm) +for i in range(4): + for j in range(i+1, 4): + symm_lower[i,j] = symm_lower[j,i] + symm_upper[j,i] = symm_upper[i,j] + +print("symmetric_lower %s" % ("OK" if (symmetric_lower(asymm) == symm_lower).all() else "FAILED")) +print("symmetric_upper %s" % ("OK" if (symmetric_upper(asymm) == symm_upper).all() else "FAILED")) diff --git a/example/eigen.ref b/example/eigen.ref index 8ccd1f472..755012f43 100644 --- a/example/eigen.ref +++ b/example/eigen.ref @@ -50,3 +50,6 @@ diagonal_n(6) OK block(2,1,3,3) OK block(1,4,4,2) OK block(1,4,3,2) OK +incr_diag OK +symmetric_lower OK +symmetric_upper OK diff --git a/include/pybind11/eigen.h b/include/pybind11/eigen.h index 16ccaa845..9c531fa2e 100644 --- a/include/pybind11/eigen.h +++ b/include/pybind11/eigen.h @@ -61,6 +61,19 @@ public: static constexpr bool value = decltype(test(std::declval()))::value; }; +// 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 class is_eigen_base { +private: + template static std::true_type test(const Eigen::EigenBase &); + static std::false_type test(...); +public: + static constexpr bool value = !is_eigen_dense::value && !is_eigen_sparse::value && + decltype(test(std::declval()))::value; +}; + template struct type_caster::value && !is_eigen_ref::value>::type> { typedef typename Type::Scalar Scalar; @@ -164,11 +177,10 @@ protected: template struct type_caster::value && is_eigen_ref::value>::type> { -private: +protected: using Derived = typename std::remove_const::Derived>::type; using DerivedCaster = type_caster; DerivedCaster derived_caster; -protected: std::unique_ptr value; 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; } @@ -182,6 +194,25 @@ public: template using cast_op_type = pybind11::detail::cast_op_type<_T>; }; +// 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. +template +struct type_caster::value && !is_eigen_ref::value>::type> { +protected: + using Matrix = Eigen::Matrix; + using MatrixCaster = type_caster; +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 PYBIND11_DESCR name() { return MatrixCaster::name(); } + + [[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>; +}; + template struct type_caster::value>::type> { typedef typename Type::Scalar Scalar;