possible to use eigenpy with custom numeric scalar types?

[ofloveandhate]

My package uses Boost.Multiprecision complex numbers. Is it possible to extend eigenpy to include Eigen matrices / vectors / etc with custom scalar types?

[jcarpent]

Good question and good feature of Eigenpy ;)

Yes, you can.
You just need to declare the type converters for that.
I've done that for Pinocchio.
Here is a code taken from Pinocchio 3x version, not yet online.

//
// Copyright (c) 2021-2022 INRIA
//

#ifndef __pinocchio_python_math_multiprecision_boost_number_hpp__
#define __pinocchio_python_math_multiprecision_boost_number_hpp__

#include "pinocchio/math/multiprecision.hpp"

#include <boost/python.hpp>
#include <boost/python/return_value_policy.hpp>
#include <eigenpy/user-type.hpp>
#include <eigenpy/ufunc.hpp>
#include <sstream>

namespace
{

  template<class Backend>
  struct get_backend_precision {};

  template <unsigned Digits10, ::boost::multiprecision::mpfr_allocation_type AllocateType>
  struct get_backend_precision<::boost::multiprecision::mpfr_float_backend<Digits10,AllocateType>>
  {
    enum { value = Digits10 };
  };

}

namespace eigenpy
{
  namespace internal
  {
    template <class Backend, ::boost::multiprecision::expression_template_option ExpressionTemplates>
    struct getitem<::boost::multiprecision::number<Backend,ExpressionTemplates>>
    {
      
      typedef ::boost::multiprecision::number<Backend,ExpressionTemplates> Scalar;

      static PyObject* run(void* data, void* /* arr */) {
        Scalar & mpfr_scalar = *static_cast<Scalar*>(data);
        Backend & backend = mpfr_scalar.backend();
       
        if(backend.data()[0]._mpfr_d == 0) // If the mpfr_scalar is not initialized, we have to init it.
        {
          mpfr_scalar = Scalar(0);
//          unsigned int precision = get_backend_precision<Backend>::value ? get_backend_precision<Backend>::value : backend.default_precision();
//          mpfr_init2(backend.data(), ::boost::multiprecision::detail::digits10_2_2(precision));
        }
        bp::object m(boost::ref(mpfr_scalar));
        Py_INCREF(m.ptr());
        return m.ptr();
      }
    };
  }

}

namespace pinocchio
{
  namespace python
  {
    namespace bp = boost::python;

