Skip to content

Commit

Permalink
Eigen: support sparse matrix types (#126)
Browse files Browse the repository at this point in the history
This commit adds support for sparse Eigen matrices along with tests.
  • Loading branch information
hmenke authored Mar 3, 2023
1 parent 67150d4 commit 8eadee2
Show file tree
Hide file tree
Showing 9 changed files with 314 additions and 18 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ jobs:
- name: Install NumPy
if: matrix.python != 'pypy3.9' && matrix.python != '3.12.0-alpha.5'
run: |
python -m pip install numpy
python -m pip install numpy scipy
- name: Configure
run: >
Expand Down
3 changes: 3 additions & 0 deletions docs/changelog.rst
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,9 @@ Version 0.2.0 (TBA)
* Added casters for dense matrix/array types from the `Eigen library
<https://eigen.tuxfamily.org/index.php?title=Main_Page>`__. (PR `#120
<https://github.com/wjakob/nanobind/pull/120>`__).
* Added casters for sparse matrix/array types from the `Eigen library
<https://eigen.tuxfamily.org/index.php?title=Main_Page>`__. (PR `#126
<https://github.com/wjakob/nanobind/pull/126>`_).
* Implemented `nb::bind_vector\<T\>() <bind_vector>` analogous to similar
functionality in pybind11. (commit `f2df8a
<https://github.com/wjakob/nanobind/commit/f2df8a90fbfb06ee03a79b0dd85fa0e266efeaa9>`__).
Expand Down
13 changes: 13 additions & 0 deletions docs/eigen.rst
Original file line number Diff line number Diff line change
Expand Up @@ -100,3 +100,16 @@ apply:
.. code-block:: cpp
void f4(nb::DRef<Eigen::MatrixXf> x) { x *= 2; }
Sparse matrices
---------------

There is also support for Eigen's sparse matrices which are mapped to
``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix``. To use it the extra header
:file:`nanobind/eigen/sparse.h` has to be included.

.. note::

There is no support for Eigen sparse vectors because there exists no
equivalent type in ``scipy.sparse``.
1 change: 1 addition & 0 deletions docs/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,7 @@ discourse:
ownership_adv
lowlevel
typeslots
eigen

.. toctree::
:caption: API Reference
Expand Down
39 changes: 23 additions & 16 deletions include/nanobind/eigen/dense.h
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@
#include <nanobind/ndarray.h>
#include <Eigen/Core>

static_assert(EIGEN_VERSION_AT_LEAST(3, 2, 7),
"Eigen matrix support in pybind11 requires Eigen >= 3.2.7");
static_assert(EIGEN_VERSION_AT_LEAST(3, 3, 1),
"Eigen matrix support in nanobind requires Eigen >= 3.3.1");

NAMESPACE_BEGIN(NB_NAMESPACE)

Expand All @@ -25,20 +25,23 @@ template <typename T> using DMap = Eigen::Map<T, 0, DStride>;

NAMESPACE_BEGIN(detail)

template <typename T>
constexpr int NumDimensions = int(T::MaxSizeAtCompileTime) == 1 ? 0 : bool(T::IsVectorAtCompileTime) ? 1 : 2;

template <typename T>
using array_for_eigen_t = ndarray<
typename T::Scalar,
numpy,
std::conditional_t<
T::NumDimensions == 1,
NumDimensions<T> == 1,
shape<(size_t) T::SizeAtCompileTime>,
shape<(size_t) T::RowsAtCompileTime,
(size_t) T::ColsAtCompileTime>>,
std::conditional_t<
T::InnerStrideAtCompileTime == Eigen::Dynamic,
any_contig,
std::conditional_t<
T::IsRowMajor || T::NumDimensions == 1,
T::IsRowMajor || NumDimensions<T> == 1,
c_contig,
f_contig
>
Expand All @@ -53,9 +56,13 @@ is_base_of_template_v<T, Eigen::EigenBase>;
template <typename T> constexpr bool is_eigen_plain_v =
is_base_of_template_v<T, Eigen::PlainObjectBase>;

/// Detect Eigen::SparseMatrix
template <typename T> constexpr bool is_eigen_sparse_v =
is_base_of_template_v<T, Eigen::SparseMatrixBase>;

/// Detects expression templates
template <typename T> constexpr bool is_eigen_xpr_v =
is_eigen_v<T> && !is_eigen_plain_v<T> &&
is_eigen_v<T> && !is_eigen_plain_v<T> && !is_eigen_sparse_v<T> &&
!std::is_base_of_v<Eigen::MapBase<T, Eigen::ReadOnlyAccessors>, T>;

template <typename T> struct type_caster<T, enable_if_t<is_eigen_plain_v<T>>> {
Expand All @@ -71,7 +78,7 @@ template <typename T> struct type_caster<T, enable_if_t<is_eigen_plain_v<T>>> {
return false;
const NDArray &array = caster.value;

if constexpr (T::NumDimensions == 1) {
if constexpr (NumDimensions<T> == 1) {
value.resize(array.shape(0));
memcpy(value.data(), array.data(),
array.shape(0) * sizeof(Scalar));
Expand All @@ -93,10 +100,10 @@ template <typename T> struct type_caster<T, enable_if_t<is_eigen_plain_v<T>>> {
}

static handle from_cpp(const T &v, rv_policy policy, cleanup_list *cleanup) noexcept {
size_t shape[T::NumDimensions];
int64_t strides[T::NumDimensions];
size_t shape[NumDimensions<T>];
int64_t strides[NumDimensions<T>];

if constexpr (T::NumDimensions == 1) {
if constexpr (NumDimensions<T> == 1) {
shape[0] = v.size();
strides[0] = v.innerStride();
} else {
Expand Down Expand Up @@ -139,7 +146,7 @@ template <typename T> struct type_caster<T, enable_if_t<is_eigen_plain_v<T>>> {
policy == rv_policy::move ? rv_policy::reference : policy;

object o = steal(NDArrayCaster::from_cpp(
NDArray(ptr, T::NumDimensions, shape, owner, strides),
NDArray(ptr, NumDimensions<T>, shape, owner, strides),
array_rv_policy, cleanup));

return o.release();
Expand All @@ -165,7 +172,7 @@ template <typename T> struct type_caster<T, enable_if_t<is_eigen_xpr_v<T>>> {

/// Caster for Eigen::Map<T>
template <typename T, int Options, typename StrideType>
struct type_caster<Eigen::Map<T, Options, StrideType>> {
struct type_caster<Eigen::Map<T, Options, StrideType>, enable_if_t<is_eigen_plain_v<T>>> {
using Map = Eigen::Map<T, Options, StrideType>;
using NDArray = array_for_eigen_t<Map>;
using NDArrayCaster = type_caster<NDArray>;
Expand All @@ -180,10 +187,10 @@ struct type_caster<Eigen::Map<T, Options, StrideType>> {
}

static handle from_cpp(const Map &v, rv_policy, cleanup_list *cleanup) noexcept {
size_t shape[T::NumDimensions];
int64_t strides[T::NumDimensions];
size_t shape[NumDimensions<T>];
int64_t strides[NumDimensions<T>];

if constexpr (T::NumDimensions == 1) {
if constexpr (NumDimensions<T> == 1) {
shape[0] = v.size();
strides[0] = v.innerStride();
} else {
Expand All @@ -194,7 +201,7 @@ struct type_caster<Eigen::Map<T, Options, StrideType>> {
}

return NDArrayCaster::from_cpp(
NDArray((void *) v.data(), T::NumDimensions, shape, handle(), strides),
NDArray((void *) v.data(), NumDimensions<T>, shape, handle(), strides),
rv_policy::reference, cleanup);
}

Expand Down Expand Up @@ -224,7 +231,7 @@ struct type_caster<Eigen::Map<T, Options, StrideType>> {

/// Caster for Eigen::Ref<T>
template <typename T, int Options, typename StrideType>
struct type_caster<Eigen::Ref<T, Options, StrideType>> {
struct type_caster<Eigen::Ref<T, Options, StrideType>, enable_if_t<is_eigen_plain_v<T>>> {
using Ref = Eigen::Ref<T, Options, StrideType>;
using Map = Eigen::Map<T, Options, StrideType>;
using MapCaster = make_caster<Map>;
Expand Down
177 changes: 177 additions & 0 deletions include/nanobind/eigen/sparse.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,177 @@
/*
nanobind/eigen/sparse.h: type casters for sparse Eigen matrices
Copyright (c) 2023 Henri Menke and Wenzel Jakob
All rights reserved. Use of this source code is governed by a
BSD-style license that can be found in the LICENSE file.
*/

#pragma once

#include <nanobind/ndarray.h>
#include <nanobind/eigen/dense.h>
#include <Eigen/SparseCore>

#include <memory>
#include <type_traits>
#include <utility>

NAMESPACE_BEGIN(NB_NAMESPACE)

NAMESPACE_BEGIN(detail)

/// Detect Eigen::SparseMatrix
template <typename T> constexpr bool is_eigen_sparse_matrix_v =
is_eigen_sparse_v<T> &&
!std::is_base_of_v<Eigen::SparseMapBase<T, Eigen::ReadOnlyAccessors>, T>;


/// Caster for Eigen::SparseMatrix
template <typename T> struct type_caster<T, enable_if_t<is_eigen_sparse_matrix_v<T>>> {
using Scalar = typename T::Scalar;
using StorageIndex = typename T::StorageIndex;
using Index = typename T::Index;
using SparseMap = Eigen::Map<T>;

static_assert(std::is_same_v<T, Eigen::SparseMatrix<Scalar, T::Options, StorageIndex>>,
"nanobind: Eigen sparse caster only implemented for matrices");

static constexpr bool row_major = T::IsRowMajor;

using ScalarNDArray = ndarray<numpy, Scalar, shape<any>>;
using StorageIndexNDArray = ndarray<numpy, StorageIndex, shape<any>>;

using ScalarCaster = make_caster<ScalarNDArray>;
using StorageIndexCaster = make_caster<StorageIndexNDArray>;

NB_TYPE_CASTER(T, const_name<row_major>("scipy.sparse.csr_matrix[",
"scipy.sparse.csc_matrix[")
+ make_caster<Scalar>::Name + const_name("]"));

ScalarCaster data_caster;
StorageIndexCaster indices_caster, indptr_caster;

bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept {
object obj = borrow(src);
try {
object matrix_type = module_::import_("scipy.sparse").attr(row_major ? "csr_matrix" : "csc_matrix");
if (!obj.type().is(matrix_type))
obj = matrix_type(obj);
} catch (const python_error &) {
return false;
}

if (object data_o = obj.attr("data"); !data_caster.from_python(data_o, flags, cleanup))
return false;
ScalarNDArray& values = data_caster.value;

if (object indices_o = obj.attr("indices"); !indices_caster.from_python(indices_o, flags, cleanup))
return false;
StorageIndexNDArray& inner_indices = indices_caster.value;

if (object indptr_o = obj.attr("indptr"); !indptr_caster.from_python(indptr_o, flags, cleanup))
return false;
StorageIndexNDArray& outer_indices = indptr_caster.value;

object shape_o = obj.attr("shape"), nnz_o = obj.attr("nnz");
Index rows, cols, nnz;
try {
if (len(shape_o) != 2)
return false;
rows = cast<Index>(shape_o[0]);
cols = cast<Index>(shape_o[1]);
nnz = cast<Index>(nnz_o);
} catch (const python_error &) {
return false;
}

value = SparseMap(rows, cols, nnz, outer_indices.data(), inner_indices.data(), values.data());

return true;
}

static handle from_cpp(T &&v, rv_policy policy, cleanup_list *cleanup) noexcept {
if (policy == rv_policy::automatic ||
policy == rv_policy::automatic_reference)
policy = rv_policy::move;

return from_cpp((const T &) v, policy, cleanup);
}

static handle from_cpp(const T &v, rv_policy policy, cleanup_list *) noexcept {
if (!v.isCompressed()) {
PyErr_SetString(PyExc_ValueError,
"nanobind: unable to return an Eigen sparse matrix that is not in a compressed format. "
"Please call `.makeCompressed()` before returning the value on the C++ end.");
return handle();
}

object matrix_type;
try {
matrix_type = module_::import_("scipy.sparse").attr(row_major ? "csr_matrix" : "csc_matrix");
} catch (python_error &e) {
e.restore();
return handle();
}

const Index rows = v.rows();
const Index cols = v.cols();
const size_t data_shape[] = { (size_t)v.nonZeros() };
const size_t outer_indices_shape[] = { (size_t)((row_major ? rows : cols) + 1) };

T *src = std::addressof(const_cast<T &>(v));
object owner;
if (policy == rv_policy::move) {
src = new T(std::move(v));
owner = capsule(src, [](void *p) noexcept { delete (T *) p; });
}

ScalarNDArray data(src->valuePtr(), 1, data_shape, owner);
StorageIndexNDArray outer_indices(src->outerIndexPtr(), 1, outer_indices_shape, owner);
StorageIndexNDArray inner_indices(src->innerIndexPtr(), 1, data_shape, owner);

try {
return matrix_type(make_tuple(
std::move(data), std::move(inner_indices), std::move(outer_indices)),
make_tuple(rows, cols))
.release();
} catch (python_error &e) {
e.restore();
return handle();
}
}
};


/// Caster for Eigen::Map<Eigen::SparseMatrix>
template <typename T>
struct type_caster<Eigen::Map<T>, enable_if_t<is_eigen_sparse_matrix_v<T>>> {
using Map = Eigen::Map<T>;
using SparseMatrixCaster = type_caster<T>;
static constexpr auto Name = SparseMatrixCaster::Name;
template <typename T_> using Cast = Map;

bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept = delete;

static handle from_cpp(const Map &v, rv_policy policy, cleanup_list *cleanup) noexcept = delete;
};


/// Caster for Eigen::Ref<Eigen::SparseMatrix>
template <typename T, int Options>
struct type_caster<Eigen::Ref<T, Options>, enable_if_t<is_eigen_sparse_matrix_v<T>>> {
using Ref = Eigen::Ref<T, Options>;
using Map = Eigen::Map<T, Options>;
using MapCaster = make_caster<Map>;
static constexpr auto Name = MapCaster::Name;
template <typename T_> using Cast = Ref;

bool from_python(handle src, uint8_t flags, cleanup_list *cleanup) noexcept = delete;

static handle from_cpp(const Ref &v, rv_policy policy, cleanup_list *cleanup) noexcept = delete;
};

NAMESPACE_END(detail)

NAMESPACE_END(NB_NAMESPACE)
2 changes: 1 addition & 1 deletion tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ nanobind_add_module(test_intrusive_ext test_intrusive.cpp object.cpp object.h ${
nanobind_add_module(test_exception_ext test_exception.cpp object.cpp object.h ${NB_EXTRA_ARGS})
nanobind_add_module(test_make_iterator_ext test_make_iterator.cpp ${NB_EXTRA_ARGS})

find_package (Eigen3 NO_MODULE)
find_package (Eigen3 3.3.1 NO_MODULE)
if (TARGET Eigen3::Eigen)
nanobind_add_module(test_eigen_ext test_eigen.cpp ${NB_EXTRA_ARGS})
target_link_libraries(test_eigen_ext PRIVATE Eigen3::Eigen)
Expand Down
24 changes: 24 additions & 0 deletions tests/test_eigen.cpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include <nanobind/eigen/dense.h>
#include <nanobind/eigen/sparse.h>

namespace nb = nanobind;

Expand Down Expand Up @@ -93,6 +94,29 @@ NB_MODULE(test_eigen_ext, m) {
m.def("updateV3i", [](Eigen::Ref<Eigen::Vector3i> a) { a[2] = 123; });
m.def("updateVXi", [](Eigen::Ref<Eigen::VectorXi> a) { a[2] = 123; });

using SparseMatrixR = Eigen::SparseMatrix<float, Eigen::RowMajor>;
using SparseMatrixC = Eigen::SparseMatrix<float>;
Eigen::MatrixXf mat(5, 6);
mat <<
0, 3, 0, 0, 0, 11,
22, 0, 0, 0, 17, 11,
7, 5, 0, 1, 0, 11,
0, 0, 0, 0, 0, 11,
0, 0, 14, 0, 8, 11;
m.def("sparse_r", [mat]() -> SparseMatrixR {
return Eigen::SparseView<Eigen::MatrixXf>(mat);
});
m.def("sparse_c", [mat]() -> SparseMatrixC {
return Eigen::SparseView<Eigen::MatrixXf>(mat);
});
m.def("sparse_copy_r", [](const SparseMatrixR &m) -> SparseMatrixR { return m; });
m.def("sparse_copy_c", [](const SparseMatrixC &m) -> SparseMatrixC { return m; });
m.def("sparse_r_uncompressed", []() -> SparseMatrixR {
SparseMatrixR m(2,2);
m.coeffRef(0,0) = 1.0f;
return m;
});

struct Buffer {
uint32_t x[30] { };

Expand Down
Loading

0 comments on commit 8eadee2

Please sign in to comment.