[PYTHON MATRIX-SIG] MLab.py for 1.0alpha1

Janne Sinkkonen janne@avocado.pc.helsinki.fi
11 Aug 1996 18:20:28 +0300


Here is a (hopefully mostly) working MLab.py for 1.0beta1. I have
tried to test it with 1D and 2D arrays, but not with higher
dimensional arrays. max() and min() are changed to mmax() and mmin()
because of name collision with the builtin functions (are mmax() and
mmin() really needed?)

diag(<a 2D matrix>) could be more effective, but I left it as it is
because the diagonal() of Numeric.py does not seem to work right now.



"""Matlab(tm) compatibility functions.

This will hopefully become a complete set of the basic functions available in
matlab.  The syntax is kept as close to the matlab syntax as possible.  One 
fundamental change is that the first index in matlab varies the fastest (as in 
FORTRAN).  That means that it will usually perform reductions over columns, 
whereas with this object the most natural reductions are over rows.  It's perfectly
possible to make this work the way it does in matlab if that's desired.
"""
# ChangeLog:
# Sun Aug 11 18:12:20 EET DST 1996:
#      Hacked for 1.0alpha1 by <janne@iki.fi>

from Numeric import *

# Elementary Matrices

# zeros is from matrixmodule in C
# ones is from Numeric.py

## This should be replaced by Paul Dubois' URNG very soon now. 
import Ranf
def rand(*args):
	"""rand(d1,...,dn, typecode='d') returns a matrix of the given dimensions
	which is initialized to random number in the range [0,1).
	"""
	return Ranf.random_sample(args)


def eye(N, M=None, k=0, typecode=None):
	"""eye(N, M=N, k=0, typecode=None) returns a N-by-M matrix where the 
	k-th diagonal is all ones, and everything else is zeros.
	"""

	if M == None: M = N
	if type(M) == type('d'): 
		typecode = M
		M = N
	m = equal(subtract.outer(arange(N), arange(M)),-k)
	if (typecode == None):
	    return m
	else:
	    return m.asType(typecode)


def tri(N, M=None, k=0, typecode=None):
	if M == None: M = N
	if type(M) == type('d'): 
		typecode = M
		M = N
	m = greaterEqual(subtract.outer(arange(N), arange(M)),-k)
	if (typecode == None):
	    return m
	else:
	    return m.asType(typecode)

# Matrix manipulation

def diag(v, k=0):
	s = v.shape
	if len(s)==1:
		n = s[0]+abs(k)
		if k > 0:
		    v = concatenate((zeros((k,)),v))
		elif k < 0:
		    v = concatenate((v,zeros((-k,))))
		return multiply(eye(n, k=k), v)
	elif len(s)==2:
	    v = add.reduce(eye(s[0], s[1], k=k)*v)
	    i0=max(0,k)
	    i1=min(k+s[0],s[1])
	    if i0>i1:
		return array([])
	    else:
		return v[i0:i1]
		

		

def fliplr(m): 
    return m[:, ::-1]

def flipud(m):
    return m[::-1]

# reshape(x, m, n) is not used, instead use reshape(x, (m, n))

def rot90(m, k=1):
	k = k % 4
	if k == 0: return m
	elif k == 1: return fliplr(transpose(m))
	elif k == 2: return fliplr(flipud(m))
	elif k == 3: return transpose(fliplr(m))

def tril(m, k=0):
	return tri(m.shape[0], m.shape[1], k=k, typecode=m.typecode())*m

def triu(m, k=0):
	return (1-tri(m.shape[0], m.shape[1], k-1, m.typecode()))*m 

# Data analysis

# Basic operations
def mmax(m):
	return maximum.reduce(m)

def mmin(m):
	return minimum.reduce(m)

# Actually from BASIS, but it fits in so naturally here...

def ptp(m):
	return max(m)-min(m)

def mean(m):
	return add.reduce(m)/len(m)

# sort is done in C but is done row-wise rather than column-wise
def msort(m):
	return transpose(sort(transpose(m)))

def median(m):
	return msort(m)[m.shape[0]/2]

def std(m):
	mu = mean(m)
	return sqrt(add.reduce(pow(m-mu,2)))/sqrt(len(m)-1)

def sum(m):
	return add.reduce(m)

def cumsum(m):
	return add.accumulate(m)

def prod(m):
	return multiply.reduce(m)

def cumprod(m):
	return multiply.accumulate(m)

def trapz(y, x=None):
	"""Integrate f using the trapezoidal rule, where y is f(x).
	"""

	if x == None: d = 1
	else: d = diff(x)
	return sum(d * (y[1:]+y[0:-1])/2)

def diff(x, n=1):
	"""Discrete difference approximation to the derivative
	"""
	if n > 1:
	    return diff(x[1:]-x[:-1], n-1)
	else:
	    return x[1:]-x[:-1]
	
def corrcoef(x, y=None):
	"""The correlation coefficients
	"""
	if y==None: y=x
	c = cov(x, y)
	return c/sqrt(multiply.outer(diag(cov(x)),diag(cov(y))))

def cov(a,b=None):
    if b == None: b=a
    ma = mean(a)
    mb = mean(b)
    sum_cov = 0.0
    
    for ai,bi in map(lambda i,j: (i,j), a,b):
	sum_cov = sum_cov+multiply.outer(ai,bi)
    return (sum_cov-len(a)*multiply.outer(ma,mb))/(len(a)-1)




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

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