First commit of MAMMULT code

6 files changed, 1132 insertions, 0 deletions

diff --git a/dynamics/ising/Makefile b/dynamics/ising/Makefile
new file mode 100644
index 0000000..1c09fa1
--- /dev/null
+++ b/dynamics/ising/Makefile
@@ -0,0 +1,11 @@
+CFLAGS="-O3"
+CC="gcc"
+MFLAG=-lm
+
+all: multiplex_ising 
+
+multiplex_ising: multiplex_ising.c 
+	$(CC) $(CFLAGS) -o multiplex_ising multiplex_ising.c utils.c iltree.c $(MFLAG)
+
+clean:
+	rm multiplex_ising 
diff --git a/dynamics/ising/iltree.c b/dynamics/ising/iltree.c
new file mode 100755
index 0000000..5a693d0
--- /dev/null
+++ b/dynamics/ising/iltree.c
@@ -0,0 +1,201 @@
+/*
+ *
+ * A simple insert-lookup static btree datatype
+ *
+ */
+
+#include <stdlib.h>
+#include "iltree.h"
+#include <stdio.h>
+
+
+void __recursive_preorder(node_t *cur, ilfunc_t *funs){
+  
+  if(cur->left){
+    __recursive_preorder(cur->left, funs);
+  }
+  funs->print(cur->info, funs->fileout);
+  if(cur->right){
+    __recursive_preorder(cur->right, funs);
+  }
+}
+
+/*
+ *
+ * Recursive push of nodes in the nodecache :-)
+ *
+ */
+
+void __recursive_destroy(node_t *cur, ilfunc_t *funs){
+  if(cur->left){
+    __recursive_destroy(cur->left, funs);
+    cur->left = NULL;
+  }
+  if(cur->right){
+    __recursive_destroy(cur->right, funs);
+    cur->right = NULL;
+  }
+}
+
+
+int __recursive_insert(node_t *cur, node_t *elem, ilfunc_t *f){
+  
+  int res ;
+  res = f->compare(cur->info, elem->info);
+  /*  printf("res: %d\n", res); */
+  if ( res > 0){
+    if (cur->left){
+      return __recursive_insert(cur->left, elem, f);
+    }
+    else{
+      cur->left = elem;
+      return 0;
+    }
+  }
+  else if (res < 0){
+    if (cur->right){
+      return __recursive_insert(cur->right, elem, f);
+    }
+    else{
+      cur->right = elem;
+      return 0;
+    }
+  }
+  printf("warning!!!!! duplicate entry!!!!!!\n\n");
+  return -1;
+}
+
+
+
+void* __recursive_lookup(node_t *cur, void *v, ilfunc_t *f){
+  
+  int res;
+
+  res = f->compare(cur->info, v);
+
+  if (res > 0){
+    if(cur->left)
+      return __recursive_lookup(cur->left, v, f);
+    else
+      return NULL;
+    
+  }
+  else if (res < 0){
+    if(cur->right)
+      return __recursive_lookup(cur->right, v, f);
+    else
+      return NULL;
+  }
+  else
+    return cur->info;
+}
+
+void __recursive_map(node_t *cur, void (*func)(void*)){
+  
+  if (cur->left)
+    __recursive_map(cur->left, func);
+  func(cur->info);
+  if (cur->right)
+    __recursive_map(cur->right, func);
+}
+
+void __recursive_map_args(node_t *cur, void (*func)(void*, void*), void *args){
+  
+  if (cur->left)
+    __recursive_map_args(cur->left, func, args);
+  func(cur->info, args);
+  if (cur->right)
+    __recursive_map_args(cur->right, func, args);
+}
+
+
+
+iltree_t iltree_create(iltree_t t){
+  if (!t) {
+    t = (iltree_t)malloc(sizeof(iltree_struct_t));
+  }
+  t->root = NULL;
+  return t;
+}
+
+
+void iltree_set_funs(iltree_t t, ilfunc_t *funs){
+  
+  t->funs = *funs;
+}
+
+
+void iltree_insert(iltree_t t, void *elem){
+  
+  node_t *n;
+  
+  n = (node_t*)malloc(sizeof(node_t));
+  n->info = t->funs.alloc(); 
+  t->funs.copy(elem, n->info);
+  n->left = n->right = NULL;
+  if (t->root == NULL){
+    t->root = n;
+  }
+  else{
+    __recursive_insert(t->root, n, & (t->funs));
+  }
+}
+
+
+void iltree_destroy(iltree_t t){
+  
+  if(t->root)
+    __recursive_destroy(t->root, & (t->funs));
+  free(t);
+}
+
+
+
+
+void iltree_view_pre(iltree_t t){
+  
+  if (t->root){
+    /*printf("----\n");*/
+    __recursive_preorder(t->root, & (t->funs));
+    /*printf("----\n");*/
+  }
+  else
+    printf("----- Empty tree!!!! -----\n");
+  
+}
+
+
+
+void* iltree_lookup(iltree_t t , void *elem){
+
+  node_t n;
+  
+  if(t->root)
+    return __recursive_lookup(t->root, elem, & (t->funs) );
+  else 
+    return NULL;
+}
+
+
+void iltree_map(iltree_t t, void (*func)(void*)){
+  
+  __recursive_map(t->root, func);
+  
+}
+
+
+void iltree_map_args(iltree_t t, void (*func)(void*, void*), void *args){
+  
+  __recursive_map_args(t->root, func, args);
+  
+}
+
+void* iltree_get_fileout(iltree_t t){
+
+  return t->funs.fileout;
+}
+
+void iltree_set_fileout(iltree_t t, void *f){
+  
+  t->funs.fileout = f;
+}
diff --git a/dynamics/ising/iltree.h b/dynamics/ising/iltree.h
new file mode 100755
index 0000000..3e835e4
--- /dev/null
+++ b/dynamics/ising/iltree.h
@@ -0,0 +1,55 @@
+#ifndef __ILTREE_H__
+#define __ILTREE_H__
+
+
+typedef struct node{
+  void* info;
+  struct node* left;
+  struct node* right;
+} node_t;
+
+typedef struct{
+  void* (*alloc)();
+  void (*dealloc)(void*);
+  void (*copy)(void *src, void *dst);
+  int (*compare)(void*, void*);
+  void (*print)(void*, void*);
+  void *fileout;
+} ilfunc_t;
+
+
+typedef struct {
+  node_t* root;
+  ilfunc_t funs;
+} iltree_struct_t;
+
+
+
+typedef iltree_struct_t* iltree_t;
+
+
+void iltree_set_funs(iltree_t, ilfunc_t *);
+
+void iltree_destroy(iltree_t);
+
+void iltree_empty(iltree_t);
+
+void iltree_insert(iltree_t, void*);
+
+void* iltree_lookup(iltree_t, void*);
+
+void iltree_view_pre(iltree_t);
+
+iltree_t iltree_create(iltree_t);
+
+void iltree_empty_cache(iltree_t);
+
+void iltree_map(iltree_t, void (*func)(void*));
+
+void iltree_map_args(iltree_t, void (*func)(void*, void*), void*);
+
+void* iltree_get_fileout(iltree_t t);
+
+void iltree_set_fileout(iltree_t t, void *f);
+
+#endif /* __ILTREE_H__*/
diff --git a/dynamics/ising/multiplex_ising.c b/dynamics/ising/multiplex_ising.c
new file mode 100644
index 0000000..8dbbf68
--- /dev/null
+++ b/dynamics/ising/multiplex_ising.c
@@ -0,0 +1,313 @@
+#include <stdio.h>

