/* Code to generate samples from the autonormal distribution, also
 * known as a Gaussian field or free field.
 * Uses the algorithm described in Wilson 1999.
 * Makes a massless field on a 2-D toroidal grid with wraparound
 * boundary conditions.
 */

/*
Parameters:
   $\sigma =$ standard deviation of normals

Somebody else does:
   $X_0 := \Normal(\mu_0,\sigma)$

We generate random variables:
   $X := (X0 - \mu_0)$
   $Y := \exp(-(X/\sigma)^2/2) * \Uniform(0,1)$
   If ($X>0$) Then $Y := 1-Y$
   $L := -\sigma\sqrt{-2\log(Y)}$
   $R :=  \sigma\sqrt{-2\log(1-Y)}$

The random map is:
   $f(\mu) = \phi(\mu,(L,R,\mu_0,X))
          = \lfloor (\mu-\mu_0+R-X)/(R-L)\rfloor (R-L) +\mu_0+X$
*/


#include <math.h>
#include <stdlib.h>
#define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
#define CALLOC(NUMBER,TYPE) ((TYPE*)Calloc((NUMBER),sizeof(TYPE)))
#define REALLOC(PTR,NUMBER,TYPE) ((TYPE*)Realloc((PTR),(NUMBER)*sizeof(TYPE)))
#define FREE(PTR) free(PTR)
void *Malloc(), *Calloc(), *Realloc(), Free();

#define pi 3.14159265358979323846264338327950288419716939937510
inline double uniform() {
  return (((((double)random()) / 2147483648e0) + (double)random()) / 2147483648e0);
}

/* Use the (basic form of) the Box-Muller (1958) to get Gaussian */
double normal () {
  return sqrt(-2*log(uniform())) * cos(2*pi*uniform());
}

/* Update a lower, middle and upper configuration of the autonormal
 * by doing a heat bath at each site in sequence.  Middle configuration
 * may be NULL, in which case it does not exist and is not updated.
 * The three (or two) configurations are updated in a coupled fashion,
 * using the layered multishift coupler for the normal distribution.
 */
HeatBathSweep (int width, int height, double force,
	       double *lower, double *middle, double *upper)
{
  int i, n, north, east, west, south;
  double x,y,l,r,mu;
  /* Four springs with given force collectively behave like 4*force */
  double sigma = sqrt(1/(4*force));
  n = width*height;
  /* site 0 tied to 0 */
  lower[0] = upper[0] = 0; if (middle) middle[0] = 0;
  /* sweep through remaining sites */
  for (i=1; i<n; ++i) {
    /* generate layered multishift coupler parameters */
    x = normal();
    y = exp(-x*x/2) * uniform();
    if (x<0) y=1-y;
    l = -sqrt(-2*log(1-y));
    r =  sqrt(-2*log(y));
    r *= sigma;
    l *= sigma;
    x *= sigma;
    
    /* compute the 4 neighboring sites */
    north = i - width; if (north<0) north += n;
    south = i + width; if (south>=n) south -= n;
    east = i+1;        if (east%width==0) east -= width;
    west = i-1;        if (i%width==0) west += width;

    /* For the three (or two) configurations, compute the average of
       the neighbors, and set site to be normal centered on this
       average, using the layered multishift coupler. */
    
    mu = (lower[north] + lower[east] + lower[west] + lower[south])/4;
    lower[i] = floor((mu+r-x)/(r-l))*(r-l) + x;
    if (middle) {
    mu = (middle[north] + middle[east] + middle[west] + middle[south])/4;
    middle[i]= floor((mu+r-x)/(r-l))*(r-l) + x;
    }
    mu = (upper[north] + upper[east] + upper[west] + upper[south])/4;
    upper[i] = floor((mu+r-x)/(r-l))*(r-l) + x;
  }
}

/* Generate a proposal state for a Metropolis-Hastings update,
 * using the tree method for picking the proposal.
 * Part of the independence sampler.
 */
MHproposal (int width, int height, double force,
	    double *proposal)
{
  int i, j;
  /* For the proposal, act as if springs in tree have only half strength */
  double sigma = sqrt(1/(force/2));
  proposal[0] = 0;
  for (i=1; i<width; ++i)
    proposal[i] = proposal[i-1] + sigma*normal();
  for (i=0; i<width; ++i)
    for (j=1; j<height; ++j)
      proposal[i+j*width] = proposal[i+(j-1)*width] + sigma*normal();
}

