csb.numeric

1 """ 2 Low level numeric / math utility functions. 3 """ 4 5 import sys 6 import math 7 import numpy 8 9 10 EXP_MIN = -308 11 EXP_MAX = +709 12 13 LOG_MIN = 1e-308 14 LOG_MAX = 1e+308 15 16 ## Euler-Mascheroni constant 17 EULER_MASCHERONI = 0.57721566490153286060651209008240243104215933593992

18 19 20 -def log(x, x_min=LOG_MIN, x_max=LOG_MAX):

21 """ 22 Safe version of log, clips argument such that overflow does not occur. 23 24 @param x: input 25 @type x: numpy array or float or int 26 27 @param x_min: lower value for clipping 28 @type x_min: float 29 30 @param x_max: upper value for clipping 31 @type x_max: float 32 """ 33 34 x_min = max(x_min, LOG_MIN) 35 x_max = min(x_max, LOG_MAX) 36 37 return numpy.log(numpy.clip(x, x_min, x_max))

38

39 -def exp(x, x_min=EXP_MIN, x_max=EXP_MAX):

40 """ 41 Safe version of exp, clips argument such that overflow does not occur. 42 43 @param x: input 44 @type x: numpy array or float or int 45 46 @param x_min: lower value for clipping 47 @type x_min: float 48 49 @param x_max: upper value for clipping 50 @type x_max: float 51 """ 52 53 x_min = max(x_min, EXP_MIN) 54 x_max = min(x_max, EXP_MAX) 55 56 return numpy.exp(numpy.clip(x, x_min, x_max))

57

58 -def sign(x):

59 """ 60 Return the sign of the input. 61 """ 62 return numpy.sign(x)

63

64 -def isreal(x, tol=1e-14):

65 """ 66 Check if input array has no imaginary part. 67 68 @param x: input array 69 @type x: numpy array 70 71 @param tol: tolerance to check for equality zero 72 @type tol: float 73 """ 74 return not hasattr(x, 'real') or abs(x.imag) < tol

75

76 -def log_sum_exp(x, axis=0):

77 """ 78 Return the logarithm of the sum of exponentials. 79 80 @type x: Numpy array 81 """ 82 xmax = x.max(axis) 83 return log(exp(x - xmax).sum(axis)) + xmax

84

85 -def log_sum_exp_accumulate(x, axis=0):

86 """ 87 Return the logarithm of the accumulated sums of exponentials. 88 89 @type x: Numpy array 90 """ 91 xmax = x.max(axis) 92 return log(numpy.add.accumulate(exp(x - xmax), axis)) + xmax

93

94 -def radian2degree(x):

95 """ 96 Convert radians angles to torsion angles. 97 98 @param x: radian angle 99 @return: torsion angle of x 100 """ 101 x = x % (2 * numpy.pi) 102 numpy.putmask(x, x > numpy.pi, x - 2 * numpy.pi) 103 return x * 180. / numpy.pi

104

105 -def degree2radian(x):

106 """ 107 Convert randian angles to torsion angles. 108 109 @param x: torsion angle 110 @return: radian angle of x 111 """ 112 numpy.putmask(x, x < 0., x + 360.) 113 return x * numpy.pi / 180.

114

115 -def euler_angles(r):

116 """ 117 Calculate the euler angles from a three dimensional rotation matrix. 118 119 @param r: 3x3 Rotation matrix 120 """ 121 a = numpy.arctan2(r[2, 1], r[2, 0]) % (2 * numpy.pi) 122 b = numpy.arctan2((r[2, 0] + r[2, 1]) / (numpy.cos(a) + numpy.sin(a)), r[2, 2]) % (2 * numpy.pi) 123 c = numpy.arctan2(r[1, 2], -r[0, 2]) % (2 * numpy.pi) 124 125 return a, b, c

126

127 -def euler(a, b, c):