    template<typename BoostNumber>
    struct BoostNumberPythonVisitor
    : public boost::python::def_visitor< BoostNumberPythonVisitor<BoostNumber> >
    {

    public:

      template<class PyClass>
      void visit(PyClass& cl) const
      {
        cl
        .def(bp::init<>("Default constructor.",bp::arg("self")))
        .def(bp::init<BoostNumber>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<bool>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<float>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<double>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<int>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<long int>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<unsigned int>("Copy constructor.",bp::args("self","value")))
//        .def(bp::init<unsigned long int>("Copy constructor.",bp::args("self","value")))
        .def(bp::init<std::string>("Constructor from a string.",bp::args("self","str_value")))
        
PINOCCHIO_COMPILER_DIAGNOSTIC_PUSH
PINOCCHIO_COMPILER_DIAGNOSTIC_IGNORED_SELF_ASSIGN_OVERLOADED
        .def(bp::self +  bp::self)
        .def(bp::self += bp::self)
        .def(bp::self -  bp::self)
        .def(bp::self -= bp::self)
        .def(bp::self *  bp::self)
        .def(bp::self *= bp::self)
        .def(bp::self /  bp::self)
        .def(bp::self /= bp::self)
        
        .def(bp::self <  bp::self)
        .def(bp::self <= bp::self)
        .def(bp::self >  bp::self)
        .def(bp::self >= bp::self)
        .def(bp::self == bp::self)
        .def(bp::self != bp::self)
        .def(bp::self_ns::pow(bp::self_ns::self,long()))
PINOCCHIO_COMPILER_DIAGNOSTIC_POP

        .def("str",&BoostNumber::str,bp::args("self","precision","scientific"))

        .def("default_precision",
             static_cast<unsigned (*)()>(BoostNumber::default_precision),
             "Get the default precision of the class.")
        .def("default_precision",
             static_cast<void (*)(unsigned)>(BoostNumber::default_precision),bp::arg("digits10"),
             "Set the default precision of the class.")
        .staticmethod("default_precision")

        .def("precision",
             static_cast<unsigned (BoostNumber::*)() const>(&BoostNumber::precision),
             bp::arg("self"),
             "Get the precision of this.")
        .def("precision",
             static_cast<void (BoostNumber::*)(unsigned)>(&BoostNumber::precision),
             bp::args("self","digits10"),
             "Set the precision of this.")
        
        .def("__float__",&cast<double>,bp::arg("self"),"Cast to float.")
        .def("__int__",&cast<int64_t>,bp::arg("self"),"Cast to int.")
        
        .def("__str__",&print,bp::arg("self"))
        .def("__repr__",&print,bp::arg("self"))
        
        .def("set_display_precision",&set_display_precision,bp::arg("digit"),
             "Set the precision when printing values.")
        .staticmethod("set_display_precision")
        
        .def("get_display_precision",&get_display_precision,
             "Get the precision when printing values.",
             bp::return_value_policy<bp::copy_non_const_reference>())
        .staticmethod("get_display_precision")

//#ifndef PINOCCHIO_PYTHON_NO_SERIALIZATION
//        .def_pickle(Pickle())
//#endif
        ;
        
      }
      
      static void expose(const std::string & type_name)
      {
        bp::class_<BoostNumber>(type_name.c_str(),
                                "",bp::no_init)
        .def(BoostNumberPythonVisitor<BoostNumber>())
        ;
        
        eigenpy::registerNewType<BoostNumber>();
        eigenpy::registerCommonUfunc<BoostNumber>();
        
#define IMPLICITLY_CONVERTIBLE(T1,T2) \
  bp::implicitly_convertible<T1,T2>();
//  bp::implicitly_convertible<T2,T1>();
        
        IMPLICITLY_CONVERTIBLE(double,BoostNumber);
        IMPLICITLY_CONVERTIBLE(float,BoostNumber);
        IMPLICITLY_CONVERTIBLE(long int,BoostNumber);
        IMPLICITLY_CONVERTIBLE(int,BoostNumber);
        IMPLICITLY_CONVERTIBLE(long,BoostNumber);
        IMPLICITLY_CONVERTIBLE(unsigned int,BoostNumber);
        IMPLICITLY_CONVERTIBLE(unsigned long int,BoostNumber);
        IMPLICITLY_CONVERTIBLE(bool,BoostNumber);
        
#undef IMPLICITLY_CONVERTIBLE

        eigenpy::registerCast<BoostNumber,double>(false);
        eigenpy::registerCast<double,BoostNumber>(true);
        eigenpy::registerCast<BoostNumber,float>(false);
        eigenpy::registerCast<float,BoostNumber>(true);
        eigenpy::registerCast<BoostNumber,long>(false);
        eigenpy::registerCast<long,BoostNumber>(true);
        eigenpy::registerCast<BoostNumber,int>(false);
        eigenpy::registerCast<int,BoostNumber>(true);;
        eigenpy::registerCast<BoostNumber,int64_t>(false);
        eigenpy::registerCast<int64_t,BoostNumber>(true);
      }
      
    private:
      
      template<typename T>
      static T cast(const BoostNumber & self)
      {
        return static_cast<T>(self);
      }
      
      static std::string print(const BoostNumber & self)
      {
        return self.str(get_display_precision(), std::ios_base::dec);
      }
      
      static void set_display_precision(const int digit)
      {
        get_display_precision() = digit;
      }
      
      static int & get_display_precision()
      {
        static int precision = BoostNumber::default_precision();
        return precision;
      }
      
    };
    
  } // namespace python
} // namespace pinocchio

#endif // ifndef __pinocchio_python_math_multiprecision_boost_number_hpp__

[ofloveandhate]

This is very helpful example code. I already have exposed Boost.Multiprecision variable precision types, so the main two lines I think I need to interact with EigenPy are:

eigenpy::registerNewType<BoostNumber>();
eigenpy::registerCommonUfunc<BoostNumber>();

and probably the casting registrations between numeric types.

Thank you.

[ofloveandhate]

I'm trying to register a complex number type and failing, because comparitors like < are not defined for complex numbers. Here's a partial traceback for my compiler error:

/opt/homebrew/include/eigenpy/ufunc.hpp:124:43: error: invalid operands to binary expression ('boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>' and 'boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>')
EIGENPY_REGISTER_BINARY_OPERATOR(greater, >)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~
/opt/homebrew/include/eigenpy/ufunc.hpp:103:15: note: expanded from macro 'EIGENPY_REGISTER_BINARY_OPERATOR'
      res = x op y;                                                          \
            ~ ^  ~
/opt/homebrew/include/eigenpy/ufunc.hpp:257:3: note: in instantiation of function template specialization 'eigenpy::internal::binary_op_greater<boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>, boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>, bool>' requested here
  EIGENPY_REGISTER_BINARY_UFUNC(greater, type_code, Scalar, Scalar, bool);
  ^
