This thread discusses one of the things highest on my wishlist for numpy. I have attached my first attempt at a code that will create a broadcastingfunction that will broadcast a function over arrays where the last N indices are assumed to be for use by the function, N=2 would be used for matrices. It is implemented for Numeric. /Jörgen Pau Gargallo skrev:
Pau, can you confirm that this is the same as the routine you're interested in?
def dot2(a,b): '''Returns dot product of last two dimensions of two 3-D arrays, threaded over first dimension.''' try: assert a.shape[1] == b.shape[2] assert a.shape[0] == b.shape[0] except AssertionError: print "incorrect input shapes" res = zeros( (a.shape[0], a.shape[1], a.shape[1]), dtype=float ) for i in range(a.shape[0]): res[i,...] = dot( a[i,...], b[i,...] ) return res
yes, that is what I would like. I would like it even with more dimensions and with all the broadcasting rules ;-) These can probably be achieved by building actual 'arrays of matrices' (an array of matrix objects) and then using the ufunc machinery. But I think that a simple dot2 function (with an other name of course) will still very useful.
pau
------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=k&kid0709&bid&3057&dat1642 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
# -*- coding: ISO-8859-1 -*- from exceptions import ArithmeticError import Numeric,LinearAlgebra,pdb from Numeric import array,NewAxis,repeat,zeros from LinearAlgebra import inverse import operator def transpose(array): """transpose list of list """ return Numeric.transpose(array).tolist() def outer(data): result=[[]] for l in data: newresult=[] for x in result: for e in l: newresult.append(x+[e]) result=newresult return result def is_squaremultiarray(arraymatrix): try: arraymatrix.shape except AttributeError: return False if len(arraymatrix.shape)>=2: return arraymatrix.shape[-1]==arraymatrix.shape[-2] else: return False def is_multiarray(arraymatrix): try: arraymatrix.shape except AttributeError: return False if len(arraymatrix.shape)>=2: return True else: return False def broad(x): """broad(x) will pick out the highest value of the list x. There can be one or two distinct values in x one of which must be one >>> broad([1,3,3]) 3 >>> broad([1]) 1 >>> broad([1]) 1 >>> broad([]) Traceback (most recent call last): ... IndexError: Incompatible indices for broadcasting >>> broad([1,2,3]) Traceback (most recent call last): ... IndexError: Incompatible indices for broadcasting """ if len(x)==1: return x[0] elif x==[1]: return 1 q=set(x) if q==set([1]): return 1 q.discard(1) if len(q)==1: return q.pop() else: raise IndexError("Incompatible indices for broadcasting") def make_index(x,broad): if x==1: return [0 for x in range(broad)] else: return range(x) def gen_iter_idx(shapes): result=[] maxidx=map(max,zip(*shapes)) for x in zip(*shapes): result.append([make_index(q,max(x)) for q in x]) return zip(*map(outer,zip(*result))) def gen_result_shape(shapes): maxl=max(map(len,shapes)) shapes=[x for x in shapes if len(x)==maxl] return tuple(map(max,zip(*shapes))) def broadcast_shapes(matrices,dims=2): shapes=[m.shape[:-dims] for m in matrices] maxlen=max([len(x) for x in shapes]) shapes=[tuple(1 for i in range(maxlen-len(x)))+tuple(x) for x in shapes] return shapes def make_flat(mat,dims=2): """make_flat(mat) takes list of arrays and changes shape of arrays to only have three indices the last two are assumed to represent a matrix convert x.shape=(n1,n2,...,nn,i,j) to x.shape=(n1*n2...*nn,i,j) return old list of old shapes shapes """ old=[] olddict={} newdict={} for m in mat: if id(m) not in olddict: #make sure shape is changed only once oldshape=m.shape olddict[id(m)]=oldshape totlen=reduce(operator.mul,m.shape[:-dims],1) shape=(totlen,)+m.shape[-dims:] newdict[id(m)]=shape else: shape=newdict[id(m)] oldshape=olddict[id(m)] old.append(oldshape) m.shape=shape return old def make_broad(mat,dims=2): """make_broad(mat,dims) change shape of arrays in the list mat so they all have the same number of dimensions. To simplify further operations that want to do broadcasting operations assuming the last <dims> are part of the object being used for operations. dims=2 would be used for operations on matrices. """ old=[] nshape=broadcast_shapes(mat,dims=dims) olddict={} newdict={} for m,s in zip(mat,nshape): if id(m) not in olddict: #make sure shape is changed only once oldshape=m.shape if s==tuple(): newshape=(1,)+m.shape[-dims:] else: newshape=s+m.shape[-dims:] newdict[id(m)]=newshape olddict[id(m)]=oldshape else: newshape=newdict[id(m)] oldshape=olddict[id(m)] m.shape=newshape old.append(oldshape) return old def unmake_flat(mat,old): for m,s in zip(mat,old): m.shape=s def broad_fun(infunc,dims=2,retvalue=None): """broad_fun(infunc,dims,retvalue) returns a new function based on infunc that does broadcasting over the n-dims dimensions of its input variables. retvalue is used to specify which returnvalue is used for infuncs that return more than one value. example: >>> from Numeric import matrixmultiply,array,NewAxis >>> mul=broad_fun(matrixmultiply,2) >>> x=array([[1.,2],[3,4]]) >>> y=repeat(repeat(x[NewAxis,NewAxis],3,0),2,1) >>> y=y*array([1,2,3])[:,NewAxis,NewAxis,NewAxis] >>> mul(x,y) array([[[[ 7., 10.], [ 15., 22.]], [[ 7., 10.], [ 15., 22.]]], [[[ 14., 20.], [ 30., 44.]], [[ 14., 20.], [ 30., 44.]]], [[[ 21., 30.], [ 45., 66.]], [[ 21., 30.], [ 45., 66.]]]]) """ if retvalue is not None: def func(*mat): return infunc(*mat)[retvalue] else: func=infunc def newfunc(*mat): oldshape=make_flat(mat) testresult=func(*[m[0] for m in mat]) unmake_flat(mat,oldshape) oldshape=make_broad(mat) shapes=[m.shape[:-dims] for m in mat] try: resultshape=gen_result_shape(shapes)+testresult.shape utshape=gen_result_shape([m[:-dims] for m in oldshape])+testresult.shape result=zeros(resultshape,testresult.typecode()) except AttributeError: resultshape=gen_result_shape(shapes) utshape=gen_result_shape([m[:-dims] for m in oldshape]) if utshape==tuple(): return testresult result=zeros(resultshape).astype(type(testresult)) indices=gen_iter_idx(shapes+[resultshape]) for i in indices: residx=i[-1] result[residx]=func(*[m[idx] for m,idx in zip(mat,i[:-1])]) unmake_flat(mat,oldshape) result.shape=utshape return result newfunc.__name__="broadcasting_%s"%infunc.__name__ return newfunc matmul=broad_fun(Numeric.matrixmultiply) det=broad_fun(LinearAlgebra.determinant) inv=broad_fun(LinearAlgebra.inverse) geninv=broad_fun(LinearAlgebra.generalized_inverse) eigenvalues=broad_fun(LinearAlgebra.eigenvalues) linear_least_squares=broad_fun(LinearAlgebra.linear_least_squares,retvalue=0) solve_linear_equations=broad_fun(LinearAlgebra.solve_linear_equations) def _test(): import doctest doctest.testmod() if __name__=="__main__": x=array([[1.,2],[3,4]]) z=array([ [[1.,2],[3,4]], [[4.,3],[2,1]]]) y=repeat(repeat(x[NewAxis,NewAxis],3,0),2,1) a=array([[1,2],[1.1,2.05],[0.98,1.97]]) b=array([[.7],[.75],[.69]]) A=repeat(repeat(a[NewAxis,NewAxis],3,0),2,1)