search

Chaste / PyChaste

A Python Wrapper for Chaste
https://chaste.github.io/pychaste/
Other
1 stars 0 forks source link

Ublas vector wrapping support is broken #11

Closed kwabenantim closed 1 year ago

kwabenantim commented 1 year ago

Describe the bug See PyChaste/issues/4

It is currently not possible to return UBlas vectors from C++ to Python, which breaks some tests and tutorials.

To Reproduce Install PyChaste from conda:

mamba install -c pychaste -c conda-forge -c bioconda chaste

In Python:

import numpy as np

import chaste
import chaste.cell_based

chaste.init()

generator = chaste.mesh.HoneycombVertexMeshGenerator(2, 2)
mesh = generator.GetMesh()

transit_type = chaste.cell_based.TransitCellProliferativeType()
cell_generator = chaste.cell_based.CellsGeneratorBiasedBernoulliTrialCellCycleModel_2()
cells = cell_generator.GenerateBasicRandom(mesh.GetNumElements(), transit_type)

cell_population = chaste.cell_based.VertexBasedCellPopulation2(mesh, cells)

point = np.array([0.0, 0.0])
normal = np.array([0.0, 1.0])
boundary_condition = chaste.cell_based.PlaneBoundaryCondition2_2(cell_population, point, normal)