128 """ 129 Calculate a three dimensional rotation matrix from the euler angles. 130 131 @param a: alpha, angle between the x-axis and the line of nodes 132 @param b: beta, angle between the z axis of the different coordinate systems 133 @param c: gamma, angle between the line of nodes and the X-axis 134 """ 135 from numpy import cos, sin, array 136 137 ca, cb, cc = cos(a), cos(b), cos(c) 138 sa, sb, sc = sin(a), sin(b), sin(c) 139 140 return array([[ cc * cb * ca - sc * sa, cc * cb * sa + sc * ca, -cc * sb], 141 [-sc * cb * ca - cc * sa, -sc * cb * sa + cc * ca, sc * sb], 142 [ sb * ca, sb * sa, cb ]])

143

144 -def rotation_matrix(axis, angle):

145 """ 146 Calculate a three dimensional rotation matrix for a rotation around 147 the given angle and axis. 148 149 @type axis: (3,) numpy array 150 @param angle: angle in radians 151 @type angle: float 152 153 @rtype: (3,3) numpy.array 154 """ 155 axis = numpy.asfarray(axis) / norm(axis) 156 assert axis.shape == (3,) 157 158 c = math.cos(angle) 159 s = math.sin(angle) 160 161 r = (1.0 - c) * numpy.outer(axis, axis) 162 r.flat[[0,4,8]] += c 163 r.flat[[5,6,1]] += s * axis 164 r.flat[[7,2,3]] -= s * axis 165 166 return r

167

168 -def axis_and_angle(r):

169 """ 170 Calculate axis and angle of rotation from a three dimensional 171 rotation matrix. 172 173 @param r: 3x3 Rotation matrix 174 175 @return: axis unit vector as (3,) numpy.array and angle in radians as float 176 @rtype: tuple 177 """ 178 from scipy.linalg import logm 179 180 B = logm(r).real 181 a = numpy.array([B[1,2], -B[0,2], B[0,1]]) 182 n = norm(a) 183 184 return a / n, n

185

186 -def norm(x):

187 """ 188 Calculate the Eucledian norm of a d-dimensional vector. 189 190 @param x: vector (i.e. rank one array) 191 @return: length of vector 192 """ 193 return numpy.linalg.norm(x)

194

195 -def reverse(array, axis=0):

196 """ 197 Reverse the order of elements in an array. 198 """ 199 from numpy import take, arange 200 return take(array, arange(array.shape[axis] - 1, -1, -1), axis)

201

202 -def polar(x):

203 """ 204 Polar coordinate representation of a d-dimensional vector. 205 206 @param x: vector (i.e. rank one array) 207 @return: polar coordinates (radius and polar angles) 208 """ 209 210 (d,) = x.shape 211 phi = numpy.zeros(d) 212 for i in reversed(range(1, d)): 213 phi[i - 1] = numpy.arctan2(x[i] / numpy.cos(phi[i]), x[i - 1]) 214 215 return numpy.array([norm(x)] + phi[:-1].tolist())

216

217 -def from_polar(x):

218 """ 219 Reconstruct d-dimensional vector from polar coordinates. 220 221 @param x: vector (i.e. rank one array) 222 @return: position in d-dimensional space 223 """ 224 225 (d,) = x.shape 226 227 c = numpy.cos(x[1:]) 228 s = numpy.sin(x[1:]) 229 r = x[0] 230 231 x = numpy.zeros(d) 232 x[0] = r 233 234 for i in range(d - 1): 235 236 x[i + 1] = x[i] * s[i] 237 x[i] *= c[i] 238 239 return x

240

241 -def polar3d(x):

242 """ 243 Polar coordinate representation of a three-dimensional vector. 244 245 @param x: vector (i.e. rank one array) 246 @return: polar coordinates (radius and polar angles) 247 """ 248 249 if x.shape != (3,): 250 raise ValueError(x) 251 252 r = norm(x) 253 theta = numpy.arccos(x[2] / r) 254 phi = numpy.arctan2(x[1], x[0]) 255 256 return numpy.array([r, theta, phi])

