matrix.h 32 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723
  1. /*
  2. pybind11/eigen/matrix.h: Transparent conversion for dense and sparse Eigen matrices
  3. Copyright (c) 2016 Wenzel Jakob <wenzel.jakob@epfl.ch>
  4. All rights reserved. Use of this source code is governed by a
  5. BSD-style license that can be found in the LICENSE file.
  6. */
  7. #pragma once
  8. #include <pybind11/numpy.h>
  9. #include "common.h"
  10. /* HINT: To suppress warnings originating from the Eigen headers, use -isystem.
  11. See also:
  12. https://stackoverflow.com/questions/2579576/i-dir-vs-isystem-dir
  13. https://stackoverflow.com/questions/1741816/isystem-for-ms-visual-studio-c-compiler
  14. */
  15. PYBIND11_WARNING_PUSH
  16. PYBIND11_WARNING_DISABLE_MSVC(5054) // https://github.com/pybind/pybind11/pull/3741
  17. // C5054: operator '&': deprecated between enumerations of different types
  18. #if defined(__MINGW32__)
  19. PYBIND11_WARNING_DISABLE_GCC("-Wmaybe-uninitialized")
  20. #endif
  21. #include <Eigen/Core>
  22. #include <Eigen/SparseCore>
  23. PYBIND11_WARNING_POP
  24. // Eigen prior to 3.2.7 doesn't have proper move constructors--but worse, some classes get implicit
  25. // move constructors that break things. We could detect this an explicitly copy, but an extra copy
  26. // of matrices seems highly undesirable.
  27. static_assert(EIGEN_VERSION_AT_LEAST(3, 2, 7),
  28. "Eigen matrix support in pybind11 requires Eigen >= 3.2.7");
  29. PYBIND11_NAMESPACE_BEGIN(PYBIND11_NAMESPACE)
  30. PYBIND11_WARNING_DISABLE_MSVC(4127)
  31. // Provide a convenience alias for easier pass-by-ref usage with fully dynamic strides:
  32. using EigenDStride = Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>;
  33. template <typename MatrixType>
  34. using EigenDRef = Eigen::Ref<MatrixType, 0, EigenDStride>;
  35. template <typename MatrixType>
  36. using EigenDMap = Eigen::Map<MatrixType, 0, EigenDStride>;
  37. PYBIND11_NAMESPACE_BEGIN(detail)
  38. #if EIGEN_VERSION_AT_LEAST(3, 3, 0)
  39. using EigenIndex = Eigen::Index;
  40. template <typename Scalar, int Flags, typename StorageIndex>
  41. using EigenMapSparseMatrix = Eigen::Map<Eigen::SparseMatrix<Scalar, Flags, StorageIndex>>;
  42. #else
  43. using EigenIndex = EIGEN_DEFAULT_DENSE_INDEX_TYPE;
  44. template <typename Scalar, int Flags, typename StorageIndex>
  45. using EigenMapSparseMatrix = Eigen::MappedSparseMatrix<Scalar, Flags, StorageIndex>;
  46. #endif
  47. // Matches Eigen::Map, Eigen::Ref, blocks, etc:
  48. template <typename T>
  49. using is_eigen_dense_map = all_of<is_template_base_of<Eigen::DenseBase, T>,
  50. std::is_base_of<Eigen::MapBase<T, Eigen::ReadOnlyAccessors>, T>>;
  51. template <typename T>
  52. using is_eigen_mutable_map = std::is_base_of<Eigen::MapBase<T, Eigen::WriteAccessors>, T>;
  53. template <typename T>
  54. using is_eigen_dense_plain
  55. = all_of<negation<is_eigen_dense_map<T>>, is_template_base_of<Eigen::PlainObjectBase, T>>;
  56. template <typename T>
  57. using is_eigen_sparse = is_template_base_of<Eigen::SparseMatrixBase, T>;
  58. // Test for objects inheriting from EigenBase<Derived> that aren't captured by the above. This
  59. // basically covers anything that can be assigned to a dense matrix but that don't have a typical
  60. // matrix data layout that can be copied from their .data(). For example, DiagonalMatrix and
  61. // SelfAdjointView fall into this category.
  62. template <typename T>
  63. using is_eigen_other
  64. = all_of<is_template_base_of<Eigen::EigenBase, T>,
  65. negation<any_of<is_eigen_dense_map<T>, is_eigen_dense_plain<T>, is_eigen_sparse<T>>>>;
  66. // Captures numpy/eigen conformability status (returned by EigenProps::conformable()):
  67. template <bool EigenRowMajor>
  68. struct EigenConformable {
  69. bool conformable = false;
  70. EigenIndex rows = 0, cols = 0;
  71. EigenDStride stride{0, 0}; // Only valid if negativestrides is false!
  72. bool negativestrides = false; // If true, do not use stride!
  73. // NOLINTNEXTLINE(google-explicit-constructor)
  74. EigenConformable(bool fits = false) : conformable{fits} {}
  75. // Matrix type:
  76. EigenConformable(EigenIndex r, EigenIndex c, EigenIndex rstride, EigenIndex cstride)
  77. : conformable{true}, rows{r}, cols{c},
  78. // TODO: when Eigen bug #747 is fixed, remove the tests for non-negativity.
  79. // http://eigen.tuxfamily.org/bz/show_bug.cgi?id=747
  80. stride{EigenRowMajor ? (rstride > 0 ? rstride : 0)
  81. : (cstride > 0 ? cstride : 0) /* outer stride */,
  82. EigenRowMajor ? (cstride > 0 ? cstride : 0)
  83. : (rstride > 0 ? rstride : 0) /* inner stride */},
  84. negativestrides{rstride < 0 || cstride < 0} {}
  85. // Vector type:
  86. EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride)
  87. : EigenConformable(r, c, r == 1 ? c * stride : stride, c == 1 ? r : r * stride) {}
  88. template <typename props>
  89. bool stride_compatible() const {
  90. // To have compatible strides, we need (on both dimensions) one of fully dynamic strides,
  91. // matching strides, or a dimension size of 1 (in which case the stride value is
  92. // irrelevant). Alternatively, if any dimension size is 0, the strides are not relevant
  93. // (and numpy ≥ 1.23 sets the strides to 0 in that case, so we need to check explicitly).
  94. if (negativestrides) {
  95. return false;
  96. }
  97. if (rows == 0 || cols == 0) {
  98. return true;
  99. }
  100. return (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner()
  101. || (EigenRowMajor ? cols : rows) == 1)
  102. && (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer()
  103. || (EigenRowMajor ? rows : cols) == 1);
  104. }
  105. // NOLINTNEXTLINE(google-explicit-constructor)
  106. operator bool() const { return conformable; }
  107. };
  108. template <typename Type>
  109. struct eigen_extract_stride {
  110. using type = Type;
  111. };
  112. template <typename PlainObjectType, int MapOptions, typename StrideType>
  113. struct eigen_extract_stride<Eigen::Map<PlainObjectType, MapOptions, StrideType>> {
  114. using type = StrideType;
  115. };
  116. template <typename PlainObjectType, int Options, typename StrideType>
  117. struct eigen_extract_stride<Eigen::Ref<PlainObjectType, Options, StrideType>> {
  118. using type = StrideType;
  119. };
  120. // Helper struct for extracting information from an Eigen type
  121. template <typename Type_>
  122. struct EigenProps {
  123. using Type = Type_;
  124. using Scalar = typename Type::Scalar;
  125. using StrideType = typename eigen_extract_stride<Type>::type;
  126. static constexpr EigenIndex rows = Type::RowsAtCompileTime, cols = Type::ColsAtCompileTime,
  127. size = Type::SizeAtCompileTime;
  128. static constexpr bool row_major = Type::IsRowMajor,
  129. vector
  130. = Type::IsVectorAtCompileTime, // At least one dimension has fixed size 1
  131. fixed_rows = rows != Eigen::Dynamic, fixed_cols = cols != Eigen::Dynamic,
  132. fixed = size != Eigen::Dynamic, // Fully-fixed size
  133. dynamic = !fixed_rows && !fixed_cols; // Fully-dynamic size
  134. template <EigenIndex i, EigenIndex ifzero>
  135. using if_zero = std::integral_constant<EigenIndex, i == 0 ? ifzero : i>;
  136. static constexpr EigenIndex inner_stride
  137. = if_zero<StrideType::InnerStrideAtCompileTime, 1>::value,
  138. outer_stride = if_zero < StrideType::OuterStrideAtCompileTime,
  139. vector ? size
  140. : row_major ? cols
  141. : rows > ::value;
  142. static constexpr bool dynamic_stride
  143. = inner_stride == Eigen::Dynamic && outer_stride == Eigen::Dynamic;
  144. static constexpr bool requires_row_major
  145. = !dynamic_stride && !vector && (row_major ? inner_stride : outer_stride) == 1;
  146. static constexpr bool requires_col_major
  147. = !dynamic_stride && !vector && (row_major ? outer_stride : inner_stride) == 1;
  148. // Takes an input array and determines whether we can make it fit into the Eigen type. If
  149. // the array is a vector, we attempt to fit it into either an Eigen 1xN or Nx1 vector
  150. // (preferring the latter if it will fit in either, i.e. for a fully dynamic matrix type).
  151. static EigenConformable<row_major> conformable(const array &a) {
  152. const auto dims = a.ndim();
  153. if (dims < 1 || dims > 2) {
  154. return false;
  155. }
  156. if (dims == 2) { // Matrix type: require exact match (or dynamic)
  157. EigenIndex np_rows = a.shape(0), np_cols = a.shape(1),
  158. np_rstride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar)),
  159. np_cstride = a.strides(1) / static_cast<ssize_t>(sizeof(Scalar));
  160. if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols)) {
  161. return false;
  162. }
  163. return {np_rows, np_cols, np_rstride, np_cstride};
  164. }
  165. // Otherwise we're storing an n-vector. Only one of the strides will be used, but
  166. // whichever is used, we want the (single) numpy stride value.
  167. const EigenIndex n = a.shape(0),
  168. stride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar));
  169. if (vector) { // Eigen type is a compile-time vector
  170. if (fixed && size != n) {
  171. return false; // Vector size mismatch
  172. }
  173. return {rows == 1 ? 1 : n, cols == 1 ? 1 : n, stride};
  174. }
  175. if (fixed) {
  176. // The type has a fixed size, but is not a vector: abort
  177. return false;
  178. }
  179. if (fixed_cols) {
  180. // Since this isn't a vector, cols must be != 1. We allow this only if it exactly
  181. // equals the number of elements (rows is Dynamic, and so 1 row is allowed).
  182. if (cols != n) {
  183. return false;
  184. }
  185. return {1, n, stride};
  186. } // Otherwise it's either fully dynamic, or column dynamic; both become a column vector
  187. if (fixed_rows && rows != n) {
  188. return false;
  189. }
  190. return {n, 1, stride};
  191. }
  192. static constexpr bool show_writeable
  193. = is_eigen_dense_map<Type>::value && is_eigen_mutable_map<Type>::value;
  194. static constexpr bool show_order = is_eigen_dense_map<Type>::value;
  195. static constexpr bool show_c_contiguous = show_order && requires_row_major;
  196. static constexpr bool show_f_contiguous
  197. = !show_c_contiguous && show_order && requires_col_major;
  198. static constexpr auto descriptor
  199. = const_name("typing.Annotated[")
  200. + io_name("numpy.typing.ArrayLike, ", "numpy.typing.NDArray[")
  201. + npy_format_descriptor<Scalar>::name + io_name("", "]") + const_name(", \"[")
  202. + const_name<fixed_rows>(const_name<(size_t) rows>(), const_name("m")) + const_name(", ")
  203. + const_name<fixed_cols>(const_name<(size_t) cols>(), const_name("n"))
  204. + const_name("]\"")
  205. // For a reference type (e.g. Ref<MatrixXd>) we have other constraints that might need to
  206. // be satisfied: writeable=True (for a mutable reference), and, depending on the map's
  207. // stride options, possibly f_contiguous or c_contiguous. We include them in the
  208. // descriptor output to provide some hint as to why a TypeError is occurring (otherwise
  209. // it can be confusing to see that a function accepts a
  210. // 'typing.Annotated[numpy.typing.NDArray[numpy.float64], "[3,2]"]' and an error message
  211. // that you *gave* a numpy.ndarray of the right type and dimensions.
  212. + const_name<show_writeable>(", \"flags.writeable\"", "")
  213. + const_name<show_c_contiguous>(", \"flags.c_contiguous\"", "")
  214. + const_name<show_f_contiguous>(", \"flags.f_contiguous\"", "") + const_name("]");
  215. };
  216. // Casts an Eigen type to numpy array. If given a base, the numpy array references the src data,
  217. // otherwise it'll make a copy. writeable lets you turn off the writeable flag for the array.
  218. template <typename props>
  219. handle
  220. eigen_array_cast(typename props::Type const &src, handle base = handle(), bool writeable = true) {
  221. constexpr ssize_t elem_size = sizeof(typename props::Scalar);
  222. array a;
  223. if (props::vector) {
  224. a = array({src.size()}, {elem_size * src.innerStride()}, src.data(), base);
  225. } else {
  226. a = array({src.rows(), src.cols()},
  227. {elem_size * src.rowStride(), elem_size * src.colStride()},
  228. src.data(),
  229. base);
  230. }
  231. if (!writeable) {
  232. array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
  233. }
  234. return a.release();
  235. }
  236. // Takes an lvalue ref to some Eigen type and a (python) base object, creating a numpy array that
  237. // reference the Eigen object's data with `base` as the python-registered base class (if omitted,
  238. // the base will be set to None, and lifetime management is up to the caller). The numpy array is
  239. // non-writeable if the given type is const.
  240. template <typename props, typename Type>
  241. handle eigen_ref_array(Type &src, handle parent = none()) {
  242. // none here is to get past array's should-we-copy detection, which currently always
  243. // copies when there is no base. Setting the base to None should be harmless.
  244. return eigen_array_cast<props>(src, parent, !std::is_const<Type>::value);
  245. }
  246. // Takes a pointer to some dense, plain Eigen type, builds a capsule around it, then returns a
  247. // numpy array that references the encapsulated data with a python-side reference to the capsule to
  248. // tie its destruction to that of any dependent python objects. Const-ness is determined by
  249. // whether or not the Type of the pointer given is const.
  250. template <typename props, typename Type, typename = enable_if_t<is_eigen_dense_plain<Type>::value>>
  251. handle eigen_encapsulate(Type *src) {
  252. capsule base(src, [](void *o) { delete static_cast<Type *>(o); });
  253. return eigen_ref_array<props>(*src, base);
  254. }
  255. // Type caster for regular, dense matrix types (e.g. MatrixXd), but not maps/refs/etc. of dense
  256. // types.
  257. template <typename Type>
  258. struct type_caster<Type, enable_if_t<is_eigen_dense_plain<Type>::value>> {
  259. using Scalar = typename Type::Scalar;
  260. static_assert(!std::is_pointer<Scalar>::value,
  261. PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED);
  262. using props = EigenProps<Type>;
  263. bool load(handle src, bool convert) {
  264. // If we're in no-convert mode, only load if given an array of the correct type
  265. if (!convert && !isinstance<array_t<Scalar>>(src)) {
  266. return false;
  267. }
  268. // Coerce into an array, but don't do type conversion yet; the copy below handles it.
  269. auto buf = array::ensure(src);
  270. if (!buf) {
  271. return false;
  272. }
  273. auto dims = buf.ndim();
  274. if (dims < 1 || dims > 2) {
  275. return false;
  276. }
  277. auto fits = props::conformable(buf);
  278. if (!fits) {
  279. return false;
  280. }
  281. PYBIND11_WARNING_PUSH
  282. PYBIND11_WARNING_DISABLE_GCC("-Wmaybe-uninitialized") // See PR #5516
  283. // Allocate the new type, then build a numpy reference into it
  284. value = Type(fits.rows, fits.cols);
  285. PYBIND11_WARNING_POP
  286. auto ref = reinterpret_steal<array>(eigen_ref_array<props>(value));
  287. if (dims == 1) {
  288. ref = ref.squeeze();
  289. } else if (ref.ndim() == 1) {
  290. buf = buf.squeeze();
  291. }
  292. int result = detail::npy_api::get().PyArray_CopyInto_(ref.ptr(), buf.ptr());
  293. if (result < 0) { // Copy failed!
  294. PyErr_Clear();
  295. return false;
  296. }
  297. return true;
  298. }
  299. private:
  300. // Cast implementation
  301. template <typename CType>
  302. static handle cast_impl(CType *src, return_value_policy policy, handle parent) {
  303. switch (policy) {
  304. case return_value_policy::take_ownership:
  305. case return_value_policy::automatic:
  306. return eigen_encapsulate<props>(src);
  307. case return_value_policy::move:
  308. return eigen_encapsulate<props>(new CType(std::move(*src)));
  309. case return_value_policy::copy:
  310. return eigen_array_cast<props>(*src);
  311. case return_value_policy::reference:
  312. case return_value_policy::automatic_reference:
  313. return eigen_ref_array<props>(*src);
  314. case return_value_policy::reference_internal:
  315. return eigen_ref_array<props>(*src, parent);
  316. default:
  317. throw cast_error("unhandled return_value_policy: should not happen!");
  318. };
  319. }
  320. public:
  321. // Normal returned non-reference, non-const value:
  322. static handle cast(Type &&src, return_value_policy /* policy */, handle parent) {
  323. return cast_impl(&src, return_value_policy::move, parent);
  324. }
  325. // If you return a non-reference const, we mark the numpy array readonly:
  326. static handle cast(const Type &&src, return_value_policy /* policy */, handle parent) {
  327. return cast_impl(&src, return_value_policy::move, parent);
  328. }
  329. // lvalue reference return; default (automatic) becomes copy
  330. static handle cast(Type &src, return_value_policy policy, handle parent) {
  331. if (policy == return_value_policy::automatic
  332. || policy == return_value_policy::automatic_reference) {
  333. policy = return_value_policy::copy;
  334. }
  335. return cast_impl(&src, policy, parent);
  336. }
  337. // const lvalue reference return; default (automatic) becomes copy
  338. static handle cast(const Type &src, return_value_policy policy, handle parent) {
  339. if (policy == return_value_policy::automatic
  340. || policy == return_value_policy::automatic_reference) {
  341. policy = return_value_policy::copy;
  342. }
  343. return cast(&src, policy, parent);
  344. }
  345. // non-const pointer return
  346. static handle cast(Type *src, return_value_policy policy, handle parent) {
  347. return cast_impl(src, policy, parent);
  348. }
  349. // const pointer return
  350. static handle cast(const Type *src, return_value_policy policy, handle parent) {
  351. return cast_impl(src, policy, parent);
  352. }
  353. static constexpr auto name = props::descriptor;
  354. // NOLINTNEXTLINE(google-explicit-constructor)
  355. operator Type *() { return &value; }
  356. // NOLINTNEXTLINE(google-explicit-constructor)
  357. operator Type &() { return value; }
  358. // NOLINTNEXTLINE(google-explicit-constructor)
  359. operator Type &&() && { return std::move(value); }
  360. template <typename T>
  361. using cast_op_type = movable_cast_op_type<T>;
  362. private:
  363. Type value;
  364. };
  365. // Base class for casting reference/map/block/etc. objects back to python.
  366. template <typename MapType>
  367. struct eigen_map_caster {
  368. static_assert(!std::is_pointer<typename MapType::Scalar>::value,
  369. PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED);
  370. private:
  371. using props = EigenProps<MapType>;
  372. public:
  373. // Directly referencing a ref/map's data is a bit dangerous (whatever the map/ref points to has
  374. // to stay around), but we'll allow it under the assumption that you know what you're doing
  375. // (and have an appropriate keep_alive in place). We return a numpy array pointing directly at
  376. // the ref's data (The numpy array ends up read-only if the ref was to a const matrix type.)
  377. // Note that this means you need to ensure you don't destroy the object in some other way (e.g.
  378. // with an appropriate keep_alive, or with a reference to a statically allocated matrix).
  379. static handle cast(const MapType &src, return_value_policy policy, handle parent) {
  380. switch (policy) {
  381. case return_value_policy::copy:
  382. return eigen_array_cast<props>(src);
  383. case return_value_policy::reference_internal:
  384. return eigen_array_cast<props>(src, parent, is_eigen_mutable_map<MapType>::value);
  385. case return_value_policy::reference:
  386. case return_value_policy::automatic:
  387. case return_value_policy::automatic_reference:
  388. return eigen_array_cast<props>(src, none(), is_eigen_mutable_map<MapType>::value);
  389. default:
  390. // move, take_ownership don't make any sense for a ref/map:
  391. pybind11_fail("Invalid return_value_policy for Eigen Map/Ref/Block type");
  392. }
  393. }
  394. // return_descr forces the use of NDArray instead of ArrayLike in args
  395. // since Ref<...> args can only accept arrays.
  396. static constexpr auto name = return_descr(props::descriptor);
  397. // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return
  398. // types but not bound arguments). We still provide them (with an explicitly delete) so that
  399. // you end up here if you try anyway.
  400. bool load(handle, bool) = delete;
  401. operator MapType() = delete;
  402. template <typename>
  403. using cast_op_type = MapType;
  404. };
  405. // We can return any map-like object (but can only load Refs, specialized next):
  406. template <typename Type>
  407. struct type_caster<Type, enable_if_t<is_eigen_dense_map<Type>::value>> : eigen_map_caster<Type> {};
  408. // Loader for Ref<...> arguments. See the documentation for info on how to make this work without
  409. // copying (it requires some extra effort in many cases).
  410. template <typename PlainObjectType, typename StrideType>
  411. struct type_caster<
  412. Eigen::Ref<PlainObjectType, 0, StrideType>,
  413. enable_if_t<is_eigen_dense_map<Eigen::Ref<PlainObjectType, 0, StrideType>>::value>>
  414. : public eigen_map_caster<Eigen::Ref<PlainObjectType, 0, StrideType>> {
  415. private:
  416. using Type = Eigen::Ref<PlainObjectType, 0, StrideType>;
  417. using props = EigenProps<Type>;
  418. using Scalar = typename props::Scalar;
  419. static_assert(!std::is_pointer<Scalar>::value,
  420. PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED);
  421. using MapType = Eigen::Map<PlainObjectType, 0, StrideType>;
  422. using Array
  423. = array_t<Scalar,
  424. array::forcecast
  425. | ((props::row_major ? props::inner_stride : props::outer_stride) == 1
  426. ? array::c_style
  427. : (props::row_major ? props::outer_stride : props::inner_stride) == 1
  428. ? array::f_style
  429. : 0)>;
  430. static constexpr bool need_writeable = is_eigen_mutable_map<Type>::value;
  431. // Delay construction (these have no default constructor)
  432. std::unique_ptr<MapType> map;
  433. std::unique_ptr<Type> ref;
  434. // Our array. When possible, this is just a numpy array pointing to the source data, but
  435. // sometimes we can't avoid copying (e.g. input is not a numpy array at all, has an
  436. // incompatible layout, or is an array of a type that needs to be converted). Using a numpy
  437. // temporary (rather than an Eigen temporary) saves an extra copy when we need both type
  438. // conversion and storage order conversion. (Note that we refuse to use this temporary copy
  439. // when loading an argument for a Ref<M> with M non-const, i.e. a read-write reference).
  440. Array copy_or_ref;
  441. public:
  442. bool load(handle src, bool convert) {
  443. // First check whether what we have is already an array of the right type. If not, we
  444. // can't avoid a copy (because the copy is also going to do type conversion).
  445. bool need_copy = !isinstance<Array>(src);
  446. EigenConformable<props::row_major> fits;
  447. if (!need_copy) {
  448. // We don't need a converting copy, but we also need to check whether the strides are
  449. // compatible with the Ref's stride requirements
  450. auto aref = reinterpret_borrow<Array>(src);
  451. if (aref && (!need_writeable || aref.writeable())) {
  452. fits = props::conformable(aref);
  453. if (!fits) {
  454. return false; // Incompatible dimensions
  455. }
  456. if (!fits.template stride_compatible<props>()) {
  457. need_copy = true;
  458. } else {
  459. copy_or_ref = std::move(aref);
  460. }
  461. } else {
  462. need_copy = true;
  463. }
  464. }
  465. if (need_copy) {
  466. // We need to copy: If we need a mutable reference, or we're not supposed to convert
  467. // (either because we're in the no-convert overload pass, or because we're explicitly
  468. // instructed not to copy (via `py::arg().noconvert()`) we have to fail loading.
  469. if (!convert || need_writeable) {
  470. return false;
  471. }
  472. Array copy = Array::ensure(src);
  473. if (!copy) {
  474. return false;
  475. }
  476. fits = props::conformable(copy);
  477. if (!fits || !fits.template stride_compatible<props>()) {
  478. return false;
  479. }
  480. copy_or_ref = std::move(copy);
  481. loader_life_support::add_patient(copy_or_ref);
  482. }
  483. ref.reset();
  484. map.reset(new MapType(data(copy_or_ref),
  485. fits.rows,
  486. fits.cols,
  487. make_stride(fits.stride.outer(), fits.stride.inner())));
  488. ref.reset(new Type(*map));
  489. return true;
  490. }
  491. // NOLINTNEXTLINE(google-explicit-constructor)
  492. operator Type *() { return ref.get(); }
  493. // NOLINTNEXTLINE(google-explicit-constructor)
  494. operator Type &() { return *ref; }
  495. template <typename _T>
  496. using cast_op_type = pybind11::detail::cast_op_type<_T>;
  497. private:
  498. template <typename T = Type, enable_if_t<is_eigen_mutable_map<T>::value, int> = 0>
  499. Scalar *data(Array &a) {
  500. return a.mutable_data();
  501. }
  502. template <typename T = Type, enable_if_t<!is_eigen_mutable_map<T>::value, int> = 0>
  503. const Scalar *data(Array &a) {
  504. return a.data();
  505. }
  506. // Attempt to figure out a constructor of `Stride` that will work.
  507. // If both strides are fixed, use a default constructor:
  508. template <typename S>
  509. using stride_ctor_default = bool_constant<S::InnerStrideAtCompileTime != Eigen::Dynamic
  510. && S::OuterStrideAtCompileTime != Eigen::Dynamic
  511. && std::is_default_constructible<S>::value>;
  512. // Otherwise, if there is a two-index constructor, assume it is (outer,inner) like
  513. // Eigen::Stride, and use it:
  514. template <typename S>
  515. using stride_ctor_dual
  516. = bool_constant<!stride_ctor_default<S>::value
  517. && std::is_constructible<S, EigenIndex, EigenIndex>::value>;
  518. // Otherwise, if there is a one-index constructor, and just one of the strides is dynamic, use
  519. // it (passing whichever stride is dynamic).
  520. template <typename S>
  521. using stride_ctor_outer
  522. = bool_constant<!any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value
  523. && S::OuterStrideAtCompileTime == Eigen::Dynamic
  524. && S::InnerStrideAtCompileTime != Eigen::Dynamic
  525. && std::is_constructible<S, EigenIndex>::value>;
  526. template <typename S>
  527. using stride_ctor_inner
  528. = bool_constant<!any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value
  529. && S::InnerStrideAtCompileTime == Eigen::Dynamic
  530. && S::OuterStrideAtCompileTime != Eigen::Dynamic
  531. && std::is_constructible<S, EigenIndex>::value>;
  532. template <typename S = StrideType, enable_if_t<stride_ctor_default<S>::value, int> = 0>
  533. static S make_stride(EigenIndex, EigenIndex) {
  534. return S();
  535. }
  536. template <typename S = StrideType, enable_if_t<stride_ctor_dual<S>::value, int> = 0>
  537. static S make_stride(EigenIndex outer, EigenIndex inner) {
  538. return S(outer, inner);
  539. }
  540. template <typename S = StrideType, enable_if_t<stride_ctor_outer<S>::value, int> = 0>
  541. static S make_stride(EigenIndex outer, EigenIndex) {
  542. return S(outer);
  543. }
  544. template <typename S = StrideType, enable_if_t<stride_ctor_inner<S>::value, int> = 0>
  545. static S make_stride(EigenIndex, EigenIndex inner) {
  546. return S(inner);
  547. }
  548. };
  549. // type_caster for special matrix types (e.g. DiagonalMatrix), which are EigenBase, but not
  550. // EigenDense (i.e. they don't have a data(), at least not with the usual matrix layout).
  551. // load() is not supported, but we can cast them into the python domain by first copying to a
  552. // regular Eigen::Matrix, then casting that.
  553. template <typename Type>
  554. struct type_caster<Type, enable_if_t<is_eigen_other<Type>::value>> {
  555. static_assert(!std::is_pointer<typename Type::Scalar>::value,
  556. PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED);
  557. protected:
  558. using Matrix
  559. = Eigen::Matrix<typename Type::Scalar, Type::RowsAtCompileTime, Type::ColsAtCompileTime>;
  560. using props = EigenProps<Matrix>;
  561. public:
  562. static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) {
  563. handle h = eigen_encapsulate<props>(new Matrix(src));
  564. return h;
  565. }
  566. static handle cast(const Type *src, return_value_policy policy, handle parent) {
  567. return cast(*src, policy, parent);
  568. }
  569. static constexpr auto name = props::descriptor;
  570. // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return
  571. // types but not bound arguments). We still provide them (with an explicitly delete) so that
  572. // you end up here if you try anyway.
  573. bool load(handle, bool) = delete;
  574. operator Type() = delete;
  575. template <typename>
  576. using cast_op_type = Type;
  577. };
  578. template <typename Type>
  579. struct type_caster<Type, enable_if_t<is_eigen_sparse<Type>::value>> {
  580. using Scalar = typename Type::Scalar;
  581. static_assert(!std::is_pointer<Scalar>::value,
  582. PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED);
  583. using StorageIndex = remove_reference_t<decltype(*std::declval<Type>().outerIndexPtr())>;
  584. using Index = typename Type::Index;
  585. static constexpr bool rowMajor = Type::IsRowMajor;
  586. bool load(handle src, bool) {
  587. if (!src) {
  588. return false;
  589. }
  590. auto obj = reinterpret_borrow<object>(src);
  591. object sparse_module = module_::import("scipy.sparse");
  592. object matrix_type = sparse_module.attr(rowMajor ? "csr_matrix" : "csc_matrix");
  593. if (!type::handle_of(obj).is(matrix_type)) {
  594. try {
  595. obj = matrix_type(obj);
  596. } catch (const error_already_set &) {
  597. return false;
  598. }
  599. }
  600. auto values = array_t<Scalar>((object) obj.attr("data"));
  601. auto innerIndices = array_t<StorageIndex>((object) obj.attr("indices"));
  602. auto outerIndices = array_t<StorageIndex>((object) obj.attr("indptr"));
  603. auto shape = pybind11::tuple((pybind11::object) obj.attr("shape"));
  604. auto nnz = obj.attr("nnz").cast<Index>();
  605. if (!values || !innerIndices || !outerIndices) {
  606. return false;
  607. }
  608. value = EigenMapSparseMatrix<Scalar,
  609. Type::Flags &(Eigen::RowMajor | Eigen::ColMajor),
  610. StorageIndex>(shape[0].cast<Index>(),
  611. shape[1].cast<Index>(),
  612. std::move(nnz),
  613. outerIndices.mutable_data(),
  614. innerIndices.mutable_data(),
  615. values.mutable_data());
  616. return true;
  617. }
  618. static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) {
  619. const_cast<Type &>(src).makeCompressed();
  620. object matrix_type
  621. = module_::import("scipy.sparse").attr(rowMajor ? "csr_matrix" : "csc_matrix");
  622. array data(src.nonZeros(), src.valuePtr());
  623. array outerIndices((rowMajor ? src.rows() : src.cols()) + 1, src.outerIndexPtr());
  624. array innerIndices(src.nonZeros(), src.innerIndexPtr());
  625. return matrix_type(pybind11::make_tuple(
  626. std::move(data), std::move(innerIndices), std::move(outerIndices)),
  627. pybind11::make_tuple(src.rows(), src.cols()))
  628. .release();
  629. }
  630. PYBIND11_TYPE_CASTER(Type,
  631. const_name<(Type::IsRowMajor) != 0>("scipy.sparse.csr_matrix[",
  632. "scipy.sparse.csc_matrix[")
  633. + npy_format_descriptor<Scalar>::name + const_name("]"));
  634. };
  635. PYBIND11_NAMESPACE_END(detail)
  636. PYBIND11_NAMESPACE_END(PYBIND11_NAMESPACE)