std::complex<T> not resolved correctly

[nils-werner]

Consider the following example

import numpy
import ctest

print ctest.square(numpy.arange(3, dtype=float))
print ctest.square(numpy.arange(3, dtype=float) * 1j)

print ctest.square(9.0)
print ctest.square(9.0 * 1j)

where ctest is my pybind11 module

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

namespace py = pybind11;

template <typename T>
T square(T x) {
    return x * x;
}

void pyexport(py::module& m) {
    m.def("square", py::vectorize(square<std::complex<double>>));
    m.def("square", py::vectorize(square<std::complex<float>>));
    m.def("square", py::vectorize(square<float>));
    m.def("square", py::vectorize(square<double>));
}

The problem is that my functions are always resolved to the std::complex type, the above Python code returning

[ 0.+0.j  1.+0.j  4.+0.j]
[ 0.+0.j -1.+0.j -4.+0.j]
(81+0j)
(-81+0j)

1. should be real valued
2. is OK
3. should be real valued
4. is OK

If I now reorder my defs like

void pyexport(py::module& m) {
    m.def("square", py::vectorize(square<float>));
    m.def("square", py::vectorize(square<double>));
    m.def("square", py::vectorize(square<std::complex<double>>));
    m.def("square", py::vectorize(square<std::complex<float>>));
}

the std::complex type is never resolved:

[ 0.  1.  4.]
test.py:8: ComplexWarning: Casting complex values to real discards the imaginary part
  print ctest.square(numpy.arange(3, dtype=float) * 1j)
[ 0.  0.  0.]
81.0
Traceback (most recent call last):
  File "test.py", line 11, in <module>
    print ctest.square(9.0 * 1j)
TypeError: Incompatible function arguments. The following argument types are supported:
    1. (array[float]) -> object
    2. (array[float]) -> object
    3. (array[complex]) -> object
    4. (array[complex]) -> object

1. is OK
2. is implicitly cast back to real valued number (not OK)
3. is Ok
4. fails to resolve

[wjakob]

First of all: there is no way to distinguish between float and single precision in Python. The first overload that you specify will always be used.

[nils-werner]

I also noticed that

print type(ctest.square(9))
print type(ctest.square(9.0))

and

void pyexport(py::module& m) {
    m.def("square", py::vectorize(square<int>));
    m.def("square", py::vectorize(square<float>));
    m.def("square", py::vectorize(square<double>));
}

will always yield long.

Did I miss that overloading m.def()s are not possible?

[wjakob]

Overloading is fine, problems appear when one type is implicitly convertible to another type. Since your use case is fairly advanced, I fear that you are on your own -- you'll have to create a single overload which takes a py::object and then write the code yourself that invokes the various specializations after inspecting the contents using the Python C API.

[wjakob]

You could also accept a py::buffer and then inspect the buffer's "format" flag -- maybe that is the easiest. Afterwards, you could simply call py::vectorize(square)(array_t(object_arg)) if it's an int, etc.

[nils-werner]

Implicit conversion makes sense, that still doesn't explain why in some cases the correct type cannot be resolved at all.

Also, why would it pick square<int> for doubles even when square<double> is defined?

[wjakob]

After thinking a bit more about it, I suppose that these sorts of overloads would be useful to have without the forceful casts. This commit should address your problems: b47a9de -- see the selective_func examples. You will then have to call py::vectorize from that.

[wjakob]

(All of this was related to the forcecast argument that was provided to NumPy by deafult. Can you confirm that it works for you with this change?)

[nils-werner]

Not yet, I am still struggling with making it compile when pip install -e installing pybind and using cppimport. :-)

[nils-werner]

Yes, that indeed fixes it. Order of definition still matters though.

[wjakob]

Yes, that always matters in pybind11. It tries one overload after the other until one works.

[nils-werner]

For future reference

import cppimport
cppimport.set_quiet(False)

import numpy
import ctest

# builtin types
for v in [
    9,
    9 * 1j,
    9.0,
    9.0 * 1j,
]:
    print type(v), type(ctest.square(v))

# numpy types
for v in [
    numpy.arange(3),
    numpy.arange(3, dtype=float),
    numpy.arange(3) * 1j,
    numpy.arange(3, dtype=float) * 1j,
]:
    print v.dtype, ctest.square(v).dtype

and

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

namespace py = pybind11;

template <typename T>
T square(T x) {
    return x * x;
}

void pyexport(py::module& m) {
    m.def("square", square<long>);
    m.def("square", square<double>);
    m.def("square", [](py::array_t<long, py::array::c_style> x) { return py::vectorize(square<long>)(x); });
    m.def("square", [](py::array_t<double, py::array::c_style> x) { return py::vectorize(square<double>)(x); });
    m.def("square", [](py::array_t<std::complex<double>, py::array::c_style> x) { return py::vectorize(square<std::complex<double>>)(x); });
}

now yield

<type 'int'> <type 'long'>
<type 'complex'> <type 'complex'>
<type 'float'> <type 'float'>
<type 'complex'> <type 'complex'>
int64 int64
float64 float64
complex128 complex128
complex128 complex128

[aldanor]

@wjakob Would you consider reverting the forcecast default (or changing it to c_style)? This is an "implicit versus explicit" kind of issue -- if your function accepts arrays of ints, it probably wants just that. If it expects a structured type and you give it a bunch of ints, it's likely a user error and you don't want those values to be copied over all fields. Also, (as with the example in this issue) the overloads may not work as you would expect them to unless you disable the forcecast; I believe it's quite counter-intuitive from C++ point of view.

It gets even more weird once you start working with structured arrays. NumPy will happily cast anything to anything according to a set of strange conventions. Chances are, you'll never want to see any of this happen accidentally, and in a rare case you do, you will knowingly and explicitly provide a forcecast flag. Just a few examples:

# cast ints to structs
>>> a = np.array([1, 2, 3], dtype='u2')
>>> a.astype([('a', 'u4'), ('b', 'u1')], casting='unsafe')
array([(1, 1), (2, 2), (3, 3)], 
      dtype=[('a', '<u4'), ('b', 'u1')])

# happily cast integer structs to floats, discarding half of the data
>>> np.rec.fromarrays([[1, 2, 3], [4, 5, 6]]).astype('f4', casting='unsafe')
rec.array([ 1.,  2.,  3.], 
          dtype=float32)

# strings to bools, no problem
>>> np.array(['-1', '0', '1']).astype('bool', casting='unsafe')
array([ True, False,  True], dtype=bool)

It's also the default NumPy behaviour for most routines to not use forcecast and only do safe casting:

Unless FORCECAST is present in flags, this call will generate an error if the data type cannot be safely obtained from the object.

[wjakob]

Yes, I am open to this. It will go into the next major version.

[wjakob]

Actually, let me take a step back before committing to this yet. The only reason why the forcecast convention was there is that NumPy will by default not convert a list like [0, 0, 1] to a 3D float vector. Do you see any way to accomplish this while avoiding over-eager casts?