inline double square (double x) {return x*x;}

enum {total, max, tree, tt_m_ht};

/* Compute the energy of associated with a proposal configuration, either
 * the total energy,
 * the energy contained in only the tree springs,
 * the upper bound on energy of any configuration that rejects the proposal,
 * or the total energy minus half the tree energy.
 */
double Energy (int width, int height, double force,
	     double *proposal,
	     int which)
{
  int i, j, n;
  double etree, eother;
  n = width*height;
  etree = 0;
  eother = square(proposal[0] - proposal[width-1])*force/2;
  for (i=1; i<width; ++i)
    etree += square(proposal[i] - proposal[i-1])*force/2;
  for (i=0; i<width; ++i)
    eother+= square(proposal[i] - proposal[i-width+n])*force/2;
  for (i=0; i<width; ++i) {
    int i1 = (i+1) % width;
    for (j=1; j<height; ++j) {
      etree += square(proposal[i+j*width] - proposal[i+(j-1)*width])*force/2;
      eother+= square(proposal[i+j*width] - proposal[i1+j*width])*force/2;
    }
  }
  switch (which) {
  case max: return 2*eother+etree;
  case tt_m_ht: return eother+etree/2;
  case total: return eother+etree;
  case tree:  return etree;
  default: exit(-1);
  }
}

/* Given the upper bound on the energy of any configuration that
 * rejects the proposal, determine upper and lower bounds on the
 * configuration after the Metropolis-Hastings update.
 */ 
ExtremalStates (int width, int height, double force,
		double *lower,
		double *upper,
		double emax)
{
  int i, j;
  for (i=0; i<width; ++i) {
    for (j=0; j<height; ++j) {
      upper[i+j*width] = sqrt(2*emax*(i+j)/force);
      lower[i+j*width] = -upper[i+j*width];
    }
  }
}

/* Return the gap between upper and lower configurations.
 * For diagnostic purposes only.
 */
double Gap (int width, int height, double *lower, double *upper)
{
  double gap=0;
  int i, j;
  for (i=0; i<width; ++i) {
    for (j=0; j<height; ++j) {
      gap += upper[i+j*width] - lower[i+j*width];
    }
  }
  return gap;
}

#include <string.h>
#include <stdio.h>
#include <time.h>
/* Apply a random composite map to a given state, and output whether
 * or not the composite map is officially coalescent.
 * First do a timing experiment.
 * Then update the given state using the independence sampler,
 * while also obtaining upper and lower bounds on the updated value
 * of any state.
 * Then update the given state and upper and lower bounds using
 * sweeps for the heatbath, for as long as specified by the timing
 * experiment.
 */
#define Diagnose fprintf(stderr,"%d,%lf ",count,Gap(width,height,lower,upper))
ApplyCompositeMap (int width, int height, double force, double *state)
{
  double *proposal, *upper, *lower;
  double emax;
  int count;
  int flag;
  time_t last, current;
  proposal = CALLOC(width*height,double);
  lower    = CALLOC(width*height,double);
  upper    = CALLOC(width*height,double);

  /* Timing experiment */
  MHproposal(width,height,force,proposal);
  emax = Energy(width,height,force,proposal,max);
  ExtremalStates(width, height, force, lower, upper, emax);
  last = time(0);
  for (count=0,Diagnose; memcmp(lower,upper,width*height*sizeof(double)); ++count) {
    HeatBathSweep(width,height,force, lower, NULL, upper);
    if ((current=time(0))-last > 1) {
      Diagnose;
      fflush(stderr);
      last=current;
    }
  }
  Diagnose;
  fprintf(stderr,"\n"); fflush(stderr);

  /* Actual updates */
  MHproposal(width,height,force,proposal);
  /* Maybe the state gets changed to the proposal. */
  if (uniform() < exp(Energy(width,height,force,proposal,tt_m_ht) -
		      Energy(width,height,force,state,tt_m_ht)))
    memcpy(state,proposal,width*height*sizeof(double));
  emax = Energy(width,height,force,proposal,max);
  ExtremalStates(width,height,force, lower, upper, emax);
  last = time(0);
  for (Diagnose; count;  --count) {
    HeatBathSweep(width,height,force, lower, state, upper);
    if ((current=time(0))-last > 1) {
      Diagnose;
      fflush(stderr);
      last=current;
    }
  }
  Diagnose;
  fprintf(stderr,"\n"); fflush(stderr);

  flag = !(memcmp(lower,upper,width*height*sizeof(double)));
  FREE(upper); FREE(lower); FREE(proposal);
  return flag;
}

