test_eigen.cpp 16.7 KB
Newer Older
1
/*
Dean Moldovan's avatar
Dean Moldovan committed
2
    tests/eigen.cpp -- automatic conversion of Eigen types
3
4
5
6
7
8
9

    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.
*/

Dean Moldovan's avatar
Dean Moldovan committed
10
#include "pybind11_tests.h"
11
#include "constructor_stats.h"
12
#include <pybind11/eigen.h>
13
#include <pybind11/stl.h>
14
15
16
17
18

#if defined(_MSC_VER)
#  pragma warning(disable: 4996) // C4996: std::unary_negation is deprecated
#endif

19
#include <Eigen/Cholesky>
20

21
using MatrixXdR = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
22
23


24

25
26
27
28
29
30
// 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;
}
31

32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
// 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)
57
double get_elem(const Eigen::Ref<const Eigen::MatrixXd> &m) { return m(2, 1); };
58
59
60
61
62

// 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);
63
64
65
    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)
66
67
68
    return ret;
}

69
70
71
72
73
74
75
76
77
struct CustomOperatorNew {
    CustomOperatorNew() = default;

    Eigen::Matrix4d a = Eigen::Matrix4d::Zero();
    Eigen::Matrix4d b = Eigen::Matrix4d::Identity();

    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};

78
79
80
81
82
83
84
85
86
87
88
TEST_SUBMODULE(eigen, 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>;
89

90
    // various tests
91
92
    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; });
93
    m.def("double_complex", [](const Eigen::VectorXcf &x) -> Eigen::VectorXcf { return 2.0f * x; });
94
95
    m.def("double_threec", [](py::EigenDRef<Eigen::Vector3f> x) { x *= 2; });
    m.def("double_threer", [](py::EigenDRef<Eigen::RowVector3f> x) { x *= 2; });
96
97
    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; });
98

99
    // test_eigen_ref_to_python
100
    // Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended
101
102
    m.def("cholesky1",
          [](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
103
104
    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(); });
105
106
107
    m.def("cholesky4", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd {
        return x.llt().matrixL();
    });
108

109
    // test_eigen_ref_mutators
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    // 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; });

luz.paz's avatar
Typo  
luz.paz committed
128
    // Return mutable references (numpy maps into eigen variables)
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    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);
162

163
164
165
166
167
168
169
170
171
172
    // 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", [](const Eigen::Ref<const Eigen::MatrixXd> &x, int start_row, int start_col, int block_rows, int block_cols) {
        return x.block(start_row, start_col, block_rows, block_cols);
    });

173
    // test_eigen_return_references, test_eigen_keepalive
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    // 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<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)
        ;

220
    // test_special_matrix_objects
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    // 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>();
    });

237
238
239
240
241
242
243
    // 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;
244

245
    // test_fixed, and various other tests
246
247
248
249
250
    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; });
251
    // test_mutator_descriptors
252
253
254
    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> &) {});
255
    // test_dense
256
257
258
259
    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; });
260
    // test_sparse, test_sparse_signature
261
262
263
264
    m.def("sparse_r", [mat]() -> SparseMatrixR {
        // NOLINTNEXTLINE(clang-analyzer-core.uninitialized.UndefReturn)
        return Eigen::SparseView<Eigen::MatrixXf>(mat);
    });
265
266
267
    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; });
268
    // test_partially_fixed
269
270
271
272
    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; });
273

274
    // test_cpp_casting
275
276
277
278
279
    // 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); });
280

281
    // [workaround(intel)] ICC 20/21 breaks with py::arg().stuff, using py::arg{}.stuff works.
282

283
    // test_nocopy_wrapper
284
285
286
287
    // 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:
288
289
290
291
    m.def(
        "get_elem_nocopy",
        [](const Eigen::Ref<const Eigen::MatrixXd> &m) -> double { return get_elem(m); },
        py::arg{}.noconvert());
292
293
    // 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); },
294
            py::arg{}.noconvert());
295

296
    // test_issue738
297
298
299
    // Issue #738: 1xN or Nx1 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.
300
301
    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());
302

303
304
    // test_issue1105
    // Issue #1105: when converting from a numpy two-dimensional (Nx1) or (1xN) value into a dense
305
    // eigen Vector or RowVector, the argument would fail to load because the numpy copy would
306
307
308
    // 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; });
309

310
    // test_named_arguments
311
    // Make sure named arguments are working properly:
312
313
314
315
316
317
318
319
320
321
    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"));
322

323
    // test_custom_operator_new
324
325
326
327
    py::class_<CustomOperatorNew>(m, "CustomOperatorNew")
        .def(py::init<>())
        .def_readonly("a", &CustomOperatorNew::a)
        .def_readonly("b", &CustomOperatorNew::b);
328
329
330
331
332

    // 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.
333
    m.def("get_elem_direct", [](const Eigen::Ref<const Eigen::VectorXd> &v) {
334
        py::module_::import("numpy").attr("ones")(10);
335
336
337
        return v(5);
    });
    m.def("get_elem_indirect", [](std::vector<Eigen::Ref<const Eigen::VectorXd>> v) {
338
        py::module_::import("numpy").attr("ones")(10);
339
340
        return v[0](5);
    });
341
}