Modifying the value of a float-like object

[Arnaud Delobelle]

Eric.Le.Bigot at spectro.jussieu.fr writes:

> Arnaud, your code is very interesting!
>
> On Apr 15, 1:00 pm, Arnaud Delobelle <arno... at googlemail.com> wrote:
>> I still don't understand why you need mutable floats.
>
> Here is why: the code that your proposed (or any code that does
> straightforward error propagation, for that matter) does not generally
> calculate uncertainties correctly:
>
>>>> a=Value(1, 0.1)
>>>> a-a
> 0.0 +- 0.14142135623823598
>
> The correct result is obviously 0.0 +- 0.0.  This is the effect of
> what I referred to in previous posts as "correlated errors".
>
> Anyway, I learned some interesting stuff about what you can do in
> Python, thanks to your code! :)

I still don't think mutable floats are necessary.  Here is an approach
below - I'll let the code speak because I have to do some shopping!

It still relies on Peter Otten's method for error calculation - which I
trust is good as my numerical analysis is to say the least very rusty!

---------- uncert2.py ----------
def uexpr(x):
    return x if isinstance(x, UBase) else UVal(x)

def uified(f):
    def decorated(*args):
        args = map(uexpr, args)
        basis = set()
        for arg in args:
            basis |= arg.basis
        uf = lambda values: f(*(x.f(values) for x in args))
        ret = UExpr(basis, uf)
        return ret if ret.err else ret.value
    return decorated

class UBase(object):
    def __repr__(self): 
        return "%r +- %r" % (self.value, self.err) 
    def __hash__(self):
        return id(self)
    __neg__ = uified(float.__neg__)
    for op in 'add', 'sub', 'mul', 'div', 'pow':
        for r in '', 'r':
            exec """__%s__ = uified(float.__%s__)""" % (r+op, r+op)
    del op, r

class UVal(UBase): 
    def __init__(self, value, err=0):
        if isinstance(value, UVal):
            value, err = value.value, value.err
        self.value = value 
        self.err = err 
        self.basis = set([self])
    def f(self, values):
        return values[self]


class UExpr(UBase):
    def __init__(self, basis, f):
        self.basis = basis
        self.f = f
        self.calc()
    def derive(self, i, eps=1e-5):
        values = dict((x, x.value) for x in self.basis)
        values[i] += eps 
        x2 = self.f(values) 
        return (x2-self.value)/eps     
    def calc(self): 
        sigma = 0
        values = dict((x, x.value) for x in self.basis)
        self.value = self.f(values)
        for i in self.basis:
            x = self.derive(i)*i.err 
            sigma += x*x 
        self.err = sigma**0.5


builtinf = type(sum)
import math
for name in dir(math):
    obj = getattr(math, name)
    if isinstance(obj, builtinf):
        setattr(math, name,  uified(obj))
----------------------------------------

Example:

marigold:junk arno$ python -i uncert2.py 
>>> a = UVal(2.0, 0.1)
>>> b = UVal(10.0, 1)
>>> a + a
4.0 +- 0.20000000000131024
>>> a*b
20.0 +- 2.2360679774310253
>>> a - a
0.0
>>> a / a
1.0
>>> from math import *
>>> pow(b, a)
100.0 +- 30.499219977998791
>>> sin(a) - tan(b)
0.26093659936659497 +- 1.4209904731243463
>>> sin(a)*sin(a) + cos(a)*cos(a)
1.0
>>> sin(a)/cos(a) - tan(a)
0.0
>>> a*b - b*a
0.0
>>> # Etc...

-- 
Arnaud