/* gcc  eoglassplot.c -o eoglassplot.x -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM  -lICE -lXext -lX11 -lm */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>
#include<sys/types.h>
#include<unistd.h>
/*#include<termios.h>*/
#include<term.h>
#include<ncurses.h>
#define min(A,B) ((A) < (B)) ? (A) : (B)
#define MAXVERTEX 10000       /* maximal number of vertices in a graph */
#define MAXLINK 20     /* maximal number of links for a vertex */
#define LENGTH 65535
#define C2 32767
#define MAXNEIGHBOR 3  /* maximal number of neighbors in hyperlink */
#define Pspin 2  

typedef unsigned short int Natural;


struct Vertex {
  short int spin;        /* value of vertex */
  Natural links;         /* # of connections */
  Natural rank;             /* position in fitness rank list */
  Natural fitness;          /* fitness of vertex */
  short int *bondcut[MAXLINK]; /* pointer to cutlist, which is +|bond| if bond 
				is satisfied, -|bond| if bond is cut */
  struct Vertex *neighbor[MAXLINK][MAXNEIGHBOR];   /* nearest neighbor connectivity vector of i-th edge to j neighbors */
  Natural index;   /* vertex' own index label, makes only sense for array */
  int updatetime;  /* time of last update */
};

/* RNG */
double rng1(long int);
/* Input and Output */
void initialize(int *,float *,float *,int *,int *,int *);
void output(Natural,float,float,int,int,int,int,int,int);
/* EO update */
void flip(struct Vertex *, int *);
/* making graphs: */
int initrandomgraph(Natural,float,int);
int initlattice(int,int,Natural *);
int intpow(int, int);
int neigh(int, int, int, int , int);
/* fixing bonds */
int initbonds(int, float , int);
void initspins(Natural, int, int, int *, int *);
void initfitness(Natural);
/* tau-distribution */
void inittau(double,double);
Natural taupick(void);
/* rank-ordering */
void initrank(Natural);
void rerank(struct Vertex *, Natural);
struct Vertex *findvertex(Natural);

/************* Plot Stuff *****************************/
void init_keyboard(void);
void close_keyboard(void);
int kbhit(void);
int readch(void);
void draw(struct Vertex *,int,int,int,float);
void plot(int,float,float,int,int,int,float,int);
void initplotwindow(void);
void closeplotwindow(void);     

static struct termios initial_settings, new_settings;
static int peek_character = -1;
static int pl_handle;
const float pl_scale = 100.0;
/************* END Plot Stuff *****************************/

struct Vertex V[MAXVERTEX];
short int cutlist[MAXLINK*MAXVERTEX];
struct Vertex *fitlist[MAXLINK][MAXVERTEX];
Natural len[MAXLINK];
Natural minfitness;
Natural pdf[LENGTH];
short int bondlist[MAXLINK*MAXVERTEX];

main(void){
  int L,d,numlines,timescale;
  int instance,maxinstance,run,maxrun,maxiter,t,mintime;
  int seed,seed1,seed2,balance,energy,energy0,minenergy,magnit;
  struct Vertex *v;
  Natural nv,rank;
  float tau,f,dx,p;
  int plotcounter=0;

  rng1((unsigned int)getpid()); /* Create a random seed from system if seeds = 0 */
  d = 2;
  initialize(&L,&f,&tau,&maxiter,&seed1,&seed2);
  numlines = initlattice(L,d,&nv);
  inittau(tau,nv);
  balance = initbonds(numlines,f,seed1);
  initspins(nv,numlines,seed2,&magnit,&energy);
      energy0 = energy;
      minenergy = energy;
      initfitness(nv);
      initrank(nv);
      initplotwindow();  
      dx = pl_scale/(float)L;
      timescale = 1+LENGTH/maxiter;/*set fade-time for plot*/
      for(t=0; t<maxiter; t++){
	rank = taupick();
        v = findvertex(rank);
	v->updatetime = t; /* remember update-time (age) of variable */
        flip(v,&energy);	
	if(plotcounter++ > L){ /* plot only every L-th update */
	  plotcounter = 0;
	  plot(L,f,tau,t,energy,minenergy,dx,timescale);
	}
        if(energy < minenergy){  /* remember best-found energy */
	  minenergy = energy;
	  mintime = t;
        }
      }
      output(nv,f,tau,numlines,energy0,energy,minenergy,maxiter,mintime);
      closeplotwindow();
}
 

