/* Metropolis-within-Gibbs algorithm for Guadeloupe sugar cane data */

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

#define MAXN 1742
#define K 30
#define infinity 999999999.0
#define PI 3.1415926536

double drand48();
int L[MAXN], U[MAXN];
double lambda[MAXN][K], D[MAXN][MAXN], sq();
int N;

/* BEGIN MAIN PROGRAM. */

int main(int argc, char **argv)
{

/* Initial declarations. */
FILE *fp;
int M, B, i, j, k, t;
char tmpstring[20];
double x[MAXN], y[MAXN];
int S6[MAXN], S10[MAXN], S14[MAXN], S19[MAXN], S23[MAXN], S30[MAXN];
double uniform(), exponential(), normal(), logpi();
void seedrand();
double curtheta, curmu, curlogpi, newlogpi, newmu, newtheta, summu, sumtheta;
int curtau[MAXN], newtau[MAXN], sumtau[MAXN];
int propmu, propth, proptau, accmu, accth, acctau;

/* Seed the random number generator. */
seedrand();

/* Read in the data. */
printf("Reading data ...\n");
if ((fp = fopen("canedata","r")) == NULL) {
    fprintf(stderr, "Unable to read file 'canedata'.\n");
    exit(1);
}
N = 0;
for (i=0; i<MAXN; i++) {
    fscanf(fp, "%s", &tmpstring);
    x[N] = atof(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    y[N] = atof(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S6[N] = atoi(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S10[N] = atoi(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S14[N] = atoi(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S19[N] = atoi(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S23[N] = atoi(tmpstring);
    fscanf(fp, "%s", &tmpstring);
    S30[N] = atoi(tmpstring);
    /* Do "cheat" of restricting to 10 x 10 grid. */
    if ( (x[N] < 10.0) && (y[N] < 10.0) ) {
      /* printf("Keeping site %d.\n", N); */
      N++;
    }
}
fclose(fp);
printf("Number of sites studied: %d\n", N);

/* OUTPUT SOME TEST VALUES -- NO.
printf("%f  %f  %d  %d  %d\n", x[2], y[2], S6[2], S10[2], S14[2]);
for (i=0; i<N; i++)
  printf("i=%d: x=%f, y=%f\n", i, x[i], y[i]);
*/

/* Determine the L_x and U_x values, etc. */
printf("Computing infection time ranges ...\n");
for (i=0; i<N; i++) {
  if (S6[i]) {
    L[i] = 0;
    U[i] = 6;
  } else if (S10[i]) {
    L[i] = 6;
    U[i] = 10;
  } else if (S14[i]) {
    L[i] = 10;
    U[i] = 14;
  } else if (S19[i]) {
    L[i] = 14;
    U[i] = 19;
  } else if (S23[i]) {
    L[i] = 19;
    U[i] = 23;
  } else if (S30[i]) {
    L[i] = 23;
    U[i] = 30;
  } else {
    L[i] = 30;
    U[i] = K+10;
  }
}

/* Compute the cane-cane distances. */
printf("Computing pairwise distances ... ");
for (i=0; i<N; i++)
  for (j=0; j<N; j++)
    D[i][j] = sqrt( sq(x[i]-x[j]) + sq(y[i]-y[j]) );
printf("done.\n");
printf("L[0]=%d\n", L[0]);

/* Test the logpi function ... */
for (i=0; i<N; i++)
  curtau[i] = U[i];
printf("Test logpi value: %f\n", logpi(2.0,3.0,curtau));

/* Initialise the Markov chain values. */
M = 110000;  /* run length */
B = 10000;  /* burn-in */
M = 11000;
B = 1000;
/* NO, TOO UNSTABLE:
curtheta = 10 * exponential();
curmu = 10 * normal();
*/
curtheta = 2.0;
curmu = 3.0;
sumtheta = summu = 0.0;
for (i=0; i<N; i++) {
  curtau[i] = iround((L[i]+U[i])/2.0);
  sumtau[i] = 0;
}
propth = propmu = proptau = accth = accmu = acctau = 0;
curlogpi = logpi(curtheta,curmu,curtau);
printf("Initial logpi value: %f\n", curlogpi);

/* Run the Markov chain! */
printf("Running chain, for %d iterations (burn-in %d) ...\n", M, B);
for (t=1; t<=M; t++) {
  printf(" t=%d", t);
  fflush(stdout);

  /* Propose new theta value. */
  propth++;
  newtheta = curtheta + 2.0 * normal();
  newlogpi = logpi(newtheta,curmu,curtau);
  if ( log(uniform()) < newlogpi - curlogpi ) {
    /* Accept proposal. */
    accth++;
    curtheta = newtheta;
    curlogpi = newlogpi;
  }

  /* Propose new mu value. */
  propmu++;
  newmu = curmu + 0.1 * normal();
  newlogpi = logpi(curtheta,newmu,curtau);
  if ( log(uniform()) < newlogpi - curlogpi ) {
    /* Accept proposal. */
    accmu++;
    curmu = newmu;
    curlogpi = newlogpi;
  }

  /* Loop through the cane sites. */
  for (i=0; i<N; i++) {
    /* printf(" i=%d", i);
    fflush(stdout); */
    if (curtau[i] <= K) {
      /* Propose new tau[i] value. */
      proptau++;
      for (j=0; j<N; j++) {
        if (j==i)
          newtau[j] = curtau[j] + pmo();
        else
          newtau[j] = curtau[j];
      }
      newlogpi = logpi(curtheta,curmu,newtau);
/* printf("t=%d, i=%d; logpi diff = %f\n", t, i, newlogpi - curlogpi); */
      if ( log(uniform()) < newlogpi - curlogpi ) {
        /* Accept proposal. */
	acctau++;
/* printf("Accepted tau proposal at i=%d; acctau=%d.\n", i, acctau); */
        curtau[i] = newtau[i];
        curlogpi = newlogpi;
      }
    }
  }

  /* Update our running sums. */
  if (t > B) {
    sumtheta = sumtheta + curtheta;
    summu = summu + curmu;
    for (j=0; j<N; j++)
      sumtau[j] = sumtau[j] + curtau[j];
  }

} /* End of Markov chain run. */

/* Output the results. */
printf("\n\nARtheta=%f,  ARmu=%f,  ARtau=%f\n",
  ((double)accth)/propth, ((double)accmu)/propmu, ((double)acctau)/proptau );
printf("Mean theta: %f\n", sumtheta/(M-B));
printf("Mean mu: %f\n", summu/(M-B));
if ((fp = fopen("tauvals","w")) == NULL) {
    fprintf(stderr, "Unable to write file 'tauvals'.\n");
    exit(1);
}
for (i=0; i<N; i++)
  fprintf(fp, "%f  %f  %f\n", x[i], y[i], ((double)sumtau[i])/(M-B));
fclose(fp);

return(0);

}  /* End of Main Program. */


/* pmo: function which returns +1 or -1, with probability 1/2 each. */
int pmo() {
  if (drand48() < 0.5)
    return(+1);
  return(-1);
}

/* Specify the target log density. */
double logpi(double thetheta, double themu, int thetau[]) {
  int ii, jj, kk;
  double tmpsum;
  /* Compute the lambda_{x,k} values. */
  for (ii=0; ii<N; ii++) {
    if ( (thetau[ii] <= L[ii]) || (thetau[ii] > U[ii]) ) {
      /* printf("out of range: ii=%d, L=%d, tau=%d, U=%d\n",
      					ii, L[ii], thetau[ii], U[ii]); */
      return(-infinity);
    }
    for (kk=0; kk<K; kk++) {
      tmpsum = themu;
      for (jj=0; jj<N; jj++) {
        if (thetau[jj] <= kk-1)
	  tmpsum = tmpsum + exp(-thetheta*D[ii][jj]);
      }
      lambda[ii][kk] = exp(tmpsum);
    }
  }
  /* Compute the sums. */
  tmpsum = -thetheta/10 - sq(themu)/200;
  for (ii=0; ii<N; ii++) {
    if (thetau[ii] <= K)
      tmpsum = tmpsum + log(1-exp(-lambda[ii][thetau[ii]]));
    for (kk=0; kk < imin(thetau[ii]-1,K); kk++)
      tmpsum = tmpsum - lambda[ii][kk];
  }
  return(tmpsum);
}


/* SEEDRAND: SEED RANDOM NUMBER GENERATOR. */
void seedrand()
{
    int i, seed;
    struct timeval tmptv;
    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) tmptv.tv_usec;
    srand48(seed);
    (void)drand48();  /* Spin it once. */
}


double sq(double xxx)
{
  return(xxx*xxx);
}


double uniform()
{
    double drand48();
    return( drand48() );
}

double exponential()
{
    double uniform();
    return( -log(uniform()) );
}

/* NORMAL:  return a standard normal random number. */
double normal()
{
    double RRR, ttt, uniform();

    RRR = - log(uniform());
    ttt = 2 * PI * uniform();

    return( sqrt(2*RRR) * cos(ttt));
}

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

int iround(double xxx)
{
    return( ifloor(xxx+0.5) );
}

int imin(int iii, int jjj)
{
    if (iii < jjj)
      return(iii);
    return(jjj);
}