1 """
2 Basic linear algebra components and functions.
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 os
25 import logging
26
27 from util import NullHandler
28 from util import numbering
29 from util import prod_scalar
30
31 _h = NullHandler()
32 _logger = logging.getLogger('linalg')
33 _logger.addHandler(_h)
34
36
37 - def __init__(self, dim_row, dim_col=None):
38 if dim_col is None:
39 dim_col = dim_row
40 self.dim_row = dim_row
41 self.dim_col = dim_col
42 self._vectors = dict()
43
50
52 assert i < self.get_dim_row(), 'row %d exceeding dimension %d' % (i, self.get_dim_row())
53 assert j < self.get_dim_col(), 'column %d exceeding dimension %d' % (j, self.get_dim_col())
54 if self._vectors.has_key(i):
55 v = self._vectors[i]
56 else:
57 v = self._create_vector()
58 return v.get_component(j)
59
66
68 table = list()
69 for i in range(self.get_dim_row()):
70 row = list()
71 for j in range(self.get_dim_col()):
72 row.append(self.get_value(i, j))
73 table.append(row)
74 return table
75
77 return (self.dim_row, self.dim_col)
78
81
84
87
89 msg = str(self.get_dimension())
90 assert i < self.get_dim_row(), 'row %d exceeding dimension %s' % (i, msg)
91 assert j < self.get_dim_col(), 'column %d exceeding dimension %s' % (j, msg)
92 if self._vectors.has_key(i):
93 v = self._vectors[i]
94 else:
95 v = self._create_vector()
96 self._vectors[i] = v
97 try:
98 v.set_component(j, float(val))
99 except TypeError, e:
100 logging.error('not a number: %s' % str(val))
101 raise e
102
104 """
105 Sub-matrix excluding the specified row and column.
106 """
107 m = Matrix(self.get_dim_row() - 1, self.get_dim_col() - 1)
108 dr = self.get_dim_row()
109 dc = self.get_dim_col()
110 for k, i in numbering(range(0, r) + range(r+1, dr)):
111 for l, j in numbering(range(0, c) + range(c+1, dc)):
112 v = self.get_value(i, j)
113 m.set_value(k, l, v)
114 return m
115
117 size = self.get_dim_row()
118 det = 0.0
119
120 if size == 1:
121 det = self.get_value(0, 0)
122 else:
123 for i in range(size):
124 minor = self.minor(0, i)
125 if i % 2 == 0:
126 sign = 1.0
127 else:
128 sign = -1.0
129 det += sign * self.get_value(0, i) * minor.determinant()
130 return det
131
143
146
156
165
166 - def scale(self, factor):
173
184
186 """
187 Full row rank pseudo inverse.
188
189 A+ = transp(A).inv(A.transp(A))
190 """
191 a = self
192 t_a = a.transpose()
193 a_t_a = a.multiply(t_a)
194 inv_a_transp_a = a_t_a.inverse()
195 return t_a.multiply(inv_a_transp_a)
196
203
205 out = '(M%dx%d)' % (self.get_dim_row(), self.get_dim_col()) + os.linesep
206 for row in range(self.get_dim_row()):
207 line = []
208 for col in range(self.get_dim_col()):
209 line.append(self.get_value(row, col))
210 out += ', '.join(map(str, line)) + os.linesep
211 return out
212
214 """
215 Makes a vector compatible with Matrix operations.
216 """
217
227
229 """
230 @todo: should somehow be merged with VectorMatrix.
231 """
232
234 self._dim = dim
235 self._values = dict()
236
239
241 """ Using specified values to initialize the vector. """
242 for n, v in numbering(values):
243 self.set_component(n, v)
244 return self
245
247 """ Using raw data (python list) to initialize the vector. """
248 self.set_values(*data)
249
251 assert i < self.get_length(), 'index %d exceeding dimension %d' % (i, self.get_length())
252 assert i >= 0, 'non positive index %d' % i
253 if not self._values.has_key(i):
254 return 0.0
255 else:
256 return self._values[i]
257
259 assert i < self.get_length(), 'index %d exceeding dimension %d' % (i, self.get_length())
260 assert i >= 0, 'non positive index %d' % i
261 self._values[i] = v
262
264 """ Raw data as built-in python list."""
265 l = list()
266 for i in range(self.get_length()):
267 l.append(self.get_component(i))
268 return l
269
270 - def sub(self, vector):
275
276 - def add(self, vector):
281
287
290
292 return math.sqrt(self.product(self))
293
297
298 - def units(self, unit_vector):
299 """
300 How many times the current vector fits in the specified units (in norm).
301
302 The sign has a meaning only if both vectors are colinear.
303 """
304 ratio = self.norm() / unit_vector.norm()
305 if self.sub(unit_vector).norm() > unit_vector.norm():
306
307 ratio = -ratio
308 return ratio
309
311 return self._values == other._values
312
314 return self._values != other._values
315
317 line = []
318 for col in range(self.get_length()):
319 line.append(self.get_component(col))
320 out = '(V)[' + ', '.join(map(str, line)) + ']'
321 return out
322
325
327 """
328 Equivalent to a Vector
329 """
333
335 line = []
336 for col in range(self.get_length()):
337 line.append(self.get_component(col))
338 out = '(P)[' + ', '.join(map(str, line)) + ']'
339 return out
340
342 matrix = Matrix(dimension)
343 for i in range(dimension):
344 matrix.set_value(i, i, 1.0)
345 return matrix
346
347 -def zero(dimension):
350
362