Re: [boost] units review

From: Kevin Lynch (krlynch_at_[hidden])
Date: 2007-04-05 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 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:	[phone]
Physics Department			Fax: [phone]
Boston University			office:	 PRB-361
590 Commonwealth Ave.			e-mail:	 krlynch_at_[hidden]
Boston, MA 02215 USA			http://budoe.bu.edu/~krlynch
-------------------------------------------------------------------------------