ComboCode.cc.ivs.units.uncertainties

1 #!/usr/bin/env python 2 3 #!! When the documentation below is updated, setup.py should be 4 # checked for consistency. 5 6 ''' 7 Calculations with full error propagation for quantities with uncertainties. 8 Derivatives can also be calculated. 9 10 Web user guide: http://packages.python.org/uncertainties/. 11 12 Example of possible calculation: (0.2 +/- 0.01)**2 = 0.04 +/- 0.004. 13 14 Correlations between expressions are correctly taken into account (for 15 instance, with x = 0.2+/-0.01, 2*x-x-x is exactly zero, as is y-x-x 16 with y = 2*x). 17 18 Examples: 19 20 import uncertainties 21 from uncertainties import ufloat 22 from uncertainties.umath import * # sin(), etc. 23 24 # Mathematical operations: 25 x = ufloat((0.20, 0.01)) # x = 0.20+/-0.01 26 x = ufloat("0.20+/-0.01") # Other representation 27 x = ufloat("0.20(1)") # Other representation 28 x = ufloat("0.20") # Implicit uncertainty of +/-1 on the last digit 29 print x**2 # Square: prints "0.04+/-0.004" 30 print sin(x**2) # Prints "0.0399...+/-0.00399..." 31 32 print x.position_in_sigmas(0.17) # Prints "-3.0": deviation of -3 sigmas 33 34 # Access to the nominal value, and to the uncertainty: 35 square = x**2 # Square 36 print square # Prints "0.04+/-0.004" 37 print square.nominal_value # Prints "0.04" 38 print square.std_dev() # Prints "0.004..." 39 40 print square.derivatives[x] # Partial derivative: 0.4 (= 2*0.20) 41 42 # Correlations: 43 u = ufloat((1, 0.05), "u variable") # Tag 44 v = ufloat((10, 0.1), "v variable") 45 sum_value = u+v 46 47 u.set_std_dev(0.1) # Standard deviations can be updated on the fly 48 print sum_value - u - v # Prints "0.0" (exact result) 49 50 # List of all sources of error: 51 print sum_value # Prints "11+/-0.1414..." 52 for (var, error) in sum_value.error_components().iteritems(): 53 print "%s: %f" % (var.tag, error) # Individual error components 54 55 # Covariance matrices: 56 cov_matrix = uncertainties.covariance_matrix([u, v, sum_value]) 57 print cov_matrix # 3x3 matrix 58 59 # Correlated variables can be constructed from a covariance matrix, if 60 # NumPy is available: 61 (u2, v2, sum2) = uncertainties.correlated_values([1, 10, 11], 62 cov_matrix) 63 print u2 # Value and uncertainty of u: correctly recovered (1+/-0.1) 64 print uncertainties.covariance_matrix([u2, v2, sum2]) # == cov_matrix 65 66 - The main function provided by this module is ufloat, which creates 67 numbers with uncertainties (Variable objects). Variable objects can 68 be used as if they were regular Python numbers. The main attributes 69 and methods of Variable objects are defined in the documentation of 70 the Variable class. 71 72 - Valid operations on numbers with uncertainties include basic 73 mathematical functions (addition, etc.). 74 75 Most operations from the standard math module (sin, etc.) can be applied 76 on numbers with uncertainties by using their generalization from the 77 uncertainties.umath module: 78 79 from uncertainties.umath import sin 80 print sin(ufloat("1+/-0.01")) # 0.841...+/-0.005... 81 print sin(1) # umath.sin() also works on floats, exactly like math.sin() 82 83 Logical operations (>, ==, etc.) are also supported. 84 85 Basic operations on NumPy arrays or matrices of numbers with 86 uncertainties can be performed: 87 88 2*numpy.array([ufloat((1, 0.01)), ufloat((2, 0.1))]) 89 90 More complex operations on NumPy arrays can be performed through the 91 dedicated uncertainties.unumpy sub-module (see its documentation). 92 93 Calculations that are performed through non-Python code (Fortran, C, 94 etc.) can handle numbers with uncertainties instead of floats through 95 the provided wrap() wrapper: 96 97 import uncertainties 98 99 # wrapped_f is a version of f that can take arguments with 100 # uncertainties, even if f only takes floats: 101 wrapped_f = uncertainties.wrap(f) 102 103 If some derivatives of the wrapped function f are known (analytically, 104 or numerically), they can be given to wrap()--see the documentation 105 for wrap(). 106 107 - Utility functions are also provided: the covariance matrix between 108 random variables can be calculated with covariance_matrix(), or used 109 as input for the definition of correlated quantities (correlated_values() 110 function--defined only if the NumPy module is available). 111 112 - Mathematical expressions involving numbers with uncertainties 113 generally return AffineScalarFunc objects, which also print as a value 114 with uncertainty. Their most useful attributes and methods are 115 described in the documentation for AffineScalarFunc. Note that 116 Variable objects are also AffineScalarFunc objects. UFloat is an 117 alias for AffineScalarFunc, provided as a convenience: testing whether 118 a value carries an uncertainty handled by this module should be done 119 with insinstance(my_value, UFloat). 120 121 - Mathematically, numbers with uncertainties are, in this package, 122 probability distributions. These probabilities are reduced to two 123 numbers: a nominal value and an uncertainty. Thus, both variables 124 (Variable objects) and the result of mathematical operations 125 (AffineScalarFunc objects) contain these two values (respectively in 126 their nominal_value attribute and through their std_dev() method). 127 128 The uncertainty of a number with uncertainty is simply defined in 129 this package as the standard deviation of the underlying probability 130 distribution. 131 132 The numbers with uncertainties manipulated by this package are assumed 133 to have a probability distribution mostly contained around their 134 nominal value, in an interval of about the size of their standard 135 deviation. This should cover most practical cases. A good choice of 136 nominal value for a number with uncertainty is thus the median of its 137 probability distribution, the location of highest probability, or the 138 average value. 139 140 - When manipulating ensembles of numbers, some of which contain 141 uncertainties, it can be useful to access the nominal value and 142 uncertainty of all numbers in a uniform manner: 143 144 x = ufloat("3+/-0.1") 145 print nominal_value(x) # Prints 3 146 print std_dev(x) # Prints 0.1 147 print nominal_value(3) # Prints 3: nominal_value works on floats 148 print std_dev(3) # Prints 0: std_dev works on floats 149 150 - Probability distributions (random variables and calculation results) 151 are printed as: 152 153 nominal value +/- standard deviation 154 155 but this does not imply any property on the nominal value (beyond the 156 fact that the nominal value is normally inside the region of high 157 probability density), or that the probability distribution of the 158 result is symmetrical (this is rarely strictly the case). 159 160 - Linear approximations of functions (around the nominal values) are 161 used for the calculation of the standard deviation of mathematical 162 expressions with this package. 163 164 The calculated standard deviations and nominal values are thus 165 meaningful approximations as long as the functions involved have 166 precise linear expansions in the region where the probability 167 distribution of their variables is the largest. It is therefore 168 important that uncertainties be small. Mathematically, this means 169 that the linear term of functions around the nominal values of their 170 variables should be much larger than the remaining higher-order terms 171 over the region of significant probability. 172 173 For instance, sin(0+/-0.01) yields a meaningful standard deviation 174 since it is quite linear over 0+/-0.01. However, cos(0+/-0.01) yields 175 an approximate standard deviation of 0 (because the cosine is not well 176 approximated by a line around 0), which might not be precise enough 177 for all applications. 178 179 - Comparison operations (>, ==, etc.) on numbers with uncertainties 180 have a pragmatic semantics, in this package: numbers with 181 uncertainties can be used wherever Python numbers are used, most of 182 the time with a result identical to the one that would be obtained 183 with their nominal value only. However, since the objects defined in 184 this module represent probability distributions and not pure numbers, 185 comparison operator are interpreted in a specific way. 186 187 The result of a comparison operation ("==", ">", etc.) is defined so as 188 to be essentially consistent with the requirement that uncertainties 189 be small: the value of a comparison operation is True only if the 190 operation yields True for all infinitesimal variations of its random 191 variables, except, possibly, for an infinitely small number of cases. 192 193 Example: 194 195 "x = 3.14; y = 3.14" is such that x == y 196 197 but 198 199 x = ufloat((3.14, 0.01)) 200 y = ufloat((3.14, 0.01)) 201 202 is not such that x == y, since x and y are independent random 203 variables that almost never give the same value. However, x == x 204 still holds. 205 206 The boolean value (bool(x), "if x...") of a number with uncertainty x 207 is the result of x != 0. 208 209 - The uncertainties package is for Python 2.5 and above. 210 211 - This package contains tests. They can be run either manually or 212 automatically with the nose unit testing framework (nosetests). 213 214 (c) 2009-2010 by Eric O. LEBIGOT (EOL) <eric.lebigot@normalesup.org>. 215 Please send feature requests, bug reports, or feedback to this address. 216 217 Please support future development by donating $5 or more through PayPal! 218 219 This software is released under a dual license. (1) The BSD license. 220 (2) Any other license, as long as it is obtained from the original 221 author.''' 222 223 # The idea behind this module is to replace the result of mathematical 224 # operations by a local approximation of the defining function. For 225 # example, sin(0.2+/-0.01) becomes the affine function 226 # (AffineScalarFunc object) whose nominal value is sin(0.2) and 227 # whose variations are given by sin(0.2+delta) = 0.98...*delta. 228 # Uncertainties can then be calculated by using this local linear 229 # approximation of the original function. 230 231 from __future__ import division # Many analytical derivatives depend on this 232 233 import re 234 import math 235 from math import sqrt, log # Optimization: no attribute look-up 236 import copy 237 238 # Numerical version: 239 __version_info__ = (1, 7, 1) 240 __version__ = '.'.join(str(num) for num in __version_info__) 241 242 __author__ = 'Eric O. LEBIGOT (EOL)' 243 244 # Attributes that are always exported (some other attributes are 245 # exported only if the NumPy module is available...): 246 __all__ = [ 247 248 # All sub-modules and packages are not imported by default, 249 # because NumPy could not be installed: 250 251 'ufloat', # Main function: returns a number with uncertainty 252 253 # Uniform access to nominal values and standard deviations: 254 'nominal_value', 255 'std_dev', 256 257 # Utility functions (more are exported if NumPy is present): 258 'covariance_matrix', 259 260 # Class for testing whether an object is a number with 261 # uncertainty. Not usually created by users (except through the 262 # Variable subclass), but possibly manipulated by external code 263 # ['derivatives()' method, etc.]. 264 'UFloat', 265 266 # Wrapper for allowing non-pure-Python function to handle 267 # quantitites with uncertainties: 268 'wrap', 269 270 # The documentation for wrap() indicates that numerical 271 # derivatives are calculated through partial_derivative(). The 272 # user might also want to change the size of the numerical 273 # differentiation step. 274 'partial_derivative' 275 276 ]

