| 1 | # This program is public domain |
|---|
| 2 | |
|---|
| 3 | """ |
|---|
| 4 | Error propogation algorithms for simple arithmetic |
|---|
| 5 | |
|---|
| 6 | Warning: like the underlying numpy library, the inplace operations |
|---|
| 7 | may return values of the wrong type if some of the arguments are |
|---|
| 8 | integers, so be sure to create them with floating point inputs. |
|---|
| 9 | """ |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | from __future__ import division # Get true division |
|---|
| 13 | import numpy |
|---|
| 14 | |
|---|
| 15 | def div(X,varX, Y,varY): |
|---|
| 16 | """Division with error propagation""" |
|---|
| 17 | # Direct algorithm: |
|---|
| 18 | # Z = X/Y |
|---|
| 19 | # varZ = (varX/X**2 + varY/Y**2) * Z**2 |
|---|
| 20 | # = (varX + varY * Z**2) / Y**2 |
|---|
| 21 | # Indirect algorithm to minimize intermediates |
|---|
| 22 | Z = X/Y # truediv => Z is a float |
|---|
| 23 | varZ = Z**2 # Z is a float => varZ is a float |
|---|
| 24 | varZ *= varY |
|---|
| 25 | varZ += varX |
|---|
| 26 | T = Y**2 # Doesn't matter if T is float or int |
|---|
| 27 | varZ /= T |
|---|
| 28 | return Z, varZ |
|---|
| 29 | |
|---|
| 30 | def mul(X,varX, Y,varY): |
|---|
| 31 | """Multiplication with error propagation""" |
|---|
| 32 | # Direct algorithm: |
|---|
| 33 | Z = X * Y |
|---|
| 34 | varZ = Y**2 * varX + X**2 * varY |
|---|
| 35 | # Indirect algorithm won't ensure floating point results |
|---|
| 36 | # varZ = Y**2 |
|---|
| 37 | # varZ *= varX |
|---|
| 38 | # Z = X**2 # Using Z to hold the temporary |
|---|
| 39 | # Z *= varY |
|---|
| 40 | # varZ += Z |
|---|
| 41 | # Z[:] = X |
|---|
| 42 | # Z *= Y |
|---|
| 43 | return Z, varZ |
|---|
| 44 | |
|---|
| 45 | def sub(X,varX, Y, varY): |
|---|
| 46 | """Subtraction with error propagation""" |
|---|
| 47 | Z = X - Y |
|---|
| 48 | varZ = varX + varY |
|---|
| 49 | return Z, varZ |
|---|
| 50 | |
|---|
| 51 | def add(X,varX, Y,varY): |
|---|
| 52 | """Addition with error propagation""" |
|---|
| 53 | Z = X + Y |
|---|
| 54 | varZ = varX + varY |
|---|
| 55 | return Z, varZ |
|---|
| 56 | |
|---|
| 57 | def exp(X,varX): |
|---|
| 58 | """Exponentiation with error propagation""" |
|---|
| 59 | Z = numpy.exp(X) |
|---|
| 60 | varZ = varX * Z**2 |
|---|
| 61 | return Z,varZ |
|---|
| 62 | |
|---|
| 63 | def log(X,varX): |
|---|
| 64 | """Logarithm with error propagation""" |
|---|
| 65 | Z = numpy.log(X) |
|---|
| 66 | varZ = varX / X**2 |
|---|
| 67 | return Z,varZ |
|---|
| 68 | |
|---|
| 69 | # Confirm this formula before using it |
|---|
| 70 | # def pow(X,varX, Y,varY): |
|---|
| 71 | # Z = X**Y |
|---|
| 72 | # varZ = (Y**2 * varX/X**2 + varY * numpy.log(X)**2) * Z**2 |
|---|
| 73 | # return Z,varZ |
|---|
| 74 | # |
|---|
| 75 | |
|---|
| 76 | def pow(X,varX,n): |
|---|
| 77 | """X**n with error propagation""" |
|---|
| 78 | # Direct algorithm |
|---|
| 79 | # Z = X**n |
|---|
| 80 | # varZ = n*n * varX/X**2 * Z**2 |
|---|
| 81 | # Indirect algorithm to minimize intermediates |
|---|
| 82 | Z = X**n |
|---|
| 83 | varZ = varX/X |
|---|
| 84 | varZ /= X |
|---|
| 85 | varZ *= Z |
|---|
| 86 | varZ *= Z |
|---|
| 87 | varZ *= n**2 |
|---|
| 88 | return Z, varZ |
|---|
| 89 | |
|---|
| 90 | def div_inplace(X,varX, Y,varY): |
|---|
| 91 | """In-place division with error propagation""" |
|---|
| 92 | # Z = X/Y |
|---|
| 93 | # varZ = (varX + varY * (X/Y)**2) / Y**2 = (varX + varY * Z**2) / Y**2 |
|---|
| 94 | X /= Y # X now has Z = X/Y |
|---|
| 95 | T = X**2 # create T with Z**2 |
|---|
| 96 | T *= varY # T now has varY * Z**2 |
|---|
| 97 | varX += T # varX now has varX + varY*Z**2 |
|---|
| 98 | del T # may want to use T[:] = Y for vectors |
|---|
| 99 | T = Y # reuse T for Y |
|---|
| 100 | T **=2 # T now has Y**2 |
|---|
| 101 | varX /= T # varX now has varZ |
|---|
| 102 | return X,varX |
|---|
| 103 | |
|---|
| 104 | def mul_inplace(X,varX, Y,varY): |
|---|
| 105 | """In-place multiplication with error propagation""" |
|---|
| 106 | # Z = X * Y |
|---|
| 107 | # varZ = Y**2 * varX + X**2 * varY |
|---|
| 108 | T = Y**2 # create T with Y**2 |
|---|
| 109 | varX *= T # varX now has Y**2 * varX |
|---|
| 110 | del T # may want to use T[:] = X for vectors |
|---|
| 111 | T = X # reuse T for X**2 * varY |
|---|
| 112 | T **=2 # T now has X**2 |
|---|
| 113 | T *= varY # T now has X**2 * varY |
|---|
| 114 | varX += T # varX now has varZ |
|---|
| 115 | X *= Y # X now has Z |
|---|
| 116 | return X,varX |
|---|
| 117 | |
|---|
| 118 | def sub_inplace(X,varX, Y, varY): |
|---|
| 119 | """In-place subtraction with error propagation""" |
|---|
| 120 | # Z = X - Y |
|---|
| 121 | # varZ = varX + varY |
|---|
| 122 | X -= Y |
|---|
| 123 | varX += varY |
|---|
| 124 | return X,varX |
|---|
| 125 | |
|---|
| 126 | def add_inplace(X,varX, Y,varY): |
|---|
| 127 | """In-place addition with error propagation""" |
|---|
| 128 | # Z = X + Y |
|---|
| 129 | # varZ = varX + varY |
|---|
| 130 | X += Y |
|---|
| 131 | varX += varY |
|---|
| 132 | return X,varX |
|---|
| 133 | |
|---|
| 134 | def pow_inplace(X,varX,n): |
|---|
| 135 | """In-place X**n with error propagation""" |
|---|
| 136 | # Direct algorithm |
|---|
| 137 | # Z = X**n |
|---|
| 138 | # varZ = abs(n) * varX/X**2 * Z**2 |
|---|
| 139 | # Indirect algorithm to minimize intermediates |
|---|
| 140 | varX /= X |
|---|
| 141 | varX /= X # varX now has varX/X**2 |
|---|
| 142 | X **= n # X now has Z = X**n |
|---|
| 143 | varX *= X |
|---|
| 144 | varX *= X # varX now has varX/X**2 * Z**2 |
|---|
| 145 | varX *= n**2 # varX now has varZ |
|---|
| 146 | return X,varX |
|---|
