[MATRIX-SIG] Making pseudo-ufuncs?

[Johann Hibschman]

Hi all,

I've been trying to write a class that wraps a normal unary or binary
function to allow it to be simply (if slowly) applied to arrays.
Hopefully, it would increase the readibility of some code.

First, is there a better way to do this?

Second, is the method I use to detect the types of the arrays a good idea? 
It would seem to fall apart for special Matrix and Vector classes.

Third, is there a good way to implement outer, accumulate, and reduce for
binary functions?

So far, I have (FuncOps is my function operations class):

class FuncBinder(FuncOps):
    def __init__(self, a_f):
	if ((type(a_f) == UfuncType) or
	    (type(a_f) == InstanceType and
	     FuncOps in a_f.__class__.__bases__)):
	    self.__call__ = a_f        # overwrite the existing method
	self.f = a_f

    def __call__(self, arg):
	"Default call routine, used for ordinary functions."
	if type(arg) == ArrayType:
	    return array_map(self.f, arg)
	else:
	    return self.f(arg)

where

def array_map(f, ar):
    "Apply an ordinary function to all values in an array."
    flat_ar = ravel(ar)
    out = zeros(len(flat_ar), typecode=flat_ar.typecode())
    for i in xrange(len(flat_ar)):
	out[i] = f(flat_ar[i])
    out.shape = ar.shape
    return out

for unary functions, and (for binary functions), I've defined methods:

def reduce(self, a, axis=0):
    result = take(a, [0], axis)
    for i in range(1, a.shape[axis]):
	result = self(result, take(a, [i], axis))
    return result

def accumulate(self, a, axis=0):
    n = len(a.shape)
    sum = take(a, [0], axis)
    out = zeros(a.shape, a.typecode())
    for i in range(1, a.shape[axis]):
	out[all_but_axis(i, axis, n)] = self(sum, take(a, [i], axis))
    return out

def all_but_axis(i, axis, num_axes):
    """
    Return a slice covering all combinations with coordinate i along
    axis.  (Effectively the hyperplane perpendicular to axis at i.)
    """
    the_slice  = ()
    for j in range(num_axes):
	if j == axis:
	    the_slice = the_slice + (i,)
	else:
	    the_slice = the_slice + (slice(None),)
    return the_slice

def outer(op, a, b):
    n_a = len(a.shape)
    n_b = len(b.shape)
    a2  = reshape(a, a.shape + (1,)*n_b)
    b2  = reshape(b, (1,)*n_a + b.shape)
    
    # duplicate each array in the appropriate directions
    a3  = a2
    for i in range(n_b):
	a3 = repeat(a3, (b.shape[i],), n_a+i)
    b3 = b2
    for i in range(n_a):
	b3 = repeat(b3, (a.shape[i],), i)

    answer = array_map_2(op, a3, b3)
    return answer

These all seem over-complex to me.  Is there are more straightforward way
to do this?

---
Johann A. Hibschman         | Grad student in Physics, working in Astronomy.
johann@physics.berkeley.edu | Probing pulsar pair production processes.


_______________
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
_______________