Error

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
[<ipython-input-12-3f9294ac1a2f>](https://sk6pcemouof-496ff2e9c6d22116-0-colab.googleusercontent.com/outputframe.html?vrz=colab-20230614-060133-RC01_540321000#) in <cell line: 14>()
     12 point = np.array([0.0, 0.0])
     13 normal = np.array([0.0, 1.0])
---> 14 boundary_condition = chaste.cell_based.PlaneBoundaryCondition2_2(cell_population, point, normal)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. chaste.cell_based._chaste_project_PyChaste_cell_based.PlaneBoundaryCondition2_2(pCellPopulation: AbstractCellPopulation<2u, 2u>, point: boost::numeric::ublas::c_vector<double, 2ul>, normal: boost::numeric::ublas::c_vector<double, 2ul>)

Invoked with: <chaste.cell_based._chaste_project_PyChaste_cell_based.VertexBasedCellPopulation2 object at 0x7f82c576bfb0>, array([0., 0.]), array([0., 1.])
kwabenantim commented 1 year ago

Also add support for zero_vector and unit_vector:

MILeach commented 1 year ago

The existing PYBIND11_CVECTOR_TYPECASTER3 macro can be used to solve this issue. It should be called in any cppwg.cpp file which contains methods using c_vector<double, 3>.

E.g, to solve the example described in this issue,

PlaneBoundaryCondition2_2.cppwg.cpp should include PythonObjectConverters.hpp

and the line

PYBIND11_CVECTOR_TYPECASTER3(); should be added to the file. This will enable automatic conversion between c_vector<double, 3> and python lists withing this file.

MILeach commented 1 year ago

The solution provided in the previous comment only appears to fix some of the failures.

Further investigation shows that the load method of the typecaster is called and successfully returns true without error. Locating the error message given in the pybind11 source code has led to identifying the cause of failure as a reference_cast_error being caught at ~line 930 of pybind11.h.

MILeach commented 1 year ago

Using the test script

import os, signal
import numpy as np
import chaste
import chaste.cell_based
import chaste.mesh

from io import StringIO
from wurlitzer import pipes, STDOUT

n = chaste.mesh.Node3(0, [0.0, 0.0, 0.0])

out = StringIO()
with pipes(stdout=out, stderr=out):
    try:
        n.AddAppliedForceContribution(np.array([0.1, 0.0, 1.0]))
    except TypeError:
        print("reference error caught")

print(out.getvalue())

The program throws a reference_cast_error if the second value of the parameter to AddAppliedForceContribution is 0.0.

I have investigated the call stack to this point with gdb, and for both zero and non-zero values, the types involved are identical. The only difference is that in the non-zero ones, the base type_caster_generic::value has been set. The next step is to investigate where & why this is/isn't set.

GDB trace for non-zero value given below:

942                     result = func.impl(call);
(gdb) s
0x00007fffeeb87485 in pybind11::cpp_function::initialize<pybind11::cpp_function::initialize<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}, void, Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(pybind11::cpp_function::initialize<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&&, void (*)(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(pybind11::detail::function_call&)#3}::_FUN(pybind11::detail::function_call&) () at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/pybind11.h:224
224         rec->impl = [](function_call &call) -> handle {
(gdb) s
0x00007fffeeb86d1b in pybind11::cpp_function::initialize<pybind11::cpp_function::initialize<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}, void, Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(pybind11::cpp_function::initialize<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&&, void (*)(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(pybind11::detail::function_call&)#3}::operator()(pybind11::detail::function_call&) const (__closure=__closure@entry=0x0, call=...) at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/pybind11.h:224
224         rec->impl = [](function_call &call) -> handle {
(gdb) s
225             cast_in args_converter;
(gdb) n
227             printf("Constructed cast_in\n");
(gdb) n
230             if (!args_converter.load_args(call)) {
(gdb) n
235             printf("Load args succeeded\n");
(gdb) n
238             process_attributes<Extra...>::precall(call);
(gdb) n
243             printf("Got data pointer\n");
(gdb) n
245             printf("Got capture pointer\n");
(gdb) n
250             printf("Set return value policy\n");
(gdb) n
257                 = cast_out::cast(std::move(args_converter).template call<Return, Guard>(cap->f),
(gdb) s
pybind11::detail::argument_loader<Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&>::call<void, pybind11::detail::void_type, pybind11::cpp_function::cpp_function<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&>(pybind11::cpp_function::cpp_function<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&) && (f=..., this=0x7fffffffd370) at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/cast.h:1417
1417            std::move(*this).template call_impl<remove_cv_t<Return>>(
(gdb) s
pybind11::detail::argument_loader<Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&>::call_impl<void, pybind11::cpp_function::cpp_function<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&, 0ul, 1ul, pybind11::detail::void_type>(pybind11::cpp_function::cpp_function<void, Node<3u>, boost::numeric::ublas::c_vector<double, 3ul> const&, pybind11::name, pybind11::is_method, pybind11::sibling, char [2], pybind11::arg>(void (Node<3u>::*)(boost::numeric::ublas::c_vector<double, 3ul> const&), pybind11::name const&, pybind11::is_method const&, pybind11::sibling const&, char const (&) [2], pybind11::arg const&)::{lambda(Node<3u>*, boost::numeric::ublas::c_vector<double, 3ul> const&)#1}&, std::integer_sequence<unsigned long, 0ul, 1ul>, pybind11::detail::void_type&&) && (f=..., this=0x7fffffffd370) at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/cast.h:1442
1442        Return call_impl(Func &&f, index_sequence<Is...>, Guard &&) && {
(gdb) s
1443            printf("In call_impl\n");
(gdb) n
233       _M_head(_Head_base& __b) noexcept { return __b._M_head_impl; }
(gdb) s
pybind11::detail::cast_op<boost::numeric::ublas::c_vector<double, 3ul> const&> (caster=...) at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/cast.h:48
48      printf("Using cast_op add rvalue reference on arg\n");
(gdb) print caster
$1 = (pybind11::detail::make_caster &&) @0x7fffffffd370: {<pybind11::detail::type_caster_base<boost::numeric::ublas::c_vector<double, 3> >> = {<pybind11::detail::type_caster_generic> = {typeinfo = 0x3, 
      cpptype = 0x3fb999999999999a, value = 0x3fb999999999999a}, static name = {text = "%"}}, <No data fields>}
(gdb) s
printf (__fmt=0x7fffeecda358 "Using cast_op add rvalue reference on arg\n") at /usr/include/x86_64-linux-gnu/bits/stdio2.h:112
112   return __printf_chk (__USE_FORTIFY_LEVEL - 1, __fmt, __va_arg_pack ());
(gdb) n
pybind11::detail::cast_op<boost::numeric::ublas::c_vector<double, 3ul> const&> (caster=...) at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/cast.h:50
50          template cast_op_type<typename std::add_rvalue_reference<T>::type>();
(gdb) s
pybind11::detail::type_caster_base<boost::numeric::ublas::c_vector<double, 3ul> >::operator boost::numeric::ublas::c_vector<double, 3ul>& (this=0x7fffffffd370)
    at /atlas/RSE/CHASTE/PyChaste/dynamic/pybind11/include/pybind11/detail/type_caster_base.h:978
978     operator itype &() {
(gdb) print value
$2 = (void *) 0x3fb999999999999a
(gdb)
kwabenantim commented 1 year ago

Possible Workaround

class PlaneBoundaryCondition2_2Factory
{
public:
  static PlaneBoundaryCondition2_2 create(AbstractCellPopulation<2, 2> * pCellPopulation, py::array_t<double> point_array, py::array_t<double> normal_array)
  {
    c_vector<double, 2> point_c_vector;
    c_vector<double, 2> normal_c_vector;

    // Do conversion from point_array to point_c_vector
    // Do conversion from normal_array to normal_c_vector

    return PlaneBoundaryCondition2_2(pCellPopulation, point_c_vector, normal_c_vector);
  }
};

...

py::class_< PlaneBoundaryCondition2_2 >(m, "PlaneBoundaryCondition2_2")
  .def(py::init(&PlaneBoundaryCondition2_2Factory::create))
kwabenantim commented 1 year ago

Simplified example with typecaster

A simplified example seems to work fine with the current typecaster:

#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <boost/numeric/ublas/vector.hpp>

#include <memory>
#include <string>
#include <vector>

using boost::numeric::ublas::c_vector;

// Typecaster for c_vector<double, 2>
namespace pybind11
{
  namespace detail
  {
    template <>
    struct type_caster<c_vector<double, 2>>
    {
    public:
      typedef c_vector<double, 2> c_vector_double_2_t;
      PYBIND11_TYPE_CASTER(c_vector_double_2_t, const_name("c_vector_double_2_t"));
      bool load(handle src, bool convert)
      {
        if (!convert && !array_t<double>::check_(src)) { return false;  }

        auto buf = array_t<double, array::c_style | array::forcecast>::ensure(src);
        if (!buf) { return false; }
        if (buf.ndim() != 1 or buf.shape()[0] != 2)  { return false;  }

        value.resize(2);
        for (int i = 0; i < 2; i++)
        {
          value[i] = buf.data()[i];
        }

        return true;
      }
      static handle cast(const c_vector<double, 2> &src, return_value_policy policy, handle parent)
      {
        std::vector<size_t> shape(1, 2);
        std::vector<size_t> strides(1, sizeof(double));
        double *data = src.size() ? const_cast<double *>(&src[0]) : static_cast<double *>(NULL);
        array a(std::move(shape), std::move(strides), data);
        return a.release();
      }
    };
  }
}

// Simplified CellPopulation
template <unsigned SPACE_DIM>
class CellPopulation
{
public:
  unsigned getDim() { return SPACE_DIM; }
};

template class CellPopulation<2>;

// Simplified BoundaryCondition
template <unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
class BoundaryCondition
{
public:
  BoundaryCondition(CellPopulation<SPACE_DIM> *population, c_vector<double, 2> point, c_vector<double, 2> normal)
  {
    _total = point[0] + point[1] + normal[0] + normal[1] + population->getDim();
  }

  double getTotal() { return _total; }

private:
  double _total;
};

template class BoundaryCondition<2, 2>;

// Pybind11 Bindings
namespace py = pybind11;
PYBIND11_MODULE(example, m)
{
  py::class_<CellPopulation<2>>(m, "CellPopulation2")
      .def(py::init<>())
      .def("getDim", &CellPopulation<2>::getDim);

  py::class_<BoundaryCondition<2, 2>>(m, "BoundaryCondition2_2")
      .def(py::init<CellPopulation<2> *, c_vector<double, 2>, c_vector<double, 2>>())
      .def("getTotal", &BoundaryCondition<2, 2>::getTotal);
}

Python:

import numpy as np
from example import CellPopulation2, BoundaryCondition2_2

pop = CellPopulation2()
print(pop.getDim())

point = np.array([0.0,0.0])
normal = np.array([0.1,0.0])
bc = BoundaryCondition2_2(pop, point, normal)
print(bc.getTotal())
kwabenantim commented 1 year ago

Passing Tests

kwabenantim commented 1 year ago

TestVertexBasedCellSimulationsPythonTutorial Runs into a segmentation fault.

Error

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

Trace

OffLatticeSimulation::Solve(...);
OffLatticeSimulation::ApplyBoundaries(...);
PlaneBoundaryCondition::ImposeBoundaryCondition(...);
if (dynamic_cast<AbstractOffLatticeCellPopulation<ELEMENT_DIM, SPACE_DIM>*>(this->mpCellPopulation) == nullptr){...}
kwabenantim commented 1 year ago

TestPetscToolsPython Needs wrapping for PETSc Vec and Mat.

Error

======================================================================
ERROR: test_some_vecs (__main__.TestPetscTools)
----------------------------------------------------------------------
TypeError: Unregistered type : _p_Vec

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/vboxuser/PyChaste/test/python/./core/TestPetscToolsPython.py", line 44, in test_some_vecs
    vec = chaste.core.PetscTools.CreateVec(local)
TypeError: Unable to convert function return value to a Python type! The signature was
    (data: List[float]) -> _p_Vec
kwabenantim commented 1 year ago

Closing this issue as the remaining failing tests are unrelated to c_vectors. The remaining problems will be addressed by #15

kwabenantim commented 1 year ago

It appears that the errors in constructors which contained c_vectors was not due to the c_vector typecasters, which seem to work fine everywhere else. Incidentally, a pointer to an AbstractOffLatticeCellPopulation was required in all the constructors where the c_vector typecasters appeared to be failing. The actual problem was with casting a VertexBasedCellPopulation to an AbstractOffLatticeCellPopulation. This was the same problem raised by TestVertexBasedCellSimulationsPythonTutorial and now seems to have been solved. It is not yet clear what changed, but the errors stopped after switching from C++14 to C++17 and fetching latest upstream updates from Chaste develop (fast-forwarded from e1fa6a6 to 3ed690e).