source: sasview/sansutil/err1d.py @ fd080830

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since fd080830 was 4025019, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Adding util and park

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