Matrix.py :  » Math » Matrix-package-for-Python » MatPy-0.4.0 » Python Open Source

Home
Python Open Source
1.3.1.2 Python
2.Ajax
3.Aspect Oriented
4.Blog
5.Build
6.Business Application
7.Chart Report
8.Content Management Systems
9.Cryptographic
10.Database
11.Development
12.Editor
13.Email
14.ERP
15.Game 2D 3D
16.GIS
17.GUI
18.IDE
19.Installer
20.IRC
21.Issue Tracker
22.Language Interface
23.Log
24.Math
25.Media Sound Audio
26.Mobile
27.Network
28.Parser
29.PDF
30.Project Management
31.RSS
32.Search
33.Security
34.Template Engines
35.Test
36.UML
37.USB Serial
38.Web Frameworks
39.Web Server
40.Web Services
41.Web Unit
42.Wiki
43.Windows
44.XML
Python Open Source » Math » Matrix package for Python 
Matrix package for Python » MatPy 0.4.0 » Matrix.py
#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <hzhu@users.sourceforge.net>.   Licence: GPL
# $Id: Matrix.py,v 1.17 2001/08/26 12:40:43 hzhu Exp $
"""
Matrix class with matlab-like interface.  Related functions.
Row vec = 1xn mat.  Col vec = nx1 mat.
1x1 mat is different from a number
Matrix multiplication follows linear algebra rules
"""

from MatPy import MatrixType
from MatPy import NaN,nan,Inf,inf

import Numeric, MLab, RandomArray, LinearAlgebra, string, __builtin__
from Numeric import Int

_eps = 2.2204e-16
_eps_approx = 1e-13
_mul = MLab.matrixmultiply
_array = MLab.array
_Axis = MLab.NewAxis

_abs = MLab.absolute
_sum = MLab.sum
_sqrt = MLab.sqrt
_pow = MLab.power
_min = MLab.min
_max = MLab.max
_minimum = Numeric.minimum
_maximum = Numeric.maximum

_trsp = MLab.transpose
_conj = MLab.conjugate
_diag = MLab.diag

_sol = LinearAlgebra.solve_linear_equations
_inv = LinearAlgebra.inverse
_eig = LinearAlgebra.eigenvalues

_ones = MLab.ones
_zeros = MLab.zeros
_eye = MLab.eye
_rand = RandomArray.random

_gt = Numeric.greater
_ge = Numeric.greater_equal
_lt = Numeric.less
_le = Numeric.less_equal
_eq = Numeric.equal

