parsimony.utils.linalgs

1 # -*- coding: utf-8 -*- 2 """ 3 Created on Wed Jul 23 19:15:22 2014 4 5 Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved. 6 7 @author: Tommy Löfstedt 8 @email: tommy.loefstedt@cea.fr 9 @license: BSD 3-clause. 10 """ 11 import abc 12 import time 13 14 import numpy as np 15 import scipy.sparse as sparse 16 17 import consts 18 19 __all__ = ["MultipartArray"] 20 21

22 -class MultipartArray(object):

23

24 - def __init__(self, parts, vertical=True):

25 26 self.parts = list(parts) 27 self.vertical = vertical 28 29 self.common_dim = 1 if vertical else 0 30 self.shape = [0, 0] 31 self.shape[self.common_dim] = self.parts[0].shape[self.common_dim] 32 self.multipart_shape = [] 33 34 for i in xrange(len(self.parts)): 35 if len(self.parts[i].shape) != 2: 36 raise ValueError("MultipartArray is only defined for 2D " \ 37 "arrays") 38 if not self.parts[0].shape[self.common_dim] \ 39 == self.parts[i].shape[self.common_dim]: 40 raise ValueError("Input vectors must have a common dimension") 41 self.multipart_shape.append( 42 self.parts[i].shape[1 - self.common_dim]) 43 self.shape[1 - self.common_dim] += self.multipart_shape[-1] 44 45 self.shape = tuple(self.shape) 46 self.multipart_shape = tuple(self.multipart_shape)

47

48 - def get_parts(self):

49 return self.parts

50

51 - def toarray(self):

52 if self.vertical: 53 return np.vstack(self.parts) 54 else: 55 return np.hstack(self.parts)

56

57 - class __ops(object):

58 add = 0 59 sub = 1 60 mul = 2 61 div = 3

62

63 - def __iop(self, other, op):

64 if np.isscalar(other): 65 for i in xrange(len(self.parts)): 66 if op == self.__ops.add: 67 self.parts[i] += other 68 elif op == self.__ops.sub: 69 self.parts[i] -= other 70 elif op == self.__ops.mul: 71 self.parts[i] *= other 72 elif op == self.__ops.div: 73 self.parts[i] *= 1.0 / other 74 else: 75 raise ValueError("Operator not yet implemented!") 76 elif isinstance(other, MultipartArray): 77 other_parts = other.get_parts() 78 for i in xrange(len(self.parts)): 79 if op == self.__ops.add: 80 self.parts[i] += other_parts[i] 81 elif op == self.__ops.sub: 82 self.parts[i] -= other_parts[i] 83 elif op == self.__ops.mul: 84 self.parts[i] *= other_parts[i] 85 else: 86 raise ValueError("Operator not yet implemented!") 87 elif self.shape == other.shape: 88 start = 0 89 end = 0 90 for i in xrange(len(self.parts)): 91 if self.vertical: 92 end += self.parts[i].shape[0] 93 if op == self.__ops.add: 94 self.parts[i] += other[start:end, :] 95 elif op == self.__ops.sub: 96 self.parts[i] -= other[start:end, :] 97 elif op == self.__ops.mul: 98 self.parts[i] *= other[start:end, :] 99 else: 100 raise ValueError("Operator not yet implemented!") 101 else: 102 end += self.parts[i].shape[1] 103 if op == self.__ops.add: 104 self.parts[i] += other[:, start:end] 105 elif op == self.__ops.sub: 106 self.parts[i] -= other[:, start:end] 107 elif op == self.__ops.mul: 108 self.parts[i] *= other[:, start:end] 109 else: 110 raise ValueError("Operator not yet implemented!") 111 start = end 112 else: 113 raise ValueError("Unknown type") 114 115 return self

116

117 - def __iadd__(self, other):

118 return self.__iop(other, self.__ops.add)

119

120 - def __isub__(self, other):

121 return self.__iop(other, self.__ops.sub)

122

123 - def __imul__(self, other):

