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
27 """
28 Null logging in order to avoid warning messages in client applications.
29 """
30 - def emit(self, record):
32
33 _h = NullHandler()
34 _logger = logging.getLogger('util')
35 _logger.addHandler(_h)
36
38 """ Maps every element of to its position."""
39 return zip(range(len(v)), v)
40
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
48
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):
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:
66 return False
67 for y2 in range(y+1, h):
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):
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):
78 m[y][x] /= c
79 return True
80
81
82
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
94 """
95 @return: the inverse of the matrix M
96 """
97
98
99 m2 = [row[:]+[int(i==j) for j in range(len(M) )] for i,row in enumerate(M) ]
100
101 return [row[len(M[0]):] for row in m2] if gauss_jordan(m2) else None
102
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