/********** END of MAIN           *** *****************************/

/******************************************************************/
/*                                                                */
/*      Input and output                                          */
/*                                                                */
/******************************************************************/

/* 
*/

void initialize(int *L, float *f, float *tau, int *maxiter, int *seed1, int *seed2){
  printf("*************************************************************\n");
  printf("*                                                           *\n");
  printf("*                Extremal Optimization                      *\n");
  printf("*                      for a                                *\n");
  printf("*                 d=2 Ising Spin-Glass                      *\n");
  printf("*                   (with periodic BC)                      *\n");
  printf("*                                                           *\n");
  printf("*                                                           *\n");
  printf("* Stefan Boettcher                            May   1, 2002 *\n");
  printf("* www.physics.emory.edu/faculty/boettcher/                  *\n");
  printf("*************************************************************\n");
  printf("\n");
  printf("Color Scheme:  Black=0 cut bonds\n");
  printf("               Blue =1 cut bonds\n");
  printf("               Green=2 cut bonds\n");
  printf("               Red  >2 cut bonds\n");
  printf("Fade-to-white: Indicates time since last update!\n");
  printf("\n");
  printf("______________BEGIN Input_________________________________\n \n");
  printf("Lattice Length: L [<100] = ");
  scanf("%d",L);
  printf("Bond Balance (prob. p of +1, 1-p of -1 bonds): p+ = ");
  scanf("%f",f);
  printf("EO-parameter: tau = ");
  scanf("%f",tau);
  printf("Update Steps per Run: T = ");
  scanf("%d",maxiter);
  printf("Random Seeds: Seed(Instance) = ");
  scanf("%d",seed1);
  printf("Random Seeds: Seed(Run) = ");
  scanf("%d",seed2); 
  printf("_____________END Input____________________________________\n");
}

void output(Natural nv,float f,float tau, int numlines, int energy0, int energy, int minenergy, int t, int mintime){
  printf("\n\nFinal EO result (after %d updates): \n",t);
  printf("========================================= \n");
  printf("Vertices = %d , Edges = %d , p+ = %4.2f, tau = %4.2f \n\n",nv,numlines,f,tau);
  printf("Initial Cuts = %d , Final Cuts = %d ",energy0,energy);
  printf("BEST CUTS = %d , obtained at update step %d \n",minenergy,mintime);
  printf("*************************************************************\n");
}

/********** END of Input/Output       *****************************/


/******************************************************************/
/*                                                                */
/*      EO update procedure                                       */
/*                                                                */
/******************************************************************/

/* "flip" inverts the state of the submitted vertex v and updates the 
system energy. Then, it recalculates fitnesses and ranks for the effected
vertices, and inverts the state of the effected bonds in cutlist INDIRECTLY
via the pointers v.bondcut.
*/

 
void flip(struct Vertex *v, int *energy){
  struct Vertex *nextV;
  Natural i,j,sign;

  v->spin *= -1;
  *energy += (v->links)-2*(v->fitness);
  rerank(v,v->links - v->fitness);
  for(i = 0; i<(v->links); i++){
    *v->bondcut[i] *= -1;
    sign = *v->bondcut[i];
    for(j = 0; j<(Pspin-1); j++){
      nextV = v->neighbor[i][j];
      rerank(nextV,nextV->fitness - sign);
    }
  }
}


/********** END EO updates            *****************************/

/******************************************************************/
/*                                                                */
/*      Various ways of making graphs                             */
/*                                                                */
/******************************************************************/