124 return self.__iop(other, self.__ops.mul)

125

126 - def __idiv__(self, other):

127 if not np.isscalar(other): 128 raise ValueError("Operator not yet implemented for type!") 129 130 return self.__iop(other, self.__ops.div)

131

132 - def __itruediv__(self, other):

133 if not np.isscalar(other): 134 raise ValueError("Operator not yet implemented for type!") 135 136 return self.__iop(float(other), self.__ops.div)

137

138 - def __op(self, other, op):

139 new_parts = [0] * len(self.parts) 140 if np.isscalar(other): 141 for i in xrange(len(self.parts)): 142 if op == self.__ops.add: 143 new_parts[i] = self.parts[i] + other 144 elif op == self.__ops.sub: 145 new_parts[i] = self.parts[i] - other 146 elif op == self.__ops.mul: 147 new_parts[i] = self.parts[i] * other 148 elif op == self.__ops.div: 149 new_parts[i] = self.parts[i] * (1.0 / other) 150 else: 151 raise ValueError("Operator not yet implemented!") 152 elif isinstance(other, MultipartArray): 153 other_parts = other.get_parts() 154 for i in xrange(len(self.parts)): 155 if op == self.__ops.add: 156 new_parts[i] = self.parts[i] + other_parts[i] 157 elif op == self.__ops.sub: 158 new_parts[i] = self.parts[i] - other_parts[i] 159 elif op == self.__ops.mul: 160 new_parts[i] = self.parts[i] * other_parts[i] 161 else: 162 raise ValueError("Operator not yet implemented!") 163 elif self.shape == other.shape: 164 start = 0 165 end = 0 166 for i in xrange(len(self.parts)): 167 if self.vertical: 168 end += self.parts[i].shape[0] 169 if op == self.__ops.add: 170 new_parts[i] = self.parts[i] + other[start:end, :] 171 elif op == self.__ops.sub: 172 new_parts[i] = self.parts[i] - other[start:end, :] 173 elif op == self.__ops.mul: 174 new_parts[i] = self.parts[i] * other[start:end, :] 175 else: 176 raise ValueError("Operator not yet implemented!") 177 else: 178 end += self.parts[i].shape[1] 179 if op == self.__ops.add: 180 new_parts[i] = self.parts[i] + other[:, start:end] 181 elif op == self.__ops.sub: 182 new_parts[i] = self.parts[i] - other[:, start:end] 183 elif op == self.__ops.mul: 184 new_parts[i] = self.parts[i] * other[:, start:end] 185 else: 186 raise ValueError("Operator not yet implemented!") 187 start = end 188 else: 189 raise ValueError("Unknown type") 190 191 return MultipartArray(new_parts, vertical=self.vertical)

192

193 - def __add__(self, other):

194 return self.__op(other, self.__ops.add)

195

196 - def __sub__(self, other):

197 return self.__op(other, self.__ops.sub)

198

199 - def __mul__(self, other):

200 return self.__op(other, self.__ops.mul)

201

202 - def __div__(self, other):

203 if not np.isscalar(other): 204 raise ValueError("Operator not yet implemented for type!") 205 206 return self.__op(other, self.__ops.div)

207

208 - def __truediv__(self, other):

209 if not np.isscalar(other): 210 raise ValueError("Operator not yet implemented for type!") 211 212 return self.__op(float(other), self.__ops.div)

213

214 - def dot(self, other):

215 if self.vertical: 216 v = [0] * len(self.parts) 217 for i in xrange(len(self.parts)): 218 v[i] = self.parts[i].dot(other) 219 v = MultipartArray(v, vertical=True) 220 else: 221 v = np.zeros((self.shape[0], 1)) 222 if isinstance(other, MultipartArray): 223 other_parts = other.get_parts() 224 for i in xrange(len(self.parts)): 225 v += self.parts[i].dot(other_parts[i]) 226 elif self.shape[1] == other.shape[0]: 227 start = 0 228 end = 0 229 for i in xrange(len(self.parts)): 230 end += self.parts[i].shape[1] 231 v += self.parts[i].dot(other[start:end, :]) 232 start = end 233 else: 234 raise ValueError("Type or shape unknown") 235 236 return v

