Package pod :: Module util
[frames] | no frames]

Source Code for Module pod.util

  1  """ 
  2  Operations on matrices and various tools. 
  3   
  4  @author: Christophe Alexandre <ch.alexandre at bluewin dot ch> 
  5   
  6  @license: 
  7  Copyright(C) 2010 Christophe Alexandre 
  8   
  9  This program is free software: you can redistribute it and/or modify 
 10  it under the terms of the GNU Lesser General Public License as published by 
 11  the Free Software Foundation, either version 3 of the License, or 
 12  (at your option) any later version. 
 13   
 14  This program is distributed in the hope that it will be useful, 
 15  but WITHOUT ANY WARRANTY; without even the implied warranty of 
 16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 17  GNU General Public License for more details. 
 18   
 19  You should have received a copy of the GNU Lesser General Public License 
 20  along with this program.  If not, see <http://www.gnu.org/licenses/lgpl.txt>. 
 21  """ 
 22   
 23  import math 
 24  import logging 
 25   
26 -class NullHandler(logging.Handler):
27 """ 28 Null logging in order to avoid warning messages in client applications. 29 """
30 - def emit(self, record):
31 pass
32 33 _h = NullHandler() 34 _logger = logging.getLogger('util') 35 _logger.addHandler(_h) 36
37 -def numbering(v):
38 """ Maps every element of to its position.""" 39 return zip(range(len(v)), v)
40
41 -def prod_scalar(v1, v2):
42 assert len(v1) == len(v2), 'input vectors must be of the same size' 43 prod = map(lambda x: x[0] * x[1], zip(v1, v2)) 44 return sum(prod)
45
46 -def norm(v):
47 return math.sqrt(prod_scalar(v, v))
48
49 -def gauss_jordan(m, eps=1E-10):
50 """ 51 Puts given matrix (2D array) into the Reduced Row Echelon Form. 52 53 NOTE: make sure all the matrix items support fractions! Int matrix will NOT work! 54 Written by Jarno Elonen in April 2005, released into Public Domain 55 56 @return: True if successful, False if 'm' is singular. 57 """ 58 (h, w) = (len(m), len(m[0])) 59 for y in range(0,h): 60 maxrow = y 61 for y2 in range(y+1, h): # Find max pivot 62 if abs(m[y2][y]) > abs(m[maxrow][y]): 63 maxrow = y2 64 (m[y], m[maxrow]) = (m[maxrow], m[y]) 65 if abs(m[y][y]) <= eps: # Singular? 66 return False 67 for y2 in range(y+1, h): # Eliminate column y 68 c = m[y2][y] / m[y][y] 69 for x in range(y, w): 70 m[y2][x] -= m[y][x] * c 71 for y in range(h-1, 0-1, -1): # Backsubstitute 72 c = m[y][y] 73 for y2 in range(0,y): 74 for x in range(w-1, y-1, -1): 75 m[y2][x] -= m[y][x] * m[y2][y] / c 76 m[y][y] /= c 77 for x in range(h, w): # Normalize row y 78 m[y][x] /= c 79 return True
80 81 # Auxiliary functions contribution by Eric Atienza 82
83 -def system_solve(M, b):
84 """ 85 solves M*x = b 86 @param M: a matrix in the form of a list of list 87 @param b: a vector in the form of a simple list of scalars 88 @return: vector x so that M*x = b 89 """ 90 m2 = [row[:]+[right] for row,right in zip(M,b) ] 91 return [row[-1] for row in m2] if gauss_jordan(m2) else None
92
93 -def mat_inverse(M):
94 """ 95 @return: the inverse of the matrix M 96 """ 97 #clone the matrix and append the identity matrix 98 # [int(i==j) for j in range_M] is nothing but the i(th row of the identity matrix 99 m2 = [row[:]+[int(i==j) for j in range(len(M) )] for i,row in enumerate(M) ] 100 # extract the appended matrix (kind of m2[m:,...] 101 return [row[len(M[0]):] for row in m2] if gauss_jordan(m2) else None
102
103 -def zeros(size, zero=0):
104 """ 105 @param size: a tuple containing dimensions of the matrix 106 @param zero: the value to use to fill the matrix (zero by default) 107 @return: matrix of dimension size 108 """ 109 return [zeros(s[1:] ) for i in range(s[0] ) ] if not len(s) else zero
110