source: sasview/src/sas/sascalc/data_util/err1d.py @ 251ef684

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 251ef684 was 9a5097c, checked in by andyfaff, 8 years ago

MAINT: import numpy as np

  • Property mode set to 100644
File size: 3.8 KB
Line 
1# This program is public domain
2"""
3Error propogation algorithms for simple arithmetic
4
5Warning: like the underlying numpy library, the inplace operations
6may return values of the wrong type if some of the arguments are
7integers, so be sure to create them with floating point inputs.
8"""
9from __future__ import division  # Get true division
10import numpy as np
11
12
13def 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
29def 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
45def sub(X, varX, Y, varY):
46    """Subtraction with error propagation"""
47    Z = X - Y
48    varZ = varX + varY
49    return Z, varZ
50
51
52def add(X, varX, Y, varY):
53    """Addition with error propagation"""
54    Z = X + Y
55    varZ = varX + varY
56    return Z, varZ
57
58
59def exp(X, varX):
60    """Exponentiation with error propagation"""
61    Z = np.exp(X)
62    varZ = varX * Z**2
63    return Z, varZ
64
65
66def log(X, varX):
67    """Logarithm with error propagation"""
68    Z = np.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 * np.log(X)**2) * Z**2
76#    return Z,varZ
77#
78
79
80def 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
95def 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
110def 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
125def 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
134def 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
143def 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
Note: See TracBrowser for help on using the repository browser.