277 278 ############################################################################### 279 280 -def set_doc(doc_string):

281 """ 282 Decorator function that sets the docstring to the given text. 283 284 It is useful for functions whose docstring is calculated 285 (including string substitutions). 286 """ 287 def set_doc_string(func): 288 func.__doc__ = doc_string 289 return func

290 return set_doc_string 291

292 ############################################################################### 293 294 # Mathematical operations with local approximations (affine scalar 295 # functions) 296 297 -class NotUpcast(Exception):

298 pass

299

300 -def to_affine_scalar(x):

301 """ 302 Transforms x into a constant affine scalar function 303 (AffineScalarFunc), unless it is already an AffineScalarFunc (in 304 which case x is returned unchanged). 305 306 Raises an exception unless 'x' belongs to some specific classes of 307 objects that are known not to depend on AffineScalarFunc objects 308 (which then cannot be considered as constants). 309 """ 310 311 if isinstance(x, AffineScalarFunc): 312 return x 313 314 #! In Python 2.6+, numbers.Number could be used instead, here: 315 if isinstance(x, (float, int, complex, long)): 316 # No variable => no derivative to define: 317 return AffineScalarFunc(x, {}) 318 319 # Case of lists, etc. 320 raise NotUpcast("%s cannot be converted to a number with" 321 " uncertainty" % type(x))

322

323 -def partial_derivative(f, param_num):

324 """ 325 Returns a function that numerically calculates the partial 326 derivative of function f with respect to its argument number 327 param_num. 328 329 The step parameter represents the shift of the parameter used in 330 the numerical approximation. 331 """ 332 333 def partial_derivative_of_f(*args): 334 """ 335 Partial derivative, calculated with the (-epsilon, +epsilon) 336 method, which is more precise than the (0, +epsilon) method. 337 """ 338 # f_nominal_value = f(*args) 339 340 shifted_args = list(args) # Copy, and conversion to a mutable 341 342 # The step is relative to the parameter being varied, so that 343 # shifting it does not suffer from finite precision: 344 step = 1e-8*abs(shifted_args[param_num]) 345 if not step: 346 step = 1e-8 # Arbitrary, but "small" with respect to 1 347 348 shifted_args[param_num] += step 349 shifted_f_plus = f(*shifted_args) 350 351 shifted_args[param_num] -= 2*step # Optimization: only 1 list copy 352 shifted_f_minus = f(*shifted_args) 353 354 return (shifted_f_plus - shifted_f_minus)/2/step

355 356 return partial_derivative_of_f 357

358 -class NumericalDerivatives(object):

359 """ 360 Sequence with the successive numerical derivatives of a function. 361 """ 362 # This is a sequence and not a list because the number of 363 # arguments of the function is not known in advance, in general. 364

365 - def __init__(self, function):

366 """ 367 'function' is the function whose derivatives can be computed. 368 """ 369 self._function = function

370

371 - def __getitem__(self, n):

372 """ 373 Returns the n-th numerical derivative of the function. 374 """ 375 return partial_derivative(self._function, n)

376

377 -def wrap(f, derivatives_funcs=None):

