| 1 | #!/usr/bin/env python |
|---|
| 2 | """Auxiliary numerical tools |
|---|
| 3 | |
|---|
| 4 | """ |
|---|
| 5 | |
|---|
| 6 | import Numeric |
|---|
| 7 | |
|---|
| 8 | def mean(x): |
|---|
| 9 | """Mean value of a vector |
|---|
| 10 | """ |
|---|
| 11 | return(float(Numeric.sum(x))/len(x)) |
|---|
| 12 | |
|---|
| 13 | |
|---|
| 14 | def cov(x, y=None): |
|---|
| 15 | """Covariance of vectors x and y. |
|---|
| 16 | |
|---|
| 17 | If y is None: return cov(x, x) |
|---|
| 18 | """ |
|---|
| 19 | |
|---|
| 20 | if y is None: |
|---|
| 21 | y = x |
|---|
| 22 | |
|---|
| 23 | assert(len(x)==len(y)) |
|---|
| 24 | N = len(x) |
|---|
| 25 | |
|---|
| 26 | cx = x - mean(x) |
|---|
| 27 | cy = y - mean(y) |
|---|
| 28 | |
|---|
| 29 | p = Numeric.innerproduct(cx,cy)/N |
|---|
| 30 | return(p) |
|---|
| 31 | |
|---|
| 32 | |
|---|
| 33 | def err(x, y=0, n=2): |
|---|
| 34 | """Relative error of ||x-y|| to ||y|| |
|---|
| 35 | n = 2: Two norm |
|---|
| 36 | n = None: Max norm |
|---|
| 37 | |
|---|
| 38 | If denominator evaluates to zero, absolute error is returned |
|---|
| 39 | """ |
|---|
| 40 | |
|---|
| 41 | if n == 2: |
|---|
| 42 | err = norm(x-y) |
|---|
| 43 | try: |
|---|
| 44 | err = err/norm(y) |
|---|
| 45 | except: |
|---|
| 46 | pass |
|---|
| 47 | |
|---|
| 48 | else: |
|---|
| 49 | err = max(abs(x-y)) |
|---|
| 50 | try: |
|---|
| 51 | err = err/max(abs(y)) |
|---|
| 52 | except: |
|---|
| 53 | pass |
|---|
| 54 | |
|---|
| 55 | return err |
|---|
| 56 | |
|---|
| 57 | |
|---|
| 58 | def norm(x): |
|---|
| 59 | """2-norm of x |
|---|
| 60 | """ |
|---|
| 61 | |
|---|
| 62 | y = Numeric.ravel(x) |
|---|
| 63 | p = Numeric.sqrt(Numeric.innerproduct(y,y)) |
|---|
| 64 | return p |
|---|
| 65 | |
|---|
| 66 | |
|---|
| 67 | def corr(x, y=None): |
|---|
| 68 | """Correlation of x and y |
|---|
| 69 | If y is None return autocorrelation of x |
|---|
| 70 | """ |
|---|
| 71 | |
|---|
| 72 | from math import sqrt |
|---|
| 73 | if y is None: |
|---|
| 74 | y = x |
|---|
| 75 | |
|---|
| 76 | varx = cov(x) |
|---|
| 77 | vary = cov(y) |
|---|
| 78 | |
|---|
| 79 | if varx == 0 or vary == 0: |
|---|
| 80 | C = 0 |
|---|
| 81 | else: |
|---|
| 82 | C = cov(x,y)/sqrt(varx * vary) |
|---|
| 83 | |
|---|
| 84 | return(C) |
|---|
| 85 | |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | def ensure_numeric(A, typecode = None): |
|---|
| 89 | """Ensure that sequence is a Numeric array. |
|---|
| 90 | Inputs: |
|---|
| 91 | A: Sequence. If A is already a Numeric array it will be returned |
|---|
| 92 | unaltered |
|---|
| 93 | If not, an attempt is made to convert it to a Numeric |
|---|
| 94 | array |
|---|
| 95 | typecode: Numeric type. If specified, use this in the conversion. |
|---|
| 96 | If not, let Numeric decide |
|---|
| 97 | |
|---|
| 98 | This function is necessary as array(A) can cause memory overflow. |
|---|
| 99 | """ |
|---|
| 100 | |
|---|
| 101 | from Numeric import ArrayType, array |
|---|
| 102 | |
|---|
| 103 | if typecode is None: |
|---|
| 104 | if type(A) == ArrayType: |
|---|
| 105 | return A |
|---|
| 106 | else: |
|---|
| 107 | return array(A) |
|---|
| 108 | else: |
|---|
| 109 | if type(A) == ArrayType: |
|---|
| 110 | if A.typecode == typecode: |
|---|
| 111 | return array(A) #FIXME: Shouldn't this just return A? |
|---|
| 112 | else: |
|---|
| 113 | return A.astype(typecode) |
|---|
| 114 | else: |
|---|
| 115 | return array(A).astype(typecode) |
|---|
| 116 | |
|---|
| 117 | |
|---|
| 118 | |
|---|
