pybind11-2.10.3/tests/test_eigen_matrix.cpp

0001 /*
0002     tests/eigen.cpp -- automatic conversion of Eigen types
0003
0004     Copyright (c) 2016 Wenzel Jakob <wenzel.jakob@epfl.ch>
0005
0006     All rights reserved. Use of this source code is governed by a
0007     BSD-style license that can be found in the LICENSE file.
0008 */
0009
0010 #include <pybind11/eigen/matrix.h>
0011 #include <pybind11/stl.h>
0012
0013 #include "constructor_stats.h"
0014 #include "pybind11_tests.h"
0015
0016 PYBIND11_WARNING_DISABLE_MSVC(4996)
0017
0018 #include <Eigen/Cholesky>
0019
0020 using MatrixXdR = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
0021
0022 // Sets/resets a testing reference matrix to have values of 10*r + c, where r and c are the
0023 // (1-based) row/column number.
0024 template <typename M>
0025 void reset_ref(M &x) {
0026     for (int i = 0; i < x.rows(); i++) {
0027         for (int j = 0; j < x.cols(); j++) {
0028             x(i, j) = 11 + 10 * i + j;
0029         }
0030     }
0031 }
0032
0033 // Returns a static, column-major matrix
0034 Eigen::MatrixXd &get_cm() {
0035     static Eigen::MatrixXd *x;
0036     if (!x) {
0037         x = new Eigen::MatrixXd(3, 3);
0038         reset_ref(*x);
0039     }
0040     return *x;
0041 }
0042 // Likewise, but row-major
0043 MatrixXdR &get_rm() {
0044     static MatrixXdR *x;
0045     if (!x) {
0046         x = new MatrixXdR(3, 3);
0047         reset_ref(*x);
0048     }
0049     return *x;
0050 }
0051 // Resets the values of the static matrices returned by get_cm()/get_rm()
0052 void reset_refs() {
0053     reset_ref(get_cm());
0054     reset_ref(get_rm());
0055 }
0056
0057 // Returns element 2,1 from a matrix (used to test copy/nocopy)
0058 double get_elem(const Eigen::Ref<const Eigen::MatrixXd> &m) { return m(2, 1); };
0059
0060 // Returns a matrix with 10*r + 100*c added to each matrix element (to help test that the matrix
0061 // reference is referencing rows/columns correctly).
0062 template <typename MatrixArgType>
0063 Eigen::MatrixXd adjust_matrix(MatrixArgType m) {
0064     Eigen::MatrixXd ret(m);
0065     for (int c = 0; c < m.cols(); c++) {
0066         for (int r = 0; r < m.rows(); r++) {
0067             ret(r, c) += 10 * r + 100 * c; // NOLINT(clang-analyzer-core.uninitialized.Assign)
0068         }
0069     }
0070     return ret;
0071 }
0072
0073 struct CustomOperatorNew {
0074     CustomOperatorNew() = default;
0075
0076     Eigen::Matrix4d a = Eigen::Matrix4d::Zero();
0077     Eigen::Matrix4d b = Eigen::Matrix4d::Identity();
0078
0079     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
0080 };
0081
0082 TEST_SUBMODULE(eigen_matrix, m) {
0083     using FixedMatrixR = Eigen::Matrix<float, 5, 6, Eigen::RowMajor>;
0084     using FixedMatrixC = Eigen::Matrix<float, 5, 6>;
0085     using DenseMatrixR = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
0086     using DenseMatrixC = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic>;
0087     using FourRowMatrixC = Eigen::Matrix<float, 4, Eigen::Dynamic>;
0088     using FourColMatrixC = Eigen::Matrix<float, Eigen::Dynamic, 4>;
0089     using FourRowMatrixR = Eigen::Matrix<float, 4, Eigen::Dynamic>;
0090     using FourColMatrixR = Eigen::Matrix<float, Eigen::Dynamic, 4>;
0091     using SparseMatrixR = Eigen::SparseMatrix<float, Eigen::RowMajor>;
0092     using SparseMatrixC = Eigen::SparseMatrix<float>;
0093
0094     // various tests
0095     m.def("double_col", [](const Eigen::VectorXf &x) -> Eigen::VectorXf { return 2.0f * x; });
0096     m.def("double_row",
0097           [](const Eigen::RowVectorXf &x) -> Eigen::RowVectorXf { return 2.0f * x; });
0098     m.def("double_complex",
0099           [](const Eigen::VectorXcf &x) -> Eigen::VectorXcf { return 2.0f * x; });
0100     m.def("double_threec", [](py::EigenDRef<Eigen::Vector3f> x) { x *= 2; });
0101     m.def("double_threer", [](py::EigenDRef<Eigen::RowVector3f> x) { x *= 2; });
0102     m.def("double_mat_cm", [](const Eigen::MatrixXf &x) -> Eigen::MatrixXf { return 2.0f * x; });
0103     m.def("double_mat_rm", [](const DenseMatrixR &x) -> DenseMatrixR { return 2.0f * x; });
0104
0105     // test_eigen_ref_to_python
0106     // Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended
0107     m.def("cholesky1",
0108           [](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
0109     m.def("cholesky2", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd {
0110         return x.llt().matrixL();
0111     });
0112     m.def("cholesky3",
0113           [](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
0114     m.def("cholesky4", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd {
0115         return x.llt().matrixL();
0116     });
0117
0118     // test_eigen_ref_mutators
0119     // Mutators: these add some value to the given element using Eigen, but Eigen should be mapping
0120     // into the numpy array data and so the result should show up there.  There are three versions:
0121     // one that works on a contiguous-row matrix (numpy's default), one for a contiguous-column
0122     // matrix, and one for any matrix.
0123     auto add_rm = [](Eigen::Ref<MatrixXdR> x, int r, int c, double v) { x(r, c) += v; };
0124     auto add_cm = [](Eigen::Ref<Eigen::MatrixXd> x, int r, int c, double v) { x(r, c) += v; };
0125
0126     // Mutators (Eigen maps into numpy variables):
0127     m.def("add_rm", add_rm); // Only takes row-contiguous
0128     m.def("add_cm", add_cm); // Only takes column-contiguous
0129     // Overloaded versions that will accept either row or column contiguous:
0130     m.def("add1", add_rm);
0131     m.def("add1", add_cm);
0132     m.def("add2", add_cm);
0133     m.def("add2", add_rm);
0134     // This one accepts a matrix of any stride:
0135     m.def("add_any",
0136           [](py::EigenDRef<Eigen::MatrixXd> x, int r, int c, double v) { x(r, c) += v; });
0137
0138     // Return mutable references (numpy maps into eigen variables)
0139     m.def("get_cm_ref", []() { return Eigen::Ref<Eigen::MatrixXd>(get_cm()); });
0140     m.def("get_rm_ref", []() { return Eigen::Ref<MatrixXdR>(get_rm()); });
0141     // The same references, but non-mutable (numpy maps into eigen variables, but is !writeable)
0142     m.def("get_cm_const_ref", []() { return Eigen::Ref<const Eigen::MatrixXd>(get_cm()); });
0143     m.def("get_rm_const_ref", []() { return Eigen::Ref<const MatrixXdR>(get_rm()); });
0144
0145     m.def("reset_refs", reset_refs); // Restores get_{cm,rm}_ref to original values
0146
0147     // Increments and returns ref to (same) matrix
0148     m.def(
0149         "incr_matrix",
0150         [](Eigen::Ref<Eigen::MatrixXd> m, double v) {
0151             m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v);
0152             return m;
0153         },
0154         py::return_value_policy::reference);
0155
0156     // Same, but accepts a matrix of any strides
0157     m.def(
0158         "incr_matrix_any",
0159         [](py::EigenDRef<Eigen::MatrixXd> m, double v) {
0160             m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v);
0161             return m;
0162         },
0163         py::return_value_policy::reference);
0164
0165     // Returns an eigen slice of even rows
0166     m.def(
0167         "even_rows",
0168         [](py::EigenDRef<Eigen::MatrixXd> m) {
0169             return py::EigenDMap<Eigen::MatrixXd>(
0170                 m.data(),
0171                 (m.rows() + 1) / 2,
0172                 m.cols(),
0173                 py::EigenDStride(m.outerStride(), 2 * m.innerStride()));
0174         },
0175         py::return_value_policy::reference);
0176
0177     // Returns an eigen slice of even columns
0178     m.def(
0179         "even_cols",
0180         [](py::EigenDRef<Eigen::MatrixXd> m) {
0181             return py::EigenDMap<Eigen::MatrixXd>(
0182                 m.data(),
0183                 m.rows(),
0184                 (m.cols() + 1) / 2,
0185                 py::EigenDStride(2 * m.outerStride(), m.innerStride()));
0186         },
0187         py::return_value_policy::reference);
0188
0189     // Returns diagonals: a vector-like object with an inner stride != 1
0190     m.def("diagonal", [](const Eigen::Ref<const Eigen::MatrixXd> &x) { return x.diagonal(); });
0191     m.def("diagonal_1",
0192           [](const Eigen::Ref<const Eigen::MatrixXd> &x) { return x.diagonal<1>(); });
0193     m.def("diagonal_n",
0194           [](const Eigen::Ref<const Eigen::MatrixXd> &x, int index) { return x.diagonal(index); });
0195
0196     // Return a block of a matrix (gives non-standard strides)
0197     m.def("block",
0198           [m](const py::object &x_obj,
0199               int start_row,
0200               int start_col,
0201               int block_rows,
0202               int block_cols) {
0203               return m.attr("_block")(x_obj, x_obj, start_row, start_col, block_rows, block_cols);
0204           });
0205
0206     m.def(
0207         "_block",
0208         [](const py::object &x_obj,
0209            const Eigen::Ref<const Eigen::MatrixXd> &x,
0210            int start_row,
0211            int start_col,
0212            int block_rows,
0213            int block_cols) {
0214             // See PR #4217 for background. This test is a bit over the top, but might be useful
0215             // as a concrete example to point to when explaining the dangling reference trap.
0216             auto i0 = py::make_tuple(0, 0);
0217             auto x0_orig = x_obj[*i0].cast<double>();
0218             if (x(0, 0) != x0_orig) {
0219                 throw std::runtime_error(
0220                     "Something in the type_caster for Eigen::Ref is terribly wrong.");
0221             }
0222             double x0_mod = x0_orig + 1;
0223             x_obj[*i0] = x0_mod;
0224             auto copy_detected = (x(0, 0) != x0_mod);
0225             x_obj[*i0] = x0_orig;
0226             if (copy_detected) {
0227                 throw std::runtime_error("type_caster for Eigen::Ref made a copy.");
0228             }
0229             return x.block(start_row, start_col, block_rows, block_cols);
0230         },
0231         py::keep_alive<0, 1>());
0232
0233     // test_eigen_return_references, test_eigen_keepalive
0234     // return value referencing/copying tests:
0235     class ReturnTester {
0236         Eigen::MatrixXd mat = create();
0237
0238     public:
0239         ReturnTester() { print_created(this); }
0240         ~ReturnTester() { print_destroyed(this); }
0241         static Eigen::MatrixXd create() { return Eigen::MatrixXd::Ones(10, 10); }
0242         // NOLINTNEXTLINE(readability-const-return-type)
0243         static const Eigen::MatrixXd createConst() { return Eigen::MatrixXd::Ones(10, 10); }
0244         Eigen::MatrixXd &get() { return mat; }
0245         Eigen::MatrixXd *getPtr() { return &mat; }
0246         const Eigen::MatrixXd &view() { return mat; }
0247         const Eigen::MatrixXd *viewPtr() { return &mat; }
0248         Eigen::Ref<Eigen::MatrixXd> ref() { return mat; }
0249         Eigen::Ref<const Eigen::MatrixXd> refConst() { return mat; }
0250         Eigen::Block<Eigen::MatrixXd> block(int r, int c, int nrow, int ncol) {
0251             return mat.block(r, c, nrow, ncol);
0252         }
0253         Eigen::Block<const Eigen::MatrixXd> blockConst(int r, int c, int nrow, int ncol) const {
0254             return mat.block(r, c, nrow, ncol);
0255         }
0256         py::EigenDMap<Eigen::Matrix2d> corners() {
0257             return py::EigenDMap<Eigen::Matrix2d>(
0258                 mat.data(),
0259                 py::EigenDStride(mat.outerStride() * (mat.outerSize() - 1),
0260                                  mat.innerStride() * (mat.innerSize() - 1)));
0261         }
0262         py::EigenDMap<const Eigen::Matrix2d> cornersConst() const {
0263             return py::EigenDMap<const Eigen::Matrix2d>(
0264                 mat.data(),
0265                 py::EigenDStride(mat.outerStride() * (mat.outerSize() - 1),
0266                                  mat.innerStride() * (mat.innerSize() - 1)));
0267         }
0268     };
0269     using rvp = py::return_value_policy;
0270     py::class_<ReturnTester>(m, "ReturnTester")
0271         .def(py::init<>())
0272         .def_static("create", &ReturnTester::create)
0273         .def_static("create_const", &ReturnTester::createConst)
0274         .def("get", &ReturnTester::get, rvp::reference_internal)
0275         .def("get_ptr", &ReturnTester::getPtr, rvp::reference_internal)
0276         .def("view", &ReturnTester::view, rvp::reference_internal)
0277         .def("view_ptr", &ReturnTester::view, rvp::reference_internal)
0278         .def("copy_get", &ReturnTester::get)       // Default rvp: copy
0279         .def("copy_view", &ReturnTester::view)     //         "
0280         .def("ref", &ReturnTester::ref)            // Default for Ref is to reference
0281         .def("ref_const", &ReturnTester::refConst) // Likewise, but const
0282         .def("ref_safe", &ReturnTester::ref, rvp::reference_internal)
0283         .def("ref_const_safe", &ReturnTester::refConst, rvp::reference_internal)
0284         .def("copy_ref", &ReturnTester::ref, rvp::copy)
0285         .def("copy_ref_const", &ReturnTester::refConst, rvp::copy)
0286         .def("block", &ReturnTester::block)
0287         .def("block_safe", &ReturnTester::block, rvp::reference_internal)
0288         .def("block_const", &ReturnTester::blockConst, rvp::reference_internal)
0289         .def("copy_block", &ReturnTester::block, rvp::copy)
0290         .def("corners", &ReturnTester::corners, rvp::reference_internal)
0291         .def("corners_const", &ReturnTester::cornersConst, rvp::reference_internal);
0292
0293     // test_special_matrix_objects
0294     // Returns a DiagonalMatrix with diagonal (1,2,3,...)
0295     m.def("incr_diag", [](int k) {
0296         Eigen::DiagonalMatrix<int, Eigen::Dynamic> m(k);
0297         for (int i = 0; i < k; i++) {
0298             m.diagonal()[i] = i + 1;
0299         }
0300         return m;
0301     });
0302
0303     // Returns a SelfAdjointView referencing the lower triangle of m
0304     m.def("symmetric_lower",
0305           [](const Eigen::MatrixXi &m) { return m.selfadjointView<Eigen::Lower>(); });
0306     // Returns a SelfAdjointView referencing the lower triangle of m
0307     m.def("symmetric_upper",
0308           [](const Eigen::MatrixXi &m) { return m.selfadjointView<Eigen::Upper>(); });
0309
0310     // Test matrix for various functions below.
0311     Eigen::MatrixXf mat(5, 6);
0312     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,
0313         0, 8, 11;
0314
0315     // test_fixed, and various other tests
0316     m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); });
0317     // Our Eigen does a hack which respects constness through the numpy writeable flag.
0318     // Therefore, the const return actually affects this type despite being an rvalue.
0319     // NOLINTNEXTLINE(readability-const-return-type)
0320     m.def("fixed_r_const", [mat]() -> const FixedMatrixR { return FixedMatrixR(mat); });
0321     m.def("fixed_c", [mat]() -> FixedMatrixC { return FixedMatrixC(mat); });
0322     m.def("fixed_copy_r", [](const FixedMatrixR &m) -> FixedMatrixR { return m; });
0323     m.def("fixed_copy_c", [](const FixedMatrixC &m) -> FixedMatrixC { return m; });
0324     // test_mutator_descriptors
0325     m.def("fixed_mutator_r", [](const Eigen::Ref<FixedMatrixR> &) {});
0326     m.def("fixed_mutator_c", [](const Eigen::Ref<FixedMatrixC> &) {});
0327     m.def("fixed_mutator_a", [](const py::EigenDRef<FixedMatrixC> &) {});
0328     // test_dense
0329     m.def("dense_r", [mat]() -> DenseMatrixR { return DenseMatrixR(mat); });
0330     m.def("dense_c", [mat]() -> DenseMatrixC { return DenseMatrixC(mat); });
0331     m.def("dense_copy_r", [](const DenseMatrixR &m) -> DenseMatrixR { return m; });
0332     m.def("dense_copy_c", [](const DenseMatrixC &m) -> DenseMatrixC { return m; });
0333     // test_sparse, test_sparse_signature
0334     m.def("sparse_r", [mat]() -> SparseMatrixR {
0335         // NOLINTNEXTLINE(clang-analyzer-core.uninitialized.UndefReturn)
0336         return Eigen::SparseView<Eigen::MatrixXf>(mat);
0337     });
0338     m.def("sparse_c",
0339           [mat]() -> SparseMatrixC { return Eigen::SparseView<Eigen::MatrixXf>(mat); });
0340     m.def("sparse_copy_r", [](const SparseMatrixR &m) -> SparseMatrixR { return m; });
0341     m.def("sparse_copy_c", [](const SparseMatrixC &m) -> SparseMatrixC { return m; });
0342     // test_partially_fixed
0343     m.def("partial_copy_four_rm_r", [](const FourRowMatrixR &m) -> FourRowMatrixR { return m; });
0344     m.def("partial_copy_four_rm_c", [](const FourColMatrixR &m) -> FourColMatrixR { return m; });
0345     m.def("partial_copy_four_cm_r", [](const FourRowMatrixC &m) -> FourRowMatrixC { return m; });
0346     m.def("partial_copy_four_cm_c", [](const FourColMatrixC &m) -> FourColMatrixC { return m; });
0347
0348     // test_cpp_casting
0349     // Test that we can cast a numpy object to a Eigen::MatrixXd explicitly
0350     m.def("cpp_copy", [](py::handle m) { return m.cast<Eigen::MatrixXd>()(1, 0); });
0351     m.def("cpp_ref_c", [](py::handle m) { return m.cast<Eigen::Ref<Eigen::MatrixXd>>()(1, 0); });
0352     m.def("cpp_ref_r", [](py::handle m) { return m.cast<Eigen::Ref<MatrixXdR>>()(1, 0); });
0353     m.def("cpp_ref_any",
0354           [](py::handle m) { return m.cast<py::EigenDRef<Eigen::MatrixXd>>()(1, 0); });
0355
0356     // [workaround(intel)] ICC 20/21 breaks with py::arg().stuff, using py::arg{}.stuff works.
0357
0358     // test_nocopy_wrapper
0359     // Test that we can prevent copying into an argument that would normally copy: First a version
0360     // that would allow copying (if types or strides don't match) for comparison:
0361     m.def("get_elem", &get_elem);
0362     // Now this alternative that calls the tells pybind to fail rather than copy:
0363     m.def(
0364         "get_elem_nocopy",
0365         [](const Eigen::Ref<const Eigen::MatrixXd> &m) -> double { return get_elem(m); },
0366         py::arg{}.noconvert());
0367     // Also test a row-major-only no-copy const ref:
0368     m.def(
0369         "get_elem_rm_nocopy",
0370         [](Eigen::Ref<const Eigen::Matrix<long, -1, -1, Eigen::RowMajor>> &m) -> long {
0371             return m(2, 1);
0372         },
0373         py::arg{}.noconvert());
0374
0375     // test_issue738, test_zero_length
0376     // Issue #738: 1×N or N×1 2D matrices were neither accepted nor properly copied with an
0377     // incompatible stride value on the length-1 dimension--but that should be allowed (without
0378     // requiring a copy!) because the stride value can be safely ignored on a size-1 dimension.
0379     // Similarly, 0×N or N×0 matrices were not accepted--again, these should be allowed since
0380     // they contain no data. This particularly affects numpy ≥ 1.23, which sets the strides to
0381     // 0 if any dimension size is 0.
0382     m.def("iss738_f1",
0383           &adjust_matrix<const Eigen::Ref<const Eigen::MatrixXd> &>,
0384           py::arg{}.noconvert());
0385     m.def("iss738_f2",
0386           &adjust_matrix<const Eigen::Ref<const Eigen::Matrix<double, -1, -1, Eigen::RowMajor>> &>,
0387           py::arg{}.noconvert());
0388
0389     // test_issue1105
0390     // Issue #1105: when converting from a numpy two-dimensional (Nx1) or (1xN) value into a dense
0391     // eigen Vector or RowVector, the argument would fail to load because the numpy copy would
0392     // fail: numpy won't broadcast a Nx1 into a 1-dimensional vector.
0393     m.def("iss1105_col", [](const Eigen::VectorXd &) { return true; });
0394     m.def("iss1105_row", [](const Eigen::RowVectorXd &) { return true; });
0395
0396     // test_named_arguments
0397     // Make sure named arguments are working properly:
0398     m.def(
0399         "matrix_multiply",
0400         [](const py::EigenDRef<const Eigen::MatrixXd> &A,
0401            const py::EigenDRef<const Eigen::MatrixXd> &B) -> Eigen::MatrixXd {
0402             if (A.cols() != B.rows()) {
0403                 throw std::domain_error("Nonconformable matrices!");
0404             }
0405             return A * B;
0406         },
0407         py::arg("A"),
0408         py::arg("B"));
0409
0410     // test_custom_operator_new
0411     py::class_<CustomOperatorNew>(m, "CustomOperatorNew")
0412         .def(py::init<>())
0413         .def_readonly("a", &CustomOperatorNew::a)
0414         .def_readonly("b", &CustomOperatorNew::b);
0415
0416     // test_eigen_ref_life_support
0417     // In case of a failure (the caster's temp array does not live long enough), creating
0418     // a new array (np.ones(10)) increases the chances that the temp array will be garbage
0419     // collected and/or that its memory will be overridden with different values.
0420     m.def("get_elem_direct", [](const Eigen::Ref<const Eigen::VectorXd> &v) {
0421         py::module_::import("numpy").attr("ones")(10);
0422         return v(5);
0423     });
0424     m.def("get_elem_indirect", [](std::vector<Eigen::Ref<const Eigen::VectorXd>> v) {
0425         py::module_::import("numpy").attr("ones")(10);
0426         return v[0](5);
0427     });
0428 }