#========================= uncertainties.py ========================= """ Calculations with uncertainties (and correlations): class Number_with_uncert. Uncertainties are treated like standard deviations. Valid operations include basic mathematical functions (addition,...), as well as operations from the math module. Logical operations (>, <, etc.) are also supported, but the calculated error is generally meaningless. Applying operations on Number_with_uncert creates 'semi-formal' expressions (see below), represented by the Semi_formal_expr class. The method used for managing uncertainties consists in replacing numbers by a kind of mathematical expressions (of type Semi_formal_expr) that can be calculated for different values of the variables on which they depend--even if the dependence is not in the form of an explicit Python function call with some parameters. These expressions are 'semi-formal' in the sense that all their parameters are bound to a value at all times, but they can nonetheless be recalculated for different values of these parameters. In particular, these semi-formal expressions always have a value. Example: >>> x = Number_with_uncert(3.14, 0.01) >>> y = 2*x >>> print x-x 0.0 >>> print y-x # The error should be exactly the error on x 3.14 +- 0.010000000000 >>> x.nominal_value = 1 >>> print y # 'y' is updated 2.0 +- 0.0200000000001 (c) Eric O. LEBIGOT (EOL), 2009. Strongly inspired by code by Arnaud Delobelle (http://groups.google.com/group/comp.lang.python/msg/b92987c7787346ec) """ __all__ = ["Number_with_uncert"] ############################################################################### class _Exact_constant(object): """ Expression that represents an exact (no error) constant function. Used in wrap_func_output() for a uniform access to expression parts (wrap_func_output() expects some attributes from expression parts). """ # This class exists so as not to uselessly create semi-formal variables # (Semi_formal_var objects), which are kept in memory. def __init__(self, value): self._variables = set() self.nominal_value = value def to_func(x): """ Coerces x into a constant expression, unless x is already an expression (Semi_formal_expr object). """ return x if isinstance(x, Semi_formal_expr) else _Exact_constant(x) def wrap_func_output(f): """ Transforms a Python function into an expression that can return the result of 'f', but generally as a Semi_formal_expr object (unless its evaluation does not involve variables [Semi_formal_var objects], in which case 'f' simply returns its usual result). """ def f_with_expr_output(*args): #! The following does not seem to go into the expression returned # by wrap_func_output(). Why? """ Version of %s(...) that returns a semi-formal expression (Semi_formal_expr object), if its result depends on variables (Semi_formal_var objects). The new version returns a simple constant, when applied to constant arguments. Original documentation: %s """ % (f.__name__, f.__doc__) # Coercion of all arguments to semi-formal expressions: sub_exprs = map(to_func, args) # We keep track of all Semi_formal_var leaves: variables = set() for sub_expr in sub_exprs: variables |= sub_expr._variables # Formal version of 'f': its value can be calculated through # the specificaton of the values of the variables used in the # arguments 'args'. def f_evaluation(): """ Returns the value of function 'f', calculated at the nominal values of its arguments sub_exprs. """ return f(*(sub_expr.nominal_value for sub_expr in sub_exprs)) if variables: # Delayed evaluation: return Semi_formal_expr_node(f_evaluation, variables) else: # Constant functions do not have to be complicated # objects: return f_evaluation() return f_with_expr_output class Semi_formal_expr(object): """ Semi-formal expression object that supports the mathematical operations of Python floats and give objects of the same type (Semi_formal_expr). Semi_formal_expr objects can thus be summed, etc. They are ssentially identical to a mathematical expression involving floats, with the difference that it effectively contains a way of recalculating its value and adapt to changes in its variables. The only 'variables' considered in this class are Semi_formal_var objects (numbers with an uncertainty). Attributes and methods defined through inheritance: - (nominal_value, error): nominal (central) value and error on the expression. Before accessing these, no calculation of the expression is performed. These quantities are dynamically calculated. - derivative_value(variable): value of the partial derivative with respect to the given variable (Semi_formal_var object), at the point currently defined by the nominal values of the variables. - _variables: semi-formal variables (Semi_formal_var objects) on which the expression depends. """ def __repr__(self): return "%s object with result %s" % (type(self), str(self)) def __str__(self): (nominal_value, error) = self.nominal_value, self.error return ("%s +- %s" % (nominal_value, error) if error else str(nominal_value)) # Conversion to float would be risky: calculations could be # performed without error handling, and without the user noticing # it. A number with an uncertainty is not a pure number. Note # that float(1j) is not allowed, for instance. def __float__(self): raise TypeError("can't convert a number with uncertainty (%s)" " to float; use x.nominal_value" % self.__class__) # Operators with no reflection: # Logical operators: warning: the resulting value cannot always be # differentiated. for operator in ('eq', 'ge', 'gt', 'le', 'lt', 'ne'): exec ("__%s__ = wrap_func_output(float.__%s__)" % (operator, operator)) # __nonzero__() is supposed to return a boolean value (it is used # by bool()). It is for instance used for converting the result # of comparison operators to a boolean, in sorted(). If we want # to be able to sort Semi_formal_expr objects, __nonzero__ cannot # return a Semi_formal_expr object. Since boolean results (such # as the result of bool()) don't have a very meaningful # uncertainty, this should not be a big deal: def __nonzero__(self): return bool(self.nominal_value) # Operators that return a numerical value: for operator in ('abs', 'neg', 'pos'): exec ("__%s__ = wrap_func_output(float.__%s__)" % (operator, operator)) # Operators with a reflection: for operator in ('add', 'div', 'divmod', 'floordiv','mod', 'mul', 'pow', 'sub', 'truediv'): for prefix in ('', 'r'): method_name = "__%s%s__" % (prefix, operator) exec("%s = wrap_func_output(float.%s)" % (method_name, method_name)) del operator, prefix, method_name class Semi_formal_var(Semi_formal_expr): """ Semi-formal variable. The variable is bound to a value+uncertainty at all times. This is a special kind of semi-formal expression (Semi_formal_expr object), with a constant value. """ def __init__(self, nominal_value, error = 0): """ 'nominal_value' is the nominal value of the semi-formal variable. 'error' is the error on the value. """ # We initialize the value in the same way as users of # instances who would set it: users use the 'nominal_value' attribute: self.nominal_value = nominal_value # Defines the (constant) expression self.error = error # Expression 'x' depends on 'x': self._variables = set([self]) def get_value(self): return self.__nominal_value def set_value(self, nominal_value): """ Since evaluations of Semi_formal_expr expression use functions that use float arguments, the value of the variables needs to be a float. """ self.__nominal_value = float(nominal_value) nominal_value = property(get_value, set_value) def derivative_value(self, variable, step = 1): return 1. if self is variable else 0. class Semi_formal_expr_node(Semi_formal_expr): """ Semi-formal expression, as defined by a 'regular' Python function. """ def __init__(self, func, variables): """ Node for a semi-formal expression (Semi_formal_expr). For the meaning of object attributes see the documentation for Semi_formal_expr. """ self._evaluate = func # Way of evaluating the function self._variables = variables #! It might be numerically more relevant to evaluate at .nominal_value+.err # ... but this is not standard. AND: I'm not sure this would yield # a more correct way of estimating the variance (does this take the # second order derivative into account??) THIS WOULD BE USEFUL # if .err is a more reasonnable value than step... def derivative_value(self, variable, step = 1e-5): """ Calculates the derivative of the expression with respect to the given variable (Semi_formal_var object), at the point defined by the nominal (central) value of its variables. 'step' is the numerical variation on 'variable' for the numerical calculation of the derivative. -log10(step) is an estimate of the number of digits lost in the derivative calculation. """ # Value of the function: central_value = self._evaluate() # We temporarily shift the value of the variable: previous_nominal_value = variable.nominal_value variable.nominal_value += step new_value = self._evaluate() variable.nominal_value = previous_nominal_value derivative_value = (new_value-central_value)/step return derivative_value @property def nominal_value(self): """ Returns the nominal (central) value of the expression, for the current nominal value of its variables. """ return self._evaluate() @property def error(self): """ Returns the error on the expression, at the current nominal value of its variables. """ # Caculation of the variance: variance = 0 for variable in self._variables: #! This is not efficient: the "central" function value is # calculated many times: (an alternative would be to # effectively calculate the derivative here) value_shift = (self.derivative_value(variable)*variable.error) variance += value_shift**2 error = variance**0.5 return error # We wrap the expressions from the math module so that they keep track # of uncertainties. BUILTIN_FUNC = type(sum) # Built-in expression type import math for name in dir(math): obj = getattr(math, name) if isinstance(obj, BUILTIN_FUNC): setattr(math, name, wrap_func_output(obj)) ############################################################################### # Values with uncertainties: # Number_with_uncert = Semi_formal_var ###############################################################################