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.
Jaime
On Fri, Aug 22, 2014 at 4:40 PM, James Crist <crist...@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...
>
> Thanks,
>
> -Jim
>
>
> On Thu, Aug 21, 2014 at 4:50 PM, Nathaniel Smith <n...@pobox.com> wrote:
>
>> On Thu, Aug 21, 2014 at 2:34 AM, James Crist <crist...@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
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
--
(\__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