#------------------------------------------------------------------
class Matrix(MatrixType):
  "class Matrix : main matrix class"
  
  def __init__(self, data):
    "__init__: Initialize to a new Matrix of given data"
    if isinstance(data, Matrix):
      self.data = data.data
    else:
      self.data = _array(data)
      m = self.data.shape
      if len(m) == 2:    pass
      elif len(m) == 1:  self.data.shape = (1, m[0])
      elif len(m) == 0:  self.data.shape = (1, 1)
      else :        raise ValueError, "Matrix of shape %s"%`m`
    self.shape = self.data.shape
    self.typecode = self.data.typecode()

  def copy(self):
    "A new matrix with a copy of self.data"
    import copy
    data = copy.copy(self.data)
    return Matrix(data)

  #------------------------------------------------------------------
  # Methods that look like attributes
  
  def __getattr__(self, name):
    "Short cut to methods T H C I eig real imag"
    if   name == 'T':    return Matrix(_trsp(self.data))
    elif name == 'H':    return Matrix(_conj(_trsp(self.data)))
    elif name == 'C':    return Matrix(_conj(self.data))
    elif name == 'I':    return Matrix(_inv(self.data))
    elif name == 'eig':    return Matrix(_eig(self.data))
    elif name == 'real':
      if self.typecode == "D":  return Matrix(self.data.real)
      else:        return self
    elif name == 'imag':
      if self.typecode == "D":  return Matrix(self.data.imag)
      else:        return zeros(self.shape)
    else:
      try:
        return self.__dict__[name]
      except KeyError:
        raise AttributeError, name

  def __coerce__(self, other): return self, other

  def checkreal(self):
    if self.typecode=='D':
      print "norms of this complex value:",
      print "%s+j%s"%(norm(self.real),norm(self.imag))
    else:
      return 1

  #------------------------------------------------------------------
  # Indices and sclices
  
  def __getitem__(self, i):
    #print "getitem", i
    if type(i) is type(1):
      if self.shape[0] == 1:    return self.data[0,i]
      elif self.shape[1] == 1:  return self.data[i,0]
      else:            return self.data[i]
    elif type(i) is type((1,)):
      result = self.data[i]
      if type(i[0]) == type(1) and type(i[1]) == type(slice(1)):
        result = result[_Axis,:]
      elif type(i[1]) == type(1) and type(i[0]) == type(slice(1)):
        result = result[:,_Axis]
    else:
      raise TypeError, i, type(i)
    return to_Matrix(result)
  
  def __setitem__(self, i, item):
    #print "setitem"
    item = to_array(item)
    if type(i) is type(1):
      if self.shape[0] == 1:    self.data[0,i] = item
      elif self.shape[1] == 1:  self.data[i,0] = item
    else:              self.data[i] = item

  def __getslice__(self, i, j):
    #print "getslice", i, j, self.shape
    if self.shape[0] == 1:    result = self.data[0,i:j]
    elif self.shape[1] == 1:  result = self.data[i:j,0][:, _Axis]
    else:            result = self.data[i:j]
    return to_Matrix(result)

  def __setslice__(self, i, j, other):
    #print "setslice"
    if isinstance(other, Matrix):  tmp = other.data
    else:              tmp = other
    if self.shape[0] == 1:    self.data[0,i:j] = tmp
    elif self.shape[1] == 1:  self.data[i:j,0] = tmp
    else:            self.data[i:j] = tmp

  def __len__(self):        return len(self.data)

  #------------------------------------------------------------------
  # Elementwise binary operations
  
  def e_add (self, other):  return Matrix(self.data + to_array(other))
  def re_add(self, other):  return Matrix(to_array(other) + self.data)
  def e_sub (self, other):  return Matrix(self.data - to_array(other))
  def re_sub(self, other):  return Matrix(to_array(other) - self.data)

  def e_mul (self, other):  return Matrix(self.data * to_array(other))
  def re_mul(self, other):  return Matrix(to_array(other) * self.data)
  def e_div (self, other):  return Matrix(self.data / to_array(other))
  def re_div(self, other):  return Matrix(to_array(other) / self.data)

  def e_pow (self, n):    return Matrix(self.data ** n)
  def re_pow(self, n):    return Matrix(n ** self.data)

  # Operator names
  __tadd__  = e_add
  __rtadd__ = re_add
  __tsub__  = e_sub
  __rtsub__ = re_sub
  __tmul__  = e_mul
  __rtmul__ = re_mul
  __tdiv__  = e_div
  __rtdiv__ = re_div
  __tpow__  = e_pow
  __rtpow__ = re_pow
  
  def __neg__(self):    return Matrix(-self.data)

  def __int__(self):    return Matrix(self.data.astype(Int))

  #------------------------------------------------------------------
  # Elementwise comparison

  def e_cmp(self, other):
    "e_cmp(self, other): elementwise comparison (other may be scalar)"
    if not isinstance(other, Matrix):
      other = ones(size(self))*other
    return mmap2(cmp, self, other)

  def gt(self, other): return Matrix(_gt(self.data, to_array(other)))
  def ge(self, other): return Matrix(_ge(self.data, to_array(other)))
  def lt(self, other): return Matrix(_lt(self.data, to_array(other)))
  def le(self, other): return Matrix(_le(self.data, to_array(other)))
  def eq(self, other): return Matrix(_eq(self.data, to_array(other)))

  #------------------------------------------------------------------
  # Matrixwise binary operations
  def checkshape(self, other, op):
    "Check dimension is both are Matrix"
    if isinstance(other, Matrix) and self.shape != other.shape:
      raise ValueError, "%s shapes %s and %s" %(
        op, self.shape, other.shape)
    
  def m_add(self, other):
    self.checkshape(other, "add")
    return Matrix(self.data + to_array(other))
  def rm_add(self, other):
    #self.checkshape(other, "radd")
    return Matrix(to_array(other) + self.data)

  def m_sub(self, other):
    self.checkshape(other, "sub")
    return Matrix(self.data - to_array(other))
  def rm_sub(self, other):
    #self.checkshape(other, "rsub")
    return Matrix(to_array(other) - self.data)

  def m_mul(self, other):
    #if isinstance(other, MatrixType):    other = other.data
    return to_Matrix(_mul(self.data,  to_array(other)))
  def rm_mul(self, other):
    "Somehow this gets transposed by MLab"
    #if isinstance(other, MatrixType):    other = other.data
    return to_Matrix(_mul(to_array(other), self.data)).T

  def m_div(self, other):
    if isinstance(other, MatrixType):
      return solve(other.T, self.T).T
    else:
      return to_Matrix(self.data / other)
  def rm_div(self, other):
    if isinstance(other, MatrixType):
      return solve(self.T, other.T).T
    else:
      return other * self.I

  def m_sol(self, other):
    return solve(self, other)
  def rm_sol(self, other):
    return solve(other, self)

  def m_pow(self, n):
    return mfuncC(lambda x,n=n:_pow(x,n), self)
  def rm_pow(self, n):
    return mfuncC(lambda x,n=n:_pow(n,x), self)

  #------------------------------------------------------------------
  # Matrixwise augmented assignments
  
  def m_iadd(self, other):
    self.checkshape(other, "iadd")
    self.data = to_array(self.data + other)
    self.shape = self.data.shape
    self.typecode = self.data.typecode()
    return self

  def m_isub(self, other):
    self.checkshape(other, "isub")
    self.data = to_array(self.data - other)
    self.shape = self.data.shape
    self.typecode = self.data.typecode()
    return self

  def m_imul(self, other):
    if isinstance(other, MatrixType):    other = other.data
    self.data = _mul(self.data, other)
    self.shape = self.data.shape
    self.typecode = self.data.typecode()
    return self

  def m_idiv(self, other):
    "Temporary hack"
    if isinstance(other, MatrixType):
      self.data = _trsp(_sol(_trsp(other.data), _trsp(self.data)))
    else:
      self.data = self.data / other
    self.shape = self.data.shape
    self.typecode = self.data.typecode()
    return self

  def m_ipow(self, n):
    self.data = _pow(self.data, n)
    self.shape = self.data.shape
    self.typecode = self.data.typecode()
    return self

  # Operator names
  __add__   = m_add
  __radd__ = rm_add
  __sub__   = m_sub
  __rsub__ = rm_sub
  __mul__   = m_mul
  __rmul__ = rm_mul
  __div__   = m_div
  __rdiv__ = rm_div
  __sol__   = m_sol
  __rsol__ = rm_sol
  __pow__   = m_pow
  __rpow__ = rm_pow

  __iadd__ = m_iadd
  __isub__ = m_isub
  __imul__ = m_imul
  __idiv__ = m_idiv
  __ipow__ = m_ipow

  def __cmp__(self, other):    return norm(self-other) > 0


