1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27 """
28 general purpose math methods
29 """
30
31 import numpy.oldnumeric as N
32 import random
33 import numpy.oldnumeric.random_array as RandomArray
34 import math
35
37 pass
38
40 """
41 cumulative sum of C{ a[0], a[0]+a[1], a[0]+a[1]+[a2], ... }
42 normalized by C{ N.sum( a ) }
43
44 @param a: array('f') or float
45 @type a: array
46
47 @return: float
48 @rtype: float
49 """
50 return N.add.accumulate( a ) / N.sum( a )
51
52
54 """
55 Variance, S{sigma}^2
56
57 @param x: data
58 @type x: array('f') or float
59 @param avg: use this average, otherwise calculated from x
60 @type avg: float OR None
61
62 @return: float
63 @rtype: float
64 """
65 if avg is None:
66 avg = N.average(x)
67
68 if len(x) == 1:
69 return 0.0
70
71 return N.sum(N.power(N.array(x) - avg, 2)) / (len(x) - 1.)
72
73
74 -def SD(x, avg = None):
75 """
76 Standard deviation, S{sigma}
77
78 @param x: data
79 @type x: array('f') or float
80 @param avg: use this average, otherwise calculated from x
81 @type avg: float OR None
82
83 @return: float
84 @rtype: float
85 """
86 return N.sqrt(variance(x, avg))
87
88
90 """
91 Weighted mean: Mean of data (x) weighted by (w).
92
93 @param x: X-D array with numbers
94 @type x: array
95 @param w: 1-D array of same length as x with weight factors
96 @type w: array
97
98 @return: array('f') or float
99 @rtype: array('f') or float
100 """
101 if w is None:
102 wx = x
103 else:
104 wx = [ x[i] * 1. * w[i] for i in range( len(x) ) ]
105
106 return N.sum(wx)/N.sum(w)
107
108
110 """
111 Variance of weighted (w) data (x).
112
113 @param x: X-D array with numbers
114 @type x: array
115 @param w: 1-D array of same length as x with weight factors
116 @type w: array
117
118 @return: array('f') or float
119 @rtype: array('f') or float
120 """
121 wm = wMean(x,w)
122 return ( N.sum(w) / ( (N.sum(w)**2-N.sum(w**2)) ) ) * N.sum(w*(x-wm)**2)
123
124
126 """
127 Standard deviation of weighted data.
128
129 @param x: X-D array with numbers
130 @type x: array
131 @param w: 1-D array of same length as x with weight factors
132 @type w: array
133
134 @return: array('f') or float
135 @rtype: array('f') or float
136 """
137 return N.sqrt( wVar(x, w) )
138
139
141 """
142 Cross product between two coordinates.
143
144 @param v: coordinate vector
145 @type v: list or array
146 @param w: coordinate vector
147 @type w: list or array
148
149 @return: array('f') or float
150 @rtype: array('f') or float
151 """
152 x = v[1]*w[2] - v[2]*w[1]
153 y = v[2]*w[0] - v[0]*w[2]
154 z = v[0]*w[1] - v[1]*w[0]
155
156 return N.array( (x, y, z ) )
157
158
160 """
161 Collect att the values above the diagonal in a square
162 matrix.
163
164 @param pw_m: symmetric square matrix
165 @type pw_m: 2-D array
166
167 @return: raveled list of 'unique' values without diagonal
168 @rtype: list
169 """
170 r = []
171 for i in range( 0, len( pw_m ) ):
172
173 for j in range( i+1, len( pw_m ) ):
174
175 r.append( pw_m[i,j] )
176
177 return r
178
179
181 """
182 Compare 2 arrays or lists of numbers for equality.
183
184 @param a: first array (multi-dimensional is supported)
185 @type a: array / list
186 @param b: second array (multi-dimensional is supported)
187 @type b: array / list
188
189 @return: 1 if array/list a equals array/list b
190 @rtype: 1|0
191 """
192 if a is None or b is None:
193 return a is b
194
195 if len(a) != len(b):
196 return 0
197
198 if type(a) is list: a = N.array( a )
199 if type(b) is list: b = N.array( b )
200
201 a = N.ravel( a )
202 b = N.ravel( b )
203
204 return N.sum( a==b ) == len(a)
205
206
208 """
209 Pairwise distances between two arrays.
210
211 @param u: first array
212 @type u: array
213 @param v: second array
214 @type v: array
215
216 @return: Numeric.array( len(u) x len(v) ) of double
217 @rtype: array
218 """
219 diag1 = N.diagonal( N.dot( u, N.transpose(u) ) )
220 diag2 = N.diagonal( N.dot( v, N.transpose(v) ) )
221 dist = -N.dot( v,N.transpose(u) )\
222 -N.transpose( N.dot( u, N.transpose(v) ) )
223 dist = N.transpose( N.asarray( map( lambda column,a:column+a, \
224 N.transpose(dist), diag1) ) )
225 return N.transpose( N.sqrt( N.asarray(
226 map( lambda row,a: row+a, dist, diag2 ) ) ))
227
228
230 """
231 Create random array of given lenght and number of ones.
232
233 @param nOnes: number of ones
234 @type nOnes: int
235 @param length: lenght of array
236 @type length: int
237
238 @return: array with ones and zeros
239 @rtype: array( 1|0 )
240 """
241 r = N.zeros( length )
242 pos = []
243
244
245 for i in range( nOnes ):
246 pos += [ int( random.random() * length ) ]
247 N.put( r, pos, 1 )
248
249
250 while nOnes != N.sum(r):
251 pos = int( random.random() * length )
252 N.put( r, pos, 1 )
253
254 return r
255
256
258 """
259 Create randomized 2D array containing ones and zeros.
260
261 @param matrix: matrix to randomize
262 @type matrix: 2D array
263 @param mask: mask OR None (default: None)
264 @type mask: list(1|0)
265 @param ranNr: number of matricies to add up (default: 1)
266 @type ranNr: integer
267
268 @return: 2D array or |ranNr| added contact matricies
269 @rtype:2D array
270
271 @raise MathUtilError: if mask does not fit matrix
272 """
273
274 a,b = N.shape( matrix )
275
276
277 if mask is not None:
278 if len(mask) == len( N.ravel(matrix) ):
279 array = N.compress( mask, N.ravel(matrix) )
280
281 if len(mask) != len( N.ravel(matrix) ):
282 raise MathUtilError(
283 'MatUtils.random2DArray - mask of incorrect length' +
284 '\tMatrix length: %i Mask length: %i'\
285 %(len( N.ravel(matrix) ), len(mask)))
286
287 if not mask:
288 array = N.ravel(matrix)
289
290
291 nOnes = int( N.sum( array ) )
292 lenArray = len( array )
293 ranArray = N.zeros( lenArray )
294
295
296 for n in range(ranNr):
297 ranArray += randomMask( nOnes, lenArray )
298
299
300 if mask is not None:
301 r = N.zeros(a*b)
302 N.put( r, N.nonzero(mask), ranArray)
303 return N.reshape( r, (a,b) )
304
305 if not mask:
306 return N.reshape( ranArray, (a,b) )
307
308
310 """
311 Running average (smoothing) over a given data window.
312
313 @param x: data
314 @type x: list of int/float
315 @param interval: window size C{ (-(interval-1)/2 to +(interval-1)/2) }
316 (default: 2)
317 @type interval: int
318 @param preserve_boundaries: shrink window at edges to keep original
319 start and end value (default: 0)
320 @type preserve_boundaries: 0|1
321
322 @return: list of floats
323 @rtype: [ float ]
324 """
325
326 if interval == 0:
327 return x
328
329 l = []
330
331 interval = int((interval-1)/2)
332
333 if not preserve_boundaries:
334
335 for i in range(len(x)):
336
337 left = max(0, i - interval)
338 right = min(len(x), i + interval + 1)
339
340 slice = x[left:right]
341
342 l.append(N.average(slice))
343 else:
344
345 for i in range( len(x) ):
346
347 left = i - interval
348 right= i + interval + 1
349
350 if left < 0:
351 right = right + left
352 left = 0
353 if right > len(x):
354 left = left + right - len(x)
355 right = len(x)
356
357 slice = x[left:right]
358
359 l.append(N.average(slice))
360
361 return N.array(l)
362
363
365 """
366 Compress sparse array of 0 and ones to list of one-positions
367 (space saving function, upack with L{unpackBinaryMatrix}).
368
369 @param cm: X by Y array of int
370 @type cm: 2D array
371
372 @return: {'shape':(X,Y), 'nonzero':[int] }
373 @rtype: dict
374 """
375 if cm == None or type( cm ) == dict:
376 return cm
377
378 result = {}
379 result['shape'] = N.shape( cm )
380 result['nonzero'] = N.nonzero( N.ravel( cm ) )
381 result['nonzero'] = result['nonzero'].tolist()
382 return result
383
384
386 """
387 Uncompress array of 0 and 1 that was compressed with L{packBinaryMatrix}.
388
389 @param pcm: {'shape':(X,Y,..), 'nonzero':[int]}
390 @type pcm: dict
391 @param raveled: return raveled (default: 0)
392 @type raveled: 1|0
393
394 @return: N.array(X by Y by ..) of int
395 @rtype: 2D array
396 """
397 if type( pcm ) != dict:
398 return pcm
399
400 s = pcm['shape']
401
402 m = N.zeros( N.cumproduct( s )[-1], N.Int)
403 pass
404 N.put( m, pcm['nonzero'], 1 )
405
406 if raveled:
407 return m
408
409 return N.reshape( m, s )
410
411
413 """
414 Convert matrix into standard python list remembering the dimensions.
415 Unpack with L{listToMatrix}.
416
417 @note: Not used right now.
418
419 @param cm: array of int
420 @type cm: 2D array
421
422 @return: {'shape':(int,..), 'lst':[..] }
423 @rtype: dict
424 """
425 if cm == None or type( cm ) == dict:
426 return cm
427
428 result = {}
429 result['shape'] = N.shape( cm )
430 result['lst'] = N.ravel( cm ).tolist()
431
432 return result
433
434
436 """
437 Convert result of L{matrixToList} back into Numeric array
438
439 @note: Not used right now.
440
441 @param lcm: {'shape':(int,..), 'lst':[..] }
442 @type lcm: dict
443
444 @return: Numeric.array
445 @rtype:
446 """
447 if type( lcm ) != dict:
448 return lcm
449
450 s = lcm['shape']
451 return N.reshape( lcm['lst'], s )
452
453
455 """
456 Builds a rotation matrix from successive rotations:
457 1. rotation about y-axis by angle alpha
458 2. rotation about x-axis by angle beta
459 3. rotation about z-axis by angle gamma
460
461 @author: Michael Habeck
462
463 @param alpha: euler angle S{alpha}
464 @type alpha: float
465 @param beta: euler angle S{beta}
466 @type beta: float
467 @param gamma: euler angle S{gamma}
468 @type gamma: float
469
470 @return: 3 x 3 array of float
471 @rtype: array
472 """
473 cos_alpha = N.cos(alpha); sin_alpha = N.sin(alpha)
474 cos_beta = N.cos(beta); sin_beta = N.sin(beta)
475 cos_gamma = N.cos(gamma); sin_gamma = N.sin(gamma)
476
477 R = N.zeros((3,3), N.Float32)
478
479 R[0][0] = cos_gamma * cos_alpha - sin_gamma * cos_beta * sin_alpha
480 R[0][1] = cos_gamma * sin_alpha + sin_gamma * cos_beta * cos_alpha
481 R[0][2] = sin_gamma * sin_beta
482
483 R[1][0] = -sin_gamma * cos_alpha - cos_gamma * cos_beta * sin_alpha
484 R[1][1] = -sin_gamma * sin_alpha + cos_gamma * cos_beta * cos_alpha
485 R[1][2] = cos_gamma * sin_beta
486
487 R[2][0] = sin_beta * sin_alpha
488 R[2][1] = -sin_beta * cos_alpha
489 R[2][2] = cos_beta
490
491 return R
492
493
495 """
496 Get random rotation matrix.
497
498 @author: Michael Habeck
499
500 @return: 3 x 3 array of float
501 @rtype: array
502 """
503 alpha = RandomArray.random() * 2 * N.pi
504 gamma = RandomArray.random() * 2 * N.pi
505 beta = N.arccos(2*(RandomArray.random() - 0.5))
506
507 return eulerRotation(alpha, beta, gamma)
508
509
511 """
512 Intersection of the two lists (i.e. all values common to both two lists).
513 C{ intersection( [any], [any] ) -> [any] }
514
515 @param a: first list
516 @type a: [any]
517 @param b: second list
518 @type b: [any]
519 @param optimize: result is sorted like the shorter of the two lists
520 (default: 1)
521 @type optimize: 1|0
522
523 @return: list
524 @rtype: [any]
525 """
526 if optimize and len(a) > len(b):
527 a, b = b, a
528
529 return [ x for x in a if x in b ]
530
531
533 """
534 All members of a list without duplicates
535 C{ noredundant( [any] ) -> [any] }
536
537 @param l: list
538 @type l: [any]
539
540 @return: list
541 @rtype: [any]
542 """
543 r = []
544 for x in l:
545 if x not in r:
546 r.append( x )
547 return r
548
549
551 """
552 Union of two lists (without duplicates)
553 C{ union( [any], [any] ) -> [any] }
554
555 @param a: first list
556 @type a: [any]
557 @param b: second list
558 @type b: [any]
559
560 @return: list
561 @rtype: [any]
562 """
563 if type( a ) is N.arraytype:
564 a = a.tolist()
565 if type( b ) is N.arraytype:
566 b = b.tolist()
567
568 return nonredundant( a + b )
569
570
572 """
573 All members of a that are not in b
574 C{ difference([any], [any]) -> [any] }
575
576 @param a: first list
577 @type a: [any]
578 @param b: second list
579 @type b: [any]
580
581 @return: list
582 @rtype: [any]
583 """
584 return [ x for x in a if x not in b ]
585
586
588 """
589 Remove all or first occurrence(s) of v from l.
590 C{ removeFromList( l, v, [all=1] ) }
591
592 @param l: list
593 @type l: [ any ]
594 @param v: remove these values
595 @type v: any or [ any ]
596 @param all: remove all occurrences (1) or only first one (0) (default: 1)
597 @type all: 0|1
598
599 @return: list
600 @rtype: [any]
601 """
602 if type( v ) != type( [] ):
603 v = [ v ]
604
605 for x in v:
606
607 while x in l:
608 l.remove( x )
609 if not all:
610 break
611
612
614 """
615 Creates a set of n unique integers randomly distributed between
616 start and stop.
617
618 @param start: minimal index
619 @type start: int
620 @param stop: 1+maximal index
621 @type stop: int
622 @param n: number of indices
623 @type n: int
624
625 @return: set of unique integers evenly distributed between start and stop
626 @rtype: [int]
627 """
628 r = []
629 while len( r ) < n:
630 x = int( round( random.uniform( start, stop ) ) )
631 if x < stop and not x in r:
632 r += [x]
633 r.sort()
634 return r
635
636
638 """
639 Calculate linear least-square fit to the points given by x and y.
640 see U{http://mathworld.wolfram.com/LeastSquaresFitting.html}
641
642 @param x: x-data
643 @type x: [ float ]
644 @param y: y-data
645 @type y: [ float ]
646
647 @return: float, float - m, n, r^2 (slope, intersection, corr. coefficient)
648 @rtype: float, float, float
649
650 @raise BiskitError: if x and y have different number of elements
651 """
652 x, y = N.array( x, N.Float64), N.array( y, N.Float64)
653 if len( x ) != len( y ):
654 raise Exception, 'linfit: x and y must have same length'
655
656 av_x = N.average( x )
657 av_y = N.average( y )
658 n = len( x )
659
660 ss_xy = N.sum( x * y ) - n * av_x * av_y
661 ss_xx = N.sum( x * x ) - n * av_x * av_x
662 ss_yy = N.sum( y * y ) - n * av_y * av_y
663
664 slope = ss_xy / ss_xx
665
666 inter = av_y - slope * av_x
667
668 corr = ss_xy**2 / ( ss_xx * ss_yy )
669
670 return slope, inter, corr
671
672
674 """
675 Convert cartesian coordinate array to polar coordinate array:
676 C{ x,y,z -> r, S{theta}, S{phi} }
677
678 @param xyz: array of cartesian coordinates (x, y, z)
679 @type xyz: array
680
681 @return: array of polar coordinates (r, theta, phi)
682 @rtype: array
683 """
684 r = N.sqrt( N.sum( xyz**2, 1 ) )
685 p = N.arccos( xyz[:,2] / r )
686
687
688 t=[]
689 for i in range(len(xyz)):
690
691 t += [math.atan2( xyz[i,1], xyz[i,0] )]
692
693 return N.transpose( N.concatenate( ([r],[t],[p]) ) )
694
695
697 """
698 Convert polar coordinate array to cartesian coordinate array:
699 C{ r, S{theta}, S{phi} -> x,y,z }
700
701 @param rtp: array of cartesian coordinates (r, theta, phi)
702 @type rtp: array
703
704 @return: array of cartesian coordinates (x, y, z)
705 @rtype: array
706 """
707 x = rtp[:,0] * N.cos( rtp[:,1] ) * N.sin( rtp[:,2] )
708 y = rtp[:,0] * N.sin( rtp[:,1] ) * N.sin( rtp[:,2] )
709 z = rtp[:,0] * N.cos( rtp[:,2] )
710
711 return N.transpose( N.concatenate( ([x],[y],[z]) ) )
712
713
715 """
716 Project the coordinates xyz on a sphere with a given radius around
717 a given center.
718
719 @param xyz: cartesian coordinates
720 @type xyz: array N x 3 of float
721 @param radius: radius of target sphere, if not provided the maximal
722 distance to center will be used (default: None)
723 @type radius: float
724 @param center: center of the sphere, if not given the average of xyz
725 will be assigned to the center (default: None)
726 @type center: array 0 x 3 of float
727
728 @return: array of cartesian coordinates (x, y, z)
729 @rtype: array
730 """
731 if center is None:
732 center = N.average( xyz )
733
734 if radius is None:
735 radius = max( N.sqrt( N.sum( N.power( xyz - center, 2 ), 1 ) ) )
736
737 rtp = cartesianToPolar( xyz - center )
738 rtp[ :, 0 ] = radius
739
740 return polarToCartesian( rtp ) + center
741
742
743
744
745
746
747 import Biskit.test as BT
748
749 -class Test(BT.BiskitTest):
750 """Test case"""
751
762
763 EXPECT = N.sum( N.array([ 2.12132034, 0.70710678, 7.07106781]) )
764
765 if __name__ == '__main__':
766
767 BT.localTest()
768