/* 
"initrandomgraph" makes an Erdos-Renyi-type Random Graph with n vertices;
any pair of vertices is connected with probability p. The rng1 is initiated 
with "seed". The result is an array of struct Vertex. 
The procedure returns the number of edges in the graph.
*/
int initrandomgraph(Natural nv, float p, int seed){
  int numlines = 0;
  Natural i,j;
  rng1(seed);
  for(i = 0; i<nv; i++){
    V[i].index = i;
    V[i].links = 0;
    for(j = 0; j<i; j++){
      if(rng1(0) < p){
	V[i].bondcut[V[i].links] = &cutlist[numlines];
	V[j].bondcut[V[j].links] = &cutlist[numlines++];
	V[i].neighbor[V[i].links++][0] = &V[j];
	V[j].neighbor[V[j].links++][0] = &V[i];
      }
    }
  }
  return(numlines);
}

/* integer power of an integer (with 0^0 = 1) */
int intpow(int base, int power){
  int i,x = 1;
  for(i = 0; i<power; i++) x *= base;
  return(x);
}

/*
"initlattice" generates a graph (in form of an array of struct vertex) 
out of a "dim"-dimensional hypercubic lattice of length L with periodic 
boundary conditions. It returns the number of edges. Note that the connection 
matrix is established only in the forward (+1) direction for each vertex,
since the backward (-1) direction is the forward direction of the preceeding 
vertex, at least on a periodic lattice!
*/

int initlattice(int L, int dim, Natural *nv){
  Natural n,nextton;
  int d;
  int numlines = 0;

  *nv = (Natural) intpow(L,dim);
  for(n = 0; n<*nv; n++){
    V[n].index = n;
    V[n].links = 0;
  }
  for(n = 0; n<*nv; n++){
    for(d = 0; d<dim; d++){
      nextton = neigh(n,d,1,L,dim);
      V[n].bondcut[V[n].links] = &cutlist[numlines];
      V[nextton].bondcut[V[nextton].links] = &cutlist[numlines++];
      V[n].neighbor[V[n].links++][0] = &V[nextton];
      V[nextton].neighbor[V[nextton].links++][0] = &V[n];
    }
  }
  return(numlines);
}


/*
"neigh" returns the label of the neighbor in the d-th dimension
in "direction" = +/-1 of the vertex labeled iv for a "dim"-dimensional
hypercubic lattice of length L with periodic boundary conditions.
It is for the label  iv = x_0+x_1*L+x_2*L^2+...  where 0 <= x_d<L and
0 <= d<dim is the d-th coordinate of the vertex iv, thus x_d = [iv div L^d] mod L.
Then neighbor(iv,d,direction) = iv+[(x_d+L+direction) mod L-x_d]*L^d
*/
int neigh(int iv, int d, int direction, int L, int dim){
  int i,x,Ld;
  if(d<dim){
    Ld = intpow(L,d);
    x = (iv/Ld)%L;
    return(iv+((x+L+direction)%L - x)*Ld);
  }
  else{ 
    printf("error in neigh %d",iv);
    return(-1);
  };
}



/******************************************************************/
/*                                                                */
/*      Fixing bonds and initializing spins                       */
/*                                                                */
/* "initbonds" makes "bondlist" with +/-1, "initspins" assigns    */
/* +/-1 spins to each vertex and marks bonds either violated (cut)*/
/* or not (state -1 or +1) in array "cutlist".                    */
/* IMPORTANT: "V[n].bondcut[i]" points to respective element in   */
/* in "cutlist" (not "bondlist"!), ie. it does NOT refer to the   */
/* weight of the bond but its STATE! Thus, in an update (see      */
/* "flip") it is easier to determine all fitnesses.               */
/*                                                                */
/******************************************************************/

