/*   Metropolis_Gaussian_Parallel.c                                 */
/*   Demonstration of integration by random numbers                 */
/*   Metropolis Algorithm                                           */
/*   Gaussian Random Numbers distribution                           */
/*   Parallel version                                               */ 
/*                                                                  */

/*   Syntax:                                                        */
/*   mpiexec -np N Metropolis_Gaussian_Parallel  -p  >  data_file   */

                                          /* Michel Vallieres */

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

#define NMAX 1000   /* initial array dimensions    */
                    /* will be modified if need be */

double x00, y00;    /* initial position of random walk */


/* function to integrate */
double func( double x, double y )
{
  return M_PI * ( x*x + y*y );
}



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



/* Metropolis algorithm to generate Random       */
/* Points in the plane with probability profile  */
/* weight(x,y)                                   */
void Metropolis ( double x[], double y[], double (*weight)(double,double),  
                  double delta, int N, int *accepted_count )
{
  double xt, yt, xx, yy, eta, ratio;
  int local_count, count;
 
  count = 0;                                           /* accepted count */

  for( local_count=0 ; local_count < N-1 ; local_count++ )
    { 
       xx = (double)rand()/(double)RAND_MAX;           /* random # in 0,1 */
       yy = (double)rand()/(double)RAND_MAX;
       xt = x[local_count] + delta * ( 2*xx - 1.0 );   /* in range -1 to 1 */
       yt = y[local_count] + delta * ( 2*yy - 1.0 );   /* trial step */
                                                       /* probability ratio */
       x[local_count+1] = x[local_count];              /* new step = old step */
       y[local_count+1] = y[local_count];
       ratio = weight( xt, yt ) / weight( x[local_count] , y[local_count] );   
       eta = (double)rand()/(double)RAND_MAX;          /* random # in 0-1  */
       if ( ratio > eta )                              /* Metropolis selection */
         {
            count++;
            x[local_count+1]  = xt;                    /* selected step */
            y[local_count+1]  = yt;
         }

       // fprintf(stderr,"%d %f %f %f %f %d %d\n", 
       //      local_count,xt,yt,ratio,eta,N,count);

    }
  *accepted_count = *accepted_count + count;
}



/* demo program to illustrate Monte-Carlo integration using */
/* Gaussian random numbers in the plane                     */
/* Parallel version                                         */
int main( int argc, char *argv[] )
{
  int    round, ir, ir_min, accepted_count, N, N_total, N_input;
  int    plot_control, first, ir_first;
  double *x, *y, *x_local, *y_local;
  double delta;
  double f, sumf, sumf2, local_sumf, local_sumf2,
         total_sumf, total_sumf2, integral, sigma;
  int    node, myid, numprocs;
  int    current_x_y_size, N_average, leftover, number, nmax;
  MPI_Status recv_status;

                                              /* join the MPI virtual machine */
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

                                              /* adjust NMAX for # nodes */
  nmax = NMAX - ( NMAX % ( numprocs-1) ); 
                                              /* current size of x & y */
                                              /* arrays - reserve memory */
  current_x_y_size = nmax;
  x = (double *)malloc( current_x_y_size * sizeof(double) );
  y = (double *)malloc( current_x_y_size * sizeof(double) );
  x_local = (double *)malloc( current_x_y_size * sizeof(double) );
  y_local = (double *)malloc( current_x_y_size * sizeof(double) );

  x00 = 0.5;                                  /* initial position in plane */
  y00 = 0.6;

                                              /* MASTER NODE */
  if ( myid == 0 )
    {
                                              /* plot control */
      plot_control = 0;
      if ( argc == 2 )
	if ( strcmp( argv[1],  "-p") == 0  )
	  plot_control = 1;

      delta = 0.2;                            /* Metropolis delta - size of jump */

      N_total = 1;                            /* total # of random numbers */
      accepted_count = 1;                     /* counter for selected numbers */  
      first = 1;

      f = func(x00,y00);                      /* Monte-Carlo sums */
      sumf = f;
      sumf2 = f*f ;

                                              /* accuracy improvement */
      fprintf( stderr, "\n Enter # of random number to use ( 0 to quit ): " );
      while ( scanf( "%d", &N_input ) != EOF && N_input > 0 )
	{

                                              /* find # per processor */  
	  N_average = N_input / ( numprocs -1 );
	  leftover =  N_input % ( numprocs -1 );
                                              /* add to total */
	  N_total = N_total + N_input - first;
                                              /* generate & send random */
                                              /* numbers to each node */
	  for ( node=1 ; node<numprocs ; node++ )
	    {
                                             /* number for each slave */
	      number = N_average;
	      if ( node <= leftover )
		number++;

	      if ( first == 0 ) number++;     /* adjust count on MC steps */
	      first=0;                        /* on subsequent calls      */

	      MPI_Ssend(&number, 1, MPI_INT, node, 121, MPI_COMM_WORLD);

                                              /* more space if needed */
	      if ( number > current_x_y_size )
		{
		  free( x );
		  free( y );
		  current_x_y_size = number;
		  x = (double *)malloc( (current_x_y_size+100) * sizeof(double) );
		  y = (double *)malloc( (current_x_y_size+100) * sizeof(double) );
		}
                                              /* Metropolis Algorithm */
	      x[0] = x00;
	      y[0] = y00;
	      Metropolis ( x, y, weight, delta, number , &accepted_count );

                                              /* send x & y to nodes */
                                              /* send x & y to nodes */
	      MPI_Ssend( x, number, MPI_DOUBLE, node, 123, MPI_COMM_WORLD );
	      MPI_Ssend( y, number, MPI_DOUBLE, node, 124, MPI_COMM_WORLD );

                                                /* plot info */
	      if ( plot_control == 1 )
		{
		  ir_first =  0;
		  if ( first == 0 ) ir_first = 1;
		  for ( ir=ir_first ; ir<number ; ir++ )
		    printf( " %f %f %f\n", x[ir], y[ir], func(x[ir], y[ir]) );
		}
                                              /* ready for next bunch of */
                                              /* random numbers */
	      x00 = x[number-1];
	      y00 = y[number-1];

	    }        /* for */
                                              /* get sums from nodes */
      local_sumf = 0.0;
      MPI_Reduce ( &local_sumf, &total_sumf, 1, MPI_DOUBLE,
                                MPI_SUM, 0, MPI_COMM_WORLD );
      local_sumf2 = 0.0;
      MPI_Reduce ( &local_sumf2, &total_sumf2, 1, MPI_DOUBLE,
                                MPI_SUM, 0, MPI_COMM_WORLD );

                                              /* add to global totals */
                                              /* add to global totals */
      sumf = sumf + total_sumf;
      sumf2 = sumf2 + total_sumf2;
                                              /* finalize Monte-Carlo calc */
      fprintf( stderr, " Sumf - Sumf2 = %f %f\n", sumf, sumf2 );
      integral = sumf/(double)N_total;
      sigma = sqrt(( sumf2/(double)N_total - integral*integral )/(double)N_total );
      fprintf( stderr, " Integral = %f  +/- sigma = %f \n", integral, sigma );

      fprintf( stderr, "\n Enter # of random number to use ( 0 to quit ): " );

    }  /* while scanf */
                                              /* siganl -> nodes to suicide */
    number = -1;
    for ( node=1 ; node<numprocs ; node++ )
      MPI_Ssend(&number, 1, MPI_INT, node, 121, MPI_COMM_WORLD);

                                             /* Acceptance ratio - should be > ~0.5 */
    fprintf( stderr, "\n Acceptance ratio: %d %d %f\n\n",
                       N_total, accepted_count, (double)accepted_count/(double)N_total );

                                             /* end main code */
    free( x );
    free( y );
    MPI_Finalize();
    exit(0);

    }    /* if my_node */
                                             /* WORKER NODES  */
  else
    {
                                             /* loop over all bunches of */
                                             /* random numbers */
      for ( ; ; )
        {
                                             /* receive signal or # of numbers */
	  MPI_Recv( &number, 1, MPI_INT, 0, 121,
                               MPI_COMM_WORLD, &recv_status);

                                             /* end code - suicide */
	  if ( number <0 )
	    {
	      free( x_local );
	      free( y_local );
	      MPI_Finalize();
	      exit(0);
	    }
	  else if ( number > 0 )            /* perform MC integral */
	    {
                                            /* more space if needed */
	      if ( number > current_x_y_size )
		{
		  free( x_local );
		  free( y_local );
		  current_x_y_size = number;
		  x_local = (double *)malloc( current_x_y_size * sizeof(double) );
		  y_local = (double *)malloc( current_x_y_size * sizeof(double) );
		}

	      MPI_Recv( x_local, number, MPI_DOUBLE, 0, 123,
                               MPI_COMM_WORLD, &recv_status);
	      MPI_Recv( y_local, number, MPI_DOUBLE, 0, 124,
                               MPI_COMM_WORLD, &recv_status);
                                             /* Monte-Carlo integral */
	      local_sumf = 0.0;
	      local_sumf2 = 0.0;
	      for ( ir=1 ; ir<number ; ir++ )
		{
		  f =  func( x_local[ir], y_local[ir] );
		  local_sumf = local_sumf + f;
		  local_sumf2 = local_sumf2 + f*f;
		}
                                             /* send partial sums to node 0 */
	      MPI_Reduce ( &local_sumf, &total_sumf, 1, MPI_DOUBLE,
                                MPI_SUM, 0, MPI_COMM_WORLD );
	      MPI_Reduce ( &local_sumf2, &total_sumf2, 1, MPI_DOUBLE,
                                MPI_SUM, 0, MPI_COMM_WORLD );
	    }
        }

    }

}