ComboCode.cc.ivs.units.uncertainties.unumpy.core

1 """ 2 Core functions used by unumpy and some of its submodules. 3 4 (c) 2010 by Eric O. LEBIGOT (EOL). 5 """ 6 7 # The functions found in this module cannot be defined in unumpy or 8 # its submodule: this creates import loops, when unumpy explicitly 9 # imports one of the submodules in order to make it available to the 10 # user. 11 12 from __future__ import division 13 14 # Standard modules: 15 import sys 16 17 # 3rd-party modules: 18 import numpy 19 20 # Local modules: 21 from cc.ivs.units import uncertainties 22 from cc.ivs.units.uncertainties import umath 23 24 from cc.ivs.units.uncertainties import __author__ 25 26 __all__ = [ 27 # Factory functions: 28 'uarray', 'umatrix', 29 30 # Utilties: 31 'nominal_values', 'std_devs', 32 33 # Classes: 34 'matrix' 35 ] 36 37 ############################################################################### 38 # Utilities: 39 40 # nominal_values() and std_devs() are defined as functions (instead of 41 # as additional methods of the unumpy.matrix class) because the user 42 # might well directly build arrays of numbers with uncertainties 43 # without going through the factory functions found in this module 44 # (uarray() and umatrix()). Thus, 45 # numpy.array([uncertainties.ufloat((1, 0.1))]) would not 46 # have a nominal_values() method. Adding such a method to, say, 47 # unumpy.matrix, would break the symmetry between NumPy arrays and 48 # matrices (no nominal_values() method), and objects defined in this 49 # module. 50 51 # ! Warning: the __doc__ is set, but help(nominal_values) does not 52 # display it, but instead displays the documentation for the type of 53 # nominal_values (i.e. the documentation of its class): 54 55 to_nominal_values = numpy.vectorize( 56 uncertainties.nominal_value, 57 otypes=[float], # Because vectorize() has side effects (dtype setting) 58 doc=("Applies uncertainties.nominal_value to the elements of" 59 " a NumPy (or unumpy) array (this includes matrices).")) 60 61 to_std_devs = numpy.vectorize( 62 uncertainties.std_dev, 63 otypes=[float], # Because vectorize() has side effects (dtype setting) 64 doc=("Returns the standard deviation of the numbers with uncertainties" 65 " contained in a NumPy array, or zero for other objects."))

66 67 -def unumpy_to_numpy_matrix(arr):

68 """ 69 If arr in a unumpy.matrix, it is converted to a numpy.matrix. 70 Otherwise, it is returned unchanged. 71 """ 72 73 return arr.view(numpy.matrix) if isinstance(arr, matrix) else arr

74

75 -def nominal_values(arr):

76 """ 77 Returns the nominal values of the numbers in NumPy array arr. 78 79 Elements that are not uncertainties.AffineScalarFunc are passed 80 through untouched (because a numpy.array can contain numbers with 81 uncertainties and pure floats simultaneously). 82 83 If arr is of type unumpy.matrix, the returned array is a 84 numpy.matrix, because the resulting matrix does not contain 85 numbers with uncertainties. 86 """ 87 88 return unumpy_to_numpy_matrix(to_nominal_values(arr))

89

90 -def std_devs(arr):

91 """ 92 Returns the standard deviations of the numbers in NumPy array arr. 93 94 Elements that are not uncertainties.AffineScalarFunc are given a 95 zero uncertainty ((because a numpy.array can contain numbers with 96 uncertainties and pure floats simultaneously).. 97 98 If arr is of type unumpy.matrix, the returned array is a 99 numpy.matrix, because the resulting matrix does not contain 100 numbers with uncertainties. 101 """ 102 103 return unumpy_to_numpy_matrix(to_std_devs(arr))

104

105 ############################################################################### 106 107 -def derivative(u, var):

108 """ 109 Returns the derivative of u along var, if u is an 110 uncertainties.AffineScalarFunc instance, and if var is one of the 111 variables on which it depends. Otherwise, return 0. 112 """ 113 if isinstance(u, uncertainties.AffineScalarFunc): 114 try: 115 return u.derivatives[var] 116 except KeyError: 117 return 0. 118 else: 119 return 0.

120

121 -def wrap_array_func(func):

