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