237

238 - def transpose(self):

239 new_parts = [0] * len(self.parts) 240 for i in xrange(len(self.parts)): 241 new_parts[i] = self.parts[i].transpose() 242 vertical = not self.vertical 243 244 return MultipartArray(new_parts, vertical=vertical)

245

246 - def _get_T(self):

247 return self.transpose()

248

249 - def _set_T(self, value):

250 raise AttributeError("attribute 'T' of 'MultipartArray' objects " \ 251 "is not writable")

252

253 - def _del_T(self):

254 raise AttributeError("attribute 'T' of 'MultipartArray' objects " \ 255 "is not writable")

256 257 T = property(_get_T, _set_T, _del_T, 'Transpose of the array.') 258

259 - def copy(self):

260 new_parts = [0] * len(self.parts) 261 for i in xrange(len(self.parts)): 262 new_parts[i] = self.parts[i].copy() 263 264 return MultipartArray(new_parts, vertical=self.vertical)

265

266 - def __str__(self):

267 string = "[" 268 if self.vertical: 269 for k in xrange(len(self.parts)): 270 for i in xrange(self.parts[k].shape[0]): 271 if i > 0 or k > 0: 272 string += ' ' 273 string += str(self.parts[k][i, :]) 274 if k < len(self.parts) - 1 \ 275 or i < self.parts[k].shape[0] - 1: 276 string += '\n' 277 if k < len(self.parts) - 1: 278 string += ' ' 279 string += '-' * (len(str(self.parts[k][i, :])) - 3) 280 string += "\n" 281 else: 282 for i in xrange(self.parts[0].shape[0]): 283 for k in xrange(len(self.parts)): 284 if k == 0 and i > 0: 285 string += ' ' 286 string += str(self.parts[k][i, :]) 287 288 if i < self.parts[len(self.parts) - 1].shape[0] - 1: 289 string += '\n' 290 291 string += "]" 292 293 return string

294

295 - def __repr__(self):

296 string = "MultipartArray(\n" + str(self.parts) 297 if self.vertical: 298 string += ")" 299 else: 300 string += ",\nvertical=" + str(self.vertical) + ")" 301 302 return string

303 304

305 -class Solver(object):

306 __metaclass__ = abc.ABCMeta 307

308 - def solve(A, b, eps=consts.TOLERANCE, max_iter=consts.MAX_ITER):

309 """Solves a linear system on the form 310 311 A.x = b, 312 313 for x. 314 315 Parameters 316 ---------- 317 A : A matrix with shape n-by-p. The coefficient matrix. 318 319 b : Numpy array, n-by-1. The right-hand-side vector. 320 """ 321 raise NotImplementedError('Abstract method "solve" must be ' 322 'specialised!')

323 324

325 -class SparseSolver(Solver):

326

327 - def solve(self, A, b, **kwargs):