378 """ 379 Wraps function f so that, when applied to numbers with 380 uncertainties (AffineScalarFunc objects) or float-like arguments, 381 f returns a local approximation of its values (in the form of an 382 object of class AffineScalarFunc). In this case, if none of the 383 arguments of f involves variables [i.e. Variable objects], f 384 simply returns its usual result. 385 386 When f is not called on AffineScalarFunc or float-like 387 arguments, the original result of f is returned. 388 389 If supplied, 'derivatives_funcs' is a sequence or iterator that 390 generally contains functions; each successive function is the 391 partial derivatives of f with respect to the corresponding 392 variable (one function for each argument of f, which takes as many 393 arguments as f). If derivatives_funcs is None, or if 394 derivatives_funcs contains a finite number of elements, then 395 missing derivatives are calculated numerically through 396 partial_derivative(). 397 398 Example: wrap(math.sin) is a sine function that can be applied to 399 numbers with uncertainties. Its derivative will be calculated 400 numerically. wrap(math.sin, [None]) would have produced the same 401 result. wrap(math.sin, [math.cos]) is the same function, but with 402 an analytically defined derivative. 403 """ 404 405 if derivatives_funcs is None: 406 derivatives_funcs = NumericalDerivatives(f) 407 else: 408 # Derivatives that are not defined are calculated numerically, 409 # if there is a finite number of them (the function lambda 410 # *args: fsum(args) has a non-defined number of arguments, as 411 # it just performs a sum... 412 try: # Is the number of derivatives fixed? 413 len(derivatives_funcs) 414 except TypeError: 415 pass 416 else: 417 derivatives_funcs = [ 418 partial_derivative(f, k) if derivative is None 419 else derivative 420 for (k, derivative) in enumerate(derivatives_funcs)] 421 422 #! Setting the doc string after "def f_with...()" does not 423 # seem to work. We define it explicitly: 424 @set_doc("""\ 425 Version of %s(...) that returns an affine approximation 426 (AffineScalarFunc object), if its result depends on variables 427 (Variable objects). Otherwise, returns a simple constant (when 428 applied to constant arguments). 429 430 Warning: arguments of the function that are not AffineScalarFunc 431 objects must not depend on uncertainties.Variable objects in any 432 way. Otherwise, the dependence of the result in 433 uncertainties.Variable objects will be incorrect. 434 435 Original documentation: 436 %s""" % (f.__name__, f.__doc__)) 437 def f_with_affine_output(*args): 438 # Can this function perform the calculation of an 439 # AffineScalarFunc (or maybe float) result? 440 try: 441 aff_funcs = map(to_affine_scalar, args) 442 443 except NotUpcast: 444 445 # This function does not now how to itself perform 446 # calculations with non-float-like arguments (as they 447 # might for instance be objects whose value really changes 448 # if some Variable objects had different values): 449 450 # Is it clear that we can't delegate the calculation? 451 452 if any(isinstance(arg, AffineScalarFunc) for arg in args): 453 # This situation arises for instance when calculating 454 # AffineScalarFunc(...)*numpy.array(...). In this 455 # case, we must let NumPy handle the multiplication 456 # (which is then performed element by element): 457 return NotImplemented 458 else: 459 # If none of the arguments is an AffineScalarFunc, we 460 # can delegate the calculation to the original 461 # function. This can be useful when it is called with 462 # only one argument (as in 463 # numpy.log10(numpy.ndarray(...)): 464 return f(*args) 465 466 ######################################## 467 # Nominal value of the constructed AffineScalarFunc: 468 args_values = [e.nominal_value for e in aff_funcs] 469 f_nominal_value = f(*args_values) 470 471 ######################################## 472 473 # List of involved variables (Variable objects): 474 variables = set() 475 for expr in aff_funcs: 476 variables |= set(expr.derivatives) 477 478 ## It is sometimes useful to only return a regular constant: 479 480 # (1) Optimization / convenience behavior: when 'f' is called 481 # on purely constant values (e.g., sin(2)), there is no need 482 # for returning a more complex AffineScalarFunc object. 483 484 # (2) Functions that do not return a "float-like" value might 485 # not have a relevant representation as an AffineScalarFunc. 486 # This includes boolean functions, since their derivatives are 487 # either 0 or are undefined: they are better represented as 488 # Python constants than as constant AffineScalarFunc functions. 489 490 if not variables or isinstance(f_nominal_value, bool): 491 return f_nominal_value 492 493 # The result of 'f' does depend on 'variables'... 494 495 ######################################## 496 497 # Calculation of the derivatives with respect to the arguments 498 # of f (aff_funcs): 499 500 # The chain rule is applied. This is because, in the case of 501 # numerical derivatives, it allows for a better-controlled 502 # numerical stability than numerically calculating the partial 503 # derivatives through '[f(x + dx, y + dy, ...) - 504 # f(x,y,...)]/da' where dx, dy,... are calculated by varying 505 # 'a'. In fact, it is numerically better to control how big 506 # (dx, dy,...) are: 'f' is a simple mathematical function and 507 # it is possible to know how precise the df/dx are (which is 508 # not possible with the numerical df/da calculation above). 509 510 # We use numerical derivatives, if we don't already have a 511 # list of derivatives: 512 513 #! Note that this test could be avoided by requiring the 514 # caller to always provide derivatives. When changing the 515 # functions of the math module, this would force this module 516 # to know about all the math functions. Another possibility 517 # would be to force derivatives_funcs to contain, say, the 518 # first 3 derivatives of f. But any of these two ideas has a 519 # chance to break, one day... (if new functions are added to 520 # the math module, or if some function has more than 3 521 # arguments). 522 523 derivatives_wrt_args = [] 524 for (arg, derivative) in zip(aff_funcs, derivatives_funcs): 525 derivatives_wrt_args.append(derivative(*args_values) 526 if arg.derivatives 527 else 0) 528 529 ######################################## 530 # Calculation of the derivative of f with respect to all the 531 # variables (Variable) involved. 532 533 # Initial value (is updated below): 534 derivatives_wrt_vars = dict((var, 0.) for var in variables) 535 536 # The chain rule is used (we already have 537 # derivatives_wrt_args): 538 539 for (func, f_derivative) in zip(aff_funcs, derivatives_wrt_args): 540 for (var, func_derivative) in func.derivatives.iteritems(): 541 derivatives_wrt_vars[var] += f_derivative * func_derivative 542 543 # The function now returns an AffineScalarFunc object: 544 return AffineScalarFunc(f_nominal_value, derivatives_wrt_vars)

545 546 # It is easier to work with f_with_affine_output, which represents 547 # a wrapped version of 'f', when it bears the same name as 'f': 548 f_with_affine_output.__name__ = f.__name__ 549 550 return f_with_affine_output 551

552 -def _force_aff_func_args(func):

553 """ 554 Takes an operator op(x, y) and wraps it. 555 556 The constructed operator returns func(x, to_affine_scalar(y)) if y 557 can be upcast with to_affine_scalar(); otherwise, it returns 558 NotImplemented. 559 560 Thus, func() is only called on two AffineScalarFunc objects, if 561 its first argument is an AffineScalarFunc. 562 """ 563 564 def op_on_upcast_args(x, y): 565 """ 566 Returns %s(self, to_affine_scalar(y)) if y can be upcast 567 through to_affine_scalar. Otherwise returns NotImplemented. 568 """ % func.__name__ 569 570 try: 571 y_with_uncert = to_affine_scalar(y) 572 except NotUpcast: 573 # This module does not know how to handle the comparison: 574 # (example: y is a NumPy array, in which case the NumPy 575 # array will decide that func() should be applied 576 # element-wise between x and all the elements of y): 577 return NotImplemented 578 else: 579 return func(x, y_with_uncert)

