/**------------------------------------------------------**

   iscouple.c: sample the Ising model by coupling

   Copyright (C) 2013 by Kevin K Lin <klin@math.arizona.edu>

   This program is free software; you can redistribute it
   and/or modify it under the terms of the GNU General
   Public License as published by the Free Software
   Foundation; either version 2 of the License, or (at your
   option) any later version.

   This program is distributed in the hope that it will be
   useful, but WITHOUT ANY WARRANTY; without even the
   implied warranty of MERCHANTABILITY or FITNESS FOR A
   PARTICULAR PURPOSE.  See the GNU General Public License
   for more details.

   You should have received a copy of the GNU General Public
   License along with this program; if not, write to the
   Free Software Foundation, Inc., 51 Franklin Street, Fifth
   Floor, Boston, MA 02110-1301 USA.

   INSTRUCTIONS

   You should be able to compile and run this on any flavor
   of Un*x with GCC.  On a Mac, do

   % gcc -O3 -ffast-math -fsched-interblock -fstrict-aliasing -funroll-loops -fPIC -DNDEBUG -faltivec -fnested-functions -o iscouple.x iscouple.c -lm

   On Linux, you may have to remove some or all of the flags
   to compile.
  
   If it compiles, doing say

   % ./iscouple.x 10 10 1.2 0.1

   will run the program to sample the 10x10 ferromagnetic
   Ising model with periodic boundaries, at 1.2 times the
   critical temperature, with an external field of h=0.1.

   IMPORTANT NOTE ABOUT RANDOM NUMBERS

   The pseudo-random number generator on your system may
   differ from mine.  You may need to adjust the constant
   RNGMAX accordingly to make sure you get random samples
   from the unit interval.

 **------------------------------------------------------**/

/** libraries **/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/** you may ignore all the code bracketed by USEKLIBS **/
#ifdef USEKLIBS
#include <solartime.h>
#endif

/** system parameters **/
int m,n;
double T, h;
double beta;
const double Tc=2.269185314213022;
const double RNGMAX=2147483647.0;
const double report_interval=60.0; /* sec */

/** functions defined in this file **/
double ising_neighbor_sum(int S[m][n], int i, int j);
int ising_update(int S[m][n], int I, int J, double U);
void run_ising(int m, int n, double T, double h);
void print_ising_state(FILE *fout, int S[m][n]);
void save_ising_state_for_matlab(FILE *fout, int S[m][n]);
int same_states(int S1[m][n], int S2[m][n]);
int mod(int i, int m);


/** main() just parses arguments **/
int main(int argc, char *argv[])
{
  if ( argc != 5 ) {
    fprintf(stderr, "usage: %s <m> <n> <T/Tc> <h>\n", argv[0]);
    exit(-1);
  }
  else {
    
    /** grab parameters **/
    sscanf(argv[1], "%d", &m);
    sscanf(argv[2], "%d", &n);

    sscanf(argv[3], "%lf", &T);
    T *= Tc;
    beta = 1.0/T;
    
    sscanf(argv[4], "%lf", &h);

    printf("%%%% size=%dx%d\n", m, n);
    printf("%%%% temperature=%g (beta=%g)\n", T, beta);
    printf("%%%% h=%g\n", h);

    /** initialize rng, then run **/
    srandomdev();
    run_ising(m, n, T, h);
  }

  return 0;
}


/**------------------------------------------------------**
 ** real work starts here
 **------------------------------------------------------**/

void run_ising(int m, int n, double T, double h)
{
  /** initialize state arrays **/
  int Smin[m][n], Smax[m][n];

  { int i,j;
    for ( i=0; i < m; i++ ) {
      for ( j=0; j < n; j++ ) {
	Smin[i][j] = -1;
	Smax[i][j] = +1;
      }
    }
  }

  /** start running **/
  { unsigned long int count;
    /** diff = sum_ij (Smax[i][j] - Smin[i][j]) **/
    int diff=2*m*n;
# ifdef USEKLIBS
    double last_time=seconds()-2*report_interval, start=seconds();
# endif

    void report(void)
    {
      printf("%%%% count=%g diff=%d\n", (double)count, diff);
      fflush(stdout);
    }

    
    for ( count=0; diff>0; count++ ) {

      /** progress report **/
#ifdef USEKLIBS
      if ( seconds() > last_time + report_interval ) {
	report();
	last_time = seconds();
      }
#endif
      
      /** random() gives a random integer in {0, ..., RNGMAX} **/
      { int I = random()%m;
	int J = random()%n;
	double U = random() / RNGMAX;

	int dmin=ising_update(Smin, I, J, U);
	int dmax=ising_update(Smax, I, J, U);

	diff += (dmax-dmin);
      }
    }

    /** sanity check **/
    if ( same_states(Smax,Smin)) {
      printf("%%%% states are equal\n");
    }
    else {
      fprintf(stderr, "%%%% unequal states:\n");
      print_ising_state(stdout, Smax);
      print_ising_state(stdout, Smin);
      exit(-2);
    }

    report();
    print_ising_state(stdout, Smax);
    save_ising_state_for_matlab(stdout, Smax);
#ifdef USEKLIBS
    printf("%%%% running time: %.3lf sec\n", seconds()-start);
#endif
  }
}

int ising_update(int S[m][n], int I, int J, double U)
{
  double p = exp(beta * ising_neighbor_sum(S, I, J));
  int oldspin = S[I][J];

  if ( (p+1.0/p)*U <= p) {
    S[I][J] = +1;
  }
  else {
    S[I][J] = -1;
  }

  return S[I][J] - oldspin;
}

double ising_neighbor_sum(int S[m][n], int i, int j)
{
  return ( S[mod(i-1,m)][j] + S[mod(i+1,m)][j] +
	   S[i][mod(j-1,n)] + S[i][mod(j+1,n)] + h );
}

int mod(int i, int m)
{
  if ( i < 0 ) {
    return i+m;
  }
  else if ( i >= m ) {
    return i-m;
  }
  else {
    return i;
  }
}

int same_states(int S1[m][n], int S2[m][n])
{
  int i,j;

  for ( i=0; i < m; i++ ) {
    for ( j=0; j < n; j++ ) {
      if ( S1[i][j] != S2[i][j] ) {
	return 0;
      }
    }
  }
  
  return 1;
}

/** print ising state as ascii **/
void print_ising_state(FILE *fout, int S[m][n])
{
  int i,j;

  for ( i=0; i < m; i++ ) {
    fprintf(fout, "%%");
    for ( j=0; j < n; j++ ) {
      if ( S[i][j] > 0 ) {
	fprintf(fout, "+");
      }
      else {
	fprintf(fout, "-");
      }
    }
    fprintf(fout, "\n");
  }
}

/** save ising state as matlab matrix **/
void save_ising_state_for_matlab(FILE *fout, int S[m][n])
{
  int i,j;

  fprintf(fout, "is%dx%d=[", m, n);
  
  for ( i=0; i < m; i++ ) {
    for ( j=0; j < n; j++ ) {
      if ( S[i][j] > 0 ) {
	fprintf(fout, " +1");
      }
      else {
	fprintf(fout, " -1");
      }
    }
    fprintf(fout, "\n");
  }

  fprintf(fout, "];\n");
}