328 """Solves linear systems on the form 329 330 A.x = d, 331 332 for x. 333 334 Parameters 335 ---------- 336 A : A sparse matrix with shape n-by-p. The coefficient matrix. 337 338 b : Numpy array, n-by-1. The right-hand-side vector. 339 340 Examples 341 -------- 342 >>> import numpy as np 343 >>> import scipy.sparse as sparse 344 >>> import parsimony.utils.linalgs as linalgs 345 >>> np.random.seed(42) 346 >>> 347 >>> n = 10 348 >>> a = np.random.rand(n); a[-1] = 0.0 349 >>> b = np.random.rand(n) 350 >>> c = np.random.rand(n); c[0] = 0.0 351 >>> A_ = np.random.rand(n, n) 352 >>> A_[A_ < 0.5] = 0.0 353 >>> A = sparse.csr_matrix(A_) 354 >>> d = np.random.rand(n, 1) 355 >>> 356 >>> solver = linalgs.SparseSolver() 357 >>> x = solver.solve(A, d) 358 >>> x_ = np.linalg.solve(A.toarray(), d) 359 >>> round(np.linalg.norm(x - x_), 15) 360 1e-15 361 >>> 362 >>> import time 363 >>> n = 100 364 >>> a = np.random.rand(n); a[-1] = 0.0 365 >>> b = np.random.rand(n) 366 >>> c = np.random.rand(n); c[0] = 0.0 367 >>> A_ = np.random.rand(n, n) 368 >>> A_[A_ < 0.5] = 0.0 369 >>> A = sparse.csr_matrix(A_) 370 >>> d = np.random.rand(n, 1) 371 >>> 372 >>> t = time.time() 373 >>> x = solver.solve(A, d) 374 >>> print "Time:", time.time() - t # doctest: +SKIP 375 >>> 376 >>> t = time.time() 377 >>> x_ = np.linalg.solve(A.toarray(), d) 378 >>> print "Time:", time.time() - t # doctest: +SKIP 379 >>> np.linalg.norm(x - x_) < 5e-14 380 True 381 >>> 382 >>> n = 1000 383 >>> a = np.random.rand(n); a[-1] = 0.0 384 >>> b = np.random.rand(n) 385 >>> c = np.random.rand(n); c[0] = 0.0 386 >>> A_ = np.random.rand(n, n) 387 >>> A_[A_ < 0.5] = 0.0 388 >>> A = sparse.csr_matrix(A_) 389 >>> d = np.random.rand(n, 1) 390 >>> 391 >>> t = time.time() 392 >>> x = solver.solve(A, d) 393 >>> print "Time:", time.time() - t # doctest: +SKIP 394 >>> 395 >>> t = time.time() 396 >>> x_ = np.linalg.solve(A.toarray(), d) 397 >>> print "Time:", time.time() - t # doctest: +SKIP 398 >>> 399 >>> np.linalg.norm(x - x_) < 5e-11 400 True 401 """ 402 n, p = A.shape 403 404 x = sparse.linalg.spsolve(A, b) 405 406 return x.reshape((n, 1))

407 408

409 -class TridiagonalSolver(Solver):

410

411 - def solve(self, A, d, **kwargs):