580 581 return op_on_upcast_args 582

583 ######################################## 584 585 # Definition of boolean operators, that assume that self and 586 # y_with_uncert are AffineScalarFunc. 587 588 # The fact that uncertainties must be smalled is used, here: the 589 # comparison functions are supposed to be constant for most values of 590 # the random variables. 591 592 # Even though uncertainties are supposed to be small, comparisons 593 # between 3+/-0.1 and 3.0 are handled (even though x == 3.0 is not a 594 # constant function in the 3+/-0.1 interval). The comparison between 595 # x and x is handled too, when x has an uncertainty. In fact, as 596 # explained in the main documentation, it is possible to give a useful 597 # meaning to the comparison operators, in these cases. 598 599 -def _eq_on_aff_funcs(self, y_with_uncert):

600 """ 601 __eq__ operator, assuming that both self and y_with_uncert are 602 AffineScalarFunc objects. 603 """ 604 difference = self - y_with_uncert 605 # Only an exact zero difference means that self and y are 606 # equal numerically: 607 return not(difference._nominal_value or difference.std_dev())

608

609 -def _ne_on_aff_funcs(self, y_with_uncert):

610 """ 611 __ne__ operator, assuming that both self and y_with_uncert are 612 AffineScalarFunc objects. 613 """ 614 615 return not _eq_on_aff_funcs(self, y_with_uncert)

616

617 -def _gt_on_aff_funcs(self, y_with_uncert):

618 """ 619 __gt__ operator, assuming that both self and y_with_uncert are 620 AffineScalarFunc objects. 621 """ 622 return self._nominal_value > y_with_uncert._nominal_value

623

624 -def _ge_on_aff_funcs(self, y_with_uncert):

625 """ 626 __ge__ operator, assuming that both self and y_with_uncert are 627 AffineScalarFunc objects. 628 """ 629 630 return (_gt_on_aff_funcs(self, y_with_uncert) 631 or _eq_on_aff_funcs(self, y_with_uncert))

632

633 -def _lt_on_aff_funcs(self, y_with_uncert):

634 """ 635 __lt__ operator, assuming that both self and y_with_uncert are 636 AffineScalarFunc objects. 637 """ 638 return self._nominal_value < y_with_uncert._nominal_value

639

640 -def _le_on_aff_funcs(self, y_with_uncert):

641 """ 642 __le__ operator, assuming that both self and y_with_uncert are 643 AffineScalarFunc objects. 644 """ 645 646 return (_lt_on_aff_funcs(self, y_with_uncert) 647 or _eq_on_aff_funcs(self, y_with_uncert))

648

649 ######################################## 650 651 -class AffineScalarFunc(object):

652 """ 653 Affine functions that support basic mathematical operations 654 (addition, etc.). Such functions can for instance be used for 655 representing the local (linear) behavior of any function. 656 657 This class is mostly meant to be used internally. 658 659 This class can also be used to represent constants. 660 661 The variables of affine scalar functions are Variable objects. 662 663 AffineScalarFunc objects include facilities for calculating the 664 'error' on the function, from the uncertainties on its variables. 665 666 Main attributes and methods: 667 668 - nominal_value, std_dev(): value at the origin / nominal value, 669 and standard deviation. 670 671 - error_components(): error_components()[x] is the error due to 672 Variable x. 673 674 - derivatives: derivatives[x] is the (value of the) derivative 675 with respect to Variable x. This attribute is a dictionary 676 whose keys are the Variable objects on which the function 677 depends. 678 679 All the Variable objects on which the function depends are in 680 'derivatives'. 681 682 - position_in_sigmas(x): position of number x with respect to the 683 nominal value, in units of the standard deviation. 684 """ 685 686 # To save memory in large arrays: 687 __slots__ = ('_nominal_value', 'derivatives') 688 689 #! The code could be modify in order to accommodate for non-float 690 # nominal values. This could for instance be done through 691 # the operator module: instead of delegating operations to 692 # float.__*__ operations, they could be delegated to 693 # operator.__*__ functions (while taking care of properly handling 694 # reverse operations: __radd__, etc.). 695

696 - def __init__(self, nominal_value, derivatives):

697 """ 698 nominal_value -- value of the function at the origin. 699 nominal_value must not depend in any way of the Variable 700 objects in 'derivatives' (the value at the origin of the 701 function being defined is a constant). 702 703 derivatives -- maps each Variable object on which the function 704 being defined depends to the value of the derivative with 705 respect to that variable, taken at the nominal value of all 706 variables. 707 708 Warning: the above constraint is not checked, and the user is 709 responsible for complying with it. 710 """ 711 712 # Defines the value at the origin: 713 714 # Only float-like values are handled. One reason is that it 715 # does not make sense for a scalar function to be affine to 716 # not yield float values. Another reason is that it would not 717 # make sense to have a complex nominal value, here (it would 718 # not be handled correctly at all): converting to float should 719 # be possible. 720 self._nominal_value = float(nominal_value) 721 self.derivatives = derivatives

722 723 # The following prevents the 'nominal_value' attribute from being 724 # modified by the user: 725 @property

726 - def nominal_value(self):

727 "Nominal value of the random number." 728 return self._nominal_value

729 730 ############################################################ 731 732 733 ### Operators: operators applied to AffineScalarFunc and/or 734 ### float-like objects only are supported. This is why methods 735 ### from float are used for implementing these operators. 736 737 # Operators with no reflection: 738 739 ######################################## 740 741 # __nonzero__() is supposed to return a boolean value (it is used 742 # by bool()). It is for instance used for converting the result 743 # of comparison operators to a boolean, in sorted(). If we want 744 # to be able to sort AffineScalarFunc objects, __nonzero__ cannot 745 # return a AffineScalarFunc object. Since boolean results (such 746 # as the result of bool()) don't have a very meaningful 747 # uncertainty unless it is zero, this behavior is fine. 748

749 - def __nonzero__(self):

750 """ 751 Equivalent to self != 0. 752 """ 753 #! This might not be relevant for AffineScalarFunc objects 754 # that contain values in a linear space which does not convert 755 # the float 0 into the null vector (see the __eq__ function: 756 # __nonzero__ works fine if subtracting the 0 float from a 757 # vector of the linear space works as if 0 were the null 758 # vector of that space): 759 return self != 0. # Uses the AffineScalarFunc.__ne__ function