#------------------------------------------------------------------
# Casting

def to_list(data):
  "to_list(data) : convert to double list"
  data = to_array(data)
  return map(list, data)

def to_Matrix(data):
  "to_Matrix(data) : covert to Matrix if possible"
  if isinstance(data, Matrix):      return Matrix(data.data)
  elif type(data) == type(_array([])):  return Matrix(data)
  elif type(data) == type([]):      return Matrix(data)
  else:                  return data

def to_array(data):
  "to_array(data) : convert to array"
  if isinstance(data, Matrix):  return data.data
  else:              return _array(data)

def to_arrayC(data):
  "to_arrayC(data) : convert to complex array"
  if isinstance(data, Matrix):  return data.data+0j
  else:              return _array(data)+0j

def to_number(data):
  "to_number(data) : conbert to number if possible"
  if hasattr(data, 'shape'):
    if data.shape == (1,1): return data[0,0]
    elif data.shape == (1): return data[0]
    else:          raise ValueError, "Can't convert to number"
  else:            return data

def Matrix_c(x):
  "Matrix_c(x) : return a matrix from a list as a column"
  if type(x) not in [type([]), type(())]:
    raise TypeError, "argument must be a sequence"
  lists = map(list, map(to_array,x))
  import operator
  return Matrix(reduce (operator.add,  lists))

