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

[Markus Neteler]

On Wed, May 16, 2001 at 04:36:16AM +0100, Glynn Clements wrote:
> 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.

Glynn, I agree. Was a little tired yesterday.

> > 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",
BTW: A C implementation would be better (maybe using s.univar). Would
be a nice exercise, but I don't have the time.

> 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()).

O.k. I think I understand the problem now. But how to deal with it?
r.mapcalc seems to generate FCELL only. Shouldn't the range printed with
float precision then? Sorry for insisting here, but above range will
be expected to be wrong by the average user not familiar with precision
issues.

And r.stats: maybe the printed precision should be dependent from the
input map type, means more decimals for DCELL maps that FCELL maps.

I am no expert here, this is just my impression,

regards

 Markus

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