+#include <stdlib.h>

+#include <time.h>

+#include <math.h>

+

+#include "utils.h"

+#include "iltree.h"

+

+

+

+typedef struct{

+  double T;

+  double J;

+  double Jcoup;

+  double p0;

+  double p1;

+  double h0;

+  double h1;

+  unsigned int num_epochs;

+} param_t;

+

+

+typedef struct{

+  unsigned int N;

+  unsigned int K;

+  unsigned int *J_slap;

+  unsigned int *r_slap;

+  int *s;

+} net_t;

+

+typedef struct{

+  double m0;

+  double m1;

+  double C;

+  double b0;

+  double b1;

+  double q0;

+  double q1;

+  double IFC0;

+  double IFC1;
+  double M;

+} stats_t;

+

+void init_spins(int *s, int N, double p){

+  

+  int i;

+  double val;

+

+  for (i=0; i<N; i ++){

+    val = 1.0 * rand() / RAND_MAX;

+    if (val < p){

+      s[i] = 1;

+    }

+    else

+      s[i] = -1;

+  }

+

+}

+
+
+void init_spins_once(param_t *sys, net_t *layers) {
+  int N, N2;
+  N = layers[0].N;

+  N2 = 2 * N;
+
+  init_spins(layers[0].s, N, sys->p0);

+  init_spins(layers[1].s, N, sys->p1);
+	}
+

+void shuffle_ids(int * v, int N){

+  

+  int tmp, val, i;

+

+  for (i=N-1; i >=0; i --){

+    val = rand() % (i+1);

+    tmp = v[i];

+    v[i] = v[val];

+    v[val] = tmp;

+  }

+}

