ComboCode.cc.ivs.sigproc.lmfit.uncertainties

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

274 275 ############################################################################### 276 277 -def set_doc(doc_string):

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

287 return set_doc_string 288 289 # Some types known to not depend on Variable objects are put in 290 # CONSTANT_TYPES. The most common types can be put in front, as this 291 # may slightly improve the execution speed. 292 CONSTANT_TYPES = (float, int, complex) # , long)

293 294 ############################################################################### 295 # Utility for issuing deprecation warnings 296 297 -def deprecation(message):

298 ''' 299 Warns the user with the given message, by issuing a 300 DeprecationWarning. 301 ''' 302 warnings.warn(message, DeprecationWarning, stacklevel=2)

303 304 305 ############################################################################### 306 307 ## Definitions that depend on the availability of NumPy: 308 309 310 try: 311 import numpy 312 except ImportError: 313 pass 314 else: 315 316 # NumPy numbers do not depend on Variable objects: 317 CONSTANT_TYPES += (numpy.number,)

318 319 # Entering variables as a block of correlated values. Only available 320 # if NumPy is installed. 321 322 #! It would be possible to dispense with NumPy, but a routine should be 323 # written for obtaining the eigenvectors of a symmetric matrix. See 324 # for instance Numerical Recipes: (1) reduction to tri-diagonal 325 # [Givens or Householder]; (2) QR / QL decomposition. 326 327 - def correlated_values(nom_values, covariance_mat, tags=None):

328 """ 329 Returns numbers with uncertainties (AffineScalarFunc objects) 330 that correctly reproduce the given covariance matrix, and have 331 the given (float) values as their nominal value. 332 333 The correlated_values_norm() function returns the same result, 334 but takes a correlation matrix instead of a covariance matrix. 335 336 The list of values and the covariance matrix must have the 337 same length, and the matrix must be a square (symmetric) one. 338 339 The numbers with uncertainties returned depend on newly 340 created, independent variables (Variable objects). 341 342 If 'tags' is not None, it must list the tag of each new 343 independent variable. 344 345 nom_values -- sequence with the nominal (real) values of the 346 numbers with uncertainties to be returned. 347 348 covariance_mat -- full covariance matrix of the returned 349 numbers with uncertainties (not the statistical correlation 350 matrix, i.e., not the normalized covariance matrix). For 351 example, the first element of this matrix is the variance of 352 the first returned number with uncertainty. 353 """ 354 355 # If no tags were given, we prepare tags for the newly created 356 # variables: 357 if tags is None: 358 tags = (None,) * len(nom_values) 359 360 # The covariance matrix is diagonalized in order to define 361 # the independent variables that model the given values: 362 363 (variances, transform) = numpy.linalg.eigh(covariance_mat) 364 365 # Numerical errors might make some variances negative: we set 366 # them to zero: 367 variances[variances < 0] = 0. 368 369 # Creation of new, independent variables: 370 371 # We use the fact that the eigenvectors in 'transform' are 372 # special: 'transform' is unitary: its inverse is its transpose: 373 374 variables = tuple( 375 # The variables represent "pure" uncertainties: 376 Variable(0, sqrt(variance), tag) 377 for (variance, tag) in zip(variances, tags)) 378 379 # Representation of the initial correlated values: 380 values_funcs = tuple( 381 AffineScalarFunc(value, dict(zip(variables, coords))) 382 for (coords, value) in zip(transform, nom_values)) 383 384 return values_funcs

385 386 __all__.append('correlated_values')

387 388 - def correlated_values_norm(values_with_std_dev, correlation_mat, 389 tags=None):

390 ''' 391 Returns correlated values like correlated_values(), but takes 392 instead as input: 393 394 - nominal (float) values along with their standard deviation, and 395 396 - a correlation matrix (i.e. a normalized covariance matrix 397 normalized with individual standard deviations). 398 399 values_with_std_dev -- sequence of (nominal value, standard 400 deviation) pairs. The returned, correlated values have these 401 nominal values and standard deviations. 402 403 correlation_mat -- correlation matrix (i.e. the normalized 404 covariance matrix, a matrix with ones on its diagonal). 405 ''' 406 407 (nominal_values, std_devs) = numpy.transpose(values_with_std_dev) 408 409 return correlated_values( 410 nominal_values, 411 correlation_mat*std_devs*std_devs[numpy.newaxis].T, 412 tags)

413 414 __all__.append('correlated_values_norm')

415 416 ############################################################################### 417 418 # Mathematical operations with local approximations (affine scalar 419 # functions) 420 421 -class NotUpcast(Exception):

422 'Raised when an object cannot be converted to a number with uncertainty'

423

424 -def to_affine_scalar(x):

425 """ 426 Transforms x into a constant affine scalar function 427 (AffineScalarFunc), unless it is already an AffineScalarFunc (in 428 which case x is returned unchanged). 429 430 Raises an exception unless 'x' belongs to some specific classes of 431 objects that are known not to depend on AffineScalarFunc objects 432 (which then cannot be considered as constants). 433 """ 434 435 if isinstance(x, AffineScalarFunc): 436 return x 437 438 #! In Python 2.6+, numbers.Number could be used instead, here: 439 if isinstance(x, CONSTANT_TYPES): 440 # No variable => no derivative to define: 441 return AffineScalarFunc(x, {}) 442 443 # Case of lists, etc. 444 raise NotUpcast("%s cannot be converted to a number with" 445 " uncertainty" % type(x))

446

447 -def partial_derivative(f, param_num):

