How to pass numpy array to compiled expression member variable?

1

Coming back to an old problem I never got to solve properly… I need to pass a numpy array to a compiled expression, to no avail. Consider the following code:

import numpy
import dolfin
mesh = dolfin.UnitSquareMesh(1,1)
fe = dolfin.FiniteElement(
    family="Quadrature",
    cell=mesh.ufl_cell(),
    degree=1,
    quad_scheme="default")
cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    const Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        const Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This compiles and runs fine but of course, since the member pointer v is const, even after assigning an array to it in the init_v function (actually I’m not even sure how/why I can do that since it is const, but anyway), it is empty when executing the print_v function, and equal to first parameter of the function when executing the eval function (so weird, right?). Now, if I make the member pointer not const:

cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        const Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This won’t compile with the error:

error: invalid conversion from ‘const dolfin::Array<double>*’ to ‘dolfin::Array<double>*’ [-fpermissive]
     v = &v_;
       ^

This makes sense. Now, if I remove the other const as well:

cpp_code = '''\
namespace dolfin
{
class MyExpr : public Expression
{
    Array<double>* v;
public:
    MyExpr(): Expression() {}
    void print_v() const
    {
        std::cout << "v = " << v->str(1) << std::endl;
    }
    void init_v(
        Array<double>& v_)
    {
        print_v();
        v = &v_;
        print_v();
    }
    void eval(
              Array<double> &expr,
        const Array<double> &pos) const
    {
        print_v();
    }
};
}'''
my_expr = dolfin.Expression(
    cppcode=cpp_code,
    element=fe)
my_expr.init_v(numpy.ones(2))
my_expr.print_v()
my_expr.eval(numpy.ones(1), numpy.ones(3))

This compiles fine, but at runtime I’m getting the following error:

my_expr.init_v(numpy.ones(2))
TypeError: in method 'MyExpr_init_v', argument 2 of type 'dolfin::Array< double > &'

So apparently, the compiled module recognizes a numpy array as const dolfin::Array< double > & , but not as dolfin::Array< double > &. Any idea on how to make this work? Thanks! Martin

2

Hey Martin,

Here’s a tester code I tried to pass a numpy array to a dolfin expression in C++ . It should help :

import dolfin as d
import numpy as np

code = """
    
    #include <pybind11/pybind11.h>
    #include <pybind11/eigen.h>
    #include <dolfin/function/Expression.h>
    
    class test_exp : public dolfin::Expression {
      public:
        
        Eigen::VectorXd arr;
        
        test_exp() : dolfin::Expression() { }

        void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const {
            
            values[0] = arr.sum();
        }

    };
    
    PYBIND11_MODULE(SIGNATURE, m) {
        pybind11::class_<test_exp, std::shared_ptr<test_exp>, dolfin::Expression>
        (m, "test_exp")
        .def(pybind11::init<>())
        .def_readwrite("arr", &test_exp::arr)
        ;
    }
    """

exp = d.CompiledExpression(d.compile_cpp_code(code).test_exp(), degree=1)
exp.arr = np.array([1., 2., 4., 8.], dtype=float)
print(exp([0.,0.,0.]))   # Prints 15

Note : Initialize the exp.arr object in python to avoid seg_faults.
Ref: test_jit.py
(Tested on 2019.1.0)

3 Likes
3

OMG thank you!! (test_jit.py is such a useful source of solution!) Just wondering: is there a similar procedure with dolfin::Array object?

4

Hey Martin,

I haven’t personally tried using the dolfin::Array class as an interface between Python and C++, but I think it should technically be possible. I can maybe suggest you this -

Interfacing NumPy objects in pybind is possible by exposing pointers to the data buffer. Maybe you can pass along the buffer as reference and construct a dolfin::Array object within C++ code and use it? (may require verifying data type compatibility and conversion)

If you do use the interface between Python and C++, I think using Eigen classes with Numpy is the easiest as in the earlier answer. You can read a little more on it here. Keep in mind that there is a proposal in FEniCS to adapt to this (Ref), although it hasn’t yet been sorted.

1 Like
5

All clear, thank you so much!

1 Like
6

One issue is:

v = np.array([1., 2., 4., 8.], dtype=float)
exp.arr = v
print(exp([0.,0.,0.]))   # Prints 15, OK
v *= 2
print(exp([0.,0.,0.]))   # Still prints 15, not 30
exp.arr = v
print(exp([0.,0.,0.]))   # Prints 30, OK

So I guess exp.arr = v actually involves some copy, right? Is there a way to make the two variables (python & c++) linked to the same data? Thanks again so much.

1 Like
7

The following modified example seems to do the trick, although there is probably a “cleaner” approach:

import dolfin as d
import numpy as np

code = """
    
    #include <pybind11/eigen.h>
    #include <dolfin/function/Expression.h>
    #include <pybind11/pybind11.h>
    #include <pybind11/numpy.h>
    namespace py = pybind11;

    typedef py::array_t<double, py::array::c_style | py::array::forcecast> 
      npArray;
    
    class test_exp : public dolfin::Expression {
      public:

        double *arr;
        int n;
        
        test_exp(npArray na) : dolfin::Expression(){
           auto nab = na.request();
           arr = (double *)nab.ptr;
           n = nab.shape[0];
        }

        void eval(Eigen::Ref<Eigen::VectorXd> values, 
                  Eigen::Ref<const Eigen::VectorXd> x) const {
            
            values[0] = 0.0;
            for(int i=0; i<n; i++)
              values[0] += arr[i];
        }

    };
    
    PYBIND11_MODULE(SIGNATURE, m) {
        pybind11::class_<test_exp, std::shared_ptr<test_exp>, 
                         dolfin::Expression>
        (m, "test_exp")
        .def(pybind11::init<npArray>())
        ;
    }
    """
arr = np.array([1., 2., 4., 8.], dtype=float)
exp = d.CompiledExpression(d.compile_cpp_code(code).test_exp(arr), degree=1)
print(exp([0.,0.,0.]))   # Prints 15
arr *= 2.0
print(exp([0.,0.,0.]))   # Prints 30
3 Likes
8

Thanks so much David!

By the way, Matthias Rambausek pointed me to this web page, which has additional info: https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html

9

Thanks Kamensky!

Just to add: Using the Eigen::Ref<> template instead of py::array_t<> , pybind11 casts numpy arrays directly and also gives access to the Eigen library methods.

i.e Change npArray definition to

typedef Eigen::Ref<Eigen::VectorXd> npArray;

arr declaration to npArray arr;

and the constructor method to

test_exp(npArray a) : dolfin::Expression(), arr(a) {}

This can help in using Eigen Methods e.g arr.sum() in the eval method.

1 Like
10

Does anyone knows how to call make “exp” call
“void eval(Eigen::RefEigen::VectorXd values,
Eigen::Ref x,
const ufc::cell& cell)” instead of
“void eval(Eigen::RefEigen::VectorXd values,
Eigen::Ref x)”
I want cell information.