[GRASS5] Re: [GRASS_DE] r.mapcalc round()

[Glynn Clements]

Markus Neteler wrote:

> > > > Does this go away if you compile r.univar with "-ffloat-store"?
> > > 
> > > Well, r.univar is only a script using r.stats and awk. However, after
> > > recompiling r.stats with above flag (adding 
> > > EXTRA_CFLAGS=-ffloat-store
> > > to the r.stats' Gmakefile) I can't see any difference.
> > 
> > That wouldn't help; it's the awk script which is doing the relevant
> > calculations. Looking at the awk script, I don't think that compiling
> > awk with -ffloat-store would help.
> > 
> > The problem is basically that none of the axioms of real arithmetic
> > quite hold true for floating-point arithmetic.
> > 
> > In this particular case, the variance is:
> > 
> > 	(SUM[x(i)^2] - SUM[x(i)]^2/N)/N
> > 
> > If all of the x(i)'s are identical, then the two sides of the
> > subtraction should be equal. However, given the error introduced by
> > each floating point operation, they will differ slightly; if the RHS
> > is slightly larger, then result will be negative, causing the sqrt()
> > to fail.
> > 
> > Whilst there probably is an algorithmic solution (a decent textbook on
> > numerical methods should list many approaches for working around the
> > problems inherent in floating-point subtraction), I'm not sure whether
> > that would be justified, or if it would be better to just include an
> > explicit check for a negative variance (indicating a variance so low
> > that it's been swamped by rounding error).
> 
> Perhaps I am wrong, but doesn't the problems in r.univar result from
> the input data?

No. Or at least this problem doesn't relate to wrong input data:

> Variance:  -2.49502e-12
> awk: cmd. line:17: (FILENAME=- FNR=18894) warning: sqrt called with negative argument -2.49502e-12
> Standarddeviation: nan

If all of the input values are identical (even if they're wrong), the
variance and SD should be zero. Also, whatever the inputs, the
variance should never be negative.

> If r.mapcalc produces wrong numbers, the r.univar shouldn't
> be better. Here, the r.univar replicates the values generated from
> r.mapcalc. So the problem may be in r.mapcalc.

The variance/SD problem is in one of the following, depending upon
your point of view:

a) the awk script "r.univar",
b) awk generally, or
c) floating-point arithmetic generally.

> The range results of r.mapcalc are quite strange with FP numbers.

Not that strange; with the "-1" flag, r.stats uses "%.10f" to print
the values; without it, it uses "%10f" then trims the result (this is
equivalent to "%.6f", at least for GNU libc).

> Even if I compile r.mapcalc with -ffloat-store the result is odd:

"-ffloat-store" was basically a hunch regarding the variance/SD issue,
which turned out to be incorrect. It will fix some FP problems
(specifically, when subtracting identical values doesn't give zero);
it won't fix all of them (in this case, subtracting values which
should be identical but aren't quite).

> r.mapcalc test=1.1
> CREATING SUPPORT FILES FOR test
> range: 1.1000000238 1.1000000238

This is to be expected; r.mapcalc also uses "%.10f", which is more
precision than a float has (float will be promoted to double when
passed as an unprototyped argument, as is the case with printf()).

-- 
Glynn Clements <glynn.clements at virgin.net>

---------------------------------------- 
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo at geog.uni-hannover.de with
subject 'unsubscribe grass5'