\ Forth tools for random lookup tables
\
\ ---------------------------------------------------
\     (c) Copyright 1998  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------
\
\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependencies:
\       Assumes independent floating point stack

FALSE [IF]

 Algorithm:

    If we need to sample a given random distribution, but have
    only a prng producing numbers y[k] evenly distributed on [0,1],
    we must solve the (generally transcendental) equation

        P(x<X[k]) = y[k]

    where P(x<X) is the cumulative distribution. In some cases
    P(x<X) is not known in closed form and must be obtained by
    numerical quadrature. Thus the inversion process can become
    a serious bottleneck in a large simulation, where many samples
    are required.

    The defining word DIST: creates a table with N entries

        x[k] = P^{-1} ( k/N ), k=0, 1, ... , N-1

    The lookup code defined by DOES> employs a prng with a long
    cycle to traverse the table in random order. This is equivalent
    to approximating the distribution function p(x) by a set of step
    functions (i.e. a histogram) of uneven width.


[THEN]


MARKER -dist

: undefined    BL WORD  FIND  NIP  0=  ;

undefined s>f   [IF] : s>f   S>D  D>F  ;    [THEN]
undefined f>s   [IF] : f>s   F>D  D>S  ;    [THEN]

undefined dxor  [IF]
    : dxor  ROT XOR >R XOR R>  ;    \ double xor
[THEN]

undefined  use( [IF]
\ Vectoring: for using function names as arguments
: use(      '       \ state-smart ' for syntactic sugar
    STATE @  IF  POSTPONE LITERAL  THEN  ;  IMMEDIATE

' NOOP  CONSTANT  'noop
: v:   CREATE  'noop  ,  DOES> PERFORM  ;   \ create dummy def'n
: 'dfa   ' >BODY  ;                         ( -- data field address)
: defines    'dfa   STATE @
             IF   POSTPONE  LITERAL    POSTPONE  !
             ELSE   !   THEN  ;  IMMEDIATE
\ end vectoring
[THEN]


include prng.f

include ansfalsi.f

\ program begins here

: sf,   \ "comma in" an IEEE 32-bit number
    HERE  SF!  1 SFLOATS  ALLOT  ;

: random_seed   ( -- d)
    TIME&DATE                   ( sec min hour day month year)
    2DROP 2DROP  60 *  +  S>D  ;

: random_ndx    ( adr -- index)
    LOCALS| adr |
    adr  CELL+                      ( -- adr+cell)
    (rand)                          ( f: -- seed')
    FDUP  F0<                       ( f: -- seed')  ( flag)
    IF    bigdiv  F+   THEN         \ normalize
    bigdiv  F/                      ( f: prng)
    adr @   s>f                     ( f: prng n)
    F*  f>s  ;                      ( index)

: )dist:    ( xt n --)              \ dist_name converts xi to X(xi)
    LOCALS| n xt |
    CREATE                          \ make entry under dist_name
        n  ,                        \ save table size
        random_seed  , ,            \ initialize and save seed
        n  s>f                      ( f: -- n)
        0e  sf,                     \ store  0.0e as first entry
        n 1 DO                      \ entries from I=1 to n-1
            I s>f   FOVER   F/      ( f: -- n I/n)
            xt EXECUTE              ( f: -- n X[I/n])
            sf,                     \ store in table
        LOOP
        FDROP                       ( f: --)

    DOES>   DUP                     ( -- base_adr base_adr)
            random_ndx              ( -- base_adr index)
            SFLOATS +   3 CELLS +   \ denormalize
            SF@   ;                 ( f: -- variate)


\ Usage:   use( inverse_dist_name  Table_size )dist:  dist_name

\ Example: inverse Poisson distribution

\ FVARIABLE xi

\ : f1   ( f: x -- 1-[1+x]/e^x-xi)
\       1e  FOVER  FLN  F-  F*  FNEGATE  1e  F+   xi F@  F-   ;

\ : invP      ( f: xi -- X[xi])
\       xi F!
\       use( f1  1e-10  1e  1e-6 )falsi    ( f: -- e^{-X})
\       FLN  FNEGATE  ;

\ Now create a table:
\   use( invP 64 )dist:   poisson
\   program to test with
\   : variates    ( N )   0 DO  poisson  F.  LOOP   ;