/* This C program uses the theory from Rosenthal (1995b) applied to the
   frog example to find a bound on the total variation distance after 200
   iterations */

#define NODES 100 /* the size of the state space */
#define M 4       /* P has edges between nodes within M of each other */
#define ITERS 200 /* the number of jumps before we take an approximate
                     sample from the stationary distribution */
#define BETA 2    /* the target distribution is proportional to
                     e^(-((distance to node zero)^BETA)) */
#define HEXP 2    /* the exponent in the calculation of h below */

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

int distance(int);

int main()
{
  int min, i, j, x, y, disti, distj, distx, disty;
  double P[NODES][NODES], pi[NODES];
  double bestbound=1, Eh0=0, Eh1=0, alph=100, epsilon, A=0, denom=0,
     acceptalpha, bound, temp;

  /* the following two loops calculate the stationary
     distribution of the chain */
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    denom = denom + exp(-pow(min, BETA));
  }
  for (j=0; j<NODES; j++)
  {
    min = distance(j);
    pi[j] = exp(-pow(min, BETA))/denom;
  }

  /* calculates a transition matrix with the desired stationary
     distribution */
  for (i=0; i<NODES; i++)
  {
    P[i][i] = 1.0/(2*M+1);
    for (j=0; j<NODES; j++)
    {
       if (pi[j] > pi[i])
          acceptalpha = 1;
       else
       {
          disti = distance(i);
          distj = distance(j);
          /* to avoid dividing zero by zero; this makes P slightly
             different than actually specified, but the difference is
             minor, and the same thing had to be done in the eigenvalue
             calculation */
          if ((pi[i]<pow(10, -317)) && (distj>disti))
            acceptalpha = 0;
          else if ((pi[i]<pow(10, -317)) && (distj<=disti))
            acceptalpha = 1;
          else
            acceptalpha = pi[j]/pi[i];
       }

       if ((0 < abs(i-j)) && (abs(i-j) < M+1))
       {
          P[i][j] = acceptalpha/(2*M+1);
          P[i][i] = P[i][i] + (1-acceptalpha)/(2*M+1);
       }
       else if (abs(i-j) >= NODES-M)
       {
          P[i][j] = acceptalpha/(2*M+1);
          P[i][i] = P[i][i] + (1-acceptalpha)/(2*M+1);
       }
       else if ((M+1 <= abs(i-j)) && (abs(i-j) < NODES-M))
          P[i][j] = 0;
    }
  }

  /* finds expected value of h after 1 iteration and alph using
     h = 1 + dist(x,0)^HEXP + dist(y,0)^HEXP and taking the low value of
     temp since inequality holds for all x, y
     covering all pairs of x and y isn't necessary with this particular h
     since x will equal y, but it's nice to have the code in case we want
     to change h */
  for (x=1; x<NODES; x++)
  {
    for (y=x; y<NODES; y++)
    {
      Eh1 = 0;
      for (i=0; i<NODES; i++)
      {
        min = distance(i);
        /* adding the two sums at the same time */
        Eh1 += pow(1.0*min, HEXP)*(P[x][i]+P[y][i]);
      }
      Eh1++;
      distx = distance(x);
      disty = distance(y);
      temp = (1+pow(distx, HEXP)+pow(disty, HEXP))/Eh1;
      if (temp < alph)
        alph = temp;
    }
  }

  epsilon = 1; /* for R = only one state */

  /* finds A; we're using k0 = 1 */
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    A += pow(min, HEXP)*P[0][i];
  }
  A = A*2 + 1;

  /* finds expected value of h for the initial distribution */
  for (i=0; i<NODES; i++)
  {
    min = distance(i);
    Eh0 += pow(min, HEXP)*pi[i];
  }
  Eh0 += 2; /* since the initial state is always zero here */

  /* tries different j's to find the best bound; in our case, the best j
     will always be 1 since epsilon is 1, but it's nice to have the code
     here if other cases are to be explored */
  for (j=1; j<=ITERS; j++)
  {
    bound = pow(1-epsilon, j) +
             Eh0*pow(A, j-1)*pow(alph, -ITERS+(j-1));
    if (bound < bestbound)
      bestbound = bound;
  }

  printf ("Bound is %g\n", bestbound);
  exit(0);
}

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