[Numpy-discussion] Best way to broadcast a function from C
Jaime Fernández del Río
jaime.frio at gmail.com
Fri Aug 22 20:50:26 EDT 2014
You can always write your own gufunc with signature '(),(),()->(a, a)', and
write a Python wrapper that always call it with an `out=` parameter of
shape (..., 3, 3), something along the lines of:
def my_wrapper(a, b, c, out=None):
if out is None:
out = np.empty(np.broadcast(a,b,c).shape + (3, 3))
if out.shape[-2:] != (3, 3):
raise ValueError("Wrong shape for 'out'")
return my_gufunc(a, b, c, out=out)
Writing your own gufunc is a little challenging, but not that hard with all
the examples now available in the numpy.linalg code base, and it is a
terribly powerful tool.
On Fri, Aug 22, 2014 at 4:40 PM, James Crist <crist042 at umn.edu> wrote:
> I suspected as much. This is actually part of my work on numerical
> evaluation in SymPy. In its current state compilation to C and autowrapping
> *works*, but I think it could definitely be more versatile/efficient. Since
> numpy seemed to have solved the broadcasting/datatype issues internally I
> hoped I could reuse this.
> I'll look into nditer (as this is for a library, it's better to do it
> *right* than do it now), but right now it's looking like I may need to
> implement this functionality myself...
> On Thu, Aug 21, 2014 at 4:50 PM, Nathaniel Smith <njs at pobox.com> wrote:
>> 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.
>> > 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
>> > 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
>> > call isn't a scalar but an array, but the effect I'm looking for is
>> > 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
>> > 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.
>> > read the "writing your own ufunc" part of the docs, but am unsure if
>> > 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
>> 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.
>> Nathaniel J. Smith
>> Postdoctoral researcher - Informatics - University of Edinburgh
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion