[Numpy-discussion] Best way to broadcast a function from C

Nathaniel Smith njs at pobox.com
Thu Aug 21 17:50:03 EDT 2014


On Thu, Aug 21, 2014 at 2:34 AM, James Crist <crist042 at umn.edu> wrote:
> All,
>
> I have a C function func that takes in scalar arguments, and an array of
> fixed dimension that is modified in place to provide the output. The
> prototype is something like:
>
> `void func(double a, double b, double c, double *arr);`
>
> I've wrapped this in Cython and called it from python with no problem. What
> I'd like to do now though is get it set up to broadcast over the input
> arguments and return a 3 dimensional array of the results. By this I mean
>
> a = array([1, 2, 3])
> b = array([2.0, 3.0, 4.0])
> c = array([3, 4, 5])
>
> func(a, b, c) -> a 3d array containing the results of func for (a, b, c) =
> (1, 2.0, 3), (2, 3.0, 4), (3, 4.0, 5)
>
> I'm not sure if this would qualify as a ufunc, as the result of one function
> call isn't a scalar but an array, but the effect I'm looking for is similar.
> Ideally it would handle datatype conversions (in the above `a` and `c`
> aren't double, but `func` takes in double). It would also be awesome to
> allow an argument to be a scalar and not an array, and have it be broadcast
> as if it were.
>
> I'm just wondering what the best way for me to hook my code up to the
> internals of numpy and get this kind of behavior in an efficient way. I've
> read the "writing your own ufunc" part of the docs, but am unsure if what
> I'm looking for qualifies. Note that I can change the inner workings of
> `func` if this is required to achieve this behavior.

I don't think it's currently possible to write a ufunc that maps
scalars to fixed-length arrays.

There's probably some not-too-terrible way to do this using the C
nditer interface, but IIRC the docs aren't for that aren't very
helpful.

The simplest approach is just to do a bit of work at the Python level
using the standard numpy API to do error-checking and set up your
arrays in the way you want, and then pass them to the C
implementation. This is not the best way to get a Right solution that
will handle all the funky corner cases that ufuncs handle, but it's by
far the fastest way to get something that's good enough for whatever
you need Right Now.

-n

-- 
Nathaniel J. Smith
Postdoctoral researcher - Informatics - University of Edinburgh
http://vorpus.org



More information about the NumPy-Discussion mailing list