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
17 EULER_MASCHERONI = 0.57721566490153286060651209008240243104215933593992
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
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
59 """
60 Return the sign of the input.
61 """
62 return numpy.sign(x)
63
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
336
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
353 """
354 Inverse digamma function
355 """
356 from numpy import exp
357 from scipy.special import digamma
358
359
360 if y < -5 / 3. - EULER_MASCHERONI:
361 x = -1 / (EULER_MASCHERONI + y)
362 else:
363 x = 0.5 + exp(y)
364
365
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
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
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
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
414 return 0.5 * numpy.dot(x[1:] - x[:-1], y[1:] + y[:-1])
415
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
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
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
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
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
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
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
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
544
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
612
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
624
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
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
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
651
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
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
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
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
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
722 if self._matrix is not None:
723 return len(self._matrix)
724 else:
725 return len(self._inverse_matrix)
726
728 self._matrix_updated = False
729
731 self._inverse_matrix_updated = False
732
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
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
774 if self._inverse_matrix_updated == False:
775 self._update_inverse_matrix()
776
777 return self._inverse_matrix.copy()
778 @inverse.setter
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
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