760 761 # Compatibility with Python 3: 762 __bool__ = __nonzero__ 763 764 ######################################## 765 766 ## Logical operators: warning: the resulting value cannot always 767 ## be differentiated. 768 769 # The boolean operations are not differentiable everywhere, but 770 # almost... 771 772 # (1) I can rely on the assumption that the user only has "small" 773 # errors on variables, as this is used in the calculation of the 774 # standard deviation (which performs linear approximations): 775 776 # (2) However, this assumption is not relevant for some 777 # operations, and does not have to hold, in some cases. This 778 # comes from the fact that logical operations (e.g. __eq__(x,y)) 779 # are not differentiable for many usual cases. For instance, it 780 # is desirable to have x == x for x = n+/-e, whatever the size of e. 781 # Furthermore, n+/-e != n+/-e', if e != e', whatever the size of e or 782 # e'. 783 784 # (3) The result of logical operators does not have to be a 785 # function with derivatives, as these derivatives are either 0 or 786 # don't exist (i.e., the user should probably not rely on 787 # derivatives for his code). 788 789 # __eq__ is used in "if data in [None, ()]", for instance. It is 790 # therefore important to be able to handle this case too, which is 791 # taken care of when _force_aff_func_args(_eq_on_aff_funcs) 792 # returns NotImplemented. 793 __eq__ = _force_aff_func_args(_eq_on_aff_funcs) 794 795 __ne__ = _force_aff_func_args(_ne_on_aff_funcs) 796 __gt__ = _force_aff_func_args(_gt_on_aff_funcs) 797 798 # __ge__ is not the opposite of __lt__ because these operators do 799 # not always yield a boolean (for instance, 0 <= numpy.arange(10) 800 # yields an array). 801 __ge__ = _force_aff_func_args(_ge_on_aff_funcs) 802 803 __lt__ = _force_aff_func_args(_lt_on_aff_funcs) 804 __le__ = _force_aff_func_args(_le_on_aff_funcs) 805 806 ######################################## 807 808 # Uncertainties handling: 809

810 - def error_components(self):

811 """ 812 Individual components of the standard deviation of the affine 813 function (in absolute value), returned as a dictionary with 814 Variable objects as keys. 815 816 This method assumes that the derivatives contained in the 817 object take scalar values (and are not a tuple, like what 818 math.frexp() returns, for instance). 819 """ 820 821 # Calculation of the variance: 822 error_components = {} 823 for (variable, derivative) in self.derivatives.iteritems(): 824 # Individual standard error due to variable: 825 error_components[variable] = abs(derivative*variable._std_dev) 826 827 return error_components

828

829 - def std_dev(self):

830 """ 831 Standard deviation of the affine function. 832 833 This method assumes that the function returns scalar results. 834 835 This returned standard deviation depends on the current 836 standard deviations [std_dev()] of the variables (Variable 837 objects) involved. 838 """ 839 #! It would be possible to not allow the user to update the 840 #std dev of Variable objects, in which case AffineScalarFunc 841 #objects could have a pre-calculated or, better, cached 842 #std_dev value (in fact, many intermediate AffineScalarFunc do 843 #not need to have their std_dev calculated: only the final 844 #AffineScalarFunc returned to the user does). 845 return sqrt(sum( 846 delta**2 for delta in self.error_components().itervalues()))

847

848 - def _general_representation(self, to_string):

849 """ 850 Uses the to_string() conversion function on both the nominal 851 value and the standard deviation, and returns a string that 852 describes them. 853 854 to_string() is typically repr() or str(). 855 """ 856 857 (nominal_value, std_dev) = (self._nominal_value, self.std_dev()) 858 859 # String representation: 860 861 # Not putting spaces around "+/-" helps with arrays of 862 # Variable, as each value with an uncertainty is a 863 # block of signs (otherwise, the standard deviation can be 864 # mistaken for another element of the array). 865 866 return ("%s+/-%s" % (to_string(nominal_value), to_string(std_dev)) 867 if std_dev 868 else to_string(nominal_value))

869

870 - def __repr__(self):

871 return self._general_representation(repr)

872

873 - def __str__(self):

874 return self._general_representation(str)

875

876 - def position_in_sigmas(self, value):

877 """ 878 Returns 'value' - nominal value, in units of the standard 879 deviation. 880 881 Raises a ValueError exception if the standard deviation is zero. 882 """ 883 try: 884 # The ._nominal_value is a float: there is no integer division, 885 # here: 886 return (value - self._nominal_value) / self.std_dev() 887 except ZeroDivisionError: 888 raise ValueError("The standard deviation is zero:" 889 " undefined result.")

890

891 - def __deepcopy__(self, memo):

892 """ 893 Hook for the standard copy module. 894 895 The returned AffineScalarFunc is a completely fresh copy, 896 which is fully independent of any variable defined so far. 897 New variables are specially created for the returned 898 AffineScalarFunc object. 899 """ 900 return AffineScalarFunc( 901 self._nominal_value, 902 dict((copy.deepcopy(var), deriv) 903 for (var, deriv) in self.derivatives.iteritems()))

904

905 - def __getstate__(self):

