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

typedef struct {
  int i;
  int s;
  int k;
  int * a;
  double h;
} cell;

typedef struct {
  cell ** c;
  int N;
} network;

cell * new_cell(int N, int i, int K, double p);
network * new_network(int N, int K, double p);
void compute_field(network * x);
void evolve(network * x, double J, double W);
double magnetization(network *x);
void randomize(network * x);
void destroy_network(network *x);

int main() {
  network * x;
  int N;
  int t, T, TRANS;
  int K;
  double J, W;
  double p; // fraction of rewiring
  double m;
  
  
  J = -8;
  W = 10;
  N = 1000;
  T = 100;
  TRANS=1000;
  K = 20;
  p = 0;

  srand48(time(NULL));
  
  
//  for (J =0; J > -10; J-=0.2 ) {
  for (p=0;p<=1; p+=0.005) {
    x = new_network(N, K, p);
    randomize(x);
    for (t=0; t<T+TRANS; t++) {
      compute_field(x);
#if 0
      for (int i=0; i<N; i++) {
        printf("%d: %d [", i, x->c[i]->s);
        for (int j=0; j<x->c[i]->k; j++) {
          printf("%d ", x->c[i]->a[j]);
        }
        printf("] -> %f\n ",  x->c[i]->h);
      }
      getchar();
#endif
      evolve(x, J, W);
      //display(x, t);
      if (t>TRANS) {
        m = magnetization(x);
//        printf("%f %f\n", J, m);
        printf("%f %f\n", p, m);
      }
    }
    destroy_network(x);
  }
}

void evolve(network * x, double J, double W) {
  int i;
  
  for (i=0; i<x->N; i++) {
    if (drand48() < 1/(1+exp(-2*(J*x->c[i]->h + W * pow(x->c[i]->h, 3))))) {
      x->c[i]->s = 1;
    } else {
      x->c[i]->s = -1;
    }
  }
}



void compute_field(network * x) {
  int i, j, k;
  
  for (i=0; i<x->N; i++) {
    x->c[i]->h = 0;
    for (j=0; j<x->c[i]->k; j++) {
      k = x->c[i]->a[j];
      x->c[i]->h += x->c[k]->s;
    }
    x->c[i]->h /= x->c[i]->k;
  }
}

double magnetization(network *x){
  int i;
  double m;
  
  m =0;
  for (i=0; i<x->N; i++) {
    m += x->c[i]->s;
  }
  return m/x->N;
}

network * new_network(int N, int K, double p) {
  network * x;
  int i;
  
  x = calloc(1, sizeof(network));
  x->N=N;
  x->c = calloc(N, sizeof(cell *));
  for (i=0; i<N; i++) {
    x->c[i] = new_cell(N, i, K, p);
  }
  return(x);
}  

void destroy_network(network *x) {
  int i;
  for (i=0; i<x->N; i++) {
    free (x->c[i]->a);
    free (x->c[i]);
  }
  free (x->c);
  free(x);
}

void randomize(network *x) {
  int i;
  
  for (i=0; i<x->N; i++) {
    x->c[i]->s = drand48() < 0.5? -1 : 1;
  }
}


cell * new_cell(int N, int i, int K, double p) {
  cell * c;
  int j;
  
  c = calloc(1, sizeof(cell));
  c->i = i;
  c->k =  K;
  c->s = drand48() < 0.5? -1 : 1;
  c->a = calloc(K, sizeof(int));
  for (j=0; j<K; j++) {
    if (drand48() < p) {
      c->a[j] = (int) (drand48() * N);
    } else {
      c->a[j] = (i + j - (K/2) + N) % N;
    }
  }
  return(c);
}