
Boost : 
From: Kevin Lynch (krlynch_at_[hidden])
Date: 20070405 14:52:22
Matthias Schabel wrote:
> I'm actually thinking of including measurement<> with the library
> itself as a way to provide both physical constants and their
> associated measurement errors (from CODATA, in particular), so your
> comments are very timely. I was planning on tweaking the
> implementation, but I had naively failed to consider the issue of
> correlated errors... I'll have to see what's possible with a
> reasonable amount of effort. I certainly don't want to try to come up
> with a general solution for dealing with error covariances  any such
> solution would be messy and complex  but the self correlation case
> is clearly important... If you have any other ideas/input, it would
> be most welcome.
>
I regret that I don't have any good ideas. Solving this is a problem
near and dear to me ... in my early forays into C++, I wrote the
equivalent of you measurement<> class, spending way too much time on it
before I realized that it certinly wasn't going to do the right thing.
The experiment I currently expend most of my effort on has an ultimate
goal of measuring the positive muon lifetime at 1ppm (hence my love for
the Fermi constant :) Keeping very careful track of errors and
correctly propagating them through calculations is critical for us.
Long story short, I've spent a lot of time thinking very hard about how
one might go about doing so "correctly" in pure C++, and I haven't found
a way to do it (note, however, that I have no real experience or
facility in the template metaprogramming ... most of my scientific
programming effort has gone in to floating point and numerical issues,
and interacting with DAQ hardware). Even ignoring the parameter
covariances (by assuming all variable errors are independent, or
equivalently, a diagonal covariance matrix), I think it's a very very
hard problem. The very simple example I posted earlier:
measurement<> x(1.,1.);
assert( 2*x == x+x );
will assert, but that's just the tip of the iceberg. A robust, correct
implementation would need to be able to track which measurement instance
it was looking at _across function calls_ and perhaps even across
translation units. None of the following should assert ... all of them
will, with either your implementation or mine:
measurement<> x(1.,1.);
measurement<> f(y+z) { return y+z; }
measurement<> y = f(x,x) + f(x,x); // not independent
assert ( y.e = 4*x.e );
measurement<> z = x; // not independent
measurement<> w = z+x; // not independent
assert ( w.e == 2*x.e );
measurement<> a(1.,1.);
extern f(measurement<>& b);
measurement<> c = a+f(a); // what happens here?
I think you would need to write the equivalent of a computer algebra
system to track what variables are "independent", and you would have to
do so across programming boundaries in a nontrivial way.
I think this is a problem akin to global compiler optimization. I'm not
sure there's even a way to do it in C++ without significant
restrictions. The only way I know to deal with this problem is to use a
computer algebra package designed to support this need, and even then,
in the packages I know, you have to supply all calculations in the
equivalent of a single translation unit. For instance, there's a
Mathematica package (the name escapes me) and open source tools like
"fussy", which I use extensively.
I'd love for someone to figure out how to do it right, but I'm not going
to hold my breath....
  Kevin Lynch voice: (617) 3536025 Physics Department Fax: (617) 3539393 Boston University office: PRB361 590 Commonwealth Ave. email: krlynch_at_[hidden] Boston, MA 02215 USA http://budoe.bu.edu/~krlynch 
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk