Skip to content

Commit

Permalink
sparse: add first working version for Eigen::SparseMatrix
Browse files Browse the repository at this point in the history
  • Loading branch information
jcarpent committed Jan 25, 2024
1 parent 64f012a commit a2484d7
Show file tree
Hide file tree
Showing 6 changed files with 566 additions and 2 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,8 @@ set(${PROJECT_NAME}_HEADERS
include/eigenpy/pickle-vector.hpp
include/eigenpy/stride.hpp
include/eigenpy/tensor/eigen-from-python.hpp
include/eigenpy/sparse/eigen-from-python.hpp
include/eigenpy/scipy-allocator.hpp
include/eigenpy/scipy-type.hpp
include/eigenpy/swig.hpp
include/eigenpy/version.hpp)
Expand Down
2 changes: 2 additions & 0 deletions include/eigenpy/eigen-from-python.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -565,4 +565,6 @@ struct EigenFromPy<const Eigen::Ref<const MatType, Options, Stride> > {
#include "eigenpy/tensor/eigen-from-python.hpp"
#endif

#include "eigenpy/sparse/eigen-from-python.hpp"

#endif // __eigenpy_eigen_from_python_hpp__
51 changes: 49 additions & 2 deletions include/eigenpy/eigen-to-python.hpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
//
// Copyright (c) 2014-2023 CNRS INRIA
// Copyright (c) 2014-2024 CNRS INRIA
//

#ifndef __eigenpy_eigen_to_python_hpp__
Expand All @@ -11,7 +11,9 @@

#include "eigenpy/eigen-allocator.hpp"
#include "eigenpy/numpy-allocator.hpp"
#include "eigenpy/scipy-allocator.hpp"
#include "eigenpy/numpy-type.hpp"
#include "eigenpy/scipy-type.hpp"
#include "eigenpy/registration.hpp"

namespace boost {
Expand Down Expand Up @@ -116,6 +118,50 @@ struct eigen_to_py_impl_matrix {
// Create an instance (either np.array or np.matrix)
return NumpyType::make(pyArray).ptr();
}

static PyTypeObject const* get_pytype() { return getPyArrayType(); }
};

template <typename MatType>
struct eigen_to_py_impl_sparse_matrix;

template <typename MatType>
struct eigen_to_py_impl<MatType, Eigen::SparseMatrixBase<MatType> >
: eigen_to_py_impl_sparse_matrix<MatType> {};

template <typename MatType>
struct eigen_to_py_impl<MatType&, Eigen::SparseMatrixBase<MatType> >
: eigen_to_py_impl_sparse_matrix<MatType&> {};

template <typename MatType>
struct eigen_to_py_impl<const MatType, const Eigen::SparseMatrixBase<MatType> >
: eigen_to_py_impl_sparse_matrix<const MatType> {};

template <typename MatType>
struct eigen_to_py_impl<const MatType&, const Eigen::SparseMatrixBase<MatType> >
: eigen_to_py_impl_sparse_matrix<const MatType&> {};

template <typename MatType>
struct eigen_to_py_impl_sparse_matrix {
enum { IsRowMajor = MatType::IsRowMajor };

static PyObject* convert(
typename boost::add_reference<
typename boost::add_const<MatType>::type>::type mat) {
typedef typename boost::remove_const<
typename boost::remove_reference<MatType>::type>::type MatrixDerived;

// Allocate and perform the copy
PyObject* pyArray =
ScipyAllocator<MatType>::allocate(const_cast<MatrixDerived&>(mat));

return pyArray;
}

static PyTypeObject const* get_pytype() {
return IsRowMajor ? ScipyType::getScipyCSRMatrixType()
: ScipyType::getScipyCSCMatrixType();
}
};

#ifdef EIGENPY_WITH_TENSOR_SUPPORT
Expand Down Expand Up @@ -149,6 +195,8 @@ struct eigen_to_py_impl_tensor {
// Create an instance (either np.array or np.matrix)
return NumpyType::make(pyArray).ptr();
}

static PyTypeObject const* get_pytype() { return getPyArrayType(); }
};
#endif

Expand All @@ -163,7 +211,6 @@ template <typename EigenType, typename _Scalar>
struct EigenToPy
#endif
: eigen_to_py_impl<EigenType> {
static PyTypeObject const* get_pytype() { return getPyArrayType(); }
};

template <typename MatType>
Expand Down
225 changes: 225 additions & 0 deletions include/eigenpy/scipy-allocator.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,225 @@
/*
* Copyright 2024 INRIA
*/

#ifndef __eigenpy_scipy_allocator_hpp__
#define __eigenpy_scipy_allocator_hpp__

#include "eigenpy/fwd.hpp"
#include "eigenpy/eigen-allocator.hpp"
#include "eigenpy/scipy-type.hpp"
#include "eigenpy/register.hpp"

#include <iostream>

