[Numpy-discussion] using numpy functions on an array of objects

[Sebastian Walter]

On Sun, Feb 1, 2009 at 12:24 AM, Robert Kern <robert.kern at gmail.com> wrote:
> On Sat, Jan 31, 2009 at 10:30, Sebastian Walter
> <sebastian.walter at gmail.com> wrote:
>> Wouldn't it be nice to have numpy a little more generic?
>> All that would be needed was a little check of the arguments.
>>
>> If I do:
>> numpy.trace(4)
>> shouldn't numpy be smart enough to regard the 4 as a 1x1 array?
>
> Why? It's not a 1x1 array. It's a scalar. If you want a 1x1 array,
> give it a 1x1 array.
>
>> numpy.sin(4) works!
>
> Yes, numpy.sin() operates on scalars in addition to arrays.
>
>> and if
>> x = my_class(4)
>>
>> wouldn't it be nice if
>>
>> numpy.trace(x)
>> would call
>> x.trace() ?
>>
>> numpy.sin(my_class(4)) works!
>>
>> Wouldn't it be nice if numpy worked a little more consistent.
>> Is this worth a ticket? Or am I missing something here?
>
> numpy.sin() is a ufunc. Unary ufuncs will call the method of the same
> name on objects in an object array (or the scalar itself if given an
> object scalar). For example:
>
> In [8]: class MyClass(object):
>   ...:     def __init__(self, x):
>   ...:         self.x = x
>   ...:     def __repr__(self):
>   ...:         return 'MyClass(%r)' % (self.x,)
>   ...:     def sin(self):
>   ...:         return MyClass(self.x+1)
>   ...:
>   ...:
>
> In [9]: sin(MyClass(4))
> Out[9]: MyClass(5)
>
> In [10]: sin([MyClass(4), MyClass(5)])
> Out[10]: array([MyClass(5), MyClass(6)], dtype=object)
>
>
> You'll notice that numpy.sin() does not try to call the list.sin()
> method when given the list. It interprets it as an object array, and
> calls the MyClass.sin() method on each of the elements.
>
> numpy.trace() is not an unary ufunc. It's just a function that
> operates on (N>=2)-D arrays. You simply couldn't apply the same rules
> as numpy.sin(). Otherwise, it would try to call the .trace() method on
> each of the objects in your container, and obviously you can't
> implement trace that way.
>
> Having numpy.trace(x) simply call x.trace() would not be making numpy
> more consistent.
>
> Now, that said, the implementation of numpy.trace(x, *args) is
> actually simply asarray(x).trace(*args). That should probably be
> asanyarray(x) in order to allow ndarray subclasses. But this only
> works because ndarray.trace() already exists. Making every function in
> numpy check for a method first is just not going to happen.

Ok I see. I understand your reasoning.
Nonetheless, I didn't suggest that trace() and sin() are the same,
because they are not. I just wanted to express that they should act
the same
if the object is of *unknown type*.

I mean, numpy.sin(MyClass(3)) works.
In the worst of all possible worlds, numpy would raise an exception
because MyClass is *not an array or a scalar*.
But it doesn't. And that is really cool! It's awesome that numpy works
on arbitrary type for the most part.

In contrast, if trace(X) encounters an unknown type it simply raises
an exception. It could as well try *in the very end*
to call X.trace().
I.e. *not* "numpy check for a method first" but *numpy check for a
method as last resort*.

That wouldn't do any harm, would it? And it is not a major effort to
add those simple checks to dot(), trace(), inv() .
I could provide a patch.

But if that is deemed "not going to happen" for whatever reasons. Is
there a good workaround?
I.e. if I do
import numpy
import mypackage

can I overwrite the functions of numpy? I mean, that is quite a hack,
but is the next *best* option.

The reason for that need is, that I am writing a Python module to
compute higher order derivatives of functions that are given as
an algorithm on matrix operations:

e.g. we want to compute the Hessian of the function

def f(X,Y,Z):
    """ X is (N,M) array, Y is (M,K) array, Z is  (K, L) array"""
     V = numpy.dot(X,Y)
     W = numpy.dot(Y,Z)
     return numpy.trace(V*W)

To do that, I generalized the real numbers to the field of truncated
Taylor polynomials of scalars and
real matrices to truncated Taylor polynomials of matrices.
The theory is explained on
http://en.wikipedia.org/wiki/Automatic_differentiation

You can have a look at the unit test, e.g. at def test_2x2dot2x2_reverse() on
http://github.com/b45ch1/algopy/blob/bd7154e2e7a7e6e1931addc0d9ec0604d488d73f/unit_tests/matrix_reverse.py
Mtc is a class of Matrix Taylor Coefficients. The class Function is
used to build the computational graph of function nodes.

I think it would be a nice addition to everyone who is doing
scientific programming with numpy because derivatives are often
required
and divided differences suck for anything but first order on small
algorithms and especially if one wants to differentiate solutions of
ODEs or PDEs.
If trace, dot, inv, etc. don't work that way, ppl would have to define
two versions of the function f to make it work.


best regards,
Sebastian Walter






>
> --
> Robert Kern
>
> "I have come to believe that the whole world is an enigma, a harmless
> enigma that is made terrible by our own mad attempt to interpret it as
> though it had an underlying truth."
>  -- Umberto Eco
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>