[Matrix-SIG] A Fixed version of MLab.py

[Travis E. Oliphant]

This is a multi-part message in MIME format.
--------------55B8511DC514F83C10E6D0E7
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Here is a fixed version of MLab.py with a couple of additions.

I've tested it a bit, but not extensively.




----------------------------------------------------
Travis Oliphant            200 First St SW          
    	                   Rochester MN 55905       
Ultrasound Research Lab	   (507) 286-5923           
Mayo Graduate School	   Oliphant.Travis@mayo.edu
--------------55B8511DC514F83C10E6D0E7
Content-Type: text/plain; charset=us-ascii; name="MLab.py"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline; filename="MLab.py"

"""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.
"""
from Numeric import *

# Elementary Matrices

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

import RandomArray
def rand(*args):
	"""rand(d1,...,dn) returns a matrix of the given dimensions
	which is initialized to random numbers from a uniform distribution
	in the range [0,1).
	"""
	return RandomArray.random(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)
	return asarray(m,typecode=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 = greater_equal(subtract.outer(arange(N), arange(M)),-k)
	return asarray(m,typecode=typecode)
	
# Matrix manipulation

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

def fliplr(m):
	m = asarray(m)
	if len(m.shape) != 2:
		raise ValueError, "Input must be 2-D."
	return m[:, ::-1]

def flipud(m):
	m = asarray(m)
	if len(m.shape) != 2:
		raise ValueError, "Input must be 2-D."
	return m[::-1]
	
# reshape(x, m, n) is not used, instead use reshape(x, (m, n))

def rot90(m, k=1):
	m = asarray(m)
	if len(m.shape) != 2:
		raise ValueError, "Input must be 2-D."
	k = k % 4
	if k == 0: return m
	elif k == 1: return transpose(fliplr(m))
	elif k == 2: return fliplr(flipud(m))
	elif k == 3: return fliplr(transpose(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 max(m):
	return maximum.reduce(m)

def min(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
	"""
	c = cov(x, y)
	d = diag(c)
	return c/sqrt(multiply.outer(d,d))

def cov(m,y=None):
	m = asarray(m)
	mu = mean(m)
	if y != None: m = concatenate((m,y))
	sum_cov = 0.0
	for v in m:
		sum_cov = sum_cov+multiply.outer(v,v)
	return (sum_cov-len(m)*multiply.outer(mu,mu))/(len(m)-1.0)


import LinearAlgebra
try:
	import cephes
except:
	print "Some functions depend on the cephes module which is not on your system."


def exit():
    import sys; sys.exit()

def squeeze(a):
    b = asarray(a.shape)
    b = compress(not_equal(b,1),b)
    a.shape = tuple(b)
    return

# This window depends on the cephes module for modified bessel function i0
def kaiser(M,beta):
    n = arange(0,M)
    alpha = (M-1)/2.0
    return cephes.i0(beta * sqrt(1-((n-alpha)/alpha)**2))/cephes.i0(beta)

def blackman(M):
    n = arange(0,M)
    return 0.42-0.5*cos(2*pi*n/M) + 0.08*cos(4*pi*n/M)

def rectangular(M):
    return ones((M),)

def bartlett(M):
    n = arange(0,M)
    return where(less_equal(n,M/2.0),2.0*n/M,2-2.0*n/M)

def hanning(M):
    n = arange(0,M)
    return 0.5-0.5*cos(2*pi*n/M)

def hamming(M):
    n = arange(0,M)
    return 0.54-0.46*cos(2*pi*n/M)

def sinc(x):
    return where(equal(x,0),1,sin(pi*x)/(pi*x))


def eig(v):
    return LinearAlgebra.eigenvectors(v)

def svd(v):
    return LinearAlgebra.singular_value_decomposition(v)










--------------55B8511DC514F83C10E6D0E7--