+

+void compute_stats (net_t *layers, stats_t *stats, int *spinold0, int *spinold1){

+  

+  int i, N;

+  

+  N = layers[0].N;

+  

+  stats->M =stats->m0 = stats->m1 = stats->C = stats->b0 = stats->b1 = stats->q0 = stats->q1 = stats->IFC0 = stats->IFC1 = 0;

+  

+  double *bubblevec0, *bubblevec1;

+  bubblevec0 = (double *)malloc(N * sizeof(double));

+  bubblevec1 = (double *)malloc(N * sizeof(double));

+	

+  

+  double bubble0, bubble1, deg;

+  int j;

+  for (i=0; i<N; i++){

+    stats->m0 += layers[0].s[i];

+    stats->m1 += layers[1].s[i];

+    stats->C += layers[0].s[i] * layers[1].s[i];

+    

+

+    bubble0=0;

+    for(j=layers[0].r_slap[i]; 

+          j< layers[0].r_slap[i + 1]; 

+          j ++){

+          bubble0+= layers[0].s[i] * layers[0].s[ layers[0].J_slap[j] ] ;

+	  stats->IFC0 += fabs(layers[0].s[i] * layers[0].s[ layers[0].J_slap[j] ]-1.)/2.;

+          }

+    deg = (layers[0].r_slap[i + 1] - layers[0].r_slap[i])*1.0;

+    bubblevec0[i]=bubble0/deg;

+    bubble1=0;

+    for(j=layers[1].r_slap[i]; 

+          j< layers[1].r_slap[i + 1]; 

+          j ++){

+          bubble1+= layers[1].s[i] * layers[1].s[ layers[1].J_slap[j] ] ;

+	  stats->IFC1 += fabs(layers[1].s[i] * layers[1].s[ layers[1].J_slap[j] ]-1.)/2.;

+          }

+    deg = (layers[1].r_slap[i + 1] - layers[1].r_slap[i])*1.0;

+    bubblevec1[i]=bubble1/deg;

+

+    stats->q0 += layers[0].s[i]*spinold0[i];

+    stats->q1 += layers[1].s[i]*spinold1[i];

+  }

+

+  stats->b0=0;

+  for (i=0; i<N; i++) {

+	stats->b0 = stats->b0 + bubblevec0[i];

+	}

+  stats->b0 /= N;

+  stats->b1=0;

+  for (i=0; i<N; i++) {

+	stats->b1 = stats->b1 + bubblevec1[i];

+	}

+  stats->b1 /= N;

+

+  stats->m0 /= N;

+  stats->m1 /= N;

+  stats->C /= N;

+  stats->q0 /= N;

+  stats->q1 /= N;

+  stats->IFC0 /= layers[0].K;

+  stats->IFC1 /= layers[1].K;
+  
+  stats->M = (fabs(stats->m0)+fabs(stats->m1))/2.0;

+ }

+

+void dump_spins(int *s1, int *s2, int N){

+  

+  int i;

+

+  for(i=0; i<N; i++){

+    printf("%d %d\n", s1[i], s2[i]);

+  }

+  

+}

+

+

