Welcome, guest | Sign In | My Account | Store | Cart
```#========================= 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
"""

__all__ = ["Number_with_uncert"]

###############################################################################

class _Exact_constant(object):
"""
Expression that represents an exact (no error) constant function.

(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
###############################################################################
```