set up proper communication with numpy

[matsen]

set up proper communication with numpy

[matsen]

I think the PYBIND11_MAKE_OPAQUE(std::vector<double>); and https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html approach is going to work, but in case it doesn't here are some code snippets:

From pybind/pybind11#1042:

pybind11::array WrapVector(std::vector<double> *v) { 
  auto capsule = py::capsule(
      &v, [](void *v) { delete reinterpret_cast<std::vector<double> *>(v); });
  return py::array(static_cast<pybind11::ssize_t>(v->size()), v->data(),
                   capsule);
}

which led to double-free problems.

This silliness also worked:

class Vector : public std::vector<double> {
 public:
  Vector(size_t size) : std::vector<double>(size){};
};

[matsen]

Here's the code I used to time the various versions:

import numpy as np
import sbn
import timeit

v1_c = sbn.make_vector()
v2_c = sbn.make_vector()
v3_c = sbn.make_vector()
a1_c = np.array(v1_c, copy=False)
a2_c = np.array(v2_c, copy=False)
a3_c = np.array(v3_c, copy=False)

print("pybind11 version: ", timeit.timeit("a3_c = a1_c + a2_c", globals=globals()))

a1 = np.random.uniform(size=len(a1_c))
a2 = np.random.uniform(size=len(a1_c))
a3 = np.random.uniform(size=len(a1_c))

print("numpy version: ", timeit.timeit("a3 = a1 + a2", globals=globals()))

class Vector {
 public:
  Vector(size_t size) : m_size(size) { m_data = new double[size]; }
  double *data() { return m_data; }
  size_t size() const { return m_size; }

 private:
  size_t m_size;
  double *m_data;
};

Vector make_vector() {
  std::random_device
      rd;  // Will be used to obtain a seed for the random number engine
  std::mt19937 gen(rd());  // Standard mersenne_twister_engine seeded with rd()
  std::uniform_real_distribution<> dis(0., 1.);
  Vector x(1000);
  for (size_t i = 0; i < x.size(); i++) {
    x.data()[i] = dis(gen);
  }
  return x;
}
...
  py::class_<Vector>(m, "Vector", py::buffer_protocol())
      .def_buffer([](Vector &v) -> py::buffer_info {
        return py::buffer_info(
            v.data(),                                /* Pointer to buffer */
            sizeof(double),                          /* Size of one scalar */
            py::format_descriptor<double>::format(), /* Python
                                                       struct-style format
                                                       descriptor */
            1,                                       /* Number of dimensions */
            {v.size()},                              /* Buffer dimensions */
            {sizeof(double)});
      });
  m.def("make_vector", &make_vector, "test");