+void make_simulation(param_t *sys, net_t *layers, stats_t *stats){

+  

+ 

+  int  *ids;

+  int e, i, j, num_flips, id, l;

+  unsigned int N, N2;

+

+  double E_old, E_new, val, exp_val;

+

+  N = layers[0].N;

+  N2 = 2 * N;

+    

+  ids = malloc(N2 * sizeof(int));

+  for (i=0; i< N2; i++){

+    ids[i] = i;

+  }

+

+  int *spinold0, *spinold1;

+  spinold0 = (int *)malloc(N * sizeof(int));

+  spinold1 = (int *)malloc(N * sizeof(int));

+

+

+  

+  for (e=0; e < sys->num_epochs; e++){

+    num_flips = 0;

+    shuffle_ids(ids, N2);

+    for (i=0; i< N2; i++){

+      id = ids[i] % N;

+      l = ids[i]/N;

+      //printf("i: %d id: %d l:%d\n", ids[i], id, l);

+      E_old = 0;

+      E_new = 0;

+      for(j=layers[l].r_slap[id]; 

+          j< layers[l].r_slap[id + 1]; 

+          j ++){

+        E_old -= layers[l].s[id] * layers[l].s[ layers[l].J_slap[j] ];

+        E_new -= -layers[l].s[id] * layers[l].s[ layers[l].J_slap[j] ];

+      }

+      E_old *= sys->J;

+      E_new *= sys->J;

+

+      if (l==0) {

+	E_old -= sys->h0 * layers[l].s[id];

+	E_new -= sys->h0 * (-layers[l].s[id]);

+	}

+      else {

+	E_old -= sys->h1 * layers[l].s[id];

+	E_new -= sys->h1 * (-layers[l].s[id]);

+	}

+

+      E_old -= sys->Jcoup * layers[l].s[id] * layers[1-l].s[id];

+      E_new -= sys->Jcoup * -(layers[l].s[id]) * layers[1-l].s[id];

+

+      E_old = 2*E_old;

+      E_new = 2*E_new;

+      if (E_new <= E_old){ /* The new configuration has smaller energy -> flip! */

+        layers[l].s[id] = - layers[l].s[id];

+      }

+      else if (sys->T > 0){ /* The new conf has higher energy -> flop with e-(\Delta E)/T */

+        val = 1.0 * rand() / RAND_MAX;

+        exp_val = exp( - (1.0*(E_new - E_old)) / sys->T);

+        if (val < exp_val){

+          layers[l].s[id] = - layers[l].s[id];

+        }

+      }

+    }

+

+  if (e==(sys->num_epochs-2)) {

+	int u;

+	for (u=0; u<N; u++) {

+		spinold0[u]=layers[0].s[u];

+		spinold1[u]=layers[1].s[u];

+		}

+	}

+

+

+

+  }

+  compute_stats(layers, stats, spinold0, spinold1);
+  

+  //free(spinold0);

+  //free(spinold1);

+  //dump_spins(layers[0].s, layers[1].s, N);

+}

+

+void dump_stats(param_t *sys, stats_t *s){

+  

+  printf("## T J gamma h0 h1 p1 p2 m1 m2 C\n");

+  printf("%g %g %g %g %g %g %g %g %g %g\n", 

+         sys->T, sys->J, sys->Jcoup, sys->h0, sys->h1, sys->p0, 

+         sys->p1, s->m0, s->m1, s->C);
+  fflush(stdout);

+}

+

+

+int main(int argc, char *argv[]){

+  

+

+  net_t layers[2];

+  param_t sys;

+  stats_t stats;

+  unsigned int *J_slap, *r_slap;

+  

+  FILE *fin;

+

+  if (argc < 10){

+    printf("Usage: %s <layer1> <layer2> <T> <J> <gamma> <h1> <h2> <p1> <p2> <num_epochs>\n", argv[0]);

+    exit(1);

+  }

+

+  sys.T = atof(argv[3]);

+  sys.J = atof(argv[4]);

+  sys.Jcoup = atof(argv[5]);

+  sys.p0 = atof(argv[8]);

+  sys.p1 = atof(argv[9]);

+  sys.h0 = atof(argv[6]);

+  sys.h1 = atof(argv[7]);

+  sys.num_epochs = atoi(argv[10]);

+  

+  srand(time(NULL));

+

+  J_slap = r_slap = NULL;

+  

+  fin = openfile_or_exit(argv[1], "r", 2);

+  read_slap(fin, &(layers[0].K), &(layers[0].N), &J_slap, &r_slap);

+  layers[0].J_slap = J_slap;

+  layers[0].r_slap = r_slap;

+  fclose(fin);

+

+  J_slap = r_slap = NULL;

+

+  fin = openfile_or_exit(argv[2], "r", 2);

+  read_slap(fin, &(layers[1].K), &(layers[1].N), &J_slap, &r_slap);

+  layers[1].J_slap = J_slap;

+  layers[1].r_slap = r_slap;

+  fclose(fin);

+

+  if (layers[0].N != layers[1].N){

+    printf("Error!!! Both layers must have the same number of nodes!!!!\n");

+    exit(3);

+  }

+

+   /* allocate space for the spins on each layer */

+  layers[0].s = malloc(layers[0].N * sizeof(double));

+  layers[1].s = malloc(layers[1].N * sizeof(double));

+  
+  /* inizialize the spins */
+  init_spins_once(&sys, layers);
+
+

+  make_simulation(&sys, layers, &stats);

+  

+ 

+  

+}

