| 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-2013 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.sigproc.lmfit import uncertainties
49
50 from cc.ivs.sigproc.lmfit.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_scalars_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 with numerical derivatives:
87 num_deriv_funcs = ['fmod', 'gamma', 'isinf', 'isnan',
88 'lgamma', 'trunc']
89
90 # Functions that do not belong in many_scalars_to_scalar_funcs, but
91 # that have a version that handles uncertainties:
92 non_std_wrapped_funcs = []
93
94 # Function that copies the relevant attributes from generalized
95 # functions from the math module:
96 wraps = functools.partial(functools.update_wrapper,
97 assigned=('__doc__', '__name__'))
98
99 ########################################
100 # Wrapping of math functions:
101
102 # Fixed formulas for the derivatives of some functions from the math
103 # module (some functions might not be present in all version of
104 # Python). Singular points are not taken into account. The user
105 # should never give "large" uncertainties: problems could only appear
106 # if this assumption does not hold.
107
108 # Functions not mentioned in _fixed_derivatives have their derivatives
109 # calculated numerically.
110
111 # Functions that have singularities (possibly at infinity) benefit
112 # from analytical calculations (instead of the default numerical
113 # calculation) because their derivatives generally change very fast.
114 # Even slowly varying functions (e.g., abs()) yield more precise
115 # results when differentiated analytically, because of the loss of
116 # precision in numerical calculations.
117
118 #def log_1arg_der(x):
119 # """
120 # Derivative of log(x) (1-argument form).
121 # """
122 # return 1/x
123
124 -def log_der0(*args):
125 """
126 Derivative of math.log() with respect to its first argument.
127
128 Works whether 1 or 2 arguments are given.
129 """
130 if len(args) == 1:
131 return 1/args[0]
132 else:
133 return 1/args[0]/math.log(args[1]) # 2-argument form
134
135 # The following version goes about as fast:
136
137 ## A 'try' is used for the most common case because it is fast when no
138 ## exception is raised:
139 #try:
140 # return log_1arg_der(*args) # Argument number check
141 #except TypeError:
142 # return 1/args[0]/math.log(args[1]) # 2-argument form
143
144 _erf_coef = 2/math.sqrt(math.pi) # Optimization for erf()
145
146 fixed_derivatives = {
147 # In alphabetical order, here:
148 'acos': [lambda x: -1/math.sqrt(1-x**2)],
149 'acosh': [lambda x: 1/math.sqrt(x**2-1)],
150 'asin': [lambda x: 1/math.sqrt(1-x**2)],
151 'asinh': [lambda x: 1/math.sqrt(1+x**2)],
152 'atan': [lambda x: 1/(1+x**2)],
153 'atan2': [lambda y, x: x/(x**2+y**2), # Correct for x == 0
154 lambda y, x: -y/(x**2+y**2)], # Correct for x == 0
155 'atanh': [lambda x: 1/(1-x**2)],
156 'ceil': [lambda x: 0],
157 'copysign': [lambda x, y: (1 if x >= 0 else -1) * math.copysign(1, y),
158 lambda x, y: 0],
159 'cos': [lambda x: -math.sin(x)],
160 'cosh': [math.sinh],
161 'degrees': [lambda x: math.degrees(1)],
162 'erf': [lambda x: exp(-x**2)*_erf_coef],
163 'erfc': [lambda x: -exp(-x**2)*_erf_coef],
164 'exp': [math.exp],
165 'expm1': [math.exp],
166 'fabs': [lambda x: 1 if x >= 0 else -1],
167 'floor': [lambda x: 0],
168 'hypot': [lambda x, y: x/math.hypot(x, y),
169 lambda x, y: y/math.hypot(x, y)],
170 'log': [log_der0,
171 lambda x, y: -math.log(x, y)/y/math.log(y)],
172 'log10': [lambda x: 1/x/math.log(10)],
173 'log1p': [lambda x: 1/(1+x)],
174 'pow': [lambda x, y: y*math.pow(x, y-1),
175 lambda x, y: math.log(x) * math.pow(x, y)],
176 'radians': [lambda x: math.radians(1)],
177 'sin': [math.cos],
178 'sinh': [math.cosh],
179 'sqrt': [lambda x: 0.5/math.sqrt(x)],
180 'tan': [lambda x: 1+math.tan(x)**2],
181 'tanh': [lambda x: 1-math.tanh(x)**2]
182 }
183
184 # Many built-in functions in the math module are wrapped with a
185 # version which is uncertainty aware:
186
187 this_module = sys.modules[__name__]
188
189 # for (name, attr) in vars(math).items():
190 for name in dir(math):
191
192 if name in fixed_derivatives: # Priority to functions in fixed_derivatives
193 derivatives = fixed_derivatives[name]
194 elif name in num_deriv_funcs:
195 # Functions whose derivatives are calculated numerically by
196 # this module fall here (isinf, fmod,...):
197 derivatives = None # Means: numerical calculation required
198 else:
199 continue # 'name' not wrapped by this module (__doc__, e, etc.)
200
201 func = getattr(math, name)
202
203 setattr(this_module, name,
204 wraps(uncertainties.wrap(func, derivatives), func))
205
206 many_scalars_to_scalar_funcs.append(name)
207
208 ###############################################################################
209
210 ########################################
211 # Special cases: some of the functions from no_std_wrapping:
212
213 ##########
214 # The math.factorial function is not converted to an uncertainty-aware
215 # function, because it does not handle non-integer arguments: it does
216 # not make sense to give it an argument with a numerical error
217 # (whereas this would be relevant for the gamma function).
218
219 ##########
220
221 # fsum takes a single argument, which cannot be differentiated.
222 # However, each of the arguments inside this single list can
223 # be a variable. We handle this in a specific way:
224
225 if sys.version_info[:2] >= (2, 6):
226
227 # For drop-in compatibility with the math module:
228 factorial = math.factorial
229 non_std_wrapped_funcs.append('factorial')
230
231
232 # We wrap math.fsum
233
234 original_func = math.fsum # For optimization purposes
235
236 # The function below exists so that temporary variables do not
237 # pollute the module namespace:
238 - def wrapped_fsum():
239 """
240 Returns an uncertainty-aware version of math.fsum, which must
241 be contained in _original_func.
242 """
243
244 # The fsum function is flattened, in order to use the
245 # wrap() wrapper:
246
247 flat_fsum = lambda *args: original_func(args)
248
249 flat_fsum_wrap = uncertainties.wrap(
250 flat_fsum, itertools.repeat(lambda *args: 1))
251
252 return wraps(lambda arg_list: flat_fsum_wrap(*arg_list),
253 original_func)
254
255 fsum = wrapped_fsum()
256 non_std_wrapped_funcs.append('fsum')
257
258 ##########
259 @uncertainties.set_doc(math.modf.__doc__)
261 """
262 Version of modf that works for numbers with uncertainty, and also
263 for regular numbers.
264 """
265
266 # The code below is inspired by uncertainties.wrap(). It is
267 # simpler because only 1 argument is given, and there is no
268 # delegation to other functions involved (as for __mul__, etc.).
269
270 aff_func = to_affine_scalar(x)
271
272 (frac_part, int_part) = math.modf(aff_func.nominal_value)
273
274 if aff_func.derivatives:
275 # The derivative of the fractional part is simply 1: the
276 # derivatives of modf(x)[0] are the derivatives of x:
277 return (AffineScalarFunc(frac_part, aff_func.derivatives), int_part)
278 else:
279 # This function was not called with an AffineScalarFunc
280 # argument: there is no need to return numbers with uncertainties:
281 return (frac_part, int_part)
282
283 many_scalars_to_scalar_funcs.append('modf')
284
285 @uncertainties.set_doc(math.ldexp.__doc__)
287 # The code below is inspired by uncertainties.wrap(). It is
288 # simpler because only 1 argument is given, and there is no
289 # delegation to other functions involved (as for __mul__, etc.).
290
291 # Another approach would be to add an additional argument to
292 # uncertainties.wrap() so that some arguments are automatically
293 # considered as constants.
294
295 aff_func = to_affine_scalar(x) # y must be an integer, for math.ldexp
296
297 if aff_func.derivatives:
298 factor = 2**y
299 return AffineScalarFunc(
300 math.ldexp(aff_func.nominal_value, y),
301 # Chain rule:
302 dict((var, factor*deriv)
303 for (var, deriv) in aff_func.derivatives.iteritems()))
304 else:
305 # This function was not called with an AffineScalarFunc
306 # argument: there is no need to return numbers with uncertainties:
307
308 # aff_func.nominal_value is not passed instead of x, because
309 # we do not have to care about the type of the return value of
310 # math.ldexp, this way (aff_func.nominal_value might be the
311 # value of x coerced to a difference type [int->float, for
312 # instance]):
313 return math.ldexp(x, y)
314 many_scalars_to_scalar_funcs.append('ldexp')
315
316 @uncertainties.set_doc(math.frexp.__doc__)
318 """
319 Version of frexp that works for numbers with uncertainty, and also
320 for regular numbers.
321 """
322
323 # The code below is inspired by uncertainties.wrap(). It is
324 # simpler because only 1 argument is given, and there is no
325 # delegation to other functions involved (as for __mul__, etc.).
326
327 aff_func = to_affine_scalar(x)
328
329 if aff_func.derivatives:
330 result = math.frexp(aff_func.nominal_value)
331 # With frexp(x) = (m, e), dm/dx = 1/(2**e):
332 factor = 1/(2**result[1])
333 return (
334 AffineScalarFunc(
335 result[0],
336 # Chain rule:
337 dict((var, factor*deriv)
338 for (var, deriv) in aff_func.derivatives.iteritems())),
339 # The exponent is an integer and is supposed to be
340 # continuous (small errors):
341 result[1])
342 else:
343 # This function was not called with an AffineScalarFunc
344 # argument: there is no need to return numbers with uncertainties:
345 return math.frexp(x)
346 non_std_wrapped_funcs.append('frexp')
347
348 ###############################################################################
349 # Exported functions:
350
351 __all__ = many_scalars_to_scalar_funcs + non_std_wrapped_funcs
352
| Home | Trees | Indices | Help |
|---|
| Generated by Epydoc 3.0.1 on Mon Nov 7 18:02:03 2016 | http://epydoc.sourceforge.net |