home *** CD-ROM | disk | FTP | other *** search
/ MacHack 2000 / MacHack 2000.toast / pc / The Hacks / MacHacksBug / Python 1.5.2c1 / Extensions / Numerical / Lib / MLab.py < prev    next >
Encoding:
Python Source  |  2000-06-23  |  7.3 KB  |  279 lines
  1. """Matlab(tm) compatibility functions.
  2.  
  3. This will hopefully become a complete set of the basic functions available in
  4. matlab.  The syntax is kept as close to the matlab syntax as possible.  One 
  5. fundamental change is that the first index in matlab varies the fastest (as in 
  6. FORTRAN).  That means that it will usually perform reductions over columns, 
  7. whereas with this object the most natural reductions are over rows.  It's perfectly
  8. possible to make this work the way it does in matlab if that's desired.
  9. """
  10. from Numeric import *
  11.  
  12. # Elementary Matrices
  13.  
  14. # zeros is from matrixmodule in C
  15. # ones is from Numeric.py
  16.  
  17. import RandomArray
  18. def rand(*args):
  19.     """rand(d1,...,dn) returns a matrix of the given dimensions
  20.     which is initialized to random numbers from a uniform distribution
  21.     in the range [0,1).
  22.     """
  23.     return RandomArray.random(args)
  24.  
  25. def eye(N, M=None, k=0, typecode=None):
  26.     """eye(N, M=N, k=0, typecode=None) returns a N-by-M matrix where the 
  27.     k-th diagonal is all ones, and everything else is zeros.
  28.     """
  29.     if M == None: M = N
  30.     if type(M) == type('d'): 
  31.         typecode = M
  32.         M = N
  33.     m = equal(subtract.outer(arange(N), arange(M)),-k)
  34.     return asarray(m,typecode=typecode)
  35.  
  36. def tri(N, M=None, k=0, typecode=None):
  37.     """tri(N, M=N, k=0, typecode=None) returns a N-by-M matrix where all
  38.     the diagonals starting from lower left corner up to the k-th are all ones.
  39.     """
  40.     if M == None: M = N
  41.     if type(M) == type('d'): 
  42.         typecode = M
  43.         M = N
  44.     m = greater_equal(subtract.outer(arange(N), arange(M)),-k)
  45.     return m.astype(typecode)
  46.     
  47. # Matrix manipulation
  48.  
  49. def diag(v, k=0):
  50.     """diag(v,k=0) returns the k-th diagonal if v is a matrix or
  51.     returns a matrix with v as the k-th diagonal if v is a vector.
  52.     """
  53.     v = asarray(v)
  54.     s = v.shape
  55.     if len(s)==1:
  56.         n = s[0]+abs(k)
  57.         if k > 0:
  58.             v = concatenate((zeros(k, v.typecode()),v))
  59.         elif k < 0:
  60.             v = concatenate((v,zeros(-k, v.typecode())))
  61.         return eye(n, k=k)*v
  62.     elif len(s)==2:
  63.         v = add.reduce(eye(s[0], s[1], k=k)*v)
  64.         if k > 0: return v[k:]
  65.         elif k < 0: return v[:k]
  66.         else: return v
  67.     else:
  68.             raise ValueError, "Input must be 1- or 2-D."
  69.     
  70.  
  71. def fliplr(m):
  72.     """fliplr(m) returns a 2-D matrix m with the rows preserved and
  73.     columns flipped in the left/right direction.  Only works with 2-D
  74.     arrays.
  75.     """
  76.     m = asarray(m)
  77.     if len(m.shape) != 2:
  78.         raise ValueError, "Input must be 2-D."
  79.     return m[:, ::-1]
  80.  
  81. def flipud(m):
  82.     """flipud(m) returns a 2-D matrix with the columns preserved and
  83.     rows flipped in the up/down direction.  Only works with 2-D arrays.
  84.     """
  85.     m = asarray(m)
  86.     if len(m.shape) != 2:
  87.         raise ValueError, "Input must be 2-D."
  88.     return m[::-1]
  89.     
  90. # reshape(x, m, n) is not used, instead use reshape(x, (m, n))
  91.  
  92. def rot90(m, k=1):
  93.     """rot90(m,k=1) returns the matrix found by rotating m by k*90 degrees
  94.     in the counterclockwise direction.
  95.     """
  96.     m = asarray(m)
  97.     if len(m.shape) != 2:
  98.         raise ValueError, "Input must be 2-D."
  99.     k = k % 4
  100.     if k == 0: return m
  101.     elif k == 1: return transpose(fliplr(m))
  102.     elif k == 2: return fliplr(flipud(m))
  103.     elif k == 3: return fliplr(transpose(m))
  104.  
  105. def tril(m, k=0):
  106.     """tril(m,k=0) returns the elements on and below the k-th diagonal of
  107.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  108.     diagonal.
  109.     """
  110.     return tri(m.shape[0], m.shape[1], k=k, typecode=m.typecode())*m
  111.  
  112. def triu(m, k=0):
  113.     """triu(m,k=0) returns the elements on and above the k-th diagonal of
  114.     m.  k=0 is the main diagonal, k > 0 is above and k < 0 is below the main
  115.     diagonal.
  116.     """    
  117.     return (1-tri(m.shape[0], m.shape[1], k-1, m.typecode()))*m 
  118.  
  119. # Data analysis
  120.  
  121. # Basic operations
  122. def max(m):
  123.     """max(m) returns the maximum along the first dimension of m.
  124.     """
  125.     return maximum.reduce(m)
  126.  
  127. def min(m):
  128.     """min(m) returns the minimum along the first dimension of m.
  129.     """
  130.     return minimum.reduce(m)
  131.  
  132. # Actually from BASIS, but it fits in so naturally here...
  133.  
  134. def ptp(m):
  135.     """ptp(m) returns the maximum - minimum along the first dimension of m.
  136.     """
  137.     return max(m)-min(m)
  138.  
  139. def mean(m):
  140.     """mean(m) returns the mean along the first dimension of m.  Note:  if m is
  141.     an integer array, integer division will occur.
  142.     """
  143.     return add.reduce(m)/len(m)
  144.  
  145. # sort is done in C but is done row-wise rather than column-wise
  146. def msort(m):
  147.     """msort(m) returns a sort along the first dimension of m as in MATLAB.
  148.     """
  149.     return transpose(sort(transpose(m)))
  150.  
  151. def median(m):
  152.     """median(m) returns a mean of m along the first dimension of m.
  153.     """
  154.     return msort(m)[m.shape[0]/2]
  155.  
  156. def std(m):
  157.     """std(m) returns the standard deviation along the first
  158.     dimension of m.  The result is unbiased meaning division by len(m)-1.
  159.     """
  160.     mu = mean(m)
  161.     return sqrt(add.reduce(pow(m-mu,2)))/sqrt(len(m)-1)
  162.  
  163. def sum(m):
  164.     """sum(m) returns the sum of the elements along the first
  165.     dimension of m.
  166.     """
  167.     return add.reduce(m)
  168.  
  169. def cumsum(m):
  170.     """cumsum(m) returns the cumulative sum of the elements along the
  171.     first dimension of m.
  172.     """
  173.     return add.accumulate(m)
  174.  
  175. def prod(m):
  176.     """prod(m) returns the product of the elements along the first
  177.     dimension of m.
  178.     """
  179.     return multiply.reduce(m)
  180.  
  181. def cumprod(m):
  182.     """cumprod(m) returns the cumulative product of the elments along the
  183.     first dimension of m.
  184.     """
  185.     return multiply.accumulate(m)
  186.  
  187. def trapz(y, x=None):
  188.     """trapz(y,x=None) integrates y = f(x) using the trapezoidal rule.
  189.     """
  190.     if x == None: d = 1
  191.     else: d = diff(x)
  192.     return sum(d * (y[1:]+y[0:-1])/2)
  193.  
  194. def diff(x, n=1):
  195.     """diff(x,n=1) calculates the first-order, discrete difference
  196.     approximation to the derivative.
  197.     """
  198.     if n > 1:
  199.         return diff(x[1:]-x[:-1], n-1)
  200.     else:
  201.         return x[1:]-x[:-1]
  202.     
  203. def corrcoef(x, y=None):
  204.     """The correlation coefficients
  205.     """
  206.     c = cov(x, y)
  207.     d = diag(c)
  208.     return c/sqrt(multiply.outer(d,d))
  209.  
  210. def cov(m,y=None):
  211.     m = asarray(m)
  212.     mu = mean(m)
  213.     if y != None: m = concatenate((m,y))
  214.     sum_cov = 0.0
  215.     for v in m:
  216.         sum_cov = sum_cov+multiply.outer(v,v)
  217.     return (sum_cov-len(m)*multiply.outer(mu,mu))/(len(m)-1.0)
  218.  
  219. # Added functions supplied by Travis Oliphant
  220. import LinearAlgebra
  221. def squeeze(a):
  222.     "squeeze(a) removes any ones from the shape of a"
  223.     b = asarray(a.shape)
  224.     reshape (a, tuple (compress (not_equal (b, 1), b)))
  225.     return
  226.  
  227. def kaiser(M,beta):
  228.     """kaiser(M, beta) returns a Kaiser window of length M with shape parameter
  229.     beta. It depends on the cephes module for the modified bessel function i0.
  230.     """
  231.     import cephes
  232.     n = arange(0,M)
  233.     alpha = (M-1)/2.0
  234.     return cephes.i0(beta * sqrt(1-((n-alpha)/alpha)**2))/cephes.i0(beta)
  235.  
  236. def blackman(M):
  237.     """blackman(M) returns the M-point Blackman window.
  238.     """
  239.     n = arange(0,M)
  240.     return 0.42-0.5*cos(2*pi*n/M) + 0.08*cos(4*pi*n/M)
  241.  
  242.  
  243. def bartlett(M):
  244.     """bartlett(M) returns the M-point Bartlett window.
  245.     """
  246.     n = arange(0,M)
  247.     return where(less_equal(n,M/2.0),2.0*n/M,2-2.0*n/M)
  248.  
  249. def hanning(M):
  250.     """hanning(M) returns the M-point Hanning window.
  251.     """
  252.     n = arange(0,M)
  253.     return 0.5-0.5*cos(2*pi*n/M)
  254.  
  255. def hamming(M):
  256.     """hamming(M) returns the M-point Hamming window.
  257.     """
  258.     n = arange(0,M)
  259.     return 0.54-0.46*cos(2*pi*n/M)
  260.  
  261. def sinc(x):
  262.     """sinc(x) returns sin(pi*x)/(pi*x) at all points of array x.
  263.     """
  264.     return where(equal(x,0.0),1,sin(pi*x)/(pi*x))
  265.  
  266.  
  267. def eig(v):
  268.     """[x,v] = eig(m) returns the the eigenvalues of m in x and the corresponding
  269.     eigenvectors in the rows of v.
  270.     """
  271.     return LinearAlgebra.eigenvectors(v)
  272.  
  273. def svd(v):
  274.     """[u,x,v] = svd(m) return the singular value decomposition of m.
  275.     """
  276.     return LinearAlgebra.singular_value_decomposition(v)
  277.  
  278.  
  279.