Package csb :: Package numeric
[frames] | no frames]

Source Code for Package 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(m) 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 return (matrix.T == matrix).all()
612
613 - def _check_unity_multiple(self, matrix):
614 615 diagonal_elements_equal = numpy.all([matrix[i][i] == matrix[i+1][i+1] 616 for i in range(len(matrix) - 1)]) 617 618 return self._check_diagonal(matrix) and diagonal_elements_equal
619
620 - def _check_square(self, matrix):
621 622 if matrix.shape[0] != matrix.shape[1]: 623 raise ValueError("Matrix " + matrix.__name__ + " must be a square matrix!")
624
625 - def _check_equal_shape(self, matrix1, matrix2):
626 627 if not numpy.all(matrix1.shape == matrix2.shape): 628 raise ValueError("Matrices " + matrix1.__name__ + " and " + matrix2.__name__ + 629 " must have equal shape!")
630
631 - def __getitem__(self, loc):
632 if self._matrix_updated == False: 633 if self.is_unity_multiple: 634 self._matrix = numpy.diag(1. / self._inverse_matrix.diagonal()) 635 else: 636 self._matrix = numpy.linalg.inv(self._inverse_matrix) 637 self._matrix_updated = True 638 639 return self._matrix[loc]
640
641 - def __setitem__(self, i, value):
642 if type(value) != numpy.ndarray: 643 raise TypeError("Arguments have to be of type numpy.ndarray!") 644 self._matrix[i] = numpy.array(value) 645 self._matrix_updated = True 646 self._inverse_matrix_updated = False 647 self._is_unity_multiple = self._check_unity_multiple(self._matrix)
648
649 - def __array__(self, dtype=None):
650 return self._matrix
651
652 - def __mul__(self, value):
653 if type(value) in (float, int): 654 v = float(value) 655 if self._matrix_updated: 656 return InvertibleMatrix(v * self._matrix) 657 else: 658 return InvertibleMatrix(inverse_matrix = self._inverse_matrix / v) 659 660 else: 661 raise ValueError("Only float and int are supported for multiplication!")
662 663 __rmul__ = __mul__ 664
665 - def __truediv__(self, value):
666 if type(value) in (float, int): 667 v = float(value) 668 if self._matrix_updated: 669 return InvertibleMatrix(self._matrix / v) 670 else: 671 return InvertibleMatrix(inverse_matrix=self._inverse_matrix * v) 672 673 else: 674 raise ValueError("Only float and int are supported for division!")
675 676 __div__ = __truediv__ 677
678 - def __imul__(self, value):
679 if type(value) in (float, int): 680 if self._matrix_updated: 681 self._matrix *= float(value) 682 self._deprecate_inverse_matrix() 683 else: 684 self._inverse_matrix /= float(value) 685 self._inverse_matrix_updated = True 686 self._deprecate_matrix() 687 688 return self 689 else: 690 raise ValueError("Only float and int are supported for multiplication!")
691
692 - def __itruediv__(self, value):
693 if type(value) in (float, int): 694 if self._matrix_updated: 695 self._matrix /= float(value) 696 self._matrix_updated = True 697 self._deprecate_inverse_matrix() 698 else: 699 self._inverse_matrix *= float(value) 700 self._inverse_matrix_updated = True 701 self._deprecate_matrix() 702 703 return self 704 else: 705 raise ValueError("Only float and int are supported for division!")
706 707 __idiv__ = __itruediv__ 708
709 - def __str__(self):
710 if self._matrix is not None and self._inverse_matrix is not None: 711 return "csb.numeric.InvertibleMatrix object holding the following numpy matrix:\n"\ 712 +self._matrix.__str__() +"\n and its inverse:\n"+self._inverse_matrix.__str__() 713 else: 714 if self._matrix is None: 715 return "csb.numeric.InvertibleMatrix object holding only the inverse matrix:\n"\ 716 +self._inverse_matrix.__str__() 717 else: 718 return "csb.numeric.InvertibleMatrix object holding the matrix:\n"\ 719 +self._matrix.__str__()
720
721 - def __len__(self):
722 if self._matrix is not None: 723 return len(self._matrix) 724 else: 725 return len(self._inverse_matrix)
726
727 - def _deprecate_matrix(self):
728 self._matrix_updated = False
729
731 self._inverse_matrix_updated = False
732
733 - def _update_matrix(self, matrix=None):
734 """ 735 Updates the _matrix field given a new matrix or by setting 736 it to the inverse of the _inverse_matrix field. 737 738 @param matrix: matrix-like object which the Matrix object 739 is supposed to represent 740 @type matrix: (n,n)-shaped numpy.ndarray or list 741 """ 742 743 if matrix is None: 744 self._matrix = numpy.linalg.inv(self._inverse_matrix) 745 else: 746 self._matrix = matrix 747 748 self._matrix_updated = True 749 self._deprecate_inverse_matrix()
750 751
752 - def _update_inverse_matrix(self, inverse=None):
753 """ 754 Updates the __inverse_matrix field given a new matrix or by setting 755 it to the inverse of the _matrix field. 756 757 @param inverse: matrix-like object which the Matrix object 758 is supposed to represent 759 @type inverse: (n,n)-shaped numpy.ndarray or list 760 """ 761 762 if inverse is None: 763 if self.is_unity_multiple: 764 self._inverse_matrix = numpy.diag(1. / self._matrix.diagonal()) 765 else: 766 self._inverse_matrix = numpy.linalg.inv(self._matrix) 767 else: 768 self._inverse_matrix = inverse 769 770 self._inverse_matrix_updated = True
771 772 @property
773 - def inverse(self):
774 if self._inverse_matrix_updated == False: 775 self._update_inverse_matrix() 776 777 return self._inverse_matrix.copy()
778 @inverse.setter
779 - def inverse(self, value):
780 if type(value) != numpy.ndarray: 781 raise TypeError("Arguments have to be of type numpy.ndarray!") 782 self._check_equal_shape(value, self._matrix) 783 self._update_inverse_matrix(numpy.array(value)) 784 self._deprecate_matrix() 785 self._is_unity_multiple = self._check_unity_multiple(self._inverse_matrix)
786 787 @property
788 - def is_unity_multiple(self):
789 """ 790 This property can be used to save computation time when drawing 791 from multivariate normal distributions with the covariance matrix 792 given by an instance of this object. By probing this property, 793 one can instead draw from normal distributions and rescale the samples 794 afterwards to avoid costly matrix inversions 795 """ 796 return self._is_unity_multiple
797