/opt/homebrew/include/eigenpy/ufunc.hpp:173:47: note: expanded from macro 'EIGENPY_REGISTER_BINARY_UFUNC'
                                    internal::binary_op_##name<T1, T2, R>,     \
                                              ^
<scratch space>:20:1: note: expanded from here
binary_op_greater
^
src/eigenpy_interaction.cpp:42:11: note: in instantiation of function template specialization 'eigenpy::registerCommonUfunc<boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>>' requested here
        eigenpy::registerCommonUfunc<bertini::mpfr_complex>();
                 ^

I see that eigenpy/ufunc.hpp has some lines in it that seem to require that user types have these comparisons implemented.

Do you have any advice?

[jcarpent]

I would suggest only using the lines in registerCommonUfunc which are aligned with your type to register.

[ofloveandhate]

I think I see what you're saying. I see the function void registerCommonUfunc() in ufunc.hpp. Are you suggesting I derive from this, and omit the lines corresponding to unimplemented functionality?

[jcarpent]

Yes, this is totally what I suggest indeed.

[ofloveandhate]

Groovy. I did it and it compiled and linked. Now I'm getting segfaults, but that's my own adventure right now. Thank you for your assistance!

[jcarpent]

With pleasure. I hope EigenPy could be useful in your research ;)

[ofloveandhate]

Do you know, when I link my shared library for Python, do I need to link against eigenpy? or numpy? I'm struggling to get past my segfault issue.

[jcarpent]

Yes, you should link to eigenpy but not to Python shared lib. Could you share your CMakeLists.txt?

[ofloveandhate]

Sadly, I'm using autoconf, for better or for worse. I have an open issue to convert to CMake. Perhaps now is the time? Ugh, it's a big move and I've been avoiding it...

[ofloveandhate]

Probing the links to libararies from my built .so file, I find that I linked to eigenpy:

otool -L /usr/local/lib/python3.11/site-packages/_pybertini.so 
/usr/local/lib/python3.11/site-packages/_pybertini.so:
	/opt/homebrew/opt/python@3.11/Frameworks/Python.framework/Versions/3.11/Python (compatibility version 3.11.0, current version 3.11.0)
	/opt/homebrew/opt/boost-python3/lib/libboost_python311-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/eigenpy/lib/libeigenpy.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/local/lib/libbertini2.1.dylib (compatibility version 2.0.0, current version 2.0.0)
	/opt/homebrew/opt/boost/lib/libboost_filesystem.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_system.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_chrono-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_regex-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_timer-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_unit_test_framework-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_serialization-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_log-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_log_setup-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost/lib/libboost_thread-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1319.100.3)
	/opt/homebrew/opt/libmpc/lib/libmpc.3.dylib (compatibility version 7.0.0, current version 7.1.0)
	/opt/homebrew/opt/mpfr/lib/libmpfr.6.dylib (compatibility version 9.0.0, current version 9.0.0)
	/opt/homebrew/opt/gmp/lib/libgmp.10.dylib (compatibility version 15.0.0, current version 15.1.0)
	/usr/lib/libc++.1.dylib (compatibility version 1.0.0, current version 1500.65.0)

---

I also find that eigenpy and my library are linked against the same versions of boost-python, so I don't think that's the issue. I'm stumped, but not giving up!

[ofloveandhate]

Here's the list of things that eigenpy is linked against:

 otool -L /opt/homebrew/opt/eigenpy/lib/libeigenpy.dylib
/opt/homebrew/opt/eigenpy/lib/libeigenpy.dylib:
	/opt/homebrew/opt/eigenpy/lib/libeigenpy.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/boost-python3/lib/libboost_python311-mt.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libc++.1.dylib (compatibility version 1.0.0, current version 1500.65.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1319.100.3)

Numpy is not linked, so I think I don't need to link. I wonder if I found incorrect version of numpy headers?

[ofloveandhate]

I finally got to a useful backtrace using lldb on Python as I try to import my library. Here's what I'm getting:

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0x50)
    frame #0: 0x0000000104965160 libeigenpy.dylib`eigenpy::Register::registerNewType(_typeobject*, std::type_info const*, int, int, _object* (*)(void*, void*), int (*)(_object*, void*, void*), unsigned char (*)(void*, void*), void (*)(void*, void*, int, void*), void (*)(void*, long, void*, long, long, int, void*), void (*)(void*, long, void*, long, void*, long, void*), int (*)(void*, long, void*), int (*)(void*, long, void*, void*)) + 80
libeigenpy.dylib`eigenpy::Register::registerNewType:

So something isn't right about my call to registerNewType.