257

258 -def from_polar3d(x):

259 """ 260 Reconstruct 3-dimensional vector from polar coordinates. 261 262 @param x: vector (i.e. rank one array) 263 @return: position in 3-dimensional space 264 """ 265 assert x.shape == (3,) 266 267 r, theta, phi = x[:] 268 s = numpy.sin(theta) 269 c = numpy.cos(theta) 270 S = numpy.sin(phi) 271 C = numpy.cos(phi) 272 273 return r * numpy.array([s * C, s * S, c])

274

275 -def dihedral_angle(a, b, c, d):

276 """ 277 Calculate the dihedral angle between 4 vectors, 278 representing 4 connected points. The angle is in range [-180, 180]. 279 280 @param a: the four points that define the dihedral angle 281 @type a: array 282 283 @return: angle in [-180, 180] 284 """ 285 286 v = b - c 287 m = numpy.cross((a - b), v) 288 m /= norm(m) 289 n = numpy.cross((d - c), v) 290 n /= norm(n) 291 292 c = numpy.dot(m, n) 293 s = numpy.dot(numpy.cross(n, m), v) / norm(v) 294 295 angle = math.degrees(math.atan2(s, c)) 296 297 if angle > 0: 298 return numpy.fmod(angle + 180, 360) - 180 299 else: 300 return numpy.fmod(angle - 180, 360) + 180

301

302 -def psi(x):

303 """ 304 Digamma function 305 """ 306 from numpy import inf, log, sum, exp 307 308 coef = [-1. / 12., 1. / 120., -1. / 252., 1. / 240., -1. / 132., 309 691. / 32760., -1. / 12.] 310 311 if x == 0.: 312 return -inf 313 elif x < 0.: 314 raise ValueError('not defined for negative values') 315 elif x < 6.: 316 return psi(x + 1) - 1. / x 317 else: 318 logx = log(x) 319 res = logx - 0.5 / x 320 res += sum([coef[i] * exp(-2 * (i + 1) * logx) for i in range(7)]) 321 return res

322

323 -def approx_psi(x):

324 from numpy import log, clip, where 325 326 if type(x) == numpy.ndarray: 327 y = 0. * x 328 y[where(x < 0.6)] = -EULER_MASCHERONI - 1. / clip(x[where(x < 0.6)], 1e-154, 1e308) 329 y[where(x >= 0.6)] = log(x[where(x >= 0.6)] - 0.5) 330 return y 331 else: 332 if x < 0.6: 333 return -EULER_MASCHERONI - 1 / clip(x, 1e-154, 1e308) 334 else: 335 return log(x - 0.5)

336

337 -def d_approx_psi(x):

338 from numpy import clip, where 339 340 if type(x) == numpy.ndarray: 341 y = 0. * x 342 y[where(x < 0.6)] = 1. / clip(x[where(x < 0.6)], 1e-154, 1e308) ** 2 343 y[where(x >= 0.6)] = 1. / (x[where(x >= 0.6)] - 0.5) 344 return y 345 else: 346 347 if x < 0.6: 348 return 1 / clip(x, 1e-154, 1e308) ** 2 349 else: 350 return 1 / (x - 0.5)

351

352 -def inv_psi(y, tol=1e-10, n_iter=100, psi=psi):

353 """ 354 Inverse digamma function 355 """ 356 from numpy import exp 357 from scipy.special import digamma 358 ## initial value 359 360 if y < -5 / 3. - EULER_MASCHERONI: 361 x = -1 / (EULER_MASCHERONI + y) 362 else: 363 x = 0.5 + exp(y) 364 365 ## Newton root finding 366 367 for dummy in range(n_iter): 368 369 z = digamma(x) - y 370 371 if abs(z) < tol: 372 break 373 374 x -= z / d_approx_psi(x) 375 376 return x

377

378 -def log_trapezoidal(log_y, x=None):

