/*

BACKMET.C:  A computer program for running the "backmet" web browsing
algorithm developed in the paper "Achieving Limiting Distributions for
Markov Chains Using Back Buttons", by A. Feuerverger and J.S. Rosenthal,
available from

	http://probability.ca/jeff/research.html


Copyright (c) 2003 by Jeffrey S. Rosenthal.

Licensed for general copying, distribution and modification according to
the GNU General Public License (http://www.gnu.org/copyleft/gpl.html).

*/



#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <sys/types.h>

#define N 1000  /* Number of states. */
#define K 6  /* Number of extra links per state. */
#define C 10  /* To make target distribution pi[i] = i+C. */
#define TSERIESLENGTH 20

#define RUNLENGTH 40 /* number of k-blocks to run each time */
#define NUMREPS 100000000
#define PRINTLENGTH 10000

#define MAXOUTLINKS 30

int outlinks[N];
int linkplace[N][MAXOUTLINKS];
double qf[N][100];
double t[N];
int allpos = 0;
int i,j,h,r;

main()

{

    double pi[N];
    double sum;
    int k;
    int timespent[N][RUNLENGTH];
    int state, proposal;
    double tv;
    double absval();

    seedrnd();

    /* Initialise outlinks[i]. */
    for (i=0; i<N; i++) {
        outlinks[i] = 0;
    }

    /* Add base outlinks. */
    for (i=0; i<N-1; i++) {
        addlink(i,i+1);
    }
    addlink(N-1,0);

    /* Add random outlinks. */
    for (h=0; h<K*N; h++) {
        i = ifloor( N * drand48() );
        j = ifloor( N * drand48() );
        addlink(i,j);
    }

    /* Initialise pi[i]. */
    sum = 0.0;
    for (i=0; i<N; i++) {
        pi[i] = ((double)i+1+C);  /* Unnormalised target distribution. */
        sum = sum + pi[i];
    }

    /* Compute qf[i][0]. */
    for (i=0; i<N; i++) {
        qf[i][0] = pi[i]/sum;  /* Normalised target distribution. */
    }

    /* Find minimal value of k. */
    k=2;
    allpos = 0;
    while (allpos==0) {
	k++;
	computet(k);
    }

    printf("\nN = %d;   k = %d\n\n", N, k);

#ifdef VERBOSE
    for (i=0; i<N; i++)
      printf("t[%d] = %.14f\n", i, t[i]);
    for (i=0; i<N; i++) {
	printf("outlinks[%d] = %d\n", i, outlinks[i]);
    }
    for (j=0; j<outlinks[0]; j++) {
	printf("linkplace[0][%d] = %d\n", j, linkplace[0][j]);
    }
    for (i=0; i<N; i++) {
	printf("qf[%d][5] = %f;   t[%d] = %f\n", i, qf[i][5], i, t[i]);
    }
#endif

    /* Initialise timespent[i][j]. */
    for (i=0; i<N; i++) {
      for (j=0; j<RUNLENGTH; j++) {
        timespent[i][j] = 0;
      }
    }

    /* Do the runs. */
    for (h=0; h<NUMREPS; h++) {
      state = 0; /* Start in lowest-pi state. */
      for (j=0; j<RUNLENGTH; j++) {
	/* Select proposal uniformly from outlinks. */
	proposal = linkplace[state][ifloor(outlinks[state]*drand48())];
	/* Accept or reject it. */
	if ( drand48() < t[proposal] / t[state] 
		* outlinks[state] / outlinks[proposal] ) {
	    /* Accept the proposal. */
	    timespent[proposal][j] = timespent[proposal][j] + k;
/* printf("k=%d; timespent[%d][%d] = %d\n", k, proposal, j, timespent[proposal][j]); */
	    state = proposal;
	} else {
	    /* Reject the proposal. */
	    timespent[proposal][j] = timespent[proposal][j] + 1;
	    timespent[state][j] = timespent[state][j] + k-1;
	}
      }
      if ( h == PRINTLENGTH * (h/PRINTLENGTH) ) {
        printf("%d repetitions of %d blocks of %d time-steps completed.\n", 
						h, RUNLENGTH, k);
      }
    }

    /* Compute the results. */
    printf("\n");
    for (j=0; j<RUNLENGTH; j++) {
      /* Compute tv after j'th block of k time-steps. */
      tv = 0.0;
      for (i=0; i<N; i++) {
	tv = tv + absval( (((double)timespent[i][j])/NUMREPS/k) - qf[i][0] );
      }
      tv = tv/2;
      printf("On %d'th block of %d time-steps, T.V. = %f\n", j+1, k, tv);
    }

#ifdef VERBOSE
    printf("\n");
    for (i=0; i<N; i++) {
        printf("timespent[%d][%d] = %d; norm = %f; stat = %f\n", 
		i, 0, timespent[i][0],
		((double)timespent[i][0])/NUMREPS/k, qf[i][0]);
    }
    printf("\n");
    for (i=0; i<N; i++) {
        printf("timespent[%d][%d] = %d; norm = %f; stat = %f\n", 
		i, RUNLENGTH-2, timespent[i][RUNLENGTH-2],
		((double)timespent[i][RUNLENGTH-2])/NUMREPS/k, qf[i][0]);
    }
#endif

}


computet(int kk)
{
    /* Begin computation of t[i]. */
    for (i=0; i<N; i++) {
        t[i] = (((double)kk)/(kk-1)) * qf[i][0];
    }

    /* Inductively compute qf[i,j] and update t[i]. */
    for (j=1; j<TSERIESLENGTH; j++) {
        for (i=0; i<N; i++) {
	    qf[i][j] = 0;
            for (h=0; h<outlinks[i]; h++) {
                r = linkplace[i][h];
                qf[i][j] = qf[i][j] + ((double)qf[r][j-1]) / outlinks[r];
	    }
	    t[i] = t[i] - (((double)kk)/pow(-(kk-1.0),j+1.0)) * qf[i][j];
	}
    }

    /* Check if all positive. */
    allpos = 1;
    for (i=0; i<N; i++) {
	if (t[i] <= 0) {
	    allpos = 0;
        }
    }
}


/* SEEDRND. */
seedrnd()
{
    int seed;
    struct timeval tmptv;
    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * ifloor ( ( (double) tmptv.tv_usec )
                                  / ( (double) 1000000.0 ) ) );
    srand48(seed);
}


/* IFLOOR */
ifloor(double x)  /* returns floor(x) as an integer */
{
    return((int)floor(x));
}


addlink(int ii, int jj)
{

    if( (outlinks[ii] >= MAXOUTLINKS) ||
        (outlinks[jj] >= MAXOUTLINKS) ) {
	fprintf(stderr, "Too many outlinks from link %d or %d.\n", ii, jj);
	exit(1);
    }

    linkplace[ii][outlinks[ii]] = jj;
#ifdef VERBOSE
  printf("ii=%d, jj=%d, outlinks[ii]=%d, outlinks[jj]=%d ... ",
				ii,jj,outlinks[ii],outlinks[jj]);
  printf("linkplace[%d][%d] = %d\n", ii, outlinks[ii], jj);
#endif
    linkplace[jj][outlinks[jj]] = ii;
    outlinks[ii]++;
    outlinks[jj]++;
}


double absval(double x)
{
    if (x<0)
	return(-x);
    return(x);
}