vectorizing a method that operates on an array

[nschloe]

I have a pybind11ified function that operates on a vector of dimension 1,

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

namespace py = pybind11;

void distill(py::array_t<double, py::array::c_style | py::array::forcecast> p) {
  auto r = p.mutable_unchecked<1>();
  for (ssize_t i = 1; i < r.shape(0); i++) {
    r(i) = 3.14 * (r(i) + r(i-1));  // whatever
  }
}

PYBIND11_MODULE(_accupy, m) {
  m.def("distill", &distill);
}

I would now like to change the function to have it operate on the first dimension of numpy arrays of arbitrary dimensionality. The first thing that comes to mind is adding an inner loop

for (ssize_t j = 0; j < r.shape(1); j++) {
}

to the above and wrapping the method in Python code à la

def distill(p):
    q = p.reshape(p.shape[0], numpy.prod(p.shape[1:]))
    out = _mymod.distill(q)
    return out.reshape(p.shape)

This however wouldn't work with arrays of one dimension anymore. Also, instead of for (ssize_t j = 0; j < r.shape(1); j++) {} it's probably better to use BLAS's multiply-add functions there; a rabbit hole I'd rather not descend into.

Is there a canonical way for vectorizing pybind11 methods that operate on an array?

[calumhrmurray]

Is this what you are looking for?

http://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html#vectorizing-functions

For me it was as simple as,

.def("function_name",py::vectorize(&class_name::function_name))

[nschloe]

This only works for functions that natively take scalar values, doesn't it?

[calumhrmurray]

I think I misunderstood your problem. You want to give an array to a vectorized c++ function? Whereas my solution was giving an array to a vectorized pybind11 c++ function?

Why not just have a double function in c++ and then vectorize it?

In my experience a py::vectorize tool works very quickly. In my particular use a comparison of vectorized c++ code versus vectorizing with pybind11, the pybind11 vectorization came out quicker. This could be that my c++ was not great though.

[nschloe]

The operation is supposed to work only in the first dimension of the array. Check out the example in the original post:

  for (ssize_t i = 1; i < r.shape(0); i++) {
    r(i) = 3.14 * (r(i) + r(i-1));
  }

[EricCousineau-TRI]

Just to make sure, are you aiming towards this functionality?
https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.ufunc.reduce.html

One option (albeit more risky) is to consider implementing your functions as a UFunc, and leverage NumPy's mechanisms, like numpy.ufunc.reduce.

I have a WIP branch here which does some basic ufunc stuff, though for me it's namely for enabling user-defined dtypes:
https://github.com/EricCousineau-TRI/pybind11/blob/31fe60a00ca6cbdd83bc28044267022d8301d6a7/tests/test_numpy_dtype_user.cpp#L222-L232
https://github.com/EricCousineau-TRI/pybind11/blob/31fe60a00ca6cbdd83bc28044267022d8301d6a7/tests/test_numpy_dtype_user.py#L79
https://github.com/EricCousineau-TRI/pybind11/blob/31fe60a00ca6cbdd83bc28044267022d8301d6a7/include/pybind11/detail/numpy_ufunc.h#L111

I can try to upstream it a little sooner if you think it'd be useful.

[mgeier]

@nschloe I'm not quite sure if I understood your problem correctly, but I think if you want to do something along the first dimension (as opposed to along an arbitrary user-specified dimension), you can change your function to take a py::iterable instead of a py::array_t<...>.
Then you can use .begin() and .end() to iterate over arrays with one less dimension.

In my case I also added a custom cast_iterator to the mix (see #360 and #462, which was sadly closed), that converts the sub-arrays to a custom C++ type.