def Matrix_r(x):
  "Matrix_r(x) : return a matrix from a list as a row"
  if type(x) not in [type([]), type(())]:
    raise TypeError, "argument must be a sequence"
  return  Matrix_c(map(lambda a:a.T, x)).T

def Matrix_cr(x):
  "Matrix_cr(x) : return a matrix from a double list as a column of rows"
  return Matrix_c(map(Matrix_r, x))

def rows(x):
  "rows(x) : split x into a list of rows"
  #print x
  return map(to_Matrix, tuple(x.data))

def cols(x):
  "cols(x) : split x into a list of columns"
  return map(lambda x:x.T, rows(x.T))

def size(data):
  "size(data) : return data.shape"
  return to_Matrix(data).shape

#------------------------------------------------------------------
# General elementwise operations
# These are used when there is no NumPy functions.

def map1(f=None, *args):
  "map1(f, *args): a costly reimplementation of map"
  return apply(map, [f] + list(args))

def map2(f=None, *args):
  "map2(f, *args): elementwise f on matrices args"
  matlist = map(to_list, args)
  return apply(map, [map] + [[f]*len(args[0])] + matlist)

def mmap2(*args):
  "mmap2(f, *args): return Matrix as elementwise f on matrices args"
  return Matrix(apply(map2, args))

def reduce1(f, vec):
  "reduce1(f, list): a costly reimplementation of reduce"
  return reduce(f, vec)

def reduce2(f, mat):
  "reduce2(f, *args): elementwise f on matrices args"
  return reduce(f, map(lambda x,f=f:reduce(f,x), to_list(mat)))

def mreduce2(f, mat):
  "mreduce2(f, *args): return Matrix as elementwise f on matrices args"
  return reduce2(f, map(list, args))

def all(mat):
  "all(x) : return row 0-1 vector, indicates each column being all true"
  #return reduce2(operator.__and__, mat)
  return efunc(MLab.alltrue, mat)

def any(mat):
  "any(x) : return row 0-1 vector, indicates each column being partly true" 
  #return reduce2(operator.__or__, mat)
  return efunc(MLab.sometrue, mat)

def find(mat):
  "find(x) : return list of indices (i,j) for which x[i,j] is true"
  m,n = mat.shape
  result = []
  for i in xrange(m):
    for j in xrange(n):
      if mat[i,j]: result.append((i,j))
  return result

#------------------------------------------------------------------
# Linear algebra
# These used in efuncs.py and mfuncs.py to general functions.

def plainsolve(a, b):
  "solve(a,b) : return solution of a*x==b"
  aT = a.T
  return to_Matrix(_sol(a.data, b.data))

def solve(a, b):
  "solve(a,b) : return LMS solution of a*x==b"
  aT = a.T
  return to_Matrix(_sol((aT*a).data, (aT*b).data))