/*
"initbonds" assigns a random +/-1 bond to each edge in the graph, with 
probability f it's +1 and (1-f) it's -1. Returned is the imbalance in the 
weights, #(+1)-#(-1).
*/
int initbonds(int numlines, float f, int seed){
  short int bond;
  int balance = 0,l = 0;
  rng1(seed);
  while(l < numlines){
      bond = 2*(rng1(0) < f)-1;
      balance += bond;
      bondlist[l++] = bond;
  }
  return(balance);
}
      

/*
"initspins" assigns a startup spin configuration to the graph. It calculates
the  startup magnitization and energy.
The procedure is slightly tricky: First, the bondlist (which contains the
quenched bond values J_i = +/-1) is put into cutlist, then cutlist is negated 
(if V.spin = -1) through V.bondcut which points to cutlist. So, by the end, cutlist[i] = J_i*s_1*s_2*...*s_k, which is = +1, if the bond is satisfied, or = -1 if the bond is cut. 
*/
void initspins(Natural nv, int numlines, int seed, int *magnit, int *energy){
  int l;
  Natural i,n;

  /* reset cutlist to the bond values */
  *magnit = 0;
  for(l = 0; l<numlines; l++) cutlist[l] = bondlist[l];
  rng1(seed);
  for(n = 0; n<nv; n++){
    V[n].updatetime = 0;
    V[n].spin = 2*(rng1(0) < .5)-1;
/*    V[n].spin = ((n/2)*2 == n) ? -1 : 1;   */
    *magnit += V[n].spin;
    /* A negative spin flips the state of each of its edges */ 
    if(V[n].spin < 0){
       for(i = 0; i<V[n].links; i++){
	 *V[n].bondcut[i] *= -1;
       }
    }
  }
  *energy = 0;
  for(l = 0; l<numlines; l++) if(cutlist[l] < 0) ++(*energy);
}


/* initfitness: sets initial fitness value for each variable between
0 <= fitness < MAXLINK. Each variable has fitness = #(cut bonds), thus 
MAXLINK must be bigger than the number of bonds of the most connected variable.
NOTE: the best "fitness" is zero while the worst is each variables 
connectivity (>0 !).
*/ 
void initfitness(Natural nv){
  Natural n,i;
  for(n = 0; n<nv; n++){
    V[n].fitness = 0;
        /* bond violated, increase fitness (bigger = worse!) */
    for(i = 0; i<V[n].links; i++) if(*V[n].bondcut[i] < 0) V[n].fitness++; 
  }
}

/********** END fixing bonds          *****************************/

/******************************************************************/
/*                                                                */
/*      tau-probability distribution                              */
/*                                                                */
/*  This provides a discrete approximation for P(k)~k^{-tau},     */
/* 1 <= k <= nv. It maps the distribution onto a large array where*/
/* entries = k occur with frequency in occordance with P(k).      */
/* Shortcoming: if "LENGTH" too small, and/or "tau" too large,    */
/* there will be only values of k<n_0<nv for some cutoff n_0, ie. */
/* certain ranks can't be reached!                                */
/*                                                                */
/******************************************************************/


void inittau(double tau, double nv){
  int i;
  double a,b;
  a = (1-pow(nv+1,1-tau))/(double) LENGTH;
  b = 1/(1-tau);
  for(i = 0; i<LENGTH; i++){
    pdf[i] = (Natural) pow(1-i*a,b);
  }
}


/* taupick: rank is selected according to the tau-probability distribution.
 */
Natural taupick(){
  return(pdf[(int) (LENGTH*rng1(0))]);
}


/*********************** END tau probability **********************/



/******************************************************************/
/*                                                                */
/*      Boxed rank-order list                                     */
/*                                                                */
/* This rank-ordering is useful when each variable only takes on  */
/* a small set of values, resulting in large degeneracies. All    */
/* variables of same fitness are put in the same fitlist and drawn at */
/* random, if the selected rank lies within the interval of ranks */
/* covered by the fitlist. Thus, the ordering is exact and the access */
/* time is merely O(1)!                                           */
/*                                                                */
/******************************************************************/

