/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/


/*---------------------------------

randomWalkDistanceRMSD.c
Computes mean and variance between corresponding
points between 2 random walks. Used for algorithms
problem.
Also computes RMSD and max to see if max is
Gumbel distributed

Additionally to randomWalkDistance.c generates optionally data on root
mean square deviation between 2 random walks on 2 dimensional lattice.

--------------------------------*/

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

// if RMSD defined compute RMSD, otherwise MSD.
#define RMSD 1

double rreal();
int randint(int);
 
main(int argc, char *argv[]) {
  int i,j;
  int M = 50000; // number of runs
  int n; // n points in walk
  int x,y,u,v;
  double distance;
  double maxDist; // max distance between success pts between both walks
        // want to see if this term dominates, so get Gumbel distribution
  double distSum; // sum of distances between successive points on walk
  double sum=0,sumSq=0; 
    // used to compute mean and variance of distances between walks
  double sumMax=0,sumSqMax=0; 
    // used to compute mean and variance of max of distances between walks

  /*--------error handling for input */
  if (argc != 2) {
     fprintf(stderr,"%s int\n",argv[0]);
     exit(1);
     }
  
  n = atoi(argv[1]);

#if RMSD
  printf("MaxDist\t%RMSD\n");
#else
  printf("MaxDist\t%AveDist\n");
#endif


  for (i=0;i<M;i++){
    
    /* generate pair of walks */
    x=y=u=v=0; //start from origin
    distSum=0;
    maxDist=0;
    for (j=0;j<n;j++){
       switch( randint(3) ) {
         case 0: x += 1; break;
         case 1: x -= 1; break;
         default: break;  // x = x
       } 
       switch( randint(3) ) {
         case 0: y += 1; break;
         case 1: y -= 1; break;
         default: break;  // y = y
       } 
       switch( randint(3) ) {
         case 0: u += 1; break;
         case 1: u -= 1; break;
         default: break;  // u = u
       } 
       switch( randint(3) ) {
         case 0: v += 1; break;
         case 1: v -= 1; break;
         default: break;  // v = v
       } 
       distance = sqrt((x-u)*(x-u) + (y-v)*(y-v));
       if (distance > maxDist) maxDist = distance;

#if RMSD
       distSum += (x-u)*(x-u) + (y-v)*(y-v);
#else
       distSum += sqrt((x-u)*(x-u) + (y-v)*(y-v));
#endif
    } /* end for loop to generate pair of walks*/

    // update auxiliary variables to compute mean,var of maxDist
    sumMax += maxDist;
    sumSqMax += maxDist*maxDist;

#if RMSD
    printf("%lf\t%lf\n",maxDist,sqrt(distSum/n));
#else
    printf("%lf\t%lf\n",maxDist,distSum/n);
#endif

#if RMSD
    /* RMSD of current pair of walks is the
       root of the average square distance */
    sum += sqrt(distSum/n);
    sumSq += distSum/n; 
#else
    sum += distSum/n;
    sumSq += (distSum/n * distSum/n); 
#endif

  } /* end of for loop to generate M pairs of walks */

  /*
    output results 
  */

#if RMSD
  printf("\n\nExpected RMSD between 2 persons,\n");
#else
  printf("\n\nExpected distance between 2 persons,\n");
#endif

  printf("both of whom perform a random walk\n");
  printf("of %d steps is:",n);
  printf(" %lf\n\n", sum/M);
  printf("\n\nVariance for this random variable is: %lf\n",
           sumSq/M - sum/M*sum/M);
  printf("Expected maxDist:%lf\t with var:%lf\n",sumMax/M,
           sumSqMax/M - sumMax/M*sumMax/M);
}


double rreal() {
  return ((double)rand())  / (RAND_MAX+1.0);
}

int randint(int m) {
  return (int) (rreal()*m);
}