Merge pull request #411 from Anatoscope/sofapython-sparse

[SofaPython] sparse matrix aliasing scipy/eigen

SofaKernel/framework/sofa/defaulttype/DataTypeInfo.h

@@ -84,7 +84,10 @@ namespace defaulttype
     generically in non-template code.
 */
 template<class TDataType>
-struct DataTypeInfo
+struct DataTypeInfo;
+
+template<class TDataType>
+struct DefaultDataTypeInfo
 {
     /// Template parameter.
     typedef TDataType DataType;
@@ -139,6 +142,8 @@ struct DataTypeInfo
     {
     }
 
+    // mtournier: wtf is this supposed to do?
+    // mtournier: wtf is this not returning &type?
     static const void* getValuePtr(const DataType& /*type*/)
     {
         return NULL;
@@ -153,6 +158,11 @@ struct DataTypeInfo
 
 };
 
+template<class TDataType>
+struct DataTypeInfo : DefaultDataTypeInfo<TDataType> { };
+
+
+
 /// Type name template: default to using DataTypeInfo::name(), but can be overriden for types with shorter typedefs
 template<class TDataType>
 struct DataTypeName : public DataTypeInfo<TDataType>

SofaKernel/modules/SofaEigen2Solver/EigenBaseSparseMatrix.h

@@ -52,8 +52,6 @@ namespace linearsolver
 
 
 
-
-
 /** Sparse matrix based on the Eigen library.
 
 An Eigen::SparseMatrix<Real, RowMajor> matrix is used to store the data in Compressed Row Storage mode.
@@ -477,6 +475,31 @@ template<> inline const char* EigenBaseSparseMatrix<float>::Name()  { return "Ei
 
 } // namespace component
 
+
+namespace defaulttype {
+
+template<class Real>
+struct DataTypeInfo< component::linearsolver::EigenBaseSparseMatrix<Real> > 
+    : DefaultDataTypeInfo<component::linearsolver::EigenBaseSparseMatrix<Real> > {
+
+    using typename DataTypeInfo::DefaultDataTypeInfo::DataType;
+
+    static const char* name() {
+        return DataType::Name();
+    }
+
+    static const void* getValuePtr(const DataType& type) {
+        return &type.compressedMatrix;
+    }
+
+    static void* getValuePtr(DataType& type) {
+        return &type.compressedMatrix;
+    }
+
+};
+
+}
+
 } // namespace sofa
 
 #endif

SofaKernel/modules/SofaEigen2Solver/EigenSparseMatrix.h

@@ -652,6 +652,15 @@ template<> inline const char* EigenSparseMatrix<defaulttype::Vec3fTypes, default
     }
 
 
+namespace defaulttype {
+
+template<class TIn, class TOut>
+struct DataTypeInfo< component::linearsolver::EigenSparseMatrix<TIn, TOut> > 
+    : DataTypeInfo< typename component::linearsolver::EigenSparseMatrix<TIn, TOut>::Inherit > {
+
+};
+
+}
 
 } // namespace sofa
 

applications/plugins/Compliant/CMakeLists.txt

@@ -250,6 +250,8 @@ if( SofaPython_FOUND )
         mapping/PythonMultiMapping.cpp
         python/_Compliant.cpp
         python/Binding_AssembledSystem.cpp
+
+        # note: need to fix assembled solver in python before removing
         python/python.cpp
         )
 

applications/plugins/Compliant/python/Compliant/sparse.py

@@ -0,0 +1,113 @@
+import Sofa
+
+import numpy as np
+import scipy as sp
+from scipy import sparse
+from contextlib import contextmanager
+
+from ctypes import *
+
+dll_path = Sofa.loadPlugin('Compliant')
+dll = CDLL(dll_path)
+
+
+# note: we don't alias eigen matrices directly as the memory layout could change
+# between versions (it did once in the past IIRC), so we use a low-level
+# exchange data structure instead
+
+class Matrix(Structure):
+    '''all the data needed to alias a sparse matrix in eigen/scipy'''
+
+    _fields_ = (('rows', c_size_t),
+                ('cols', c_size_t),
+                ('outer_index', POINTER(c_int)),
+                ('inner_nonzero', POINTER(c_int)),
+                ('values', POINTER(c_double)),
+                ('indices', POINTER(c_int)),
+                ('size', c_size_t))
+
+
+    @staticmethod
+    def from_scipy(s):
+        data = Matrix()
+        values, inner_indices, outer_index = s.data, s.indices, s.indptr
+
+        data.rows, data.cols = s.shape
+
+        data.outer_index = outer_index.ctypes.data_as(POINTER(c_int))
+        data.inner_nonzero = None
+
+        data.values = values.ctypes.data_as(POINTER(c_double))
+        data.indices = inner_indices.ctypes.data_as(POINTER(c_int))
+        data.size = values.size
+
+        return data
+
+    @staticmethod
+    def from_eigen(ptr):
+        data = Matrix()
+        dll.eigen_to_scipy(byref(data), ptr)
+        return data
+
+
+    def to_eigen(self, ptr):
+        return dll.eigen_from_scipy(ptr, byref(self))
+
+
+    def to_scipy(self):
+        '''warning: if the scipy view reallocates, it will no longer alias the sparse
+        matrix.
+        '''
+
+        data = self
+
+        # needed: outer_index, data.values, data.size, data.indices, outer_size, inner_size
+        outer_index = np.ctypeslib.as_array( as_buffer(data.outer_index, data.rows + 1) )
+
+        shape = (data.rows, data.cols)
+
+        if not data.values:
+            return sp.sparse.csr_matrix( shape )
+
+        values = np.ctypeslib.as_array( as_buffer(data.values, data.size) )
+        inner_indices = np.ctypeslib.as_array( as_buffer(data.indices, data.size) )
+
+        return sp.sparse.csr_matrix( (values, inner_indices, outer_index), shape)
+
+
+    @staticmethod
+    @contextmanager
+    def view(ptr):
+        view = Matrix.from_eigen(ptr).to_scipy()
+        data = view.data.ctypes.data
+
+        try:
+            yield view
+        finally:
+            new = view.data.ctypes.data
+            if new != data:
+                # data pointer changed: rebuild view and assign back to eigen
+                Matrix.from_scipy(view).to_eigen(ptr)
+
+
+# sparse <- eigen        
+dll.eigen_to_scipy.restype = None
+dll.eigen_to_scipy.argtypes = (POINTER(Matrix), c_void_p)
+
+# eigen <- sparse
+dll.eigen_from_scipy.restype = None
+dll.eigen_from_scipy.argtypes = (c_void_p, POINTER(Matrix))
+
+
+def as_buffer(ptr, *size):
+    '''cast a ctypes pointer to a multidimensional array of given sizes'''
+
+    addr = addressof(ptr.contents)
+
+    buffer_type = type(ptr.contents)
+
+    for s in reversed(size):
+        buffer_type = buffer_type * s
+
+    return buffer_type.from_address(addr)
+

applications/plugins/SofaPython/CMakeLists.txt

@@ -124,6 +124,7 @@ set(SOURCE_FILES
     ScriptEvent.cpp
     ScriptFunction.cpp
     initSofaPython.cpp
+    ctypes.cpp
 )
 
 set(PYTHON_FILES

applications/plugins/SofaPython/ctypes.cpp

@@ -0,0 +1,154 @@
+// a bunch of ctypes-bound functions
+
+#include <Eigen/Sparse>
+#include <cstddef>
+#include <type_traits>
+
+template<class U>
+struct scipy_csr_matrix {
+
+    // SparseMatrix
+    std::size_t rows, cols;
+    int* outer_index;
+    int* inner_nonzero;
+
+    // CompressedStorage
+    struct storage_type {
+        U* values;
+        int* indices;
+        std::size_t size;
+    } storage;
+
+};
+
+
+template<class U>
+struct eigen_csr_matrix : Eigen::SparseMatrix<U, Eigen::RowMajor> {
+
+    // this is a dummy class that provides glue between scipy matrices and eigen
+    // matrices (mostly adapting pointers/shapes)
+
+
+private:
+    // that's right: use placement new to make sure the destructor is never
+    // called since memory is owned by scipy.
+    ~eigen_csr_matrix() { }
+
+public:
+
+    eigen_csr_matrix( const scipy_csr_matrix<U>* source ) {
+        this->m_outerSize = source->rows;
+        this->m_innerSize = source->cols;
+
+        this->m_outerIndex = source->outer_index;
+        this->m_innerNonZeros = source->inner_nonzero;
+
+        // same here
+        new (&this->m_data) eigen_compressed_storage(source->storage);
+    }
+
+
+    using eigen_compressed_storage_base = typename eigen_csr_matrix::Storage;
+
+    struct eigen_compressed_storage : eigen_compressed_storage_base {
+
+        eigen_compressed_storage(const typename scipy_csr_matrix<U>::storage_type& source) {
+            this->m_values = source.values;
+            this->m_indices = source.indices;
+            this->m_size = source.size;
+        }
+
+        typename scipy_csr_matrix<U>::storage_type to_scipy() const {
+            typename scipy_csr_matrix<U>::storage_type res;
+
+            res.size = this->m_size;
+
+            if( res.size ) {
+                res.values = this->m_values;
+                res.indices = this->m_indices;
+            } else {
+                // this keeps ctypes for yelling a warning
+                res.values = nullptr;
+                res.indices = nullptr;
+            }
+
+            return res;
+        }
+
+    };
+
+    scipy_csr_matrix<U> to_scipy() const {
+        scipy_csr_matrix<U> res;
+
+        res.rows = this->m_outerSize;
+        res.cols = this->m_innerSize;
+
+        res.outer_index = this->m_outerIndex;
+        res.inner_nonzero = this->m_innerNonZeros;
+
+        res.storage = static_cast<const eigen_compressed_storage&>(this->m_data).to_scipy();
+
+        return res;
+    }
+};
+
+
+template<class U>
+static void eigen_from_scipy_impl(eigen_csr_matrix<U>* lvalue,
+                                  const scipy_csr_matrix<U>* rvalue) {
+    // note: we use placement new to make sure destructor is never called
+    // since memory is owned by scipy
+
+    // note: damn you clang-3.4 you're supposed to be a c++11 compiler ffs
+    // typename std::aligned_union<0, eigen_csr_matrix<U> >::type storage;
+
+    union storage_type {
+        eigen_csr_matrix<U> matrix;
+        char bytes[0]; // ye olde c trick ahoy
+        storage_type() { }
+        ~storage_type() { }
+    } storage;
+
+    const eigen_csr_matrix<U>* alias = new (storage.bytes) eigen_csr_matrix<U>(rvalue);
+
+    *lvalue = *alias;
+    lvalue->makeCompressed();
+}
+
+
+
+extern "C" {
+
+    std::size_t eigen_sizeof_double() {
+        return sizeof(eigen_csr_matrix<double>);
+    }
+
+    void eigen_to_scipy_double(scipy_csr_matrix<double>* dst, const eigen_csr_matrix<double>* src) {
+        *dst = src->to_scipy();
+    }
+
+    void eigen_from_scipy_double(eigen_csr_matrix<double>* lvalue,
+                                 const scipy_csr_matrix<double>* rvalue) {    
+        eigen_from_scipy_impl<double>(lvalue, rvalue);
+    }
+
+
+    std::size_t eigen_sizeof_float() {
+        return sizeof(eigen_csr_matrix<float>);
+    }
+
+
+    void eigen_to_scipy_float(scipy_csr_matrix<float>* dst, const eigen_csr_matrix<float>* src) {
+        *dst = src->to_scipy();
+    }
+
+    void eigen_from_scipy_float(eigen_csr_matrix<float>* lvalue,
+                                const scipy_csr_matrix<float>* rvalue) {    
+        eigen_from_scipy_impl<float>(lvalue, rvalue);
+    }
+
+}
+
+
+
+