/* Initiates the rank list in terms of a ordered set of boxes. Since there
 * is only a finite set of distinct fitnesses, fitness = 0,...,MAXLINK, there
 * is a large degeneracy between vertex-fitnesses. fitlist[i][j] contains all
 * pointers-to-vertex for the j-th vertex of fitness = i, where 0 <= j <= len[i]
 * is the size of the i-th fitlist. "minfitness" is the fitlist corresponding to
 * the worst fitness (remember: worst = highest!) that is non-empty, 
 * ie len[minfitness]>0 and len[i] = 0 f.a. i>minfitness.
 */
void initrank(Natural nv){
  Natural n,i;
  minfitness = 0;
  for(i = 0; i<MAXLINK; i++) len[i] = 0;
  for(n = 0; n<nv; n++){
    i = V[n].fitness;
    if(i > minfitness) minfitness = i;
    V[n].rank = len[i];
    fitlist[i][len[i]++] = &V[n];
  }
}
 
/* rerank: removes vertex v out of the "entry" position in its current fitlist 
 * "old", moves the last vertex in fitlist[old] from position len[old] to the
 * "entry" position (to keep fitlist[old] compact!), then makes v the new, 
 * last element in fitlist[new].
 * Also, "minfitness" may have to be adjusted, up or down.
 */ 
void rerank(struct Vertex *v, Natural new){
  Natural entry,old;
  old = v->fitness;
      /* no need to change if old=new */
  if(old == new) return;
  entry = v->rank;
      /* move top-of-list into emptied spot, decrement list */
  fitlist[old][entry] = fitlist[old][--len[old]]; 
      /* update location key for moved item */
  fitlist[old][entry]->rank = entry; 
  v->fitness = new;
      /* add new entry on top of new list */
  fitlist[new][len[new]] = v;
      /* update location key for entered item */
  v->rank = len[new]++;
      /* reset minfitness to worst non-zero fitlist */
  if(new > minfitness) 
     minfitness = new;
  else
     while(len[minfitness] == 0) minfitness--;
}


/* findvertex: returns the pointer to the k-th ranked vertex.
 * If the k-th ranked vertex is to be chosen, we look through the fitlist
 * starting with fitlist[minfitness] in order of increasing fitness until the
 * fitlist "i" which contains the k-th worst-fit vertex is located. One element 
 * of fitlist[i] is selected at random. This procedure approximates the 
 * distribution well enough.
 */
struct Vertex *findvertex(Natural rank){
  Natural i = minfitness;
  while(rank > len[i]) rank -= len[i--];
  return(fitlist[i][(Natural) (len[i]*rng1(0))]);
}

/*********************** END fitlist ranking ***********************/



/******************************************************************/
/*                                                                */
/*           Random Number Generator                              */
/*                                                                */
/* This one does the multiplicative LCG of Lewis, Goodman and     */
/* Miller, using the Schrage modulus trick                        */
/*                                                                */
/******************************************************************/

double rng1(long int n){
  static const long int ia = 16807;
  static const long int RMAX = 2147483647;
  static const double conv1 = .4656612875245797e-9; /* =1/RMAX */
  static const long int iq = 127773;
  static const long int ir = 2836;
  static long int irandom;
  long int l;

  if(n > 0) irandom = n;
  l = irandom/iq;
  irandom = (irandom-iq*l)*ia - ir*l;
  if(irandom < 0) irandom = irandom+RMAX;
  return(conv1*irandom);
}

/************  END Random Number Generator ***********************/

  /* set up a plotter window            */
  /* an animation with double buffering */
void initplotwindow(){
  pl_parampl ("BITMAPSIZE", "700x700");
  pl_parampl ("VANISH_ON_DELETE", "yes");
  pl_parampl ("USE_DOUBLE_BUFFERING", "fast"); 
  pl_handle = pl_newpl ("X", stdin, stdout, stderr);
  pl_selectpl (pl_handle);
  pl_openpl ();
  pl_fspace (0.0,-0.2*pl_scale,pl_scale,pl_scale);
  pl_flinewidth (0.001*pl_scale);

  /* run over-ridden keyboard ala' DOS */
  /*init_keyboard();*/
}