sol = solve

def eig(x):
  """eig(x) : return eigenvalues as an array, eigenvectors as a Matrix
  Definition: x * u.T = u.T * v
  """
  s = x.shape
  assert s[0] == s[1]
  v, u =  MLab.eig(x.data)
  return v, to_Matrix(u)

def svd(a):
  "svd(x) : return Matrix u, array s, Matrix v as svd decomposition"
  u, s, v =  MLab.svd(a.data)
  return Matrix(u), s, Matrix(v)

def efunc(f, *args):
  "efunc(f, *args) : elementwise function f of args"
  return to_Matrix(apply(f, map(to_array, args)))

def efuncC(f, *args):
  """efuncC(f, *args) : elementwise function f of args with possible complex elements"""
  y = to_Matrix(apply(f, map(to_arrayC, args)))
  return approx_real(y)

def efuncRC(*args):
  """efuncRC(*args) : elementwise function of matrices args,
  allow complex if one arg is complex"""
  for x in args[1:]:
    if x.typecode == "D":
      return apply(efuncC, args)
  else:
    return apply(efunc, args)

def mfunc(f, x):
  """mfunc(f, x) : matrix function f of matrix x.
     Note: Numeric defines (v,u) = eig(x) => x*u.T = u.T * Diag(v)
  """
  (v, u) = eig(x)
  uT = u.T
  V = Matrix(_diag(f(v)))
  y = uT * V * uT.I
  return approx_real(y)

def mfuncC(f, x):
  """mfuncC(f, x) : matrix function with possibly complex eigenvalues.
     Note: Numeric defines (v,u) = eig(x) => x*u.T = u.T * Diag(v)
  """
  (v, u) = eig(x)
  uT = u.T
  V = Matrix(_diag(f(v+0j)))
  y = uT * V * uT.I
  return approx_real(y)

def mfunc_p(f, x):
  "mfunc_p(f, x) : matrix function with positive eigenvalues"
  (v, u) = eig(x)
  if v.typecode() == 'D':
    if norm(Matrix(v.imag)) > norm(Matrix(v.real))*1e-14:
      raise TypeError, "%s only applicable to real numbers" % `f`
    else:
      v = v.real
  print v, v.typecode()
  
  if _min(v) + _eps < 0:
    raise TypeError,  "%s only applicable to positive numbers" % `f`
  else:
    uT = u.T
    V = diag(Matrix(f(__builtin__.max(v,0))))
    y = u.T * V * uT.I
    return approx_real(y)

def approx_real(x):
  "approx_real(x) : return x.real if |x.imag| < |x.real| * _eps_approx"
  try:
    xtype = x.typecode
  except AttributeError:
    if abs(x.imag) <= abs(x.real) * _eps_approx:
      return x.real
    else:
      return x
    
  if xtype == 'D' and norm(x.imag) <= norm(x.real) * _eps_approx:
    return x.real
  else:
    return x

#------------------------------------------------------------------
# New matrices

def Mrange(*i):
  "Mrange(*i) : row vector of range(i) - deprecated"
  return to_Matrix(apply(range,i))

def r_range(*i):
  "r_range(*i) : row vector of range(i)"
  return to_Matrix(apply(range,i))

def c_range(*i):
  "c_range(*i) : col vector of range(i)"
  return to_Matrix(apply(range,i)).T

def zeros(x, typecode='d'):
  "zeros(size, *typecode) : zero matrix of given size and type"
  return to_Matrix(_zeros(x, typecode))

def ones(x):
  "ones(size) : matrix of ones of given size"
  return to_Matrix(_ones(x))

def eye(x):
  "eye(size) : identity matrix of given size"
  return to_Matrix(_eye(x))

def rand(*s):
  "rand(size) : standard uniformly distributed random matrix of given size"
  return to_Matrix(apply(_rand, s))

