/* pybind11/eigen/matrix.h: Transparent conversion for dense and sparse Eigen matrices Copyright (c) 2016 Wenzel Jakob All rights reserved. Use of this source code is governed by a BSD-style license that can be found in the LICENSE file. */ #pragma once #include "../numpy.h" #include "common.h" /* HINT: To suppress warnings originating from the Eigen headers, use -isystem. See also: https://stackoverflow.com/questions/2579576/i-dir-vs-isystem-dir https://stackoverflow.com/questions/1741816/isystem-for-ms-visual-studio-c-compiler */ PYBIND11_WARNING_PUSH PYBIND11_WARNING_DISABLE_MSVC(5054) // https://github.com/pybind/pybind11/pull/3741 // C5054: operator '&': deprecated between enumerations of different types #if defined(__MINGW32__) PYBIND11_WARNING_DISABLE_GCC("-Wmaybe-uninitialized") #endif #include #include PYBIND11_WARNING_POP // Eigen prior to 3.2.7 doesn't have proper move constructors--but worse, some classes get implicit // move constructors that break things. We could detect this an explicitly copy, but an extra copy // of matrices seems highly undesirable. static_assert(EIGEN_VERSION_AT_LEAST(3, 2, 7), "Eigen matrix support in pybind11 requires Eigen >= 3.2.7"); PYBIND11_NAMESPACE_BEGIN(PYBIND11_NAMESPACE) PYBIND11_WARNING_DISABLE_MSVC(4127) // Provide a convenience alias for easier pass-by-ref usage with fully dynamic strides: using EigenDStride = Eigen::Stride; template using EigenDRef = Eigen::Ref; template using EigenDMap = Eigen::Map; PYBIND11_NAMESPACE_BEGIN(detail) #if EIGEN_VERSION_AT_LEAST(3, 3, 0) using EigenIndex = Eigen::Index; template using EigenMapSparseMatrix = Eigen::Map>; #else using EigenIndex = EIGEN_DEFAULT_DENSE_INDEX_TYPE; template using EigenMapSparseMatrix = Eigen::MappedSparseMatrix; #endif // Matches Eigen::Map, Eigen::Ref, blocks, etc: template using is_eigen_dense_map = all_of, std::is_base_of, T>>; template using is_eigen_mutable_map = std::is_base_of, T>; template using is_eigen_dense_plain = all_of>, is_template_base_of>; template using is_eigen_sparse = is_template_base_of; // 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 using is_eigen_other = all_of, negation, is_eigen_dense_plain, is_eigen_sparse>>>; // Captures numpy/eigen conformability status (returned by EigenProps::conformable()): template struct EigenConformable { bool conformable = false; EigenIndex rows = 0, cols = 0; EigenDStride stride{0, 0}; // Only valid if negativestrides is false! bool negativestrides = false; // If true, do not use stride! // NOLINTNEXTLINE(google-explicit-constructor) EigenConformable(bool fits = false) : conformable{fits} {} // Matrix type: EigenConformable(EigenIndex r, EigenIndex c, EigenIndex rstride, EigenIndex cstride) : conformable{true}, rows{r}, cols{c}, // TODO: when Eigen bug #747 is fixed, remove the tests for non-negativity. // http://eigen.tuxfamily.org/bz/show_bug.cgi?id=747 stride{EigenRowMajor ? (rstride > 0 ? rstride : 0) : (cstride > 0 ? cstride : 0) /* outer stride */, EigenRowMajor ? (cstride > 0 ? cstride : 0) : (rstride > 0 ? rstride : 0) /* inner stride */}, negativestrides{rstride < 0 || cstride < 0} {} // Vector type: EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride) : EigenConformable(r, c, r == 1 ? c * stride : stride, c == 1 ? r : r * stride) {} template bool stride_compatible() const { // To have compatible strides, we need (on both dimensions) one of fully dynamic strides, // matching strides, or a dimension size of 1 (in which case the stride value is // irrelevant). Alternatively, if any dimension size is 0, the strides are not relevant // (and numpy ≥ 1.23 sets the strides to 0 in that case, so we need to check explicitly). if (negativestrides) { return false; } if (rows == 0 || cols == 0) { return true; } return (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner() || (EigenRowMajor ? cols : rows) == 1) && (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer() || (EigenRowMajor ? rows : cols) == 1); } // NOLINTNEXTLINE(google-explicit-constructor) operator bool() const { return conformable; } }; template struct eigen_extract_stride { using type = Type; }; template struct eigen_extract_stride> { using type = StrideType; }; template struct eigen_extract_stride> { using type = StrideType; }; // Helper struct for extracting information from an Eigen type template struct EigenProps { using Type = Type_; using Scalar = typename Type::Scalar; using StrideType = typename eigen_extract_stride::type; static constexpr EigenIndex rows = Type::RowsAtCompileTime, cols = Type::ColsAtCompileTime, size = Type::SizeAtCompileTime; static constexpr bool row_major = Type::IsRowMajor, vector = Type::IsVectorAtCompileTime, // At least one dimension has fixed size 1 fixed_rows = rows != Eigen::Dynamic, fixed_cols = cols != Eigen::Dynamic, fixed = size != Eigen::Dynamic, // Fully-fixed size dynamic = !fixed_rows && !fixed_cols; // Fully-dynamic size template using if_zero = std::integral_constant; static constexpr EigenIndex inner_stride = if_zero::value, outer_stride = if_zero < StrideType::OuterStrideAtCompileTime, vector ? size : row_major ? cols : rows > ::value; static constexpr bool dynamic_stride = inner_stride == Eigen::Dynamic && outer_stride == Eigen::Dynamic; static constexpr bool requires_row_major = !dynamic_stride && !vector && (row_major ? inner_stride : outer_stride) == 1; static constexpr bool requires_col_major = !dynamic_stride && !vector && (row_major ? outer_stride : inner_stride) == 1; // Takes an input array and determines whether we can make it fit into the Eigen type. If // the array is a vector, we attempt to fit it into either an Eigen 1xN or Nx1 vector // (preferring the latter if it will fit in either, i.e. for a fully dynamic matrix type). static EigenConformable conformable(const array &a) { const auto dims = a.ndim(); if (dims < 1 || dims > 2) { return false; } if (dims == 2) { // Matrix type: require exact match (or dynamic) EigenIndex np_rows = a.shape(0), np_cols = a.shape(1), np_rstride = a.strides(0) / static_cast(sizeof(Scalar)), np_cstride = a.strides(1) / static_cast(sizeof(Scalar)); if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols)) { return false; } return {np_rows, np_cols, np_rstride, np_cstride}; } // Otherwise we're storing an n-vector. Only one of the strides will be used, but // whichever is used, we want the (single) numpy stride value. const EigenIndex n = a.shape(0), stride = a.strides(0) / static_cast(sizeof(Scalar)); if (vector) { // Eigen type is a compile-time vector if (fixed && size != n) { return false; // Vector size mismatch } return {rows == 1 ? 1 : n, cols == 1 ? 1 : n, stride}; } if (fixed) { // The type has a fixed size, but is not a vector: abort return false; } if (fixed_cols) { // Since this isn't a vector, cols must be != 1. We allow this only if it exactly // equals the number of elements (rows is Dynamic, and so 1 row is allowed). if (cols != n) { return false; } return {1, n, stride}; } // Otherwise it's either fully dynamic, or column dynamic; both become a column vector if (fixed_rows && rows != n) { return false; } return {n, 1, stride}; } static constexpr bool show_writeable = is_eigen_dense_map::value && is_eigen_mutable_map::value; static constexpr bool show_order = is_eigen_dense_map::value; static constexpr bool show_c_contiguous = show_order && requires_row_major; static constexpr bool show_f_contiguous = !show_c_contiguous && show_order && requires_col_major; static constexpr auto descriptor = const_name("numpy.ndarray[") + npy_format_descriptor::name + const_name("[") + const_name(const_name<(size_t) rows>(), const_name("m")) + const_name(", ") + const_name(const_name<(size_t) cols>(), const_name("n")) + const_name("]") + // For a reference type (e.g. Ref) we have other constraints that might need to // be satisfied: writeable=True (for a mutable reference), and, depending on the map's // stride options, possibly f_contiguous or c_contiguous. We include them in the // descriptor output to provide some hint as to why a TypeError is occurring (otherwise // it can be confusing to see that a function accepts a 'numpy.ndarray[float64[3,2]]' and // an error message that you *gave* a numpy.ndarray of the right type and dimensions. const_name(", flags.writeable", "") + const_name(", flags.c_contiguous", "") + const_name(", flags.f_contiguous", "") + const_name("]"); }; // Casts an Eigen type to numpy array. If given a base, the numpy array references the src data, // otherwise it'll make a copy. writeable lets you turn off the writeable flag for the array. template handle eigen_array_cast(typename props::Type const &src, handle base = handle(), bool writeable = true) { constexpr ssize_t elem_size = sizeof(typename props::Scalar); array a; if (props::vector) { a = array({src.size()}, {elem_size * src.innerStride()}, src.data(), base); } else { a = array({src.rows(), src.cols()}, {elem_size * src.rowStride(), elem_size * src.colStride()}, src.data(), base); } if (!writeable) { array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; } return a.release(); } // Takes an lvalue ref to some Eigen type and a (python) base object, creating a numpy array that // reference the Eigen object's data with `base` as the python-registered base class (if omitted, // the base will be set to None, and lifetime management is up to the caller). The numpy array is // non-writeable if the given type is const. template handle eigen_ref_array(Type &src, handle parent = none()) { // none here is to get past array's should-we-copy detection, which currently always // copies when there is no base. Setting the base to None should be harmless. return eigen_array_cast(src, parent, !std::is_const::value); } // Takes a pointer to some dense, plain Eigen type, builds a capsule around it, then returns a // numpy array that references the encapsulated data with a python-side reference to the capsule to // tie its destruction to that of any dependent python objects. Const-ness is determined by // whether or not the Type of the pointer given is const. template ::value>> handle eigen_encapsulate(Type *src) { capsule base(src, [](void *o) { delete static_cast(o); }); return eigen_ref_array(*src, base); } // Type caster for regular, dense matrix types (e.g. MatrixXd), but not maps/refs/etc. of dense // types. template struct type_caster::value>> { using Scalar = typename Type::Scalar; static_assert(!std::is_pointer::value, PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED); using props = EigenProps; bool load(handle src, bool convert) { // If we're in no-convert mode, only load if given an array of the correct type if (!convert && !isinstance>(src)) { return false; } // Coerce into an array, but don't do type conversion yet; the copy below handles it. auto buf = array::ensure(src); if (!buf) { return false; } auto dims = buf.ndim(); if (dims < 1 || dims > 2) { return false; } auto fits = props::conformable(buf); if (!fits) { return false; } // Allocate the new type, then build a numpy reference into it value = Type(fits.rows, fits.cols); auto ref = reinterpret_steal(eigen_ref_array(value)); if (dims == 1) { ref = ref.squeeze(); } else if (ref.ndim() == 1) { buf = buf.squeeze(); } int result = detail::npy_api::get().PyArray_CopyInto_(ref.ptr(), buf.ptr()); if (result < 0) { // Copy failed! PyErr_Clear(); return false; } return true; } private: // Cast implementation template static handle cast_impl(CType *src, return_value_policy policy, handle parent) { switch (policy) { case return_value_policy::take_ownership: case return_value_policy::automatic: return eigen_encapsulate(src); case return_value_policy::move: return eigen_encapsulate(new CType(std::move(*src))); case return_value_policy::copy: return eigen_array_cast(*src); case return_value_policy::reference: case return_value_policy::automatic_reference: return eigen_ref_array(*src); case return_value_policy::reference_internal: return eigen_ref_array(*src, parent); default: throw cast_error("unhandled return_value_policy: should not happen!"); }; } public: // Normal returned non-reference, non-const value: static handle cast(Type &&src, return_value_policy /* policy */, handle parent) { return cast_impl(&src, return_value_policy::move, parent); } // If you return a non-reference const, we mark the numpy array readonly: static handle cast(const Type &&src, return_value_policy /* policy */, handle parent) { return cast_impl(&src, return_value_policy::move, parent); } // lvalue reference return; default (automatic) becomes copy static handle cast(Type &src, return_value_policy policy, handle parent) { if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference) { policy = return_value_policy::copy; } return cast_impl(&src, policy, parent); } // const lvalue reference return; default (automatic) becomes copy static handle cast(const Type &src, return_value_policy policy, handle parent) { if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference) { policy = return_value_policy::copy; } return cast(&src, policy, parent); } // non-const pointer return static handle cast(Type *src, return_value_policy policy, handle parent) { return cast_impl(src, policy, parent); } // const pointer return static handle cast(const Type *src, return_value_policy policy, handle parent) { return cast_impl(src, policy, parent); } static constexpr auto name = props::descriptor; // NOLINTNEXTLINE(google-explicit-constructor) operator Type *() { return &value; } // NOLINTNEXTLINE(google-explicit-constructor) operator Type &() { return value; } // NOLINTNEXTLINE(google-explicit-constructor) operator Type &&() && { return std::move(value); } template using cast_op_type = movable_cast_op_type; private: Type value; }; // Base class for casting reference/map/block/etc. objects back to python. template struct eigen_map_caster { static_assert(!std::is_pointer::value, PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED); private: using props = EigenProps; public: // Directly referencing a ref/map's data is a bit dangerous (whatever the map/ref points to has // to stay around), but we'll allow it under the assumption that you know what you're doing // (and have an appropriate keep_alive in place). We return a numpy array pointing directly at // the ref's data (The numpy array ends up read-only if the ref was to a const matrix type.) // Note that this means you need to ensure you don't destroy the object in some other way (e.g. // with an appropriate keep_alive, or with a reference to a statically allocated matrix). static handle cast(const MapType &src, return_value_policy policy, handle parent) { switch (policy) { case return_value_policy::copy: return eigen_array_cast(src); case return_value_policy::reference_internal: return eigen_array_cast(src, parent, is_eigen_mutable_map::value); case return_value_policy::reference: case return_value_policy::automatic: case return_value_policy::automatic_reference: return eigen_array_cast(src, none(), is_eigen_mutable_map::value); default: // move, take_ownership don't make any sense for a ref/map: pybind11_fail("Invalid return_value_policy for Eigen Map/Ref/Block type"); } } static constexpr auto name = props::descriptor; // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return // types but not bound arguments). We still provide them (with an explicitly delete) so that // you end up here if you try anyway. bool load(handle, bool) = delete; operator MapType() = delete; template using cast_op_type = MapType; }; // We can return any map-like object (but can only load Refs, specialized next): template struct type_caster::value>> : eigen_map_caster {}; // Loader for Ref<...> arguments. See the documentation for info on how to make this work without // copying (it requires some extra effort in many cases). template struct type_caster< Eigen::Ref, enable_if_t>::value>> : public eigen_map_caster> { private: using Type = Eigen::Ref; using props = EigenProps; using Scalar = typename props::Scalar; static_assert(!std::is_pointer::value, PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED); using MapType = Eigen::Map; using Array = array_t; static constexpr bool need_writeable = is_eigen_mutable_map::value; // Delay construction (these have no default constructor) std::unique_ptr map; std::unique_ptr ref; // Our array. When possible, this is just a numpy array pointing to the source data, but // sometimes we can't avoid copying (e.g. input is not a numpy array at all, has an // incompatible layout, or is an array of a type that needs to be converted). Using a numpy // temporary (rather than an Eigen temporary) saves an extra copy when we need both type // conversion and storage order conversion. (Note that we refuse to use this temporary copy // when loading an argument for a Ref with M non-const, i.e. a read-write reference). Array copy_or_ref; public: bool load(handle src, bool convert) { // First check whether what we have is already an array of the right type. If not, we // can't avoid a copy (because the copy is also going to do type conversion). bool need_copy = !isinstance(src); EigenConformable fits; if (!need_copy) { // We don't need a converting copy, but we also need to check whether the strides are // compatible with the Ref's stride requirements auto aref = reinterpret_borrow(src); if (aref && (!need_writeable || aref.writeable())) { fits = props::conformable(aref); if (!fits) { return false; // Incompatible dimensions } if (!fits.template stride_compatible()) { need_copy = true; } else { copy_or_ref = std::move(aref); } } else { need_copy = true; } } if (need_copy) { // We need to copy: If we need a mutable reference, or we're not supposed to convert // (either because we're in the no-convert overload pass, or because we're explicitly // instructed not to copy (via `py::arg().noconvert()`) we have to fail loading. if (!convert || need_writeable) { return false; } Array copy = Array::ensure(src); if (!copy) { return false; } fits = props::conformable(copy); if (!fits || !fits.template stride_compatible()) { return false; } copy_or_ref = std::move(copy); loader_life_support::add_patient(copy_or_ref); } ref.reset(); map.reset(new MapType(data(copy_or_ref), fits.rows, fits.cols, make_stride(fits.stride.outer(), fits.stride.inner()))); ref.reset(new Type(*map)); return true; } // NOLINTNEXTLINE(google-explicit-constructor) operator Type *() { return ref.get(); } // NOLINTNEXTLINE(google-explicit-constructor) operator Type &() { return *ref; } template using cast_op_type = pybind11::detail::cast_op_type<_T>; private: template ::value, int> = 0> Scalar *data(Array &a) { return a.mutable_data(); } template ::value, int> = 0> const Scalar *data(Array &a) { return a.data(); } // Attempt to figure out a constructor of `Stride` that will work. // If both strides are fixed, use a default constructor: template using stride_ctor_default = bool_constant::value>; // Otherwise, if there is a two-index constructor, assume it is (outer,inner) like // Eigen::Stride, and use it: template using stride_ctor_dual = bool_constant::value && std::is_constructible::value>; // Otherwise, if there is a one-index constructor, and just one of the strides is dynamic, use // it (passing whichever stride is dynamic). template using stride_ctor_outer = bool_constant, stride_ctor_dual>::value && S::OuterStrideAtCompileTime == Eigen::Dynamic && S::InnerStrideAtCompileTime != Eigen::Dynamic && std::is_constructible::value>; template using stride_ctor_inner = bool_constant, stride_ctor_dual>::value && S::InnerStrideAtCompileTime == Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic && std::is_constructible::value>; template ::value, int> = 0> static S make_stride(EigenIndex, EigenIndex) { return S(); } template ::value, int> = 0> static S make_stride(EigenIndex outer, EigenIndex inner) { return S(outer, inner); } template ::value, int> = 0> static S make_stride(EigenIndex outer, EigenIndex) { return S(outer); } template ::value, int> = 0> static S make_stride(EigenIndex, EigenIndex inner) { return S(inner); } }; // type_caster for special matrix types (e.g. DiagonalMatrix), which are EigenBase, but not // EigenDense (i.e. they don't have a data(), at least not with the usual matrix layout). // 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>> { static_assert(!std::is_pointer::value, PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED); protected: using Matrix = Eigen::Matrix; using props = EigenProps; public: static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) { handle h = eigen_encapsulate(new Matrix(src)); return h; } static handle cast(const Type *src, return_value_policy policy, handle parent) { return cast(*src, policy, parent); } static constexpr auto name = props::descriptor; // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return // types but not bound arguments). We still provide them (with an explicitly delete) so that // you end up here if you try anyway. bool load(handle, bool) = delete; operator Type() = delete; template using cast_op_type = Type; }; template struct type_caster::value>> { using Scalar = typename Type::Scalar; static_assert(!std::is_pointer::value, PYBIND11_EIGEN_MESSAGE_POINTER_TYPES_ARE_NOT_SUPPORTED); using StorageIndex = remove_reference_t().outerIndexPtr())>; using Index = typename Type::Index; static constexpr bool rowMajor = Type::IsRowMajor; bool load(handle src, bool) { if (!src) { return false; } auto obj = reinterpret_borrow(src); object sparse_module = module_::import("scipy.sparse"); object matrix_type = sparse_module.attr(rowMajor ? "csr_matrix" : "csc_matrix"); if (!type::handle_of(obj).is(matrix_type)) { try { obj = matrix_type(obj); } catch (const error_already_set &) { return false; } } auto values = array_t((object) obj.attr("data")); auto innerIndices = array_t((object) obj.attr("indices")); auto outerIndices = array_t((object) obj.attr("indptr")); auto shape = pybind11::tuple((pybind11::object) obj.attr("shape")); auto nnz = obj.attr("nnz").cast(); if (!values || !innerIndices || !outerIndices) { return false; } value = EigenMapSparseMatrix(shape[0].cast(), shape[1].cast(), std::move(nnz), outerIndices.mutable_data(), innerIndices.mutable_data(), values.mutable_data()); return true; } static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) { const_cast(src).makeCompressed(); object matrix_type = module_::import("scipy.sparse").attr(rowMajor ? "csr_matrix" : "csc_matrix"); array data(src.nonZeros(), src.valuePtr()); array outerIndices((rowMajor ? src.rows() : src.cols()) + 1, src.outerIndexPtr()); array innerIndices(src.nonZeros(), src.innerIndexPtr()); return matrix_type(pybind11::make_tuple( std::move(data), std::move(innerIndices), std::move(outerIndices)), pybind11::make_tuple(src.rows(), src.cols())) .release(); } PYBIND11_TYPE_CASTER(Type, const_name<(Type::IsRowMajor) != 0>("scipy.sparse.csr_matrix[", "scipy.sparse.csc_matrix[") + npy_format_descriptor::name + const_name("]")); }; PYBIND11_NAMESPACE_END(detail) PYBIND11_NAMESPACE_END(PYBIND11_NAMESPACE)