Poisson Gap Sampling

Poisson Gap Sampling Scheduling

more ~/poisson_1d_c.txt

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

        //////////////////////////////////////////////////
        // algorithm poisson random number (Knuth):     //
        // init:                                        //
        //    Let L ← exp(-lambda), k ← 0 and p ← 1.    //
        // do:                                          //
        //    k ← k + 1.                                //
        //    Generate uniform random number u in [0,1] //
        //    and let p ← p × u.                        //
        // while p ≥ L.                                 //
        // return k − 1.                                //
        //////////////////////////////////////////////////

int     poisson ( double lmbd )
{
        double  L = exp( -lmbd );
        int     k = 0;
        double  p = 1;

        do {
                double  u = drand48();
                p *= u;
                k += 1;
        } while ( p >= L );

        return( k-1 );
}

int     main ( int argc, char** argv )
{

        int     i;
        float   s = atof( argv[1] );  //  initial seed     //
        int     p = atoi( argv[2] );  //  sampling points  //
        int     z = atoi( argv[3] );  //  total size       //
        int     *v;                   //  vector of sampling points //
        float   ld = (float) z / (float) p;
        float   w = 2.0;              //  inital weight    //
        int     k = 0;                //  currently found sampling points //

        v = ( int* ) malloc( z*sizeof( z ) );

        srand48( s );                 //  initialize seed //

        do {
                i = 0; k = 0;         //  N.B. first point always acqured //
                                      //  we use the nomenclature of first point == 0 //

                while ( i < z ) {
                        v[k] = i;     //  store point to be acquired //
                        i += 1;       //  put pointer at next point  //
                        k += 1;       //  next index of acqured point //

                        //  Now make a gap : //

                        i += poisson( (ld-1.0)*w*sin((float)i/(float)(z+1)*M_PI) );
                }

//  if more points T.B. acquired found than than wanted: try again with weight 2% larger //
// if fewer points T.B. acquired found than than wanted: try again with weight 2% smaller //

                if ( k > p ) w *= 1.02;
                if ( k < p ) w /= 1.02;

        } while ( k != p ); // try until correct number points to be acquired found //

//  print the data on standard output //

        for ( k = 0 ; k < p ; k++ ) printf( "%d\n", v[k] );

}