122 """ 123 Returns a version of the function func() that works even when 124 func() is given a NumPy array that contains numbers with 125 uncertainties. 126 127 func() is supposed to return a NumPy array. 128 129 This wrapper is similar to uncertainties.wrap(), except that it 130 handles an array argument instead of float arguments. 131 132 func -- version that takes and returns a single NumPy array. 133 """ 134 135 @uncertainties.set_doc("""\ 136 Version of %s(...) that works even when its first argument is a NumPy 137 array that contains numbers with uncertainties. 138 139 Warning: elements of the first argument array that are not 140 AffineScalarFunc objects must not depend on uncertainties.Variable 141 objects in any way. Otherwise, the dependence of the result in 142 uncertainties.Variable objects will be incorrect. 143 144 Original documentation: 145 %s""" % (func.__name__, func.__doc__)) 146 def wrapped_func(arr, *args): 147 # Nominal value: 148 arr_nominal_value = nominal_values(arr) 149 func_nominal_value = func(arr_nominal_value, *args) 150 151 # The algorithm consists in numerically calculating the derivatives 152 # of func: 153 154 # Variables on which the array depends are collected: 155 variables = set() 156 for element in arr.flat: 157 # floats, etc. might be present 158 if isinstance(element, uncertainties.AffineScalarFunc): 159 variables |= set(element.derivatives.iterkeys()) 160 161 # If the matrix has no variables, then the function value can be 162 # directly returned: 163 if not variables: 164 return func_nominal_value 165 166 # Calculation of the derivatives of each element with respect 167 # to the variables. Each element must be independent of the 168 # others. The derivatives have the same shape as the output 169 # array (which might differ from the shape of the input array, 170 # in the case of the pseudo-inverse). 171 derivatives = numpy.vectorize(lambda _: {})(func_nominal_value) 172 for var in variables: 173 174 # A basic assumption of this package is that the user 175 # guarantees that uncertainties cover a zone where 176 # evaluated functions are linear enough. Thus, numerical 177 # estimates of the derivative should be good over the 178 # standard deviation interval. This is true for the 179 # common case of a non-zero standard deviation of var. If 180 # the standard deviation of var is zero, then var has no 181 # impact on the uncertainty of the function func being 182 # calculated: an incorrect derivative has no impact. One 183 # scenario can give incorrect results, however, but it 184 # should be extremely uncommon: the user defines a 185 # variable x with 0 standard deviation, sets y = func(x) 186 # through this routine, changes the standard deviation of 187 # x, and prints y; in this case, the uncertainty on y 188 # might be incorrect, because this program had no idea of 189 # the scale on which func() is linear, when it calculated 190 # the numerical derivative. 191 192 # The standard deviation might be numerically too small 193 # for the evaluation of the derivative, though: we set the 194 # minimum variable shift. 195 196 shift_var = max(var._std_dev/1e5, 1e-8*abs(var._nominal_value)) 197 # An exceptional case is that of var being exactly zero. 198 # In this case, an arbitrary shift is used for the 199 # numerical calculation of the derivative. The resulting 200 # derivative value might be quite incorrect, but this does 201 # not matter as long as the uncertainty of var remains 0, 202 # since it is, in this case, a constant. 203 if not shift_var: 204 shift_var = 1e-8 205 206 # Shift of all the elements of arr when var changes by shift_var: 207 shift_arr = array_derivative(arr, var)*shift_var 208 209 # Origin value of array arr when var is shifted by shift_var: 210 shifted_arr_values = arr_nominal_value + shift_arr 211 func_shifted = func(shifted_arr_values, *args) 212 numerical_deriv = (func_shifted-func_nominal_value)/shift_var 213 214 # Update of the list of variables and associated 215 # derivatives, for each element: 216 for (derivative_dict, derivative_value) \ 217 in zip(derivatives.flat, numerical_deriv.flat): 218 if derivative_value: 219 derivative_dict[var] = derivative_value 220 221 # numbers with uncertainties are build from the result: 222 return numpy.vectorize(uncertainties.AffineScalarFunc)( 223 func_nominal_value, derivatives)

