Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ublas vector wrapping support is broken #11

Closed
kwabenantim opened this issue Jun 8, 2023 · 11 comments · Fixed by #18
Closed

Ublas vector wrapping support is broken #11

kwabenantim opened this issue Jun 8, 2023 · 11 comments · Fixed by #18
Assignees
Labels
bug Something isn't working

Comments

@kwabenantim
Copy link
Member

kwabenantim commented Jun 8, 2023

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 kwabenantim self-assigned this Jun 8, 2023
@kwabenantim kwabenantim added the bug Something isn't working label Jun 8, 2023
@kwabenantim
Copy link
Member Author

Also add support for zero_vector and unit_vector:

  • Cylindrical2dNodesOnlyMesh::SetUpBoxCollection has a unit_vector argument
  • PeriodicNodesOnlyMesh::SetUpBoxCollection has a zero_vector argument

@kwabenantim kwabenantim changed the title Fix c_vector support Ublas vector wrapping support is broken Jun 17, 2023
@MILeach
Copy link
Collaborator

MILeach commented Jun 20, 2023

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.

This was referenced Jun 20, 2023
@kwabenantim kwabenantim linked a pull request Jun 21, 2023 that will close this issue
@MILeach
Copy link
Collaborator

MILeach commented Jun 21, 2023

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
Copy link
Collaborator

MILeach commented Jun 27, 2023

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
Copy link
Member Author

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
Copy link
Member Author

kwabenantim commented Jun 30, 2023

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
Copy link
Member Author

kwabenantim commented Jul 11, 2023

Passing Tests

  • ode/TestOdeSystemPython.py
  • cell_based/tutorials/TestScratchAssayTutorial.py
  • cell_based/tutorials/TestMeshBasedCellSimulationsPythonTutorial.py
  • cell_based/tutorials/TestPottsBasedCellSimulationsPythonTutorial.py
  • cell_based/tutorials/TestSpheroidTutorial.py
  • cell_based/tutorials/TestCellSortingTutorial.py
  • cell_based/TestVertexBasedCellPopulationPython.py
  • cell_based/TestMeshBasedCellPopulationPython.py
  • cell_based/TestCaBasedCellPopulationPython.py
  • mesh/TestPottsMeshPython.py
  • core/TestFileFinderPython.py
  • core/TestRandomNumberGeneratorPython.py
  • core/TestTimerPython.py
  • core/TestVersionPython.py
  • core/TestOutputFileHandlerPython.py
  • cell_based/tutorials/TestVertexBasedCellSimulationsPythonTutorial.py
  • cell_based/tutorials/TestTensileTestTutorial.py
  • cell_based/tutorials/TestNodeBasedCellSimulationsPythonTutorial.py
  • cell_based/TestNodeBasedCellPopulationPython.py
  • mesh/TestChastePointPython.py
  • core/TestPetscToolsPython.py

@kwabenantim
Copy link
Member Author

kwabenantim commented Jul 12, 2023

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
Copy link
Member Author

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
Copy link
Member Author

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

@kwabenantim
Copy link
Member Author

kwabenantim commented Jul 18, 2023

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).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants