/*   Metropolis_Gaussian_Random.c                                   */
/*   Demonstration of random numbers                                */
/*   Metropolis Algorithm                                           */
/*   Gaussian Random Numbers distribution                           */ 
/*                                                                  */

/*   Syntax:                                                        */
/*   Metropolis_Gaussian_Demo 10000 > data_file                     */

                                                        /* Michel Vallieres */

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


/* probability profile of the random numbers */
double weight( double x, double y )
{
  return exp( - ( x*x + y*y ) );
}


/* Metropolis algorithm to generate Random       */
/* Numbers in the plane with probability profile */
/* weight(x,y)                                   */
void metropolis ( double *x, double *y, double (*weight)(double,double), double delta )
{
  double xt, yt, xx, yy, eta, ratio;
  
  xx = (double)rand()/(double)RAND_MAX;          /* new random numbers 0 to 1 */
  yy = (double)rand()/(double)RAND_MAX;
  xt = *x + delta * ( 2*xx - 1.0 );              /* in range -1 to 1 */
  yt = *y + delta * ( 2*yy - 1.0 );
  ratio = weight( xt, yt ) / weight( *x, *y );   /* ration of probability profiles */
  eta = (double)rand()/(double)RAND_MAX;         /* random # in 0-1  */
  if ( ratio > eta )                             /* Metropolis selection */
    {
      *x = xt;                                   /* selected numbers */
      *y = yt;
    }
}


/* demo program to illustrate Gaussina random numbers in the plane */
int main( int argc, char *argv[] )
{
  int    i, count, N;
  double x, y, x_old, y_old, delta;
                                                 /* argument list -> N */
  if ( argc != 2 )
    {
      printf( "\n" );
      printf( "Syntax: Metropolis_Gaussian_Random  N \n" );
      exit(1);
    }
  N = atoi( argv[1] );
                                                /*  Metropolis delta - size of jump */
  delta = 0.2;
  x = 0.5;                                      /* initila position in plane */
  y = 0.6;
  count = 0;                                    /* counter for selected numbers */
  printf( " %f %f\n", x, y );


  for ( i=0 ; i<N ; i++ )
    {
      x_old = x;
      y_old = y;
      metropolis ( &x, &y, weight, delta );     /* Metropolis Algorithm */
      if ( x != x_old && y != y_old ) 
          count++;
      printf( " %f %f\n", x, y );
    }
                                                /* Acceptance ratio - should be > ~0.5 */
  fprintf( stderr, "\n Acceptance ratio: %f\n\n", (double)count/(double)N );

}