/* This is a sample C program to perform Wang-Landau sampling
 * for the two-dimensional nearest-neighbor Ising model on a 
 * square lattice with periodic boundary conditions (no external 
 * magnetic field). This program was written as part of a 
 * graduate teaching exercise and is not intended to be in any 
 * way optimized. The random number generator used in this 
 * program is a simple congruential one that is part of the
 * standard C library. High resolution simulations often require 
 * better generators that can produce random numbers with better 
 * statistical properties and longer sequences than generated by 
 * a congruential method. 
 *
 * Reference paper: " A new approach to Monte Carlo simulations 
 * in statistical physics: Wang-Landau sampling", D. P. Landau,
 * Shan-Ho Tsai, and M. Exler, Amer. J. Phys. 2004.
 *
 */

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

/* Macro for transformation energy level number -> energy */
#define energy(b) (4.0*(b)-2.0*L*L)

/* Macro for reading parameters */
#define MAX_LEN 256
#define read_variable(string,type,variable)   \
 {                       \
  char input[MAX_LEN];            \
  fputs(string,stdout);            \
  fgets(input,MAX_LEN,stdin);         \
  sscanf(input,type,&variable);        \
 }

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

 int ** s;  /* Field of Ising spins (dynamic memory) */
 int L;   /* System size (L by L sites) */
 int top_b; /* Number of energy levels */
 double flat_thres;  /* Threshold for flat histogram */
 double * g; /* Array for g-values (density of states) */
 int * hist; /* Array for histogram */
 int b;   /* Current energy level number */
 int b_new; /* Energy Level number of proposed move */
 int seed;  /* Random number seed */
 int skip;  /* Number of moves to skip between checks of histogram */
 int min_steps;  /* Minimum number of moves for each f */
 double f;  /* Factor to multiply g(E) after move is accepted */
 double min_f;  /* f to start with */
 double dec_pow; /* Exponent of power law to decrease f after each run */
 double prob;/* Probability for acceptance of move */
 int i,j,n,mc_steps; /* Counters */
 char filename[MAX_LEN];
 FILE * out_file;

 /* Read system size */
 read_variable("System size: ","%d",L);
 if(L<2) {
  fputs("System size has to be larger than 1\n",stderr);
  exit(1);
 }

 /* Read starting f */
 read_variable("Starting f: ","%lf",f);
 if(!(f>1.0)) {
  fputs("Starting f has to be larger than 1\n",stderr);
  exit(1);
 }

 /* Read minimum f */
 read_variable("Minimum f: 1+","%lf",min_f);
 min_f+=1.0;
 if( !( (min_f<f) && (min_f>1.0))) {
  fputs("Minimum f has to be smaller than starting f\n",stderr);
  exit(1);
 }

 /* Read exponent for decrease of f */
 read_variable("Exponent to decrease f: 1/","%lf",dec_pow);
 dec_pow=1.0/dec_pow;
 if( !( (dec_pow<1.0) && (dec_pow>0.0))) {
  fputs("Decrease exponent has to be in (0.0,1.0)\n",stderr);
  exit(1);
 }

 /* Read criterion for flat histogram */
 read_variable("Criterion for flat histogram: ","%lf",flat_thres);
 if( !( (flat_thres<1.0) && (flat_thres>0.0))) {
  fputs("Fraction has to be in (0.0,1.0)\n",stderr);
  exit(1);
 }

 /* Read minimum number of MC steps before histogram is checked */
 read_variable("Minimum number of MC steps for each f: ","%d",min_steps);
 if( !(min_steps>0)) {
  fputs("Number of steps has to be positive\n",stderr);
  exit(1);
 }

 /* Read number of MC steps between subsequent histogram checks */
 read_variable("Number of MC steps between histogram checks: ","%d",skip);
 if( !(skip>=0)) {
  fputs("Number of steps has to be non-negative\n",stderr);
  exit(1);
 }

 /* Allocate memory */
 g=(double *)calloc(L*L+1,sizeof(double));
 hist=(int *)calloc(L*L+1,sizeof(int));
 s=(int **)calloc(L,sizeof(int *));

 if(!(g&&hist&&s)) {
  fputs("Error allocating memory\n",stderr);
  exit(1);
 }

 for(i=0;i<L;i++) {
  s[i]=(int *)calloc(L,sizeof(int));
  if(!s[i]) {
   fputs("Error allocating memory\n",stderr);
   exit(1);
  }
 }

 /* Read in initial random number generator seed */
 read_variable("Random number seed: ","%d",seed);
 srand(seed);
 
 /* Get filename for output */
 filename[0]=0;
 read_variable("Output filename (leave empty for stdout): ","%s",filename[0]);
 if( strlen(filename)==0 ) {
  out_file=stdout;
 }
 else {
  if(!(out_file=fopen(filename,"w"))) {
   fputs("Error opening file for writing\n",stderr);
   exit(1);
  }
 }

 /* Initialize g distribution (logarithm of g is used!) */
 for(i=0;i<=L*L;i++) {
  g[i]=0.0;
 }

 /* Form initial state */
 for(i=0;i<L;i++)
  for(j=0;j<L;j++)
   s[i][j]=+1;

 /* set upper limit for energy level number; depends on if L is odd or even */
 if(L%2)
  top_b=L*(L-1)+1;
 else
  top_b=L*L+1;

 /* Set energy level number of initial state */
 b=0;

 /* Repeat loop until f reaches minimum value */
 while(f>min_f) {
  /* Only use logarithms of f and g for the calculation */
  double lnf=log(f);
  int c,cont=1;

  /* Start with flat histogram */
  for(i=0;i<=L*L;i++) {
   hist[i]=0;
  }

  /* Reset counter for number of steps */
  n=0;
  c=skip+1;
  mc_steps=0;

  do {
   
   /* Perform one Monte Carlo sweep: randomly pick N=L*L sites */
   for(i=0;i<L*L;i++) {
    
    int si,sj,neighsum;

    /* Choose random site for flipping */
    int site=(int)(L*L*(rand()/(RAND_MAX+1.0)));
    si=site/L;
    sj=site%L;
 
    /* Calculate sum over neighbor spins */
    neighsum= s[(si!=L-1)?si+1:0][sj];
    neighsum+=s[(si!=0)?si-1:L-1][sj];
    neighsum+=s[si][(sj!=L-1)?sj+1:0];
    neighsum+=s[si][(sj!=0)?sj-1:L-1];

    /* Number of new energy level */
    b_new=b+s[si][sj]*neighsum/2;

    /* Check if energy level is inside specified range */
    if(b_new<top_b) {

     /* Calculate probability for acceptance */
     prob=exp(g[b]-g[b_new]);

     /* Accept if proposed move has lower g (prob>=1.0) */
     /* or if random number [0.0,1.0] is less than prob */
     if((prob>=1.0) || ((double)rand()/RAND_MAX < prob)) {
      b=b_new;
      s[si][sj]*=-1;
     }

     /* Multiply g[b] by f: add logarithms */
     g[b]+=lnf;

     /* Update counters */
     ++hist[b];
     n++;
    }

   }

   /* Update counters for Monte Carlo steps */
   mc_steps++;
   c++;


   /* Check for flat histogram after skipping skip steps */
   if((mc_steps>=min_steps) && (c>=skip)) {

    int div;

    /* Reset skip counter */
    c=0;

    /* Check for flat histogram */
    /* Exclude "empty" energy levels, energy jumps by 8J from groundstate to first excited state */
    div=top_b>L*L-1?top_b-2:top_b-1;

    /* Look into each energy level and continue if fraction is larger than given threshold */
    for(i=0;
      (i<top_b) &&
       ((i==1) || (i==L*L-1) || ((double)hist[i]/n*div > flat_thres));
      i++);

    /* If loop made it to the end, histogram is flat */
    if(i==top_b) cont=0;
    else cont=1;
   }
   else cont=1;
  } while(cont);

  /* Normalize (logarithmic) g values, so g(0)=1 */
  for(i=1;i<=top_b;i++) {
   g[i]-=g[0];
  }
  g[0]=0.0;

  /* Give status message with current f */
  printf("f-1=%.3e\tlog(f)=%g\tmc_steps=%d\n",f-1.,lnf,mc_steps);

  /* Decrease factor f by a power law */
  f=pow(f,dec_pow);
 }

 /* Print g distribution and final histogram */
 /* Output of g is normalized to the number of ground states */ 
 /* We have two ground states in the Ising model, so g(0)=2 */ 
 fputs("b\tE(b)\tlog(g(E))\tH(E)\n",out_file);
 fputs("-----------------------------------------\n",out_file);
 for(i=0;i<top_b;i++) {
  if((i!=1)&&(i!=L*L-1))
   fprintf(out_file,"%d\t%+6.2f\t%g\t%d\n",
       i,energy(i),g[i]+log(2.0),hist[i]);
 }

 /* Close file if we are not using stdout */
 if(out_file!=stdout)
  fclose(out_file);

 return 0;
}