448 """ 449 Returns a function that numerically calculates the partial 450 derivative of function f with respect to its argument number 451 param_num. 452 453 The step parameter represents the shift of the parameter used in 454 the numerical approximation. 455 """ 456 457 def partial_derivative_of_f(*args, **kws): 458 """ 459 Partial derivative, calculated with the (-epsilon, +epsilon) 460 method, which is more precise than the (0, +epsilon) method. 461 """ 462 # f_nominal_value = f(*args) 463 param_kw = None 464 if '__param__kw__' in kws: 465 param_kw = kws.pop('__param__kw__') 466 shifted_args = list(args) # Copy, and conversion to a mutable 467 shifted_kws = {} 468 for k, v in kws.items(): 469 shifted_kws[k] = v 470 step = 1.e-8 471 if param_kw in shifted_kws: 472 step = step*abs(shifted_kws[param_kw]) 473 elif param_num < len(shifted_args): 474 # The step is relative to the parameter being varied, so that 475 # shsifting it does not suffer from finite precision: 476 step = step*abs(shifted_args[param_num]) 477 478 if param_kw in shifted_kws: 479 shifted_kws[param_kw] += step 480 elif param_num < len(shifted_args): 481 shifted_args[param_num] += step 482 483 shifted_f_plus = f(*shifted_args, **shifted_kws) 484 485 if param_kw in shifted_kws: 486 shifted_kws[param_kw] -= 2*step 487 elif param_num < len(shifted_args): 488 shifted_args[param_num] -= 2*step 489 shifted_f_minus = f(*shifted_args, **shifted_kws) 490 491 return (shifted_f_plus - shifted_f_minus)/2/step

492 493 return partial_derivative_of_f 494

495 -class NumericalDerivatives(object):

496 """ 497 Convenient access to the partial derivatives of a function, 498 calculated numerically. 499 """ 500 # This is not a list because the number of arguments of the 501 # function is not known in advance, in general. 502

503 - def __init__(self, function):

504 """ 505 'function' is the function whose derivatives can be computed. 506 """ 507 self._function = function

508

509 - def __getitem__(self, n):

510 """ 511 Returns the n-th numerical derivative of the function. 512 """ 513 return partial_derivative(self._function, n)

514

515 -def wrap(f, derivatives_iter=None):

516 """ 517 Wraps a function f into a function that also accepts numbers with 518 uncertainties (UFloat objects) and returns a number with 519 uncertainties. Doing so may be necessary when function f cannot 520 be expressed analytically (with uncertainties-compatible operators 521 and functions like +, *, umath.sin(), etc.). 522 523 f must return a scalar (not a list, etc.). 524 525 In the wrapped function, the standard Python scalar arguments of f 526 (float, int, etc.) can be replaced by numbers with 527 uncertainties. The result will contain the appropriate 528 uncertainty. 529 530 If no argument to the wrapped function has an uncertainty, f 531 simply returns its usual, scalar result. 532 533 If supplied, derivatives_iter can be an iterable that generally 534 contains functions; each successive function is the partial 535 derivative of f with respect to the corresponding variable (one 536 function for each argument of f, which takes as many arguments as 537 f). If instead of a function, an element of derivatives_iter 538 contains None, then it is automatically replaced by the relevant 539 numerical derivative; this can be used for non-scalar arguments of 540 f (like string arguments). 541 542 If derivatives_iter is None, or if derivatives_iter contains a 543 fixed (and finite) number of elements, then any missing derivative 544 is calculated numerically. 545 546 An infinite number of derivatives can be specified by having 547 derivatives_iter be an infinite iterator; this can for instance 548 be used for specifying the derivatives of functions with a 549 undefined number of argument (like sum(), whose partial 550 derivatives all return 1). 551 552 Example (for illustration purposes only, as 553 uncertainties.umath.sin() runs faster than the examples that 554 follow): wrap(math.sin) is a sine function that can be applied to 555 numbers with uncertainties. Its derivative will be calculated 556 numerically. wrap(math.sin, [None]) would have produced the same 557 result. wrap(math.sin, [math.cos]) is the same function, but with 558 an analytically defined derivative. 559 """ 560 561 if derivatives_iter is None: 562 derivatives_iter = NumericalDerivatives(f) 563 else: 564 # Derivatives that are not defined are calculated numerically, 565 # if there is a finite number of them (the function lambda 566 # *args: fsum(args) has a non-defined number of arguments, as 567 # it just performs a sum): 568 try: # Is the number of derivatives fixed? 569 len(derivatives_iter) 570 except TypeError: 571 pass 572 else: 573 derivatives_iter = [ 574 partial_derivative(f, k) if derivative is None 575 else derivative 576 for (k, derivative) in enumerate(derivatives_iter)] 577 578 #! Setting the doc string after "def f_with...()" does not 579 # seem to work. We define it explicitly: 580 @set_doc("""\ 581 Version of %s(...) that returns an affine approximation 582 (AffineScalarFunc object), if its result depends on variables 583 (Variable objects). Otherwise, returns a simple constant (when 584 applied to constant arguments). 585 586 Warning: arguments of the function that are not AffineScalarFunc 587 objects must not depend on uncertainties.Variable objects in any 588 way. Otherwise, the dependence of the result in 589 uncertainties.Variable objects will be incorrect. 590 591 Original documentation: 592 %s""" % (f.__name__, f.__doc__)) 593 def f_with_affine_output(*args, **kwargs): 594 # Can this function perform the calculation of an 595 # AffineScalarFunc (or maybe float) result? 596 try: 597 old_funcs = map(to_affine_scalar, args) 598 aff_funcs = [to_affine_scalar(a) for a in args] 599 aff_kws = kwargs 600 aff_varkws = [] 601 for key, val in kwargs.items(): 602 if isinstance(val, Variable): 603 aff_kws[key] = to_affine_scalar(val) 604 aff_varkws.append(key) 605 606 except NotUpcast: 607 608 # This function does not know how to itself perform 609 # calculations with non-float-like arguments (as they 610 # might for instance be objects whose value really changes 611 # if some Variable objects had different values): 612 613 # Is it clear that we can't delegate the calculation? 614 615 if any(isinstance(arg, AffineScalarFunc) for arg in args): 616 # This situation arises for instance when calculating 617 # AffineScalarFunc(...)*numpy.array(...). In this 618 # case, we must let NumPy handle the multiplication 619 # (which is then performed element by element): 620 return NotImplemented 621 else: 622 # If none of the arguments is an AffineScalarFunc, we 623 # can delegate the calculation to the original 624 # function. This can be useful when it is called with 625 # only one argument (as in 626 # numpy.log10(numpy.ndarray(...)): 627 return f(*args, **kwargs) 628 629 ######################################## 630 # Nominal value of the constructed AffineScalarFunc: 631 args_values = [e.nominal_value for e in aff_funcs] 632 kw_values = {} 633 for key, val in aff_kws.items(): 634 kw_values[key] = val 635 if key in aff_varkws: 636 kw_values[key] = val.nominal_value 637 f_nominal_value = f(*args_values, **kw_values) 638 639 ######################################## 640 641 # List of involved variables (Variable objects): 642 variables = set() 643 for expr in aff_funcs: 644 variables |= set(expr.derivatives) 645 for vname in aff_varkws: 646 variables |= set(aff_kws[vname].derivatives) 647 ## It is sometimes useful to only return a regular constant: 648 649 # (1) Optimization / convenience behavior: when 'f' is called 650 # on purely constant values (e.g., sin(2)), there is no need 651 # for returning a more complex AffineScalarFunc object. 652 653 # (2) Functions that do not return a "float-like" value might 654 # not have a relevant representation as an AffineScalarFunc. 655 # This includes boolean functions, since their derivatives are 656 # either 0 or are undefined: they are better represented as 657 # Python constants than as constant AffineScalarFunc functions. 658 659 if not variables or isinstance(f_nominal_value, bool): 660 return f_nominal_value 661 662 # The result of 'f' does depend on 'variables'... 663 664 ######################################## 665 666 # Calculation of the derivatives with respect to the arguments 667 # of f (aff_funcs): 668 669 # The chain rule is applied. This is because, in the case of 670 # numerical derivatives, it allows for a better-controlled 671 # numerical stability than numerically calculating the partial 672 # derivatives through '[f(x + dx, y + dy, ...) - 673 # f(x,y,...)]/da' where dx, dy,... are calculated by varying 674 # 'a'. In fact, it is numerically better to control how big 675 # (dx, dy,...) are: 'f' is a simple mathematical function and 676 # it is possible to know how precise the df/dx are (which is 677 # not possible with the numerical df/da calculation above). 678 679 # We use numerical derivatives, if we don't already have a 680 # list of derivatives: 681 682 #! Note that this test could be avoided by requiring the 683 # caller to always provide derivatives. When changing the 684 # functions of the math module, this would force this module 685 # to know about all the math functions. Another possibility 686 # would be to force derivatives_iter to contain, say, the 687 # first 3 derivatives of f. But any of these two ideas has a 688 # chance to break, one day... (if new functions are added to 689 # the math module, or if some function has more than 3 690 # arguments). 691 692 derivatives_wrt_args = [] 693 for (arg, derivative) in zip(aff_funcs, derivatives_iter): 694 derivatives_wrt_args.append(derivative(*args_values, **aff_kws) 695 if arg.derivatives 696 else 0) 697 698 699 kws_values = [] 700 for vname in aff_varkws: 701 kws_values.append( aff_kws[vname].nominal_value) 702 for (vname, derivative) in zip(aff_varkws, derivatives_iter): 703 derivatives_wrt_args.append(derivative(__param__kw__=vname, 704 **kw_values) 705 if aff_kws[vname].derivatives 706 else 0) 707 708 ######################################## 709 # Calculation of the derivative of f with respect to all the 710 # variables (Variable) involved. 711 712 # Initial value (is updated below): 713 derivatives_wrt_vars = dict((var, 0.) for var in variables) 714 715 # The chain rule is used (we already have 716 # derivatives_wrt_args): 717 718 for (func, f_derivative) in zip(aff_funcs, derivatives_wrt_args): 719 for (var, func_derivative) in func.derivatives.items(): 720 derivatives_wrt_vars[var] += f_derivative * func_derivative 721 722 for (vname, f_derivative) in zip(aff_varkws, derivatives_wrt_args): 723 func = aff_kws[vname] 724 for (var, func_derivative) in func.derivatives.items(): 725 derivatives_wrt_vars[var] += f_derivative * func_derivative 726 727 # The function now returns an AffineScalarFunc object: 728 return AffineScalarFunc(f_nominal_value, derivatives_wrt_vars)