/* Do Read-Once Coupling from the Past using the composite maps.
 */
/* Structure for creating a stream of autonormal configurations
 * with given width, height, and force.
 */
typedef struct {
  int width, height;
  double force;
  double *state, *oldstate;
} AutonormalStream;

/* Initialize Read-Once Coupling from the Past.
 */
AutonormalStream AutonormalInit (int width, int height, double force)
{
  AutonormalStream str;
  str.width = width;
  str.height = height;
  str.force = force;
  str.oldstate = CALLOC(width*height,double);
  str.state = CALLOC(width*height,double);
  /* Keep going until a coalescent map is found. */
  while (!ApplyCompositeMap(width, height, force, str.state));
  return str;
}

/* Generate the next sample using Read-Once Coupling from the Past.
 * If called more than once, then previously returned values are overwritten.
 */
double *NextAutonormal (AutonormalStream str)
{
  while (1) {
    memcpy(str.oldstate,str.state,str.width*str.height*sizeof(double));
    if (ApplyCompositeMap(str.width, str.height, str.force, str.state))
      return str.oldstate;
  }
}

main (int argc, char *argv[]) {
  int width, height, num_iter, seed;
  double force;
  AutonormalStream str;
  if (argc!=6) {
    printf("%s <width> <height> <force> <iterations> <seed>\n",argv[0]);
    exit(-1);
  }
  width = atoi(argv[1]);
  height = atoi(argv[2]);
  force = atof(argv[3]);
  num_iter = atoi(argv[4]);
  seed = atoi(argv[5]);
  if (width<=0 || height<=0 || force<=0 || num_iter<=0) {
    printf("Width, height, force, and iterations must be positive.\n");
    exit(-1);
  }
  printf("Generating %d autonormal's on %dX%d grid with coupling constant %lf\nusing seed %d.\n",
	 num_iter,width,height,force,seed);
  srandom(seed);
  str = AutonormalInit(width, height, force);
  for (; num_iter; --num_iter)
    ShowAutonormal(width,height,NextAutonormal(str));
}

ShowAutonormal (int width, int height, double *state)
{
  int i,j;
  for (j=0; j<height; ++j) {
    for (i=0; i<width; ++i) {
      printf("%+lf ",state[i+j*width]);
    }
    printf("\n");
  }
  printf("\n");
}
  
test (int argc, char *argv[]) {
  double a, b, c;
  a = atof(argv[1]);
  b = atof(argv[2]);
  c = a/b;
  printf("%s / %s\n",argv[1],argv[2]);
  printf("%lf / %lf = %lf\n",a,b,c);
}

void *Malloc (n)
     int n;
{
  void *ptr;
/*  printf("Request for %d bytes\n",n);*/
  ptr = (void*)malloc(n);
/*  printf("Ptr = %d\n",ptr);*/
  if (ptr) return ptr;
  fprintf(stderr,"Memory shortage!!  Allocating %d bytes.\n",n);
  exit(-1);
}

void *Calloc (n,elsize)
     int n,elsize;
{
  void *ptr;
/*  printf("Request for %d bytes\n",n,elsize);*/
  ptr = (void*)calloc(n,elsize);
/*  printf("Ptr = %d\n",ptr);*/
  if (ptr) return ptr;
  fprintf(stderr,"Memory shortage!!  Callocing %d bytes.\n",n*elsize);
  exit(-1);
}

void *Realloc (ptr,n)
     void *ptr;
     int n;
{
  ptr = (void*)realloc(ptr,n);
  if (ptr) return ptr;
  fprintf(stderr,"Memory shortage!!  Reallocating %d bytes.\n",n);
  exit(-1);
}