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(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