ComboCode.cc.ivs.sigproc.lmfit.uncertainties.unumpy.core

1 """ 2 Core functions used by unumpy and some of its submodules. 3 4 (c) 2010-2013 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 from numpy.core import numeric 20 21 # Local modules: 22 from cc.ivs.sigproc.lmfit import uncertainties 23 from cc.ivs.sigproc.lmfit.uncertainties import umath 24 25 from cc.ivs.sigproc.lmfit.uncertainties import __author__ 26 27 __all__ = [ 28 # Factory functions: 29 'uarray', 'umatrix', 30 31 # Utilities: 32 'nominal_values', 'std_devs', 33 34 # Classes: 35 'matrix' 36 ] 37 38 ############################################################################### 39 # Utilities: 40 41 # nominal_values() and std_devs() are defined as functions (instead of 42 # as additional methods of the unumpy.matrix class) because the user 43 # might well directly build arrays of numbers with uncertainties 44 # without going through the factory functions found in this module 45 # (uarray() and umatrix()). Thus, 46 # numpy.array([uncertainties.ufloat((1, 0.1))]) would not 47 # have a nominal_values() method. Adding such a method to, say, 48 # unumpy.matrix, would break the symmetry between NumPy arrays and 49 # matrices (no nominal_values() method), and objects defined in this 50 # module. 51 52 # ! Warning: the __doc__ is set, but help(nominal_values) does not 53 # display it, but instead displays the documentation for the type of 54 # nominal_values (i.e. the documentation of its class): 55 56 to_nominal_values = numpy.vectorize( 57 uncertainties.nominal_value, 58 otypes=[float], # Because vectorize() has side effects (dtype setting) 59 doc=("Applies uncertainties.nominal_value to the elements of" 60 " a NumPy (or unumpy) array (this includes matrices).")) 61 62 to_std_devs = numpy.vectorize( 63 uncertainties.std_dev, 64 otypes=[float], # Because vectorize() has side effects (dtype setting) 65 doc=("Returns the standard deviation of the numbers with uncertainties" 66 " contained in a NumPy array, or zero for other objects."))

67 68 -def unumpy_to_numpy_matrix(arr):

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

75

76 -def nominal_values(arr):

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

90

91 -def std_devs(arr):

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

105

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

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

121

122 -def wrap_array_func(func):

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

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

244 245 -def uarray((values, std_devs)):

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

257

258 ############################################################################### 259 260 -def array_derivative(array_like, var):

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

279

280 -def func_with_deriv_to_uncert_func(func_with_derivatives):

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

387 388 return wrapped_func 389

390 ########## Matrix inverse 391 392 -def inv_with_derivatives(arr, input_type, derivatives):

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

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

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

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

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

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

501 return _pinv_with_uncert(array_like, rcond)

502

503 ########## Matrix class 504 505 -class matrix(numpy.matrix):

506 # The name of this class is the same as NumPy's, which is why it 507 # does not follow PEP 8. 508 """ 509 Class equivalent to numpy.matrix, but that behaves better when the 510 matrix contains numbers with uncertainties. 511 """ 512

513 - def __rmul__(self, other):

514 # ! NumPy's matrix __rmul__ uses an apparently a restrictive 515 # dot() function that cannot handle the multiplication of a 516 # scalar and of a matrix containing objects (when the 517 # arguments are given in this order). We go around this 518 # limitation: 519 if numeric.isscalar(other): 520 return numeric.dot(self, other) 521 else: 522 return numeric.dot(other, self) # The order is important

523 524 # The NumPy doc for getI is empty: 525 # @uncertainties.set_doc(numpy.matrix.getI.__doc__)

526 - def getI(self):

527 "Matrix inverse of pseudo-inverse" 528 529 # numpy.matrix.getI is OK too, but the rest of the code assumes that 530 # numpy.matrix.I is a property object anyway: 531 532 M, N = self.shape 533 if M == N: 534 func = _inv 535 else: 536 func = _pinv 537 return func(self)

538 539 540 # ! In Python >= 2.6, this could be simplified as: 541 # I = numpy.matrix.I.getter(__matrix_inverse) 542 I = property(getI, numpy.matrix.I.fset, numpy.matrix.I.fdel, 543 numpy.matrix.I.__doc__) 544 545 @property

546 - def nominal_values(self):

547 """ 548 Nominal value of all the elements of the matrix. 549 """ 550 return nominal_values(self)

551 552 std_devs = std_devs

553

554 -def umatrix(*args):

555 """ 556 Constructs a matrix that contains numbers with uncertainties. 557 558 The input data is the same as for uarray(...): a tuple with the 559 nominal values, and the standard deviations. 560 561 The returned matrix can be inverted, thanks to the fact that it is 562 a unumpy.matrix object instead of a numpy.matrix one. 563 """ 564 565 return uarray(*args).view(matrix)

566

567 ############################################################################### 568 569 -def define_vectorized_funcs():

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

611 612 define_vectorized_funcs() 613

Source Code for Module ComboCode.cc.ivs.sigproc.lmfit.uncertainties.unumpy.core