[PYTHON MATRIX-SIG] Array.py

Hinsen Konrad hinsenk@ere.umontreal.ca
Tue, 26 Sep 1995 21:39:11 -0400


# J-style array class

# Arrays are represented by a scalar or a list, possibly containing
# other lists in case of higher-rank arrays. Array rank is limited
# only by the user's patience.

# Send comments to Konrad Hinsen <hinsenk@ere.umontreal.ca>


import math, string


######################################################################
# Error type

ArrayError = 'ArrayError'

######################################################################
# Various functions that do the real work. Classes follow.

# Construct string representation of array

def _toString(data, dimension):
    s = ''
    if dimension == 0:
	s = s + `data`
    elif dimension == 1:
	for e in data:
	    s = s + `e` + ' '
    else:
	separator = (dimension-1)*'\n'
	for e in data:
	    s = s + _toString(e,dimension-1) + separator
    return string.strip(s)


# Find the shape of an array and check for consistency

def _shape(data):
    if type(data) == type([]):
	if data and type(data[0]) == type([]):
	    shapes = map(lambda x:_shape(x),data)
	    for i in range(1,len(shapes)):
		if shapes[i] != shapes[0]:
		    raise ArrayError, 'Inconsistent shapes'
	    shape = [len(data)]
	    shape = shape + shapes[0]
	    return shape
	else:
	    return [len(data)]
    else:
	return []


# Construct a one-dimensional list of all array elements

def __ravel(data):
    if type(data) == type([]):
	if type(data[0]) == type([]):
	    return reduce(lambda a,b: a+b,
			  map(lambda x: __ravel(x), data))
	else:
	    return data
    else:
	return [data]

def _ravel(array):
    return Array(__ravel(array._data), reduce(lambda a,b: a*b, array._shape, 1))


# Reshape an array

def _reshape(array, shape):
    array = ravel(array)
    if len(shape._data) == 0:
	return take(array,0)
    else:
	array = array._data
	shape = shape._data
	n = reduce(lambda a,b: a*b, shape)
	if n > len(array):
	    nr = (n+len(array)-1)/len(array)
	    array = (nr*array)[:n]
	elif n < len(array):
	    array = array[:n]
	for i in range(len(shape)-1, 0, -1):
	    d = shape[i]
	    n = n/d
	    for j in range(n):
		array[j:j+d] = [array[j:j+d]]
	return Array(array,shape)


# Map a function on the first dimensions of an array

def _extract(a, index, dimension):
    if len(a[1]) < dimension:
	return a
    else:
	return (a[0][index], a[1][1:], a[2])

def _map(function, arglist, max_frame, scalar_flag):
    result = []
    if len(max_frame) == 0:
	if scalar_flag:
	    result = apply(function,tuple(map(lambda a: a[0], arglist)))
	else:
	    result = apply(function,tuple(map(lambda a: Array(a[0],a[2]),
					      arglist)))._data
    else:
	for index in range(max_frame[0]):
	    result.append(_map(function,
			       map(lambda a,i=index,d=len(max_frame):
				   _extract(a,i,d), arglist),
			       max_frame[1:], scalar_flag))
    return result


# Reduce an array with a given binary function

def _reduce(function, array):
    function = function[0]
    array = array[0]
    if len(array._shape) == 0:
	return array
    elif array._shape[0] == 0:
	return reshape(function._neutral, array._shape[1:])
    else:
	result = Array(array._data[0], array._shape[1:])
	for i in range(1,array._shape[0]):
	    result = function(result, Array(array._data[i], array._shape[1:]))
	return result


# Find the higher of two ranks

def _maxrank(a,b):
    if a == None or b == None:
	return None
    else:
	return max(a,b)

######################################################################
# Array class definition

class Array:

    def __init__(self, scalar_or_list, shape = None):
	self._data = scalar_or_list
	if shape == None:
	    self._shape = _shape(self._data)
	else:
	    self._shape = shape

    def __str__(self):
	return _toString(self._data,len(self._shape))

    __repr__ = __str__

    def __len__(self):
	if type(self._data) == type([]):
	    return len(self._data)
	else:
	    return 1

    def __getitem__(self, index):
	return take(self, index)

    def __getslice__(self, i, j):
	return take(self, range(i,j))

    def __add__(self, other):
	return sum(self, other)
    __radd__ = __add__

    def __sub__(self, other):
	return difference(self, other)
    def __rsub__(self, other):
	return difference(other, self)

    def __mul__(self, other):
	return product(self, other)
    __rmul__ = __mul__

    def __div__(self, other):
	return quotient(self, other)
    def __rdiv__(self, other):
	return quotient(other, self)

    def __pow__(self,other):
	return power(self, other)
    def __rpow__(self,other):
	return power(other, self)

    def __neg__(self):
	return 0-self


# Check for arrayness