void plot(int L,float f,float tau,int t,int energy,int minenergy,float dx,int timescale){     
  char sbuff[10];
  int n,lx,ly,color1,color2,color0;
  float x,y,x1,y1;
   
  pl_erase();
  pl_fmove(0.0,-0.15*pl_scale);
  pl_label("L = ");
  pl_fmove(0.04*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(L,3,sbuff));
  pl_fmove(0.09*pl_scale,-0.15*pl_scale);
  pl_label("tau = ");
  pl_fmove(0.14*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(tau,3,sbuff));

  pl_fmove(0.21*pl_scale,-0.15*pl_scale);
  pl_label("p+ = ");
  pl_fmove(0.27*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(f,3,sbuff));

   pl_fmove(0.36*pl_scale,-0.15*pl_scale);
  pl_label("E = ");
  pl_fmove(0.4*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(energy,6,sbuff));
  pl_fmove(0.5*pl_scale,-0.15*pl_scale);
  pl_label("E_0 = ");
  pl_fmove(0.58*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(minenergy,6,sbuff));
  pl_fmove(0.75*pl_scale,-0.15*pl_scale);
  pl_label("t = ");
  pl_fmove(0.8*pl_scale,-0.15*pl_scale);
  pl_label(gcvt(t,8,sbuff));
/*
  pl_filltype(0);
  pl_fbox(0.0,0.0,pl_scale,pl_scale);
*/
  y1 = 0;
  n = 0;
  for(ly=0; ly<L; ly++){
    x1 = 0;
    y = y1;
    y1 += dx;
    for(lx=0; lx<L; lx++){
      x = x1;
      x1 += dx;
      color0 = 0;
      color1 = LENGTH;
      /* color2 makes sites fade from black through gray to white according to the last time they had been update (white=old) */
      color2 = abs(timescale*(t-V[n].updatetime));
      color2 = min(color2,LENGTH);
    switch (V[n].fitness){
      case 0: 
	pl_fillcolor(color2,color2,color2);
	break;
      case 1:
	pl_fillcolor(color0,color2,color1);
	break;
      case 2:
	pl_fillcolor(color2,color1,color0);
	break;
      default:
	pl_fillcolor(color1,color0,color2);
	break;
      }
/*    pl_color(100*n,100*n+2891,10*n+21841); */
      pl_filltype(1);
      pl_fbox(x,y,x1,y1);
      ++n;
    }
  }
}

void closeplotwindow(void){
      pl_closepl ();
      pl_selectpl (0);
      pl_deletepl (pl_handle);
      /* restore keyboard defaults */
      /*      close_keyboard(); */
      exit(0);   
}


void init_keyboard(void){
    tcgetattr(0, &initial_settings);
    new_settings = initial_settings;
    new_settings.c_lflag &= ~ICANON; 
    new_settings.c_lflag &= ~ECHO;
    new_settings.c_lflag &= ~ISIG; 
    new_settings.c_cc[VMIN] = 1;
    new_settings.c_cc[VTIME] = 0; 
    tcsetattr(0, TCSANOW,&new_settings);   
}

void close_keyboard(void){
  tcsetattr(0, TCSANOW,&initial_settings); 
}

int kbhit(void){
    char ch;
    int nread;
    if(peek_character != -1)
        return 1;
  new_settings.c_cc[VMIN] = 0;  
  tcsetattr(0, TCSANOW,&new_settings); 
  nread=read(0,&ch,1);
   new_settings.c_cc[VMIN] = 1;  
  tcsetattr(0, TCSANOW,&new_settings);  
  if(nread == 1){
      peek_character=ch;
      return 1;
  }

  return 0;
}

int readch(void){
    char ch;
    if(peek_character != -1){
        ch=peek_character;
        peek_character = -1;
            return ch;
    }
    read(0,&ch,1);
    return ch;
}