Zengtudor
/
pybind11
mirror of https://github.com/pybind/pybind11.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
pybind11/tests/test_eigen_matrix.cpp

444 lines
20 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/*
tests/eigen.cpp -- automatic conversion of Eigen types
Copyright (c) 2016 Wenzel Jakob <wenzel.jakob@epfl.ch>
All rights reserved. Use of this source code is governed by a
BSD-style license that can be found in the LICENSE file.
*/
#include <pybind11/eigen/matrix.h>
#include <pybind11/stl.h>
#include "constructor_stats.h"
#include "pybind11_tests.h"
PYBIND11_WARNING_DISABLE_MSVC(4996)
#include <Eigen/Cholesky>
using MatrixXdR = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
// 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 <typename M>
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;
}
}
}
// 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(const Eigen::Ref<const Eigen::MatrixXd> &m) { return m(2, 1); };
// Returns a matrix with 10*r + 100*c added to each matrix element (to help test that the matrix
// reference is referencing rows/columns correctly).
template <typename MatrixArgType>
Eigen::MatrixXd adjust_matrix(MatrixArgType m) {
Eigen::MatrixXd ret(m);
for (int c = 0; c < m.cols(); c++) {
for (int r = 0; r < m.rows(); r++) {
ret(r, c) += 10 * r + 100 * c; // NOLINT(clang-analyzer-core.uninitialized.Assign)
}
}
return ret;
}
struct CustomOperatorNew {
CustomOperatorNew() = default;
Eigen::Matrix4d a = Eigen::Matrix4d::Zero();
Eigen::Matrix4d b = Eigen::Matrix4d::Identity();
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
TEST_SUBMODULE(eigen_matrix, m) {
using FixedMatrixR = Eigen::Matrix<float, 5, 6, Eigen::RowMajor>;
using FixedMatrixC = Eigen::Matrix<float, 5, 6>;
using DenseMatrixR = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using DenseMatrixC = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic>;
using FourRowMatrixC = Eigen::Matrix<float, 4, Eigen::Dynamic>;
using FourColMatrixC = Eigen::Matrix<float, Eigen::Dynamic, 4>;
using FourRowMatrixR = Eigen::Matrix<float, 4, Eigen::Dynamic>;
using FourColMatrixR = Eigen::Matrix<float, Eigen::Dynamic, 4>;
using SparseMatrixR = Eigen::SparseMatrix<float, Eigen::RowMajor>;
using SparseMatrixC = Eigen::SparseMatrix<float>;
// various tests
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_complex",
[](const Eigen::VectorXcf &x) -> Eigen::VectorXcf { return 2.0f * x; });
m.def("double_threec", [](py::EigenDRef<Eigen::Vector3f> x) { x *= 2; });
m.def("double_threer", [](py::EigenDRef<Eigen::RowVector3f> x) { x *= 2; });
m.def("double_mat_cm", [](const Eigen::MatrixXf &x) -> Eigen::MatrixXf { return 2.0f * x; });
m.def("double_mat_rm", [](const DenseMatrixR &x) -> DenseMatrixR { return 2.0f * x; });
// test_eigen_ref_to_python
// Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended
m.def("cholesky1",
[](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
m.def("cholesky2", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd {
return x.llt().matrixL();
});
m.def("cholesky3",
[](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
m.def("cholesky4", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd {
return x.llt().matrixL();
});
// test_eigen_ref_mutators
// 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<MatrixXdR> x, int r, int c, double v) { x(r, c) += v; };
auto add_cm = [](Eigen::Ref<Eigen::MatrixXd> 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<Eigen::MatrixXd> x, int r, int c, double v) { x(r, c) += v; });
// Return mutable references (numpy maps into eigen variables)
m.def("get_cm_ref", []() { return Eigen::Ref<Eigen::MatrixXd>(get_cm()); });
m.def("get_rm_ref", []() { return Eigen::Ref<MatrixXdR>(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<const Eigen::MatrixXd>(get_cm()); });
m.def("get_rm_const_ref", []() { return Eigen::Ref<const MatrixXdR>(get_rm()); });
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<Eigen::MatrixXd> 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<Eigen::MatrixXd> 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<Eigen::MatrixXd> m) {
return py::EigenDMap<Eigen::MatrixXd>(
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<Eigen::MatrixXd> m) {
return py::EigenDMap<Eigen::MatrixXd>(
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<const Eigen::MatrixXd> &x) { return x.diagonal(); });
m.def("diagonal_1",
[](const Eigen::Ref<const Eigen::MatrixXd> &x) { return x.diagonal<1>(); });
m.def("diagonal_n",
[](const Eigen::Ref<const Eigen::MatrixXd> &x, int index) { return x.diagonal(index); });
// Return a block of a matrix (gives non-standard strides)
m.def("block",
[m](const py::object &x_obj,
int start_row,
int start_col,
int block_rows,
int block_cols) {
return m.attr("_block")(x_obj, x_obj, start_row, start_col, block_rows, block_cols);
});
m.def(
"_block",
[](const py::object &x_obj,
const Eigen::Ref<const Eigen::MatrixXd> &x,
int start_row,
int start_col,
int block_rows,
int block_cols) {
// See PR #4217 for background. This test is a bit over the top, but might be useful
// as a concrete example to point to when explaining the dangling reference trap.
auto i0 = py::make_tuple(0, 0);
auto x0_orig = x_obj[*i0].cast<double>();
if (x(0, 0) != x0_orig) {
throw std::runtime_error(
"Something in the type_caster for Eigen::Ref is terribly wrong.");
}
double x0_mod = x0_orig + 1;
x_obj[*i0] = x0_mod;
auto copy_detected = (x(0, 0) != x0_mod);
x_obj[*i0] = x0_orig;
if (copy_detected) {
throw std::runtime_error("type_caster for Eigen::Ref made a copy.");
}
return x.block(start_row, start_col, block_rows, block_cols);
},
py::keep_alive<0, 1>());
// test_eigen_return_references, test_eigen_keepalive
// 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); }
// NOLINTNEXTLINE(readability-const-return-type)
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<Eigen::MatrixXd> ref() { return mat; }
Eigen::Ref<const Eigen::MatrixXd> refConst() { return mat; }
Eigen::Block<Eigen::MatrixXd> block(int r, int c, int nrow, int ncol) {
return mat.block(r, c, nrow, ncol);
}
Eigen::Block<const Eigen::MatrixXd> blockConst(int r, int c, int nrow, int ncol) const {
return mat.block(r, c, nrow, ncol);
}
py::EigenDMap<Eigen::Matrix2d> corners() {
return py::EigenDMap<Eigen::Matrix2d>(
mat.data(),
py::EigenDStride(mat.outerStride() * (mat.outerSize() - 1),
mat.innerStride() * (mat.innerSize() - 1)));
}
py::EigenDMap<const Eigen::Matrix2d> cornersConst() const {
return py::EigenDMap<const Eigen::Matrix2d>(
mat.data(),
py::EigenDStride(mat.outerStride() * (mat.outerSize() - 1),
mat.innerStride() * (mat.innerSize() - 1)));
}
};
using rvp = py::return_value_policy;
py::class_<ReturnTester>(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);
// test_special_matrix_objects
// Returns a DiagonalMatrix with diagonal (1,2,3,...)
m.def("incr_diag", [](int k) {
Eigen::DiagonalMatrix<int, Eigen::Dynamic> 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<Eigen::Lower>(); });
// Returns a SelfAdjointView referencing the lower triangle of m
m.def("symmetric_upper",
[](const Eigen::MatrixXi &m) { return m.selfadjointView<Eigen::Upper>(); });
// 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;
// test_fixed, and various other tests
m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); });
// Our Eigen does a hack which respects constness through the numpy writeable flag.
// Therefore, the const return actually affects this type despite being an rvalue.
// NOLINTNEXTLINE(readability-const-return-type)
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; });
// test_mutator_descriptors
m.def("fixed_mutator_r", [](const Eigen::Ref<FixedMatrixR> &) {});
m.def("fixed_mutator_c", [](const Eigen::Ref<FixedMatrixC> &) {});
m.def("fixed_mutator_a", [](const py::EigenDRef<FixedMatrixC> &) {});
// test_dense
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; });
// test_defaults
bool have_numpy = true;
try {
py::module_::import("numpy");
} catch (const py::error_already_set &) {
have_numpy = false;
}
if (have_numpy) {
py::module_::import("numpy");
Eigen::Matrix<double, 3, 3> defaultMatrix = Eigen::Matrix3d::Identity();
m.def("defaults_mat", [](const Eigen::Matrix3d &) {}, py::arg("mat") = defaultMatrix);
Eigen::VectorXd defaultVector = Eigen::VectorXd::Ones(32);
m.def("defaults_vec", [](const Eigen::VectorXd &) {}, py::arg("vec") = defaultMatrix);
}
// test_sparse, test_sparse_signature
m.def("sparse_r", [mat]() -> SparseMatrixR {
// NOLINTNEXTLINE(clang-analyzer-core.uninitialized.UndefReturn)
return Eigen::SparseView<Eigen::MatrixXf>(mat);
});
m.def("sparse_c",
[mat]() -> SparseMatrixC { return Eigen::SparseView<Eigen::MatrixXf>(mat); });
m.def("sparse_copy_r", [](const SparseMatrixR &m) -> SparseMatrixR { return m; });
m.def("sparse_copy_c", [](const SparseMatrixC &m) -> SparseMatrixC { return m; });
// test_partially_fixed
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; });
// test_cpp_casting
// Test that we can cast a numpy object to a Eigen::MatrixXd explicitly
m.def("cpp_copy", [](py::handle m) { return m.cast<Eigen::MatrixXd>()(1, 0); });
m.def("cpp_ref_c", [](py::handle m) { return m.cast<Eigen::Ref<Eigen::MatrixXd>>()(1, 0); });
m.def("cpp_ref_r", [](py::handle m) { return m.cast<Eigen::Ref<MatrixXdR>>()(1, 0); });
m.def("cpp_ref_any",
[](py::handle m) { return m.cast<py::EigenDRef<Eigen::MatrixXd>>()(1, 0); });
// [workaround(intel)] ICC 20/21 breaks with py::arg().stuff, using py::arg{}.stuff works.
// test_nocopy_wrapper
// 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",
[](const Eigen::Ref<const Eigen::MatrixXd> &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<const Eigen::Matrix<long, -1, -1, Eigen::RowMajor>> &m) -> long {
return m(2, 1);
},
py::arg{}.noconvert());
// test_issue738, test_zero_length
// Issue #738: 1×N or N×1 2D matrices were neither accepted nor properly copied with an
// incompatible stride value on the length-1 dimension--but that should be allowed (without
// requiring a copy!) because the stride value can be safely ignored on a size-1 dimension.
// Similarly, 0×N or N×0 matrices were not accepted--again, these should be allowed since
// they contain no data. This particularly affects numpy ≥ 1.23, which sets the strides to
// 0 if any dimension size is 0.
m.def("iss738_f1",
&adjust_matrix<const Eigen::Ref<const Eigen::MatrixXd> &>,
py::arg{}.noconvert());
m.def("iss738_f2",
&adjust_matrix<const Eigen::Ref<const Eigen::Matrix<double, -1, -1, Eigen::RowMajor>> &>,
py::arg{}.noconvert());
// test_issue1105
// Issue #1105: when converting from a numpy two-dimensional (Nx1) or (1xN) value into a dense
// eigen Vector or RowVector, the argument would fail to load because the numpy copy would
// fail: numpy won't broadcast a Nx1 into a 1-dimensional vector.
m.def("iss1105_col", [](const Eigen::VectorXd &) { return true; });
m.def("iss1105_row", [](const Eigen::RowVectorXd &) { return true; });
// test_named_arguments
// Make sure named arguments are working properly:
m.def(
"matrix_multiply",
[](const py::EigenDRef<const Eigen::MatrixXd> &A,
const py::EigenDRef<const Eigen::MatrixXd> &B) -> Eigen::MatrixXd {
if (A.cols() != B.rows()) {
throw std::domain_error("Nonconformable matrices!");
}
return A * B;
},
py::arg("A"),
py::arg("B"));
// test_custom_operator_new
py::class_<CustomOperatorNew>(m, "CustomOperatorNew")
.def(py::init<>())
.def_readonly("a", &CustomOperatorNew::a)
.def_readonly("b", &CustomOperatorNew::b);
// test_eigen_ref_life_support
// In case of a failure (the caster's temp array does not live long enough), creating
// a new array (np.ones(10)) increases the chances that the temp array will be garbage
// collected and/or that its memory will be overridden with different values.
m.def("get_elem_direct", [](const Eigen::Ref<const Eigen::VectorXd> &v) {
py::module_::import("numpy").attr("ones")(10);
return v(5);
});
m.def("get_elem_indirect", [](std::vector<Eigen::Ref<const Eigen::VectorXd>> v) {
py::module_::import("numpy").attr("ones")(10);
return v[0](5);
});
}
Reference in New Issue View Git Blame Copy Permalink