diff --git a/include/pybind11/detail/eigen_matrix.h b/include/pybind11/detail/eigen_matrix.h new file mode 100644 index 00000000000..c4445ccb522 --- /dev/null +++ b/include/pybind11/detail/eigen_matrix.h @@ -0,0 +1,708 @@ + +/* + pybind11/eigen.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 + +/* 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 +*/ + +#include "../numpy.h" + +// The C4127 suppression was introduced for Eigen 3.4.0. In theory we could +// make it version specific, or even remove it later, but considering that +// 1. C4127 is generally far more distracting than useful for modern template code, and +// 2. we definitely want to ignore any MSVC warnings originating from Eigen code, +// it is probably best to keep this around indefinitely. +#if defined(_MSC_VER) +# pragma warning(push) +# pragma warning(disable : 4127) // C4127: conditional expression is constant +# pragma warning(disable : 5054) // https://github.com/pybind/pybind11/pull/3741 +// C5054: operator '&': deprecated between enumerations of different types +#elif defined(__MINGW32__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" +#endif + +#include +#include + +#if defined(_MSC_VER) +# pragma warning(pop) +#elif defined(__MINGW32__) +# pragma GCC diagnostic pop +#endif + +PYBIND11_NAMESPACE_BEGIN(PYBIND11_NAMESPACE) + +// 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 ((PYBIND11_SILENCE_MSVC_C4127(fixed_rows) && np_rows != rows) + || (PYBIND11_SILENCE_MSVC_C4127(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 (PYBIND11_SILENCE_MSVC_C4127(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 (PYBIND11_SILENCE_MSVC_C4127(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; + 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 { +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; + 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>> { +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; + 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) \ No newline at end of file diff --git a/include/pybind11/detail/eigen_tensor.h b/include/pybind11/detail/eigen_tensor.h new file mode 100644 index 00000000000..c3042101ce0 --- /dev/null +++ b/include/pybind11/detail/eigen_tensor.h @@ -0,0 +1,383 @@ + +/* + pybind11/eigen_tensor.h: Transparent conversion for Eigen tensors + + 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 + +/* 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 +*/ + +#include "../numpy.h" + +// The C4127 suppression was introduced for Eigen 3.4.0. In theory we could +// make it version specific, or even remove it later, but considering that +// 1. C4127 is generally far more distracting than useful for modern template code, and +// 2. we definitely want to ignore any MSVC warnings originating from Eigen code, +// it is probably best to keep this around indefinitely. +#if defined(_MSC_VER) +# pragma warning(push) +# pragma warning(disable : 4554) // Tensor.h warning +// C5054: operator '&': deprecated between enumerations of different types +#elif defined(__MINGW32__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" +#endif + +#include + +#if defined(_MSC_VER) +# pragma warning(pop) +#elif defined(__MINGW32__) +# pragma GCC diagnostic pop +#endif + +PYBIND11_NAMESPACE_BEGIN(PYBIND11_NAMESPACE) + +PYBIND11_NAMESPACE_BEGIN(detail) + +template +constexpr int compute_array_flag_from_tensor() { + static_assert(((int) T::Layout == (int) Eigen::RowMajor) + || ((int) T::Layout == (int) Eigen::ColMajor), + "Layout must be row or column major"); + return ((int) T::Layout == (int) Eigen::RowMajor) ? array::c_style : array::f_style; +} + +template +struct eigen_tensor_helper {}; + +template +struct eigen_tensor_helper> { + using T = Eigen::Tensor; + using ValidType = void; + + static std::array get_shape(const T &f) { + return f.dimensions(); + } + + static constexpr bool + is_correct_shape(const std::array & /*shape*/) { + return true; + } + + template + struct helper {}; + + template + struct helper> { + static constexpr auto value = concat(const_name(((void) Is, "?"))...); + }; + + static constexpr auto dimensions_descriptor + = helper())>::value; +}; + +template +struct eigen_tensor_helper< + Eigen::TensorFixedSize, Options_, IndexType>> { + using T = Eigen::TensorFixedSize, Options_, IndexType>; + using ValidType = void; + + static constexpr std::array get_shape(const T & /*f*/) { + return get_shape(); + } + + static constexpr std::array get_shape() { + return {{Indices...}}; + } + + static bool is_correct_shape(const std::array &shape) { + return get_shape() == shape; + } + + static constexpr auto dimensions_descriptor = concat(const_name()...); +}; + +template +struct get_tensor_descriptor { + static constexpr auto value + = const_name("numpy.ndarray[") + npy_format_descriptor::name + + const_name("[") + eigen_tensor_helper::dimensions_descriptor + + const_name("], flags.writeable, ") + + const_name<(int) T::Layout == (int) Eigen::RowMajor>("flags.c_contiguous", + "flags.f_contiguous"); +}; + +template +struct type_caster::ValidType> { + using H = eigen_tensor_helper; + PYBIND11_TYPE_CASTER(Type, get_tensor_descriptor::value); + + bool load(handle src, bool /*convert*/) { + array_t()> a( + reinterpret_borrow(src)); + + if (a.ndim() != Type::NumIndices) { + return false; + } + + std::array shape; + std::copy(a.shape(), a.shape() + Type::NumIndices, shape.begin()); + + if (!H::is_correct_shape(shape)) { + return false; + } + + value = Eigen::TensorMap(const_cast(a.data()), shape); + + return true; + } + + static handle cast(Type &&src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::reference + || policy == return_value_policy::reference_internal) { + pybind11_fail("Cannot use a reference return value policy for an rvalue"); + } + return cast_impl(&src, return_value_policy::move, parent); + } + + static handle cast(const Type &&src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::reference + || policy == return_value_policy::reference_internal) { + pybind11_fail("Cannot use a reference return value policy for an rvalue"); + } + return cast_impl(&src, return_value_policy::move, parent); + } + + 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); + } + + 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); + } + + static handle cast(Type *src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic) { + policy = return_value_policy::take_ownership; + } else if (policy == return_value_policy::automatic_reference) { + policy = return_value_policy::reference; + } + return cast_impl(src, policy, parent); + } + + static handle cast(const Type *src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic) { + policy = return_value_policy::take_ownership; + } else if (policy == return_value_policy::automatic_reference) { + policy = return_value_policy::reference; + } + return cast_impl(src, policy, parent); + } + + template + static handle cast_impl(C *src, return_value_policy policy, handle parent) { + object parent_object; + bool writeable = false; + switch (policy) { + case return_value_policy::move: + if (std::is_const::value) { + pybind11_fail("Cannot move from a constant reference"); + } + { + Eigen::aligned_allocator allocator; + Type *copy = ::new (allocator.allocate(1)) Type(std::move(*src)); + src = copy; + } + + parent_object = capsule(src, [](void *ptr) { + Eigen::aligned_allocator allocator; + Type *copy = (Type *) ptr; + copy->~Type(); + allocator.deallocate(copy, 1); + }); + writeable = true; + break; + + case return_value_policy::take_ownership: + if (std::is_const::value) { + pybind11_fail("Cannot take ownership of a const reference"); + } + parent_object = capsule(src, [](void *ptr) { delete (Type *) ptr; }); + writeable = true; + break; + + case return_value_policy::copy: + parent_object = {}; + writeable = true; + break; + + case return_value_policy::reference: + parent_object = none(); + writeable = !std::is_const::value; + break; + + case return_value_policy::reference_internal: + // Default should do the right thing + parent_object = reinterpret_borrow(parent); + writeable = !std::is_const::value; + break; + + default: + pybind11_fail("pybind11 bug in eigen.h, please file a bug report"); + } + + handle result = array_t()>( + H::get_shape(*src), src->data(), parent_object) + .release(); + + if (!writeable) { + array_proxy(result.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; + } + + return result; + } +}; + +template +struct type_caster, typename eigen_tensor_helper::ValidType> { + using H = eigen_tensor_helper; + + bool load(handle src, bool /*convert*/) { + // Note that we have a lot more checks here as we want to make sure to avoid copies + auto a = reinterpret_borrow(src); + if ((a.flags() & compute_array_flag_from_tensor()) == 0) { + return false; + } + + if (!a.dtype().is(dtype::of())) { + return false; + } + + if (a.ndim() != Type::NumIndices) { + return false; + } + + std::array shape; + std::copy(a.shape(), a.shape() + Type::NumIndices, shape.begin()); + + if (!H::is_correct_shape(shape)) { + return false; + } + + value.reset(new Eigen::TensorMap( + reinterpret_cast(a.mutable_data()), shape)); + + return true; + } + + static handle cast(Eigen::TensorMap &&src, return_value_policy policy, handle parent) { + return cast_impl(&src, policy, parent); + } + + static handle + cast(const Eigen::TensorMap &&src, return_value_policy policy, handle parent) { + return cast_impl(&src, policy, parent); + } + + static handle cast(Eigen::TensorMap &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); + } + + static handle + cast(const Eigen::TensorMap &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); + } + + static handle cast(Eigen::TensorMap *src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic) { + policy = return_value_policy::take_ownership; + } else if (policy == return_value_policy::automatic_reference) { + policy = return_value_policy::reference; + } + return cast_impl(src, policy, parent); + } + + static handle + cast(const Eigen::TensorMap *src, return_value_policy policy, handle parent) { + if (policy == return_value_policy::automatic) { + policy = return_value_policy::take_ownership; + } else if (policy == return_value_policy::automatic_reference) { + policy = return_value_policy::reference; + } + return cast_impl(src, policy, parent); + } + + template + static handle cast_impl(C *src, return_value_policy policy, handle parent) { + object parent_object; + constexpr bool writeable = !std::is_const::value; + switch (policy) { + case return_value_policy::reference: + parent_object = none(); + break; + + case return_value_policy::reference_internal: + // Default should do the right thing + parent_object = reinterpret_borrow(parent); + break; + + default: + // move, take_ownership don't make any sense for a ref/map: + pybind11_fail("Invalid return_value_policy for Eigen Map type, must be either " + "reference or reference_internal"); + } + + handle result = array_t()>( + H::get_shape(*src), src->data(), parent_object) + .release(); + + if (!writeable) { + array_proxy(result.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; + } + + return result; + } + +protected: + // TODO: Move to std::optional once std::optional has more support + std::unique_ptr> value; + +public: + static constexpr auto name = get_tensor_descriptor::value; + explicit operator Eigen::TensorMap *() { + return value.get(); + } /* NOLINT(bugprone-macro-parentheses) */ + explicit operator Eigen::TensorMap &() { + return *value; + } /* NOLINT(bugprone-macro-parentheses) */ + explicit operator Eigen::TensorMap &&() && { + return std::move(*value); + } /* NOLINT(bugprone-macro-parentheses) */ + + template + using cast_op_type = ::pybind11::detail::movable_cast_op_type; +}; + +PYBIND11_NAMESPACE_END(detail) +PYBIND11_NAMESPACE_END(PYBIND11_NAMESPACE) \ No newline at end of file diff --git a/include/pybind11/eigen.h b/include/pybind11/eigen.h index a6d3a522118..69fb528f991 100644 --- a/include/pybind11/eigen.h +++ b/include/pybind11/eigen.h @@ -9,1044 +9,11 @@ #pragma once -/* 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 -*/ - -#include "numpy.h" - -// The C4127 suppression was introduced for Eigen 3.4.0. In theory we could -// make it version specific, or even remove it later, but considering that -// 1. C4127 is generally far more distracting than useful for modern template code, and -// 2. we definitely want to ignore any MSVC warnings originating from Eigen code, -// it is probably best to keep this around indefinitely. -#if defined(_MSC_VER) -# pragma warning(push) -# pragma warning(disable : 4554) // Tensor.h warning -# pragma warning(disable : 4127) // C4127: conditional expression is constant -# pragma warning(disable : 5054) // https://github.com/pybind/pybind11/pull/3741 -// C5054: operator '&': deprecated between enumerations of different types -#elif defined(__MINGW32__) -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" -#endif - -#include - -#include -#include -#include - -#if defined(_MSC_VER) -# pragma warning(pop) -#elif defined(__MINGW32__) -# pragma GCC diagnostic pop -#endif +#include "detail/eigen_matrix.h" +#include "detail/eigen_tensor.h" // 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 support in pybind11 requires Eigen >= 3.2.7"); - -PYBIND11_NAMESPACE_BEGIN(PYBIND11_NAMESPACE) - -// 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 ((PYBIND11_SILENCE_MSVC_C4127(fixed_rows) && np_rows != rows) - || (PYBIND11_SILENCE_MSVC_C4127(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 (PYBIND11_SILENCE_MSVC_C4127(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 (PYBIND11_SILENCE_MSVC_C4127(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; - 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 { -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; - 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>> { -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 -constexpr int compute_array_flag_from_tensor() { - static_assert(((int) T::Layout == (int) Eigen::RowMajor) - || ((int) T::Layout == (int) Eigen::ColMajor), - "Layout must be row or column major"); - return ((int) T::Layout == (int) Eigen::RowMajor) ? array::c_style : array::f_style; -} - -template -struct eigen_tensor_helper {}; - -template -struct eigen_tensor_helper> { - using T = Eigen::Tensor; - using ValidType = void; - - static std::array get_shape(const T &f) { - return f.dimensions(); - } - - static constexpr bool - is_correct_shape(const std::array & /*shape*/) { - return true; - } - - template - struct helper {}; - - template - struct helper> { - static constexpr auto value = concat(const_name(((void) Is, "?"))...); - }; - - static constexpr auto dimensions_descriptor - = helper())>::value; -}; - -template -struct eigen_tensor_helper< - Eigen::TensorFixedSize, Options_, IndexType>> { - using T = Eigen::TensorFixedSize, Options_, IndexType>; - using ValidType = void; - - static constexpr std::array get_shape(const T & /*f*/) { - return get_shape(); - } - - static constexpr std::array get_shape() { - return {{Indices...}}; - } - - static bool is_correct_shape(const std::array &shape) { - return get_shape() == shape; - } - - static constexpr auto dimensions_descriptor = concat(const_name()...); -}; - -template -struct get_tensor_descriptor { - static constexpr auto value - = const_name("numpy.ndarray[") + npy_format_descriptor::name - + const_name("[") + eigen_tensor_helper::dimensions_descriptor - + const_name("], flags.writeable, ") - + const_name<(int) T::Layout == (int) Eigen::RowMajor>("flags.c_contiguous", - "flags.f_contiguous"); -}; - -template -struct type_caster::ValidType> { - using H = eigen_tensor_helper; - PYBIND11_TYPE_CASTER(Type, get_tensor_descriptor::value); - - bool load(handle src, bool /*convert*/) { - array_t()> a( - reinterpret_borrow(src)); - - if (a.ndim() != Type::NumIndices) { - return false; - } - - std::array shape; - std::copy(a.shape(), a.shape() + Type::NumIndices, shape.begin()); - - if (!H::is_correct_shape(shape)) { - return false; - } - - value = Eigen::TensorMap(const_cast(a.data()), shape); - - return true; - } - - static handle cast(Type &&src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::reference - || policy == return_value_policy::reference_internal) { - pybind11_fail("Cannot use a reference return value policy for an rvalue"); - } - return cast_impl(&src, return_value_policy::move, parent); - } - - static handle cast(const Type &&src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::reference - || policy == return_value_policy::reference_internal) { - pybind11_fail("Cannot use a reference return value policy for an rvalue"); - } - return cast_impl(&src, return_value_policy::move, parent); - } - - 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); - } - - 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); - } - - static handle cast(Type *src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::automatic) { - policy = return_value_policy::take_ownership; - } else if (policy == return_value_policy::automatic_reference) { - policy = return_value_policy::reference; - } - return cast_impl(src, policy, parent); - } - - static handle cast(const Type *src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::automatic) { - policy = return_value_policy::take_ownership; - } else if (policy == return_value_policy::automatic_reference) { - policy = return_value_policy::reference; - } - return cast_impl(src, policy, parent); - } - - template - static handle cast_impl(C *src, return_value_policy policy, handle parent) { - object parent_object; - bool writeable = false; - switch (policy) { - case return_value_policy::move: - if (std::is_const::value) { - pybind11_fail("Cannot move from a constant reference"); - } - { - Eigen::aligned_allocator allocator; - Type *copy = ::new (allocator.allocate(1)) Type(std::move(*src)); - src = copy; - } - - parent_object = capsule(src, [](void *ptr) { - Eigen::aligned_allocator allocator; - Type *copy = (Type *) ptr; - copy->~Type(); - allocator.deallocate(copy, 1); - }); - writeable = true; - break; - - case return_value_policy::take_ownership: - if (std::is_const::value) { - pybind11_fail("Cannot take ownership of a const reference"); - } - parent_object = capsule(src, [](void *ptr) { delete (Type *) ptr; }); - writeable = true; - break; - - case return_value_policy::copy: - parent_object = {}; - writeable = true; - break; - - case return_value_policy::reference: - parent_object = none(); - writeable = !std::is_const::value; - break; - - case return_value_policy::reference_internal: - // Default should do the right thing - parent_object = reinterpret_borrow(parent); - writeable = !std::is_const::value; - break; - - default: - pybind11_fail("pybind11 bug in eigen.h, please file a bug report"); - } - - handle result = array_t()>( - H::get_shape(*src), src->data(), parent_object) - .release(); - - if (!writeable) { - array_proxy(result.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; - } - - return result; - } -}; - -template -struct type_caster, typename eigen_tensor_helper::ValidType> { - using H = eigen_tensor_helper; - - bool load(handle src, bool /*convert*/) { - // Note that we have a lot more checks here as we want to make sure to avoid copies - auto a = reinterpret_borrow(src); - if ((a.flags() & compute_array_flag_from_tensor()) == 0) { - return false; - } - - if (!a.dtype().is(dtype::of())) { - return false; - } - - if (a.ndim() != Type::NumIndices) { - return false; - } - - std::array shape; - std::copy(a.shape(), a.shape() + Type::NumIndices, shape.begin()); - - if (!H::is_correct_shape(shape)) { - return false; - } - - value.reset(new Eigen::TensorMap( - reinterpret_cast(a.mutable_data()), shape)); - - return true; - } - - static handle cast(Eigen::TensorMap &&src, return_value_policy policy, handle parent) { - return cast_impl(&src, policy, parent); - } - - static handle - cast(const Eigen::TensorMap &&src, return_value_policy policy, handle parent) { - return cast_impl(&src, policy, parent); - } - - static handle cast(Eigen::TensorMap &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); - } - - static handle - cast(const Eigen::TensorMap &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); - } - - static handle cast(Eigen::TensorMap *src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::automatic) { - policy = return_value_policy::take_ownership; - } else if (policy == return_value_policy::automatic_reference) { - policy = return_value_policy::reference; - } - return cast_impl(src, policy, parent); - } - - static handle - cast(const Eigen::TensorMap *src, return_value_policy policy, handle parent) { - if (policy == return_value_policy::automatic) { - policy = return_value_policy::take_ownership; - } else if (policy == return_value_policy::automatic_reference) { - policy = return_value_policy::reference; - } - return cast_impl(src, policy, parent); - } - - template - static handle cast_impl(C *src, return_value_policy policy, handle parent) { - object parent_object; - constexpr bool writeable = !std::is_const::value; - switch (policy) { - case return_value_policy::reference: - parent_object = none(); - break; - - case return_value_policy::reference_internal: - // Default should do the right thing - parent_object = reinterpret_borrow(parent); - break; - - default: - // move, take_ownership don't make any sense for a ref/map: - pybind11_fail("Invalid return_value_policy for Eigen Map type, must be either " - "reference or reference_internal"); - } - - handle result = array_t()>( - H::get_shape(*src), src->data(), parent_object) - .release(); - - if (!writeable) { - array_proxy(result.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_; - } - - return result; - } - -protected: - // TODO: Move to std::optional once std::optional has more support - std::unique_ptr> value; - -public: - static constexpr auto name = get_tensor_descriptor::value; - explicit operator Eigen::TensorMap *() { - return value.get(); - } /* NOLINT(bugprone-macro-parentheses) */ - explicit operator Eigen::TensorMap &() { - return *value; - } /* NOLINT(bugprone-macro-parentheses) */ - explicit operator Eigen::TensorMap &&() && { - return std::move(*value); - } /* NOLINT(bugprone-macro-parentheses) */ - - template - using cast_op_type = ::pybind11::detail::movable_cast_op_type; -}; - -template -struct type_caster::value>> { - using Scalar = typename Type::Scalar; - 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) + "Eigen support in pybind11 requires Eigen >= 3.2.7"); \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7296cd1b812..9adac392e08 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -128,7 +128,8 @@ set(PYBIND11_TEST_FILES test_custom_type_casters test_custom_type_setup test_docstring_options - test_eigen + test_eigen_tensor + test_eigen_matrix test_enum test_eval test_exceptions @@ -233,7 +234,10 @@ list(GET PYBIND11_EIGEN_VERSION_AND_HASH 1 PYBIND11_EIGEN_VERSION_HASH) # Check if Eigen is available; if not, remove from PYBIND11_TEST_FILES (but # keep it in PYBIND11_PYTEST_FILES, so that we get the "eigen is not installed" # skip message). -list(FIND PYBIND11_TEST_FILES test_eigen.cpp PYBIND11_TEST_FILES_EIGEN_I) +list(FIND PYBIND11_TEST_FILES test_eigen_tensor.cpp PYBIND11_TEST_FILES_EIGEN_I) +if(PYBIND11_TEST_FILES_EIGEN_I EQUAL -1) + list(FIND PYBIND11_TEST_FILES test_eigen_matrix.cpp PYBIND11_TEST_FILES_EIGEN_I) +endif() if(PYBIND11_TEST_FILES_EIGEN_I GREATER -1) # Try loading via newer Eigen's Eigen3Config first (bypassing tools/FindEigen3.cmake). # Eigen 3.3.1+ exports a cmake 3.0+ target for handling dependency requirements, but also diff --git a/tests/test_eigen.cpp b/tests/test_eigen_matrix.cpp similarity index 89% rename from tests/test_eigen.cpp rename to tests/test_eigen_matrix.cpp index e0bbabc21ce..63d558eb297 100644 --- a/tests/test_eigen.cpp +++ b/tests/test_eigen_matrix.cpp @@ -58,31 +58,6 @@ MatrixXdR &get_rm() { return *x; } -Eigen::Tensor &get_tensor() { - static Eigen::Tensor *x; - - if (!x) { - x = new Eigen::Tensor(3, 1, 2); - reset_tensor(*x); - } - - return *x; -} - -Eigen::TensorFixedSize> &get_fixed_tensor() { - static Eigen::TensorFixedSize> *x; - - if (!x) { - Eigen::aligned_allocator>> allocator; - x = new (allocator.allocate(1)) Eigen::TensorFixedSize>(); - reset_tensor(*x); - } - - return *x; -} - -const Eigen::Tensor &get_const_tensor() { return get_tensor(); } - // Resets the values of the static matrices returned by get_cm()/get_rm() void reset_refs() { reset_ref(get_cm()); @@ -114,7 +89,7 @@ struct CustomOperatorNew { EIGEN_MAKE_ALIGNED_OPERATOR_NEW; }; -TEST_SUBMODULE(eigen, m) { +TEST_SUBMODULE(eigen_matrix, m) { using FixedMatrixR = Eigen::Matrix; using FixedMatrixC = Eigen::Matrix; using DenseMatrixR = Eigen::Matrix; @@ -431,46 +406,4 @@ TEST_SUBMODULE(eigen, m) { py::module_::import("numpy").attr("ones")(10); return v[0](5); }); - - m.def("copy_fixed_global_tensor", []() { return get_fixed_tensor(); }); - - m.def("copy_global_tensor", []() { return get_tensor(); }); - - m.def("copy_const_global_tensor", []() { return get_const_tensor(); }); - - m.def( - "reference_global_tensor", - []() { return &get_tensor(); }, - py::return_value_policy::reference); - - m.def( - "reference_const_global_tensor", - []() { return &get_const_tensor(); }, - py::return_value_policy::reference); - - m.def( - "reference_view_of_global_tensor", - []() { return Eigen::TensorMap>(get_tensor()); }, - py::return_value_policy::reference); - - m.def( - "reference_view_of_fixed_global_tensor", - []() { - return Eigen::TensorMap>>( - get_fixed_tensor()); - }, - py::return_value_policy::reference); - - m.def("round_trip_tensor", [](const py::array_t &foo) { - auto blah = foo.cast>(); - return blah; - }); - - m.def( - "round_trip_view_tensor", - [](const py::array_t &foo) { - auto view = foo.cast>>(); - return view; - }, - py::return_value_policy::reference); } diff --git a/tests/test_eigen.py b/tests/test_eigen_matrix.py similarity index 94% rename from tests/test_eigen.py rename to tests/test_eigen_matrix.py index 0e2bef69671..97f877aaa58 100644 --- a/tests/test_eigen.py +++ b/tests/test_eigen_matrix.py @@ -3,7 +3,7 @@ from pybind11_tests import ConstructorStats np = pytest.importorskip("numpy") -m = pytest.importorskip("pybind11_tests.eigen") +m = pytest.importorskip("pybind11_tests.eigen_matrix") ref = np.array( @@ -772,63 +772,4 @@ def test_custom_operator_new(): o = m.CustomOperatorNew() np.testing.assert_allclose(o.a, 0.0) - np.testing.assert_allclose(o.b.diagonal(), 1.0) - - -tensor_ref = np.array( - [ - [[0, 3]], - [[1, 4]], - [[2, 5]], - ] -) - - -def assert_equal_tensor_ref(mat, writeable=True, modified=0): - assert mat.flags.writeable == writeable - - if modified != 0: - tensor_ref[0, 0, 0] = modified - - np.testing.assert_array_equal(mat, tensor_ref) - - tensor_ref[0, 0, 0] = 0 - - -def test_convert_tensor_to_py(): - assert_equal_tensor_ref(m.copy_global_tensor()) - assert_equal_tensor_ref(m.copy_fixed_global_tensor()) - assert_equal_tensor_ref(m.copy_const_global_tensor()) - - assert_equal_tensor_ref(m.reference_global_tensor()) - assert_equal_tensor_ref(m.reference_view_of_global_tensor()) - assert_equal_tensor_ref(m.reference_view_of_fixed_global_tensor()) - assert_equal_tensor_ref(m.reference_const_global_tensor(), writeable=False) - - -def test_references_actually_refer(): - a = m.reference_global_tensor() - a[0, 0, 0] = 100 - assert_equal_tensor_ref(m.copy_const_global_tensor(), modified=100) - a[0, 0, 0] = 0 - assert_equal_tensor_ref(m.copy_const_global_tensor()) - - a = m.reference_view_of_global_tensor() - a[0, 0, 0] = 100 - assert_equal_tensor_ref(m.copy_const_global_tensor(), modified=100) - a[0, 0, 0] = 0 - assert_equal_tensor_ref(m.copy_const_global_tensor()) - - -def test_round_trip(): - assert_equal_tensor_ref(m.round_trip_tensor(tensor_ref)) - - -def test_round_trip_references_actually_refer(): - # Need to create a copy that matches the type on the C side - copy = np.array(tensor_ref, dtype=np.float64, order="F") - a = m.round_trip_view_tensor(copy) - a[0, 0, 0] = 100 - assert_equal_tensor_ref(copy, modified=100) - a[0, 0, 0] = 0 - assert_equal_tensor_ref(copy) + np.testing.assert_allclose(o.b.diagonal(), 1.0) \ No newline at end of file diff --git a/tests/test_eigen_tensor.cpp b/tests/test_eigen_tensor.cpp new file mode 100644 index 00000000000..c3f174c3ad0 --- /dev/null +++ b/tests/test_eigen_tensor.cpp @@ -0,0 +1,88 @@ +/* + tests/eigen_tensor.cpp -- automatic conversion of Eigen Tensor + + 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. +*/ + +#include + +#include "pybind11_tests.h" + +template +void reset_tensor(M &x) { + for (int i = 0; i < x.size(); i++) { + x(i) = i; + } +} + +Eigen::Tensor &get_tensor() { + static Eigen::Tensor *x; + + if (!x) { + x = new Eigen::Tensor(3, 1, 2); + reset_tensor(*x); + } + + return *x; +} + +Eigen::TensorFixedSize> &get_fixed_tensor() { + static Eigen::TensorFixedSize> *x; + + if (!x) { + Eigen::aligned_allocator>> allocator; + x = new (allocator.allocate(1)) Eigen::TensorFixedSize>(); + reset_tensor(*x); + } + + return *x; +} + +const Eigen::Tensor &get_const_tensor() { return get_tensor(); } + +TEST_SUBMODULE(eigen_tensor, m) { + m.def("copy_fixed_global_tensor", []() { return get_fixed_tensor(); }); + + m.def("copy_global_tensor", []() { return get_tensor(); }); + + m.def("copy_const_global_tensor", []() { return get_const_tensor(); }); + + m.def( + "reference_global_tensor", + []() { return &get_tensor(); }, + py::return_value_policy::reference); + + m.def( + "reference_const_global_tensor", + []() { return &get_const_tensor(); }, + py::return_value_policy::reference); + + m.def( + "reference_view_of_global_tensor", + []() { return Eigen::TensorMap>(get_tensor()); }, + py::return_value_policy::reference); + + m.def( + "reference_view_of_fixed_global_tensor", + []() { + return Eigen::TensorMap>>( + get_fixed_tensor()); + }, + py::return_value_policy::reference); + + m.def("round_trip_tensor", [](const py::array_t &foo) { + auto blah = foo.cast>(); + return blah; + }); + + m.def( + "round_trip_view_tensor", + [](const py::array_t &foo) { + auto view = foo.cast>>(); + return view; + }, + py::return_value_policy::reference); +} diff --git a/tests/test_eigen_tensor.py b/tests/test_eigen_tensor.py new file mode 100644 index 00000000000..932104b3d77 --- /dev/null +++ b/tests/test_eigen_tensor.py @@ -0,0 +1,64 @@ +import pytest + +from pybind11_tests import ConstructorStats + +np = pytest.importorskip("numpy") +m = pytest.importorskip("pybind11_tests.eigen_tensor") + +tensor_ref = np.array( + [ + [[0, 3]], + [[1, 4]], + [[2, 5]], + ] +) + + +def assert_equal_tensor_ref(mat, writeable=True, modified=0): + assert mat.flags.writeable == writeable + + if modified != 0: + tensor_ref[0, 0, 0] = modified + + np.testing.assert_array_equal(mat, tensor_ref) + + tensor_ref[0, 0, 0] = 0 + + +def test_convert_tensor_to_py(): + assert_equal_tensor_ref(m.copy_global_tensor()) + assert_equal_tensor_ref(m.copy_fixed_global_tensor()) + assert_equal_tensor_ref(m.copy_const_global_tensor()) + + assert_equal_tensor_ref(m.reference_global_tensor()) + assert_equal_tensor_ref(m.reference_view_of_global_tensor()) + assert_equal_tensor_ref(m.reference_view_of_fixed_global_tensor()) + assert_equal_tensor_ref(m.reference_const_global_tensor(), writeable=False) + + +def test_references_actually_refer(): + a = m.reference_global_tensor() + a[0, 0, 0] = 100 + assert_equal_tensor_ref(m.copy_const_global_tensor(), modified=100) + a[0, 0, 0] = 0 + assert_equal_tensor_ref(m.copy_const_global_tensor()) + + a = m.reference_view_of_global_tensor() + a[0, 0, 0] = 100 + assert_equal_tensor_ref(m.copy_const_global_tensor(), modified=100) + a[0, 0, 0] = 0 + assert_equal_tensor_ref(m.copy_const_global_tensor()) + + +def test_round_trip(): + assert_equal_tensor_ref(m.round_trip_tensor(tensor_ref)) + + +def test_round_trip_references_actually_refer(): + # Need to create a copy that matches the type on the C side + copy = np.array(tensor_ref, dtype=np.float64, order="F") + a = m.round_trip_view_tensor(copy) + a[0, 0, 0] = 100 + assert_equal_tensor_ref(copy, modified=100) + a[0, 0, 0] = 0 + assert_equal_tensor_ref(copy)