Merge pull request #185 from jcarpent/devel

Start exposition of complex Scalar types for autodiff

CMakeLists.txt

@@ -120,6 +120,9 @@ SET(${PROJECT_NAME}_HEADERS
   include/eigenpy/registration.hpp
   include/eigenpy/angle-axis.hpp
   include/eigenpy/quaternion.hpp
+  include/eigenpy/user-type.hpp
+  include/eigenpy/ufunc.hpp
+  include/eigenpy/register.hpp
   include/eigenpy/stride.hpp
   include/eigenpy/version.hpp
 )
@@ -150,6 +153,7 @@ SET(${PROJECT_NAME}_SOURCES
   src/matrix-float.cpp
   src/matrix-complex-float.cpp
   src/matrix-complex-double.cpp
+  src/register.cpp
   src/matrix-double.cpp
   src/matrix-long-double.cpp
   src/matrix-complex-long-double.cpp

cmake @ caad0fd6

-Subproject commit a61ae61479a1d50d1ae3c988d1d9aaf20841173d
+Subproject commit caad0fd63c5824c5cc352cc856b666d8842df6d1

include/eigenpy/eigen-allocator.hpp

@@ -7,6 +7,7 @@
 
 #include "eigenpy/fwd.hpp"
 #include "eigenpy/numpy-map.hpp"
+#include "eigenpy/register.hpp"
 #include "eigenpy/scalar-conversion.hpp"
 #include "eigenpy/utils/is-aligned.hpp"
 
@@ -114,14 +115,15 @@ namespace eigenpy
       Type * mat_ptr = details::init_matrix_or_array<Type>::run(pyArray,raw_ptr);
       Type & mat = *mat_ptr;
       
-      const int pyArray_Type = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
-      if(pyArray_Type == NumpyEquivalentType<Scalar>::type_code)
+      const int pyArray_type_code = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
+      const int Scalar_type_code = Register::getTypeCode<Scalar>();
+      if(pyArray_type_code == Scalar_type_code)
       {
         mat = NumpyMap<MatType,Scalar>::map(pyArray); // avoid useless cast
         return;
       }
       
-      switch(pyArray_Type)
+      switch(pyArray_type_code)
       {
         case NPY_INT:
           EIGENPY_CAST_FROM_PYARRAY_TO_EIGEN_MATRIX(MatType,int,Scalar,pyArray,mat);
@@ -158,11 +160,12 @@ namespace eigenpy
                      PyArrayObject * pyArray)
     {
       const MatrixDerived & mat = const_cast<const MatrixDerived &>(mat_.derived());
-      const int pyArray_Type = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
+      const int pyArray_type_code = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
+      const int Scalar_type_code = Register::getTypeCode<Scalar>();
       
       typedef typename NumpyMap<MatType,Scalar>::EigenMap MapType;
       
-      if(pyArray_Type == NumpyEquivalentType<Scalar>::type_code) // no cast needed
+      if(pyArray_type_code == Scalar_type_code) // no cast needed
       {
         MapType map_pyArray = NumpyMap<MatType,Scalar>::map(pyArray);
         if(mat.rows() == map_pyArray.rows())
@@ -172,7 +175,7 @@ namespace eigenpy
         return;
       }
       
-      switch(pyArray_Type)
+      switch(pyArray_type_code)
       {
         case NPY_INT:
           EIGENPY_CAST_FROM_EIGEN_MATRIX_TO_PYARRAY(MatType,Scalar,int,mat,pyArray);
@@ -219,8 +222,9 @@ namespace eigenpy
       typedef typename StrideType<MatType,Eigen::internal::traits<RefType>::StrideType::InnerStrideAtCompileTime, Eigen::internal::traits<RefType>::StrideType::OuterStrideAtCompileTime >::type NumpyMapStride;
 
       bool need_to_allocate = false;
-      const int pyArray_Type = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
-      if(pyArray_Type != NumpyEquivalentType<Scalar>::type_code)
+      const int pyArray_type_code = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
+      const int Scalar_type_code = Register::getTypeCode<Scalar>();
+      if(pyArray_type_code != Scalar_type_code)
         need_to_allocate |= true;
       if(    (MatType::IsRowMajor && (PyArray_IS_C_CONTIGUOUS(pyArray) && !PyArray_IS_F_CONTIGUOUS(pyArray)))
           || (!MatType::IsRowMajor && (PyArray_IS_F_CONTIGUOUS(pyArray) && !PyArray_IS_C_CONTIGUOUS(pyArray)))
@@ -246,13 +250,13 @@ namespace eigenpy
         new (raw_ptr) StorageType(mat_ref,pyArray,mat_ptr);
         
         RefType & mat = *reinterpret_cast<RefType*>(raw_ptr);
-        if(pyArray_Type == NumpyEquivalentType<Scalar>::type_code)
+        if(pyArray_type_code == Scalar_type_code)
         {
           mat = NumpyMap<MatType,Scalar>::map(pyArray); // avoid useless cast
           return;
         }
         
-        switch(pyArray_Type)
+        switch(pyArray_type_code)
         {
           case NPY_INT:
             EIGENPY_CAST_FROM_PYARRAY_TO_EIGEN_MATRIX(MatType,int,Scalar,pyArray,mat);
@@ -284,7 +288,7 @@ namespace eigenpy
       }
       else
       {
-        assert(pyArray_Type == NumpyEquivalentType<Scalar>::type_code);
+        assert(pyArray_type_code == Scalar_type_code);
         typename NumpyMap<MatType,Scalar,Options,NumpyMapStride>::EigenMap numpyMap = NumpyMap<MatType,Scalar,Options,NumpyMapStride>::map(pyArray);
         RefType mat_ref(numpyMap);
         new (raw_ptr) StorageType(mat_ref,pyArray);
@@ -311,8 +315,10 @@ namespace eigenpy
       typedef typename StrideType<MatType,Eigen::internal::traits<RefType>::StrideType::InnerStrideAtCompileTime, Eigen::internal::traits<RefType>::StrideType::OuterStrideAtCompileTime >::type NumpyMapStride;
 
       bool need_to_allocate = false;
-      const int pyArray_Type = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
-      if(pyArray_Type != NumpyEquivalentType<Scalar>::type_code)
+      const int pyArray_type_code = EIGENPY_GET_PY_ARRAY_TYPE(pyArray);
+      const int Scalar_type_code = Register::getTypeCode<Scalar>();
+      
+      if(pyArray_type_code != Scalar_type_code)
         need_to_allocate |= true;
       if(    (MatType::IsRowMajor && (PyArray_IS_C_CONTIGUOUS(pyArray) && !PyArray_IS_F_CONTIGUOUS(pyArray)))
           || (!MatType::IsRowMajor && (PyArray_IS_F_CONTIGUOUS(pyArray) && !PyArray_IS_C_CONTIGUOUS(pyArray)))
@@ -338,13 +344,13 @@ namespace eigenpy
         new (raw_ptr) StorageType(mat_ref,pyArray,mat_ptr);
         
         MatType & mat = *mat_ptr;
-        if(pyArray_Type == NumpyEquivalentType<Scalar>::type_code)
+        if(pyArray_type_code == Scalar_type_code)
         {
           mat = NumpyMap<MatType,Scalar>::map(pyArray); // avoid useless cast
           return;
         }
         
-        switch(pyArray_Type)
+        switch(pyArray_type_code)
         {
           case NPY_INT:
             EIGENPY_CAST_FROM_PYARRAY_TO_EIGEN_MATRIX(MatType,int,Scalar,pyArray,mat);
@@ -376,7 +382,7 @@ namespace eigenpy
       }
       else
       {
-        assert(pyArray_Type == NumpyEquivalentType<Scalar>::type_code);
+        assert(pyArray_type_code == Scalar_type_code);
         typename NumpyMap<MatType,Scalar,Options,NumpyMapStride>::EigenMap numpyMap = NumpyMap<MatType,Scalar,Options,NumpyMapStride>::map(pyArray);
         RefType mat_ref(numpyMap);
         new (raw_ptr) StorageType(mat_ref,pyArray);

include/eigenpy/eigen-to-python.hpp

@@ -25,7 +25,7 @@ namespace boost { namespace python {
     
 #ifndef BOOST_PYTHON_NO_PY_SIGNATURES
     inline PyTypeObject const*
-    get_pytype()const
+    get_pytype() const
     {
       return converter::registered_pytype<MatrixRef>::get_pytype();
     }

include/eigenpy/numpy-allocator.hpp

@@ -9,8 +9,11 @@
 #include "eigenpy/numpy-type.hpp"
 #include "eigenpy/eigen-allocator.hpp"
 
+#include "eigenpy/register.hpp"
+
 namespace eigenpy
 {
+
   template<typename MatType>
   struct NumpyAllocator
   {
@@ -20,8 +23,8 @@ namespace eigenpy
     {
       typedef typename SimilarMatrixType::Scalar Scalar;
       
-      PyArrayObject * pyArray = (PyArrayObject*) call_PyArray_SimpleNew(static_cast<int>(nd), shape,
-                                                                        NumpyEquivalentType<Scalar>::type_code);
+      const int code = Register::getTypeCode<Scalar>();
+      PyArrayObject * pyArray = (PyArrayObject*) call_PyArray_SimpleNew(static_cast<int>(nd), shape, code);
       
       // Copy data
       EigenAllocator<SimilarMatrixType>::copy(mat,pyArray);
@@ -42,8 +45,11 @@ namespace eigenpy
       
       if(NumpyType::sharedMemory())
       {
-        PyArrayObject * pyArray = (PyArrayObject*) call_PyArray_New(static_cast<int>(nd), shape,
-                                                                    NumpyEquivalentType<Scalar>::type_code,
+        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);
         
@@ -77,8 +83,11 @@ namespace eigenpy
       
       if(NumpyType::sharedMemory())
       {
-        PyArrayObject * pyArray = (PyArrayObject*) call_PyArray_New(static_cast<int>(nd), shape,
-                                                                    NumpyEquivalentType<Scalar>::type_code,
+        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<SimilarMatrixType &>(mat.derived()).data(),
                                                                     NPY_ARRAY_MEMORY_CONTIGUOUS_RO | NPY_ARRAY_ALIGNED);
                                                                     

include/eigenpy/numpy-type.hpp

@@ -9,12 +9,16 @@
 #include "eigenpy/scalar-conversion.hpp"
 
 #include <patchlevel.h> // For PY_MAJOR_VERSION
+#include <stdexcept>
+#include <typeinfo>
+#include <sstream>
 
 namespace eigenpy
 {
   namespace bp = boost::python;
 
-  template <typename SCALAR>  struct NumpyEquivalentType {};
+  // By default, the Scalar is considered as a Python object
+  template <typename Scalar> struct NumpyEquivalentType { enum  { type_code = NPY_USERDEF };};
 
   template <> struct NumpyEquivalentType<float>   { enum { type_code = NPY_FLOAT  };};
   template <> struct NumpyEquivalentType< std::complex<float> >   { enum { type_code = NPY_CFLOAT  };};
@@ -22,9 +26,18 @@ namespace eigenpy
   template <> struct NumpyEquivalentType< std::complex<double> >  { enum { type_code = NPY_CDOUBLE };};
   template <> struct NumpyEquivalentType<long double>  { enum { type_code = NPY_LONGDOUBLE };};
   template <> struct NumpyEquivalentType< std::complex<long double> >  { enum { type_code = NPY_CLONGDOUBLE };};
+  template <> struct NumpyEquivalentType<bool>    { enum { type_code = NPY_BOOL  };};
   template <> struct NumpyEquivalentType<int>     { enum { type_code = NPY_INT    };};
   template <> struct NumpyEquivalentType<long>    { enum { type_code = NPY_LONG    };};
 
+  template<typename Scalar>
+  bool isNumpyNativeType()
+  {
+    if(NumpyEquivalentType<Scalar>::type_code == NPY_USERDEF)
+      return false;
+    return true;
+  }
+
   template<typename Scalar>
   bool np_type_is_convertible_into_scalar(const int np_type)
   {

include/eigenpy/numpy.hpp

@@ -18,6 +18,7 @@
 #endif
 
 #include <numpy/noprefix.h>
+#include <numpy/ufuncobject.h>
 
 #if defined _WIN32 || defined __CYGWIN__
   #define EIGENPY_GET_PY_ARRAY_TYPE(array) \
@@ -30,21 +31,34 @@
 namespace eigenpy
 {
   void EIGENPY_DLLEXPORT import_numpy();
+  int EIGENPY_DLLEXPORT PyArray_TypeNum(PyTypeObject * type);
 }
 
 #if defined _WIN32 || defined __CYGWIN__
 namespace eigenpy
 {
-  EIGENPY_DLLEXPORT PyObject*  call_PyArray_SimpleNew(int nd, npy_intp * shape, int np_type);
+  EIGENPY_DLLEXPORT PyObject* call_PyArray_SimpleNew(int nd, npy_intp * shape, int np_type);
 
-  EIGENPY_DLLEXPORT PyObject* call_PyArray_New(int nd, npy_intp * shape, int np_type, void * data_ptr, int options);
+  EIGENPY_DLLEXPORT PyObject* call_PyArray_New(PyTypeObject * py_type_ptr, int nd, npy_intp * shape, int np_type, void * data_ptr, int options);
 
   EIGENPY_DLLEXPORT int call_PyArray_ObjectType(PyObject *, int);
+
+  EIGENPY_DLLEXPORT PyTypeObject * getPyArrayType();
+
+  EIGENPY_DLLEXPORT PyArray_Descr * call_PyArray_DescrFromType(int typenum);
+
+  EIGENPY_DLLEXPORT void call_PyArray_InitArrFuncs(PyArray_ArrFuncs * funcs);
+
+  EIGENPY_DLLEXPORT int call_PyArray_RegisterDataType(PyArray_Descr * dtype);
 }
 #else
   #define call_PyArray_SimpleNew PyArray_SimpleNew
-  #define call_PyArray_New(nd,shape,np_type,data_ptr,options) \
-    PyArray_New(&PyArray_Type,nd,shape,np_type,NULL,data_ptr,0,options,NULL)
+  #define call_PyArray_New(py_type_ptr,nd,shape,np_type,data_ptr,options) \
+    PyArray_New(py_type_ptr,nd,shape,np_type,NULL,data_ptr,0,options,NULL)
+  #define getPyArrayType() &PyArray_Type
+  #define call_PyArray_DescrFromType(typenum) PyArray_DescrFromType(typenum)
+  #define call_PyArray_InitArrFuncs(funcs) PyArray_InitArrFuncs(funcs)
+ #define call_PyArray_RegisterDataType(dtype) PyArray_RegisterDataType(dtype)
 #endif
 
 #endif // ifndef __eigenpy_numpy_hpp__

include/eigenpy/register.hpp

+//
+// Copyright (c) 2020 INRIA
+//
+
+#ifndef __eigenpy_register_hpp__
+#define __eigenpy_register_hpp__
+
+#include "eigenpy/fwd.hpp"
+#include "eigenpy/numpy-type.hpp"
+#include "eigenpy/exception.hpp"
+
+#include <algorithm>
+#include <map>
+#include <typeinfo>
+#include <string>
+
+namespace eigenpy
+{
+
+  /// \brief Structure collecting all the types registers in Numpy via EigenPy
+  struct EIGENPY_DLLEXPORT Register
+  {
+    
+    static PyArray_Descr * getPyArrayDescr(PyTypeObject * py_type_ptr);
+    
+    template<typename Scalar>
+    static bool isRegistered()
+    {
+      return isRegistered(Register::getPyType<Scalar>());
+    }
+    
+    static bool isRegistered(PyTypeObject * py_type_ptr);
+    
+    static int getTypeCode(PyTypeObject * py_type_ptr);
+  
+    template<typename Scalar>
+    static PyTypeObject * getPyType()
+    {
+      if(!isNumpyNativeType<Scalar>())
+      {
+        const PyTypeObject * const_py_type_ptr = bp::converter::registered_pytype<Scalar>::get_pytype();
+        if(const_py_type_ptr == NULL)
+        {
+          std::stringstream ss;
+          ss << "The type " << typeid(Scalar).name() << " does not have a registered converter inside Boot.Python." << std::endl;
+          throw std::invalid_argument(ss.str());
+        }
+        PyTypeObject * py_type_ptr = const_cast<PyTypeObject *>(const_py_type_ptr);
+        return py_type_ptr;
+      }
+      else
+      {
+        PyArray_Descr * new_descr = call_PyArray_DescrFromType(NumpyEquivalentType<Scalar>::type_code);
+        return new_descr->typeobj;
+      }
+    }
+    
+    template<typename Scalar>
+    static int getTypeCode()
+    {
+      if(isNumpyNativeType<Scalar>())
+        return NumpyEquivalentType<Scalar>::type_code;
+      else
+      {
+        const std::type_info & info = typeid(Scalar);
+        if(instance().type_to_py_type_bindings.find(&info) != instance().type_to_py_type_bindings.end())
+        {
+          PyTypeObject * py_type = instance().type_to_py_type_bindings[&info];
+          int code = instance().py_array_code_bindings[py_type];
+  
+          return code;
+        }
+        else
+          return -1; // type not registered
+      }
+    }
+    
+    static int registerNewType(PyTypeObject * py_type_ptr,
+                               const std::type_info * type_info_ptr,
+                               const int type_size,
+                               PyArray_GetItemFunc * getitem,
+                               PyArray_SetItemFunc * setitem,
+                               PyArray_NonzeroFunc * nonzero,
+                               PyArray_CopySwapFunc * copyswap,
+                               PyArray_CopySwapNFunc * copyswapn,
+                               PyArray_DotFunc * dotfunc);
+    
+    static Register & instance();
+    
+  private:
+    
+    Register() {};
+    
+    struct Compare_PyTypeObject
+    {
+      bool operator()(const PyTypeObject * a, const PyTypeObject * b) const
+      {
+        return std::string(a->tp_name) < std::string(b->tp_name);
+      }
+    };
+  
+    struct Compare_TypeInfo
+    {
+      bool operator()(const std::type_info * a, const std::type_info * b) const
+      {
+        return std::string(a->name()) < std::string(b->name());
+      }
+    };
+  
+    typedef std::map<const std::type_info *,PyTypeObject *,Compare_TypeInfo> MapInfo;
+    MapInfo type_to_py_type_bindings;
+    
+    typedef std::map<PyTypeObject *,PyArray_Descr *,Compare_PyTypeObject> MapDescr;
+    MapDescr py_array_descr_bindings;
+    
+    typedef std::map<PyTypeObject *,int,Compare_PyTypeObject> MapCode;
+    MapCode py_array_code_bindings;
+  
+  };
+  
+} // namespace eigenpy
+
+#endif // __eigenpy_register_hpp__

include/eigenpy/ufunc.hpp

+//
+// Copyright (c) 2020 INRIA
+//
+
+#ifndef __eigenpy_ufunc_hpp__
+#define __eigenpy_ufunc_hpp__
+
+#include "eigenpy/register.hpp"
+
+namespace eigenpy
+{
+  namespace internal
+  {
+  
+#define EIGENPY_REGISTER_BINARY_OPERATOR(name,op) \
+    template<typename T1, typename T2, typename R> \
+    void binary_op_##name(char** args, npy_intp * dimensions, npy_intp * steps, void * /*data*/) \
+    { \
+      npy_intp is0 = steps[0], is1 = steps[1], \
+      os = steps[2], n = *dimensions; \
+      char * i0 = args[0], *i1 = args[1], *o = args[2]; \
+      int k; \
+      for (k = 0; k < n; k++) \
+      { \
+        T1 & x = *static_cast<T1*>(static_cast<void*>(i0)); \
+        T2 & y = *static_cast<T2*>(static_cast<void*>(i1)); \
+        R & res = *static_cast<R*>(static_cast<void*>(o)); \
+        res = x op y; \
+        i0 += is0; i1 += is1; o += os; \
+      } \
+    } \
+    \
+    template<typename T> \
+    void binary_op_##name(char** args, npy_intp * dimensions, npy_intp * steps, void * data) \
+    { \
+      binary_op_##name<T,T,T>(args,dimensions,steps,data); \
+    }
+  
+    EIGENPY_REGISTER_BINARY_OPERATOR(add,+)
+    EIGENPY_REGISTER_BINARY_OPERATOR(subtract,-)
+    EIGENPY_REGISTER_BINARY_OPERATOR(multiply,*)
+    EIGENPY_REGISTER_BINARY_OPERATOR(divide,/)
+    EIGENPY_REGISTER_BINARY_OPERATOR(equal,==)
+    EIGENPY_REGISTER_BINARY_OPERATOR(not_equal,!=)
+    EIGENPY_REGISTER_BINARY_OPERATOR(less,<)
+    EIGENPY_REGISTER_BINARY_OPERATOR(greater,>)
+    EIGENPY_REGISTER_BINARY_OPERATOR(less_equal,<=)
+    EIGENPY_REGISTER_BINARY_OPERATOR(greater_equal,>=)
+  
+  #define EIGENPY_REGISTER_UNARY_OPERATOR(name,op) \
+    template<typename T, typename R> \
+    void unary_op_##name(char** args, npy_intp * dimensions, npy_intp * steps, void * /*data*/) \
+    { \
+      npy_intp is = steps[0], \
+      os = steps[1], n = *dimensions; \
+      char * i = args[0], *o = args[1]; \
+      int k; \
+      for (k = 0; k < n; k++) \
+      { \
+        T & x = *static_cast<T*>(static_cast<void*>(i)); \
+        R & res = *static_cast<R*>(static_cast<void*>(o)); \
+        res = op x; \
+        i += is; o += os; \
+      } \
+    } \
+    \
+    template<typename T> \
+    void unary_op_##name(char** args, npy_intp * dimensions, npy_intp * steps, void * data) \
+    { \
+      unary_op_##name<T,T>(args,dimensions,steps,data); \
+    }
+  
+    EIGENPY_REGISTER_UNARY_OPERATOR(negative,-)
+  
+  } // namespace internal
+  
+#define EIGENPY_REGISTER_BINARY_UFUNC(name,code,T1,T2,R) { \
+   PyUFuncObject* ufunc = \
+       (PyUFuncObject*)PyObject_GetAttrString(numpy, #name); \
+   int _types[3] = { Register::getTypeCode<T1>(), Register::getTypeCode<T2>(), Register::getTypeCode<R>()}; \
+   if (!ufunc) { \
+       /*goto fail; \*/ \
+   } \
+   if (sizeof(_types)/sizeof(int)!=ufunc->nargs) { \
+       PyErr_Format(PyExc_AssertionError, \
+                    "ufunc %s takes %d arguments, our loop takes %lu", \
+                    #name, ufunc->nargs, (unsigned long) \
+                    (sizeof(_types)/sizeof(int))); \
+       Py_DECREF(ufunc); \
+   } \
+   if (PyUFunc_RegisterLoopForType((PyUFuncObject*)ufunc, code, \
+        internal::binary_op_##name<T1,T2,R>, _types, 0) < 0) { \
+       /*Py_DECREF(ufunc);*/ \
+       /*goto fail; \*/ \
+   } \
+   Py_DECREF(ufunc); \
+}
+  
+#define EIGENPY_REGISTER_UNARY_UFUNC(name,code,T,R) { \
+   PyUFuncObject* ufunc = \
+       (PyUFuncObject*)PyObject_GetAttrString(numpy, #name); \
+   int _types[2] = { Register::getTypeCode<T>(), Register::getTypeCode<R>()}; \
+   if (!ufunc) { \
+       /*goto fail; \*/ \
+   } \
+   if (sizeof(_types)/sizeof(int)!=ufunc->nargs) { \
+       PyErr_Format(PyExc_AssertionError, \
+                    "ufunc %s takes %d arguments, our loop takes %lu", \
+                    #name, ufunc->nargs, (unsigned long) \
+                    (sizeof(_types)/sizeof(int))); \
+       Py_DECREF(ufunc); \
+   } \
+   if (PyUFunc_RegisterLoopForType((PyUFuncObject*)ufunc, code, \
+        internal::unary_op_##name<T,R>, _types, 0) < 0) { \
+       /*Py_DECREF(ufunc);*/ \
+       /*goto fail; \*/ \
+   } \
+   Py_DECREF(ufunc); \
+}
+
+  template<typename Scalar>
+  void registerCommonUfunc()
+  {
+    const int code = Register::getTypeCode<Scalar>();
+  
+    PyObject* numpy_str;
+#if PY_MAJOR_VERSION >= 3
+    numpy_str = PyUnicode_FromString("numpy");
+#else
+    numpy_str = PyString_FromString("numpy");
+#endif
+    PyObject* numpy;
+    numpy = PyImport_Import(numpy_str);
+    Py_DECREF(numpy_str);
+    
+    import_ufunc();
+
+    // Binary operators
+    EIGENPY_REGISTER_BINARY_UFUNC(add,code,Scalar,Scalar,Scalar);
+    EIGENPY_REGISTER_BINARY_UFUNC(subtract,code,Scalar,Scalar,Scalar);
+    EIGENPY_REGISTER_BINARY_UFUNC(multiply,code,Scalar,Scalar,Scalar);
+    EIGENPY_REGISTER_BINARY_UFUNC(divide,code,Scalar,Scalar,Scalar);
+  
+    // Comparison operators
+    EIGENPY_REGISTER_BINARY_UFUNC(equal,code,Scalar,Scalar,bool);
+    EIGENPY_REGISTER_BINARY_UFUNC(not_equal,code,Scalar,Scalar,bool);
+    EIGENPY_REGISTER_BINARY_UFUNC(greater,code,Scalar,Scalar,bool);
+    EIGENPY_REGISTER_BINARY_UFUNC(less,code,Scalar,Scalar,bool);
+    EIGENPY_REGISTER_BINARY_UFUNC(greater_equal,code,Scalar,Scalar,bool);
+    EIGENPY_REGISTER_BINARY_UFUNC(less_equal,code,Scalar,Scalar,bool);
+  
+    // Unary operators
+    EIGENPY_REGISTER_UNARY_UFUNC(negative,code,Scalar,Scalar);
+
+    Py_DECREF(numpy);
+  }
+  
+} // namespace eigenpy
+
+#endif // __eigenpy_ufunc_hpp__

include/eigenpy/user-type.hpp

+//
+// Copyright (c) 2020 INRIA
+//
+
+#ifndef __eigenpy_user_type_hpp__
+#define __eigenpy_user_type_hpp__
+
+#include "eigenpy/fwd.hpp"
+#include "eigenpy/numpy-type.hpp"
+#include "eigenpy/register.hpp"
+
+namespace eigenpy
+{
+  namespace internal
+  {
+    template<typename T, int type_code = NumpyEquivalentType<T>::type_code>
+    struct SpecialMethods
+    {
+      static void copyswap(void * /*dst*/, void * /*src*/, int /*swap*/, void * /*arr*/) {};
+      static PyObject * getitem(void * /*ip*/, void * /*ap*/) { return NULL; };
+      static int setitem(PyObject * /*op*/, void * /*ov*/, void * /*ap*/) { return -1; }
+      static void copyswapn(void * /*dest*/, long /*dstride*/, void * /*src*/,
+                            long /*sstride*/, long /*n*/, int /*swap*/, void * /*arr*/) {};
+      static npy_bool nonzero(void * /*ip*/, void * /*array*/) { return (npy_bool)false; };
+      static void dotfunc(void * /*ip0_*/, npy_intp /*is0*/, void * /*ip1_*/, npy_intp /*is1*/,
+                          void * /*op*/, npy_intp /*n*/, void * /*arr*/);
+//      static void cast(void * /*from*/, void * /*to*/, npy_intp /*n*/, void * /*fromarr*/, void * /*toarr*/) {};
+    };
+  
+    template<typename T>
+    struct SpecialMethods<T,NPY_USERDEF>
+    {
+      static void copyswap(void * dst, void * src, int swap, void * /*arr*/)
+      {
+//        std::cout << "copyswap" << std::endl;