diff --git a/dynamics/ising/utils.c b/dynamics/ising/utils.c
new file mode 100755
index 0000000..952fcd7
--- /dev/null
+++ b/dynamics/ising/utils.c
@@ -0,0 +1,494 @@
+#include "iltree.h"
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "utils.h"
+
+void* alloc_double(){
+  return malloc(sizeof(long double));
+}
+
+void dealloc_double(void *elem){
+  free(elem);
+}
+
+void copy_double(void *elem1, void *elem2){
+  *((long double*)elem2) = *((long double*)elem1);
+}
+
+int compare_double(void *elem1, void *elem2){
+  return *((long double*)elem1) - *((long double*)elem2);
+}
+
+void print_double(void *elem, void *fileout){
+  
+  long double k, i, j;
+  long double x;
+
+  k = *((long double*)elem);
+
+  x = (1 + sqrtl(1 + 8 * (k-1))) / 2;
+  i = floorl(x) + 1;
+  j = k - ( (i-1)*1.0 * (i-2) ) /2;
+  //printf("x: %Lf\n i: %0.0Lf j: %0.0Lf\n", x, i, j);
+  fprintf((FILE*)fileout, "%0.0Lf %0.0Lf\n", i-1, j-1);
+}
+
+iltree_t init_tree(iltree_t t, void *fileout){
+  
+  ilfunc_t funs= {
+    .alloc = alloc_double,
+    .dealloc = dealloc_double, 
+    .copy = copy_double,
+    .compare = compare_double, 
+    .print = print_double,
+    .fileout = fileout
+  };
+  
+  t = iltree_create(t);
+  iltree_set_funs(t, &funs);
+  return t;
+}
+
+
+/* @@@@ CAMBIARE IL PRIMO PARAMETRO IN UN FILE* PER RENDERLA COERENTE
+   ALLE ALTRE FUNZIONI DI READ!!!!*/
+int read_deg_distr(FILE *filein, unsigned int **degs, unsigned int **Nk, double **p){
+  
+  int n_degs = 0;
+  int size = 10;
+  char *line;
+  char buff[256];
+  int k_i, num_i;
+  double p_i;
+  char *ptr;
+
+  line = NULL;
+
+  *degs = realloc(*degs, size*sizeof(unsigned int));
+  *Nk = realloc(*Nk, size*sizeof(unsigned int));
+  *p = realloc(*p, size*sizeof(double));
+  
+
+  while(fgets(buff, 256, filein)){
+    ptr = strtok(buff, " ");
+    if (ptr[0] == '#')
+      continue;
+    k_i = atoi(ptr);
+    ptr = strtok(NULL, " " );
+    num_i = atoi(ptr);
+    ptr = strtok(NULL, " \n");
+    p_i = atof(ptr);
+    if (n_degs == size){
+      size += 10;
+      *degs = realloc(*degs, size*sizeof(unsigned int));
+      *Nk = realloc(*Nk, size*sizeof(unsigned int));
+      *p = realloc(*p, size*sizeof(double));
+    }
+    (*degs)[n_degs] = k_i;
+    (*Nk)[n_degs] = num_i;
+    (*p)[n_degs] = p_i;
+    n_degs += 1;
+  }
+  *degs = realloc(*degs, n_degs*sizeof(unsigned int));
+  *Nk = realloc(*Nk, n_degs*sizeof(unsigned int));
+  *p = realloc(*p, n_degs*sizeof(double));
+  return n_degs;
+}
+
+
+int read_deg_seq(FILE *filein, unsigned int **nodes){
+
+  int size, N, k;
+  char buff[256];
+  char *ptr;
+
+  N = 0;
+  size = 10;
+  
+  *nodes = (unsigned int*)malloc(size * sizeof(unsigned int));
+
+  while(fgets(buff, 256, filein)){
+    ptr = strtok(buff, " ");
+    if (ptr[0] == '#')
+      continue;
+    k = atoi(ptr);
+    
+    if (N == size){
+      size += 10;
+      *nodes = realloc(*nodes, size*sizeof(unsigned int));
+    }
+    (*nodes)[N] = k;
+    N += 1;
+  }
+  *nodes = realloc(*nodes, N * sizeof(unsigned int));
+  return N;
+}
+
+int read_stubs(FILE *filein, unsigned int **S){
+  
+  int size, K;
+  char buff[256];
+  char *ptr;
+  
+  K=0;
+  size = 20;
+  *S = malloc(size * sizeof(unsigned int));
+  
+  while(fgets(buff, 256, filein)){
+    if (K == size){
+      size += 20;
+      *S = realloc(*S, size*sizeof(unsigned int));
+    }
+    ptr = strtok(buff, " "); /* read the first node */
+    (*S)[K++] = atoi(ptr);
+    ptr = strtok(NULL, " "); /* read the second node */
+    (*S)[K++] = atoi(ptr);
+  }
+  *S = realloc(*S, K * sizeof(unsigned int));
+  return K;
+}
+
+/*
+ * Read a file in ij format 
+ */
+int read_ij(FILE *filein, unsigned int **I, unsigned int **J){
+  
+  unsigned int size, K;
+  char buff[256];
+  char *ptr;
+  
+  size = 20;
+  K = 0;
+  
+  *I = malloc(size * sizeof(unsigned int));
+  *J = malloc(size * sizeof(unsigned int));
+  while(fgets(buff, 256, filein)){
+    if (buff[0] == '#')
+      continue;
+    if (K == size){
+      size += 20;
+      *I = realloc(*I, size*sizeof(unsigned int));
+      *J = realloc(*J, size*sizeof(unsigned int));
+    }
+    ptr = strtok(buff, " "); /* read the first node */
+    (*I)[K] = atoi(ptr);
+    ptr = strtok(NULL, " "); /* read the second node */
+    (*J)[K] = atoi(ptr);
+    K += 1;
+  }
+  
+  *I = realloc(*I, K * sizeof(unsigned int));
+  *J = realloc(*J, K * sizeof(unsigned int));
+  return K;
+}
+
+/*funzione pesata di moreno*/
+
+int read_ij_w(FILE *filein, unsigned int **I, unsigned int **J , double **W){
+  
+  unsigned int size, K;
+  char buff[256];
+  char *ptr;
+  
+  size = 20;
+  K = 0;
+  
+  *I = malloc(size * sizeof(unsigned int));
+  *J = malloc(size * sizeof(unsigned int));
+  *W = malloc(size * sizeof(double));
+  
+  while(fgets(buff, 256, filein)){
+    if (buff[0] == '#')
+      continue;
+    if (K == size){
+      size += 20;
+      *I = realloc(*I, size*sizeof(unsigned int));
+      *J = realloc(*J, size*sizeof(unsigned int));
+      *W = realloc(*W, size*sizeof(double));
+      if (*W==NULL) {
+	printf ("Errore");
+	exit(-1);
+	}
+    }
+    ptr = strtok(buff, " "); /* read the first node */
+    (*I)[K] = atoi(ptr);
+    ptr = strtok(NULL, " "); /* read the second node */
+    (*J)[K] = atoi(ptr);
+    ptr = strtok(NULL, " "); /* read the weight  */
+    (*W)[K] = atof(ptr);    
+    K += 1;
+  }
+  
+  *I = realloc(*I, K * sizeof(unsigned int));
+  *J = realloc(*J, K * sizeof(unsigned int));
+  *W = realloc(*W, K * sizeof(double));
+  return K;
+}
+
+
+
+void read_slap(FILE *filein, unsigned int *K, unsigned int *N, 
+              unsigned int **J_slap, unsigned int **r_slap){
+  
+  unsigned int *I=NULL, *J=NULL;
+  unsigned int i, k;
+  
+  k = read_ij(filein, &I, &J);
+  *K = 2 * k;
+  I = realloc(I, 2*k * sizeof(unsigned int));
+  J = realloc(J, 2*k * sizeof(unsigned int));
+  for (i=k; i<2*k; i ++){
+    I[i] = J[i-k];
+    J[i] = I[i-k];
+  }
+  
+  *N = convert_ij2slap(I, J, 2*k, r_slap, J_slap);
+  free(I);
+  free(J);
+  return;
+}
+
+/*funzione pesata di moreno*/
+
+void read_slap_w(FILE *filein, unsigned int *K, unsigned int *N, 
+              unsigned int **J_slap, unsigned int **r_slap, double **w_slap){
+  
+  unsigned int *I=NULL, *J=NULL;
+  double *W=NULL;
+  unsigned int i, k;
+  
+  k = read_ij_w(filein, &I, &J, &W);
+  *K = 2 * k;
+  I = realloc(I, 2*k * sizeof(unsigned int));
+  J = realloc(J, 2*k * sizeof(unsigned int));
+  W = realloc(W, 2*k * sizeof(double));
+  
+  for (i=k; i<2*k; i ++){
+    I[i] = J[i-k];
+    J[i] = I[i-k];
+    W[i] = W[i-k];
+  }
+  
+  *N = convert_ij2slap_w(I, J, W, 2*k, r_slap, J_slap, w_slap);
+  free(I);
+  free(J);
+  free(W);
+  return;
+}
+
+unsigned int find_max(unsigned int *v, unsigned int N){
+
+  unsigned int i, max;
+  
+  max = v[0];
+  i = 0;
+  while(++i < N){
+    if (v[i] > max)
+      max = v[i];
+  }
+  return max;
+}
+
+
+int convert_ij2slap(unsigned int *I, unsigned int *J, unsigned int K, 
+                    unsigned int ** r_slap, unsigned int **J_slap){
+  
+  unsigned int tmp, max;
+  unsigned int N;
+  unsigned int i, pos;
+  unsigned int *p;
+
+  max = find_max(I, K) + 1;
+  tmp = find_max(J, K) + 1;
+  if (tmp > max){
+    max = tmp ;
+  }
+  
+  *r_slap = malloc( (max+1) * sizeof(unsigned int));
+  p = malloc(max * sizeof(unsigned int));
+  
+  *J_slap = malloc(K * sizeof(unsigned int));
+  
+  memset(*r_slap, 0, (max+1) * sizeof(unsigned int));
+  for(i=0; i<max + 1; i++)
+    (*r_slap)[i] = 0;
+  memset(p, 0, max * sizeof(unsigned int));
+  (*r_slap)[0] = 0;
+  //fprintf(stderr, "WARNING!!!! R_SLAP[0] NOW IS SET TO ZERO!!!!!\n");
+  for(i=0; i<K; i++){
+    (*r_slap)[ I[i] + 1] += 1;
+  }
+  for(i=1; i<=max; i++){
+    (*r_slap)[i] += (*r_slap)[i-1];
+  }
+  for(i=0; i<K; i++){
+    pos = (*r_slap) [ I[i] ] + p[ I[i] ];
+    (*J_slap)[pos] = J[i];
+    p[ I[i] ] += 1;
+  }
+  free(p);
+  return max;
+}
+
+/*funzione pesata di moreno*/
+int convert_ij2slap_w(unsigned int *I, unsigned int *J, double *W, unsigned int K, 
+                    unsigned int ** r_slap, unsigned int **J_slap,double **w_slap){
+  
+  unsigned int tmp, max;
+  unsigned int N;
+  unsigned int i, pos;
+  unsigned int *p;
+
+  max = find_max(I, K) + 1;
+  tmp = find_max(J, K) + 1;
+  if (tmp > max){
+    max = tmp ;
+  }
+  
+  *r_slap = malloc( (max+1) * sizeof(unsigned int));
+  p = malloc(max * sizeof(unsigned int));
+  
+  *J_slap = malloc(K * sizeof(unsigned int));
+  *w_slap = malloc(K * sizeof(double));
+  
+  memset(*r_slap, 0, (max+1) * sizeof(unsigned int));
+  for(i=0; i<max + 1; i++)
+    (*r_slap)[i] = 0;
+  memset(p, 0, max * sizeof(unsigned int));
+  (*r_slap)[0] = 0;
+  //fprintf(stderr, "WARNING!!!! R_SLAP[0] NOW IS SET TO ZERO!!!!!\n");
+  for(i=0; i<K; i++){
+    (*r_slap)[ I[i] + 1] += 1;
+  }
+  for(i=1; i<=max; i++){
+    (*r_slap)[i] += (*r_slap)[i-1];
+  }
+  for(i=0; i<K; i++){
+    pos = (*r_slap) [ I[i] ] + p[ I[i] ];
+    (*J_slap)[pos] = J[i];
+    
+    (*w_slap)[pos] = W[i];
+    
+    p[ I[i] ] += 1;
+  }
+  free(p);
+  return max;
+}
+
+
+
+int convert_ij2slap_N(unsigned int *I, unsigned int *J, unsigned int K,   
+                      unsigned int N, unsigned int ** r_slap, 
+                      unsigned int **J_slap){
+  
+  unsigned int tmp, max;
+  unsigned int i, pos;
+  unsigned int *p;
+
+  max = N;
+  
+  *r_slap = malloc( (max+1) * sizeof(unsigned int));
+  p = malloc(max * sizeof(unsigned int));
+  
+  *J_slap = malloc(K * sizeof(unsigned int));
+  
+  memset(*r_slap, 0, (max+1) * sizeof(unsigned int));
+  for(i=0; i<max + 1; i++)
+    (*r_slap)[i] = 0;
+  memset(p, 0, max * sizeof(unsigned int));
+  (*r_slap)[0] = 0;
+  //fprintf(stderr, "WARNING!!!! R_SLAP[0] NOW IS SET TO ZERO!!!!!\n");
+  for(i=0; i<K; i++){
+    (*r_slap)[ I[i] + 1] += 1;
+  }
+  for(i=1; i<=max; i++){
+    (*r_slap)[i] += (*r_slap)[i-1];
+  }
+  for(i=0; i<K; i++){
+    pos = (*r_slap) [ I[i] ] + p[ I[i] ];
+    (*J_slap)[pos] = J[i];
+    p[ I[i] ] += 1;
+  }
+  free(p);
+  return max;
+}
+
+
+
+/* RIVEDERE QUESTA  FUNZIONE...... PASSARE UN FILE COME ARGOMENTO E
+   USARE fprintf */
+void dump_deg_distr(unsigned int *degs, double *p, int n){
+
+  int i;
+
+  for(i=0; i<n; i++){
+    printf("%d %2.6f\n", degs[i], p[i]);
+  }
+}
+
+
+
+/* RIVEDERE QUESTA  FUNZIONE...... PASSARE UN FILE COME ARGOMENTO E
+   USARE fprintf */
+void dump_deg_seq(unsigned int *nodes, int N){
+
+  int i;
+  for(i=0; i<N; i++){
+    printf("%d: %d\n", i, nodes[i]);
+  }
+}
+
+void dump_edges(iltree_t t){
+
+  iltree_view_pre(t);
+}
+
+FILE* openfile_or_exit(char *filename, char *mode, int exitcode){
+  
+  FILE *fileout;
+  char error_str[256];
+  
+  fileout = fopen(filename, mode);
+  if (!fileout){
+    sprintf(error_str, "Error opening file %s", filename);
+    perror(error_str);
+    exit(exitcode);
+  }
+  return fileout;
+}
+
+inline int compare_int(const void *x1, const void *x2){
+  return *((unsigned int*)x1) - *((unsigned int*)x2);
+}
+
+void write_edges(FILE *fileout, unsigned int *J_slap, 
+                 unsigned int *r_slap, unsigned int N){
+  
+  unsigned int i, j;
+  
+  for(i=0; i<N; i++){
+    for (j=r_slap[i]; j<r_slap[i+1]; j++){
+      if (J_slap[j] > i){
+        fprintf(fileout, "%d %d\n", i, J_slap[j]);
+      }
+    }
+  }
+}
+  
+
+/* Check if j is a neighbour of i */
+int is_neigh(unsigned int *J_slap, unsigned int *r_slap, unsigned int N, 
+             unsigned int i, unsigned int j){
+  
+  unsigned int l;
+  unsigned int count;
+  count = 0;
+  for(l=r_slap[i]; l<r_slap[i+1]; l++){
+    if (J_slap[l] == j)
+      count ++;
+  }
+  return count;
+}
diff --git a/dynamics/ising/utils.h b/dynamics/ising/utils.h
new file mode 100755
index 0000000..c844248
--- /dev/null
+++ b/dynamics/ising/utils.h
@@ -0,0 +1,58 @@
+#ifndef __UTILS_H__
+#define __UTILS_H__
+
+#include "iltree.h"
+
+iltree_t init_tree(iltree_t t, void *fileout);
+
+int read_deg_distr(FILE *filein, unsigned int **degs, unsigned int **Nk, double **p);
+
+int read_deg_seq(FILE *filein, unsigned int **nodes);
+
+int read_stubs(FILE *filein, unsigned int **S);
+
+int read_ij(FILE *filein, unsigned int **i, unsigned int **j);
+
+/*funzione pesata di moreno*/
+int read_ij_w(FILE *filein, unsigned int **i, unsigned int **j, double **w);
+
+void read_slap(FILE *filein, unsigned int *K, unsigned int *N, 
+               unsigned int **J_slap, unsigned int **r_slap);
+
+/*funzione pesata di moreno*/
+void read_slap_w(FILE *filein, unsigned int *K, unsigned int *N, 
+               unsigned int **J_slap, unsigned int **r_slap, double **w_slap);
+
+int convert_ij2slap(unsigned int *I, unsigned int *J, unsigned int K, 
+                    unsigned int ** r_slap, unsigned int **J_slap);
+
+/*funzione pesata di moreno*/
+int convert_ij2slap_w(unsigned int *I, unsigned int *J, double *W, unsigned int K, 
+                    unsigned int ** r_slap, unsigned int **J_slap,double **w_slap);
+
+int convert_ij2slap_N(unsigned int *I, unsigned int *J, unsigned int K,   
+                      unsigned int N, unsigned int ** r_slap, 
+                      unsigned int **J_slap);
+  
+
+void write_edges(FILE *fileout, unsigned int *J_slap, 
+                 unsigned int *r_slap, unsigned int N);
+  
+
+void dump_deg_distr(unsigned int *degs, double *p, int n);
+
+void dump_deg_seq(unsigned int *nodes, int N);
+
+void dump_edges(iltree_t t);
+
+FILE* openfile_or_exit(char *filename, char *mode, int exitcode);
+
+int compare_int(const void *x1, const void *x2);
+
+unsigned int find_max(unsigned int *, unsigned int);
+
+int is_neigh(unsigned int *J_slap, unsigned int *r_slap, unsigned int N, 
+             unsigned int i, unsigned int j);
+
+
+#endif /*__UTILS_H__*/