729 730 # It is easier to work with f_with_affine_output, which represents 731 # a wrapped version of 'f', when it bears the same name as 'f': 732 f_with_affine_output.__name__ = f.__name__ 733 734 return f_with_affine_output 735

736 -def _force_aff_func_args(func):

737 """ 738 Takes an operator op(x, y) and wraps it. 739 740 The constructed operator returns func(x, to_affine_scalar(y)) if y 741 can be upcast with to_affine_scalar(); otherwise, it returns 742 NotImplemented. 743 744 Thus, func() is only called on two AffineScalarFunc objects, if 745 its first argument is an AffineScalarFunc. 746 """ 747 748 def op_on_upcast_args(x, y): 749 """ 750 Returns %s(self, to_affine_scalar(y)) if y can be upcast 751 through to_affine_scalar. Otherwise returns NotImplemented. 752 """ % func.__name__ 753 754 try: 755 y_with_uncert = to_affine_scalar(y) 756 except NotUpcast: 757 # This module does not know how to handle the comparison: 758 # (example: y is a NumPy array, in which case the NumPy 759 # array will decide that func() should be applied 760 # element-wise between x and all the elements of y): 761 return NotImplemented 762 else: 763 return func(x, y_with_uncert)

764 765 return op_on_upcast_args 766

767 ######################################## 768 769 # Definition of boolean operators, that assume that self and 770 # y_with_uncert are AffineScalarFunc. 771 772 # The fact that uncertainties must be smalled is used, here: the 773 # comparison functions are supposed to be constant for most values of 774 # the random variables. 775 776 # Even though uncertainties are supposed to be small, comparisons 777 # between 3+/-0.1 and 3.0 are handled (even though x == 3.0 is not a 778 # constant function in the 3+/-0.1 interval). The comparison between 779 # x and x is handled too, when x has an uncertainty. In fact, as 780 # explained in the main documentation, it is possible to give a useful 781 # meaning to the comparison operators, in these cases. 782 783 -def _eq_on_aff_funcs(self, y_with_uncert):

784 """ 785 __eq__ operator, assuming that both self and y_with_uncert are 786 AffineScalarFunc objects. 787 """ 788 difference = self - y_with_uncert 789 # Only an exact zero difference means that self and y are 790 # equal numerically: 791 return not(difference._nominal_value or difference.std_dev())

792

793 -def _ne_on_aff_funcs(self, y_with_uncert):

794 """ 795 __ne__ operator, assuming that both self and y_with_uncert are 796 AffineScalarFunc objects. 797 """ 798 799 return not _eq_on_aff_funcs(self, y_with_uncert)

800

