/*

GALINEX.C: A computer program for running some of the simulations in
the chapter "Optimal Proposal Distributions and Adaptive MCMC", by
J.S. Rosenthal, available from:

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

Copyright (c) 2008 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).

*/

#define TARGSCALING true

#define N 10
#define sigma 0.7
#define numreps 100000

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

#define PI 3.1415926536

double M[N][N], Sigmainv[N][N];

main()
{

int i,j,k,n;
double x[N], y[N], normvec[N], sq(), normal(), logpi(), h();
int numacc;
double hsum;

#ifdef TARGSCALING
  printf("Yes targscaling ... sigma=%f\n", sigma);
#else
  printf("No targscaling ... sigma=%f\n", sigma);
#endif

  /* Define the matrix M. */
  printf("M =\n");
  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      if (i==j)
#ifdef TARGSCALING
        M[i][j] = 1.0 / (i+1);
#else
        M[i][j] = 1.0;
#endif
      else
#ifdef TARGSCALING
        M[i][j] = 0.0;
#else
        M[i][j] = (1.0*(i+1)*(j+1))/N/N;
#endif
      printf(" %f ", M[i][j]);
    }
  printf("\n");
  }
printf("\n");

  /* Let Sigmainv = M^2. */
  printf("Sigmainv =\n");
  for (i=0; i<N; i++) {
    for (j=0; j<N; j++) {
      Sigmainv[i][j] = 0.0;
      for (k=0; k<N; k++) {
        Sigmainv[i][j] = Sigmainv[i][j] + M[i][k]*M[k][j];
      }
      printf(" %f ", Sigmainv[i][j]);
    }
    printf("\n");
  }
  printf("\n");

  /* Initialise values. */
  x[0] = 1.0;
  for (i=1; i<N; i++)
    x[i] = 0.0;
  hsum = 0.0;
  numacc = 0;
  seedrnd();

for (n=0; n<numreps; n++) {
  for (i=0; i<N; i++)
    /* Choose the proposal. */
#ifdef TARGSCALING
    y[i] = x[i] + sigma * (i+1) * normal();
#else
    y[i] = x[i] + sigma * normal();
#endif
  if ( log(drand48()) < logpi(y) - logpi(x) ) {
    /* Accept the proposal. */
    for (i=0; i<N; i++)
      x[i] = y[i];
    numacc++;
  }
  hsum = hsum + h(x);
/*
  printf("x = (");
  for (i=0; i<N; i++) {
    printf(" %f ", x[i]);
  }
  printf(")\n");
  printf("h(x) = %f\n", h(x));
*/
}

/* printf("hsum = %f\n", hsum); */
printf("acceptance rate = %f\n", 1.0*numacc/numreps);
printf("average functional value = %f\n", hsum / numreps);

}


double logpi(double z[]) {
  int i,j;
  double logval = 0.0;
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
      logval = logval - z[i]*Sigmainv[i][j]*z[j]/2;
  return(logval);
}

double h(double z[]) {
  double sq();
  /* FIRST EXAMPLE: return(sq(z[0])); */
  return(sq(z[9]));
}


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


/* NORMAL:  return a standard normal random number. */
double normal()
{
    double R, theta, drand48();

    R = - log(drand48());
    theta = 2 * PI * drand48();

    return( sqrt(2*R) * cos(theta));
}


/* SEEDRND. */
int 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);
    return(0);
}


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