vectorize dispatch std::complex issue

[nbecker]

#include "pybind11/numpy.h"     // vectorize
#include "pybind11/pybind11.h"
#include "pybind11/operators.h"
#include "numpy/arrayobject.h"
#include <complex>

namespace py = pybind11;


PYBIND11_PLUGIN (bug1) {
  if (_import_array() < 0) {
    PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
    return nullptr;
  }

  py::module m("bug1", "pybind11 example plugin");

  m.def("mag_sqr_ufunc", py::vectorize([](double x) { return x * x; }));
  m.def("mag_sqr_ufunc", py::vectorize([](std::complex<double> x) { return real(x) * real(x) + imag(x) * imag(x); }));

  return m.ptr();
}

import numpy as np
from bug1 import mag_sqr_ufunc
u = np.ones (1000000, dtype=complex)
mag_sqr_ufunc(u)

produces:

In [3]: /home/nbecker/.local/bin/ipython3:4: ComplexWarning: Casting complex values to real discards the imaginary part
  import re

This leads me to believe that the wrong ufunc was dispatched, when given complex data type, the scalar overload was erroneously chosen.
This is confirmed by commenting out the scalar overload, in which case there is no message and a correct result is returned.

[nbecker]

Similarly, adding an overload for 'float':

#include "pybind11/numpy.h"     // vectorize
#include "pybind11/pybind11.h"
#include "pybind11/operators.h"
#include "numpy/arrayobject.h"
//#include "ndarray/pybind11.h"
#include <complex>

//namespace nd=ndarray;
namespace py = pybind11;


PYBIND11_PLUGIN (bug1) {
  if (_import_array() < 0) {
    PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
    return nullptr;
  }

  py::module m("bug1", "pybind11 example plugin");

  m.def("mag_sqr_ufunc", py::vectorize([](std::complex<double> x) { return real(x) * real(x) + imag(x) * imag(x); }));
  m.def("mag_sqr_ufunc", py::vectorize([](double x) { return x * x; }));
  m.def("mag_sqr_ufunc", py::vectorize([](float x) { return x * x; }));

  return m.ptr();
}

import numpy as np
from bug1 import mag_sqr_ufunc
u = np.ones (1000000, dtype=complex)
v = mag_sqr_ufunc(u)
w = mag_sqr_ufunc(u.real)
x = mag_sqr_ufunc(u.real.astype(np.float32))

I see x.dtype is

u.real.astype(np.float32).dtype
Out[52]: 
dtype('float32')

In [53]: x.dtype
Out[53]: 
dtype('float64')

which tells me that an array(float) was converted to array(double)
and then dispatched to the (double) overload.

[jagerman]

Pybind's default overload dispatch tries in FIFO order, and stops when one works. That means the float/double overload is never going to work (it's always going to call the first one registered). As to the std::complex one, what happens if you swap the order of the definitions (that'll certainly fix the problem for actual complex, but what I'm unsure of is whether regular python numeric value arguments will also end up there).

[nbecker]

I followed this example:
#563

to do my own dispatch:

  m.def("mag_sqr_ufunc", [](py::array a) {
      if (py::isinstance<py::array_t<complex_t>>(a))
        return py::vectorize([](complex_t x) { return real(x) * real(x) + imag(x) * imag(x); })(a);
      if (py::isinstance<py::array_t<double_t>>(a))
        return py::vectorize([](double_t x) { return x * x; })(a);
      if (py::isinstance<py::array_t<complex64_t>>(a))
        return py::vectorize([](complex64_t x) { return real(x) * real(x) + imag(x) * imag(x); })(a);
      if (py::isinstance<py::array_t<float>>(a))
        return py::vectorize([](float x) { return x * x; })(a);
      else
        throw py::type_error("mag_sqr unhandled type");
    });

I believe it is working

[nbecker]

I think it's unfortunate that function overloads of
f(py::array_t)
will match when T is not matched exactly. Glancing at the code, I'm thinking it's actually a conversion to array(no _t) which is matched, with

   PYBIND11_OBJECT_CVT(array, buffer, detail::npy_api::get().PyArray_Check_, raw_array)

I think it would be better if f(py::array) would allow conversion, while f(py::array_t) would not, although I don't know how that could be implemented.

[wjakob]

This is a long-standing topic of discussions, see #338 and #634.

[jagerman]

@nbecker - can you confirm whether this works as expected with current master? #634 and #643 should have taken care of it.

[nbecker]

It seems to be, but it reveals a small bug in printing
Traceback (most recent call last):
File "./test_mag_sqr.py", line 9, in
y = mag_sqr(2.0)
TypeError: mag_sqr(): incompatible function arguments. The following argument types are supported:
1. (in: numpy.ndarray[float]) -> object << this was py::array_t
2. (in: numpy.ndarray[float]) -> object << this was py::array_t
3. (in: numpy.ndarray[complex]) -> object << this was py::array_t<std::complex>
4. (in: numpy.ndarray[complex]) -> object << this was py::array_t<std::complex>

[dean0x7d]

I believe the signature there is intentionally a generic object because calling the function with a single value return a scalar, not an array.

[nbecker]

Sorry, that message wasn't rendered correctly, what I meant to say was

Traceback (most recent call last):
File "./test_mag_sqr.py", line 9, in 
y = mag_sqr(2.0)
TypeError: mag_sqr(): incompatible function arguments. The following argument types are supported:
1. (in: numpy.ndarray[float]) -> object << this was py::array_t<double>
2. (in: numpy.ndarray[float]) -> object << this was py::array_t<float>
3. (in: numpy.ndarray[complex]) -> object << this was py::array_t<std::complex<double>>
4. (in: numpy.ndarray[complex]) -> object << this was py::array_t<std::complex<float>>

The issue is that messages 1 and 2 both say 'float', while the actual c++ types are 'double' and 'float', and similarly for messages 3 and 4.

[dean0x7d]

Ah, sorry, I thought it was about the return type. The arguments types are indeed confusing. I guess the best thing would be to have it say float32, float64, complex64 and complex128 in keeping with numpy dtype naming.

[dean0x7d]

I think that between #634, #643 and #724 this should be resolved now.