/* The following C program studies the speed of convergence of a matrix to
   its stationary distribution using the Metropolis-Hastings algorithm.
   The 'frog example' is used again */

#define NODES 20      /* size of the state space */
#define ITERATIONS 40 /* number of jumps before we take an approximate
                         sample from the stationary distribution */
#define RUNS 500000   /* number of trials that are averaged together to
                         estimate the stationary distribution */
#define M 5           /* all nodes within M of each other have an edge
                         between them */
#define BETA 2        /* target distribution is proportional to (1 +
                         distance to node zero)^-BETA */

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

void simulate();
int distance (int);

int main()
{
  srand((unsigned) time(NULL));
  simulate();
  exit(0);
}

void simulate()
{
  int edges[NODES][NODES], total[NODES], sum[NODES];
  int current, temp, partsum, proposal, c, i, j, rejects, min;
  double sdist[NODES], target[NODES];
  double rejsum=0, vardist=0, l2distpi=0, l2distpiinv=0, denom=0;
  double randomnum, randomnum2, alpha, rejperc;
  FILE *fp;
  fp = fopen("outputfile", "w");

  /* setting all values in arrays to zero so we can add them */
  for (i=0; i<NODES; i++)
  {
    for (j=0; j<NODES; j++)
      edges[i][j] = 0;
    total[i] = 0;
    sum[i] = 0;
  }

  /* setting all the edges with equal weight */
  for (i=0; i<NODES; i++)
    for (j=0; j<=M; j++)
    {
      temp = (i+j) % NODES;
      edges[i][temp] = 1;
      edges[temp][i] = 1;
    }

  /* finding the number of edges from each node (actually just 2M+1 in
     this case, but written more generally to allow changes to code) */
  for (i=0; i<NODES; i++)
    for(j=0; j<NODES; j++)
      sum[i] += edges[i][j];

  /* calculating the given target distribution for NODES and BETA
     we don't actually need pi - only the ratios, but it's good to have it */
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    denom += pow(min+1, -BETA);
  }
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    target[i] = pow(min+1, -BETA)/denom;
  }

  /* the actual simulation */
  for (i=0; i<RUNS; i++)
  {
    rejects = 0; /* keeps track of rejections to find what percentage of
                    proposals are rejected (this hopefully leads to a
                    connection to 'optimal scaling'), resets often to
                    avoid overflow */
    current = 0; /* starting the simulation at node 0 here */
    for (j=0; j<ITERATIONS; j++)
    {
      partsum = 0;
      randomnum = rand()/(RAND_MAX+1.0)*sum[current];
      /* starting the loop M spots away from the current node maximizes
         efficiency since this is where all the edges are */
      c = (current + NODES - M) % NODES;
      while (partsum < randomnum)
      {
        partsum += edges[current][c];
        proposal = c;
        c = (c+1) % NODES;
      }

      randomnum2 = rand()/(RAND_MAX+1.0);
      /* sum is q's reciprocal in this case */
      alpha = target[proposal]*sum[current]/(target[current]*sum[proposal]);
      if (randomnum2 < alpha) /* accepting the proposed jump */
        current = proposal;
      else
        rejects++;
    }
    rejsum += ((double)rejects)/ITERATIONS;
    total[current]++;
  }

  /* calculating the different types of norms to quantitatively study how
     well we approximated the target distribution */
  for (i=0; i<NODES; i++)
  {
    sdist[i] = ((double)total[i])/RUNS;
    vardist += fabs(sdist[i]-target[i]);
    l2distpi += pow(sdist[i]-target[i], 2.0)*target[i];
    l2distpiinv += pow(sdist[i]-target[i], 2.0)/target[i];
  }
  vardist = 0.5*vardist;
  l2distpi = sqrt(l2distpi);
  l2distpiinv = sqrt(l2distpiinv);
  rejperc = 100*rejsum/RUNS;
  fprintf (fp, "Variation distance is %g\n", vardist);
  fprintf (fp, "L2 dist(pi) is %g\n", l2distpi);
  fprintf (fp, "L2 dist(1/pi) is %g\n", l2distpiinv);
  fprintf (fp, "Rejection percentage is %g\n", rejperc);
  fclose(fp);
}

int distance(int i)
{
  if (NODES-i < i)
    return NODES-i;
  return i;
}