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

#define idx(x, y, N) ((((y)+N)%N)*N+((x)+N)%N) // conversion between x, y and index i
int neigh[4][2] = {{1,0}, {0,1}, {-1, 0}, {0, -1}}; //neighbour distance in a square lattice
#define b2s(b, r) (2*(((b)>>(r))&1)-1) // takes bit r of b and converts it to a spin (-1,1)

typedef struct _lattice {
  int N; 
  int R; // number of replicas
  int * s; // spins (multispin) 
  double ** c; // coulpings
} lattice;

lattice * new_lattice(int N, int R, double d, double a); // degree of disorder and asymmetry
void evolve(lattice * l, double K, double H);
double overlap (lattice * l);
int bitcount (unsigned x);
void visualize(FILE * gp, lattice *l, unsigned mask, char * title);
void visualize_connections(FILE * gp, lattice *l);

#define VIS

int main() {
  lattice * l;
  int t, T, TRANS;
  double o;
  double a; // asymmetry
  double d; // disorder (anti-ferromagnetism)
  int N;
  int R; // number of replicas
  double H;
  double K; // iverse of temperature
  char title[100];

  srand48(time(NULL));
  
  T = 100;
  TRANS = 20;
  R = 4;
  a = 0.2;
  d = 1;
  N = 10;
  l = new_lattice(N, R, d, a); 
  H=3; // magnetic field for TRANS
  K = 100;
  
  

#ifdef VIS
  FILE *gp;
  gp = popen("gnuplot", "w");
  fprintf(gp, "unset label; set size square; \n");
  fprintf(gp, "set xrange  [-0.5:%f]\n", N-0.5);
  fprintf(gp, "set yrange  [-0.5:%f]\n", N-0.5);
  fprintf(gp, "set cbrange  [0:%d]\n", (1<< l->R) -1);
  visualize_connections(gp, l);

  visualize(gp, l, 3, "time=0");
  printf("premi return per cominciare\n");
  getchar();
#endif  
  
  
  

  for (t=0; t<T+TRANS; t++) {
    evolve(l, K, H);
    if (t>TRANS) {
      H=0;
      o = overlap(l);
      printf("%d %f\n", t-TRANS, o);
    }
    sprintf(title, "time=%d", t);
#ifdef VIS
    visualize(gp, l, 3, title);
#endif
  }
}

void evolve(lattice *l, double K, double H) {
  int x,y;
  int i,j,r;
  double h, deltaE;
  int tmp;
    
  for (r=0; r< l->R; r++) {// REPLICAS
    for (i=0; i<l->N*l->N; i++) { // step Monte-Carlo
      x = (int) (drand48() * l->N);
      y = (int) (drand48() * l->N);
      h=H; // here one could put the local magnetization of replicas....
      for (j=0; j<4; j++) {
        tmp = l->s[idx(x+neigh[j][0],y+neigh[j][1],l->N)];
        h += (int) b2s(tmp, r) * l->c[idx(x,y, l->N)][j];
      }
      deltaE = 2 * b2s(l->s[idx(x,y,l->N)],r) * h;
      if (deltaE <= 0 || drand48() < exp(-deltaE*K)) {
        l->s[idx(x,y,l->N)] ^= 1<<r; // flips bit 
      }
    }
  }
}
      

double overlap (lattice * l) {
  int couples;
  int i, n, R;
  double ov;
  
  couples = 0;
  R = l->R;
  for (i=0; i<l->N*l->N; i++) {
    n = bitcount(l->s[i]);
    couples += (n*(n-1) + (R-n)*(R-n-1))/2;
  }
  ov = (double) couples/(l->N*l->N)/(R*(R-1)/2);
  return(ov);
  
}


int bitcount (unsigned x) {
  int tot;
  
  tot = 0;
  while(x) {
    tot += x & 1;
    x >>= 1;
  }
  return tot;
}

lattice * new_lattice(int N, int R, double d, double a) {
  lattice * l;
  int i,j;
  int x, y;
  
  l = calloc(1, sizeof(lattice));
  l->N = N;
  l->R = R;
  l->s = calloc(N*N, sizeof(int));
  for (i=0; i<N*N; i++) {
    for (j=0; j<R; j++) {
      l->s[i] |= (drand48() < 0.5)<<j;
    }
  }
  l->c = calloc(N*N, sizeof(double *));
  for (i=0; i<N*N; i++) {
      l->c[i] = calloc(4, sizeof(double));
  }
  // first the east and north bonds
  for (y=0; y<N; y++) {
    for (x=0; x<N; x++) {
      for (j=0; j<2; j++) {
        l->c[idx(x,y,l->N)][j] = 2*(drand48() < d) - 1;  // d=0.5 => 50% +1 50% -1, d=1 all +1 , d=0 all -1 
      }
    }
  }
  // then west and south
  for (y=0; y<N; y++) {
    for (x=0; x<N; x++) {
      for (j=0; j<2; j++) { // a=0 : symmetric, a=1 : antisymmetric
        l->c[idx(x,y,l->N)][j+2] = (1-2*(drand48() < a)) * l->c[idx(x+neigh[j+2][0],y+neigh[j+2][1],l->N)][j];
     }
    }
  }
  return l;
}

void visualize(FILE * gp, lattice *l, unsigned mask, char * title) {
  int data;
  
#if 0
  fprintf(gp, "set size square; unset key; set cbrange [0:%d]\n", mask);
  fprintf(gp, "plot '-' binary array=%dx%d w image \n", l->N,l->N);
  for (int i=0; i<l->N*l->N; i++) {
    data = l->s[i] & mask;
    fwrite(&data, sizeof(int),1, gp);
    printf("%d ", data);
  }
#endif

  fprintf(gp, "set title '%s'\n", title);
  fprintf(gp, "plot '-' matrix w image\n"); 
  for (int y=0; y<l->N; y++) {
    for (int x=0; x<l->N; x++) {
      data = l->s[idx(x,y,l->N)] & mask;
      fprintf(gp, "%d ", data);
    }
    fprintf(gp, "\n");
  }
  fprintf(gp,"e\n");

  fflush(gp);
  
#ifdef VISTEXT
  // textual
  for (int y=l->N-1; y>=0; y--) {
    for (int x=0; x<l->N; x++) {
      data = l->s[idx(x,y,l->N)] & mask;
      printf("%d ", data);
    }
    printf("\n");
  }
#endif

}

void visualize_connections(FILE * gp, lattice *l) {
  fprintf(gp, "unset label; \n");
  for (int x=0; x<l->N; x++) {
    for (int y=0; y<l->N; y++) {
      for (int j=0; j<4; j++) {
        fprintf(gp, "set label front '%c' at %f,%f textcolor 'blue' center\n", l->c[idx(x,y,l->N)][j]>0? '+':'-', 
          x+(double)neigh[j][0]/4, y+(double)neigh[j][1]/4);
      }
    }
  }

  fflush(gp);
}