def unit(n, m):
  "unit(n, m) : unit vector of dimension n, pointing at mth direction"
  x = zeros(n).T
  x[m] = 1.
  return x

# Matrices of the same shape

def sort(x):
  "sort(x) : return columnwise sorted"
  return efunc(MLab.msort, x)

def diff(x):
  "diff(x) : return columnwise diff"
  return efunc(MLab.diff, x.T).T

def diff0(x):
  "diff(x) : return columnwise diff, keeping first row.  Inverse of cumsum"
  return Matrix_c((a[0,:], diff(a)))

def trapz(x):
  "trapz(x) : return columnwise integration using trapzoidal rule"
  return efunc(MLab.trapz, x)
  
def cumsum(x):
  "cumsum(x) : return columnwise cumulative sum"
  return efunc(MLab.cumsum, x)

def cumprod(x):
  "cumprod(x) : return columnwise cumulative prod"
  return efunc(MLab.cumprod, x)

def flipud(x):
  "flipud(x) : return up-down flipped"
  return efunc(MLab.flipud, x)

def fliplr(x):
  "flipud(x) : return left-right flipped"
  return efunc(MLab.fliplr, x)

# Reduce to row vector

def norm1(x):
  "norm1(x) : return columnwise norms"
  return _sqrt(_sum(to_array(x)**2))

def sum(x):
  "sum(x) : return columnwise sum"
  return efunc(_sum, x)

def prod(x):
  "prod(x) : return columnwise product"
  return efunc(MLab.prod, x)

def mean(x):
  "mean(x) : return columnwise mean"
  return efunc(MLab.mean, x)

def median(x):
  "median(x) : return columnwise median"
  return efunc(MLab.median, x)

def std(x):
  "std(x) : return columnwise standard deviation"
  return efunc(MLab.std, x)

def ptp(x):
  "ptp(x) : return columnwise ptp?"
  return efunc(MLab.ptp, x)

def min(*args):
  """min(x) : return columnwise minimum if one argument,
  return elementwise min across multiple arguments"""
  if len(args) == 1:
    return efunc(_min, args[0])
  else:
    return to_Matrix(reduce(_minimum, map(to_array,args)))

def max(*args):
  """max(x) : return columnwise maximum if one argument,
  return elementwise max across multiple arguements"""
  if len(args) == 1:
    return efunc(_max, args[0])
  else:
    return to_Matrix(reduce(_maximum, map(to_array,args)))


# Different objects

def sum2(x):
  "sum2(x) : return sum of all elements"
  return _sum(_sum(to_array(x)))

def min2(x):
  "min2(x) : return min of all elements"
  return _min(_min(to_array(x)))

def max2(x):
  "max2(x) : return max of all elements"
  return _max(_max(to_array(x)))

def cov(x):
  "std(x) : return covariance matrix for columns of x"
  return efunc(MLab.cov, x)

def norm(x):
  "norm(x) : return Frobenius norm of x"
  return _sqrt(sum2(_abs(to_array(x))**2))

def mnorm(x):
  "norm(x) : return Hilbert-Schmidt (?) norm of x"
  return _sqrt(sum2(_abs(to_array(x.H * x))))

def Diag(x):
  "Diag(x) : return diag matrix if x is row or col vector"
  if x.shape[1] == 1:    x = x.T
  a = to_array(x)
  if a.shape[0] == 1:
    return  Matrix(_diag(a[0]))
  else:
    raise ValueError, "x must be a vector"

def diag(x):
  "diag(x) : return diagonal of x as row vector"
  if (x.shape[1] == 1 or x.shape[0] == 1) and x.shape[1]*x.shape[0] != 1:
    raise ValueError, "x cannot be a vector" + `x.shape`
  a = to_array(x)
  return  Matrix(_diag(a))
  
www.java2java.com | Contact Us
Copyright 2009 - 12 Demo Source and Support. All rights reserved.
All other trademarks are property of their respective owners.