906 """ 907 Hook for the pickle module. 908 """ 909 obj_slot_values = dict((k, getattr(self, k)) for k in 910 # self.__slots__ would not work when 911 # self is an instance of a subclass: 912 AffineScalarFunc.__slots__) 913 return obj_slot_values

914

915 - def __setstate__(self, data_dict):

916 """ 917 Hook for the pickle module. 918 """ 919 for (name, value) in data_dict.iteritems(): 920 setattr(self, name, value)

921 922 # Nicer name, for users: isinstance(ufloat(...), UFloat) is True: 923 UFloat = AffineScalarFunc

924 925 -def get_ops_with_reflection():

926 927 """ 928 Returns operators with a reflection, along with their derivatives 929 (for float operands). 930 """ 931 932 # Operators with a reflection: 933 934 # We do not include divmod(). This operator could be included, by 935 # allowing its result (a tuple) to be differentiated, in 936 # derivative_value(). However, a similar result can be achieved 937 # by the user by calculating separately the division and the 938 # result. 939 940 # {operator(x, y): (derivative wrt x, derivative wrt y)}: 941 942 # Note that unknown partial derivatives can be numerically 943 # calculated by expressing them as something like 944 # "partial_derivative(float.__...__, 1)(x, y)": 945 946 # String expressions are used, so that reversed operators are easy 947 # to code, and execute relatively efficiently: 948 949 derivatives_list = { 950 'add': ("1.", "1."), 951 # 'div' is the '/' operator when __future__.division is not in 952 # effect. Since '/' is applied to 953 # AffineScalarFunc._nominal_value numbers, it is applied on 954 # floats, and is therefore the "usual" mathematical division. 955 'div': ("1/y", "-x/y**2"), 956 'floordiv': ("0.", "0."), # Non exact: there is a discontinuities 957 # The derivative wrt the 2nd arguments is something like (..., x//y), 958 # but it is calculated numerically, for convenience: 959 'mod': ("1.", "partial_derivative(float.__mod__, 1)(x, y)"), 960 'mul': ("y", "x"), 961 'pow': ("y*x**(y-1)", "log(x)*x**y"), 962 'sub': ("1.", "-1."), 963 'truediv': ("1/y", "-x/y**2") 964 } 965 966 # Conversion to Python functions: 967 ops_with_reflection = {} 968 for (op, derivatives) in derivatives_list.iteritems(): 969 ops_with_reflection[op] = [ 970 eval("lambda x, y: %s" % expr) for expr in derivatives ] 971 972 ops_with_reflection["r"+op] = [ 973 eval("lambda y, x: %s" % expr) for expr in reversed(derivatives)] 974 975 return ops_with_reflection

976 977 # Operators that have a reflexion, along with their derivatives: 978 _ops_with_reflection = get_ops_with_reflection() 979 980 # Some effectively modified operators (for the automated tests): 981 _modified_operators = []

982 983 -def add_operators_to_AffineScalarFunc():

984 """ 985 Adds many operators (__add__, etc.) to the AffineScalarFunc class. 986 """ 987 988 ######################################## 989 990 #! Derivatives are set to return floats. For one thing, 991 # uncertainties generally involve floats, as they are based on 992 # small variations of the parameters. It is also better to 993 # protect the user from unexpected integer result that behave 994 # badly with the division. 995 996 ## Operators that return a numerical value: 997 998 # Single-argument operators that should be adapted from floats to 999 # AffineScalarFunc objects: 1000 simple_numerical_operators_derivatives = { 1001 'abs': lambda x: 1. if x>=0 else -1., 1002 'neg': lambda x: -1., 1003 'pos': lambda x: 1., 1004 'trunc': lambda x: 0. 1005 } 1006 1007 for (op, derivative) in \ 1008 simple_numerical_operators_derivatives.iteritems(): 1009 1010 attribute_name = "__%s__" % op 1011 # float objects don't exactly have the same attributes between 1012 # different versions of Python (for instance, __trunc__ was 1013 # introduced with Python 2.6): 1014 if attribute_name in dir(float): 1015 setattr(AffineScalarFunc, attribute_name, 1016 wrap(getattr(float, attribute_name), 1017 [derivative])) 1018 _modified_operators.append(op) 1019 1020 ######################################## 1021 1022 # Reversed versions (useful for float*AffineScalarFunc, for instance): 1023 for (op, derivatives) in _ops_with_reflection.iteritems(): 1024 attribute_name = '__%s__' % op 1025 setattr(AffineScalarFunc, attribute_name, 1026 wrap(getattr(float, attribute_name), derivatives)) 1027 1028 ######################################## 1029 # Conversions to pure numbers are meaningless. Note that the 1030 # behavior of float(1j) is similar. 1031 for coercion_type in ('complex', 'int', 'long', 'float'): 1032 def raise_error(self): 1033 raise TypeError("can't convert an affine function (%s)" 1034 ' to %s; use x.nominal_value' 1035 # In case AffineScalarFunc is sub-classed: 1036 % (self.__class__, coercion_type))

1037 1038 setattr(AffineScalarFunc, '__%s__' % coercion_type, raise_error) 1039 1040 add_operators_to_AffineScalarFunc() # Actual addition of class attributes

1041 1042 -class Variable(AffineScalarFunc):

1043 """ 1044 Representation of a float-like scalar random variable, along with 1045 its uncertainty. 1046 """ 1047 1048 # To save memory in large arrays: 1049 __slots__ = ('_std_dev', 'tag') 1050

1051 - def __init__(self, value, std_dev, tag=None):

1052 """ 1053 The nominal value and the standard deviation of the variable 1054 are set. These values must be scalars. 1055 1056 'tag' is a tag that the user can associate to the variable. This 1057 is useful for tracing variables. 1058 1059 The meaning of the nominal value is described in the main 1060 module documentation. 1061 """ 1062 1063 #! The value, std_dev, and tag are assumed by __copy__() not to 1064 # be copied. Either this should be guaranteed here, or __copy__ 1065 # should be updated. 1066 1067 # Only float-like values are handled. One reason is that the 1068 # division operator on integers would not produce a 1069 # differentiable functions: for instance, Variable(3, 0.1)/2 1070 # has a nominal value of 3/2 = 1, but a "shifted" value 1071 # of 3.1/2 = 1.55. 1072 value = float(value) 1073 1074 # If the variable changes by dx, then the value of the affine 1075 # function that gives its value changes by 1*dx: 1076 1077 # ! Memory cycles are created. However, they are garbage 1078 # collected, if possible. Using a weakref.WeakKeyDictionary 1079 # takes much more memory. Thus, this implementation chooses 1080 # more cycles and a smaller memory footprint instead of no 1081 # cycles and a larger memory footprint. 1082 1083 # ! Using AffineScalarFunc instead of super() results only in 1084 # a 3 % speed loss (Python 2.6, Mac OS X): 1085 super(Variable, self).__init__(value, {self: 1.}) 1086 1087 # We force the error to be float-like. Since it is considered 1088 # as a Gaussian standard deviation, it is semantically 1089 # positive (even though there would be no problem defining it 1090 # as a sigma, where sigma can be negative and still define a 1091 # Gaussian): 1092 1093 assert std_dev >= 0, "the error must be a positive number" 1094 # Since AffineScalarFunc.std_dev is a property, we cannot do 1095 # "self.std_dev = ...": 1096 self._std_dev = std_dev 1097 1098 self.tag = tag

1099 1100 # Standard deviations can be modified (this is a feature). 1101 # AffineScalarFunc objects that depend on the Variable have their 1102 # std_dev() automatically modified (recalculated with the new 1103 # std_dev of their Variables):

1104 - def set_std_dev(self, value):

1105 """ 1106 Updates the standard deviation of the variable to a new value. 1107 """ 1108 1109 # A zero variance is accepted. Thus, it is possible to 1110 # conveniently use infinitely precise variables, for instance 1111 # to study special cases. 1112 1113 self._std_dev = value

1114 1115 # The following method is overridden so that we can represent the tag:

1116 - def _general_representation(self, to_string):

1117 """ 1118 Uses the to_string() conversion function on both the nominal 1119 value and standard deviation and returns a string that 1120 describes the number. 1121 1122 to_string() is typically repr() or str(). 1123 """ 1124 num_repr = super(Variable, self)._general_representation(to_string) 1125 1126 # Optional tag: only full representations (to_string == repr) 1127 # contain the tag, as the tag is required in order to recreate 1128 # the variable. Outputting the tag for regular string ("print 1129 # x") would be too heavy and produce an unusual representation 1130 # of a number with uncertainty. 1131 return (num_repr if ((self.tag is None) or (to_string != repr)) 1132 else "< %s = %s >" % (self.tag, num_repr))

1133

1134 - def __copy__(self):

1135 """ 1136 Hook for the standard copy module. 1137 """ 1138 1139 # This copy implicitly takes care of the reference of the 1140 # variable to itself (in self.derivatives): the new Variable 1141 # object points to itself, not to the original Variable. 1142 1143 # Reference: http://www.doughellmann.com/PyMOTW/copy/index.html 1144 1145 #! The following assumes that the arguments to Variable are 1146 # *not* copied upon construction, since __copy__ is not supposed 1147 # to copy "inside" information: 1148 return Variable(self.nominal_value, self.std_dev(), self.tag)

1149

1150 - def __deepcopy__(self, memo):

1151 """ 1152 Hook for the standard copy module. 1153 1154 A new variable is created. 1155 """ 1156 1157 # This deep copy implicitly takes care of the reference of the 1158 # variable to itself (in self.derivatives): the new Variable 1159 # object points to itself, not to the original Variable. 1160 1161 # Reference: http://www.doughellmann.com/PyMOTW/copy/index.html 1162 1163 return self.__copy__()

1164

1165 - def __getstate__(self):

1166 """ 1167 Hook for the standard pickle module. 1168 """ 1169 obj_slot_values = dict((k, getattr(self, k)) for k in self.__slots__) 1170 obj_slot_values.update(AffineScalarFunc.__getstate__(self)) 1171 # Conversion to a usual dictionary: 1172 return obj_slot_values

1173

1174 - def __setstate__(self, data_dict):

1175 """ 1176 Hook for the standard pickle module. 1177 """ 1178 for (name, value) in data_dict.iteritems(): 1179 setattr(self, name, value)

1180

1181 ############################################################################### 1182 1183 # Utilities 1184 1185 -def nominal_value(x):

1186 """ 1187 Returns the nominal value of x if it is a quantity with 1188 uncertainty (i.e., an AffineScalarFunc object); otherwise, returns 1189 x unchanged. 1190 1191 This utility function is useful for transforming a series of 1192 numbers, when only some of them generally carry an uncertainty. 1193 """ 1194 1195 return x.nominal_value if isinstance(x, AffineScalarFunc) else x

1196

1197 -def std_dev(x):

1198 """ 1199 Returns the standard deviation of x if it is a quantity with 1200 uncertainty (i.e., an AffineScalarFunc object); otherwise, returns 1201 the float 0. 1202 1203 This utility function is useful for transforming a series of 1204 numbers, when only some of them generally carry an uncertainty. 1205 """ 1206 1207 return x.std_dev() if isinstance(x, AffineScalarFunc) else 0.

1208

1209 -def covariance_matrix(functions):

1210 """ 1211 Returns a matrix that contains the covariances between the given 1212 sequence of numbers with uncertainties (AffineScalarFunc objects). 1213 The resulting matrix implicitly depends on their ordering in 1214 'functions'. 1215 1216 The covariances are floats (never int objects). 1217 1218 The returned covariance matrix is the exact linear approximation 1219 result, if the nominal values of the functions and of their 1220 variables are their mean. Otherwise, the returned covariance 1221 matrix should be close to it linear approximation value. 1222 """ 1223 # See PSI.411. 1224 1225 covariance_matrix = [] 1226 for (i1, expr1) in enumerate(functions): 1227 derivatives1 = expr1.derivatives # Optimization 1228 vars1 = set(derivatives1) 1229 coefs_expr1 = [] 1230 for (i2, expr2) in enumerate(functions[:i1+1]): 1231 derivatives2 = expr2.derivatives # Optimization 1232 coef = 0. 1233 for var in vars1.intersection(derivatives2): 1234 # var is a variable common to both functions: 1235 coef += (derivatives1[var]*derivatives2[var]*var._std_dev**2) 1236 coefs_expr1.append(coef) 1237 covariance_matrix.append(coefs_expr1) 1238 1239 # We symmetrize the matrix: 1240 for (i, covariance_coefs) in enumerate(covariance_matrix): 1241 covariance_coefs.extend(covariance_matrix[j][i] 1242 for j in range(i+1, len(covariance_matrix))) 1243 1244 return covariance_matrix

1245 1246 # Entering variables as a block of correlated values. Only available 1247 # if NumPy is installed. 1248 1249 #! It would be possible to dispense with NumPy, but a routine should be 1250 # written for obtaining the eigenvectors of a symmetric matrix. See 1251 # for instance Numerical Recipes: (1) reduction to tri-diagonal 1252 # [Givens or Householder]; (2) QR / QL decomposition. 1253 1254 try: 1255 import numpy 1256 except ImportError: 1257 pass 1258 else:

1259 1260 - def correlated_values(values, covariance_mat, tags=None):

1261 """ 1262 Returns numbers with uncertainties (AffineScalarFunc objects) 1263 that correctly reproduce the given covariance matrix, and have 1264 the given values as their nominal value. 1265 1266 The list of values and the covariance matrix must have the 1267 same length, and the matrix must be a square (symmetric) one. 1268 1269 The affine functions returned depend on newly created, 1270 independent variables (Variable objects). 1271 1272 If 'tags' is not None, it must list the tag of each new 1273 independent variable. 1274 """ 1275 1276 # If no tags were given, we prepare tags for the newly created 1277 # variables: 1278 if tags is None: 1279 tags = (None,) * len(values) 1280 1281 # The covariance matrix is diagonalized in order to define 1282 # the independent variables that model the given values: 1283 1284 (variances, transform) = numpy.linalg.eigh(covariance_mat) 1285 1286 # Numerical errors might make some variances negative: we set 1287 # them to zero: 1288 variances[variances < 0] = 0. 1289 1290 # Creation of new, independent variables: 1291 1292 # We use the fact that the eigenvectors in 'transform' are 1293 # special: 'transform' is unitary: its inverse is its transpose: 1294 1295 variables = tuple( 1296 # The variables represent uncertainties only: 1297 Variable(0, sqrt(variance), tag) 1298 for (variance, tag) in zip(variances, tags)) 1299 1300 # Representation of the initial correlated values: 1301 values_funcs = tuple( 1302 AffineScalarFunc(value, dict(zip(variables, coords))) 1303 for (coords, value) in zip(transform, values)) 1304 1305 return values_funcs

1306 1307 __all__.append('correlated_values') 1308 1309 ############################################################################### 1310 # Parsing of values with uncertainties: 1311 1312 POSITIVE_DECIMAL_UNSIGNED = r'(\d+)(\.\d*)?' 1313 1314 # Regexp for a number with uncertainty (e.g., "-1.234(2)e-6"), where the 1315 # uncertainty is optional (in which case the uncertainty is implicit): 1316 NUMBER_WITH_UNCERT_RE_STR = ''' 1317 ([+-])? # Sign 1318 %s # Main number 1319 (?:$%s$)? # Optional uncertainty 1320 ([eE][+-]?\d+)? # Optional exponent 1321 ''' % (POSITIVE_DECIMAL_UNSIGNED, POSITIVE_DECIMAL_UNSIGNED) 1322 1323 NUMBER_WITH_UNCERT_RE = re.compile( 1324 "^%s$" % NUMBER_WITH_UNCERT_RE_STR, re.VERBOSE)

1325 1326 -def parse_error_in_parentheses(representation):

1327 """ 1328 Returns (value, error) from a string representing a number with 1329 uncertainty like 12.34(5), 12.34(142), 12.5(3.4) or 12.3(4.2)e3. 1330 If no parenthesis is given, an uncertainty of one on the last 1331 digit is assumed. 1332 1333 Raises ValueError if the string cannot be parsed. 1334 """ 1335 1336 match = NUMBER_WITH_UNCERT_RE.search(representation) 1337 1338 if match: 1339 # The 'main' part is the nominal value, with 'int'eger part, and 1340 # 'dec'imal part. The 'uncert'ainty is similarly broken into its 1341 # integer and decimal parts. 1342 (sign, main_int, main_dec, uncert_int, uncert_dec, 1343 exponent) = match.groups() 1344 else: 1345 raise ValueError("Unparsable number representation: '%s'." 1346 " Was expecting a string of the form 1.23(4)" 1347 " or 1.234" % representation) 1348 1349 # The value of the number is its nominal value: 1350 value = float(''.join((sign or '', 1351 main_int, 1352 main_dec or '.0', 1353 exponent or ''))) 1354 1355 if uncert_int is None: 1356 # No uncertainty was found: an uncertainty of 1 on the last 1357 # digit is assumed: 1358 uncert_int = '1' 1359 # No uncertainty was found: an uncertainty of 0 on the last 1360 # digit is assumed: 1361 uncert_int = '0' 1362 1363 # Do we have a fully explicit uncertainty? 1364 if uncert_dec is not None: 1365 uncert = float("%s%s" % (uncert_int, uncert_dec or '')) 1366 else: 1367 # uncert_int represents an uncertainty on the last digits: 1368 1369 # The number of digits after the period defines the power of 1370 # 10 than must be applied to the provided uncertainty: 1371 num_digits_after_period = (0 if main_dec is None 1372 else len(main_dec)-1) 1373 uncert = int(uncert_int)/10**num_digits_after_period 1374 1375 # We apply the exponent to the uncertainty as well: 1376 uncert *= float("1%s" % (exponent or '')) 1377 1378 return (value, uncert)

1379

1380 1381 # The following function is not exposed because it can in effect be 1382 # obtained by doing x = ufloat(representation) and 1383 # x.nominal_value and x.std_dev(): 1384 -def str_to_number_with_uncert(representation):

1385 """ 1386 Given a string that represents a number with uncertainty, returns the 1387 nominal value and the uncertainty. 1388 1389 The string can be of the form: 1390 - 124.5+/-0.15 1391 - 124.50(15) 1392 - 124.50(123) 1393 - 124.5 1394 1395 When no numerical error is given, an uncertainty of 1 on the last 1396 digit is implied. 1397 1398 Raises ValueError if the string cannot be parsed. 1399 """ 1400 1401 try: 1402 # Simple form 1234.45+/-1.2: 1403 (value, uncert) = representation.split('+/-') 1404 except ValueError: 1405 # Form with parentheses or no uncertainty: 1406 parsed_value = parse_error_in_parentheses(representation) 1407 else: 1408 try: 1409 parsed_value = (float(value), float(uncert)) 1410 except ValueError: 1411 raise ValueError('Cannot parse %s: was expecting a number' 1412 ' like 1.23+/-0.1' % representation) 1413 1414 return parsed_value

1415

1416 -def ufloat(representation, tag=None):

1417 """ 1418 Returns a random variable (Variable object). 1419 1420 Converts the representation of a number into a number with 1421 uncertainty (a random variable, defined by a nominal value and 1422 a standard deviation). 1423 1424 The representation can be a (value, standard deviation) sequence, 1425 or a string. 1426 1427 Strings of the form '12.345+/-0.015', '12.345(15)', or '12.3' are 1428 recognized (see full list below). In the last case, an 1429 uncertainty of +/-1 is assigned to the last digit. 1430 1431 'tag' is an optional string tag for the variable. Variables 1432 don't have to have distinct tags. Tags are useful for tracing 1433 what values (and errors) enter in a given result (through the 1434 error_components() method). 1435 1436 Examples of valid string representations: 1437 1438 -1.23(3.4) 1439 -1.34(5) 1440 1(6) 1441 3(4.2) 1442 -9(2) 1443 1234567(1.2) 1444 12.345(15) 1445 -12.3456(78)e-6 1446 0.29 1447 31. 1448 -31. 1449 31 1450 -3.1e10 1451 169.0(7) 1452 169.1(15) 1453 """ 1454 1455 # This function is somewhat optimized so as to help with the 1456 # creation of lots of Variable objects (through unumpy.uarray, for 1457 # instance). 1458 1459 # representations is "normalized" so as to be a valid sequence of 1460 # 2 arguments for Variable(). 1461 1462 #! Accepting strings and any kind of sequence slows down the code 1463 # by about 5 %. On the other hand, massive initializations of 1464 # numbers with uncertainties are likely to be performed with 1465 # unumpy.uarray, which does not support parsing from strings and 1466 # thus does not have any overhead. 1467 1468 #! Different, in Python 3: 1469 if isinstance(representation, basestring): 1470 representation = str_to_number_with_uncert(representation) 1471 1472 #! The tag is forced to be a string, so that the user does not 1473 # create a Variable(2.5, 0.5) in order to represent 2.5 1474 # +/- 0.5. Forcing 'tag' to be a string prevents errors from being 1475 # considered as tags, here: 1476 1477 #! 'unicode' is removed in Python3: 1478 if tag is not None: 1479 assert ((type(tag) is str) or (type(tag) is unicode)), \ 1480 "The tag can only be a string." 1481 1482 #! init_args must contain all arguments, here: 1483 return Variable(*representation, **{'tag': tag})

1484

1485 ############################################################################### 1486 # Support for legacy code (will be removed in the future): 1487 1488 -def NumberWithUncert(*args):

1489 """ 1490 Wrapper for legacy code. Obsolete: do not use. Use ufloat 1491 instead. 1492 """ 1493 import warnings 1494 warnings.warn("NumberWithUncert is obsolete." 1495 " Use ufloat instead.", DeprecationWarning, 1496 stacklevel=2) 1497 return ufloat(*args)

1498

1499 -def num_with_uncert(*args):

1500 """ 1501 Wrapper for legacy code. Obsolete: do not use. Use ufloat 1502 instead. 1503 """ 1504 import warnings 1505 warnings.warn("num_with_uncert is obsolete." 1506 " Use ufloat instead.", DeprecationWarning, 1507 stacklevel=2) 1508 return ufloat(*args)

1509

1510 -def array_u(*args):

1511 """ 1512 Wrapper for legacy code. Obsolete: do not use. Use 1513 unumpy.uarray instead. 1514 """ 1515 import warnings 1516 warnings.warn("uncertainties.array_u is obsolete." 1517 " Use uncertainties.unumpy.uarray instead.", 1518 DeprecationWarning, 1519 stacklevel=2) 1520 import uncertainties.unumpy 1521 return uncertainties.unumpy.uarray(*args)

1522

1523 -def nominal_values(*args):

1524 """ 1525 Wrapper for legacy code. Obsolete: do not use. Use 1526 unumpy.nominal_values instead. 1527 """ 1528 import warnings 1529 warnings.warn("uncertainties.nominal_values is obsolete." 1530 " Use uncertainties.unumpy.nominal_values instead.", 1531 DeprecationWarning, 1532 stacklevel=2) 1533 import uncertainties.unumpy 1534 return uncertainties.unumpy.nominal_values(*args)

1535

1536 -def std_devs(*args):

1537 """ 1538 Wrapper for legacy code. Obsolete: do not use. Use ufloat 1539 instead. 1540 """ 1541 import warnings 1542 warnings.warn("uncertainties.std_devs is obsolete." 1543 " Use uncertainties.unumpy.std_devs instead.", 1544 DeprecationWarning, 1545 stacklevel=2) 1546 import uncertainties.unumpy 1547 return uncertainties.unumpy.std_devs(*args)

1548

Source Code for Package ComboCode.cc.ivs.units.uncertainties