224 225 # It is easier to work with wrapped_func, which represents a 226 # wrapped version of 'func', when it bears the same name as 227 # 'func' (the name is used by repr(wrapped_func)). 228 wrapped_func.__name__ = func.__name__ 229 230 return wrapped_func 231 232 ############################################################################### 233 # Arrays 234 235 # Vectorized creation of an array of variables: 236 237 # ! Looking up uncertainties.Variable beforehand through '_Variable = 238 # uncertainties.Variable' does not result in a significant speed up: 239 240 _uarray = numpy.vectorize(lambda v, s: uncertainties.Variable(v, s), 241 otypes=[object])

242 243 -def uarray((values, std_devs)):

244 """ 245 Returns a NumPy array of numbers with uncertainties 246 initialized with the given nominal values and standard 247 deviations. 248 249 values, std_devs -- valid arguments for numpy.array, with 250 identical shapes (list of numbers, list of lists, numpy.ndarray, 251 etc.). 252 """ 253 254 return _uarray(values, std_devs)

255

256 ############################################################################### 257 258 -def array_derivative(array_like, var):

259 """ 260 Returns the derivative of the given array with respect to the 261 given variable. 262 263 The returned derivative is a Numpy ndarray of the same shape as 264 array_like, that contains floats. 265 266 array_like -- array-like object (list, etc.) that contains 267 scalars or numbers with uncertainties. 268 269 var -- Variable object. 270 """ 271 return numpy.vectorize(lambda u: derivative(u, var), 272 # The type is set because an 273 # integer derivative should not 274 # set the output type of the 275 # array: 276 otypes=[float])(array_like)

277

278 -def func_with_deriv_to_uncert_func(func_with_derivatives):

279 """ 280 Returns a function that can be applied to array-like objects that 281 contain numbers with uncertainties (lists, lists of lists, Numpy 282 arrays, etc.). 283 284 func_with_derivatives -- defines a function that takes array-like 285 objects containing scalars and returns an array. Both the value 286 and the derivatives of this function with respect to multiple 287 scalar parameters are calculated by func_with_derivatives(). 288 289 func_with_derivatives(arr, input_type, derivatives, *args) returns 290 an iterator. The first element is the value of the function at 291 point 'arr' (with the correct type). The following elements are 292 arrays that represent the derivative of the function for each 293 derivative array from the iterator 'derivatives'. 294 295 func_with_derivatives takes the following arguments: 296 297 arr -- Numpy ndarray of scalars where the function must be 298 evaluated. 299 300 input_type -- type of the input array-like object. This type is 301 used for determining the type that the function should return. 302 303 derivatives -- iterator that returns the derivatives of the 304 argument of the function with respect to multiple scalar 305 variables. func_with_derivatives() returns the derivatives of 306 the defined function with respect to these variables. 307 308 args -- additional arguments that define the result (example: 309 for the pseudo-inverse numpy.linalg.pinv: numerical cutoff). 310 311 Examples of func_with_derivatives: inv_with_derivatives(). 312 """ 313 314 def wrapped_func(array_like, *args): 315 """ 316 array_like -- array-like object that contains numbers with 317 uncertainties (list, Numpy ndarray or matrix, etc.). 318 319 args -- additional arguments that are passed directly to 320 func_with_derivatives. 321 """ 322 323 # So that .flat works even if array_like is a list. Later 324 # useful for faster code: 325 array_version = numpy.asarray(array_like) 326 327 # Variables on which the array depends are collected: 328 variables = set() 329 for element in array_version.flat: 330 # floats, etc. might be present 331 if isinstance(element, uncertainties.AffineScalarFunc): 332 variables |= set(element.derivatives.iterkeys()) 333 334 array_nominal = nominal_values(array_version) 335 # Function value, and derivatives at array_nominal (the 336 # derivatives are with respect to the variables contained in 337 # array_like): 338 func_and_derivs = func_with_derivatives( 339 array_nominal, 340 type(array_like), 341 (array_derivative(array_version, var) for var in variables), 342 *args) 343 344 func_nominal_value = func_and_derivs.next() 345 346 if not variables: 347 return func_nominal_value 348 349 # The result is built progressively, with the contribution of 350 # each variable added in turn: 351 352 # Calculation of the derivatives of the result with respect to 353 # the variables. 354 derivatives = numpy.array( 355 [{} for _ in xrange(func_nominal_value.size)], dtype=object) 356 derivatives.resize(func_nominal_value.shape) 357 358 # Memory-efficient approach. A memory-hungry approach would 359 # be to calculate the matrix derivatives will respect to all 360 # variables and then combine them into a matrix of 361 # AffineScalarFunc objects. The approach followed here is to 362 # progressively build the matrix of derivatives, by 363 # progressively adding the derivatives with respect to 364 # successive variables. 365 for (var, deriv_wrt_var) in zip(variables, func_and_derivs): 366 367 # Update of the list of variables and associated 368 # derivatives, for each element: 369 for (derivative_dict, derivative_value) in zip( 370 derivatives.flat, deriv_wrt_var.flat): 371 if derivative_value: 372 derivative_dict[var] = derivative_value 373 374 # An array of numbers with uncertainties are built from the 375 # result: 376 result = numpy.vectorize(uncertainties.AffineScalarFunc)( 377 func_nominal_value, derivatives) 378 379 # Numpy matrices that contain numbers with uncertainties are 380 # better as unumpy matrices: 381 if isinstance(result, numpy.matrix): 382 result = result.view(matrix) 383 384 return result

385 386 return wrapped_func 387

388 ########## Matrix inverse 389 390 -def inv_with_derivatives(arr, input_type, derivatives):

391 """ 392 Defines the matrix inverse and its derivatives. 393 394 See the definition of func_with_deriv_to_uncert_func() for its 395 detailed semantics. 396 """ 397 398 inverse = numpy.linalg.inv(arr) 399 # The inverse of a numpy.matrix is a numpy.matrix. It is assumed 400 # that numpy.linalg.inv is such that other types yield 401 # numpy.ndarrays: 402 if issubclass(input_type, numpy.matrix): 403 inverse = inverse.view(numpy.matrix) 404 yield inverse 405 406 # It is mathematically convenient to work with matrices: 407 inverse_mat = numpy.asmatrix(inverse) 408 409 # Successive derivatives of the inverse: 410 for derivative in derivatives: 411 derivative_mat = numpy.asmatrix(derivative) 412 yield -inverse_mat * derivative_mat * inverse_mat

413 414 _inv = func_with_deriv_to_uncert_func(inv_with_derivatives) 415 _inv.__doc__ = """\ 416 Version of numpy.linalg.inv that works with array-like objects 417 that contain numbers with uncertainties. 418 419 The result is a unumpy.matrix if numpy.linalg.pinv would return a 420 matrix for the array of nominal values. 421 422 Analytical formulas are used. 423 424 Original documentation: 425 %s 426 """ % numpy.linalg.inv.__doc__

427 428 ########## Matrix pseudo-inverse 429 430 -def pinv_with_derivatives(arr, input_type, derivatives, rcond):

431 """ 432 Defines the matrix pseudo-inverse and its derivatives. 433 434 Works with real or complex matrices. 435 436 See the definition of func_with_deriv_to_uncert_func() for its 437 detailed semantics. 438 """ 439 440 inverse = numpy.linalg.pinv(arr, rcond) 441 # The pseudo-inverse of a numpy.matrix is a numpy.matrix. It is 442 # assumed that numpy.linalg.pinv is such that other types yield 443 # numpy.ndarrays: 444 if issubclass(input_type, numpy.matrix): 445 inverse = inverse.view(numpy.matrix) 446 yield inverse 447 448 # It is mathematically convenient to work with matrices: 449 inverse_mat = numpy.asmatrix(inverse) 450 451 # Formula (4.12) from The Differentiation of Pseudo-Inverses and 452 # Nonlinear Least Squares Problems Whose Variables 453 # Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM 454 # Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), 455 # pp. 413-432 456 457 # See also 458 # http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse 459 460 # Shortcuts. All the following factors should be numpy.matrix objects: 461 PA = arr*inverse_mat 462 AP = inverse_mat*arr 463 factor21 = inverse_mat*inverse_mat.H 464 factor22 = numpy.eye(arr.shape[0])-PA 465 factor31 = numpy.eye(arr.shape[1])-AP 466 factor32 = inverse_mat.H*inverse_mat 467 468 # Successive derivatives of the inverse: 469 for derivative in derivatives: 470 derivative_mat = numpy.asmatrix(derivative) 471 term1 = -inverse_mat*derivative_mat*inverse_mat 472 derivative_mat_H = derivative_mat.H 473 term2 = factor21*derivative_mat_H*factor22 474 term3 = factor31*derivative_mat_H*factor32 475 yield term1+term2+term3

476 477 # Default rcond argument for the generalization of numpy.linalg.pinv: 478 try: 479 # Python 2.6+: 480 _pinv_default = numpy.linalg.pinv.__defaults__[0] 481 except AttributeError: 482 _pinv_default = 1e-15 483 484 _pinv_with_uncert = func_with_deriv_to_uncert_func(pinv_with_derivatives) 485 486 @uncertainties.set_doc(""" 487 Version of numpy.linalg.pinv that works with array-like objects 488 that contain numbers with uncertainties. 489 490 The result is a unumpy.matrix if numpy.linalg.pinv would return a 491 matrix for the array of nominal values. 492 493 Analytical formulas are used. 494 495 Original documentation: 496 %s 497 """ % numpy.linalg.pinv.__doc__)

498 -def _pinv(array_like, rcond=_pinv_default):

499 return _pinv_with_uncert(array_like, rcond)

500

501 ########## Matrix class 502 503 -class matrix(numpy.matrix):

504 # The name of this class is the same as NumPy's, which is why it 505 # does not follow PEP 8. 506 """ 507 Class equivalent to numpy.matrix, but that behaves better when the 508 matrix contains numbers with uncertainties. 509 """ 510 511 # The NumPy doc for getI is empty: 512 # @uncertainties.set_doc(numpy.matrix.getI.__doc__)

513 - def getI(self):

514 "Matrix inverse of pseudo-inverse" 515 516 # numpy.matrix.getI is OK too, but the rest of the code assumes that 517 # numpy.matrix.I is a property object anyway: 518 519 M, N = self.shape 520 if M == N: 521 func = _inv 522 else: 523 func = _pinv 524 return func(self)

525 526 527 # ! In Python >= 2.6, this could be simplified as: 528 # I = numpy.matrix.I.getter(__matrix_inverse) 529 I = property(getI, numpy.matrix.I.fset, numpy.matrix.I.fdel, 530 numpy.matrix.I.__doc__) 531 532 @property

533 - def nominal_values(self):

534 """ 535 Nominal value of all the elements of the matrix. 536 """ 537 return nominal_values(self)

538 539 std_devs = std_devs

540

541 -def umatrix(*args):

542 """ 543 Constructs a matrix that contains numbers with uncertainties. 544 545 The input data is the same as for uarray(...): a tuple with the 546 nominal values, and the standard deviations. 547 548 The returned matrix can be inverted, thanks to the fact that it is 549 a unumpy.matrix object instead of a numpy.matrix one. 550 """ 551 552 return uarray(*args).view(matrix)

553

554 ############################################################################### 555 556 -def define_vectorized_funcs():

557 """ 558 Defines vectorized versions of functions from uncertainties.umath. 559 560 Some functions have their name translated, so as to follow NumPy's 561 convention (example: math.acos -> numpy.arccos). 562 """ 563 564 this_module = sys.modules[__name__] 565 # NumPy does not always use the same function names as the math 566 # module: 567 func_name_translations = dict( 568 (f_name, 'arc'+f_name[1:]) 569 for f_name in ['acos', 'acosh', 'asin', 'atan', 'atan2', 'atanh']) 570 571 new_func_names = [ 572 func_name_translations[function_name] 573 if function_name in func_name_translations 574 else function_name 575 for function_name in umath.many_scalar_to_scalar_funcs] 576 577 for (function_name, unumpy_name) in \ 578 zip(umath.many_scalar_to_scalar_funcs, new_func_names): 579 580 # ! The newly defined functions (uncertainties.unumpy.cos, etc.) 581 # do not behave exactly like their NumPy equivalent (numpy.cos, 582 # etc.): cos(0) gives an array() and not a 583 # numpy.float... (equality tests succeed, though). 584 func = getattr(umath, function_name) 585 setattr( 586 this_module, unumpy_name, 587 numpy.vectorize(func, 588 # If by any chance a function returns, 589 # in a particular case, an integer, 590 # side-effects in vectorize() would 591 # fix the resulting dtype to integer, 592 # which is not what is wanted: 593 otypes=[object], 594 doc="""\ 595 Vectorized version of umath.%s. 596 597 Original documentation: 598 %s""" % (function_name, func.__doc__))) 599 600 __all__.append(unumpy_name)

601 602 define_vectorized_funcs() 603

Source Code for Module ComboCode.cc.ivs.units.uncertainties.unumpy.core