379 """ 380 Compute the logarithm of the 1D integral of x, using trepezoidal approximation. 381 Assumes x is monotonically increasing. 382 """ 383 if x is not None: 384 log_x_diff = log(x[1:] - x[:-1]) 385 else: 386 log_x_diff = 0. 387 388 log_y_add = log_sum_exp(numpy.vstack((log_y[:-1], log_y[1:])), 0) 389 390 return log_sum_exp(log_y_add + log_x_diff) - log(2)

391

392 -def log_midpoint_rule_2d(log_f, x, y):

393 x_delta = x[:-1] - x[1:] 394 y_delta = y[:-1] - y[1:] 395 396 z = numpy.array([log_f[:, 1:] , log_f[:, :-1]]) 397 398 y_hat = log_sum_exp(z.reshape((2, -1)), 0) 399 y_hat = numpy.reshape(y_hat, (len(x), len(y) - 1)) 400 y_hat += log(y_delta) - log(2) 401 402 return log_sum_exp(y_hat + log(x_delta)) - log(2.)

403

404 -def log_trapezoidal_2d(log_f, x=None, y=None):

405 """ 406 Compute the logarithm of the 1D integral of x, using trepezoidal approximation. 407 Assumes x and y is monotonically increasing. 408 """ 409 int_y = numpy.array([log_trapezoidal(log_f[i, :], y) for i in range(len(y))]) 410 411 return log_trapezoidal(int_y, x)

412

413 -def trapezoidal(x, y):

414 return 0.5 * numpy.dot(x[1:] - x[:-1], y[1:] + y[:-1])

415

416 -def trapezoidal_2d(f):

417 """ 418 Approximate the integral of f from a to b in two dimensions, 419 using trepezoidal approximation. 420 421 @param f: 2D numpy array of function values at equally spaces positions 422 @return: approximation of the definit integral 423 """ 424 425 I = f[0, 0] + f[-1, -1] + f[0, -1] + f[-1, 0] 426 I += 2 * (f[1:-1, (0, -1)].sum() + f[(0, -1), 1:-1].sum()) 427 I += 4 * f[1:-1, 1:-1].sum() 428 429 return I / 4.

430

431 -def simpson_2d(f):

432 """ 433 Approximate the integral of f from a to b in two dimensions, 434 using Composite Simpson's rule. 435 436 @param f: 2D numpy array of function values 437 @return: approximation of the definit integral 438 """ 439 440 n = int((f.shape[0] - 1) / 2) 441 i = 2 * numpy.arange(1, n + 1) - 1 442 j = 2 * numpy.arange(1, n) 443 444 I = f[0, 0] + f[-1, -1] + f[0, -1] + f[-1, 0] 445 I += 4 * (f[0, i].sum() + f[-1, i].sum() + f[0, j].sum() + f[-1, j].sum()) 446 I += 4 * (f[i, 0].sum() + f[i, -1].sum() + f[j, 0].sum() + f[j, -1].sum()) 447 I += 16 * f[i][:, i].sum() + 8 * (f[i][:, j].sum() + f[j][:, i].sum()) 448 I += 4 * f[j][:, j].sum() 449 450 return I / 9.

451

452 -def pad(x, s):

453 """ 454 Add layers of zeros around grid. 455 """ 456 s = numpy.array(s) - 1 457 y = numpy.zeros(numpy.array(x.shape) + s) 458 s /= 2 459 slices = [slice(s[i], -s[i]) for i in range(len(s))] 460 y[slices] = x 461 462 return y

463

464 -def trim(x, s):

465 """ 466 Remove additional layers. 467 """ 468 s = numpy.array(s) - 1 469 s /= 2 470 471 slices = [slice(s[i], -s[i]) for i in range(len(s))] 472 473 return x[slices]

474

475 -def zerofill(x, s):

476 477 y = numpy.zeros(s) 478 slices = [slice(-x.shape[i], None) for i in range(len(s))] 479 y[slices] = x 480 481 return y