801 -def _gt_on_aff_funcs(self, y_with_uncert):

802 """ 803 __gt__ operator, assuming that both self and y_with_uncert are 804 AffineScalarFunc objects. 805 """ 806 return self._nominal_value > y_with_uncert._nominal_value

807

808 -def _ge_on_aff_funcs(self, y_with_uncert):

809 """ 810 __ge__ operator, assuming that both self and y_with_uncert are 811 AffineScalarFunc objects. 812 """ 813 814 return (_gt_on_aff_funcs(self, y_with_uncert) 815 or _eq_on_aff_funcs(self, y_with_uncert))

816

817 -def _lt_on_aff_funcs(self, y_with_uncert):

818 """ 819 __lt__ operator, assuming that both self and y_with_uncert are 820 AffineScalarFunc objects. 821 """ 822 return self._nominal_value < y_with_uncert._nominal_value

823

824 -def _le_on_aff_funcs(self, y_with_uncert):

825 """ 826 __le__ operator, assuming that both self and y_with_uncert are 827 AffineScalarFunc objects. 828 """ 829 830 return (_lt_on_aff_funcs(self, y_with_uncert) 831 or _eq_on_aff_funcs(self, y_with_uncert))

832

833 ######################################## 834 835 -class AffineScalarFunc(object):

836 """ 837 Affine functions that support basic mathematical operations 838 (addition, etc.). Such functions can for instance be used for 839 representing the local (linear) behavior of any function. 840 841 This class is mostly meant to be used internally. 842 843 This class can also be used to represent constants. 844 845 The variables of affine scalar functions are Variable objects. 846 847 AffineScalarFunc objects include facilities for calculating the 848 'error' on the function, from the uncertainties on its variables. 849 850 Main attributes and methods: 851 852 - nominal_value, std_dev(): value at the origin / nominal value, 853 and standard deviation. 854 855 - error_components(): error_components()[x] is the error due to 856 Variable x. 857 858 - derivatives: derivatives[x] is the (value of the) derivative 859 with respect to Variable x. This attribute is a dictionary 860 whose keys are the Variable objects on which the function 861 depends. 862 863 All the Variable objects on which the function depends are in 864 'derivatives'. 865 866 - std_score(x): position of number x with respect to the 867 nominal value, in units of the standard deviation. 868 """ 869 870 # To save memory in large arrays: 871 __slots__ = ('_nominal_value', 'derivatives') 872 873 #! The code could be modify in order to accommodate for non-float 874 # nominal values. This could for instance be done through 875 # the operator module: instead of delegating operations to 876 # float.__*__ operations, they could be delegated to 877 # operator.__*__ functions (while taking care of properly handling 878 # reverse operations: __radd__, etc.). 879

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

881 """ 882 nominal_value -- value of the function at the origin. 883 nominal_value must not depend in any way of the Variable 884 objects in 'derivatives' (the value at the origin of the 885 function being defined is a constant). 886 887 derivatives -- maps each Variable object on which the function 888 being defined depends to the value of the derivative with 889 respect to that variable, taken at the nominal value of all 890 variables. 891 892 Warning: the above constraint is not checked, and the user is 893 responsible for complying with it. 894 """ 895 896 # Defines the value at the origin: 897 898 # Only float-like values are handled. One reason is that it 899 # does not make sense for a scalar function to be affine to 900 # not yield float values. Another reason is that it would not 901 # make sense to have a complex nominal value, here (it would 902 # not be handled correctly at all): converting to float should 903 # be possible. 904 self._nominal_value = float(nominal_value) 905 self.derivatives = derivatives

906 907 # The following prevents the 'nominal_value' attribute from being 908 # modified by the user: 909 @property

910 - def nominal_value(self):

911 "Nominal value of the random number." 912 return self._nominal_value

913 914 ############################################################ 915 916 917 ### Operators: operators applied to AffineScalarFunc and/or 918 ### float-like objects only are supported. This is why methods 919 ### from float are used for implementing these operators. 920 921 # Operators with no reflection: 922 923 ######################################## 924 925 # __nonzero__() is supposed to return a boolean value (it is used 926 # by bool()). It is for instance used for converting the result 927 # of comparison operators to a boolean, in sorted(). If we want 928 # to be able to sort AffineScalarFunc objects, __nonzero__ cannot 929 # return a AffineScalarFunc object. Since boolean results (such 930 # as the result of bool()) don't have a very meaningful 931 # uncertainty unless it is zero, this behavior is fine. 932

933 - def __nonzero__(self):

934 """ 935 Equivalent to self != 0. 936 """ 937 #! This might not be relevant for AffineScalarFunc objects 938 # that contain values in a linear space which does not convert 939 # the float 0 into the null vector (see the __eq__ function: 940 # __nonzero__ works fine if subtracting the 0 float from a 941 # vector of the linear space works as if 0 were the null 942 # vector of that space): 943 return self != 0. # Uses the AffineScalarFunc.__ne__ function

944 945 ######################################## 946 947 ## Logical operators: warning: the resulting value cannot always 948 ## be differentiated. 949 950 # The boolean operations are not differentiable everywhere, but 951 # almost... 952 953 # (1) I can rely on the assumption that the user only has "small" 954 # errors on variables, as this is used in the calculation of the 955 # standard deviation (which performs linear approximations): 956 957 # (2) However, this assumption is not relevant for some 958 # operations, and does not have to hold, in some cases. This 959 # comes from the fact that logical operations (e.g. __eq__(x,y)) 960 # are not differentiable for many usual cases. For instance, it 961 # is desirable to have x == x for x = n+/-e, whatever the size of e. 962 # Furthermore, n+/-e != n+/-e', if e != e', whatever the size of e or 963 # e'. 964 965 # (3) The result of logical operators does not have to be a 966 # function with derivatives, as these derivatives are either 0 or 967 # don't exist (i.e., the user should probably not rely on 968 # derivatives for his code). 969 970 # __eq__ is used in "if data in [None, ()]", for instance. It is 971 # therefore important to be able to handle this case too, which is 972 # taken care of when _force_aff_func_args(_eq_on_aff_funcs) 973 # returns NotImplemented. 974 __eq__ = _force_aff_func_args(_eq_on_aff_funcs) 975 976 __ne__ = _force_aff_func_args(_ne_on_aff_funcs) 977 __gt__ = _force_aff_func_args(_gt_on_aff_funcs) 978 979 # __ge__ is not the opposite of __lt__ because these operators do 980 # not always yield a boolean (for instance, 0 <= numpy.arange(10) 981 # yields an array). 982 __ge__ = _force_aff_func_args(_ge_on_aff_funcs) 983 984 __lt__ = _force_aff_func_args(_lt_on_aff_funcs) 985 __le__ = _force_aff_func_args(_le_on_aff_funcs) 986 987 ######################################## 988 989 # Uncertainties handling: 990

