eric.le.bi...@spectro.jussieu.fr writes:
> Arnaud, your code is very interesting!
>
> On Apr 15, 1:00 pm, Arnaud Delobelle <arno...@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
--
http://mail.python.org/mailman/listinfo/python-list