482

483 -def convolve(x, f):

484 485 from numpy import fft, all 486 487 sx = numpy.array(x.shape) 488 sf = numpy.array(f.shape) 489 490 if not all(sx >= sf): return convolve(f, x) 491 492 y = fft.ifftn(fft.fftn(x) * fft.fftn(f, sx)).real 493 slices = [slice(sf[i] - 1, sx[i]) for i in range(len(sf))] 494 495 return y[slices]

496

497 -def correlate(x, y):

498 499 from numpy import fft 500 501 sx = numpy.array(x.shape) 502 sy = numpy.array(y.shape) 503 504 if (sx >= sy).sum(): 505 506 slices = [slice(None, sx[i] - sy[i] + 1) for i in range(len(sx))] 507 508 X = fft.fftn(x) 509 Y = fft.fftn(zerofill(y, sx)) 510 511 else: 512 513 sf = sx + sy - 1 514 slices = [slice(None, sf[i]) for i in range(len(sf))] 515 516 X = fft.fftn(x, sf) 517 Y = fft.fftn(zerofill(y, sf), sf) 518 519 return fft.ifftn(X.conjugate() * Y)[slices].real

520

521 522 -def gower_matrix(X):

523 """ 524 Gower, J.C. (1966). Some distance properties of latent root 525 and vector methods used in multivariate analysis. 526 Biometrika 53: 325-338 527 528 @param X: ensemble coordinates 529 @type X: (m,n,k) numpy.array 530 531 @return: symmetric dissimilarity matrix 532 @rtype: (n,n) numpy.array 533 """ 534 X = numpy.asarray(X) 535 536 B = sum(numpy.dot(x, x.T) for x in X) / float(len(X)) 537 b = B.mean(1) 538 bb = b.mean() 539 540 return (B - numpy.add.outer(b, b)) + bb

541

542 -class MatrixInitError(Exception):

543 pass

544

545 -class InvertibleMatrix(object):

546 """ 547 Matrix object which is intended to save time in MCMC sampling algorithms involving 548 repeated integration of Hamiltonian equations of motion and frequent draws from 549 multivariate normal distributions involving mass matrices as covariance matrices. 550 It can be initialized either with the matrix one wants to use or its inverse. 551 The main feature compared to standard numpy matrices / arrays is that it has also 552 a property "inverse", which gives the inverse matrix. If the matrix (its inverse) is 553 changed, the inverse (regular matrix) is calculated only when needed. This avoids costly 554 matrix inversions. 555 556 @param matrix: matrix-like object with whose values the Matrix object is supposed to 557 hold 558 @type matrix: invertible (n,n)-shaped numpy.ndarray 559 560 @param inverse_matrix: matrix-like object with whose inverse the Matrix object is supposed 561 to hold 562 @type inverse_matrix: invertible (n,n)-shaped numpy.ndarray 563 """ 564

565 - def __init__(self, matrix=None, inverse_matrix=None):

566 567 if (matrix is None and inverse_matrix is None): 568 raise MatrixInitError("At least one matrix argument has to be specified") 569 570 self._matrix = None 571 self._inverse_matrix = None 572 573 self._matrix_updated = False 574 self._inverse_matrix_updated = False 575 576 self._is_unity_multiple = False 577 578 if matrix is not None and inverse_matrix is not None: 579 if type(matrix) != numpy.ndarray or type(inverse_matrix) != numpy.ndarray: 580 raise TypeError("Arguments have to be of type numpy.ndarray!") 581 matrix = matrix.copy() 582 inverse_matrix = inverse_matrix.copy() 583 self._check_equal_shape(matrix, inverse_matrix) 584 self._matrix = matrix 585 self._inverse_matrix = inverse_matrix 586 self._is_unity_multiple = self._check_unity_multiple(self._matrix) 587 self._matrix_updated = True 588 self._inverse_matrix_updated = True 589 else: 590 if matrix is not None: 591 if type(matrix) != numpy.ndarray: 592 raise TypeError("Arguments have to be of type numpy.ndarray!") 593 matrix = matrix.copy() 594 self._check_square(matrix) 595 self._matrix = matrix 596 self._matrix_updated = True 597 self._inverse_matrix_updated = False 598 self._is_unity_multiple = self._check_unity_multiple(self._matrix) 599 else: 600 if type(inverse_matrix) != numpy.ndarray: 601 raise TypeError("Arguments have to be of type numpy.ndarray!") 602 inverse_matrix = inverse_matrix.copy() 603 self._check_square(inverse_matrix) 604 self._inverse_matrix = inverse_matrix 605 self._matrix_updated = False 606 self._inverse_matrix_updated = True 607 self._is_unity_multiple = self._check_unity_multiple(self._inverse_matrix)

