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 non-trivial 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) 353-6025 Physics Department Fax: (617) 353-9393 Boston University office: PRB-361 590 Commonwealth Ave. e-mail: krlynch@bu.edu Boston, MA 02215 USA http://budoe.bu.edu/~krlynch -------------------------------------------------------------------------------