991 - def error_components(self):

992 """ 993 Individual components of the standard deviation of the affine 994 function (in absolute value), returned as a dictionary with 995 Variable objects as keys. 996 997 This method assumes that the derivatives contained in the 998 object take scalar values (and are not a tuple, like what 999 math.frexp() returns, for instance). 1000 """ 1001 1002 # Calculation of the variance: 1003 error_components = {} 1004 for (variable, derivative) in self.derivatives.items(): 1005 # Individual standard error due to variable: 1006 error_components[variable] = abs(derivative*variable._std_dev) 1007 1008 return error_components

1009

1010 - def std_dev(self):

1011 """ 1012 Standard deviation of the affine function. 1013 1014 This method assumes that the function returns scalar results. 1015 1016 This returned standard deviation depends on the current 1017 standard deviations [std_dev()] of the variables (Variable 1018 objects) involved. 1019 """ 1020 #! It would be possible to not allow the user to update the 1021 #std dev of Variable objects, in which case AffineScalarFunc 1022 #objects could have a pre-calculated or, better, cached 1023 #std_dev value (in fact, many intermediate AffineScalarFunc do 1024 #not need to have their std_dev calculated: only the final 1025 #AffineScalarFunc returned to the user does). 1026 return sqrt(sum( 1027 delta**2 for delta in self.error_components().itervalues()))

1028

1029 - def _general_representation(self, to_string):

1030 """ 1031 Uses the to_string() conversion function on both the nominal 1032 value and the standard deviation, and returns a string that 1033 describes them. 1034 1035 to_string() is typically repr() or str(). 1036 """ 1037 1038 (nominal_value, std_dev) = (self._nominal_value, self.std_dev()) 1039 1040 # String representation: 1041 1042 # Not putting spaces around "+/-" helps with arrays of 1043 # Variable, as each value with an uncertainty is a 1044 # block of signs (otherwise, the standard deviation can be 1045 # mistaken for another element of the array). 1046 1047 return ("%s+/-%s" % (to_string(nominal_value), to_string(std_dev)) 1048 if std_dev 1049 else to_string(nominal_value))

1050

1051 - def __repr__(self):

1052 return self._general_representation(repr)

1053

1054 - def __str__(self):

1055 return self._general_representation(str)

1056

1057 - def std_score(self, value):

1058 """ 1059 Returns 'value' - nominal value, in units of the standard 1060 deviation. 1061 1062 Raises a ValueError exception if the standard deviation is zero. 1063 """ 1064 try: 1065 # The ._nominal_value is a float: there is no integer division, 1066 # here: 1067 return (value - self._nominal_value) / self.std_dev() 1068 except ZeroDivisionError: 1069 raise ValueError("The standard deviation is zero:" 1070 " undefined result.")

1071

1072 - def __deepcopy__(self, memo):

1073 """ 1074 Hook for the standard copy module. 1075 1076 The returned AffineScalarFunc is a completely fresh copy, 1077 which is fully independent of any variable defined so far. 1078 New variables are specially created for the returned 1079 AffineScalarFunc object. 1080 """ 1081 return AffineScalarFunc( 1082 self._nominal_value, 1083 dict((copy.deepcopy(var), deriv) 1084 for (var, deriv) in self.derivatives.items()))

1085

1086 - def __getstate__(self):