namespace eigenpy {

template <typename EigenType, typename BaseType>
struct scipy_allocator_impl;

template <typename EigenType>
struct scipy_allocator_impl_sparse_matrix;

template <typename MatType>
struct scipy_allocator_impl<
MatType,
Eigen::SparseMatrixBase<typename remove_const_reference<MatType>::type> >
: scipy_allocator_impl_sparse_matrix<MatType> {};

template <typename MatType>
struct scipy_allocator_impl<
const MatType, const Eigen::SparseMatrixBase<
typename remove_const_reference<MatType>::type> >
: scipy_allocator_impl_sparse_matrix<const MatType> {};

// template <typename MatType>
// struct scipy_allocator_impl<MatType &, Eigen::MatrixBase<MatType> > :
// scipy_allocator_impl_sparse_matrix<MatType &>
//{};

template <typename MatType>
struct scipy_allocator_impl<const MatType &,
const Eigen::SparseMatrixBase<MatType> >
: scipy_allocator_impl_sparse_matrix<const MatType &> {};

template <typename EigenType,
typename BaseType = typename get_eigen_base_type<EigenType>::type>
struct ScipyAllocator : scipy_allocator_impl<EigenType, BaseType> {};

template <typename MatType>
struct scipy_allocator_impl_sparse_matrix {
template <typename SimilarMatrixType>
static PyObject *allocate(
const Eigen::SparseCompressedBase<SimilarMatrixType> &mat_,
bool copy = false) {
EIGENPY_UNUSED_VARIABLE(copy);
typedef typename SimilarMatrixType::Scalar Scalar;
typedef typename SimilarMatrixType::StorageIndex StorageIndex;

enum { IsRowMajor = SimilarMatrixType::IsRowMajor };

typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> DataVector;
typedef const Eigen::Map<const DataVector> MapDataVector;
typedef Eigen::Matrix<StorageIndex, Eigen::Dynamic, 1> StorageIndexVector;
typedef Eigen::Matrix<int32_t, Eigen::Dynamic, 1> ScipyStorageIndexVector;
typedef const Eigen::Map<const StorageIndexVector> MapStorageIndexVector;

SimilarMatrixType &mat = mat_.const_cast_derived();
bp::object scipy_sparse_matrix_type =
ScipyType::get_pytype_object<SimilarMatrixType>();

MapDataVector data(mat.valuePtr(), mat.nonZeros());
MapStorageIndexVector outer_indices(
mat.outerIndexPtr(), (IsRowMajor ? mat.rows() : mat.cols()) + 1);
MapStorageIndexVector inner_indices(mat.innerIndexPtr(), mat.nonZeros());

bp::object scipy_sparse_matrix = scipy_sparse_matrix_type(bp::make_tuple(
DataVector(data),
ScipyStorageIndexVector(inner_indices.template cast<int32_t>()),
ScipyStorageIndexVector(outer_indices.template cast<int32_t>()))); //,
// bp::make_tuple(mat.rows(),
// mat.cols())));

Py_INCREF(scipy_sparse_matrix.ptr());
return scipy_sparse_matrix.ptr();
}
};

// template <typename MatType>
// struct scipy_allocator_impl_sparse_matrix<MatType &> {
// template <typename SimilarMatrixType>
// static PyArrayObject *allocate(Eigen::PlainObjectBase<SimilarMatrixType>
// &mat,
// npy_intp nd, npy_intp *shape) {
// typedef typename SimilarMatrixType::Scalar Scalar;
// enum {
// NPY_ARRAY_MEMORY_CONTIGUOUS =
// SimilarMatrixType::IsRowMajor ? NPY_ARRAY_CARRAY : NPY_ARRAY_FARRAY
// };
//
// if (NumpyType::sharedMemory()) {
// const int Scalar_type_code = Register::getTypeCode<Scalar>();
// PyArrayObject *pyArray = (PyArrayObject *)call_PyArray_New(
// getPyArrayType(), static_cast<int>(nd), shape, Scalar_type_code,
// mat.data(), NPY_ARRAY_MEMORY_CONTIGUOUS | NPY_ARRAY_ALIGNED);
//
// return pyArray;
// } else {
// return NumpyAllocator<MatType>::allocate(mat, nd, shape);
// }
// }
// };

#if EIGEN_VERSION_AT_LEAST(3, 2, 0)

// template <typename MatType, int Options, typename Stride>
// struct scipy_allocator_impl_sparse_matrix<Eigen::Ref<MatType, Options,
// Stride> > {
// typedef Eigen::Ref<MatType, Options, Stride> RefType;
//
// static PyArrayObject *allocate(RefType &mat, npy_intp nd, npy_intp *shape)
// {
// typedef typename RefType::Scalar Scalar;
// enum {
// NPY_ARRAY_MEMORY_CONTIGUOUS =
// RefType::IsRowMajor ? NPY_ARRAY_CARRAY : NPY_ARRAY_FARRAY
// };
//
// if (NumpyType::sharedMemory()) {
// const int Scalar_type_code = Register::getTypeCode<Scalar>();
// const bool reverse_strides = MatType::IsRowMajor || (mat.rows() == 1);
// Eigen::DenseIndex inner_stride = reverse_strides ? mat.outerStride()
// : mat.innerStride(),
// outer_stride = reverse_strides ? mat.innerStride()
// : mat.outerStride();
//
// const int elsize =
// call_PyArray_DescrFromType(Scalar_type_code)->elsize; npy_intp
// strides[2] = {elsize * inner_stride, elsize * outer_stride};
//
// PyArrayObject *pyArray = (PyArrayObject *)call_PyArray_New(
// getPyArrayType(), static_cast<int>(nd), shape, Scalar_type_code,
// strides, mat.data(), NPY_ARRAY_MEMORY_CONTIGUOUS |
// NPY_ARRAY_ALIGNED);
//
// return pyArray;
// } else {
// return NumpyAllocator<MatType>::allocate(mat, nd, shape);
// }
// }
// };

#endif

// template <typename MatType>
// struct scipy_allocator_impl_sparse_matrix<const MatType &> {
// template <typename SimilarMatrixType>
// static PyArrayObject *allocate(
// const Eigen::PlainObjectBase<SimilarMatrixType> &mat, npy_intp nd,
// npy_intp *shape) {
// typedef typename SimilarMatrixType::Scalar Scalar;
// enum {
// NPY_ARRAY_MEMORY_CONTIGUOUS_RO = SimilarMatrixType::IsRowMajor
// ? NPY_ARRAY_CARRAY_RO
// : NPY_ARRAY_FARRAY_RO
// };
//
// if (NumpyType::sharedMemory()) {
// const int Scalar_type_code = Register::getTypeCode<Scalar>();
// PyArrayObject *pyArray = (PyArrayObject *)call_PyArray_New(
// getPyArrayType(), static_cast<int>(nd), shape, Scalar_type_code,
// const_cast<Scalar *>(mat.data()),
// NPY_ARRAY_MEMORY_CONTIGUOUS_RO | NPY_ARRAY_ALIGNED);
//
// return pyArray;
// } else {
// return NumpyAllocator<MatType>::allocate(mat, nd, shape);
// }
// }
// };

#if EIGEN_VERSION_AT_LEAST(3, 2, 0)

// template <typename MatType, int Options, typename Stride>
// struct scipy_allocator_impl_sparse_matrix<
// const Eigen::Ref<const MatType, Options, Stride> > {
// typedef const Eigen::Ref<const MatType, Options, Stride> RefType;
//
// static PyArrayObject *allocate(RefType &mat, npy_intp nd, npy_intp *shape)
// {
// typedef typename RefType::Scalar Scalar;
// enum {
// NPY_ARRAY_MEMORY_CONTIGUOUS_RO =
// RefType::IsRowMajor ? NPY_ARRAY_CARRAY_RO : NPY_ARRAY_FARRAY_RO
// };
//
// if (NumpyType::sharedMemory()) {
// const int Scalar_type_code = Register::getTypeCode<Scalar>();
//
// const bool reverse_strides = MatType::IsRowMajor || (mat.rows() == 1);
// Eigen::DenseIndex inner_stride = reverse_strides ? mat.outerStride()
// : mat.innerStride(),
// outer_stride = reverse_strides ? mat.innerStride()
// : mat.outerStride();
//
// const int elsize =
// call_PyArray_DescrFromType(Scalar_type_code)->elsize; npy_intp
// strides[2] = {elsize * inner_stride, elsize * outer_stride};
//
// PyArrayObject *pyArray = (PyArrayObject *)call_PyArray_New(
// getPyArrayType(), static_cast<int>(nd), shape, Scalar_type_code,
// strides, const_cast<Scalar *>(mat.data()),
// NPY_ARRAY_MEMORY_CONTIGUOUS_RO | NPY_ARRAY_ALIGNED);
//
// return pyArray;
// } else {
// return NumpyAllocator<MatType>::allocate(mat, nd, shape);
// }
// }
// };

#endif

} // namespace eigenpy

#endif // ifndef __eigenpy_scipy_allocator_hpp__
21 changes: 21 additions & 0 deletions include/eigenpy/scipy-type.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,27 @@ struct EIGENPY_DLLAPI ScipyType {
: getInstance().csc_matrix_obj;
}

template <typename SparseMatrix>
static PyTypeObject const* get_pytype(
const Eigen::SparseMatrixBase<SparseMatrix>* ptr = nullptr) {
EIGENPY_UNUSED_VARIABLE(ptr);
return SparseMatrix::IsRowMajor ? getInstance().csr_matrix_type
: getInstance().csc_matrix_type;
}

static int get_numpy_type_num(const bp::object& obj) {
const PyTypeObject* type = Py_TYPE(obj.ptr());
EIGENPY_USED_VARIABLE_ONLY_IN_DEBUG_MODE(type);
assert(type == getInstance().csr_matrix_type ||
type == getInstance().csc_matrix_type);

bp::object dtype = obj.attr("dtype");

const PyArray_Descr* npy_type =
reinterpret_cast<PyArray_Descr*>(dtype.ptr());
return npy_type->type_num;
}

protected:
ScipyType();

Expand Down
Loading

0 comments on commit a2484d7

Please sign in to comment.