1 /*
2     Copyright (c) 2016 Wenzel Jakob <wenzel.jakob@epfl.ch>
3 
4     All rights reserved. Use of this source code is governed by a
5     BSD-style license that can be found in the LICENSE file.
6 */
7 
8 
9 
10 
11 #pragma once
12 
13 #include "numpy.h"
14 
15 #if defined(__INTEL_COMPILER)
16 #  pragma warning(disable: 1682)
17 #elif defined(__GNUG__) || defined(__clang__)
18 #  pragma GCC diagnostic push
19 #  pragma GCC diagnostic ignored "-Wconversion"
20 #  pragma GCC diagnostic ignored "-Wdeprecated-declarations"
21 #  ifdef __clang__
22 
23 
24 #    pragma GCC diagnostic ignored "-Wdeprecated"
25 #  endif
26 #  if __GNUC__ >= 7
27 #    pragma GCC diagnostic ignored "-Wint-in-bool-context"
28 #  endif
29 #endif
30 
31 #if defined(_MSC_VER)
32 #  pragma warning(push)
33 #  pragma warning(disable: 4127)
34 #  pragma warning(disable: 4996)
35 #endif
36 
37 #include <Eigen/Core>
38 #include <Eigen/SparseCore>
39 
40 
41 
42 
43 static_assert(EIGEN_VERSION_AT_LEAST(3,2,7), "Eigen support in pybind11 requires Eigen >= 3.2.7");
44 
45 NAMESPACE_BEGIN(PYBIND11_NAMESPACE)
46 
47 
48 using EigenDStride = Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>;
49 template <typename MatrixType> using EigenDRef = Eigen::Ref<MatrixType, 0, EigenDStride>;
50 template <typename MatrixType> using EigenDMap = Eigen::Map<MatrixType, 0, EigenDStride>;
51 
52 NAMESPACE_BEGIN(detail)
53 
54 #if EIGEN_VERSION_AT_LEAST(3,3,0)
55 using EigenIndex = Eigen::Index;
56 #else
57 using EigenIndex = EIGEN_DEFAULT_DENSE_INDEX_TYPE;
58 #endif
59 
60 
61 template <typename T> using is_eigen_dense_map = all_of<is_template_base_of<Eigen::DenseBase, T>, std::is_base_of<Eigen::MapBase<T, Eigen::ReadOnlyAccessors>, T>>;
62 template <typename T> using is_eigen_mutable_map = std::is_base_of<Eigen::MapBase<T, Eigen::WriteAccessors>, T>;
63 template <typename T> using is_eigen_dense_plain = all_of<negation<is_eigen_dense_map<T>>, is_template_base_of<Eigen::PlainObjectBase, T>>;
64 template <typename T> using is_eigen_sparse = is_template_base_of<Eigen::SparseMatrixBase, T>;
65 
66 
67 
68 
69 template <typename T> using is_eigen_other = all_of<
70     is_template_base_of<Eigen::EigenBase, T>,
71     negation<any_of<is_eigen_dense_map<T>, is_eigen_dense_plain<T>, is_eigen_sparse<T>>>
72 >;
73 
74 
75 template <bool EigenRowMajor> struct EigenConformable {
76     bool conformable = false;
77     EigenIndex rows = 0, cols = 0;
78     EigenDStride stride{0, 0};
79     bool negativestrides = false;
80 
81     EigenConformable(bool fits = false) : conformable{fits} {}
82 
EigenConformableEigenConformable83     EigenConformable(EigenIndex r, EigenIndex c,
84             EigenIndex rstride, EigenIndex cstride) :
85         conformable{true}, rows{r}, cols{c} {
86 
87         if (rstride < 0 || cstride < 0) {
88             negativestrides = true;
89         } else {
90             stride = {EigenRowMajor ? rstride : cstride  ,
91                       EigenRowMajor ? cstride : rstride   };
92         }
93     }
94 
EigenConformableEigenConformable95     EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride)
96         : EigenConformable(r, c, r == 1 ? c*stride : stride, c == 1 ? r : r*stride) {}
97 
stride_compatibleEigenConformable98     template <typename props> bool stride_compatible() const {
99 
100 
101         return
102             !negativestrides &&
103             (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner() ||
104                 (EigenRowMajor ? cols : rows) == 1) &&
105             (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer() ||
106                 (EigenRowMajor ? rows : cols) == 1);
107     }
108     operator bool() const { return conformable; }
109 };
110 
111 template <typename Type> struct eigen_extract_stride { using type = Type; };
112 template <typename PlainObjectType, int MapOptions, typename StrideType>
113 struct eigen_extract_stride<Eigen::Map<PlainObjectType, MapOptions, StrideType>> { using type = StrideType; };
114 template <typename PlainObjectType, int Options, typename StrideType>
115 struct eigen_extract_stride<Eigen::Ref<PlainObjectType, Options, StrideType>> { using type = StrideType; };
116 
117 
118 template <typename Type_> struct EigenProps {
119     using Type = Type_;
120     using Scalar = typename Type::Scalar;
121     using StrideType = typename eigen_extract_stride<Type>::type;
122     static constexpr EigenIndex
123         rows = Type::RowsAtCompileTime,
124         cols = Type::ColsAtCompileTime,
125         size = Type::SizeAtCompileTime;
126     static constexpr bool
127         row_major = Type::IsRowMajor,
128         vector = Type::IsVectorAtCompileTime,
129         fixed_rows = rows != Eigen::Dynamic,
130         fixed_cols = cols != Eigen::Dynamic,
131         fixed = size != Eigen::Dynamic,
132         dynamic = !fixed_rows && !fixed_cols;
133 
134     template <EigenIndex i, EigenIndex ifzero> using if_zero = std::integral_constant<EigenIndex, i == 0 ? ifzero : i>;
135     static constexpr EigenIndex inner_stride = if_zero<StrideType::InnerStrideAtCompileTime, 1>::value,
136                                 outer_stride = if_zero<StrideType::OuterStrideAtCompileTime,
137                                                        vector ? size : row_major ? cols : rows>::value;
138     static constexpr bool dynamic_stride = inner_stride == Eigen::Dynamic && outer_stride == Eigen::Dynamic;
139     static constexpr bool requires_row_major = !dynamic_stride && !vector && (row_major ? inner_stride : outer_stride) == 1;
140     static constexpr bool requires_col_major = !dynamic_stride && !vector && (row_major ? outer_stride : inner_stride) == 1;
141 
142 
143 
144 
145     static EigenConformable<row_major> conformable(const array &a) {
146         const auto dims = a.ndim();
147         if (dims < 1 || dims > 2)
148             return false;
149 
150         if (dims == 2) {
151 
152             EigenIndex
153                 np_rows = a.shape(0),
154                 np_cols = a.shape(1),
155                 np_rstride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar)),
156                 np_cstride = a.strides(1) / static_cast<ssize_t>(sizeof(Scalar));
157             if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols))
158                 return false;
159 
160             return {np_rows, np_cols, np_rstride, np_cstride};
161         }
162 
163 
164 
165         const EigenIndex n = a.shape(0),
166               stride = a.strides(0) / static_cast<ssize_t>(sizeof(Scalar));
167 
168         if (vector) {
169             if (fixed && size != n)
170                 return false;
171             return {rows == 1 ? 1 : n, cols == 1 ? 1 : n, stride};
172         }
173         else if (fixed) {
174 
175             return false;
176         }
177         else if (fixed_cols) {
178 
179 
180             if (cols != n) return false;
181             return {1, n, stride};
182         }
183         else {
184 
185             if (fixed_rows && rows != n) return false;
186             return {n, 1, stride};
187         }
188     }
189 
190     static constexpr bool show_writeable = is_eigen_dense_map<Type>::value && is_eigen_mutable_map<Type>::value;
191     static constexpr bool show_order = is_eigen_dense_map<Type>::value;
192     static constexpr bool show_c_contiguous = show_order && requires_row_major;
193     static constexpr bool show_f_contiguous = !show_c_contiguous && show_order && requires_col_major;
194 
195     static constexpr auto descriptor =
196         _("numpy.ndarray[") + npy_format_descriptor<Scalar>::name +
197         _("[")  + _<fixed_rows>(_<(size_t) rows>(), _("m")) +
198         _(", ") + _<fixed_cols>(_<(size_t) cols>(), _("n")) +
199         _("]") +
200 
201 
202 
203 
204 
205 
206         _<show_writeable>(", flags.writeable", "") +
207         _<show_c_contiguous>(", flags.c_contiguous", "") +
208         _<show_f_contiguous>(", flags.f_contiguous", "") +
209         _("]");
210 };
211 
212 
213 
214 template <typename props> handle eigen_array_cast(typename props::Type const &src, handle base = handle(), bool writeable = true) {
215     constexpr ssize_t elem_size = sizeof(typename props::Scalar);
216     array a;
217     if (props::vector)
218         a = array({ src.size() }, { elem_size * src.innerStride() }, src.data(), base);
219     else
220         a = array({ src.rows(), src.cols() }, { elem_size * src.rowStride(), elem_size * src.colStride() },
221                   src.data(), base);
222 
223     if (!writeable)
224         array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
225 
226     return a.release();
227 }
228 
229 
230 
231 
232 
233 template <typename props, typename Type>
234 handle eigen_ref_array(Type &src, handle parent = none()) {
235 
236 
237     return eigen_array_cast<props>(src, parent, !std::is_const<Type>::value);
238 }
239 
240 
241 
242 
243 
244 template <typename props, typename Type, typename = enable_if_t<is_eigen_dense_plain<Type>::value>>
245 handle eigen_encapsulate(Type *src) {
246     capsule base(src, [](void *o) { delete static_cast<Type *>(o); });
247     return eigen_ref_array<props>(*src, base);
248 }
249 
250 
251 
252 template<typename Type>
253 struct type_caster<Type, enable_if_t<is_eigen_dense_plain<Type>::value>> {
254     using Scalar = typename Type::Scalar;
255     using props = EigenProps<Type>;
256 
257     bool load(handle src, bool convert) {
258 
259         if (!convert && !isinstance<array_t<Scalar>>(src))
260             return false;
261 
262 
263         auto buf = array::ensure(src);
264 
265         if (!buf)
266             return false;
267 
268         auto dims = buf.ndim();
269         if (dims < 1 || dims > 2)
270             return false;
271 
272         auto fits = props::conformable(buf);
273         if (!fits)
274             return false;
275 
276 
277         value = Type(fits.rows, fits.cols);
278         auto ref = reinterpret_steal<array>(eigen_ref_array<props>(value));
279         if (dims == 1) ref = ref.squeeze();
280         else if (ref.ndim() == 1) buf = buf.squeeze();
281 
282         int result = detail::npy_api::get().PyArray_CopyInto_(ref.ptr(), buf.ptr());
283 
284         if (result < 0) {
285             PyErr_Clear();
286             return false;
287         }
288 
289         return true;
290     }
291 
292 private:
293 
294 
295     template <typename CType>
296     static handle cast_impl(CType *src, return_value_policy policy, handle parent) {
297         switch (policy) {
298             case return_value_policy::take_ownership:
299             case return_value_policy::automatic:
300                 return eigen_encapsulate<props>(src);
301             case return_value_policy::move:
302                 return eigen_encapsulate<props>(new CType(std::move(*src)));
303             case return_value_policy::copy:
304                 return eigen_array_cast<props>(*src);
305             case return_value_policy::reference:
306             case return_value_policy::automatic_reference:
307                 return eigen_ref_array<props>(*src);
308             case return_value_policy::reference_internal:
309                 return eigen_ref_array<props>(*src, parent);
310             default:
311                 throw cast_error("unhandled return_value_policy: should not happen!");
312         };
313     }
314 
315 public:
316 
317 
318     static handle cast(Type &&src, return_value_policy  , handle parent) {
319         return cast_impl(&src, return_value_policy::move, parent);
320     }
321 
322     static handle cast(const Type &&src, return_value_policy  , handle parent) {
323         return cast_impl(&src, return_value_policy::move, parent);
324     }
325 
326     static handle cast(Type &src, return_value_policy policy, handle parent) {
327         if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference)
328             policy = return_value_policy::copy;
329         return cast_impl(&src, policy, parent);
330     }
331 
332     static handle cast(const Type &src, return_value_policy policy, handle parent) {
333         if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference)
334             policy = return_value_policy::copy;
335         return cast(&src, policy, parent);
336     }
337 
338     static handle cast(Type *src, return_value_policy policy, handle parent) {
339         return cast_impl(src, policy, parent);
340     }
341 
342     static handle cast(const Type *src, return_value_policy policy, handle parent) {
343         return cast_impl(src, policy, parent);
344     }
345 
346     static constexpr auto name = props::descriptor;
347 
348     operator Type*() { return &value; }
349     operator Type&() { return value; }
350     operator Type&&() && { return std::move(value); }
351     template <typename T> using cast_op_type = movable_cast_op_type<T>;
352 
353 private:
354     Type value;
355 };
356 
357 
358 template <typename MapType> struct eigen_map_caster {
359 private:
360     using props = EigenProps<MapType>;
361 
362 public:
363 
364 
365 
366 
367 
368 
369 
370     static handle cast(const MapType &src, return_value_policy policy, handle parent) {
371         switch (policy) {
372             case return_value_policy::copy:
373                 return eigen_array_cast<props>(src);
374             case return_value_policy::reference_internal:
375                 return eigen_array_cast<props>(src, parent, is_eigen_mutable_map<MapType>::value);
376             case return_value_policy::reference:
377             case return_value_policy::automatic:
378             case return_value_policy::automatic_reference:
379                 return eigen_array_cast<props>(src, none(), is_eigen_mutable_map<MapType>::value);
380             default:
381 
382                 pybind11_fail("Invalid return_value_policy for Eigen Map/Ref/Block type");
383         }
384     }
385 
386     static constexpr auto name = props::descriptor;
387 
388 
389 
390 
391     bool load(handle, bool) = delete;
392     operator MapType() = delete;
393     template <typename> using cast_op_type = MapType;
394 };
395 
396 
397 template <typename Type> struct type_caster<Type, enable_if_t<is_eigen_dense_map<Type>::value>>
398     : eigen_map_caster<Type> {};
399 
400 
401 
402 template <typename PlainObjectType, typename StrideType>
403 struct type_caster<
404     Eigen::Ref<PlainObjectType, 0, StrideType>,
405     enable_if_t<is_eigen_dense_map<Eigen::Ref<PlainObjectType, 0, StrideType>>::value>
406 > : public eigen_map_caster<Eigen::Ref<PlainObjectType, 0, StrideType>> {
407 private:
408     using Type = Eigen::Ref<PlainObjectType, 0, StrideType>;
409     using props = EigenProps<Type>;
410     using Scalar = typename props::Scalar;
411     using MapType = Eigen::Map<PlainObjectType, 0, StrideType>;
412     using Array = array_t<Scalar, array::forcecast |
413                 ((props::row_major ? props::inner_stride : props::outer_stride) == 1 ? array::c_style :
414                  (props::row_major ? props::outer_stride : props::inner_stride) == 1 ? array::f_style : 0)>;
415     static constexpr bool need_writeable = is_eigen_mutable_map<Type>::value;
416 
417     std::unique_ptr<MapType> map;
418     std::unique_ptr<Type> ref;
419 
420 
421 
422 
423 
424 
425     Array copy_or_ref;
426 public:
427     bool load(handle src, bool convert) {
428 
429 
430         bool need_copy = !isinstance<Array>(src);
431 
432         EigenConformable<props::row_major> fits;
433         if (!need_copy) {
434 
435 
436             Array aref = reinterpret_borrow<Array>(src);
437 
438             if (aref && (!need_writeable || aref.writeable())) {
439                 fits = props::conformable(aref);
440                 if (!fits) return false;
441                 if (!fits.template stride_compatible<props>())
442                     need_copy = true;
443                 else
444                     copy_or_ref = std::move(aref);
445             }
446             else {
447                 need_copy = true;
448             }
449         }
450 
451         if (need_copy) {
452 
453 
454 
455             if (!convert || need_writeable) return false;
456 
457             Array copy = Array::ensure(src);
458             if (!copy) return false;
459             fits = props::conformable(copy);
460             if (!fits || !fits.template stride_compatible<props>())
461                 return false;
462             copy_or_ref = std::move(copy);
463             loader_life_support::add_patient(copy_or_ref);
464         }
465 
466         ref.reset();
467         map.reset(new MapType(data(copy_or_ref), fits.rows, fits.cols, make_stride(fits.stride.outer(), fits.stride.inner())));
468         ref.reset(new Type(*map));
469 
470         return true;
471     }
472 
473     operator Type*() { return ref.get(); }
474     operator Type&() { return *ref; }
475     template <typename _T> using cast_op_type = pybind11::detail::cast_op_type<_T>;
476 
477 private:
478     template <typename T = Type, enable_if_t<is_eigen_mutable_map<T>::value, int> = 0>
479     Scalar *data(Array &a) { return a.mutable_data(); }
480 
481     template <typename T = Type, enable_if_t<!is_eigen_mutable_map<T>::value, int> = 0>
482     const Scalar *data(Array &a) { return a.data(); }
483 
484 
485 
486     template <typename S> using stride_ctor_default = bool_constant<
487         S::InnerStrideAtCompileTime != Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic &&
488         std::is_default_constructible<S>::value>;
489 
490 
491     template <typename S> using stride_ctor_dual = bool_constant<
492         !stride_ctor_default<S>::value && std::is_constructible<S, EigenIndex, EigenIndex>::value>;
493 
494 
495     template <typename S> using stride_ctor_outer = bool_constant<
496         !any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value &&
497         S::OuterStrideAtCompileTime == Eigen::Dynamic && S::InnerStrideAtCompileTime != Eigen::Dynamic &&
498         std::is_constructible<S, EigenIndex>::value>;
499     template <typename S> using stride_ctor_inner = bool_constant<
500         !any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value &&
501         S::InnerStrideAtCompileTime == Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic &&
502         std::is_constructible<S, EigenIndex>::value>;
503 
504     template <typename S = StrideType, enable_if_t<stride_ctor_default<S>::value, int> = 0>
505     static S make_stride(EigenIndex, EigenIndex) { return S(); }
506     template <typename S = StrideType, enable_if_t<stride_ctor_dual<S>::value, int> = 0>
507     static S make_stride(EigenIndex outer, EigenIndex inner) { return S(outer, inner); }
508     template <typename S = StrideType, enable_if_t<stride_ctor_outer<S>::value, int> = 0>
509     static S make_stride(EigenIndex outer, EigenIndex) { return S(outer); }
510     template <typename S = StrideType, enable_if_t<stride_ctor_inner<S>::value, int> = 0>
511     static S make_stride(EigenIndex, EigenIndex inner) { return S(inner); }
512 
513 };
514 
515 
516 
517 
518 
519 template <typename Type>
520 struct type_caster<Type, enable_if_t<is_eigen_other<Type>::value>> {
521 protected:
522     using Matrix = Eigen::Matrix<typename Type::Scalar, Type::RowsAtCompileTime, Type::ColsAtCompileTime>;
523     using props = EigenProps<Matrix>;
524 public:
525     static handle cast(const Type &src, return_value_policy  , handle  ) {
526         handle h = eigen_encapsulate<props>(new Matrix(src));
527         return h;
528     }
529     static handle cast(const Type *src, return_value_policy policy, handle parent) { return cast(*src, policy, parent); }
530 
531     static constexpr auto name = props::descriptor;
532 
533 
534 
535 
536     bool load(handle, bool) = delete;
537     operator Type() = delete;
538     template <typename> using cast_op_type = Type;
539 };
540 
541 template<typename Type>
542 struct type_caster<Type, enable_if_t<is_eigen_sparse<Type>::value>> {
543     typedef typename Type::Scalar Scalar;
544     typedef remove_reference_t<decltype(*std::declval<Type>().outerIndexPtr())> StorageIndex;
545     typedef typename Type::Index Index;
546     static constexpr bool rowMajor = Type::IsRowMajor;
547 
548     bool load(handle src, bool) {
549         if (!src)
550             return false;
551 
552         auto obj = reinterpret_borrow<object>(src);
553         object sparse_module = module::import("scipy.sparse");
554         object matrix_type = sparse_module.attr(
555             rowMajor ? "csr_matrix" : "csc_matrix");
556 
557         if (!obj.get_type().is(matrix_type)) {
558             try {
559                 obj = matrix_type(obj);
560             } catch (const error_already_set &) {
561                 return false;
562             }
563         }
564 
565         auto values = array_t<Scalar>((object) obj.attr("data"));
566         auto innerIndices = array_t<StorageIndex>((object) obj.attr("indices"));
567         auto outerIndices = array_t<StorageIndex>((object) obj.attr("indptr"));
568         auto shape = pybind11::tuple((pybind11::object) obj.attr("shape"));
569         auto nnz = obj.attr("nnz").cast<Index>();
570 
571         if (!values || !innerIndices || !outerIndices)
572             return false;
573 
574         value = Eigen::MappedSparseMatrix<Scalar, Type::Flags, StorageIndex>(
575             shape[0].cast<Index>(), shape[1].cast<Index>(), nnz,
576             outerIndices.mutable_data(), innerIndices.mutable_data(), values.mutable_data());
577 
578         return true;
579     }
580 
581     static handle cast(const Type &src, return_value_policy  , handle  ) {
582         const_cast<Type&>(src).makeCompressed();
583 
584         object matrix_type = module::import("scipy.sparse").attr(
585             rowMajor ? "csr_matrix" : "csc_matrix");
586 
587         array data(src.nonZeros(), src.valuePtr());
588         array outerIndices((rowMajor ? src.rows() : src.cols()) + 1, src.outerIndexPtr());
589         array innerIndices(src.nonZeros(), src.innerIndexPtr());
590 
591         return matrix_type(
592             std::make_tuple(data, innerIndices, outerIndices),
593             std::make_pair(src.rows(), src.cols())
594         ).release();
595     }
596 
597     PYBIND11_TYPE_CASTER(Type, _<(Type::IsRowMajor) != 0>("scipy.sparse.csr_matrix[", "scipy.sparse.csc_matrix[")
598             + npy_format_descriptor<Scalar>::name + _("]"));
599 };
600 
601 NAMESPACE_END(detail)
602 NAMESPACE_END(PYBIND11_NAMESPACE)
603 
604 #if defined(__GNUG__) || defined(__clang__)
605 #  pragma GCC diagnostic pop
606 #elif defined(_MSC_VER)
607 #  pragma warning(pop)
608 #endif
609