def isArray(x):
    return hasattr(x,'_shape')


######################################################################
# Array function class

class ArrayFunction:

    def __init__(self, function, ranks, intrinsic_ranks=None):
	self._function = function
	if isArray(ranks):
	    self._ranks = ranks._data
	elif type(ranks) == type([]):
	    self._ranks = ranks
	else:
	    self._ranks = [ranks]
	if intrinsic_ranks == None:
	    self._intrinsic_ranks = self._ranks
	else:
	    self._intrinsic_ranks = intrinsic_ranks
	if len(self._ranks)  == 1:
	    self._ranks = len(self._intrinsic_ranks)*self._ranks
	    

    def __call__(self, *args):
	if len(self._ranks) != len(args):
	    raise ArrayError, 'Wrong number of arguments for an array function'
	arglist = []
	framelist = []
	shapelist = []
	for i in range(len(args)):
	    if isArray(args[i]):
		arglist.append(args[i])
	    else:
		arglist.append(Array(args[i]))
	    shape = arglist[i]._shape
	    rank = self._ranks[i]
	    intrinsic_rank = self._intrinsic_ranks[i]
	    if rank == None:
		cell = 0
	    elif rank < 0:
		cell = min(-rank,len(shape))
	    else:
		cell = max(0,len(shape)-rank)
	    if intrinsic_rank != None:
		cell = max(cell,len(shape)-intrinsic_rank)
	    framelist.append(shape[:cell])
	    shapelist.append(shape[cell:])
	max_frame = []
	for frame in framelist:
	    if len(frame) > len(max_frame):
		max_frame = frame
	for i in range(len(framelist)):
	    if framelist[i] != max_frame[len(max_frame)-len(framelist[i]):]:
		raise ArrayError, 'Incompatible arguments'
	scalar_function = reduce(lambda a,b:_maxrank(a,b),
				 self._intrinsic_ranks) == 0
	return Array(_map(self._function, map(lambda a,b,c: (a._data,b,c),
					      arglist, framelist, shapelist),
			  max_frame, scalar_function))

    def __getitem__(self, ranks):
	return ArrayFunction(self._function,ranks,self._intrinsic_ranks)


class BinaryArrayFunction(ArrayFunction):

    def __init__(self, function, neutral_element, ranks, intrinsic_ranks=None):
	ArrayFunction.__init__(self, function, ranks, intrinsic_ranks)
	self._neutral = neutral_element
	self.over = ArrayFunction(ArrayOperator(_reduce, [self]), [None])

    def __getitem__(self, ranks):
	return BinaryArrayFunction(self._function, self._neutral,
				   ranks, self._intrinsic_ranks)


class ArrayOperator:

    def __init__(self, operator, function_list):
	self._operator = operator
	self._functions = function_list

    def __call__(self, *args):
	return apply(self._operator, (self._functions, args))


######################################################################
# Array functions

# Structural functions
shape =   ArrayFunction(lambda a: Array(a._shape), [None])
reshape = ArrayFunction(_reshape, [None, 1])
ravel =   ArrayFunction(_ravel, [None])
take =    ArrayFunction(lambda a,i: Array(a._data[i._data], a._shape[1:]),
			[None, 0])

# Elementwise binary functions
_sum =        ArrayFunction(lambda a,b: a+b, [0, 0])
_difference = ArrayFunction(lambda a,b: a-b, [0, 0])
_product =    ArrayFunction(lambda a,b: a*b, [0, 0])
_quotient =   ArrayFunction(lambda a,b: a/b, [0, 0])
_power =      ArrayFunction(pow, [0, 0])
sum  =        BinaryArrayFunction(_sum, 0, [None, None])
difference  = BinaryArrayFunction(_difference, 0, [None, None])
product  =    BinaryArrayFunction(_product, 1, [None, None])
quotient  =   BinaryArrayFunction(_quotient, 1, [None, None])
power =       BinaryArrayFunction(_power, 1, [None, None])

# Scalar functions of one variable
sqrt  = ArrayFunction(math.sqrt, [0])
exp   = ArrayFunction(math.exp, [0])
log   = ArrayFunction(math.log, [0])
log10 = ArrayFunction(math.log10, [0])
sin   = ArrayFunction(math.sin, [0])
cos   = ArrayFunction(math.cos, [0])
tan   = ArrayFunction(math.tan, [0])
asin  = ArrayFunction(math.asin, [0])
acos  = ArrayFunction(math.acos, [0])
atan  = ArrayFunction(math.atan, [0])
sinh  = ArrayFunction(math.sinh, [0])
cosh  = ArrayFunction(math.cosh, [0])
tanh  = ArrayFunction(math.tanh, [0])
floor = ArrayFunction(math.floor, [0])
ceil  = ArrayFunction(math.ceil, [0])

=================
MATRIX-SIG  - SIG on Matrix Math for Python

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