608

609 - def _check_diagonal(self, matrix):

610 611 i, j = numpy.nonzero(matrix) 612 613 return numpy.array_equal(i, j)

614

615 - def _check_unity_multiple(self, matrix):

616 617 diagonal_elements_equal = numpy.all([matrix[i][i] == matrix[i+1][i+1] 618 for i in range(len(matrix) - 1)]) 619 620 return self._check_diagonal(matrix) and diagonal_elements_equal

621

622 - def _check_square(self, matrix):

623 624 if matrix.shape[0] != matrix.shape[1]: 625 raise ValueError("Matrix " + matrix.__name__ + " must be a square matrix!")

626

627 - def _check_equal_shape(self, matrix1, matrix2):

628 629 if not numpy.all(matrix1.shape == matrix2.shape): 630 raise ValueError("Matrices " + matrix1.__name__ + " and " + matrix2.__name__ + 631 " must have equal shape!")

632

633 - def __getitem__(self, loc):

634 if self._matrix_updated == False: 635 if self.is_unity_multiple: 636 self._matrix = numpy.diag(1. / self._inverse_matrix.diagonal()) 637 else: 638 self._matrix = numpy.linalg.inv(self._inverse_matrix) 639 self._matrix_updated = True 640 641 return self._matrix[loc]

642

643 - def __setitem__(self, i, value):

644 if type(value) != numpy.ndarray: 645 raise TypeError("Arguments have to be of type numpy.ndarray!") 646 self._matrix[i] = numpy.array(value) 647 self._matrix_updated = True 648 self._inverse_matrix_updated = False 649 self._is_unity_multiple = self._check_unity_multiple(self._matrix)

650

651 - def __array__(self, dtype=None):

652 return self._matrix

653

654 - def __mul__(self, value):

655 if type(value) in (float, int): 656 v = float(value) 657 if self._matrix_updated: 658 return InvertibleMatrix(v * self._matrix) 659 else: 660 return InvertibleMatrix(inverse_matrix = self._inverse_matrix / v) 661 662 else: 663 raise ValueError("Only float and int are supported for multiplication!")

664 665 __rmul__ = __mul__ 666

667 - def __truediv__(self, value):

668 if type(value) in (float, int): 669 v = float(value) 670 if self._matrix_updated: 671 return InvertibleMatrix(self._matrix / v) 672 else: 673 return InvertibleMatrix(inverse_matrix=self._inverse_matrix * v) 674 675 else: 676 raise ValueError("Only float and int are supported for division!")

677 678 __div__ = __truediv__ 679

680 - def __imul__(self, value):

681 if type(value) in (float, int): 682 if self._matrix_updated: 683 self._matrix *= float(value) 684 self._deprecate_inverse_matrix() 685 else: 686 self._inverse_matrix /= float(value) 687 self._inverse_matrix_updated = True 688 self._deprecate_matrix() 689 690 return self 691 else: 692 raise ValueError("Only float and int are supported for multiplication!")

693

694 - def __itruediv__(self, value):

