[GRASS-dev] random issues and bug
glynn at gclements.plus.com
Fri Aug 24 06:49:22 CEST 2007
> sorry this is a bit disorganized, but there are several issues.
> But there is a common theme. :)
> r.malcalc rand(a,b) for integer values of a,b
> For a cellular automata seed I was trying to make a random 0/1 base map.
> G63> r.mapcalc "rand.boolean = rand(0,1)"
> G63> r.info -r rand.boolean
> G63> r.mapcalc "rand.boolean = rand(0,2)"
> G63> r.info -r rand.boolean
> I assume what is happening here is that the range for rand(0,2) is 0.0 to
> 1.999999 and the conversion to integer (CELL map) is done by a cast which
> throws away the least significant bits. But to the casual user this seems
> very counter intuitive.
There is no cast; the rand() function has separate CELL/FCELL/DCELL
> r.mapcalc manual says:
> rand(a,b) random value x : a < x < b
> ... but that's only true if a,b are floats. For integer values it is more
> like: a <= x <= b-1
The manual page is wrong. The code (intentionally) calculates values x
such that "a <= x < b".
I've fixed this in CVS.
> My first attempt was:
> G63> r.mapcalc "rand.boolean = int( rand(0,1) + 0.5 )"
> but that made an all 0s map as rand(n,n+1) == n as mentioned above.
> (I didn't expect the cast to int)
"0" and "1" are integer literals, so the above will use the integer
version of rand().
> In hindsight I can get that method to work by passing rand(,) two floats:
> G63> r.mapcalc "rand.boolean = int( rand(0.0,1.0) + 0.5 )"
> G63> r.info -r rand.boolean
That's inefficient, but it will work.
Ultimately, if you want 0 or 1, use rand(0,2).
> Bug: r.mapcalc rand() not acting randomly. (in at least 5.4 - 6.3cvs)
> #simple_xy location
> export GRASS_WIDTH=$BOXSIZE
> export GRASS_HEIGHT=$BOXSIZE
> d.mon x0
> g.region n=$BOXSIZE s=0 w=0 e=$BOXSIZE res=1
> r.mapcalc "rand.boolean = rand(0,2)" # make a map of 0s and 1s
> d.rast rand.boolean
> screenshot: (1:1 region rows/cols to image rows/cols)
> vertical bands, the rows have a repeating pattern to them.
> weird-- in GRASS 6.0.2 each row is one cell shorter and the pattern
> drifts to the left (box is 1:1 so ends up in opp corner).
> In 5.4.0 it's the same as both 6.3cvs and 6.2.2, ie vertical bands.
> My first guess was that this arrived with the addition of $GRASS_RND_SEED
> and the seed was reset to the same value at the start of each row when it
> should really only be set once at the start of the module.
It is only set once at the start of the module.
> After short reflection that's not the case (it's not a smear) but there
> certainly is a repeating pattern. Aliasing effect? Resizing the xmon
> doesn't get rid of the pattern, it's in the cells. By setting the region
> to be one bigger than 512x512 the effect goes away:
> g.region n=$((1+ $BOXSIZE)) s=0 w=0 e=$((1 + $BOXSIZE)) res=1
> r.mapcalc "rand.boolean = rand(0,2)"
> eh? power of 2 issue?
Yep; it's an artifact of the PRNG.
If you want "very" random numbers, use a Mersenne Twister, but it's
quite a bit slower than an LCG.
> to back up a little, my first attempt to make the 0/1 map was with
> r.random, but that requires a base map and is for placing points.
> Wrong module.
> Next I tried r.mapcalc and ran into above issues.
> Next I tried r.surf.random. It creates a DCELL map with data in the range
> of two integers.
> One thing which struck me as odd is that r.surf.random creates a DCELL map
> from two integer command line options. If you create a map with floating
> points numbers, why not allow min and max to be any real number??
Probably overlooked when r.surf.random was extended from CELL to
> I've just now added an -i flag to r.surf.random to write out a CELL map.
> I made that be in the range of [min <= x <= max] for int, and for DCELL
> kept (min < x < max) [or whatever rand() gives]. As this doesn't match
> r.mapcalc behaviour (see above) I think we need to change either one or
> the other module for consistency's sake, then document the choice well.
> My first attempt to add the -i flag was like:
> *(row_out_C + col_count) = (CELL) floor(rand1(2742)*(max-min)+min +0.5);
> but that under-represents first and last bin.
> e.g. for min=0 max=10 the 0 and 10 bins only have half the cells of 1-9.
> G63> r.surf.random.MOD rsrand.i -i min=0 max=10
> G63> d.histogram rsrand.i
> G63> r.stats -c rsrand.i
> 0 13032
> 1 26190
> 2 26142
> 3 26262
> 4 26336
> 5 26250
> 6 26171
> 7 26513
> 8 26262
> 9 25999
> 10 12987
> as the "1" bin draws from the range of 0.5-1.49999 and the 0 bin only
> has half that space (0.0-0.49999) to draw from.
> so I made the -i version as follows:
> *(row_out_C + col_count) = (CELL) (rand1(2742)*(max+1-min)+min);
> but that goes back to the question of should min=0, max=10 create 11 bins
> or 10?
10. 0 to 9, inclusive.
> For the integer case, I'd say include the end numbers; but that
> probably means changing r.mapcalc's rand(int,int) to be the same. ?????
r.mapcalc's behaviour is correct.
PRNGs invariably use start-closed/end-open intervals, i.e. a <= x < b.
More generally, if you partition the interval a-c into a-b and b-c,
you normally want b to be in one or the other but not both.
> As for Numerical Recipes code in the module (or not), it uses rand1()
> from the gmath library. lib/gmath/rand1.c
rand1() uses either rand() or drand48() from the standard C library.
There is no Numerical Recipes code involved here.
> How to cite that in the r.surf.random.html help page? I put:
> "It uses the random number generator drand48() or rand(), depending on
> the user's platform."
> but I'd much rather be specific and say those come from the GNU Project
> or as defined in POSIX 1.1 or C99 std, or libc6 or whatnot, but don't
> know the correct answer.
rand() is C89, drand48() is POSIX.
> as for drand48, the man page says:
> These functions are declared obsolete by SVID 3, which states
> that rand(3) should be used instead.
> So should we remove the "#if defined(HAVE_DRAND48)" parts from
ANSI C doesn't dictate the specific PRNG algorithm or the number of
bits, only that RAND_MAX is at least 32767. The end result is that
some rand() implementations are quite poor (I've even seen
implementations which take the *least* significant bits of an LCG,
which produces alternating even/odd values).
If lrand48() et al exist, they will always use a 48-bit state, which
puts a lower bound on the quality of the PRNG.
OTOH, if a system has lrand48(), it's likely that rand() will either
be an alias for lrand48() or something at least as good.
Glynn Clements <glynn at gclements.plus.com>
More information about the grass-dev