| Home | Trees | Indices | Help |
|---|
|
|
1 '''
2 Mathematical operations that generalize many operations from the
3 standard math module so that they also work on numbers with
4 uncertainties.
5
6 Examples:
7
8 from umath import sin
9
10 # Manipulation of numbers with uncertainties:
11 x = uncertainties.ufloat((3, 0.1))
12 print sin(x) # prints 0.141120008...+/-0.098999...
13
14 # The umath functions also work on regular Python floats:
15 print sin(3) # prints 0.141120008... This is a Python float.
16
17 Importing all the functions from this module into the global namespace
18 is possible. This is encouraged when using a Python shell as a
19 calculator. Example:
20
21 import uncertainties
22 from uncertainties.umath import * # Imports tan(), etc.
23
24 x = uncertainties.ufloat((3, 0.1))
25 print tan(x) # tan() is the uncertainties.umath.tan function
26
27 The numbers with uncertainties handled by this module are objects from
28 the uncertainties module, from either the Variable or the
29 AffineScalarFunc class.
30
31 (c) 2009-2010 by Eric O. LEBIGOT (EOL) <eric.lebigot@normalesup.org>.
32 Please send feature requests, bug reports, or feedback to this address.
33
34 This software is released under a dual license. (1) The BSD license.
35 (2) Any other license, as long as it is obtained from the original
36 author.'''
37
38 from __future__ import division # Many analytical derivatives depend on this
39
40 # Standard modules
41 import math
42 import sys
43 import itertools
44 import functools
45 import inspect
46
47 # Local modules
48 from cc.ivs.units import uncertainties
49
50 from cc.ivs.units.uncertainties import __author__, to_affine_scalar, AffineScalarFunc
51
52 ###############################################################################
53
54 # We wrap the functions from the math module so that they keep track of
55 # uncertainties by returning a AffineScalarFunc object.
56
57 # Some functions from the math module cannot be adapted in a standard
58 # way so to work with AffineScalarFunc objects (either as their result
59 # or as their arguments):
60
61 # (1) Some functions return a result of a type whose value and
62 # variations (uncertainties) cannot be represented by AffineScalarFunc
63 # (e.g., math.frexp, which returns a tuple). The exception raised
64 # when not wrapping them with wrap() is more obvious than the
65 # one obtained when wrapping them (in fact, the wrapped functions
66 # attempts operations that are not supported, such as calculation a
67 # subtraction on a result of type tuple).
68
69 # (2) Some functions don't take continuous scalar arguments (which can
70 # be varied during differentiation): math.fsum, math.factorial...
71 # Such functions can either be:
72
73 # - wrapped in a special way.
74
75 # - excluded from standard wrapping by adding their name to
76 # no_std_wrapping
77
78 # Math functions that have a standard interface: they take
79 # one or more float arguments, and return a scalar:
80 many_scalar_to_scalar_funcs = []
81
82 # Some functions require a specific treatment and must therefore be
83 # excluded from standard wrapping. Functions
84 no_std_wrapping = ['modf', 'frexp', 'ldexp', 'fsum', 'factorial']
85
86 # Functions that do not belong in many_scalar_to_scalar_funcs, but
87 # that have a version that handles uncertainties:
88 non_std_wrapped_funcs = []
89
90 # Function that copies the relevant attributes from generalized
91 # functions from the math module:
92 wraps = functools.partial(functools.update_wrapper,
93 assigned=('__doc__', '__name__'))
94
95 ########################################
96 # Wrapping of built-in math functions not in no_std_wrapping:
97
98 # Fixed formulas for the derivatives of some functions from the math
99 # module (some functions might not be present in all version of
100 # Python). Singular points are not taken into account. The user
101 # should never give "large" uncertainties: problems could only appear
102 # if this assumption does not hold.
103
104 # Functions not mentioned in _fixed_derivatives have their derivatives
105 # calculated numerically.
106
107 # Functions that have singularities (possibly at infinity) benefit
108 # from analytical calculations (instead of the default numerical
109 # calculation). Even slowly varying functions (e.g., abs()) yield
110 # more precise results when differentiated analytically, because of
111 # the loss of precision in numerical calculations.
112
113 #def log_1arg_der(x):
114 # """
115 # Derivative of log(x) (1-argument form).
116 # """
117 # return 1/x
118
119 -def log_der0(*args):
120 """
121 Derivative of math.log() with respect to its first argument.
122
123 Works whether 1 or 2 arguments are given.
124 """
125 if len(args) == 1:
126 return 1/args[0]
127 else:
128 return 1/args[0]/math.log(args[1]) # 2-argument form
129
130 # The following version goes about as fast:
131
132 ## A 'try' is used for the most common case because it is fast when no
133 ## exception is raised:
134 #try:
135 # return log_1arg_der(*args) # Argument number check
136 #except TypeError:
137 # return 1/args[0]/math.log(args[1]) # 2-argument form
138
139 fixed_derivatives = {
140 # In alphabetical order, here:
141 'acos': [lambda x: -1/math.sqrt(1-x**2)],
142 'acosh': [lambda x: 1/math.sqrt(x**2-1)],
143 'asin': [lambda x: 1/math.sqrt(1-x**2)],
144 'asinh': [lambda x: 1/math.sqrt(1+x**2)],
145 'atan': [lambda x: 1/(1+x**2)],
146 'atan2': [lambda y, x: x/(x**2+y**2), # Correct for x == 0
147 lambda y, x: -y/(x**2+y**2)], # Correct for x == 0
148 'atanh': [lambda x: 1/(1-x**2)],
149 'ceil': [lambda x: 0],
150 'copysign': [lambda x, y: (1 if x >= 0 else -1) * math.copysign(1, y),
151 lambda x, y: 0],
152 'cos': [lambda x: -math.sin(x)],
153 'cosh': [math.sinh],
154 'degrees': [lambda x: math.degrees(1)],
155 'exp': [math.exp],
156 'fabs': [lambda x: 1 if x >= 0 else -1],
157 'floor': [lambda x: 0],
158 'hypot': [lambda x, y: x/math.hypot(x, y),
159 lambda x, y: y/math.hypot(x, y)],
160 'log': [log_der0,
161 lambda x, y: -math.log(x, y)/y/math.log(y)],
162 'log10': [lambda x: 1/x/math.log(10)],
163 'log1p': [lambda x: 1/(1+x)],
164 'pow': [lambda x, y: y*math.pow(x, y-1),
165 lambda x, y: math.log(x) * math.pow(x, y)],
166 'radians': [lambda x: math.radians(1)],
167 'sin': [math.cos],
168 'sinh': [math.cosh],
169 'sqrt': [lambda x: 0.5/math.sqrt(x)],
170 'tan': [lambda x: 1+math.tan(x)**2],
171 'tanh': [lambda x: 1-math.tanh(x)**2]
172 }
173
174 # Many built-in functions in the math module are wrapped with a
175 # version which is uncertainty aware:
176
177 this_module = sys.modules[__name__]
178
179 # We do not want to wrap module attributes such as __doc__, etc.:
180 for (name, func) in inspect.getmembers(math, inspect.isbuiltin):
181
182 if name in no_std_wrapping:
183 continue
184
185 if name in fixed_derivatives:
186 derivatives = fixed_derivatives[name]
187 else:
188 # Functions whose derivatives are calculated numerically by
189 # this module fall here (isinf, fmod,...):
190 derivatives = None # Means: numerical calculation required
191 setattr(this_module, name,
192 wraps(uncertainties.wrap(func, derivatives), func))
193 many_scalar_to_scalar_funcs.append(name)
194
195 ###############################################################################
196
197 ########################################
198 # Special cases: some of the functions from no_std_wrapping:
199
200 ##########
201 # The math.factorial function is not converted to an uncertainty-aware
202 # function, because it does not handle non-integer arguments: it does
203 # not make sense to give it an argument with a numerical error
204 # (whereas this would be relevant for the gamma function).
205
206 ##########
207
208 # fsum takes a single argument, which cannot be differentiated.
209 # However, each of the arguments inside this single list can
210 # be a variable. We handle this in a specific way:
211
212 if sys.version_info[:2] >= (2, 6):
213
214 # For drop-in compatibility with the math module:
215 factorial = math.factorial
216 non_std_wrapped_funcs.append('factorial')
217
218
219 # We wrap math.fsum
220
221 original_func = math.fsum # For optimization purposes
222
223 # The function below exists so that temporary variables do not
224 # pollute the module namespace:
225 - def wrapped_fsum():
226 """
227 Returns an uncertainty-aware version of math.fsum, which must
228 be contained in _original_func.
229 """
230
231 # The fsum function is flattened, in order to use the
232 # wrap() wrapper:
233
234 flat_fsum = lambda *args: original_func(args)
235
236 flat_fsum_wrap = uncertainties.wrap(
237 flat_fsum, itertools.repeat(lambda *args: 1))
238
239 return wraps(lambda arg_list: flat_fsum_wrap(*arg_list),
240 original_func)
241
242 fsum = wrapped_fsum()
243 non_std_wrapped_funcs.append('fsum')
244
245 ##########
246 @uncertainties.set_doc(math.modf.__doc__)
248 """
249 Version of modf that works for numbers with uncertainty, and also
250 for regular numbers.
251 """
252
253 # The code below is inspired by uncertainties.wrap(). It is
254 # simpler because only 1 argument is given, and there is no
255 # delegation to other functions involved (as for __mul__, etc.).
256
257 aff_func = to_affine_scalar(x)
258
259 (frac_part, int_part) = math.modf(aff_func.nominal_value)
260
261 if aff_func.derivatives:
262 # The derivative of the fractional part is simply 1: the
263 # derivatives of modf(x)[0] are the derivatives of x:
264 return (AffineScalarFunc(frac_part, aff_func.derivatives), int_part)
265 else:
266 # This function was not called with an AffineScalarFunc
267 # argument: there is no need to return numbers with uncertainties:
268 return (frac_part, int_part)
269
270 many_scalar_to_scalar_funcs.append('modf')
271
272 @uncertainties.set_doc(math.ldexp.__doc__)
274 # The code below is inspired by uncertainties.wrap(). It is
275 # simpler because only 1 argument is given, and there is no
276 # delegation to other functions involved (as for __mul__, etc.).
277
278 # Another approach would be to add an additional argument to
279 # uncertainties.wrap() so that some arguments are automatically
280 # considered as constants.
281
282 aff_func = to_affine_scalar(x) # y must be an integer, for math.ldexp
283
284 if aff_func.derivatives:
285 factor = 2**y
286 return AffineScalarFunc(
287 math.ldexp(aff_func.nominal_value, y),
288 # Chain rule:
289 dict((var, factor*deriv)
290 for (var, deriv) in aff_func.derivatives.iteritems()))
291 else:
292 # This function was not called with an AffineScalarFunc
293 # argument: there is no need to return numbers with uncertainties:
294
295 # aff_func.nominal_value is not passed instead of x, because
296 # we do not have to care about the type of the return value of
297 # math.ldexp, this way (aff_func.nominal_value might be the
298 # value of x coerced to a difference type [int->float, for
299 # instance]):
300 return math.ldexp(x, y)
301 many_scalar_to_scalar_funcs.append('ldexp')
302
303 @uncertainties.set_doc(math.frexp.__doc__)
305 """
306 Version of frexp that works for numbers with uncertainty, and also
307 for regular numbers.
308 """
309
310 # The code below is inspired by uncertainties.wrap(). It is
311 # simpler because only 1 argument is given, and there is no
312 # delegation to other functions involved (as for __mul__, etc.).
313
314 aff_func = to_affine_scalar(x)
315
316 if aff_func.derivatives:
317 result = math.frexp(aff_func.nominal_value)
318 # With frexp(x) = (m, e), dm/dx = 1/(2**e):
319 factor = 1/(2**result[1])
320 return (
321 AffineScalarFunc(
322 result[0],
323 # Chain rule:
324 dict((var, factor*deriv)
325 for (var, deriv) in aff_func.derivatives.iteritems())),
326 # The exponent is an integer and is supposed to be
327 # continuous (small errors):
328 result[1])
329 else:
330 # This function was not called with an AffineScalarFunc
331 # argument: there is no need to return numbers with uncertainties:
332 return math.frexp(x)
333 non_std_wrapped_funcs.append('frexp')
334
335 ###############################################################################
336 # Exported functions:
337
338 __all__ = many_scalar_to_scalar_funcs + non_std_wrapped_funcs
339
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Mon Nov 7 18:02:02 2016 | http://epydoc.sourceforge.net |