695 if type(value) in (float, int): 696 if self._matrix_updated: 697 self._matrix /= float(value) 698 self._matrix_updated = True 699 self._deprecate_inverse_matrix() 700 else: 701 self._inverse_matrix *= float(value) 702 self._inverse_matrix_updated = True 703 self._deprecate_matrix() 704 705 return self 706 else: 707 raise ValueError("Only float and int are supported for division!")

708 709 __idiv__ = __itruediv__ 710

711 - def __eq__(self, other):

712 713 if self._matrix is not None and other._matrix is not None: 714 715 return self._matrix == other._matrix 716 717 if self._inverse_matrix is not None and other._inverse_matrix is not None: 718 719 return self._inverse_matrix == other._inverse_matrix

720

721 - def __str__(self):

722 if self._matrix is not None and self._inverse_matrix is not None: 723 return "csb.numeric.InvertibleMatrix object holding the following numpy matrix:\n"\ 724 +self._matrix.__str__() +"\n and its inverse:\n"+self._inverse_matrix.__str__() 725 else: 726 if self._matrix is None: 727 return "csb.numeric.InvertibleMatrix object holding only the inverse matrix:\n"\ 728 +self._inverse_matrix.__str__() 729 else: 730 return "csb.numeric.InvertibleMatrix object holding the matrix:\n"\ 731 +self._matrix.__str__()

732

733 - def __len__(self):

734 if self._matrix is not None: 735 return len(self._matrix) 736 else: 737 return len(self._inverse_matrix)

738

739 - def _deprecate_matrix(self):

740 self._matrix_updated = False

741

742 - def _deprecate_inverse_matrix(self):

743 self._inverse_matrix_updated = False

744

745 - def _update_matrix(self, matrix=None):

746 """ 747 Updates the _matrix field given a new matrix or by setting 748 it to the inverse of the _inverse_matrix field. 749 750 @param matrix: matrix-like object which the Matrix object 751 is supposed to represent 752 @type matrix: (n,n)-shaped numpy.ndarray or list 753 """ 754 755 if matrix is None: 756 self._matrix = numpy.linalg.inv(self._inverse_matrix) 757 else: 758 self._matrix = matrix 759 760 self._matrix_updated = True 761 self._deprecate_inverse_matrix()

762 763

764 - def _update_inverse_matrix(self, inverse=None):

765 """ 766 Updates the __inverse_matrix field given a new matrix or by setting 767 it to the inverse of the _matrix field. 768 769 @param inverse: matrix-like object which the Matrix object 770 is supposed to represent 771 @type inverse: (n,n)-shaped numpy.ndarray or list 772 """ 773 774 if inverse is None: 775 if self.is_unity_multiple: 776 self._inverse_matrix = numpy.diag(1. / self._matrix.diagonal()) 777 else: 778 self._inverse_matrix = numpy.linalg.inv(self._matrix) 779 else: 780 self._inverse_matrix = inverse 781 782 self._inverse_matrix_updated = True

783 784 @property

785 - def inverse(self):

786 if self._inverse_matrix_updated == False: 787 self._update_inverse_matrix() 788 789 return self._inverse_matrix.copy()

790 @inverse.setter

791 - def inverse(self, value):

792 if type(value) != numpy.ndarray: 793 raise TypeError("Arguments have to be of type numpy.ndarray!") 794 self._check_equal_shape(value, self._matrix) 795 self._update_inverse_matrix(numpy.array(value)) 796 self._deprecate_matrix() 797 self._is_unity_multiple = self._check_unity_multiple(self._inverse_matrix)

798 799 @property

800 - def is_unity_multiple(self):

801 """ 802 This property can be used to save computation time when drawing 803 from multivariate normal distributions with the covariance matrix 804 given by an instance of this object. By probing this property, 805 one can instead draw from normal distributions and rescale the samples 806 afterwards to avoid costly matrix inversions 807 """ 808 return self._is_unity_multiple

809

Source Code for Package csb.numeric