[Numpy-discussion] Speeding up numarray -- questions on its design

[David M. Cooke]

Robert Kern <rkern at ucsd.edu> writes:

> Travis Oliphant wrote:
>
>>>> 4) Object arrays must be supported.  This was a bad oversight and an
>>>> important feature of Numeric arrays.
>>>>
>>> The current implementation does support them (though in a different
>>> way, and generally not as efficiently, though Todd is more up on the
>>> details here). What aspect of object arrays are you finding lacking?
>>> C-api?
>> I did not see such support when I looked at it, but given the
>> previous comment, I could easily have missed where that support is
>> provided.  I'm mainly following up on Konrad's comment that his
>> Automatic differentiation does not work with Numarray because of the
>> missing support for object arrays.  There are other applications for
>> object arrays as well.   Most of the support needs to come from the
>> ufunc side.
>
> It's tucked away in numarray.objects. Unfortunately for Konrad's
> application, numarray ufuncs don't recognize that it's being passed an
> object with the special methods defined, and they won't automatically
> create 0-D object "arrays". 0-D object arrays will work just fine when
> using operators (x+y works), but not when explicitly calling the
> ufuncs (add(x,y) does not work). Both methods work fine for 0-D
> numerical arrays.

Are the 0-D object arrays necessary for this? The behaviour that
Konrad needs is this (highly abstracted):

class A:
    def __add__(self, other):
        return 0.1
    def sin(self):
        return 0.5

Then:
>>> a = A()
>>> a + a
0.10000000000000001
>>> Numeric.add(a,a)
0.10000000000000001
>>> Numeric.sin(a)
0.5

The Numeric ufuncs, if the argument isn't an array, look for a method
of the right name (here, sin) on the object, and call that.

You could define a delegate class that does this with something like

class MathFunctionDelegate:
    def __init__(self, fallback=Numeric):
        self._fallback = fallback
    def add(self, a, b):
        try:
            return a + b
        except TypeError:
            return self._fallback.add(a, b)
    def sin(self, x):
        sin = getattr(x, 'sin', None)
        if sin is None:
            return self._fallback.sin(x)
        else:
            return sin(x)

    ... etc. ...

(This could be a module, too. This just allows parameterisation.)

In ScientificPython, FirstDerivatives.py has a method of the DerivVar
class that looks like this:

    def sin(self):
        v = Numeric.sin(self.value)
        d = Numeric.cos(self.value)
        return DerivVar(v, map(lambda x,f=d: f*x, self.deriv))

Add something like this to the __init__:
        self._mathfuncs = MathFunctionDelegate(Numeric)
and that sin method becomes
    def sin(self):
        v = self._mathfuncs.sin(self.value)
        d = self._mathfuncs.cos(self.value)
        return DerivVar(v, map(lambda x,f=d: f*x, self.deriv))

That's not quite perfect, as the user has to use a mathfuncs object
also; that's why having Numeric or numarray do the delegation
automatically is nice.

This would work equally well with numarray (or the math or cmath
modules!) replacing Numeric. You could get fancy and be polymorphic:
choose the right module to use depending on the type of the argument
(Numeric arrays use Numeric, floats use math, etc.).

If this was a module instead, you could have registration of types.
I'll call this module numpy. Here's a possible (low-level) usage:

import numpy
import Numeric, numarray, math, cmath
from Scientific.Functions import Derivatives
numpy.register_type(Numeric.arraytype, Numeric)
numpy.register_type(numarray.NumArray, numarray)
numpy.register_type(float, math)
numpy.register_type(complex, cmath)
numpy.register_type(Derivatives.DerivVar, Derivates.derivate_math)
numpy.default_constructor(numarray.array)

a = numpy.array([1,2,3])        # makes a numarray
b = Numeric.array([1,2,3])      # Numeric array

print numpy.sin(a), numpy.sin(b)


Things to consider with this would be:
* how to handle a + b
* where should the registering of types be done? (Probably by the
  packages themselves)
* more complex predicates for registering handlers? (to handle
  subclasses, etc.)

etc.

Ok, I hope that's not too rambling. But the idea is that neither
Numeric nor numarray need to provide the delegation ability.

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|cookedm at physics.mcmaster.ca