412 """Solves linear systems with a tridiagonal coefficient matrix. 413 414 A solver that uses the Thomas algorithm (the Tridiagonal matrix 415 algorithm) for systems on the form 416 417 0 418 [b c ] [x] [d] 419 A.x = [a b c ] [x] = [d] = d. 420 [ a b c] [x] [d] 421 [ a b] [x] [d] 422 0 423 424 Parameters 425 ---------- 426 A : A sparse diagonal matrix (dia format) with shape n-by-p. The 427 coefficient matrix. 428 429 b : Numpy array, n-by-1. The right-hand-side vector. 430 431 Examples 432 -------- 433 >>> import numpy as np 434 >>> import scipy.sparse as sparse 435 >>> import parsimony.utils.linalgs as linalgs 436 >>> np.random.seed(42) 437 >>> 438 >>> n = 10 439 >>> a = np.random.rand(n); a[-1] = 0.0 440 >>> b = np.random.rand(n) 441 >>> c = np.random.rand(n); c[0] = 0.0 442 >>> abc = np.vstack((a, b, c)) 443 >>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n)) 444 >>> d = np.random.rand(n, 1) 445 >>> 446 >>> solver = linalgs.TridiagonalSolver() 447 >>> x = solver.solve(A, d) 448 >>> print x 449 [[ -1.84339326] 450 [ 4.62737333] 451 [-12.41571989] 452 [ 16.38029815] 453 [ 14.38143172] 454 [-14.58969243] 455 [ 6.21233944] 456 [ 1.34271395] 457 [ -1.63358708] 458 [ 4.88318651]] 459 >>> x_ = np.linalg.solve(A.toarray(), d) 460 >>> print x_ 461 [[ -1.84339326] 462 [ 4.62737333] 463 [-12.41571989] 464 [ 16.38029815] 465 [ 14.38143172] 466 [-14.58969243] 467 [ 6.21233944] 468 [ 1.34271395] 469 [ -1.63358708] 470 [ 4.88318651]] 471 >>> np.linalg.norm(x - x_) < 5e-14 472 True 473 >>> 474 >>> import time 475 >>> n = 100 476 >>> a = np.random.rand(n); a[-1] = 0.0 477 >>> b = np.random.rand(n) 478 >>> c = np.random.rand(n); c[0] = 0.0 479 >>> abc = np.vstack((a, b, c)) 480 >>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n)) 481 >>> d = np.random.rand(n, 1) 482 >>> 483 >>> t = time.time() 484 >>> x = solver.solve(A, d) 485 >>> print "Time:", time.time() - t # doctest: +SKIP 486 >>> 487 >>> t = time.time() 488 >>> x_ = np.linalg.solve(A.toarray(), d) 489 >>> print "Time:", time.time() - t # doctest: +SKIP 490 >>> 491 >>> np.linalg.norm(x - x_) < 5e-12 492 True 493 >>> 494 >>> n = 1000 495 >>> a = np.random.rand(n); a[-1] = 0.0 496 >>> b = np.random.rand(n) 497 >>> c = np.random.rand(n); c[0] = 0.0 498 >>> abc = np.vstack((a, b, c)) 499 >>> A = sparse.dia_matrix((abc, [-1, 0, 1]), shape=(n, n)) 500 >>> d = np.random.rand(n, 1) 501 >>> 502 >>> t = time.time() 503 >>> x = solver.solve(A, d) 504 >>> print "Time:", time.time() - t # doctest: +SKIP 505 >>> 506 >>> t = time.time() 507 >>> x_ = np.linalg.solve(A.toarray(), d) 508 >>> print "Time:", time.time() - t # doctest: +SKIP 509 >>> 510 >>> np.linalg.norm(x - x_) < 5e-9 511 True 512 """ 513 # TODO: Put in compiled code for speed. 514 515 if not sparse.isspmatrix_dia(A): 516 A = A.todia() 517 518 abc = A.data 519 a = abc[0, :] 520 b = abc[1, :] 521 c = abc[2, :] 522 523 if abc.dtype != np.float: 524 a = np.asarray(a, np.float) 525 b = np.asarray(b, np.float) 526 c = np.asarray(c, np.float) 527 528 n = len(a) 529 x = np.zeros(n) 530 531 # Decomposition and forward substitution. 532 c_ = np.zeros(n) 533 d_ = np.zeros(n) 534 i = 0 535 if abs(b[i]) < consts.TOLERANCE: 536 # TODO: Do this instead: In this case x0 is found trivially and we 537 # recurse to a problem of order n-1. 538 solver = SparseSolver() 539 return solver.solve(A, d) 540 c_[i + 1] = c[i + 1] / b[i] 541 d_[i] = d[i] / b[i] 542 for i in range(1, n - 1): 543 i_1 = i - 1 544 den = (b[i] - a[i_1] * c_[i]) 545 if abs(den) < consts.TOLERANCE: # We cannot handle this case! 546 # TODO: Use algorithm for banded matrices instead! 547 solver = SparseSolver() 548 return solver.solve(A, d) 549 c_[i + 1] = c[i + 1] / den 550 d_[i] = (d[i] - a[i_1] * d_[i_1]) / den 551 i = n - 1 552 d_[i] = (d[i] - a[i - 1] * d_[i - 1]) / (b[i] - a[i - 1] * c_[i]) 553 554 # Back substitution. 555 i = n - 1 556 x[i] = d_[i] 557 for i in reversed(xrange(n - 1)): 558 x[i] = d_[i] - c_[i + 1] * x[i + 1] 559 560 return x.reshape((n, 1))

561 562 if __name__ == "__main__": 563 import doctest 564 doctest.testmod() 565

Source Code for Module parsimony.utils.linalgs