1087 """ 1088 Hook for the pickle module. 1089 """ 1090 obj_slot_values = dict((k, getattr(self, k)) for k in 1091 # self.__slots__ would not work when 1092 # self is an instance of a subclass: 1093 AffineScalarFunc.__slots__) 1094 return obj_slot_values

1095

1096 - def __setstate__(self, data_dict):

1097 """ 1098 Hook for the pickle module. 1099 """ 1100 for (name, value) in data_dict.items(): 1101 setattr(self, name, value)

1102 1103 # Nicer name, for users: isinstance(ufloat(...), UFloat) is True: 1104 UFloat = AffineScalarFunc

1105 1106 -def get_ops_with_reflection():

1107 1108 """ 1109 Returns operators with a reflection, along with their derivatives 1110 (for float operands). 1111 """ 1112 1113 # Operators with a reflection: 1114 1115 # We do not include divmod(). This operator could be included, by 1116 # allowing its result (a tuple) to be differentiated, in 1117 # derivative_value(). However, a similar result can be achieved 1118 # by the user by calculating separately the division and the 1119 # result. 1120 1121 # {operator(x, y): (derivative wrt x, derivative wrt y)}: 1122 1123 # Note that unknown partial derivatives can be numerically 1124 # calculated by expressing them as something like 1125 # "partial_derivative(float.__...__, 1)(x, y)": 1126 1127 # String expressions are used, so that reversed operators are easy 1128 # to code, and execute relatively efficiently: 1129 1130 derivatives_list = { 1131 'add': ("1.", "1."), 1132 # 'div' is the '/' operator when __future__.division is not in 1133 # effect. Since '/' is applied to 1134 # AffineScalarFunc._nominal_value numbers, it is applied on 1135 # floats, and is therefore the "usual" mathematical division. 1136 'div': ("1/y", "-x/y**2"), 1137 'floordiv': ("0.", "0."), # Non exact: there is a discontinuities 1138 # The derivative wrt the 2nd arguments is something like (..., x//y), 1139 # but it is calculated numerically, for convenience: 1140 'mod': ("1.", "partial_derivative(float.__mod__, 1)(x, y)"), 1141 'mul': ("y", "x"), 1142 'pow': ("y*x**(y-1)", "log(x)*x**y"), 1143 'sub': ("1.", "-1."), 1144 'truediv': ("1/y", "-x/y**2") 1145 } 1146 1147 # Conversion to Python functions: 1148 ops_with_reflection = {} 1149 for (op, derivatives) in derivatives_list.items(): 1150 ops_with_reflection[op] = [ 1151 eval("lambda x, y: %s" % expr) for expr in derivatives ] 1152 1153 ops_with_reflection["r"+op] = [ 1154 eval("lambda y, x: %s" % expr) for expr in reversed(derivatives)] 1155 1156 return ops_with_reflection

1157 1158 # Operators that have a reflection, along with their derivatives: 1159 _ops_with_reflection = get_ops_with_reflection() 1160 1161 # Some effectively modified operators (for the automated tests): 1162 _modified_operators = [] 1163 _modified_ops_with_reflection = []

1164 1165 -def add_operators_to_AffineScalarFunc():

1166 """ 1167 Adds many operators (__add__, etc.) to the AffineScalarFunc class. 1168 """ 1169 1170 ######################################## 1171 1172 #! Derivatives are set to return floats. For one thing, 1173 # uncertainties generally involve floats, as they are based on 1174 # small variations of the parameters. It is also better to 1175 # protect the user from unexpected integer result that behave 1176 # badly with the division. 1177 1178 ## Operators that return a numerical value: 1179 1180 # Single-argument operators that should be adapted from floats to 1181 # AffineScalarFunc objects, associated to their derivative: 1182 simple_numerical_operators_derivatives = { 1183 'abs': lambda x: 1. if x>=0 else -1., 1184 'neg': lambda x: -1., 1185 'pos': lambda x: 1., 1186 'trunc': lambda x: 0. 1187 } 1188 1189 for (op, derivative) in ( 1190 simple_numerical_operators_derivatives.items()): 1191 1192 attribute_name = "__%s__" % op 1193 # float objects don't exactly have the same attributes between 1194 # different versions of Python (for instance, __trunc__ was 1195 # introduced with Python 2.6): 1196 try: 1197 setattr(AffineScalarFunc, attribute_name, 1198 wrap(getattr(float, attribute_name), 1199 [derivative])) 1200 except AttributeError: 1201 pass 1202 else: 1203 _modified_operators.append(op) 1204 1205 ######################################## 1206 1207 # Reversed versions (useful for float*AffineScalarFunc, for instance): 1208 for (op, derivatives) in _ops_with_reflection.items(): 1209 attribute_name = '__%s__' % op 1210 # float objects don't exactly have the same attributes between 1211 # different versions of Python (for instance, __div__ and 1212 # __rdiv__ were removed, in Python 3): 1213 try: 1214 setattr(AffineScalarFunc, attribute_name, 1215 wrap(getattr(float, attribute_name), derivatives)) 1216 except AttributeError: 1217 pass 1218 else: 1219 _modified_ops_with_reflection.append(op) 1220 1221 ######################################## 1222 # Conversions to pure numbers are meaningless. Note that the 1223 # behavior of float(1j) is similar. 1224 for coercion_type in ('complex', 'int', 'long', 'float'): 1225 def raise_error(self): 1226 raise TypeError("can't convert an affine function (%s)" 1227 ' to %s; use x.nominal_value' 1228 # In case AffineScalarFunc is sub-classed: 1229 % (self.__class__, coercion_type))

1230 1231 setattr(AffineScalarFunc, '__%s__' % coercion_type, raise_error) 1232 1233 add_operators_to_AffineScalarFunc() # Actual addition of class attributes

1234 1235 -class Variable(AffineScalarFunc):

1236 """ 1237 Representation of a float-like scalar random variable, along with 1238 its uncertainty. 1239 1240 Objects are meant to represent variables that are independent from 1241 each other (correlations are handled through the AffineScalarFunc 1242 class). 1243 """ 1244 1245 # To save memory in large arrays: 1246 __slots__ = ('_std_dev', 'tag') 1247

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

1249 """ 1250 The nominal value and the standard deviation of the variable 1251 are set. These values must be scalars. 1252 1253 'tag' is a tag that the user can associate to the variable. This 1254 is useful for tracing variables. 1255 1256 The meaning of the nominal value is described in the main 1257 module documentation. 1258 """ 1259 1260 #! The value, std_dev, and tag are assumed by __copy__() not to 1261 # be copied. Either this should be guaranteed here, or __copy__ 1262 # should be updated. 1263 1264 # Only float-like values are handled. One reason is that the 1265 # division operator on integers would not produce a 1266 # differentiable functions: for instance, Variable(3, 0.1)/2 1267 # has a nominal value of 3/2 = 1, but a "shifted" value 1268 # of 3.1/2 = 1.55. 1269 value = float(value) 1270 1271 # If the variable changes by dx, then the value of the affine 1272 # function that gives its value changes by 1*dx: 1273 1274 # ! Memory cycles are created. However, they are garbage 1275 # collected, if possible. Using a weakref.WeakKeyDictionary 1276 # takes much more memory. Thus, this implementation chooses 1277 # more cycles and a smaller memory footprint instead of no 1278 # cycles and a larger memory footprint. 1279 1280 # ! Using AffineScalarFunc instead of super() results only in 1281 # a 3 % speed loss (Python 2.6, Mac OS X): 1282 super(Variable, self).__init__(value, {self: 1.}) 1283 1284 # We force the error to be float-like. Since it is considered 1285 # as a Gaussian standard deviation, it is semantically 1286 # positive (even though there would be no problem defining it 1287 # as a sigma, where sigma can be negative and still define a 1288 # Gaussian): 1289 1290 assert std_dev >= 0, "the error must be a positive number" 1291 # Since AffineScalarFunc.std_dev is a property, we cannot do 1292 # "self.std_dev = ...": 1293 self._std_dev = std_dev 1294 1295 self.tag = tag

1296 1297 # Standard deviations can be modified (this is a feature). 1298 # AffineScalarFunc objects that depend on the Variable have their 1299 # std_dev() automatically modified (recalculated with the new 1300 # std_dev of their Variables):

1301 - def set_std_dev(self, value):

1302 """ 1303 Updates the standard deviation of the variable to a new value. 1304 """ 1305 1306 # A zero variance is accepted. Thus, it is possible to 1307 # conveniently use infinitely precise variables, for instance 1308 # to study special cases. 1309 1310 self._std_dev = value

1311 1312 # The following method is overridden so that we can represent the tag:

1313 - def _general_representation(self, to_string):

1314 """ 1315 Uses the to_string() conversion function on both the nominal 1316 value and standard deviation and returns a string that 1317 describes the number. 1318 1319 to_string() is typically repr() or str(). 1320 """ 1321 num_repr = super(Variable, self)._general_representation(to_string) 1322 1323 # Optional tag: only full representations (to_string == repr) 1324 # contain the tag, as the tag is required in order to recreate 1325 # the variable. Outputting the tag for regular string ("print 1326 # x") would be too heavy and produce an unusual representation 1327 # of a number with uncertainty. 1328 return (num_repr if ((self.tag is None) or (to_string != repr)) 1329 else "< %s = %s >" % (self.tag, num_repr))

1330

1331 - def __hash__(self):

1332 # All Variable objects are by definition independent 1333 # variables, so they never compare equal; therefore, their 1334 # id() are therefore allowed to differ 1335 # (http://docs.python.org/reference/datamodel.html#object.__hash__): 1336 return id(self)

1337

1338 - def __copy__(self):

1339 """ 1340 Hook for the standard copy module. 1341 """ 1342 1343 # This copy implicitly takes care of the reference of the 1344 # variable to itself (in self.derivatives): the new Variable 1345 # object points to itself, not to the original Variable. 1346 1347 # Reference: http://www.doughellmann.com/PyMOTW/copy/index.html 1348 1349 #! The following assumes that the arguments to Variable are 1350 # *not* copied upon construction, since __copy__ is not supposed 1351 # to copy "inside" information: 1352 return Variable(self.nominal_value, self.std_dev(), self.tag)

1353

1354 - def __deepcopy__(self, memo):

1355 """ 1356 Hook for the standard copy module. 1357 1358 A new variable is created. 1359 """ 1360 1361 # This deep copy implicitly takes care of the reference of the 1362 # variable to itself (in self.derivatives): the new Variable 1363 # object points to itself, not to the original Variable. 1364 1365 # Reference: http://www.doughellmann.com/PyMOTW/copy/index.html 1366 1367 return self.__copy__()

1368

1369 - def __getstate__(self):

1370 """ 1371 Hook for the standard pickle module. 1372 """ 1373 obj_slot_values = dict((k, getattr(self, k)) for k in self.__slots__) 1374 obj_slot_values.update(AffineScalarFunc.__getstate__(self)) 1375 # Conversion to a usual dictionary: 1376 return obj_slot_values

1377

1378 - def __setstate__(self, data_dict):

1379 """ 1380 Hook for the standard pickle module. 1381 """ 1382 for (name, value) in data_dict.items(): 1383 setattr(self, name, value)

1384

1385 ############################################################################### 1386 1387 # Utilities 1388 1389 -def nominal_value(x):

1390 """ 1391 Returns the nominal value of x if it is a quantity with 1392 uncertainty (i.e., an AffineScalarFunc object); otherwise, returns 1393 x unchanged. 1394 1395 This utility function is useful for transforming a series of 1396 numbers, when only some of them generally carry an uncertainty. 1397 """ 1398 1399 return x.nominal_value if isinstance(x, AffineScalarFunc) else x

1400

1401 -def std_dev(x):

1402 """ 1403 Returns the standard deviation of x if it is a quantity with 1404 uncertainty (i.e., an AffineScalarFunc object); otherwise, returns 1405 the float 0. 1406 1407 This utility function is useful for transforming a series of 1408 numbers, when only some of them generally carry an uncertainty. 1409 """ 1410 1411 return x.std_dev() if isinstance(x, AffineScalarFunc) else 0.

1412

1413 -def covariance_matrix(nums_with_uncert):

1414 """ 1415 Returns a matrix that contains the covariances between the given 1416 sequence of numbers with uncertainties (AffineScalarFunc objects). 1417 The resulting matrix implicitly depends on their ordering in 1418 'nums_with_uncert'. 1419 1420 The covariances are floats (never int objects). 1421 1422 The returned covariance matrix is the exact linear approximation 1423 result, if the nominal values of the numbers with uncertainties 1424 and of their variables are their mean. Otherwise, the returned 1425 covariance matrix should be close to its linear approximation 1426 value. 1427 1428 The returned matrix is a list of lists. 1429 """ 1430 # See PSI.411 in EOL's notes. 1431 1432 covariance_matrix = [] 1433 for (i1, expr1) in enumerate(nums_with_uncert): 1434 derivatives1 = expr1.derivatives # Optimization 1435 vars1 = set(derivatives1) 1436 coefs_expr1 = [] 1437 for (i2, expr2) in enumerate(nums_with_uncert[:i1+1]): 1438 derivatives2 = expr2.derivatives # Optimization 1439 coef = 0. 1440 for var in vars1.intersection(derivatives2): 1441 # var is a variable common to both numbers with 1442 # uncertainties: 1443 coef += (derivatives1[var]*derivatives2[var]*var._std_dev**2) 1444 coefs_expr1.append(coef) 1445 covariance_matrix.append(coefs_expr1) 1446 1447 # We symmetrize the matrix: 1448 for (i, covariance_coefs) in enumerate(covariance_matrix): 1449 covariance_coefs.extend(covariance_matrix[j][i] 1450 for j in range(i+1, len(covariance_matrix))) 1451 1452 return covariance_matrix

1453 1454 try: 1455 import numpy 1456 except ImportError: 1457 pass 1458 else:

1459 - def correlation_matrix(nums_with_uncert):

1460 ''' 1461 Returns the correlation matrix of the given sequence of 1462 numbers with uncertainties, as a NumPy array of floats. 1463 ''' 1464 1465 cov_mat = numpy.array(covariance_matrix(nums_with_uncert)) 1466 1467 std_devs = numpy.sqrt(cov_mat.diagonal()) 1468 1469 return cov_mat/std_devs/std_devs[numpy.newaxis].T

1470 1471 __all__.append('correlation_matrix') 1472 1473 ############################################################################### 1474 # Parsing of values with uncertainties: 1475 1476 POSITIVE_DECIMAL_UNSIGNED = r'(\d+)(\.\d*)?' 1477 1478 # Regexp for a number with uncertainty (e.g., "-1.234(2)e-6"), where the 1479 # uncertainty is optional (in which case the uncertainty is implicit): 1480 NUMBER_WITH_UNCERT_RE_STR = ''' 1481 ([+-])? # Sign 1482 %s # Main number 1483 (?:$%s$)? # Optional uncertainty 1484 ([eE][+-]?\d+)? # Optional exponent 1485 ''' % (POSITIVE_DECIMAL_UNSIGNED, POSITIVE_DECIMAL_UNSIGNED) 1486 1487 NUMBER_WITH_UNCERT_RE = re.compile( 1488 "^%s$" % NUMBER_WITH_UNCERT_RE_STR, re.VERBOSE)

1489 1490 -def parse_error_in_parentheses(representation):

1491 """ 1492 Returns (value, error) from a string representing a number with 1493 uncertainty like 12.34(5), 12.34(142), 12.5(3.4) or 12.3(4.2)e3. 1494 If no parenthesis is given, an uncertainty of one on the last 1495 digit is assumed. 1496 1497 Raises ValueError if the string cannot be parsed. 1498 """ 1499 1500 match = NUMBER_WITH_UNCERT_RE.search(representation) 1501 1502 if match: 1503 # The 'main' part is the nominal value, with 'int'eger part, and 1504 # 'dec'imal part. The 'uncert'ainty is similarly broken into its 1505 # integer and decimal parts. 1506 (sign, main_int, main_dec, uncert_int, uncert_dec, 1507 exponent) = match.groups() 1508 else: 1509 raise ValueError("Unparsable number representation: '%s'." 1510 " Was expecting a string of the form 1.23(4)" 1511 " or 1.234" % representation) 1512 1513 # The value of the number is its nominal value: 1514 value = float(''.join((sign or '', 1515 main_int, 1516 main_dec or '.0', 1517 exponent or ''))) 1518 1519 if uncert_int is None: 1520 # No uncertainty was found: an uncertainty of 1 on the last 1521 # digit is assumed: 1522 uncert_int = '1' 1523 1524 # Do we have a fully explicit uncertainty? 1525 if uncert_dec is not None: 1526 uncert = float("%s%s" % (uncert_int, uncert_dec or '')) 1527 else: 1528 # uncert_int represents an uncertainty on the last digits: 1529 1530 # The number of digits after the period defines the power of 1531 # 10 than must be applied to the provided uncertainty: 1532 num_digits_after_period = (0 if main_dec is None 1533 else len(main_dec)-1) 1534 uncert = int(uncert_int)/10**num_digits_after_period 1535 1536 # We apply the exponent to the uncertainty as well: 1537 uncert *= float("1%s" % (exponent or '')) 1538 1539 return (value, uncert)

1540

1541 1542 # The following function is not exposed because it can in effect be 1543 # obtained by doing x = ufloat(representation) and 1544 # x.nominal_value and x.std_dev(): 1545 -def str_to_number_with_uncert(representation):

1546 """ 1547 Given a string that represents a number with uncertainty, returns the 1548 nominal value and the uncertainty. 1549 1550 The string can be of the form: 1551 - 124.5+/-0.15 1552 - 124.50(15) 1553 - 124.50(123) 1554 - 124.5 1555 1556 When no numerical error is given, an uncertainty of 1 on the last 1557 digit is implied. 1558 1559 Raises ValueError if the string cannot be parsed. 1560 """ 1561 1562 try: 1563 # Simple form 1234.45+/-1.2: 1564 (value, uncert) = representation.split('+/-') 1565 except ValueError: 1566 # Form with parentheses or no uncertainty: 1567 parsed_value = parse_error_in_parentheses(representation) 1568 else: 1569 try: 1570 parsed_value = (float(value), float(uncert)) 1571 except ValueError: 1572 raise ValueError('Cannot parse %s: was expecting a number' 1573 ' like 1.23+/-0.1' % representation) 1574 1575 return parsed_value

1576

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

1578 """ 1579 Returns a random variable (Variable object). 1580 1581 Converts the representation of a number into a number with 1582 uncertainty (a random variable, defined by a nominal value and 1583 a standard deviation). 1584 1585 The representation can be a (value, standard deviation) sequence, 1586 or a string. 1587 1588 Strings of the form '12.345+/-0.015', '12.345(15)', or '12.3' are 1589 recognized (see full list below). In the last case, an 1590 uncertainty of +/-1 is assigned to the last digit. 1591 1592 'tag' is an optional string tag for the variable. Variables 1593 don't have to have distinct tags. Tags are useful for tracing 1594 what values (and errors) enter in a given result (through the 1595 error_components() method). 1596 1597 Examples of valid string representations: 1598 1599 -1.23(3.4) 1600 -1.34(5) 1601 1(6) 1602 3(4.2) 1603 -9(2) 1604 1234567(1.2) 1605 12.345(15) 1606 -12.3456(78)e-6 1607 12.3(0.4)e-5 1608 0.29 1609 31. 1610 -31. 1611 31 1612 -3.1e10 1613 169.0(7) 1614 169.1(15) 1615 """ 1616 1617 # This function is somewhat optimized so as to help with the 1618 # creation of lots of Variable objects (through unumpy.uarray, for 1619 # instance). 1620 1621 # representations is "normalized" so as to be a valid sequence of 1622 # 2 arguments for Variable(). 1623 1624 #! Accepting strings and any kind of sequence slows down the code 1625 # by about 5 %. On the other hand, massive initializations of 1626 # numbers with uncertainties are likely to be performed with 1627 # unumpy.uarray, which does not support parsing from strings and 1628 # thus does not have any overhead. 1629 1630 #! Different, in Python 3: 1631 if isinstance(representation, basestring): 1632 representation = str_to_number_with_uncert(representation) 1633 1634 #! The tag is forced to be a string, so that the user does not 1635 # create a Variable(2.5, 0.5) in order to represent 2.5 +/- 0.5. 1636 # Forcing 'tag' to be a string prevents numerical uncertainties 1637 # from being considered as tags, here: 1638 if tag is not None: 1639 #! 'unicode' is removed in Python3: 1640 assert isinstance(tag, (str, unicode)), "The tag can only be a string." 1641 1642 #! The special ** syntax is for Python 2.5 and before (Python 2.6+ 1643 # understands tag=tag): 1644 return Variable(*representation, **{'tag': tag})

1645

Source Code for Package ComboCode.cc.ivs.sigproc.lmfit.uncertainties