summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--dynamics/ising/Makefile11
-rwxr-xr-xdynamics/ising/iltree.c201
-rwxr-xr-xdynamics/ising/iltree.h55
-rw-r--r--dynamics/ising/multiplex_ising.c313
-rwxr-xr-xdynamics/ising/utils.c494
-rwxr-xr-xdynamics/ising/utils.h58
-rw-r--r--dynamics/randomwalks/Makefile20
-rwxr-xr-xdynamics/randomwalks/entropyrate2add.c147
-rwxr-xr-xdynamics/randomwalks/entropyrate2int.c151
-rwxr-xr-xdynamics/randomwalks/entropyrate2mult.c149
-rwxr-xr-xdynamics/randomwalks/iltree.c201
-rwxr-xr-xdynamics/randomwalks/iltree.h55
-rwxr-xr-xdynamics/randomwalks/statdistr2.c231
-rwxr-xr-xdynamics/randomwalks/utils.c494
-rwxr-xr-xdynamics/randomwalks/utils.h58
-rw-r--r--models/correlations/Makefile15
-rw-r--r--models/correlations/fit_utils.c382
-rw-r--r--models/correlations/fit_utils.h25
-rw-r--r--models/correlations/rank_utils.c217
-rw-r--r--models/correlations/rank_utils.h44
-rw-r--r--models/correlations/tune_qnn_adaptive.c208
-rw-r--r--models/correlations/tune_rho.c142
-rw-r--r--models/growth/Makefile26
-rw-r--r--models/growth/nibilab_linear_delay.c414
-rw-r--r--models/growth/nibilab_linear_delay_mix.c457
-rw-r--r--models/growth/nibilab_linear_delta.c403
-rw-r--r--models/growth/nibilab_linear_random_times.c417
-rw-r--r--models/growth/nibilab_nonlinear.c493
-rw-r--r--models/growth/node_deg_over_time.py82
-rw-r--r--models/nullmodels/model_MDM.py66
-rw-r--r--models/nullmodels/model_MSM.py65
-rw-r--r--models/nullmodels/model_hypergeometric.py56
-rw-r--r--models/nullmodels/model_layer_growth.py121
-rw-r--r--structure/activity/degs_to_activity_overlap.py36
-rw-r--r--structure/activity/degs_to_binary.py45
-rw-r--r--structure/activity/hamming_dist.py41
-rw-r--r--structure/activity/layer_activity.py27
-rw-r--r--structure/activity/layer_activity_vectors.py44
-rw-r--r--structure/activity/multiplexity.py41
-rw-r--r--structure/activity/node_activity.py51
-rw-r--r--structure/activity/node_activity_vectors.py68
-rw-r--r--structure/activity/node_degree_vectors.py68
-rw-r--r--structure/correlations/Makefile15
-rw-r--r--structure/correlations/compute_pearson.py33
-rw-r--r--structure/correlations/compute_rho.py118
-rw-r--r--structure/correlations/compute_tau.py133
-rw-r--r--structure/correlations/dump_k_q.c50
-rw-r--r--structure/correlations/fit_knn.c59
-rw-r--r--structure/correlations/fit_utils.c382
-rw-r--r--structure/correlations/fit_utils.h25
-rw-r--r--structure/correlations/knn_q_from_degrees.py49
-rw-r--r--structure/correlations/knn_q_from_layers.py98
-rw-r--r--structure/correlations/rank_nodes.py74
-rw-r--r--structure/correlations/rank_nodes_thresh.py87
-rw-r--r--structure/correlations/rank_occurrence.py45
-rw-r--r--structure/correlations/rank_utils.c217
-rw-r--r--structure/correlations/rank_utils.h44
-rw-r--r--structure/metrics/aggregate_layers_w.py37
-rw-r--r--structure/metrics/avg_edge_overlap.py47
-rw-r--r--structure/metrics/cartography_from_columns.py44
-rw-r--r--structure/metrics/cartography_from_deg_vectors.py37
-rw-r--r--structure/metrics/cartography_from_layers.py54
-rw-r--r--structure/metrics/edge_overlap.py43
-rw-r--r--structure/metrics/intersect_layers.py47
-rw-r--r--structure/metrics/overlap_degree.py47
-rw-r--r--structure/reinforcement/reinforcement.py61
66 files changed, 8538 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__*/
diff --git a/dynamics/randomwalks/Makefile b/dynamics/randomwalks/Makefile
new file mode 100644
index 0000000..8dc9763
--- /dev/null
+++ b/dynamics/randomwalks/Makefile
@@ -0,0 +1,20 @@
+CFLAGS="-O3"
+CC="gcc"
+MFLAG=-lm
+
+all: statdistr2 entropyrate2add entropyrate2mult entropyrate2int
+
+statdistr2: statdistr2.c
+ $(CC) $(CFLAGS) -o statdistr2 statdistr2.c utils.c iltree.c $(MFLAG)
+
+entropyrate2add: entropyrate2add.c
+ $(CC) $(CFLAGS) -o entropyrate2add entropyrate2add.c utils.c iltree.c $(MFLAG)
+
+entropyrate2mult: entropyrate2mult.c
+ $(CC) $(CFLAGS) -o entropyrate2mult entropyrate2mult.c utils.c iltree.c $(MFLAG)
+
+entropyrate2int: entropyrate2int.c
+ $(CC) $(CFLAGS) -o entropyrate2int entropyrate2int.c utils.c iltree.c $(MFLAG)
+
+clean:
+ rm statdistr2 entropyrate2add entropyrate2mult entropyrate2int
diff --git a/dynamics/randomwalks/entropyrate2add.c b/dynamics/randomwalks/entropyrate2add.c
new file mode 100755
index 0000000..fae4fc8
--- /dev/null
+++ b/dynamics/randomwalks/entropyrate2add.c
@@ -0,0 +1,147 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "utils.h"
+
+
+int main(int argc, char *argv[]){
+
+ if (argc < 6){
+ printf("Usage: %s <layer1> <layer2> <overlapping network> <N nodes> <b1> <b2>\n", argv[0]);
+ exit(1);
+ }
+
+
+ FILE *filein0,*filein1, *filein, *fileout;
+ unsigned int N0, K0,N1, K1, N, K;
+ unsigned int *J_slap0, *r_slap0, *J_slap1, *r_slap1, *J_slap, *r_slap;
+ double *w_slap;
+
+
+ int i, j;
+ double c_i, d_i, f_i, f_i_2;
+ double alpha = (atof(argv[5]));
+ double beta = (atof(argv[6]));
+
+
+ int ov;
+ int deg0, deg1;
+ double degM, part, f_j, f_j_2;
+ double degMrid, maxdegM = 200.0;
+ double num1, num2, den, h;
+ int number_nodes=(atoi(argv[4]));
+ double M=2.0;
+
+ filein0 = openfile_or_exit(argv[1], "r", 2);
+ read_slap(filein0, &K0, &N0, &J_slap0, &r_slap0);
+
+ filein1 = openfile_or_exit(argv[2], "r", 2);
+ read_slap(filein1, &K1, &N1, &J_slap1, &r_slap1);
+
+
+ filein = openfile_or_exit(argv[3], "r", 2);
+ read_slap_w(filein, &K, &N, &J_slap, &r_slap,&w_slap);
+
+ int r_slap0_n[N+1],r_slap1_n[N+1];
+ for (i=0; i<=N; i++) {
+ if (i<=N0) {
+ r_slap0_n[i]=r_slap0[i];
+ }
+ else {
+ r_slap0_n[i]=r_slap0[N0];
+ }
+ if (i<=N1) {
+ r_slap1_n[i]=r_slap1[i];
+ }
+ else {
+ r_slap1_n[i]=r_slap1[N1];
+ }
+
+ }
+
+
+
+ double c_i_vec[N];
+ double d_i_vec[N];
+ double f_i_vec[N];
+
+ for (i=0; i<N; i++) {
+ c_i=0;
+ d_i=0;
+
+ for (j=r_slap[i]; j<r_slap[i+1]; j++) {
+ ov = w_slap[j];
+
+ deg0=r_slap0_n[J_slap[j]+1]-r_slap0_n[J_slap[j]];
+
+ deg1=r_slap1_n[J_slap[j]+1]-r_slap1_n[J_slap[j]];
+
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+ c_i+=ov*(f_j+f_j_2);
+
+ d_i+=ov*(f_j+f_j_2)*log((ov*(f_j+f_j_2)));
+
+
+
+
+ }
+ c_i_vec[i]=c_i;
+ d_i_vec[i]=d_i;
+ deg0=r_slap0_n[i+1]-r_slap0_n[i];
+ deg1=r_slap1_n[i+1]-r_slap1_n[i];
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+
+ if (deg0>0.0000000001) {
+ f_i = pow (deg0, alpha);
+ }
+ else {
+ f_i = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_i_2 = pow (deg1, beta);
+ }
+ else {
+ f_i_2=0;
+ }
+ f_i_vec[i]=f_i+f_i_2;
+
+ }
+ num1=0;
+ num2=0;
+ den=0;
+ for (i=0; i<N; i++) {
+ if (c_i_vec[i]>0.0) {
+ num1+=f_i_vec[i]*c_i_vec[i]*log(c_i_vec[i]);
+ }
+
+
+ num2=num2+f_i_vec[i]*d_i_vec[i];
+ den=den+f_i_vec[i]*c_i_vec[i];
+
+
+
+
+ }
+
+ h=(num1-num2)/den;
+ printf("%f %f %f\n", h, alpha, beta);
+
+}
diff --git a/dynamics/randomwalks/entropyrate2int.c b/dynamics/randomwalks/entropyrate2int.c
new file mode 100755
index 0000000..757c7b7
--- /dev/null
+++ b/dynamics/randomwalks/entropyrate2int.c
@@ -0,0 +1,151 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "utils.h"
+
+
+int main(int argc, char *argv[]){
+
+ if (argc < 6){
+ printf("Usage: %s <layer1> <layer2> <overlapping network> <N nodes> <b_p> <b_o>\n", argv[0]);
+ exit(1);
+ }
+
+
+ FILE *filein0,*filein1,*filein2,*filein3,*filein4,*filein5,*filein6,*filein7,*filein8,*filein9, *filein, *fileout;
+ unsigned int N0, K0,N1, K1,N2, K2,N3, K3,N4, K4,N5, K5,N6, K6,N7, K7,N8, K8,N9, K9, N, K;
+ unsigned int *J_slap0, *r_slap0, *J_slap1, *r_slap1,*J_slap2, *r_slap2, *J_slap3, *r_slap3,*J_slap4, *r_slap4, *J_slap5, *r_slap5,*J_slap6, *r_slap6, *J_slap7, *r_slap7,*J_slap8, *r_slap8, *J_slap9, *r_slap9, *J_slap, *r_slap;
+ double *w_slap;
+
+
+ int i, j;
+ double c_i, d_i, f_i, f_i_2;
+ double alpha = (atof(argv[5]));
+ double beta = (atof(argv[6]));
+
+ int ov;
+ int deg0, deg1,deg2, deg3,deg4, deg5,deg6, deg7,deg8, deg9;
+ double degM, part, f_j, f_j_2;
+ double degMrid, maxdegM = 200.0;
+ double num1, num2, den, h;
+ int number_nodes=(atoi(argv[4]));
+ double M=2.0;
+
+ filein0 = openfile_or_exit(argv[1], "r", 2);
+ read_slap(filein0, &K0, &N0, &J_slap0, &r_slap0);
+
+ filein1 = openfile_or_exit(argv[2], "r", 2);
+ read_slap(filein1, &K1, &N1, &J_slap1, &r_slap1);
+
+ filein = openfile_or_exit(argv[3], "r", 2);
+ read_slap_w(filein, &K, &N, &J_slap, &r_slap,&w_slap);
+
+ int r_slap0_n[N+1],r_slap1_n[N+1],r_slap2_n[N+1],r_slap3_n[N+1],r_slap4_n[N+1],r_slap5_n[N+1],r_slap6_n[N+1],r_slap7_n[N+1],r_slap8_n[N+1],r_slap9_n[N+1];
+ for (i=0; i<=N; i++) {
+ if (i<=N0) {
+ r_slap0_n[i]=r_slap0[i];
+ }
+ else {
+ r_slap0_n[i]=r_slap0[N0];
+ }
+ if (i<=N1) {
+ r_slap1_n[i]=r_slap1[i];
+ }
+ else {
+ r_slap1_n[i]=r_slap1[N1];
+ }
+
+
+ }
+
+
+
+ double c_i_vec[N];
+ double d_i_vec[N];
+ double f_i_vec[N];
+
+
+ for (i=0; i<N; i++) {
+ c_i=0;
+ d_i=0;
+
+ for (j=r_slap[i]; j<r_slap[i+1]; j++) {
+ ov = w_slap[j];
+
+ deg0=r_slap0_n[J_slap[j]+1]-r_slap0_n[J_slap[j]];
+ deg1=r_slap1_n[J_slap[j]+1]-r_slap1_n[J_slap[j]];
+
+
+
+
+
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+ if (part>0.0000000001) {
+ f_j = pow (part, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ f_j_2 = pow (degM, beta);
+ c_i+=ov*f_j*f_j_2;
+ if (part>0.0000000001) {
+ d_i+=ov*f_j*f_j_2*log((ov*f_j*f_j_2));
+ }
+
+
+
+ }
+ c_i_vec[i]=c_i;
+ d_i_vec[i]=d_i;
+
+ deg0=r_slap0_n[i+1]-r_slap0_n[i];
+ deg1=r_slap1_n[i+1]-r_slap1_n[i];
+
+
+
+
+
+
+
+
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)));
+
+
+
+ if (part>0.0000000001) {
+ f_i = pow (part, alpha);
+ }
+ else {
+ f_i = 0;
+ }
+ f_i_2 = pow (degM, beta);
+
+ f_i_vec[i]=f_i*f_i_2;
+
+ }
+ num1=0;
+ num2=0;
+ den=0;
+ for (i=0; i<N; i++) {
+ if (c_i_vec[i]>0.0) {
+
+ num1+=f_i_vec[i]*c_i_vec[i]*log(c_i_vec[i]);
+ }
+
+
+ num2=num2+f_i_vec[i]*d_i_vec[i];
+ den=den+f_i_vec[i]*c_i_vec[i];
+
+
+ }
+
+ h=(num1-num2)/den*1.0;
+ printf("%f %f %f\n", h, alpha, beta);
+
+}
diff --git a/dynamics/randomwalks/entropyrate2mult.c b/dynamics/randomwalks/entropyrate2mult.c
new file mode 100755
index 0000000..cded4eb
--- /dev/null
+++ b/dynamics/randomwalks/entropyrate2mult.c
@@ -0,0 +1,149 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "utils.h"
+
+
+int main(int argc, char *argv[]){
+
+
+ if (argc < 6){
+ printf("Usage: %s <layer1> <layer2> <overlapping network> <N nodes> <b1> <b2>\n", argv[0]);
+ exit(1);
+ }
+
+
+ FILE *filein0,*filein1, *filein, *fileout;
+ unsigned int N0, K0,N1, K1, N, K;
+ unsigned int *J_slap0, *r_slap0, *J_slap1, *r_slap1, *J_slap, *r_slap;
+ double *w_slap;
+
+
+ int i, j;
+ double c_i, d_i, f_i, f_i_2;
+ double alpha = (atof(argv[5]));
+ double beta = (atof(argv[6]));
+
+
+ int ov;
+ int deg0, deg1;
+ double degM, part, f_j, f_j_2;
+ double degMrid, maxdegM = 200.0;
+ double num1, num2, den, h;
+ int number_nodes=(atoi(argv[4]));
+ double M=2.0;
+
+ filein0 = openfile_or_exit(argv[1], "r", 2);
+ read_slap(filein0, &K0, &N0, &J_slap0, &r_slap0);
+
+ filein1 = openfile_or_exit(argv[2], "r", 2);
+ read_slap(filein1, &K1, &N1, &J_slap1, &r_slap1);
+
+
+ filein = openfile_or_exit(argv[3], "r", 2);
+ read_slap_w(filein, &K, &N, &J_slap, &r_slap,&w_slap);
+
+ int r_slap0_n[N+1],r_slap1_n[N+1];
+ for (i=0; i<=N; i++) {
+ if (i<=N0) {
+ r_slap0_n[i]=r_slap0[i];
+ }
+ else {
+ r_slap0_n[i]=r_slap0[N0];
+ }
+ if (i<=N1) {
+ r_slap1_n[i]=r_slap1[i];
+ }
+ else {
+ r_slap1_n[i]=r_slap1[N1];
+ }
+
+ }
+
+
+
+ double c_i_vec[N];
+ double d_i_vec[N];
+ double f_i_vec[N];
+
+
+ for (i=0; i<N; i++) {
+ c_i=0;
+ d_i=0;
+
+ for (j=r_slap[i]; j<r_slap[i+1]; j++) {
+ ov = w_slap[j];
+
+ deg0=r_slap0_n[J_slap[j]+1]-r_slap0_n[J_slap[j]];
+
+ deg1=r_slap1_n[J_slap[j]+1]-r_slap1_n[J_slap[j]];
+
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+ c_i+=ov*f_j*f_j_2;
+ if (deg0>0.0000000001 && deg1>0.0000000001) {
+ d_i+=ov*f_j*f_j_2*log((ov*f_j*f_j_2));
+ }
+
+
+ /*chiudo il for*/
+ }
+ c_i_vec[i]=c_i;
+ d_i_vec[i]=d_i;
+ deg0=r_slap0_n[i+1]-r_slap0_n[i];
+ deg1=r_slap1_n[i+1]-r_slap1_n[i];
+
+ degM=(deg0+deg1)*1.0;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+
+ if (deg0>0.0000000001) {
+ f_i = pow (deg0, alpha);
+ }
+ else {
+ f_i = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_i_2 = pow (deg1, beta);
+ }
+ else {
+ f_i_2=0;
+ }
+ f_i_vec[i]=f_i*f_i_2;
+
+ }
+ num1=0;
+ num2=0;
+ den=0;
+ for (i=0; i<N; i++) {
+ if (c_i_vec[i]>0.0) {
+ num1+=f_i_vec[i]*c_i_vec[i]*log(c_i_vec[i]);
+ }
+
+
+ num2=num2+f_i_vec[i]*d_i_vec[i];
+ den=den+f_i_vec[i]*c_i_vec[i];
+
+
+
+
+ }
+
+ h=(num1-num2)/den;
+ printf("%f %f %f\n", h, alpha, beta);
+
+}
diff --git a/dynamics/randomwalks/iltree.c b/dynamics/randomwalks/iltree.c
new file mode 100755
index 0000000..5a693d0
--- /dev/null
+++ b/dynamics/randomwalks/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/randomwalks/iltree.h b/dynamics/randomwalks/iltree.h
new file mode 100755
index 0000000..3e835e4
--- /dev/null
+++ b/dynamics/randomwalks/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/randomwalks/statdistr2.c b/dynamics/randomwalks/statdistr2.c
new file mode 100755
index 0000000..10528da
--- /dev/null
+++ b/dynamics/randomwalks/statdistr2.c
@@ -0,0 +1,231 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "utils.h"
+
+
+int main(int argc, char *argv[]){
+
+
+ if (argc < 6){
+ printf("Usage: %s <layer1> <layer2> <overlapping network> <N nodes> <b1> <b2>\n", argv[0]);
+ exit(1);
+ }
+
+
+ /*dichiarazioni enzo*/
+ FILE *filein0,*filein1, *filein, *fileout;
+ unsigned int N0, K0,N1, K1, N, K;
+ unsigned int *J_slap0, *r_slap0, *J_slap1, *r_slap1, *J_slap, *r_slap;
+ double *w_slap;
+
+ /*dichiarazioni mie*/
+ /*r_slap e' lungo N+1, da 0 a N compresi*/
+ int i, j;
+ double f_j, f_j_2;
+ double alpha = (atof(argv[5]));
+ double beta = (atof(argv[6]));
+
+ /*il bias*/
+ int ov;
+ int deg0, deg1;
+ double degM, prodM, part, intM;
+ double degMrid, maxdegM = 200.0;
+
+ int number_nodes=(atoi(argv[4]));
+ double M=2.0;
+
+ filein0 = openfile_or_exit(argv[1], "r", 2);
+ read_slap(filein0, &K0, &N0, &J_slap0, &r_slap0);
+
+ filein1 = openfile_or_exit(argv[2], "r", 2);
+ read_slap(filein1, &K1, &N1, &J_slap1, &r_slap1);
+
+
+ filein = openfile_or_exit(argv[3], "r", 2);
+ read_slap_w(filein, &K, &N, &J_slap, &r_slap,&w_slap);
+
+ int r_slap0_n[N+1],r_slap1_n[N+1];
+ for (i=0; i<=N; i++) {
+ if (i<=N0) {
+ r_slap0_n[i]=r_slap0[i];
+ }
+ else {
+ r_slap0_n[i]=r_slap0[N0];
+ }
+ if (i<=N1) {
+ r_slap1_n[i]=r_slap1[i];
+ }
+ else {
+ r_slap1_n[i]=r_slap1[N1];
+ }
+
+ }
+
+
+
+ double cf_i_vec_add[N];
+ double cf_i_vec_mult[N];
+ double cf_i_vec_part[N];
+ double cf_i_vec_int[N];
+
+ double tot_add=0, tot_mult=0, tot_part=0, tot_int=0;
+ double c_i_add, c_i_mult, c_i_part, c_i_int;
+ double f_add, f_mult, f_int;
+ /*ciclo sui nodi dell'aggregato*/
+ for (i=0; i<N; i++) {
+ c_i_add=0;
+ c_i_mult=0;
+ c_i_part=0;
+ c_i_int=0;
+
+ /*ciclo sui primi vicini di i*/
+ for (j=r_slap[i]; j<r_slap[i+1]; j++) {
+ ov = w_slap[j];
+
+ deg0=r_slap0_n[J_slap[j]+1]-r_slap0_n[J_slap[j]];
+
+ deg1=r_slap1_n[J_slap[j]+1]-r_slap1_n[J_slap[j]];
+
+
+ degM=(deg0+deg1)*1.0;
+ //prodM=(deg0*deg1)*1.0;
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+ //intM=degM*part;
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+
+ c_i_add+=ov*(f_j+f_j_2);
+
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+
+ c_i_mult+=ov*(f_j*f_j_2);
+ //c_i_part+=ov*part;
+
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+ if (part>0.0000000001) {
+ f_j = pow (part, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ f_j_2 = pow (degM, beta);
+
+ c_i_int+=ov*(f_j*f_j_2);
+
+ /*chiudo il for*/
+ }
+
+ deg0=r_slap0_n[i+1]-r_slap0_n[i];
+ deg1=r_slap1_n[i+1]-r_slap1_n[i];
+
+ degM=(deg0+deg1)*1.0;
+ prodM=(deg0*deg1)*1.0;
+ part = (M/(M-1))* (1-(pow((deg0/degM),2))-(pow((deg1/degM),2)) );
+
+ intM=degM*part;
+
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+
+
+ cf_i_vec_add[i]=c_i_add*(f_j+f_j_2);
+
+
+ if (deg0>0.0000000001) {
+ f_j = pow (deg0, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ if (deg1>0.0000000001) {
+ f_j_2 = pow (deg1, beta);
+ }
+ else {
+ f_j_2=0;
+ }
+
+
+
+ cf_i_vec_mult[i]=c_i_mult*(f_j*f_j_2);
+ //cf_i_vec_part[i]=c_i_part*part;
+ if (part>0.0000000001) {
+ f_j = pow (part, alpha);
+ }
+ else {
+ f_j = 0;
+ }
+ f_j_2 = pow (degM, beta);
+
+
+
+ cf_i_vec_int[i]=c_i_int*(f_j*f_j_2);
+
+ tot_add+=cf_i_vec_add[i];
+ tot_mult+=cf_i_vec_mult[i];
+ tot_part+=cf_i_vec_part[i];
+ tot_int+=cf_i_vec_int[i];
+ }
+
+
+ double vec_add[N];
+ double vec_mult[N];
+ double vec_part[N];
+ double vec_int[N];
+ double tot_add_rid=0, tot_mult_rid=0, tot_part_rid=0, tot_int_rid=0;
+
+ for (i=0; i<N; i++) {
+ vec_add[i]=cf_i_vec_add[i]/tot_add;
+ vec_mult[i]=cf_i_vec_mult[i]/tot_mult;
+ vec_part[i]=cf_i_vec_part[i]/tot_part;
+ vec_int[i]=cf_i_vec_int[i]/tot_int;
+ tot_add_rid+=vec_add[i];
+ tot_mult_rid+=vec_mult[i];
+ tot_part_rid+=vec_part[i];
+ tot_int_rid+=vec_int[i];
+
+ }
+
+ //sigma delle distr
+ double average_add, variance_add, std_deviation_add, sum_add = 0, sum1_add = 0;
+ double average_mult, variance_mult, std_deviation_mult, sum_mult = 0, sum1_mult = 0;
+
+ double average_int, variance_int, std_deviation_int, sum_int = 0, sum1_int = 0;
+ double sigma_norm_add, sigma_norm_mult, sigma_norm_int;
+ for (i=0; i<N; i++) {
+ printf("%d %g %g %g %g %g\n", i, vec_add[i],vec_mult[i],vec_int[i], alpha, beta);
+
+ }
+
+}
diff --git a/dynamics/randomwalks/utils.c b/dynamics/randomwalks/utils.c
new file mode 100755
index 0000000..952fcd7
--- /dev/null
+++ b/dynamics/randomwalks/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/randomwalks/utils.h b/dynamics/randomwalks/utils.h
new file mode 100755
index 0000000..c844248
--- /dev/null
+++ b/dynamics/randomwalks/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__*/
diff --git a/models/correlations/Makefile b/models/correlations/Makefile
new file mode 100644
index 0000000..b3828aa
--- /dev/null
+++ b/models/correlations/Makefile
@@ -0,0 +1,15 @@
+CFLAGS="-O3"
+CC="gcc"
+GSL_FLAGS=-lgsl -lgslcblas -lm
+MFLAG=-lm
+
+all: tune_qnn_adaptive tune_rho
+
+tune_qnn_adaptive: tune_qnn_adaptive.c rank_utils.c fit_utils.c
+ $(CC) $(CFLAGS) -o tune_qnn_adaptive tune_qnn_adaptive.c rank_utils.c fit_utils.c $(GSL_FLAGS)
+
+tune_rho: tune_rho.c rank_utils.c
+ $(CC) $(CFLAGS) -o tune_rho tune_rho.c rank_utils.c $(MFLAG)
+
+clean:
+ rm tune_rho tune_qnn_adaptive
diff --git a/models/correlations/fit_utils.c b/models/correlations/fit_utils.c
new file mode 100644
index 0000000..e9e3954
--- /dev/null
+++ b/models/correlations/fit_utils.c
@@ -0,0 +1,382 @@
+/**
+ *
+ * Fit a given sequence with a power-law function
+ *
+ * a*x^{b}
+ *
+ * The fit is actually performed as a linear fit on the
+ * exponential-binned log-log distribution
+ *
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <gsl/gsl_fit.h>
+
+/**
+ *
+ * Load a sequence from a file, which contains one element on each line
+ *
+ */
+
+void load_sequence(char *fname, double **v, int *N){
+
+ int size;
+ char buff[256];
+
+ FILE *f;
+
+ f = fopen(fname, "r");
+ if (!f){
+ fprintf(stderr, "Error opening file %s!!!! Exiting...\n", fname);
+ exit(1);
+ }
+
+ *N =0;
+ size = 10;
+ if (*v == NULL)
+ *v = malloc(size * sizeof(double));
+ else{
+ *v = realloc(*v, size * sizeof(double));
+ }
+
+ while (fgets(buff, 255, f)){
+ (*v)[*N] = atof(buff);
+ *N += 1;
+ if (*N == size){
+ size += 10;
+ *v = realloc(*v, size * sizeof(double));
+ }
+ }
+ *v = realloc(*v, (*N) * sizeof(double));
+ fclose(f);
+}
+
+
+/**
+ *
+ * Load a sequence, getting the "col"-th column of the input file
+ *
+ */
+
+void load_sequence_col(char *fname, double **v, int *N, int col){
+
+ int size, n;
+ char buff[256];
+ char *ptr;
+
+ FILE *f;
+
+ f = fopen(fname, "r");
+ if (!f){
+ fprintf(stderr, "Error opening file %s!!!! Exiting...\n", fname);
+ exit(1);
+ }
+
+ *N =0;
+ size = 10;
+ if (*v == NULL)
+ *v = malloc(size * sizeof(double));
+ else{
+ *v = realloc(*v, size * sizeof(double));
+ }
+
+ while (fgets(buff, 255, f)){
+ ptr = strtok(buff, " ");
+ if (col > 0){
+ n = 0;
+ while (n<col){
+ ptr = strtok(NULL, " ");
+ n += 1;
+ }
+ }
+ (*v)[*N] = atof(ptr);
+ *N += 1;
+
+ if (*N == size){
+ size += 10;
+ *v = realloc(*v, size * sizeof(double));
+ }
+ }
+ *v = realloc(*v, (*N) * sizeof(double));
+ fclose(f);
+}
+
+
+/**
+ *
+ * Make the exponential binning of a distribution, with a giving
+ * exponent alpha. The value of y[i] is the number of elements of v
+ * whose value lies between x[i-1] and x[i]...
+ *
+ */
+
+void exp_bin_cnt(double *v, int N, double alpha, double **x, double **y, int *num){
+
+ double min_v, max_v, val, last_val;
+ int i, j, size, num_x;
+ double last_size, new_size;
+
+ min_v = max_v = v[0];
+
+ for (i=1; i<N; i ++){
+ if (v[i] > max_v)
+ max_v = v[i];
+ else if (v[i] > 0 && v[i] < min_v)
+ min_v = v[i];
+ }
+
+ size = 10;
+ if (*x == NULL){
+ *x = malloc(size * sizeof(double));
+ }
+ else{
+ *x = realloc(*x, size * sizeof(double));
+ }
+
+ val = min_v;
+ last_size = min_v;
+ (*x)[0] = min_v;
+ num_x = 1;
+
+ while(val < max_v){
+ new_size = last_size * alpha;
+ val = last_size + new_size;
+ last_size = new_size;
+ last_val = val;
+ (*x)[num_x] = val;
+ num_x += 1;
+ if (num_x == size){
+ size += 10;
+ *x = realloc(*x, size * sizeof(double));
+ }
+ }
+
+ if (*y == NULL){
+ *y = malloc(num_x * sizeof(double));
+ }
+ else{
+ *y = realloc(*y, num_x * sizeof(double));
+ }
+ for (i=0; i < num_x; i++){
+ (*y)[i] = 0;
+ }
+
+
+
+ for(i=0; i <N; i ++){
+ j = 0;
+ while(v[i] > (*x)[j]){
+ j ++;
+ }
+ (*y)[j] += 1;
+ }
+ *num = num_x;
+}
+
+/**
+ *
+ * Make the exponential binning of a distribution, with a giving
+ * exponent alpha. The value of y[i] is the average of the values in
+ * the vector "w" whose corresponding v lies between x[i-1] and x[i]...
+ *
+ */
+
+
+void exp_bin_avg(double *v, double *w, int N, double alpha, double **x, double **y, int *num){
+
+ double min_v, max_v, val, last_val;
+ int i, j, size, num_x;
+ double last_size, new_size;
+ int *cnt;
+
+
+ min_v = max_v = v[0];
+
+ for (i=1; i<N; i ++){
+ if (v[i] > max_v)
+ max_v = v[i];
+ else if (v[i] > 0 && v[i] < min_v)
+ min_v = v[i];
+ }
+
+ size = 10;
+ if (*x == NULL){
+ *x = malloc(size * sizeof(double));
+ }
+ else{
+ *x = realloc(*x, size * sizeof(double));
+ }
+
+ val = min_v;
+ last_size = min_v;
+ (*x)[0] = min_v;
+ num_x = 1;
+
+ while(val < max_v){
+ new_size = last_size * alpha;
+ val = last_size + new_size;
+ last_size = new_size;
+ last_val = val;
+ (*x)[num_x] = val;
+ num_x += 1;
+ if (num_x == size){
+ size += 10;
+ *x = realloc(*x, size * sizeof(double));
+ }
+ }
+
+
+ cnt = malloc(num_x * sizeof(int));
+
+ if (*y == NULL){
+ *y = malloc(num_x * sizeof(double));
+ }
+ else{
+ *y = realloc(*y, num_x * sizeof(double));
+ }
+ for (i=0; i < num_x; i++){
+ (*y)[i] = 0;
+ cnt[i] = 0;
+ }
+
+ for(i=0; i <N; i ++){
+ j = 0;
+ while(j < num_x && v[i] > (*x)[j]){
+ j ++;
+ }
+ if(j == num_x){
+ printf("Error!!!!! Trying to assing a non-existing bin!!! -- fit_utils.c:exp_bin_avg!!!\n");
+ exit(37);
+ }
+ (*y)[j] += w[i];
+ cnt[j] += 1;
+ }
+ *num = num_x;
+
+ for(i=0; i< num_x; i++){
+ if (cnt[i] > 0){
+ (*y)[i] = (*y)[i] / cnt[i];
+ }
+ }
+ free(cnt);
+}
+
+
+/**
+ *
+ * Print a distribution on stdout
+ *
+ */
+
+void dump_distr(double *x, double *y, int N){
+ int i;
+
+ for(i=0; i<N; i ++){
+ printf("%g %g\n", x[i], y[i]);
+ }
+
+}
+
+
+/**
+ * Compact a distribution, leaving only the pairs (x_i, y_i) for which
+ * y_i > 0
+ *
+ */
+
+void compact_distr(double *x, double *y, int *num){
+
+ int i, j;
+
+ i = j = 0;
+ while(j < *num){
+ while(j < *num && y[j] == 0){
+ j ++;
+ }
+ if (j==*num){
+ break;
+ }
+ x[i] = x[j];
+ y[i] = y[j];
+ i ++;
+ j ++;
+ }
+ *num = i;
+}
+
+
+/**
+ *
+ * Apply the function "f" on all the elemnts of a vector, in-place
+ *
+ */
+
+void map_vec(double *v, int N, double (*f)(double)){
+ int i;
+
+ for (i=0; i<N; i ++){
+ v[i] = f(v[i]);
+ }
+}
+
+
+/**
+ *
+ * Normalize a distribution, dividing each y[i] for the width of the
+ * corresponding bin (i.e., x[i] - x[i-1])
+ *
+ */
+void normalize_distr(double *x, double *y, int num){
+
+ int i;
+
+ for(i=1; i<num; i ++){
+ y[i] /= (x[i] - x[i-1]);
+ }
+}
+
+/**
+ *
+ * take two vectors (k and q) and a pairing, and compute the best
+ * power-law fit "a*k^{mu}" of qnn(k)
+ *
+ */
+
+void fit_current_trend(double *k, double *q, int N, int *pairing, double *mu, double *a,
+ double *corr){
+
+ static int *num;
+ static double *q_pair, *x, *y;
+
+ int i;
+ int fit_num;
+ double c0, c1, cov00, cov01, cov11, sqsum;
+
+ if (q_pair == NULL){
+ q_pair = malloc(N * sizeof(double));
+ }
+
+ for(i=0; i<N; i++){
+ q_pair[i] = q[pairing[i]];
+ }
+
+ exp_bin_avg(k, q_pair, N, 1.3, &x, &y, &fit_num);
+ //normalize_distr(x,y, fit_num);
+ compact_distr(x, y, &fit_num);
+
+ map_vec(x, fit_num, log);
+ map_vec(y, fit_num, log);
+ gsl_fit_linear(x, 1, y, 1, fit_num, &c0, &c1, &cov00, &cov01, &cov11, &sqsum);
+
+ //fprintf(stderr,"cov00: %g cov01: %g cov11: %g\n", cov00, cov01, cov11);
+ //fprintf(stderr, "corr: %g ", cov01 / sqrt(cov00 * cov11));
+ *a = c0;
+ *mu = c1;
+ *corr = cov01/sqrt(cov00 * cov11);
+}
+
diff --git a/models/correlations/fit_utils.h b/models/correlations/fit_utils.h
new file mode 100644
index 0000000..183f47c
--- /dev/null
+++ b/models/correlations/fit_utils.h
@@ -0,0 +1,25 @@
+#ifndef __FIT_UTILS_H__
+#define __FIT_UTILS_H__
+
+
+void load_sequence(char *fname, double **v, int *N);
+
+void load_sequence_col(char *fname, double **v, int *N, int col);
+
+void exp_bin_cnt(double *v, int N, double alpha, double **x, double **y, int *num);
+
+void exp_bin_avg(double *v, double *w, int N, double alpha, double **x, double **y, int *num);
+
+void dump_distr(double *x, double *y, int N);
+
+void compact_distr(double *x, double *y, int *num);
+
+void map_vec(double *v, int N, double (*f)(double));
+
+void normalize_distr(double *x, double *y, int num);
+
+void fit_current_trend(double *k, double *q, int N, int *pairing,
+ double *mu, double *a, double *corr);
+
+
+#endif /* __FIT_UTILS_H__ */
diff --git a/models/correlations/rank_utils.c b/models/correlations/rank_utils.c
new file mode 100644
index 0000000..4b8ef41
--- /dev/null
+++ b/models/correlations/rank_utils.c
@@ -0,0 +1,217 @@
+/**
+ *
+ * Some utility functions to be used with tune_rho, compute_rho and
+ * the like
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+
+#include "rank_utils.h"
+
+void load_ranking(char *fname, int *N, double **R){
+
+ char buff[256];
+ int size = 10;
+ FILE *f;
+
+ if (*R == NULL){
+ *R = malloc(size * sizeof (double));
+ }
+ f = fopen(fname, "r");
+ if (!f){
+ printf("Unable to open file: %s!!! Exiting\n");
+ exit(1);
+ }
+
+ *N = 0;
+
+ while (fgets(buff, 255, f)){
+ if (* N == size){
+ size += 10;
+ *R = realloc(*R, size * sizeof(double));
+ }
+ (*R)[*N] = atof(buff);
+ *N += 1;
+ }
+}
+
+double avg_array(double *v, int N){
+
+ double sum = 0.0;
+ int i;
+
+ for (i = 0; i < N; i ++){
+ sum += v[i];
+ }
+ return sum/N;
+}
+
+double compute_C(double *R1, double *R2, int N){
+ double mu1, mu2, sum1, sum2;
+
+ mu1 = avg_array(R1, N);
+ mu2 = avg_array(R2, N);
+
+ sum1 = mu1 * N;
+ sum2 = mu2 * N;
+
+ return N * mu1 * mu2 - mu2 * sum1 - mu1 * sum2;
+}
+
+
+double compute_D(double *R1, double *R2, int N){
+
+ double mu1, mu2, s1, s2;
+ int i;
+
+ mu1 = avg_array(R1, N);
+ mu2 = avg_array(R1, N);
+
+ s1 = s2 = 0.0;
+
+ for (i=0 ; i < N; i ++){
+ s1 += pow((R1[i] - mu1), 2);
+ s2 += pow((R2[i] - mu2), 2);
+ }
+
+ return sqrt(s1 * s2);
+}
+
+
+double compute_rho(double *R1, double *R2, int N, int *pairing){
+
+ double rho = 0;
+ int i;
+
+ for (i=0; i < N; i ++){
+ rho += R1[i] * R2[ pairing[i] ];
+ }
+
+ rho = (rho + compute_C(R1, R2, N))/ compute_D(R1, R2, N);
+ return rho;
+}
+
+void dump_ranking(double *R, int N){
+ int i;
+
+ for (i=0; i < N; i ++){
+ printf("%d: %f\n", i, R[i] );
+ }
+}
+
+
+void init_pairing_natural(int *pairing, int N){
+ int i;
+
+ for (i = 0; i< N; i ++){
+ pairing[i] = i;
+ }
+}
+
+void init_pairing_inverse(int *pairing, int N){
+ int i;
+
+ for (i = 0; i< N; i ++){
+ pairing[i] = N-i-1;
+ }
+}
+
+void select_pairing(int *pairing, int N, int argc, char *argv[], int pos){
+
+ if (argc < pos + 1 || !strncasecmp("rnd", argv[pos], 3)){
+ init_pairing_random(pairing, N);
+ }
+ else if (!strncasecmp("nat", argv[pos], 3)){
+ init_pairing_natural(pairing, N);
+ }
+ else if (!strncasecmp("inv", argv[pos], 3)){
+ init_pairing_inverse(pairing, N);
+ }
+ else{
+ printf ("Pairing strategy \"%s\" unknown!!! Exiting...\n", argv[pos]);
+ exit(1);
+ }
+
+}
+
+
+void shuffle_sequence(int *s, int N){
+
+ int i, j, tmp;
+
+ for (i=N-1; i>=0; i--){
+ j = rand() % (i+1);
+ tmp = s[j];
+ s[j] = s[i];
+ s[i] = tmp;
+ }
+}
+
+
+
+void init_pairing_random(int *pairing, int N){
+
+ init_pairing_natural(pairing, N);
+ shuffle_sequence(pairing, N);
+
+}
+
+/* Loads a pairing from a file, in the format:
+ *
+ * rank1 rank2
+ * ...........
+ */
+void load_pairing(int **pairing, int N, char *fname){
+
+ FILE *f;
+ int i, j, num;
+ char buff[256];
+ char *ptr;
+
+ f = fopen(fname, "r");
+ if (!f){
+ printf("Error opening file \"%s\"!!!! Exiting....\n", fname);
+ exit(2);
+ }
+
+ if (*pairing == NULL){
+ *pairing = malloc(N * sizeof(int));
+ init_pairing_natural(*pairing, N);
+ }
+
+ num = 0;
+ while(num < N){
+ fgets(buff, 255, f);
+ if (buff[0] == '#')
+ continue;
+ ptr = strtok(buff, " "); /* read the first node */
+ i = atoi(ptr);
+ ptr = strtok(NULL, " "); /* read the second node */
+ j = atoi(ptr);
+ (*pairing)[i] = j;
+ num += 1;
+ }
+}
+
+
+void dump_pairing(int *pairing, int N){
+ int i;
+
+ for(i=0; i< N; i ++){
+ printf("%d %d\n", i, pairing[i]);
+ }
+}
+
+void copy_pairing(int *p1, int *p2, int N){
+ int i;
+
+ for (i=0 ; i < N; i ++){
+ p2[i] = p1[i];
+ }
+}
diff --git a/models/correlations/rank_utils.h b/models/correlations/rank_utils.h
new file mode 100644
index 0000000..521582a
--- /dev/null
+++ b/models/correlations/rank_utils.h
@@ -0,0 +1,44 @@
+#ifndef __RANK_UTILS_H__
+#define __RANK_UTILS_H__
+
+/* Load a ranking */
+void load_ranking(char *fname, int *N, double **R);
+
+/* Compute the average of an array */
+double avg_array(double *v, int N);
+
+/* Compute the term "C" in the rho correlation coefficient */
+double compute_C(double *R1, double *R2, int N);
+
+/* Compute the term "D" in the rho correlation coefficient */
+double compute_D(double *R1, double *R2, int N);
+
+/* Compute the Spearman's rank correlation coefficient \rho */
+double compute_rho(double *R1, double *R2, int N, int *pairing);
+
+void dump_ranking(double *R, int N);
+
+void init_pairing_natural(int *pairing, int N);
+
+void init_pairing_inverse(int *pairing, int N);
+
+void init_pairing_random(int *pairing, int N);
+
+void select_pairing(int *pairing, int N, int argc, char *argv[], int pos);
+
+void shuffle_sequence(int *s, int N);
+
+void dump_pairing(int *pairing, int N);
+
+void copy_pairing(int *pairing1, int *pairing2, int N);
+
+
+#endif /* __RANK_UTILS_H__ */
+
+
+
+
+
+
+
+
diff --git a/models/correlations/tune_qnn_adaptive.c b/models/correlations/tune_qnn_adaptive.c
new file mode 100644
index 0000000..71643a1
--- /dev/null
+++ b/models/correlations/tune_qnn_adaptive.c
@@ -0,0 +1,208 @@
+/**
+ *
+ * Tune the inter-layer degree correlation function \bar{q}(k) of a
+ * two-layer multiplex, in order to make it as similar as possible to
+ * a power law with exponent $\mu$, where $\mu$ is given as input
+ *
+ * This version of the program is (or, better, "shoud be") able to
+ * tune also the value of the pre-factor "a"
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <strings.h>
+#include <time.h>
+
+#include "rank_utils.h"
+
+
+
+
+inline double compute_delta(double q, double k, double mu, double a){
+
+ return fabs(log(q) - mu * log(k) - log(a));
+ //return fabs (q - a*pow(k,mu));
+}
+
+
+
+
+
+
+void tune_qnn_adaptive(double *R1, double *R2, int N, int *pairing, double eps,
+ double beta, double mu_target){
+
+ double act_mu, test_mu, F, F_new, F_old;
+ double delta1_old, delta2_old, delta1_new, delta2_new;
+ double val, prob;
+ double mu, a, err, dummy_a;
+ int *new_pairing;
+ int p1, p2, tmp_val;
+ int tot;
+ char swap = 0;
+
+ new_pairing = malloc(N * sizeof(int));
+ copy_pairing(pairing, new_pairing, N);
+
+ a = 1.0;
+
+ F = 10000;
+
+ //fprintf("%f %f %f %f %f\n", eps, beta, mu_target, act_mu, F);
+
+ tot = 0;
+ while (F > eps){
+ /* sample two positions of "pairing" and shuffle them in "new_pairing" */
+
+ p1 = rand() % N;
+ p2 = rand() % N;
+ while (p2 == p1){
+ p2 = rand() % N;
+ }
+ tmp_val = new_pairing[p1];
+ new_pairing[p1] = new_pairing[p2];
+ new_pairing[p2] = tmp_val;
+
+ delta1_old = compute_delta(R2[pairing[p1]], R1[p1], mu_target, a);
+ delta2_old = compute_delta(R2[pairing[p2]], R1[p2], mu_target, a);
+
+ delta1_new = compute_delta(R2[new_pairing[p1]], R1[p1], mu_target, a);
+ delta2_new = compute_delta(R2[new_pairing[p2]], R1[p2], mu_target, a);
+
+ //F_new = log(delta1_new * delta2_new);
+ //F_old = log(delta1_old * delta2_old);
+
+ F_new = delta1_new + delta2_new;
+ F_old = delta1_old + delta2_old;
+
+
+ if (F_new <= F_old){/* Accept this swap with probability 1 */
+ tmp_val = pairing[p1];
+ pairing[p1] = pairing[p2];
+ pairing[p2] = tmp_val;
+ //act_mu = test_mu;
+ swap = 1;
+ }
+ else{/* Accept the swap with a certain probability */
+ val = 1.0 * rand() / RAND_MAX;
+
+ //prob = pow(fabs(F - F_new)/(F+F_new), beta);
+ prob = exp(-(F_new - F_old)/beta);
+ //fprintf(stderr, "-- %g %g %g -> %f \n", F_new, F_old, F_new - F_old, prob);
+ if (val < prob){ /* Accept the swap */
+ tmp_val = pairing[p1];
+ pairing[p1] = pairing[p2];
+ pairing[p2] = tmp_val;
+ //act_mu = test_mu;
+ swap = 1;
+ }
+ else{ /* Rollback the swap */
+ tmp_val = new_pairing[p1];
+ new_pairing[p1] = new_pairing[p2];
+ new_pairing[p2] = tmp_val;
+ swap = 0;
+ }
+
+ }
+
+
+ ///F = fabs(act_mu - mu_target);
+ ///if (tot % 200 == 0){
+ //fprintf(stderr, "%d %g\n", tot, act_mu);
+
+ //fprintf(stderr, "%d: %f %f ---- %f %f ---- %f %f ---- %d \n",
+ //tot, delta1_old, delta2_old, delta1_new, delta2_new, F_old, F_new, swap);
+
+ ///}
+ tot += 1;
+ if (tot %200 == 0){
+ fit_current_trend(R1, R2, N, pairing, &mu, &a, &err);
+ fprintf(stderr, "mu: %g a: %g corr: %g\n", mu, a, err);
+ //a = a - 0.01 *(a - exp(a));
+ a = exp(a);
+ }
+ fit_current_trend(R1, R2, N, pairing, &mu, &dummy_a, &err);
+ F = fabs(mu - mu_target);
+ }
+}
+
+void dump_qnn(double *R1, double *R2, int N, int *pairing){
+
+ int i;
+ int *qnn, *num;
+
+ qnn = malloc(N * sizeof(int));
+ num = malloc(N * sizeof(int));
+
+ for (i=0; i<N; i++){
+ qnn[i]=0;
+ num[i]=0;
+ }
+
+ for (i=0; i<N; i++){
+ qnn[(int)(R1[i])] += R2[pairing[i]];
+ num[(int)(R1[i])] += 1;
+ }
+
+ for(i=0; i<N; i++){
+ if (num[i] >0){
+ printf("%d %f\n", i, 1.0*qnn[i]/num[i]);
+ }
+ }
+ free(num);
+ free(qnn);
+}
+
+
+void dump_degs(double *R1, double *R2, int N, int *pairing){
+
+ int i;
+
+ for(i=0; i<N; i++){
+ printf("%g %g\n", R1[i], R2[pairing[i]]);
+ }
+
+}
+
+int main (int argc, char *argv[]){
+
+ int N1, N2;
+ double *R1 = NULL;
+ double *R2 = NULL;
+ double eps, beta, mu_target;
+
+ int *pairing = NULL;
+
+ srand(time(NULL));
+
+ if (argc < 6){
+ printf("Usage: %s <degs1> <degs2> <mu> <eps> <beta> [RND|NAT|INV]\n", argv[0]);
+ exit(1);
+ }
+
+ mu_target = atof(argv[3]);
+ eps = atof(argv[4]);
+ beta = atof(argv[5]);
+
+ load_ranking(argv[1], &N1, &R1);
+ load_ranking(argv[2], &N2, &R2);
+
+ if (N1 != N2){
+ printf("Error!!!! The two files must have the same number of lines!!!! Exiting...\n");
+ exit(1);
+ }
+
+ pairing = malloc(N1 * sizeof(int));
+
+ select_pairing(pairing, N1, argc, argv, 6);
+
+ tune_qnn_adaptive(R1, R2, N1, pairing, eps, beta, mu_target);
+
+ dump_pairing(pairing, N1);
+ //dump_qnn(R1, R2, N1, pairing);
+ //dump_degs(R1, R2, N1, pairing);
+
+}
diff --git a/models/correlations/tune_rho.c b/models/correlations/tune_rho.c
new file mode 100644
index 0000000..4495b0e
--- /dev/null
+++ b/models/correlations/tune_rho.c
@@ -0,0 +1,142 @@
+/**
+ *
+ * Tune the Spearman's \rho correlation coefficient between two rankings
+ * given as input, using a simulated-anneling procedure
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <strings.h>
+#include <time.h>
+
+#include "rank_utils.h"
+
+//void select_pairing(int *pairing, int N, int argc, char *argv[]){
+//
+// if (argc < 7 || !strncasecmp("rnd", argv[6], 3)){
+// init_pairing_random(pairing, N);
+// }
+// else if (!strncasecmp("nat", argv[6], 3)){
+// init_pairing_natural(pairing, N);
+// }
+// else if (!strncasecmp("inv", argv[6], 3)){
+// init_pairing_inverse(pairing, N);
+// }
+// else{
+// printf ("Pairing strategy \"%s\" unknown!!! Exiting...\n", argv[6]);
+// exit(1);
+// }
+//
+//}
+
+void tune_rho(double *R1, double *R2, int N, int *pairing, double eps,
+ double beta, double rho_target){
+
+ double act_rho, test_rho, F, F_new;
+ double val, prob;
+ int *new_pairing;
+ int p1, p2, tmp_val;
+ int tot;
+
+ new_pairing = malloc(N * sizeof(int));
+ copy_pairing(pairing, new_pairing, N);
+
+ act_rho = compute_rho(R1, R2, N, pairing);
+
+ F = fabs(act_rho - rho_target);
+
+ //fprintf("%f %f %f %f %f\n", eps, beta, rho_target, act_rho, F);
+
+ tot = 0;
+ while (F > eps){
+ /* sample two positions of "pairing" and shuffle them in "new_pairing" */
+ p1 = rand() % N;
+ p2 = rand() % N;
+ while (p2 == p1){
+ p2 = rand() % N;
+ }
+ tmp_val = new_pairing[p1];
+ new_pairing[p1] = new_pairing[p2];
+ new_pairing[p2] = tmp_val;
+ test_rho = compute_rho(R1, R2, N, new_pairing);
+ F_new = fabs(test_rho - rho_target);
+ if (F_new < F){/* Accept this swap with probability 1 */
+ tmp_val = pairing[p1];
+ pairing[p1] = pairing[p2];
+ pairing[p2] = tmp_val;
+ act_rho = test_rho;
+ }
+ else{/* Accept the swap with a certain probability */
+ val = 1.0 * rand() / RAND_MAX;
+
+ //prob = pow(fabs(F - F_new)/(F+F_new), beta);
+ prob = exp(-(F_new - F)/beta);
+ //fprintf(stderr, "-- %f\n", prob);
+ if (val < prob){ /* Accept the swap */
+ tmp_val = pairing[p1];
+ pairing[p1] = pairing[p2];
+ pairing[p2] = tmp_val;
+ act_rho = test_rho;
+ }
+ else{ /* Rollback the swap */
+ tmp_val = new_pairing[p1];
+ new_pairing[p1] = new_pairing[p2];
+ new_pairing[p2] = tmp_val;
+ }
+
+ }
+ F = fabs(act_rho - rho_target);
+ if (tot % 200 == 0){
+ fprintf(stderr, "%d %g\n", tot, act_rho);
+ }
+ tot += 1;
+ }
+}
+
+
+
+
+int main (int argc, char *argv[]){
+
+ int N1, N2;
+ double *R1 = NULL;
+ double *R2 = NULL;
+ double eps, beta, rho, rho_target;
+
+ int *pairing = NULL;
+
+ srand(time(NULL));
+
+ if (argc < 6){
+ printf("Usage: %s <rank1> <rank2> <rho> <eps> <beta> [RND|NAT|INV]\n", argv[0]);
+ exit(1);
+ }
+
+ rho_target = atof(argv[3]);
+ eps = atof(argv[4]);
+ beta = atof(argv[5]);
+
+ load_ranking(argv[1], &N1, &R1);
+ load_ranking(argv[2], &N2, &R2);
+
+ if (N1 != N2){
+ printf("Error!!!! The two files must have the same number of lines!!!! Exiting...\n");
+ exit(1);
+ }
+
+ pairing = malloc(N1 * sizeof(int));
+
+ select_pairing(pairing, N1, argc, argv, 6);
+
+ rho = compute_rho(R1, R2, N1, pairing);
+
+ fprintf(stderr, "%g\n", rho);
+
+ tune_rho(R1, R2, N1, pairing, eps, beta, rho_target);
+
+ dump_pairing(pairing, N1);
+
+}
diff --git a/models/growth/Makefile b/models/growth/Makefile
new file mode 100644
index 0000000..d0b461d
--- /dev/null
+++ b/models/growth/Makefile
@@ -0,0 +1,26 @@
+CFLAGS="-O3"
+CC="gcc"
+MFLAG=-lm
+
+all: nibilab_linear_delta nibilab_linear_delay nibilab_linear_delay_mix \
+ nibilab_linear_random_times nibilab_nonlinear
+
+nibilab_linear_delta: nibilab_linear_delta.c
+ $(CC) $(CFLAGS) -o nibilab_linear_delta nibilab_linear_delta.c $(MFLAG)
+
+nibilab_linear_delay: nibilab_linear_delay.c
+ $(CC) $(CFLAGS) -o nibilab_linear_delay nibilab_linear_delay.c $(MFLAG)
+
+nibilab_linear_delay_mix: nibilab_linear_delay_mix.c
+ $(CC) $(CFLAGS) -o nibilab_linear_delay_mix nibilab_linear_delay_mix.c $(MFLAG)
+
+nibilab_linear_random_times: nibilab_linear_random_times.c
+ $(CC) $(CFLAGS) -o nibilab_linear_random_times nibilab_linear_random_times.c $(MFLAG)
+
+nibilab_nonlinear: nibilab_nonlinear.c
+ $(CC) $(CFLAGS) -o nibilab_nonlinear nibilab_nonlinear.c $(MFLAG)
+
+
+clean:
+ rm nibilab_linear_delta nibilab_linear_delay nibilab_linear_delay_mix \
+ nibilab_linear_random_times nibilab_nonlinear
diff --git a/models/growth/nibilab_linear_delay.c b/models/growth/nibilab_linear_delay.c
new file mode 100644
index 0000000..518b09d
--- /dev/null
+++ b/models/growth/nibilab_linear_delay.c
@@ -0,0 +1,414 @@
+/**
+ *
+ * NiBiLaB model (Nicosia, Bianconi, Latora, Barthelemy) of multiplex
+ * linear preferential attachment
+ *
+ * A new node arrives on layer 1 at each time step, but its replica on
+ * layer 2 arrives after a power-law distributed delay \tau
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+
+typedef struct{
+ int size;
+ int N;
+ int *val;
+} int_array_t;
+
+
+typedef struct{
+ int K;
+ int N;
+ int size;
+ int *i;
+ int *j;
+ int *degrees;
+ int *arrived;
+ double *cumul;
+ int_array_t *times;
+} ij_net_t;
+
+
+
+typedef struct{
+ double a;
+ double b;
+ double c;
+ double d;
+} coupling_t;
+
+
+typedef struct{
+ int N;
+ double x_min;
+ double *distr;
+ double gamma;
+} distr_t;
+
+
+int init_network(ij_net_t *G, int m0){
+
+ int n;
+
+ for(n=0; n<m0; n++){
+ G->i[n] = n;
+ G->j[n] = n+1 % m0;
+ G->degrees[n] = 2;
+ G->arrived[n] = n;
+ }
+ G->K = m0;
+ G->N = m0;
+ return n;
+}
+
+
+
+void dump_network(ij_net_t *G){
+
+ int i;
+ for(i=0; i< G->K; i ++){
+ printf("%d %d\n",G->i[i], G->j[i]);
+ }
+}
+
+
+void dump_network_to_file(ij_net_t *G, char *fname){
+
+
+ FILE *f;
+ int i;
+
+
+ f = fopen(fname, "w+");
+
+ for(i=0; i< G->K; i ++){
+ fprintf(f, "%d %d\n",G->i[i], G->j[i]);
+ }
+ fclose(f);
+}
+
+
+
+double f1(double v1, double v2, coupling_t *matr){
+ return v1 * matr->a + v2 * matr->b;
+}
+
+double f2(double v1, double v2, coupling_t *matr){
+ return v1 * matr->c + v2 * matr->d;
+}
+
+
+void compute_cumul(ij_net_t *G1, ij_net_t *G2, coupling_t *matr){
+
+ int i;
+ double val;
+
+ G1->cumul[0] = f1(G1->degrees[ G1->arrived[0] ], G2->degrees[ G1->arrived[0] ], matr) ;
+ G2->cumul[0] = f2(G1->degrees[ G2->arrived[0] ], G2->degrees[ G2->arrived[0] ], matr) ;
+
+ for(i=1; i<G1->N; i ++){
+ val = f1(G1->degrees[ G1->arrived[i] ], G2->degrees[ G1->arrived[i] ],matr) ;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ val = f2(G1->degrees[ G2->arrived[i] ], G2->degrees[ G2->arrived[i] ], matr) ;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+void dump_cumul(ij_net_t *G){
+
+ int i;
+
+ for (i=0; i<G->N; i ++){
+ printf("%d %2.6f\n", G->times[i], G->cumul[i]);
+ }
+
+}
+
+
+
+int select_neigh_cumul(ij_net_t *G){
+
+ double r;
+ int j;
+
+ r = (rand() * 1.0 / RAND_MAX) * G->cumul[ G->N-1 ];
+ //printf(" r: %f ", r);
+
+ j = 0;
+ /* Change with a binary search ???!?!?!?! */
+ while (r > G->cumul[j]){
+ j ++;
+ }
+ return G->arrived[j];
+}
+
+int already_neighbour(ij_net_t *G, int j, int d, int offset){
+
+ int i;
+
+ for(i=G-> K + offset; i< G->K + offset + j; i ++){
+ if (G->j[i] == d)
+ return 1;
+ }
+ return 0;
+}
+
+
+void sample_edges(ij_net_t *G, int n, int m, int offset){
+
+ int j;
+ int d;
+
+ //printf("----");
+ for(j=0; j<m; j++){
+ G->i[G->K + offset + j] = n;
+ d = select_neigh_cumul(G);
+ while(already_neighbour(G, j, d, offset)){
+ d = select_neigh_cumul(G);
+ }
+ //printf(" link %d--%d\n", n, d);
+ G->j[G->K + offset + j] = d;
+ G->degrees[d] += 1;
+ }
+ //printf("\n");
+}
+
+void create_distr(double *v, double gamma, int x_min, int x_max){
+ int n, i;
+ double sum, new_sum;
+
+ n = x_max-x_min + 1;
+
+ sum = 0;
+ for(i=0; i<n; i++){
+ v[i] = pow((x_min + i), gamma);
+ sum += v[i];
+ }
+ new_sum = 0;
+ /* Now we normalize the array*/
+ for(i=0; i<n; i++){
+ v[i]/= sum;
+ new_sum += v[i];
+ }
+ //printf("New_sum: %lf\n", new_sum);
+ /* Now we compute the cumulative distribution*/
+ for(i=1; i<n; i++){
+ v[i] += v[i-1];
+ }
+}
+
+
+inline int find_degree(double *v, int dim, double xi){
+ int i;
+
+ i=0;
+ while(xi > v[i])
+ i++;
+ return i;
+
+}
+
+
+int sample_power_law(distr_t *d){
+
+ double xi, q, val;
+ int k;
+
+ while(1){
+ xi = rand()* 1.0 / RAND_MAX;
+ k = find_degree(d->distr, d->N, xi);
+ val = rand()*1.0/RAND_MAX;
+ q = d->x_min + xi * d->N;
+ q = q / (floor(q) + 1);
+ q = pow(q, d->gamma);
+ if (val <= q){
+ return k;
+ }
+ }
+}
+
+
+
+
+int grow_multi_net_delay(ij_net_t *G1, ij_net_t *G2, int N, int m, int m0, coupling_t *matr){
+
+ int i, j;
+ int d;
+ int i2;
+
+ i2 = m0;
+
+ for(i=m0; i<N; i++){
+ compute_cumul(G1, G2, matr);
+ //dump_cumul(G1);
+ //dump_cumul(G2);
+ //n = m0 + i; /* This is the id of the new node */
+ G1->arrived[i] = i;
+ sample_edges(G1, G1->arrived[i], m, 0);
+ G1->degrees[G1->arrived[i]] += m;
+
+ for (j=0; j < G2->times[i].N; j++){
+ G2->arrived[i2] = G2->times[i].val[j];
+ //printf("%d\n", G2->arrived[i2]);
+ sample_edges(G2, G2->arrived[i2], m, m *j);
+ G2->degrees[G2->arrived[i2]] += m;
+ i2 += 1;
+ }
+ G1->N += 1;
+ G1->K += m;
+ G2->N += G2->times[i].N;
+ G2->K += m * G2->times[i].N;
+ }
+ return 0;
+}
+
+void init_structure(ij_net_t *G, int N){
+
+ int i;
+
+ G->i = malloc(G->size * sizeof(int));
+ G->j = malloc(G->size * sizeof(int));
+ G->degrees = malloc(N * sizeof(int));
+ G->arrived = malloc(N * sizeof(int));
+ G->cumul = malloc(N * sizeof(double));
+ G->times = malloc(N * sizeof(int_array_t));
+ for (i=0; i<N; i++){
+ G->times[i].size = 5;
+ G->times[i].N = 0;
+ G->times[i].val = malloc(5 * sizeof(int));
+ }
+ memset(G->degrees, 0, N*sizeof(int));
+ G->N = 0;
+
+}
+
+void add_node_to_time(ij_net_t *G, int node, int t){
+
+ if (G->times[t].size == G->times[t].N){
+ G->times[t].size += 5;
+ G->times[t].val = realloc(G->times[t].val, G->times[t].size);
+ }
+ G->times[t].val[G->times[t].N] = node;
+ G->times[t].N += 1;
+}
+
+
+void init_times_delta(ij_net_t *G, int N){
+
+ int i;
+
+ for (i=0; i<N; i ++){
+ add_node_to_time(G,i,i);
+ }
+
+}
+
+
+
+void init_times_delay(ij_net_t *G, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ //printf(" %d", t);
+ if (i+t < N){
+ add_node_to_time(G, i, t+i);
+ }
+ }
+ //printf("\n");
+}
+
+
+
+void dump_times(ij_net_t *G, int N){
+
+ int i, j;
+
+ for (i=0; i<N; i++){
+ for (j=0; j< G->times[i].N; j++){
+ printf("%d %d\n", i,G->times[i].val[j]);
+ }
+ //printf("\n");
+ }
+
+}
+
+
+int main(int argc, char *argv[]){
+
+ ij_net_t G1, G2;
+
+ int N, m, m0, k_max;
+ coupling_t matr;
+ double gamma;
+ distr_t delay_distr;
+
+ char str[256];
+
+ if (argc < 10){
+ printf("Usage: %s <N> <m> <m0> <outfile> <a> <b> <c> <d> <gamma>\n", argv[0]);
+ exit(1);
+ }
+
+ srand(time(NULL));
+
+ /* Diagonal coupling */
+ matr.a = atof(argv[5]);
+ matr.b = atof(argv[6]);
+ matr.c = atof(argv[7]);
+ matr.d = atof(argv[8]);
+ gamma = atof(argv[9]);
+
+ N = atoi(argv[1]);
+ m = atoi(argv[2]);
+ m0 = atoi(argv[3]);
+
+ G1.size = (N+m0) * m;
+ G2.size = (N+m0) * m;
+
+ init_structure(&G1, N);
+ init_structure(&G2, N);
+
+
+ G1.K = init_network(&G1, m0);
+ G2.K = init_network(&G2, m0);
+
+ delay_distr.N = N;
+ delay_distr.gamma = gamma;
+ delay_distr.x_min = 1;
+ delay_distr.distr = malloc(N * sizeof(double));
+
+ create_distr(delay_distr.distr, delay_distr.gamma, delay_distr.x_min, N);
+
+ init_times_delay(&G2, &delay_distr, m0, N);
+
+ dump_times(&G2, N);
+
+ fprintf(stderr, "Init finished!\n");
+
+ grow_multi_net_delay(&G1, &G2, N, m, m0, &matr);
+
+ //printf("### G1\n");
+ sprintf(str, "%s_layer1.txt", argv[4]);
+ dump_network_to_file(&G1, str);
+
+ //printf("### G2\n");
+ sprintf(str, "%s_layer2.txt", argv[4]);
+ dump_network_to_file(&G2, str);
+
+ /* dump_network_to_file(S, S_num, argv[4]); */
+ /* printf("Network dumped!\n"); */
+ /* k_max = degree_distr(S, S_num, &distr); */
+ /* printf("k_max is: %d\n", k_max); */
+ /* dump_distr_to_file(distr, k_max, argv[5]); */
+
+}
diff --git a/models/growth/nibilab_linear_delay_mix.c b/models/growth/nibilab_linear_delay_mix.c
new file mode 100644
index 0000000..e48d16b
--- /dev/null
+++ b/models/growth/nibilab_linear_delay_mix.c
@@ -0,0 +1,457 @@
+/**
+ *
+ * NiBiLaB model (Nicosia, Bianconi, Latora, Barthelemy) of multiplex
+ * linear preferential attachment
+ *
+ * At each time step, a new node arrives on one of the two layers
+ * (chosen uniformly at random), but its replica on the other layer
+ * arrives after a power-law distributed delay \tau
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+
+typedef struct{
+ int size;
+ int N;
+ int *val;
+} int_array_t;
+
+
+typedef struct{
+ int K;
+ int N;
+ int size;
+ int *i;
+ int *j;
+ int *degrees;
+ int *arrived;
+ double *cumul;
+ int_array_t *times;
+} ij_net_t;
+
+
+
+typedef struct{
+ double a;
+ double b;
+ double c;
+ double d;
+} coupling_t;
+
+
+typedef struct{
+ int N;
+ double x_min;
+ double *distr;
+ double gamma;
+} distr_t;
+
+
+int init_network(ij_net_t *G, int m0){
+
+ int n;
+
+ for(n=0; n<m0; n++){
+ G->i[n] = n;
+ G->j[n] = (n+1) % m0;
+ G->degrees[n] = 2;
+ G->arrived[n] = n;
+ }
+ G->K = m0;
+ G->N = m0;
+ return n;
+}
+
+
+
+void dump_network(ij_net_t *G){
+
+ int i;
+ for(i=0; i< G->K; i ++){
+ printf("%d %d\n",G->i[i], G->j[i]);
+ }
+}
+
+
+void dump_network_to_file(ij_net_t *G, char *fname){
+
+
+ FILE *f;
+ int i;
+
+
+ f = fopen(fname, "w+");
+
+ for(i=0; i< G->K; i ++){
+ fprintf(f, "%d %d\n",G->i[i], G->j[i]);
+ }
+ fclose(f);
+}
+
+
+
+double f1(double v1, double v2, coupling_t *matr){
+ return v1 * matr->a + v2 * matr->b;
+}
+
+double f2(double v1, double v2, coupling_t *matr){
+ return v1 * matr->c + v2 * matr->d;
+}
+
+
+void compute_cumul(ij_net_t *G1, ij_net_t *G2, coupling_t *matr){
+
+ int i;
+ double val;
+
+ G1->cumul[0] = f1(G1->degrees[ G1->arrived[0] ], G2->degrees[ G1->arrived[0] ], matr) ;
+ G2->cumul[0] = f2(G1->degrees[ G2->arrived[0] ], G2->degrees[ G2->arrived[0] ], matr) ;
+
+ for(i=1; i<G1->N; i ++){
+ val = f1(G1->degrees[ G1->arrived[i] ], G2->degrees[ G1->arrived[i] ],matr) ;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ val = f2(G1->degrees[ G2->arrived[i] ], G2->degrees[ G2->arrived[i] ], matr) ;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+void dump_cumul(ij_net_t *G){
+
+ int i;
+
+ for (i=0; i<G->N; i ++){
+ printf("%d %2.6f\n", G->times[i], G->cumul[i]);
+ }
+
+}
+
+
+
+int select_neigh_cumul(ij_net_t *G){
+
+ double r;
+ int j;
+
+ r = (rand() * 1.0 / RAND_MAX) * G->cumul[ G->N-1 ];
+ //printf(" r: %f ", r);
+
+ j = 0;
+ /* Change with a binary search ???!?!?!?! */
+ while (r > G->cumul[j]){
+ j ++;
+ }
+ return G->arrived[j];
+}
+
+int already_neighbour(ij_net_t *G, int j, int d, int offset){
+
+ int i;
+
+ for(i=G-> K + offset; i< G->K + offset + j; i ++){
+ if (G->j[i] == d)
+ return 1;
+ }
+ return 0;
+}
+
+
+void sample_edges(ij_net_t *G, int n, int m, int offset){
+
+ int j;
+ int d;
+
+ //printf("----");
+ for(j=0; j<m; j++){
+ G->i[G->K + offset + j] = n;
+ d = select_neigh_cumul(G);
+ while(already_neighbour(G, j, d, offset)){
+ d = select_neigh_cumul(G);
+ }
+ //printf(" link %d--%d\n", n, d);
+ G->j[G->K + offset + j] = d;
+ G->degrees[d] += 1;
+ }
+ //printf("\n");
+}
+
+void create_distr(double *v, double gamma, int x_min, int x_max){
+ int n, i;
+ double sum, new_sum;
+
+ n = x_max-x_min + 1;
+
+ sum = 0;
+ for(i=0; i<n; i++){
+ v[i] = pow((x_min + i), gamma);
+ sum += v[i];
+ }
+ new_sum = 0;
+ /* Now we normalize the array*/
+ for(i=0; i<n; i++){
+ v[i]/= sum;
+ new_sum += v[i];
+ }
+ //printf("New_sum: %lf\n", new_sum);
+ /* Now we compute the cumulative distribution*/
+ for(i=1; i<n; i++){
+ v[i] += v[i-1];
+ }
+}
+
+
+inline int find_degree(double *v, int dim, double xi){
+ int i;
+
+ i=0;
+ while(xi > v[i])
+ i++;
+ return i;
+
+}
+
+
+int sample_power_law(distr_t *d){
+
+ double xi, q, val;
+ int k;
+
+ while(1){
+ xi = rand()* 1.0 / RAND_MAX;
+ k = find_degree(d->distr, d->N, xi);
+ val = rand()*1.0/RAND_MAX;
+ q = d->x_min + xi * d->N;
+ q = q / (floor(q) + 1);
+ q = pow(q, d->gamma);
+ if (val <= q){
+ return k;
+ }
+ }
+}
+
+
+
+
+int grow_multi_net_delay_mixed(ij_net_t *G1, ij_net_t *G2, int N, int m, int m0, coupling_t *matr){
+
+ int i, j;
+ int d;
+ int i1, i2;
+
+ i1 = m0;
+ i2 = m0;
+
+ for(i=m0; i<N; i++){
+ compute_cumul(G1, G2, matr);
+ //dump_cumul(G1);
+ //dump_cumul(G2);
+ //n = m0 + i; /* This is the id of the new node */
+ for (j=0; j < G1->times[i].N; j++){
+ G1->arrived[i1] = G1->times[i].val[j];
+ sample_edges(G1, G1->arrived[i1], m, m * j);
+ G1->degrees[G1->arrived[i1]] += m;
+ i1 += 1;
+ }
+ for (j=0; j < G2->times[i].N; j++){
+ G2->arrived[i2] = G2->times[i].val[j];
+ //printf("%d\n", G2->arrived[i2]);
+ sample_edges(G2, G2->arrived[i2], m, m *j);
+ G2->degrees[G2->arrived[i2]] += m;
+ i2 += 1;
+ }
+ G1->N += G1->times[i].N;
+ G1->K += m * G1->times[i].N;
+ G2->N += G2->times[i].N;
+ G2->K += m * G2->times[i].N;
+ }
+ return 0;
+}
+
+void init_structure(ij_net_t *G, int N){
+
+ int i;
+
+ G->i = malloc(G->size * sizeof(int));
+ G->j = malloc(G->size * sizeof(int));
+ G->degrees = malloc(N * sizeof(int));
+ G->arrived = malloc(N * sizeof(int));
+ G->cumul = malloc(N * sizeof(double));
+ G->times = malloc(N * sizeof(int_array_t));
+ for (i=0; i<N; i++){
+ G->times[i].size = 5;
+ G->times[i].N = 0;
+ G->times[i].val = malloc(5 * sizeof(int));
+ }
+ memset(G->degrees, 0, N*sizeof(int));
+ G->N = 0;
+}
+
+void add_node_to_time(ij_net_t *G, int node, int t){
+
+ if (G->times[t].size == G->times[t].N){
+ G->times[t].size += 5;
+ G->times[t].val = realloc(G->times[t].val, G->times[t].size);
+ }
+ G->times[t].val[G->times[t].N] = node;
+ G->times[t].N += 1;
+}
+
+
+void init_times_delta(ij_net_t *G, int N){
+
+ int i;
+
+ for (i=0; i<N; i ++){
+ add_node_to_time(G,i,i);
+ }
+
+}
+
+
+
+void init_times_delay(ij_net_t *G, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ //printf(" %d", t);
+ if (i+t < N){
+ add_node_to_time(G, i, t+i);
+ }
+ }
+ //printf("\n");
+}
+
+/**
+ * Set the arrival time of nodes on each layer. A node $i$ selects its
+ * master layer uniformly at random, and then arrives in the other
+ * layer after a power-lay distributed delay
+ *
+ */
+
+void init_times_delay_mixed(ij_net_t *G1, ij_net_t *G2, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ double val;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ //printf(" %d", t);
+ val = rand() * 1.0 / RAND_MAX;
+ if (val <=0.5){ /* select layer 1 as the master layer*/
+ add_node_to_time(G1, i, i);
+ if (i+t < N){
+ add_node_to_time(G2, i, t+i);
+ }
+ }
+ else{/* select layer 2 as the master layer */
+ add_node_to_time(G2, i, i);
+ if (i+t < N){
+ add_node_to_time(G1, i, t+i);
+ }
+ }
+ }
+ //printf("\n");
+}
+
+
+
+
+void dump_times(ij_net_t *G, int N){
+
+ int i, j;
+
+ for (i=0; i<N; i++){
+ for (j=0; j< G->times[i].N; j++){
+ printf("%d %d\n", i,G->times[i].val[j]);
+ }
+ //printf("\n");
+ }
+
+}
+
+
+int main(int argc, char *argv[]){
+
+ ij_net_t G1, G2;
+
+ int N, m, m0, k_max;
+ coupling_t matr;
+ double gamma;
+ distr_t delay_distr;
+
+ char str[256];
+
+ if (argc < 10){
+ printf("Usage: %s <N> <m> <m0> <outfile> <a> <b> <c> <d> <gamma>\n", argv[0]);
+ exit(1);
+ }
+
+ srand(time(NULL));
+
+ /* Diagonal coupling */
+ matr.a = atof(argv[5]);
+ matr.b = atof(argv[6]);
+ matr.c = atof(argv[7]);
+ matr.d = atof(argv[8]);
+ gamma = atof(argv[9]);
+
+ N = atoi(argv[1]);
+ m = atoi(argv[2]);
+ m0 = atoi(argv[3]);
+
+ G1.size = (N+m0) * m;
+ G2.size = (N+m0) * m;
+
+ init_structure(&G1, N);
+ init_structure(&G2, N);
+
+
+ G1.K = init_network(&G1, m0);
+ G2.K = init_network(&G2, m0);
+
+ delay_distr.N = N;
+ delay_distr.gamma = gamma;
+ delay_distr.x_min = 1;
+ delay_distr.distr = malloc(N * sizeof(double));
+
+ create_distr(delay_distr.distr, delay_distr.gamma, delay_distr.x_min, N);
+
+ //init_times_delay(&G2, &delay_distr, m0, N);
+ init_times_delay_mixed(&G1, &G2, &delay_distr, m0, N);
+
+ printf("----- G1 times ----- \n");
+ dump_times(&G1, N);
+
+ printf("----- G2 times ----- \n");
+ dump_times(&G2, N);
+ //exit(2);
+
+ fprintf(stderr, "Init finished!\n");
+
+ grow_multi_net_delay_mixed(&G1, &G2, N, m, m0, &matr);
+
+ //printf("### G1\n");
+ sprintf(str, "%s_layer1.txt", argv[4]);
+ dump_network_to_file(&G1, str);
+
+ //printf("### G2\n");
+ sprintf(str, "%s_layer2.txt", argv[4]);
+ dump_network_to_file(&G2, str);
+
+ /* dump_network_to_file(S, S_num, argv[4]); */
+ /* printf("Network dumped!\n"); */
+ /* k_max = degree_distr(S, S_num, &distr); */
+ /* printf("k_max is: %d\n", k_max); */
+ /* dump_distr_to_file(distr, k_max, argv[5]); */
+
+}
diff --git a/models/growth/nibilab_linear_delta.c b/models/growth/nibilab_linear_delta.c
new file mode 100644
index 0000000..b1f7e90
--- /dev/null
+++ b/models/growth/nibilab_linear_delta.c
@@ -0,0 +1,403 @@
+/**
+ *
+ * NiBiLaB model (Nicosia, Bianconi, Latora, Barthelemy) of multiplex
+ * linear preferential attachment
+ *
+ * Node arrive at the same time on both layers.
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+
+typedef struct{
+ int size;
+ int N;
+ int *val;
+} int_array_t;
+
+
+typedef struct{
+ int K;
+ int N;
+ int size;
+ int *i;
+ int *j;
+ int *degrees;
+ int *arrived;
+ double *cumul;
+ int_array_t *times;
+} ij_net_t;
+
+
+
+typedef struct{
+ double a;
+ double b;
+ double c;
+ double d;
+} coupling_t;
+
+
+typedef struct{
+ int N;
+ double x_min;
+ double *distr;
+ double gamma;
+} distr_t;
+
+
+int init_network(ij_net_t *G, int m0){
+
+ int n;
+
+ for(n=0; n<m0; n++){
+ G->i[n] = n;
+ G->j[n] = n+1 % m0;
+ G->degrees[n] = 2;
+ G->arrived[n] = n;
+ }
+ G->K = m0;
+ G->N = m0;
+ return n;
+}
+
+
+
+void dump_network(ij_net_t *G){
+
+ int i;
+ for(i=0; i< G->K; i ++){
+ printf("%d %d\n",G->i[i], G->j[i]);
+ }
+}
+
+
+void dump_network_to_file(ij_net_t *G, char *fname){
+
+
+ FILE *f;
+ int i;
+
+
+ f = fopen(fname, "w+");
+
+ for(i=0; i< G->K; i ++){
+ fprintf(f, "%d %d\n",G->i[i], G->j[i]);
+ }
+ fclose(f);
+}
+
+
+
+double f1(double v1, double v2, coupling_t *matr){
+ return v1 * matr->a + v2 * matr->b;
+}
+
+double f2(double v1, double v2, coupling_t *matr){
+ return v1 * matr->c + v2 * matr->d;
+}
+
+
+void compute_cumul(ij_net_t *G1, ij_net_t *G2, coupling_t *matr){
+
+ int i;
+ double val;
+
+ G1->cumul[0] = f1(G1->degrees[ G1->arrived[0] ], G2->degrees[ G1->arrived[0] ], matr) ;
+ G2->cumul[0] = f2(G1->degrees[ G2->arrived[0] ], G2->degrees[ G2->arrived[0] ], matr) ;
+
+ for(i=1; i<G1->N; i ++){
+ val = f1(G1->degrees[ G1->arrived[i] ], G2->degrees[ G1->arrived[i] ],matr) ;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ val = f2(G1->degrees[ G2->arrived[i] ], G2->degrees[ G2->arrived[i] ], matr) ;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+void dump_cumul(ij_net_t *G){
+
+ int i;
+
+ for (i=0; i<G->N; i ++){
+ printf("%d %2.6f\n", G->times[i], G->cumul[i]);
+ }
+
+}
+
+
+
+int select_neigh_cumul(ij_net_t *G){
+
+ double r;
+ int j;
+
+ r = (rand() * 1.0 / RAND_MAX) * G->cumul[ G->N-1 ];
+ //printf(" r: %f ", r);
+
+ j = 0;
+ /* Change with a binary search ???!?!?!?! */
+ while (r > G->cumul[j]){
+ j ++;
+ }
+ return G->arrived[j];
+}
+
+int already_neighbour(ij_net_t *G, int j, int d, int offset){
+
+ int i;
+
+ for(i=G-> K + offset; i< G->K + offset + j; i ++){
+ if (G->j[i] == d)
+ return 1;
+ }
+ return 0;
+}
+
+
+void sample_edges(ij_net_t *G, int n, int m, int offset){
+
+ int j;
+ int d;
+
+ //printf("----");
+ for(j=0; j<m; j++){
+ G->i[G->K + offset + j] = n;
+ d = select_neigh_cumul(G);
+ while(already_neighbour(G, j, d, offset)){
+ d = select_neigh_cumul(G);
+ }
+ //printf(" link %d--%d\n", n, d);
+ G->j[G->K + offset + j] = d;
+ G->degrees[d] += 1;
+ }
+ //printf("\n");
+}
+
+void create_distr(double *v, double gamma, int x_min, int x_max){
+ int n, i;
+ double sum, new_sum;
+
+ n = x_max-x_min + 1;
+
+ sum = 0;
+ for(i=0; i<n; i++){
+ v[i] = pow((x_min + i), gamma);
+ sum += v[i];
+ }
+ new_sum = 0;
+ /* Now we normalize the array*/
+ for(i=0; i<n; i++){
+ v[i]/= sum;
+ new_sum += v[i];
+ }
+ //printf("New_sum: %lf\n", new_sum);
+ /* Now we compute the cumulative distribution*/
+ for(i=1; i<n; i++){
+ v[i] += v[i-1];
+ }
+}
+
+
+inline int find_degree(double *v, int dim, double xi){
+ int i;
+
+ i=0;
+ while(xi > v[i])
+ i++;
+ return i;
+
+}
+
+
+int sample_power_law(distr_t *d){
+
+ double xi, q, val;
+ int k;
+
+ while(1){
+ xi = rand()* 1.0 / RAND_MAX;
+ k = find_degree(d->distr, d->N, xi);
+ val = rand()*1.0/RAND_MAX;
+ q = d->x_min + xi * d->N;
+ q = q / (floor(q) + 1);
+ q = pow(q, d->gamma);
+ if (val <= q){
+ return k;
+ }
+ }
+}
+
+
+
+
+int grow_multi_net_delta(ij_net_t *G1, ij_net_t *G2, int N, int m, int m0, coupling_t *matr){
+
+ int i, j;
+ int d;
+ int i2;
+
+
+ for(i=m0; i<N; i++){
+ compute_cumul(G1, G2, matr);
+ //dump_cumul(G1);
+ //dump_cumul(G2);
+ //n = m0 + i; /* This is the id of the new node */
+ G1->arrived[i] = i;
+ sample_edges(G1, G1->arrived[i], m, 0);
+ G1->degrees[G1->arrived[i]] += m;
+
+ G2->arrived[i] = i;
+ //printf("%d\n", G2->arrived[i2]);
+ sample_edges(G2, G2->arrived[i], m, 0);
+ G2->degrees[G2->arrived[i]] += m;
+ G1->N += 1;
+ G1->K += m;
+ G2->N += 1;
+ G2->K += m;
+ }
+ return 0;
+}
+
+void init_structure(ij_net_t *G, int N){
+
+ int i;
+
+ G->i = malloc(G->size * sizeof(int));
+ G->j = malloc(G->size * sizeof(int));
+ G->degrees = malloc(N * sizeof(int));
+ G->arrived = malloc(N * sizeof(int));
+ G->cumul = malloc(N * sizeof(double));
+ G->times = malloc(N * sizeof(int_array_t));
+ for (i=0; i<N; i++){
+ G->times[i].size = 5;
+ G->times[i].N = 0;
+ G->times[i].val = malloc(5 * sizeof(int));
+ }
+ memset(G->degrees, 0, N*sizeof(int));
+ G->N = 0;
+
+}
+
+void add_node_to_time(ij_net_t *G, int node, int t){
+
+ if (G->times[t].size == G->times[t].N){
+ G->times[t].size += 5;
+ G->times[t].val = realloc(G->times[t].val, G->times[t].size);
+ }
+ G->times[t].val[G->times[t].N] = node;
+ G->times[t].N += 1;
+}
+
+
+void init_times_delta(ij_net_t *G, int N){
+
+ int i;
+
+ for (i=0; i<N; i ++){
+ add_node_to_time(G,i,i);
+ }
+
+}
+
+
+
+void init_times_delay(ij_net_t *G, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ printf(" %d", t);
+ if (i+t < N){
+ add_node_to_time(G, i, t+i);
+ }
+ }
+ printf("\n");
+}
+
+
+
+void dump_times(ij_net_t *G, int N){
+
+ int i, j;
+
+ for (i=0; i<N; i++){
+ printf("%d: ", i);
+ for (j=0; j< G->times[i].N; j++){
+ printf(" %d", G->times[i].val[j]);
+ }
+ printf("\n");
+ }
+
+}
+
+
+int main(int argc, char *argv[]){
+
+ ij_net_t G1, G2;
+
+ int N, m, m0, k_max;
+ coupling_t matr;
+ double gamma;
+ distr_t delay_distr;
+
+ char str[256];
+
+ if (argc < 9){
+ printf("Usage: %s <N> <m> <m0> <outfile> <a> <b> <c> <d>\n", argv[0]);
+ exit(1);
+ }
+
+ srand(time(NULL));
+
+ /* Diagonal coupling */
+ matr.a = atof(argv[5]);
+ matr.b = atof(argv[6]);
+ matr.c = atof(argv[7]);
+ matr.d = atof(argv[8]);
+
+ N = atoi(argv[1]);
+ m = atoi(argv[2]);
+ m0 = atoi(argv[3]);
+
+ G1.size = (N+m0) * m;
+ G2.size = (N+m0) * m;
+
+ init_structure(&G1, N);
+ init_structure(&G2, N);
+
+
+ G1.K = init_network(&G1, m0);
+ G2.K = init_network(&G2, m0);
+
+
+ init_times_delta(&G2, N);
+
+ //dump_times(&G2, N);
+
+ fprintf(stderr, "Init finished!\n");
+
+ grow_multi_net_delta(&G1, &G2, N, m, m0, &matr);
+
+ //printf("### G1\n");
+ sprintf(str, "%s_layer1.txt", argv[4]);
+ dump_network_to_file(&G1, str);
+
+ //printf("### G2\n");
+ sprintf(str, "%s_layer2.txt", argv[4]);
+ dump_network_to_file(&G2, str);
+
+ /* dump_network_to_file(S, S_num, argv[4]); */
+ /* printf("Network dumped!\n"); */
+ /* k_max = degree_distr(S, S_num, &distr); */
+ /* printf("k_max is: %d\n", k_max); */
+ /* dump_distr_to_file(distr, k_max, argv[5]); */
+
+}
diff --git a/models/growth/nibilab_linear_random_times.c b/models/growth/nibilab_linear_random_times.c
new file mode 100644
index 0000000..48929aa
--- /dev/null
+++ b/models/growth/nibilab_linear_random_times.c
@@ -0,0 +1,417 @@
+/**
+ *
+ * NiBiLaB model (Nicosia, Bianconi, Latora, Barthelemy) of multiplex
+ * linear preferential attachment
+ *
+ * Nodes arrive sequentially on the first layer, but their replicas on
+ * the second layer arrive at uniformly distributed times.
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+
+typedef struct{
+ int size;
+ int N;
+ int *val;
+} int_array_t;
+
+
+typedef struct{
+ int K;
+ int N;
+ int size;
+ int *i;
+ int *j;
+ int *degrees;
+ int *arrived;
+ double *cumul;
+ int_array_t *times;
+} ij_net_t;
+
+
+
+typedef struct{
+ double a;
+ double b;
+ double c;
+ double d;
+} coupling_t;
+
+
+typedef struct{
+ int N;
+ double x_min;
+ double *distr;
+ double gamma;
+} distr_t;
+
+
+int init_network(ij_net_t *G, int m0){
+
+ int n;
+
+ for(n=0; n<m0; n++){
+ G->i[n] = n;
+ G->j[n] = n+1 % m0;
+ G->degrees[n] = 2;
+ G->arrived[n] = n;
+ }
+ G->K = m0;
+ G->N = m0;
+ return n;
+}
+
+
+
+void dump_network(ij_net_t *G){
+
+ int i;
+ for(i=0; i< G->K; i ++){
+ printf("%d %d\n",G->i[i], G->j[i]);
+ }
+}
+
+
+void dump_network_to_file(ij_net_t *G, char *fname){
+
+
+ FILE *f;
+ int i;
+
+
+ f = fopen(fname, "w+");
+
+ for(i=0; i< G->K; i ++){
+ fprintf(f, "%d %d\n",G->i[i], G->j[i]);
+ }
+ fclose(f);
+}
+
+
+
+double f1(double v1, double v2, coupling_t *matr){
+ return v1 * matr->a + v2 * matr->b;
+}
+
+double f2(double v1, double v2, coupling_t *matr){
+ return v1 * matr->c + v2 * matr->d;
+}
+
+
+void compute_cumul(ij_net_t *G1, ij_net_t *G2, coupling_t *matr){
+
+ int i;
+ double val;
+
+ G1->cumul[0] = f1(G1->degrees[ G1->arrived[0] ], G2->degrees[ G1->arrived[0] ], matr) ;
+ G2->cumul[0] = f2(G1->degrees[ G2->arrived[0] ], G2->degrees[ G2->arrived[0] ], matr) ;
+
+ for(i=1; i<G1->N; i ++){
+ val = f1(G1->degrees[ G1->arrived[i] ], G2->degrees[ G1->arrived[i] ],matr) ;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ val = f2(G1->degrees[ G2->arrived[i] ], G2->degrees[ G2->arrived[i] ], matr) ;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+void dump_cumul(ij_net_t *G){
+
+ int i;
+
+ for (i=0; i<G->N; i ++){
+ printf("%d %2.6f\n", G->times[i], G->cumul[i]);
+ }
+
+}
+
+
+
+int select_neigh_cumul(ij_net_t *G){
+
+ double r;
+ int j;
+
+ r = (rand() * 1.0 / RAND_MAX) * G->cumul[ G->N-1 ];
+ //printf(" r: %f ", r);
+
+ j = 0;
+ /* Change with a binary search ???!?!?!?! */
+ while (r > G->cumul[j]){
+ j ++;
+ }
+ return G->arrived[j];
+}
+
+int already_neighbour(ij_net_t *G, int j, int d, int offset){
+
+ int i;
+
+ for(i=G-> K + offset; i< G->K + offset + j; i ++){
+ if (G->j[i] == d)
+ return 1;
+ }
+ return 0;
+}
+
+
+void sample_edges(ij_net_t *G, int n, int m, int offset){
+
+ int j;
+ int d;
+
+ //printf("----");
+ for(j=0; j<m; j++){
+ G->i[G->K + offset + j] = n;
+ d = select_neigh_cumul(G);
+ while(already_neighbour(G, j, d, offset)){
+ d = select_neigh_cumul(G);
+ }
+ //printf(" link %d--%d\n", n, d);
+ G->j[G->K + offset + j] = d;
+ G->degrees[d] += 1;
+ }
+ //printf("\n");
+}
+
+void create_distr(double *v, double gamma, int x_min, int x_max){
+ int n, i;
+ double sum, new_sum;
+
+ n = x_max-x_min + 1;
+
+ sum = 0;
+ for(i=0; i<n; i++){
+ v[i] = pow((x_min + i), gamma);
+ sum += v[i];
+ }
+ new_sum = 0;
+ /* Now we normalize the array*/
+ for(i=0; i<n; i++){
+ v[i]/= sum;
+ new_sum += v[i];
+ }
+ //printf("New_sum: %lf\n", new_sum);
+ /* Now we compute the cumulative distribution*/
+ for(i=1; i<n; i++){
+ v[i] += v[i-1];
+ }
+}
+
+
+inline int find_degree(double *v, int dim, double xi){
+ int i;
+
+ i=0;
+ while(xi > v[i])
+ i++;
+ return i;
+
+}
+
+
+int sample_power_law(distr_t *d){
+
+ double xi, q, val;
+ int k;
+
+ while(1){
+ xi = rand()* 1.0 / RAND_MAX;
+ k = find_degree(d->distr, d->N, xi);
+ val = rand()*1.0/RAND_MAX;
+ q = d->x_min + xi * d->N;
+ q = q / (floor(q) + 1);
+ q = pow(q, d->gamma);
+ if (val <= q){
+ return k;
+ }
+ }
+}
+
+
+
+
+int grow_multi_net_delay(ij_net_t *G1, ij_net_t *G2, int N, int m, int m0, coupling_t *matr){
+
+ int i, j;
+ int d;
+ int i2;
+
+ i2 = m0;
+
+ for(i=m0; i<N; i++){
+ compute_cumul(G1, G2, matr);
+ //dump_cumul(G1);
+ //dump_cumul(G2);
+ //n = m0 + i; /* This is the id of the new node */
+ G1->arrived[i] = i;
+ sample_edges(G1, G1->arrived[i], m, 0);
+ G1->degrees[G1->arrived[i]] += m;
+
+ for (j=0; j < G2->times[i].N; j++){
+ G2->arrived[i2] = G2->times[i].val[j];
+ //printf("%d\n", G2->arrived[i2]);
+ sample_edges(G2, G2->arrived[i2], m, m *j);
+ G2->degrees[G2->arrived[i2]] += m;
+ i2 += 1;
+ }
+ G1->N += 1;
+ G1->K += m;
+ G2->N += G2->times[i].N;
+ G2->K += m * G2->times[i].N;
+ }
+ return 0;
+}
+
+void init_structure(ij_net_t *G, int N){
+
+ int i;
+
+ G->i = malloc(G->size * sizeof(int));
+ G->j = malloc(G->size * sizeof(int));
+ G->degrees = malloc(N * sizeof(int));
+ G->arrived = malloc(N * sizeof(int));
+ G->cumul = malloc(N * sizeof(double));
+ G->times = malloc(N * sizeof(int_array_t));
+ for (i=0; i<N; i++){
+ G->times[i].size = 5;
+ G->times[i].N = 0;
+ G->times[i].val = malloc(5 * sizeof(int));
+ }
+ memset(G->degrees, 0, N*sizeof(int));
+ G->N = 0;
+
+}
+
+void add_node_to_time(ij_net_t *G, int node, int t){
+
+ if (G->times[t].size == G->times[t].N){
+ G->times[t].size += 5;
+ G->times[t].val = realloc(G->times[t].val, G->times[t].size);
+ }
+ G->times[t].val[G->times[t].N] = node;
+ G->times[t].N += 1;
+}
+
+
+void init_times_delta(ij_net_t *G, int N){
+
+ int i;
+
+ for (i=0; i<N; i ++){
+ add_node_to_time(G,i,i);
+ }
+
+}
+
+
+
+void init_times_delay(ij_net_t *G, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ printf(" %d", t);
+ if (i+t < N){
+ add_node_to_time(G, i, t+i);
+ }
+ }
+ printf("\n");
+}
+
+void init_times_random(ij_net_t *G, int N){
+
+ int i,j;
+
+ for (i=0; i<N; i ++){
+ j = rand() % N;
+ add_node_to_time(G, i, j);
+ }
+}
+
+
+void dump_times(ij_net_t *G, int N){
+
+ int i, j;
+
+ for (i=0; i<N; i++){
+ printf("%d: ", i);
+ for (j=0; j< G->times[i].N; j++){
+ printf(" %d", G->times[i].val[j]);
+ }
+ printf("\n");
+ }
+
+}
+
+
+int main(int argc, char *argv[]){
+
+ ij_net_t G1, G2;
+
+ int N, m, m0, k_max;
+ coupling_t matr;
+ double gamma;
+ distr_t delay_distr;
+
+ char str[256];
+
+ if (argc < 9){
+ printf("Usage: %s <N> <m> <m0> <outfile> <a> <b> <c> <d>\n", argv[0]);
+ exit(1);
+ }
+
+ srand(time(NULL));
+
+ /* Diagonal coupling */
+ matr.a = atof(argv[5]);
+ matr.b = atof(argv[6]);
+ matr.c = atof(argv[7]);
+ matr.d = atof(argv[8]);
+
+ N = atoi(argv[1]);
+ m = atoi(argv[2]);
+ m0 = atoi(argv[3]);
+
+ G1.size = (N+m0) * m;
+ G2.size = (N+m0) * m;
+
+ init_structure(&G1, N);
+ init_structure(&G2, N);
+
+
+ G1.K = init_network(&G1, m0);
+ G2.K = init_network(&G2, m0);
+
+
+ init_times_random(&G2, N);
+
+ //dump_times(&G2, N);
+
+ fprintf(stderr, "Init finished!\n");
+
+ grow_multi_net_delay(&G1, &G2, N, m, m0, &matr);
+
+ //printf("### G1\n");
+ sprintf(str, "%s_layer1.txt", argv[4]);
+ dump_network_to_file(&G1, str);
+
+ //printf("### G2\n");
+ sprintf(str, "%s_layer2.txt", argv[4]);
+ dump_network_to_file(&G2, str);
+
+ /* dump_network_to_file(S, S_num, argv[4]); */
+ /* printf("Network dumped!\n"); */
+ /* k_max = degree_distr(S, S_num, &distr); */
+ /* printf("k_max is: %d\n", k_max); */
+ /* dump_distr_to_file(distr, k_max, argv[5]); */
+
+}
diff --git a/models/growth/nibilab_nonlinear.c b/models/growth/nibilab_nonlinear.c
new file mode 100644
index 0000000..3ad9ce8
--- /dev/null
+++ b/models/growth/nibilab_nonlinear.c
@@ -0,0 +1,493 @@
+/**
+ *
+ * NiBiLaB model (Nicosia, Bianconi, Latora, Barthelemy) of multiplex
+ * non-linear preferential attachment.
+ *
+ * The probability for a newly arrived node $i$ to create a link to an
+ * existing node $j$ is
+ *
+ * p_{i->j} \propto k\lay{1}^{\alpha}/k\lay{2}^{\beta}
+ *
+ * where alpha and beta are two parameters.
+ *
+ * Nodes arrive at the same time on both layers.
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+
+
+typedef struct{
+ int size;
+ int N;
+ int *val;
+} int_array_t;
+
+
+typedef struct{
+ int K;
+ int N;
+ int size;
+ int *i;
+ int *j;
+ int *degrees;
+ int *arrived;
+ double *cumul;
+ int_array_t *times;
+ double *eta;
+} ij_net_t;
+
+
+
+typedef struct{
+ double a;
+ double b;
+ double c;
+ double d;
+} coupling_t;
+
+
+typedef struct{
+ int N;
+ double x_min;
+ double *distr;
+ double gamma;
+} distr_t;
+
+
+int init_network(ij_net_t *G, int m0){
+
+ int n;
+
+ for(n=0; n<m0; n++){
+ G->i[n] = n;
+ G->j[n] = n+1 % m0;
+ G->degrees[n] = 2;
+ G->arrived[n] = n;
+ }
+ G->K = m0;
+ G->N = m0;
+ return n;
+}
+
+
+
+void dump_network(ij_net_t *G){
+
+ int i;
+ for(i=0; i< G->K; i ++){
+ printf("%d %d\n",G->i[i], G->j[i]);
+ }
+}
+
+
+void dump_network_to_file(ij_net_t *G, char *fname){
+
+
+ FILE *f;
+ int i;
+
+
+ f = fopen(fname, "w+");
+
+ for(i=0; i< G->K; i ++){
+ fprintf(f, "%d %d\n",G->i[i], G->j[i]);
+ }
+ fclose(f);
+}
+
+
+
+double f1(double v1, double v2, coupling_t *matr){
+ return v1 * matr->a + v2 * matr->b;
+}
+
+double f2(double v1, double v2, coupling_t *matr){
+ return v1 * matr->c + v2 * matr->d;
+}
+
+
+void compute_cumul_nlin(ij_net_t *G1, ij_net_t *G2, double alpha){
+
+ int i;
+ double val;
+
+ /* The bias at layer 1 is k1/k2^2 */
+ if (G1->degrees[G1->arrived[0]] > 0 && G2->degrees[ G1->arrived[0] ])
+ G1->cumul[0] = (G1->degrees[G1->arrived[0]]) * 1.0 / pow(G2->degrees[ G1->arrived[0] ],alpha) ;
+ else
+ G1->cumul[0] = 0.0;
+
+ if (G2->degrees[G2->arrived[0]] > 0 && G1->degrees[ G2->arrived[0] ] > 0)
+ G2->cumul[0] = ( G2->degrees[G2->arrived[0]]) * 1.0/ pow(G1->degrees[ G2->arrived[0] ],alpha);
+ else
+ G2->cumul[0] = 0.0;
+
+ for(i=1; i<G1->N; i ++){
+ if (G1->degrees[G1->arrived[i]] > 0 && G2->degrees[ G1->arrived[i] ] > 0)
+ val = (G1->degrees[G1->arrived[i]]) * 1.0 / pow(G2->degrees[ G1->arrived[i] ],alpha) ;
+ else
+ val = 0;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ if (G2->degrees[G2->arrived[i]] > 0 && G1->degrees[ G2->arrived[i] ] )
+ val = (G2->degrees[G2->arrived[i]]) * 1.0/ pow(G1->degrees[ G2->arrived[i] ], alpha);
+ else
+ val = 0.0;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+
+void compute_cumul_nlin_2(ij_net_t *G1, ij_net_t *G2, double alpha, double beta){
+
+ int i;
+ double val;
+
+ /* The bias at layer 1 is k1/k2^2 */
+ if (G1->degrees[G1->arrived[0]] > 0 && G2->degrees[ G1->arrived[0] ])
+ G1->cumul[0] = pow(G1->degrees[G1->arrived[0]], alpha) / pow(G2->degrees[ G1->arrived[0] ],beta) ;
+ else
+ G1->cumul[0] = 0.0;
+
+ if (G2->degrees[G2->arrived[0]] > 0 && G1->degrees[ G2->arrived[0] ] > 0)
+ G2->cumul[0] = pow(G2->degrees[G2->arrived[0]], alpha)/pow(G1->degrees[ G2->arrived[0] ], beta);
+ else
+ G2->cumul[0] = 0.0;
+
+ for(i=1; i<G1->N; i ++){
+ if (G1->degrees[G1->arrived[i]] > 0 && G2->degrees[ G1->arrived[i] ] > 0)
+ val = pow(G1->degrees[G1->arrived[i]], alpha)/pow(G2->degrees[ G1->arrived[i] ],beta) ;
+ else
+ val = 0;
+ G1->cumul[i] = G1->cumul[i-1] + val;
+ }
+ for(i=1; i<G2->N; i ++){
+ if (G2->degrees[G2->arrived[i]] > 0 && G1->degrees[ G2->arrived[i] ] )
+ val = pow(G2->degrees[G2->arrived[i]],alpha)/pow(G1->degrees[ G2->arrived[i] ], beta);
+ else
+ val = 0.0;
+ G2->cumul[i] = G2->cumul[i-1] + val;
+ }
+}
+
+
+void dump_cumul(ij_net_t *G){
+
+ int i;
+
+ for (i=0; i<G->N; i ++){
+ printf("%d %2.6f\n", G->times[i], G->cumul[i]);
+ }
+
+}
+
+
+
+int select_neigh_cumul(ij_net_t *G){
+
+ double r;
+ int j;
+
+ r = (rand() * 1.0 / RAND_MAX) * G->cumul[ G->N-1 ];
+ //printf(" r: %f ", r);
+
+ j = 0;
+ /* Change with a binary search ???!?!?!?! */
+ while (r > G->cumul[j]){
+ j ++;
+ }
+ return G->arrived[j];
+}
+
+int already_neighbour(ij_net_t *G, int j, int d, int offset){
+
+ int i;
+
+ for(i=G-> K + offset; i< G->K + offset + j; i ++){
+ if (G->j[i] == d)
+ return 1;
+ }
+ return 0;
+}
+
+
+void sample_edges(ij_net_t *G, int n, int m, int offset){
+
+ int j;
+ int d;
+
+ //printf("----");
+ for(j=0; j<m; j++){
+ G->i[G->K + offset + j] = n;
+ d = select_neigh_cumul(G);
+ while(already_neighbour(G, j, d, offset)){
+ d = select_neigh_cumul(G);
+ }
+ //printf(" link %d--%d\n", n, d);
+ G->j[G->K + offset + j] = d;
+ G->degrees[d] += 1;
+ }
+ //printf("\n");
+}
+
+void create_distr(double *v, double gamma, int x_min, int x_max){
+ int n, i;
+ double sum, new_sum;
+
+ n = x_max-x_min + 1;
+
+ sum = 0;
+ for(i=0; i<n; i++){
+ v[i] = pow((x_min + i), gamma);
+ sum += v[i];
+ }
+ new_sum = 0;
+ /* Now we normalize the array*/
+ for(i=0; i<n; i++){
+ v[i]/= sum;
+ new_sum += v[i];
+ }
+ //printf("New_sum: %lf\n", new_sum);
+ /* Now we compute the cumulative distribution*/
+ for(i=1; i<n; i++){
+ v[i] += v[i-1];
+ }
+}
+
+
+inline int find_degree(double *v, int dim, double xi){
+ int i;
+
+ i=0;
+ while(xi > v[i])
+ i++;
+ return i;
+
+}
+
+
+int sample_power_law(distr_t *d){
+
+ double xi, q, val;
+ int k;
+
+ while(1){
+ xi = rand()* 1.0 / RAND_MAX;
+ k = find_degree(d->distr, d->N, xi);
+ val = rand()*1.0/RAND_MAX;
+ q = d->x_min + xi * d->N;
+ q = q / (floor(q) + 1);
+ q = pow(q, d->gamma);
+ if (val <= q){
+ return k;
+ }
+ }
+}
+
+
+
+
+int grow_multi_net_nlin_2(ij_net_t *G1, ij_net_t *G2, int N, int m,
+ int m0, double alpha, double beta){
+
+ int i, j;
+ int d;
+ int i2;
+
+
+ for(i=m0; i<N; i++){
+ compute_cumul_nlin_2(G1, G2, alpha, beta);
+ //dump_cumul(G1);
+ //dump_cumul(G2);
+ //n = m0 + i; /* This is the id of the new node */
+ G1->arrived[i] = i;
+ sample_edges(G1, G1->arrived[i], m, 0);
+ G1->degrees[G1->arrived[i]] += m;
+
+ G2->arrived[i] = i;
+ sample_edges(G2, G2->arrived[i], m, 0);
+ G2->degrees[G2->arrived[i]] += m;
+ G1->N += 1;
+ G1->K += m;
+ G2->N += 1;
+ G2->K += m ;
+ }
+ return 0;
+}
+
+void init_structure(ij_net_t *G, int N){
+
+ int i;
+
+ G->i = malloc(G->size * sizeof(int));
+ G->j = malloc(G->size * sizeof(int));
+ G->degrees = malloc(N * sizeof(int));
+ G->arrived = malloc(N * sizeof(int));
+ G->cumul = malloc(N * sizeof(double));
+ G->eta = malloc(N * sizeof(double));
+ G->times = malloc(N * sizeof(int_array_t));
+ for (i=0; i<N; i++){
+ G->times[i].size = 5;
+ G->times[i].N = 0;
+ G->times[i].val = malloc(5 * sizeof(int));
+ }
+ memset(G->degrees, 0, N*sizeof(int));
+ G->N = 0;
+
+}
+
+void add_node_to_time(ij_net_t *G, int node, int t){
+
+ if (G->times[t].size == G->times[t].N){
+ G->times[t].size += 5;
+ G->times[t].val = realloc(G->times[t].val, G->times[t].size);
+ }
+ G->times[t].val[G->times[t].N] = node;
+ G->times[t].N += 1;
+}
+
+
+void init_times_delta(ij_net_t *G, int N){
+
+ int i;
+
+ for (i=0; i<N; i ++){
+ add_node_to_time(G,i,i);
+ }
+
+}
+
+
+
+void init_times_delay(ij_net_t *G, distr_t *d, int m0, int N){
+ int i;
+ int t;
+
+ for (i=m0; i<N; i ++){
+ t = sample_power_law(d);/* So t is in the interval [1,N] */
+ printf(" %d", t);
+ if (i+t < N){
+ add_node_to_time(G, i, t+i);
+ }
+ }
+ printf("\n");
+}
+
+
+
+void dump_times(ij_net_t *G, int N){
+
+ int i, j;
+
+ for (i=0; i<N; i++){
+ printf("%d: ", i);
+ for (j=0; j< G->times[i].N; j++){
+ printf(" %d", G->times[i].val[j]);
+ }
+ printf("\n");
+ }
+
+}
+
+
+void init_eta_random(ij_net_t *G, int N){
+
+ int i;
+ double val;
+
+ for (i=0; i<N; i ++){
+ val = rand() * 1.0/RAND_MAX;
+ G->eta[i] = val;
+ }
+}
+
+void init_eta_ass(ij_net_t *G2, ij_net_t *G1, int N){
+ int i;
+
+ for (i=0; i<N; i ++){
+ G2->eta[i] = G1->eta[i];
+ }
+
+}
+
+void init_eta_dis(ij_net_t *G2, ij_net_t *G1, int N){
+ int i;
+
+ for (i=0; i<N; i ++){
+ G2->eta[i] = 1 - G1->eta[i];
+ }
+
+}
+
+
+int main(int argc, char *argv[]){
+
+ ij_net_t G1, G2;
+
+ int N, m, m0, k_max;
+ coupling_t matr;
+ double alpha, beta;
+ distr_t delay_distr;
+
+ char str[256];
+
+ if (argc < 6){
+ printf("Usage: %s <N> <m> <m0> <outfile> <alpha> <beta>\n", argv[0]);
+ exit(1);
+ }
+
+ srand(time(NULL));
+
+ /* Diagonal coupling */
+
+ N = atoi(argv[1]);
+ m = atoi(argv[2]);
+ m0 = atoi(argv[3]);
+ alpha=atof(argv[5]);
+ beta = atof(argv[6]);
+
+ G1.size = (N+m0) * m;
+ G2.size = (N+m0) * m;
+
+ init_structure(&G1, N);
+ init_structure(&G2, N);
+
+
+ G1.K = init_network(&G1, m0);
+ G2.K = init_network(&G2, m0);
+
+ //init_eta_random(&G1, N);
+ //init_eta_random(&G2, N);
+
+
+
+ //init_times_delta(&G1, N);
+ //init_times_delta(&G2, N);
+
+ //dump_times(&G2, N);
+
+ fprintf(stderr, "Init finished!\n");
+
+ grow_multi_net_nlin_2(&G1, &G2, N, m, m0, alpha, beta);
+
+ //printf("### G1\n");
+ sprintf(str, "%s_layer1.txt", argv[4]);
+ dump_network_to_file(&G1, str);
+
+ //printf("### G2\n");
+ sprintf(str, "%s_layer2.txt", argv[4]);
+ dump_network_to_file(&G2, str);
+
+ /* dump_network_to_file(S, S_num, argv[4]); */
+ /* printf("Network dumped!\n"); */
+ /* k_max = degree_distr(S, S_num, &distr); */
+ /* printf("k_max is: %d\n", k_max); */
+ /* dump_distr_to_file(distr, k_max, argv[5]); */
+
+}
diff --git a/models/growth/node_deg_over_time.py b/models/growth/node_deg_over_time.py
new file mode 100644
index 0000000..a71e86c
--- /dev/null
+++ b/models/growth/node_deg_over_time.py
@@ -0,0 +1,82 @@
+####
+##
+## Get an edgelist, a file pof arrival times and a node ID and return
+## the degree of that node as a function of time (we suppose that IDs
+## are sequential)
+##
+##
+
+import sys
+
+
+if len(sys.argv) < 4:
+ print "Usage: %s <netfile> <timefile> <nodeid1> [<nodeid2> <nodeid3>...]" % sys.argv[0]
+ sys.exit(1)
+
+node_id = int(sys.argv[3])
+
+lines = open(sys.argv[2], "r").readlines()
+
+arrival_time = {}
+
+#### WARNING!!!! THIS WORKS ONLY FOR M0=3
+arrival_time[0] = 0
+arrival_time[1] = 1
+arrival_time[2] = 2
+
+
+neigh_by_time = {}
+
+max_t = -1
+
+for l in lines:
+ if l[0] == "#":
+ continue
+ t,node = [int(x) for x in l.strip(" \n").split(" ")]
+ arrival_time[node] = t
+ if t > max_t :
+ max_t = t
+
+
+lines = open(sys.argv[1], "r").readlines()
+
+
+for l in lines:
+ if l[0] == "#":
+ continue
+ n1,n2 = [int(x) for x in l.strip(" \n").split(" ")]
+
+
+ if neigh_by_time.has_key(n1):
+ neigh_by_time[n1].append(arrival_time[n2])
+ else:
+ neigh_by_time[n1] = [arrival_time[n2]]
+
+ if neigh_by_time.has_key(n2):
+ neigh_by_time[n2].append(arrival_time[n1])
+ else:
+ neigh_by_time[n2] = [arrival_time[n1]]
+
+
+
+#print neigh_by_time[node_id]
+
+
+for node_id in sys.argv[3:]:
+ node_id = int(node_id)
+ neigh_by_time[node_id].sort()
+ last_time = neigh_by_time[node_id][0]
+#### changed here
+ k = 1
+ print "#### ", node_id
+ for t in neigh_by_time[node_id][1:]:
+ if t != last_time:
+ if last_time < arrival_time[node_id]:
+ print arrival_time[node_id], k
+ else:
+ print last_time, k
+ last_time = t
+ k += 1
+ print max_t, k-1
+ print
+ print
diff --git a/models/nullmodels/model_MDM.py b/models/nullmodels/model_MDM.py
new file mode 100644
index 0000000..5b0a969
--- /dev/null
+++ b/models/nullmodels/model_MDM.py
@@ -0,0 +1,66 @@
+####
+##
+##
+## This is the vertical participation model. For each node i, we use
+## exactly the same value of B_i as in the original network, but we
+## choose at random the layers in which node i will be active. This
+## breaks down intra-layer correlations.
+##
+## We get as input a file which reports, for each value of B_i, the
+## number of nodes in the original network which have that value, i the format:
+##
+## B_i N(B_i)
+##
+##
+##
+## The output is the obtained distribution of bit-strings.
+##
+##
+
+import sys
+import random
+
+
+def to_binary(l):
+ s = 0
+ e = 0
+ for v in l:
+ s += v * pow(2,e)
+ e +=1
+ return s
+
+
+if len(sys.argv) < 3:
+ print "Usage: %s <Bi_file> <M>" % sys.argv[0]
+ sys.exit(1)
+
+M = int(sys.argv[2])
+
+layers = range(M)
+
+distr = {}
+
+with open(sys.argv[1], "r") as f:
+ for l in f:
+ if l[0] == "#":
+ continue
+ val, num = [int(x) for x in l.strip(" \n").split(" ")]
+ for j in range(num):
+ node_layers = random.sample(layers, val)
+ node_bitstring = [0 for x in range(M)]
+ #print node_bitstring, node_layers
+ for i in node_layers:
+ #print i,
+ node_bitstring[i] = 1
+ #print node_bitstring
+
+ bs = to_binary(node_bitstring)
+ if bs in distr:
+ distr[bs] += 1
+ else:
+ distr[bs] = 1
+
+for k in distr:
+ bin_list = bin(k)
+ bin_num = sum([int(x) if x=='1' else 0 for x in bin_list[2:]])
+ sys.stderr.write("%d %0175s %d \n" % (bin_num, bin_list[2:], distr[k]))
diff --git a/models/nullmodels/model_MSM.py b/models/nullmodels/model_MSM.py
new file mode 100644
index 0000000..2c30df5
--- /dev/null
+++ b/models/nullmodels/model_MSM.py
@@ -0,0 +1,65 @@
+####
+##
+##
+## Create a synthetic multiplex network in which a node $i$ appears at
+## each layer $\alpha$ with a probability equal to $B_i$, which is the
+## fraction of layers in which node $i$ participate in the original
+## multiplex.
+##
+## Take a file of node binary participation indices, and sample a
+## multiplex compatible with that distribution
+##
+##
+## The input file is in the format:
+##
+## nodeID_i B_i
+##
+## The output file is a node-layer participation file, in the format
+##
+## node_id1 layer_id1
+## node_id2 layer_id2
+## .....
+##
+
+import sys
+
+import random
+
+if len(sys.argv) < 3:
+ print "Usage: %s <filein> <M>" % sys.argv[0]
+ sys.exit(1)
+
+M = int(sys.argv[2])
+
+bin_index = {}
+node_ids = []
+
+lines = open(sys.argv[1]).readlines()
+
+
+for l in lines:
+ if l[0] == "#":
+ continue
+ elems = [int(x) for x in l.strip(" \n").split(" ")]
+ bin_index[elems[0]] = 1.0 * elems[1]/M
+ node_ids.append(elems[0])
+
+N = len(node_ids)
+
+node_layers = {}
+
+for alpha in range (M):
+ for i in node_ids:
+ val = random.random()
+ if val < bin_index[i]:
+ if node_layers.has_key(i):
+ node_layers[i].append(alpha)
+ else:
+ node_layers[i] = [alpha]
+
+
+for i in node_ids:
+ if node_layers.has_key(i):
+ for j in range(len(node_layers[i])):
+ print i, node_layers[i][j]
+
diff --git a/models/nullmodels/model_hypergeometric.py b/models/nullmodels/model_hypergeometric.py
new file mode 100644
index 0000000..0a36237
--- /dev/null
+++ b/models/nullmodels/model_hypergeometric.py
@@ -0,0 +1,56 @@
+####
+##
+## This is the hypergeometric model. Each layer has the same number of
+## non-isolated nodes as the initial graph, but the nodes are
+## activated at random. The input is a file of number of non-isolated
+## nodes in each layer, and the total number of nodes in the multiplex.
+##
+## The output file is a node-layer participation file, in the format
+##
+## node_id1 layer_id1
+## node_id2 layer_id2
+## .....
+##
+
+import sys
+import random
+
+if len(sys.argv) < 3:
+ print "Usage: %s <layer_N_file> <N>" % sys.argv[0]
+ sys.exit(1)
+
+N = int(sys.argv[2])
+
+layer_degs = []
+node_layers = {}
+
+lines = open(sys.argv[1]).readlines()
+
+M = 0
+
+
+for l in lines:
+ if l[0] == "#":
+ continue
+ n = [int(x) for x in l.strip(" \n").split(" ")][0]
+ layer_degs.append(n)
+ M += 1
+
+for i in range(M):
+ num = layer_degs[i]
+ added = []
+ n = 0
+ while n < num:
+ j = random.choice(range(N))
+ if j not in added:
+ n += 1
+ added.append(j)
+ if node_layers.has_key(j):
+ node_layers[j].append(i)
+ else:
+ node_layers[j] = [i]
+
+for n in node_layers.keys():
+ for i in node_layers[n]:
+ print n,i
+
diff --git a/models/nullmodels/model_layer_growth.py b/models/nullmodels/model_layer_growth.py
new file mode 100644
index 0000000..46b4c0f
--- /dev/null
+++ b/models/nullmodels/model_layer_growth.py
@@ -0,0 +1,121 @@
+####
+##
+## layer-by-layer multiplex growth
+##
+## - We start from a multiplex with M_0 layers, with a certain number of
+## active nodes each
+##
+## - Each new layer arrives with a certain number N\lay{\alpha} of nodes
+## to be activated (this number is sampled from the observed distribution
+## of N\lay{\alpha}, like in the airlines multiplex)
+##
+## - Each node $i$ is activated with a probability proportional to the
+## number of existing layers in which it is already active, added to an
+## attractivity A :
+##
+## P_i(t) \propto A + B_i(t)
+##
+## - the hope is that A might tune the exponent of the resulting distribution
+## of B_i.....
+##
+##
+## This script takes as input a file which contains the degrees of the
+## layers, the total number of nodes in the multiplex, the initial
+## number M0 of layers in the multiplex and the attractiveness
+## parameter A. If "RND" is specified as a third parameter, then the
+## sequence of N\lay{\alpha} is shuffled
+##
+## Gives as output a list of node-layer participation
+##
+
+import sys
+import random
+
+if len(sys.argv) < 5:
+ print "Usage: %s <layers_N> <N> <M0> <A> [RND]" % sys.argv[0]
+ sys.exit(1)
+
+N = int(sys.argv[2])
+M0 = int(sys.argv[3])
+A = int(sys.argv[4])
+
+layer_degs = []
+
+
+if len(sys.argv) > 5 and sys.argv[5] == "RND":
+ randomize = True
+else:
+ randomize = False
+
+lines = open(sys.argv[1]).readlines()
+
+M = 0
+
+
+for l in lines:
+ if l[0] == "#":
+ continue
+ n = [int(x) for x in l.strip(" \n").split(" ")][0]
+ layer_degs.append(n)
+ M += 1
+
+
+if randomize:
+ random.shuffle(layer_degs)
+
+## the list node_Bi contains, at each time, the attachment
+## probabilities, i.e. it is a list which contains each node $i$
+## exactly $A + B_i$ times
+
+node_Bi = []
+
+#
+# initialize the distribution of attachment proibabilities, giving to
+# all nodes an attachment probability equal to A
+#
+
+for i in range(N):
+ for j in range(A):
+ node_Bi.append(i)
+
+layers = []
+
+
+for i in range(M0):
+ N_alpha = layer_degs.pop()
+ node_list = []
+ num_nodes = 0
+ while num_nodes < N_alpha:
+ val = random.choice(range(N))
+ if val not in node_list:
+ node_list.append(val)
+ num_nodes += 1
+ for n in node_list:
+ node_Bi.append(n)
+ layers.append(node_list)
+ i += 1
+
+
+#sys.exit(1)
+
+while i < M:
+ N_alpha = layer_degs.pop()
+ node_list = []
+ num_nodes = 0
+ while num_nodes < N_alpha:
+ val = random.choice(node_Bi)
+ if val not in node_list:
+ node_list.append(val)
+ num_nodes += 1
+ for n in node_list:
+ node_Bi.append(n)
+ layers.append(node_list)
+ i += 1
+
+#print layers
+
+for i in range(M):
+ node_list = layers[i]
+ for n in node_list:
+ print n, i
+
diff --git a/structure/activity/degs_to_activity_overlap.py b/structure/activity/degs_to_activity_overlap.py
new file mode 100644
index 0000000..c96959b
--- /dev/null
+++ b/structure/activity/degs_to_activity_overlap.py
@@ -0,0 +1,36 @@
+####
+##
+## Take a file which contains, on the n-th line, the degrees at each
+## layer of the n-th node, and print on output the corresponding node
+## multi-activity (i.e., the number of layers in which the node is
+## active) and the overlapping degree (i.e., the total number of edges
+## incident on the node)
+##
+##
+
+
+import sys
+
+def to_binary(l):
+ s = 0
+ e = 0
+ for v in l:
+ s += v * pow(2,e)
+ e +=1
+ return s
+
+if len(sys.argv) < 2:
+ print "Usage: %s <filein>" % sys.argv[0]
+ sys.exit(1)
+
+distr = {}
+
+with open(sys.argv[1]) as f:
+ for l in f:
+ elems = [int(x) for x in l.strip(" \n").split(" ")]
+ ov = sum(elems)
+ new_list = [1 if x>0 else 0 for x in elems]
+ multi_act = sum(new_list)
+ if multi_act and ov:
+ print ov, multi_act
+
diff --git a/structure/activity/degs_to_binary.py b/structure/activity/degs_to_binary.py
new file mode 100644
index 0000000..cb7aef5
--- /dev/null
+++ b/structure/activity/degs_to_binary.py
@@ -0,0 +1,45 @@
+####
+##
+## Take a file which contains, on the n-th line, the degrees at each
+## layer of the n-th node, and print on output the corresponding node
+## participation bit-strings, i.e. the string which contains "1" if on
+## that layer the node is connected, and zero otherwise
+##
+## on the stderr, we also dump the distribution of each bit-string
+##
+
+
+import sys
+
+def to_binary(l):
+ s = 0
+ e = 0
+ for v in l:
+ s += v * pow(2,e)
+ e +=1
+ return s
+
+if len(sys.argv) < 2:
+ print "Usage: %s <filein>" % sys.argv[0]
+ sys.exit(1)
+
+distr = {}
+
+with open(sys.argv[1]) as f:
+ for l in f:
+ elems = [int(x) for x in l.strip(" \n").split(" ")]
+ new_list = [1 if x>0 else 0 for x in elems]
+ val = to_binary(new_list)
+ if val in distr:
+ distr[val] += 1
+ else:
+ distr[val] = 1
+ for i in new_list:
+ print i,
+ print
+
+for k in distr:
+ bin_list = bin(k)
+ bin_num = sum([int(x) if x=='1' else 0 for x in bin_list[2:]])
+ sys.stderr.write("%d %028s %d \n" % (bin_num, bin_list[2:], distr[k]))
+
diff --git a/structure/activity/hamming_dist.py b/structure/activity/hamming_dist.py
new file mode 100644
index 0000000..2a98e64
--- /dev/null
+++ b/structure/activity/hamming_dist.py
@@ -0,0 +1,41 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the activity of the n-th
+## layer.
+##
+##
+
+
+import sys
+
+if len(sys.argv) < 3:
+ print "Usage: %s <layer1> <layer2> [<layer3>...]" % sys.argv[0]
+ sys.exit(1)
+
+max_N = -1
+
+layers = []
+
+
+for layer in sys.argv[1:]:
+ active = []
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+ active.extend([s,d])
+ active = set(active)
+ layers.append(active)
+
+for i in range(len(layers)):
+ for j in range(i+1, len(layers)):
+ s = layer[i] ^ layer[j]
+ print i, j, len(s)*1.0 / min(len(layer[i]) + len(layer[j]), max_N+1)
+
diff --git a/structure/activity/layer_activity.py b/structure/activity/layer_activity.py
new file mode 100644
index 0000000..bea0416
--- /dev/null
+++ b/structure/activity/layer_activity.py
@@ -0,0 +1,27 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the activity of the n-th
+## layer.
+##
+##
+
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+
+for layer in sys.argv[1:]:
+ active = []
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ active.extend([s,d])
+ active = set(active)
+ print len(active)
diff --git a/structure/activity/layer_activity_vectors.py b/structure/activity/layer_activity_vectors.py
new file mode 100644
index 0000000..f32418d
--- /dev/null
+++ b/structure/activity/layer_activity_vectors.py
@@ -0,0 +1,44 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the activity of the n-th
+## layer.
+##
+##
+
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+max_N = -1
+
+layers = []
+
+
+for layer in sys.argv[1:]:
+ active = []
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+ active.extend([s,d])
+ active = set(active)
+ layers.append(active)
+
+for l in layers:
+ for n in range(max_N+1):
+ if n in l:
+ print 1,
+ else:
+ print 0,
+ print
+
diff --git a/structure/activity/multiplexity.py b/structure/activity/multiplexity.py
new file mode 100644
index 0000000..9e038c9
--- /dev/null
+++ b/structure/activity/multiplexity.py
@@ -0,0 +1,41 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the activity of the n-th
+## layer.
+##
+##
+
+
+import sys
+
+if len(sys.argv) < 3:
+ print "Usage: %s <layer1> <layer2> [<layer3>...]" % sys.argv[0]
+ sys.exit(1)
+
+max_N = -1
+
+layers = []
+
+
+for layer in sys.argv[1:]:
+ active = []
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+ active.extend([s,d])
+ active = set(active)
+ layers.append(active)
+
+for i in range(len(layers)):
+ for j in range(i+1, len(layers)):
+ s = layer[i] & layer[j]
+ print i, j, len(s)*1.0 / (max_N+1)
+
diff --git a/structure/activity/node_activity.py b/structure/activity/node_activity.py
new file mode 100644
index 0000000..6410a68
--- /dev/null
+++ b/structure/activity/node_activity.py
@@ -0,0 +1,51 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the activity of the n-th node. We
+## assume that nodes are numbered sequentially, starting from 0, with
+## no gaps (i.e., missing nodes are treated as isolated nodes)
+##
+##
+##
+
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+node_activity = {}
+
+max_N = -1
+
+for layer in sys.argv[1:]:
+ active = []
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ active.extend([s,d])
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+ active = set(active)
+ for n in active:
+ if n in node_activity:
+ node_activity[n] += 1
+ else:
+ node_activity[n] = 1
+
+
+for n in range(max_N+1):
+ if n in node_activity:
+ print node_activity[n]
+ else:
+ print 0
+
+
+
+
diff --git a/structure/activity/node_activity_vectors.py b/structure/activity/node_activity_vectors.py
new file mode 100644
index 0000000..9839334
--- /dev/null
+++ b/structure/activity/node_activity_vectors.py
@@ -0,0 +1,68 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the degrees of the n-th node at
+## each layer, separated by a space, in the format:
+##
+## node1_deglay1 node1_deglay2 .... node1_deglayM
+## node2_deglay1 node2_deglay2 .... node2_deglayM
+## ..............................................
+## nodeN_deglay1 nodeN_deglay2 .... nodeN_deglayM
+##
+##
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+node_degrees = {}
+
+max_N = -1
+
+num_layer = 0
+
+for layer in sys.argv[1:]:
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+
+ if s in node_degrees:
+ if num_layer in node_degrees[s]:
+ node_degrees[s][num_layer] += 1
+ else:
+ node_degrees[s][num_layer] = 1
+ else:
+ node_degrees[s] = {}
+ node_degrees[s][num_layer] = 1
+
+ if d in node_degrees:
+ if num_layer in node_degrees[d]:
+ node_degrees[d][num_layer] += 1
+ else:
+ node_degrees[d][num_layer] = 1
+ else:
+ node_degrees[d] = {}
+ node_degrees[d][num_layer] = 1
+ num_layer += 1
+
+
+for n in range(max_N+1):
+ for i in range(num_layer):
+ if n in node_degrees:
+ if i in node_degrees[n]:
+ print 1,
+ else:
+ print 0,
+ else:
+ print 0,
+ print
diff --git a/structure/activity/node_degree_vectors.py b/structure/activity/node_degree_vectors.py
new file mode 100644
index 0000000..8863daa
--- /dev/null
+++ b/structure/activity/node_degree_vectors.py
@@ -0,0 +1,68 @@
+####
+##
+## Take as input the layers of a multiplex, and provide as output a
+## file where the n-th line contains the degrees of the n-th node at
+## each layer, separated by a space, in the format:
+##
+## node1_deglay1 node1_deglay2 .... node1_deglayM
+## node2_deglay1 node2_deglay2 .... node2_deglayM
+## ..............................................
+## nodeN_deglay1 nodeN_deglay2 .... nodeN_deglayM
+##
+##
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+node_degrees = {}
+
+max_N = -1
+
+num_layer = 0
+
+for layer in sys.argv[1:]:
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+
+ if s > max_N:
+ max_N = s
+ if d > max_N:
+ max_N = d
+
+ if s in node_degrees:
+ if num_layer in node_degrees[s]:
+ node_degrees[s][num_layer] += 1
+ else:
+ node_degrees[s][num_layer] = 1
+ else:
+ node_degrees[s] = {}
+ node_degrees[s][num_layer] = 1
+
+ if d in node_degrees:
+ if num_layer in node_degrees[d]:
+ node_degrees[d][num_layer] += 1
+ else:
+ node_degrees[d][num_layer] = 1
+ else:
+ node_degrees[d] = {}
+ node_degrees[d][num_layer] = 1
+ num_layer += 1
+
+
+for n in range(max_N+1):
+ for i in range(num_layer):
+ if n in node_degrees:
+ if i in node_degrees[n]:
+ print node_degrees[n][i],
+ else:
+ print 0,
+ else:
+ print 0,
+ print
diff --git a/structure/correlations/Makefile b/structure/correlations/Makefile
new file mode 100644
index 0000000..58301cb
--- /dev/null
+++ b/structure/correlations/Makefile
@@ -0,0 +1,15 @@
+CFLAGS="-O3"
+CC="gcc"
+GSL_FLAGS=-lgsl -lgslcblas -lm
+MFLAG=-lm
+
+all: dump_k_q fit_knn
+
+fit_knn: fit_knn.c fit_utils.c
+ $(CC) $(CFLAGS) -o fit_knn fit_knn.c fit_utils.c $(GSL_FLAGS)
+
+dump_k_q: dump_k_q.c rank_utils.c
+ $(CC) $(CFLAGS) -o dump_k_q dump_k_q.c rank_utils.c $(MFLAG)
+
+clean:
+ rm dump_k_q fit_knn
diff --git a/structure/correlations/compute_pearson.py b/structure/correlations/compute_pearson.py
new file mode 100644
index 0000000..e2c9055
--- /dev/null
+++ b/structure/correlations/compute_pearson.py
@@ -0,0 +1,33 @@
+####
+##
+## Compute the pearson correlation coefficient between the values of
+## node properties included in the two files provided as input.
+##
+
+import sys
+import numpy
+import scipy.stats
+import math
+
+if len(sys.argv) < 3:
+ print "Usage %s <file1> <file2>" % sys.argv[0]
+ sys.exit(1)
+
+
+x1 = []
+
+with open(sys.argv[1], "r") as lines:
+ for l in lines:
+ elem = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split()][0]
+ x1.append(elem)
+
+x2 = []
+
+with open(sys.argv[2], "r") as lines:
+ for l in lines:
+ elem = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split()][0]
+ x2.append(elem)
+
+
+r2 =numpy.corrcoef(x1,x2)
+print r2[0][1]
diff --git a/structure/correlations/compute_rho.py b/structure/correlations/compute_rho.py
new file mode 100644
index 0000000..8fa4174
--- /dev/null
+++ b/structure/correlations/compute_rho.py
@@ -0,0 +1,118 @@
+####
+##
+## Get two rankings and a pairing, and compute the corresponding
+## Spearman's rho rank correlation coefficient
+##
+##
+##
+##
+
+import sys
+import random
+import math
+#import scipy.stats
+
+##
+## Compute the constant C of the Spearman's rho correlation coefficient in the draft
+##
+def compute_C(rank1, rank2):
+ N = len(rank1)
+ [mu1, mu2] = [1.0 * sum(x.values()) / len(x.values()) for x in [rank1, rank2]]
+ [sum1, sum2] = [1.0 * sum(x.values()) for x in [rank1, rank2]]
+ C = N * mu1 * mu2 - mu2 * sum1 - mu1 * sum2
+
+ #print mu1, mu2, sum1, sum2, C
+
+ return C
+
+
+##
+## Compute the constant D of the Spearman's rho correlation coefficient in the draft
+##
+def compute_D(rank1, rank2):
+ [mu1, mu2] = [1.0 * sum(x.values()) / len(x.values()) for x in [rank1, rank2]]
+ s1 = sum([math.pow((x-mu1), 2) for x in rank1.values()])
+ s2 = sum([math.pow((x-mu2), 2) for x in rank2.values()])
+
+ D = math.sqrt(s1 * s2)
+ return D
+
+def compute_rho(rank1, rank2, pairing, C, D):
+
+ rho = 0
+ for s,t in pairing:
+ rho += rank1[s] * rank2[t]
+ rho = (rho + C) / D
+ return rho
+
+if len(sys.argv) < 3:
+ print "Usage: %s <rank1> <rank2> [<pairing>]" % sys.argv[0]
+ sys.exit(1)
+
+rank1 = {}
+rank2 = {}
+
+lines = open(sys.argv[1], "r").readlines()
+
+i = 0
+for l in lines:
+ if l[0] == "#" or l.strip(" \n").split(" ") == []: ## Strip comments and empty lines
+ continue
+ r = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split(" ")][0]
+ rank1[i] = r
+ i += 1
+
+
+lines = open(sys.argv[2], "r").readlines()
+
+i = 0
+for l in lines:
+ if l[0] == "#" or l.strip(" \n").split(" ") == []: ## Strip comments and empty lines
+ continue
+ r = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split(" ")][0]
+ rank2[i] = r
+ i += 1
+
+
+N1 = len(rank1)
+N2 = len(rank2)
+
+if (N1 != N2):
+ print "The two files must have the same number of nodes!!!!!"
+ sys.exit(2)
+
+
+
+C = compute_C(rank1, rank2)
+D = compute_D(rank1, rank2)
+
+
+## We start from a random pairing, and compute the corresponding value of rho
+
+pairing = []
+
+if len(sys.argv) > 3:
+ lines = open(sys.argv[3], "r").readlines()
+
+
+
+ for l in lines:
+ if l[0] == "#" or l.strip(" \n").split(" ") == []:
+ continue
+ p1, p2 = [int(x) for x in l.strip(" \n").split(" ")][:2]
+ pairing.append((p1, p2))
+
+else:
+ for i in range (N1):
+ pairing.append((i,i))
+
+
+if len(pairing) != N1:
+ print "ERROR !!! The pairing should be of the same length of the ranking files!!!"
+ sys.exit(1)
+
+rho = compute_rho(rank1, rank2, pairing, C, D)
+
+print rho
+
+
diff --git a/structure/correlations/compute_tau.py b/structure/correlations/compute_tau.py
new file mode 100644
index 0000000..3d1f804
--- /dev/null
+++ b/structure/correlations/compute_tau.py
@@ -0,0 +1,133 @@
+####
+##
+## Take as input two files, whose n^th line contains the ranking of
+## element n, and compute the Kendall's \tau_b rank correlation
+## coefficient
+##
+##
+
+import sys
+from numpy import *
+
+
+def kendalltau(x,y):
+ initial_sort_with_lexsort = True # if True, ~30% slower (but faster under profiler!) but with better worst case (O(n log(n)) than (quick)sort (O(n^2))
+ n = len(x)
+ temp = range(n) # support structure used by mergesort
+ # this closure recursively sorts sections of perm[] by comparing
+ # elements of y[perm[]] using temp[] as support
+ # returns the number of swaps required by an equivalent bubble sort
+ def mergesort(offs, length):
+ exchcnt = 0
+ if length == 1:
+ return 0
+ if length == 2:
+ if y[perm[offs]] <= y[perm[offs+1]]:
+ return 0
+ t = perm[offs]
+ perm[offs] = perm[offs+1]
+ perm[offs+1] = t
+ return 1
+ length0 = length / 2
+ length1 = length - length0
+ middle = offs + length0
+ exchcnt += mergesort(offs, length0)
+ exchcnt += mergesort(middle, length1)
+ if y[perm[middle - 1]] < y[perm[middle]]:
+ return exchcnt
+ # merging
+ i = j = k = 0
+ while j < length0 or k < length1:
+ if k >= length1 or (j < length0 and y[perm[offs + j]] <= y[perm[middle + k]]):
+ temp[i] = perm[offs + j]
+ d = i - j
+ j += 1
+ else:
+ temp[i] = perm[middle + k]
+ d = (offs + i) - (middle + k)
+ k += 1
+ if d > 0:
+ exchcnt += d;
+ i += 1
+ perm[offs:offs+length] = temp[0:length]
+ return exchcnt
+
+ # initial sort on values of x and, if tied, on values of y
+ if initial_sort_with_lexsort:
+ # sort implemented as mergesort, worst case: O(n log(n))
+ perm = lexsort((y, x))
+ else:
+ # sort implemented as quicksort, 30% faster but with worst case: O(n^2)
+ perm = range(n)
+ perm.sort(lambda a,b: cmp(x[a],x[b]) or cmp(y[a],y[b]))
+
+ # compute joint ties
+ first = 0
+ t = 0
+ for i in xrange(1,n):
+ if x[perm[first]] != x[perm[i]] or y[perm[first]] != y[perm[i]]:
+ t += ((i - first) * (i - first - 1)) / 2
+ first = i
+ t += ((n - first) * (n - first - 1)) / 2
+
+ # compute ties in x
+ first = 0
+ u = 0
+ for i in xrange(1,n):
+ if x[perm[first]] != x[perm[i]]:
+ u += ((i - first) * (i - first - 1)) / 2
+ first = i
+ u += ((n - first) * (n - first - 1)) / 2
+
+ # count exchanges
+ exchanges = mergesort(0, n)
+ # compute ties in y after mergesort with counting
+ first = 0
+ v = 0
+ for i in xrange(1,n):
+ if y[perm[first]] != y[perm[i]]:
+ v += ((i - first) * (i - first - 1)) / 2
+ first = i
+ v += ((n - first) * (n - first - 1)) / 2
+
+ tot = (n * (n - 1)) / 2
+ if tot == u and tot == v:
+ return 1 # Special case for all ties in both ranks
+
+ tau = ((tot-(v+u-t)) - 2.0 * exchanges) / (sqrt(float(( tot - u )) * float( tot - v )))
+
+ # what follows reproduces ending of Gary Strangman's original stats.kendalltau() in SciPy
+ svar = (4.0*n+10.0) / (9.0*n*(n-1))
+ z = tau / sqrt(svar)
+ ##prob = erfc(abs(z)/1.4142136)
+ ##return tau, prob
+ return tau
+
+def main():
+
+ if len(sys.argv) < 3:
+ print "Usage: %s <file1> <file2>" % sys.argv[0]
+ sys.exit(1)
+
+ x1 = []
+ x2= []
+
+ lines = open(sys.argv[1]).readlines()
+
+ for l in lines:
+ elem = [float(x) if "e" in x or "." in x else int(x) for x in l.strip(" \n").split()][0]
+ x1.append(elem)
+
+ lines = open(sys.argv[2]).readlines()
+
+ for l in lines:
+ elem = [float(x) if "e" in x or "." in x else int(x) for x in l.strip(" \n").split()][0]
+ x2.append(elem)
+
+
+ tau = kendalltau(x1,x2)
+ print tau
+
+
+if __name__ == "__main__":
+ main()
diff --git a/structure/correlations/dump_k_q.c b/structure/correlations/dump_k_q.c
new file mode 100644
index 0000000..2b6ba79
--- /dev/null
+++ b/structure/correlations/dump_k_q.c
@@ -0,0 +1,50 @@
+/**
+ *
+ * Get the degree distributions of two layers and a pairing, and dump
+ * on output the pairs k-q
+ *
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+void dump_k_q(double *R1, double *R2, int N, int *pairing){
+
+ int i;
+
+ for (i=0; i<N; i++){
+ printf("%g %g\n", R1[i], R2[pairing[i]]);
+ }
+}
+
+
+int main(int argc, char *argv[]){
+
+ int N1, N2;
+ double *R1 = NULL;
+ double *R2 = NULL;
+
+ int *pairing = NULL;
+
+
+ if (argc < 4){
+ printf("Usage: %s <degs1> <degs2> <pairing>\n", argv[0]);
+ exit(1);
+ }
+
+ load_ranking(argv[1], &N1, &R1);
+ load_ranking(argv[2], &N2, &R2);
+
+ if (N1 != N2){
+ printf("Error!!!! The two files must have the same number of lines!!!! Exiting...\n");
+ exit(1);
+ }
+
+ pairing = malloc(N1 * sizeof(int));
+
+ load_pairing(&pairing, N1, argv[3]);
+
+ dump_k_q(R1, R2, N1, pairing);
+
+}
diff --git a/structure/correlations/fit_knn.c b/structure/correlations/fit_knn.c
new file mode 100644
index 0000000..64fc4c3
--- /dev/null
+++ b/structure/correlations/fit_knn.c
@@ -0,0 +1,59 @@
+/**
+ *
+ * Take a file which contains pairs in the form:
+ *
+ * k q
+ *
+ * compute qnn(k) and fit this function with a power-law
+ *
+ * a*x^{b}
+ *
+ * The fit is actually performed as a linear fit on the
+ * exponential-binned log-log distribution
+ *
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <gsl/gsl_fit.h>
+#include "fit_utils.h"
+
+
+
+int main(int argc, char *argv[]){
+
+ double *k, *q, *x, *y;
+ int N, num;
+ double alpha;
+
+ double c0, c1, cov00, cov01, cov11, sqsum;
+
+ k=q=x=y=NULL;
+
+ if (argc < 3){
+ printf("%s <filein> <alpha>\n", argv[0]);
+ exit(1);
+ }
+
+ alpha=atof(argv[2]);
+
+ fprintf(stderr,"Alpha: %g\n", alpha);
+
+ load_sequence_col(argv[1], &k, &N, 0);
+ load_sequence_col(argv[1], &q, &N, 1);
+ //exit(1);
+ exp_bin_avg(k, q, N, alpha, &x, &y, &num);
+ //normalize_distr(x, y, num);
+ //dump_distr(x, y, num);
+ compact_distr(x, y, &num);
+ // dump_distr(x,y,num);
+ map_vec(x, num, log);
+ map_vec(y, num, log);
+ //dump_distr(x,y,num);
+
+ gsl_fit_linear(x, 1, y, 1, num, &c0, &c1, &cov00, &cov01, &cov11, &sqsum);
+ fprintf(stderr, "a: %g b: %g\n", exp(c0), c1);
+}
diff --git a/structure/correlations/fit_utils.c b/structure/correlations/fit_utils.c
new file mode 100644
index 0000000..e9e3954
--- /dev/null
+++ b/structure/correlations/fit_utils.c
@@ -0,0 +1,382 @@
+/**
+ *
+ * Fit a given sequence with a power-law function
+ *
+ * a*x^{b}
+ *
+ * The fit is actually performed as a linear fit on the
+ * exponential-binned log-log distribution
+ *
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <gsl/gsl_fit.h>
+
+/**
+ *
+ * Load a sequence from a file, which contains one element on each line
+ *
+ */
+
+void load_sequence(char *fname, double **v, int *N){
+
+ int size;
+ char buff[256];
+
+ FILE *f;
+
+ f = fopen(fname, "r");
+ if (!f){
+ fprintf(stderr, "Error opening file %s!!!! Exiting...\n", fname);
+ exit(1);
+ }
+
+ *N =0;
+ size = 10;
+ if (*v == NULL)
+ *v = malloc(size * sizeof(double));
+ else{
+ *v = realloc(*v, size * sizeof(double));
+ }
+
+ while (fgets(buff, 255, f)){
+ (*v)[*N] = atof(buff);
+ *N += 1;
+ if (*N == size){
+ size += 10;
+ *v = realloc(*v, size * sizeof(double));
+ }
+ }
+ *v = realloc(*v, (*N) * sizeof(double));
+ fclose(f);
+}
+
+
+/**
+ *
+ * Load a sequence, getting the "col"-th column of the input file
+ *
+ */
+
+void load_sequence_col(char *fname, double **v, int *N, int col){
+
+ int size, n;
+ char buff[256];
+ char *ptr;
+
+ FILE *f;
+
+ f = fopen(fname, "r");
+ if (!f){
+ fprintf(stderr, "Error opening file %s!!!! Exiting...\n", fname);
+ exit(1);
+ }
+
+ *N =0;
+ size = 10;
+ if (*v == NULL)
+ *v = malloc(size * sizeof(double));
+ else{
+ *v = realloc(*v, size * sizeof(double));
+ }
+
+ while (fgets(buff, 255, f)){
+ ptr = strtok(buff, " ");
+ if (col > 0){
+ n = 0;
+ while (n<col){
+ ptr = strtok(NULL, " ");
+ n += 1;
+ }
+ }
+ (*v)[*N] = atof(ptr);
+ *N += 1;
+
+ if (*N == size){
+ size += 10;
+ *v = realloc(*v, size * sizeof(double));
+ }
+ }
+ *v = realloc(*v, (*N) * sizeof(double));
+ fclose(f);
+}
+
+
+/**
+ *
+ * Make the exponential binning of a distribution, with a giving
+ * exponent alpha. The value of y[i] is the number of elements of v
+ * whose value lies between x[i-1] and x[i]...
+ *
+ */
+
+void exp_bin_cnt(double *v, int N, double alpha, double **x, double **y, int *num){
+
+ double min_v, max_v, val, last_val;
+ int i, j, size, num_x;
+ double last_size, new_size;
+
+ min_v = max_v = v[0];
+
+ for (i=1; i<N; i ++){
+ if (v[i] > max_v)
+ max_v = v[i];
+ else if (v[i] > 0 && v[i] < min_v)
+ min_v = v[i];
+ }
+
+ size = 10;
+ if (*x == NULL){
+ *x = malloc(size * sizeof(double));
+ }
+ else{
+ *x = realloc(*x, size * sizeof(double));
+ }
+
+ val = min_v;
+ last_size = min_v;
+ (*x)[0] = min_v;
+ num_x = 1;
+
+ while(val < max_v){
+ new_size = last_size * alpha;
+ val = last_size + new_size;
+ last_size = new_size;
+ last_val = val;
+ (*x)[num_x] = val;
+ num_x += 1;
+ if (num_x == size){
+ size += 10;
+ *x = realloc(*x, size * sizeof(double));
+ }
+ }
+
+ if (*y == NULL){
+ *y = malloc(num_x * sizeof(double));
+ }
+ else{
+ *y = realloc(*y, num_x * sizeof(double));
+ }
+ for (i=0; i < num_x; i++){
+ (*y)[i] = 0;
+ }
+
+
+
+ for(i=0; i <N; i ++){
+ j = 0;
+ while(v[i] > (*x)[j]){
+ j ++;
+ }
+ (*y)[j] += 1;
+ }
+ *num = num_x;
+}
+
+/**
+ *
+ * Make the exponential binning of a distribution, with a giving
+ * exponent alpha. The value of y[i] is the average of the values in
+ * the vector "w" whose corresponding v lies between x[i-1] and x[i]...
+ *
+ */
+
+
+void exp_bin_avg(double *v, double *w, int N, double alpha, double **x, double **y, int *num){
+
+ double min_v, max_v, val, last_val;
+ int i, j, size, num_x;
+ double last_size, new_size;
+ int *cnt;
+
+
+ min_v = max_v = v[0];
+
+ for (i=1; i<N; i ++){
+ if (v[i] > max_v)
+ max_v = v[i];
+ else if (v[i] > 0 && v[i] < min_v)
+ min_v = v[i];
+ }
+
+ size = 10;
+ if (*x == NULL){
+ *x = malloc(size * sizeof(double));
+ }
+ else{
+ *x = realloc(*x, size * sizeof(double));
+ }
+
+ val = min_v;
+ last_size = min_v;
+ (*x)[0] = min_v;
+ num_x = 1;
+
+ while(val < max_v){
+ new_size = last_size * alpha;
+ val = last_size + new_size;
+ last_size = new_size;
+ last_val = val;
+ (*x)[num_x] = val;
+ num_x += 1;
+ if (num_x == size){
+ size += 10;
+ *x = realloc(*x, size * sizeof(double));
+ }
+ }
+
+
+ cnt = malloc(num_x * sizeof(int));
+
+ if (*y == NULL){
+ *y = malloc(num_x * sizeof(double));
+ }
+ else{
+ *y = realloc(*y, num_x * sizeof(double));
+ }
+ for (i=0; i < num_x; i++){
+ (*y)[i] = 0;
+ cnt[i] = 0;
+ }
+
+ for(i=0; i <N; i ++){
+ j = 0;
+ while(j < num_x && v[i] > (*x)[j]){
+ j ++;
+ }
+ if(j == num_x){
+ printf("Error!!!!! Trying to assing a non-existing bin!!! -- fit_utils.c:exp_bin_avg!!!\n");
+ exit(37);
+ }
+ (*y)[j] += w[i];
+ cnt[j] += 1;
+ }
+ *num = num_x;
+
+ for(i=0; i< num_x; i++){
+ if (cnt[i] > 0){
+ (*y)[i] = (*y)[i] / cnt[i];
+ }
+ }
+ free(cnt);
+}
+
+
+/**
+ *
+ * Print a distribution on stdout
+ *
+ */
+
+void dump_distr(double *x, double *y, int N){
+ int i;
+
+ for(i=0; i<N; i ++){
+ printf("%g %g\n", x[i], y[i]);
+ }
+
+}
+
+
+/**
+ * Compact a distribution, leaving only the pairs (x_i, y_i) for which
+ * y_i > 0
+ *
+ */
+
+void compact_distr(double *x, double *y, int *num){
+
+ int i, j;
+
+ i = j = 0;
+ while(j < *num){
+ while(j < *num && y[j] == 0){
+ j ++;
+ }
+ if (j==*num){
+ break;
+ }
+ x[i] = x[j];
+ y[i] = y[j];
+ i ++;
+ j ++;
+ }
+ *num = i;
+}
+
+
+/**
+ *
+ * Apply the function "f" on all the elemnts of a vector, in-place
+ *
+ */
+
+void map_vec(double *v, int N, double (*f)(double)){
+ int i;
+
+ for (i=0; i<N; i ++){
+ v[i] = f(v[i]);
+ }
+}
+
+
+/**
+ *
+ * Normalize a distribution, dividing each y[i] for the width of the
+ * corresponding bin (i.e., x[i] - x[i-1])
+ *
+ */
+void normalize_distr(double *x, double *y, int num){
+
+ int i;
+
+ for(i=1; i<num; i ++){
+ y[i] /= (x[i] - x[i-1]);
+ }
+}
+
+/**
+ *
+ * take two vectors (k and q) and a pairing, and compute the best
+ * power-law fit "a*k^{mu}" of qnn(k)
+ *
+ */
+
+void fit_current_trend(double *k, double *q, int N, int *pairing, double *mu, double *a,
+ double *corr){
+
+ static int *num;
+ static double *q_pair, *x, *y;
+
+ int i;
+ int fit_num;
+ double c0, c1, cov00, cov01, cov11, sqsum;
+
+ if (q_pair == NULL){
+ q_pair = malloc(N * sizeof(double));
+ }
+
+ for(i=0; i<N; i++){
+ q_pair[i] = q[pairing[i]];
+ }
+
+ exp_bin_avg(k, q_pair, N, 1.3, &x, &y, &fit_num);
+ //normalize_distr(x,y, fit_num);
+ compact_distr(x, y, &fit_num);
+
+ map_vec(x, fit_num, log);
+ map_vec(y, fit_num, log);
+ gsl_fit_linear(x, 1, y, 1, fit_num, &c0, &c1, &cov00, &cov01, &cov11, &sqsum);
+
+ //fprintf(stderr,"cov00: %g cov01: %g cov11: %g\n", cov00, cov01, cov11);
+ //fprintf(stderr, "corr: %g ", cov01 / sqrt(cov00 * cov11));
+ *a = c0;
+ *mu = c1;
+ *corr = cov01/sqrt(cov00 * cov11);
+}
+
diff --git a/structure/correlations/fit_utils.h b/structure/correlations/fit_utils.h
new file mode 100644
index 0000000..183f47c
--- /dev/null
+++ b/structure/correlations/fit_utils.h
@@ -0,0 +1,25 @@
+#ifndef __FIT_UTILS_H__
+#define __FIT_UTILS_H__
+
+
+void load_sequence(char *fname, double **v, int *N);
+
+void load_sequence_col(char *fname, double **v, int *N, int col);
+
+void exp_bin_cnt(double *v, int N, double alpha, double **x, double **y, int *num);
+
+void exp_bin_avg(double *v, double *w, int N, double alpha, double **x, double **y, int *num);
+
+void dump_distr(double *x, double *y, int N);
+
+void compact_distr(double *x, double *y, int *num);
+
+void map_vec(double *v, int N, double (*f)(double));
+
+void normalize_distr(double *x, double *y, int num);
+
+void fit_current_trend(double *k, double *q, int N, int *pairing,
+ double *mu, double *a, double *corr);
+
+
+#endif /* __FIT_UTILS_H__ */
diff --git a/structure/correlations/knn_q_from_degrees.py b/structure/correlations/knn_q_from_degrees.py
new file mode 100644
index 0000000..5e10bfa
--- /dev/null
+++ b/structure/correlations/knn_q_from_degrees.py
@@ -0,0 +1,49 @@
+####
+##
+##
+## Take a file with the degree sequence of nodes which are active on
+## both layers, and compute \overline{k}(q) and \overline{q}(k),
+## i.e. the two inter-layer correlation functions, where we call "k"
+## the degree on the first column, and "q" the degree on the second
+## column
+##
+##
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <filein>" % sys.argv[0]
+ sys.exit(1)
+
+qnn_k = {}
+knn_q = {}
+
+with open(sys.argv[1]) as f:
+ for l in f:
+ k, q = [int(x) for x in l.strip(" \n").split(" ")]
+ if qnn_k.has_key(k):
+ qnn_k[k].append(q)
+ else:
+ qnn_k[k] = [q]
+ if knn_q.has_key(q):
+ knn_q[q].append(k)
+ else:
+ knn_q[q] = [k]
+
+
+
+
+keys= qnn_k.keys()
+keys.sort()
+for k in keys:
+ print k, 1.0 * sum(qnn_k[k])/len(qnn_k[k])
+
+
+
+
+keys= knn_q.keys()
+keys.sort()
+for q in keys:
+ sys.stderr.write("%d %f\n" % (q, 1.0 * sum(knn_q[q])/len(knn_q[q])))
+
+
diff --git a/structure/correlations/knn_q_from_layers.py b/structure/correlations/knn_q_from_layers.py
new file mode 100644
index 0000000..c108c18
--- /dev/null
+++ b/structure/correlations/knn_q_from_layers.py
@@ -0,0 +1,98 @@
+####
+##
+##
+## Compute the degree-degree correlations of a multiplex graph, namely:
+##
+## <k_1>(k_2) and <k_2>(k_1)
+##
+## Takes as input the two lists of edges corresponding to each layer
+##
+
+import sys
+import numpy as np
+import networkx as net
+
+
+def knn(G, n):
+ neigh = G.neighbors(n)
+ l = G.degree(neigh).values()
+ return 1.0 * sum(l) / len(l)
+
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> <layer2>" % sys.argv[0]
+ sys.exit(1)
+
+
+G1 = net.read_edgelist(sys.argv[1])
+G2 = net.read_edgelist(sys.argv[2])
+
+k1_k1 = {} ## Intraleyer knn (k1)
+k2_k2 = {} ## Intralayer knn (k2)
+k1_k2 = {} ## Interlayer average degree at layer 1 of a node having degree k_2 in layer 2
+k2_k1 = {} ## Interlayer average degree at layer 2 of a node having degree k_1 in layer 1
+
+
+for n in G1.nodes():
+ k1 = G1.degree(n)
+
+
+ ##print k1,k2
+
+ knn1 = knn(G1, n)
+ if n in G2.nodes():
+ k2 = G2.degree(n)
+ knn2 = knn(G2, n)
+ else:
+ k2 = 0
+ knn2 = 0
+
+ if k1_k1.has_key(k1):
+ k1_k1[k1].append(knn1)
+ else:
+ k1_k1[k1] = [knn1]
+
+ if k2_k2.has_key(k2):
+ k2_k2[k2].append(knn2)
+ else:
+ k2_k2[k2] = [knn2]
+
+
+ if k1_k2.has_key(k2):
+ k1_k2[k2].append(k1)
+ else:
+ k1_k2[k2] = [k1]
+
+ if k2_k1.has_key(k1):
+ k2_k1[k1].append(k2)
+ else:
+ k2_k1[k1] = [k2]
+
+
+k1_keys = k1_k1.keys()
+k1_keys.sort()
+k2_keys = k2_k2.keys()
+k2_keys.sort()
+
+
+f = open("%s_%s_k1" % (sys.argv[1], sys.argv[2]), "w+")
+
+for n in k1_keys:
+ avg_knn = np.mean(k1_k1[n])
+ std_knn = np.std(k1_k1[n])
+ avg_k2 = np.mean(k2_k1[n])
+ std_k2 = np.std(k2_k1[n])
+ f.write("%d %f %f %f %f\n" % (n, avg_knn, std_knn, avg_k2, std_k2))
+
+f.close()
+
+f = open("%s_%s_k2" % (sys.argv[1], sys.argv[2]), "w+")
+
+for n in k2_keys:
+ avg_knn = np.mean(k2_k2[n])
+ std_knn = np.std(k2_k2[n])
+ avg_k1 = np.mean(k1_k2[n])
+ std_k1 = np.std(k1_k2[n])
+ f.write("%d %f %f %f %f\n" % (n, avg_knn, std_knn, avg_k1, std_k1))
+f.close()
+
diff --git a/structure/correlations/rank_nodes.py b/structure/correlations/rank_nodes.py
new file mode 100644
index 0000000..6985e88
--- /dev/null
+++ b/structure/correlations/rank_nodes.py
@@ -0,0 +1,74 @@
+####
+##
+##
+## Get a file as input, whose n-th line corresponds to the value of a
+## certain property of node n, and rank nodes according to their
+## properties, taking into account ranking ties properly.
+##
+## The output is a file whose n-th line is the "ranking" of the n-th
+## node according to the given property. (notice that rankings could
+## be fractional, due to the tie removal algorithm)
+##
+##
+
+
+import sys
+import math
+
+
+if len(sys.argv) < 2:
+ print "Usage: %s <filein>" % sys.argv[0]
+ sys.exit(1)
+
+
+
+lines = open(sys.argv[1], "r").readlines()
+
+ranking = []
+
+n=0
+for l in lines:
+ if l[0] == "#" or l.strip(" \n").split(" ") == []:
+ continue
+ v = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split(" ")][0]
+ ranking.append((v,n))
+ n +=1
+
+ranking.sort(reverse=True)
+#print ranking
+new_ranking = {}
+
+v0, n0 = ranking[0]
+
+old_value = v0
+tot_rankings = 1
+
+stack = [n0]
+l=1.0
+for v,n in ranking[1:]:
+ l += 1
+ ##print stack, tot_rankings
+ if v != old_value: ### There is a new rank
+ # We first compute the rank for all the nodes in the stack and then we set it
+ new_rank_value = 1.0 * tot_rankings / len(stack)
+ ##print new_rank_value
+ for j in stack:
+ new_ranking[j] = new_rank_value
+ old_value = v
+ tot_rankings = l
+ stack = [n]
+ else: # One more value with the same rank, keep it for the future
+ stack.append(n)
+ tot_rankings += l
+
+new_rank_value = 1.0 * tot_rankings / len(stack)
+#print new_rank_value
+for j in stack:
+ new_ranking[j] = new_rank_value
+
+#print new_ranking
+
+keys = new_ranking.keys()
+keys.sort()
+for k in keys:
+ print new_ranking[k]
diff --git a/structure/correlations/rank_nodes_thresh.py b/structure/correlations/rank_nodes_thresh.py
new file mode 100644
index 0000000..44b565a
--- /dev/null
+++ b/structure/correlations/rank_nodes_thresh.py
@@ -0,0 +1,87 @@
+####
+##
+##
+## Get a file as input, whose n-th line corresponds to the value of a
+## certain property of node n, and rank nodes according to their
+## properties, taking into account ranking ties properly.
+##
+## The output is a file whose n-th line is the "ranking" of the n-th
+## node according to the given property. (notice that rankings could
+## be fractional, due to the tie removal algorithm)
+##
+## The rank of a node is set to "0" (ZERO) if the corresponding
+## property is smaller than a value given as second parameter
+##
+
+
+import sys
+import math
+
+
+if len(sys.argv) < 3:
+ print "Usage: %s <filein> <thresh>" % sys.argv[0]
+ sys.exit(1)
+
+
+thresh = float(sys.argv[2])
+
+lines = open(sys.argv[1], "r").readlines()
+
+ranking = []
+
+n=0
+for l in lines:
+ if l[0] == "#" or l.strip(" \n").split(" ") == []:
+ continue
+ v = [float(x) if "." in x or "e" in x else int(x) for x in l.strip(" \n").split(" ")][0]
+ if v >= thresh:
+ ranking.append((v,n))
+ else:
+ ranking.append((0,n))
+ n +=1
+
+ranking.sort(reverse=True)
+#print ranking
+new_ranking = {}
+
+v0, n0 = ranking[0]
+
+
+old_value = v0
+tot_rankings = 1
+
+stack = [n0]
+l=1.0
+for v,n in ranking[1:]:
+ l += 1
+ ##print stack, tot_rankings
+ if v != old_value: ### There is a new rank
+ # We first compute the rank for all the nodes in the stack and then we set it
+ if old_value == 0:
+ new_rank_value = 0
+ else:
+ new_rank_value = 1.0 * tot_rankings / len(stack)
+ ##print new_rank_value
+ for j in stack:
+ new_ranking[j] = new_rank_value
+ old_value = v
+ tot_rankings = l
+ stack = [n]
+ else: # One more value with the same rank, keep it for the future
+ stack.append(n)
+ tot_rankings += l
+
+if v == 0 :
+ new_rank_value = 0
+else:
+ new_rank_value = 1.0 * tot_rankings / len(stack)
+#print new_rank_value
+for j in stack:
+ new_ranking[j] = new_rank_value
+
+#print new_ranking
+
+keys = new_ranking.keys()
+keys.sort()
+for k in keys:
+ print new_ranking[k]
diff --git a/structure/correlations/rank_occurrence.py b/structure/correlations/rank_occurrence.py
new file mode 100644
index 0000000..6339d5d
--- /dev/null
+++ b/structure/correlations/rank_occurrence.py
@@ -0,0 +1,45 @@
+####
+##
+## Get two rankings and compute the size of the k-intersection,
+## i.e. the number of elements which are present in the first k
+## positions of both rankings, as a function of k
+##
+##
+
+
+import sys
+
+
+if len(sys.argv)< 4:
+ print "Usage: %s <file1> <file2> <increment>" % sys.argv[0]
+ sys.exit(1)
+
+incr = int(sys.argv[3])
+
+rank1 = []
+rank2 = []
+
+lines = open(sys.argv[1], "r").readlines()
+
+for l in lines:
+ n= l.strip(" \n").split(" ")[0]
+ rank1.append(n)
+
+lines = open(sys.argv[2], "r").readlines()
+
+for l in lines:
+ n= l.strip(" \n").split(" ")[0]
+ rank2.append(n)
+
+
+N = len(rank1)
+
+i = incr
+
+while i < N+incr:
+ l = len(set(rank1[:i]) & set(rank2[:i]))
+ print i, l
+ i += incr
+
+
+
diff --git a/structure/correlations/rank_utils.c b/structure/correlations/rank_utils.c
new file mode 100644
index 0000000..4b8ef41
--- /dev/null
+++ b/structure/correlations/rank_utils.c
@@ -0,0 +1,217 @@
+/**
+ *
+ * Some utility functions to be used with tune_rho, compute_rho and
+ * the like
+ *
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <time.h>
+
+#include "rank_utils.h"
+
+void load_ranking(char *fname, int *N, double **R){
+
+ char buff[256];
+ int size = 10;
+ FILE *f;
+
+ if (*R == NULL){
+ *R = malloc(size * sizeof (double));
+ }
+ f = fopen(fname, "r");
+ if (!f){
+ printf("Unable to open file: %s!!! Exiting\n");
+ exit(1);
+ }
+
+ *N = 0;
+
+ while (fgets(buff, 255, f)){
+ if (* N == size){
+ size += 10;
+ *R = realloc(*R, size * sizeof(double));
+ }
+ (*R)[*N] = atof(buff);
+ *N += 1;
+ }
+}
+
+double avg_array(double *v, int N){
+
+ double sum = 0.0;
+ int i;
+
+ for (i = 0; i < N; i ++){
+ sum += v[i];
+ }
+ return sum/N;
+}
+
+double compute_C(double *R1, double *R2, int N){
+ double mu1, mu2, sum1, sum2;
+
+ mu1 = avg_array(R1, N);
+ mu2 = avg_array(R2, N);
+
+ sum1 = mu1 * N;
+ sum2 = mu2 * N;
+
+ return N * mu1 * mu2 - mu2 * sum1 - mu1 * sum2;
+}
+
+
+double compute_D(double *R1, double *R2, int N){
+
+ double mu1, mu2, s1, s2;
+ int i;
+
+ mu1 = avg_array(R1, N);
+ mu2 = avg_array(R1, N);
+
+ s1 = s2 = 0.0;
+
+ for (i=0 ; i < N; i ++){
+ s1 += pow((R1[i] - mu1), 2);
+ s2 += pow((R2[i] - mu2), 2);
+ }
+
+ return sqrt(s1 * s2);
+}
+
+
+double compute_rho(double *R1, double *R2, int N, int *pairing){
+
+ double rho = 0;
+ int i;
+
+ for (i=0; i < N; i ++){
+ rho += R1[i] * R2[ pairing[i] ];
+ }
+
+ rho = (rho + compute_C(R1, R2, N))/ compute_D(R1, R2, N);
+ return rho;
+}
+
+void dump_ranking(double *R, int N){
+ int i;
+
+ for (i=0; i < N; i ++){
+ printf("%d: %f\n", i, R[i] );
+ }
+}
+
+
+void init_pairing_natural(int *pairing, int N){
+ int i;
+
+ for (i = 0; i< N; i ++){
+ pairing[i] = i;
+ }
+}
+
+void init_pairing_inverse(int *pairing, int N){
+ int i;
+
+ for (i = 0; i< N; i ++){
+ pairing[i] = N-i-1;
+ }
+}
+
+void select_pairing(int *pairing, int N, int argc, char *argv[], int pos){
+
+ if (argc < pos + 1 || !strncasecmp("rnd", argv[pos], 3)){
+ init_pairing_random(pairing, N);
+ }
+ else if (!strncasecmp("nat", argv[pos], 3)){
+ init_pairing_natural(pairing, N);
+ }
+ else if (!strncasecmp("inv", argv[pos], 3)){
+ init_pairing_inverse(pairing, N);
+ }
+ else{
+ printf ("Pairing strategy \"%s\" unknown!!! Exiting...\n", argv[pos]);
+ exit(1);
+ }
+
+}
+
+
+void shuffle_sequence(int *s, int N){
+
+ int i, j, tmp;
+
+ for (i=N-1; i>=0; i--){
+ j = rand() % (i+1);
+ tmp = s[j];
+ s[j] = s[i];
+ s[i] = tmp;
+ }
+}
+
+
+
+void init_pairing_random(int *pairing, int N){
+
+ init_pairing_natural(pairing, N);
+ shuffle_sequence(pairing, N);
+
+}
+
+/* Loads a pairing from a file, in the format:
+ *
+ * rank1 rank2
+ * ...........
+ */
+void load_pairing(int **pairing, int N, char *fname){
+
+ FILE *f;
+ int i, j, num;
+ char buff[256];
+ char *ptr;
+
+ f = fopen(fname, "r");
+ if (!f){
+ printf("Error opening file \"%s\"!!!! Exiting....\n", fname);
+ exit(2);
+ }
+
+ if (*pairing == NULL){
+ *pairing = malloc(N * sizeof(int));
+ init_pairing_natural(*pairing, N);
+ }
+
+ num = 0;
+ while(num < N){
+ fgets(buff, 255, f);
+ if (buff[0] == '#')
+ continue;
+ ptr = strtok(buff, " "); /* read the first node */
+ i = atoi(ptr);
+ ptr = strtok(NULL, " "); /* read the second node */
+ j = atoi(ptr);
+ (*pairing)[i] = j;
+ num += 1;
+ }
+}
+
+
+void dump_pairing(int *pairing, int N){
+ int i;
+
+ for(i=0; i< N; i ++){
+ printf("%d %d\n", i, pairing[i]);
+ }
+}
+
+void copy_pairing(int *p1, int *p2, int N){
+ int i;
+
+ for (i=0 ; i < N; i ++){
+ p2[i] = p1[i];
+ }
+}
diff --git a/structure/correlations/rank_utils.h b/structure/correlations/rank_utils.h
new file mode 100644
index 0000000..521582a
--- /dev/null
+++ b/structure/correlations/rank_utils.h
@@ -0,0 +1,44 @@
+#ifndef __RANK_UTILS_H__
+#define __RANK_UTILS_H__
+
+/* Load a ranking */
+void load_ranking(char *fname, int *N, double **R);
+
+/* Compute the average of an array */
+double avg_array(double *v, int N);
+
+/* Compute the term "C" in the rho correlation coefficient */
+double compute_C(double *R1, double *R2, int N);
+
+/* Compute the term "D" in the rho correlation coefficient */
+double compute_D(double *R1, double *R2, int N);
+
+/* Compute the Spearman's rank correlation coefficient \rho */
+double compute_rho(double *R1, double *R2, int N, int *pairing);
+
+void dump_ranking(double *R, int N);
+
+void init_pairing_natural(int *pairing, int N);
+
+void init_pairing_inverse(int *pairing, int N);
+
+void init_pairing_random(int *pairing, int N);
+
+void select_pairing(int *pairing, int N, int argc, char *argv[], int pos);
+
+void shuffle_sequence(int *s, int N);
+
+void dump_pairing(int *pairing, int N);
+
+void copy_pairing(int *pairing1, int *pairing2, int N);
+
+
+#endif /* __RANK_UTILS_H__ */
+
+
+
+
+
+
+
+
diff --git a/structure/metrics/aggregate_layers_w.py b/structure/metrics/aggregate_layers_w.py
new file mode 100644
index 0000000..cd5ea1d
--- /dev/null
+++ b/structure/metrics/aggregate_layers_w.py
@@ -0,0 +1,37 @@
+####
+##
+## Aggregate the layers of a multiplex, weigthing each edge according
+## to the number of times it occurs in the different layers
+##
+
+import sys
+import networkx as net
+
+
+if len(sys.argv) < 3:
+ print "Usage: %s <file1> <file2> [<file3>....]" % sys.argv(0)
+ sys.exit(1)
+
+G = net.Graph()
+
+lines = open(sys.argv[1], "r").readlines()
+
+for l in lines:
+ s,d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ G.add_edge(s,d)
+ G[s][d]['weigth'] = 1
+
+for f in sys.argv[2:]:
+ lines = open(f, "r").readlines()
+ for l in lines:
+ s,d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if (G.has_edge(s,d)):
+ G[s][d]['weigth'] += 1
+ else:
+ G.add_edge(s,d)
+ G[s][d]['weigth'] = 1
+
+for s,d in G.edges():
+ print s,d, G[s][d]['weigth']
+
+
diff --git a/structure/metrics/avg_edge_overlap.py b/structure/metrics/avg_edge_overlap.py
new file mode 100644
index 0000000..e68fd8b
--- /dev/null
+++ b/structure/metrics/avg_edge_overlap.py
@@ -0,0 +1,47 @@
+####
+##
+## Compute the average edge overlap of a multiplex, i.e. the average
+## number of layers in which an edge is present
+##
+##
+
+import sys
+
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+max_N = -1
+
+all_edges = {}
+
+layer_ID = -1
+
+for layer in sys.argv[1:]:
+ layer_ID += 1
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if s > d:
+ tmp = s
+ s = d
+ d = tmp
+ if (s,d) in all_edges:
+ all_edges[(s,d)].append(layer_ID)
+ else:
+ all_edges[(s,d)] = [layer_ID]
+
+K = len(all_edges.keys())
+M = len(sys.argv) - 1
+
+numer = 0
+
+for k in all_edges.keys():
+ numer += len(set(all_edges[(k)]))
+
+
+print 1.0 * numer / K, 1.0 * numer / (K * M)
+
diff --git a/structure/metrics/cartography_from_columns.py b/structure/metrics/cartography_from_columns.py
new file mode 100644
index 0000000..a2343ed
--- /dev/null
+++ b/structure/metrics/cartography_from_columns.py
@@ -0,0 +1,44 @@
+####
+##
+## Make a cartography (i.e., sum and participation coefficient) based
+## on the values found at the given column numbers of the file given
+## as input, e.g.:
+##
+## cartography_cols.py FILEIN 1 5 8 14
+##
+## makes the assumption that the system has 4 layers, and that the
+## quantities involved in the cartography are at the 2nd, 6th, 9th and
+## 15th column of FILEIN
+##
+##
+##
+
+import sys
+
+if len(sys.argv) < 4:
+ print "Usage: %s <filein> <col1> <col2> [<col3> <col4>...]" % sys.argv[0]
+ sys.exit(1)
+
+filein=sys.argv[1]
+cols = [int(x) for x in sys.argv[2:]]
+
+M = len(cols)
+
+with open(filein,"r") as f:
+ for l in f:
+ if l[0] == "#":
+ continue
+ elems = [float(x) if "e" in x or "." in x else int(x) for x in l.strip(" \n").split(" ")]
+ sum_elems = 0
+ part = 0
+ for c in cols:
+ val = elems[c]
+ sum_elems += val
+ part += val*val
+ if sum_elems > 0:
+ part = M * 1.0 / (M -1) * (1 - part * 1.0 / (sum_elems * sum_elems))
+ else:
+ part = 0.0
+ print elems[0], sum_elems, part
+
+
diff --git a/structure/metrics/cartography_from_deg_vectors.py b/structure/metrics/cartography_from_deg_vectors.py
new file mode 100644
index 0000000..c40d701
--- /dev/null
+++ b/structure/metrics/cartography_from_deg_vectors.py
@@ -0,0 +1,37 @@
+####
+##
+## Take as input a file containing, on each line, the degree vector of
+## a node of the multiplex, and compute the multiplex cartography
+## diagram
+##
+##
+
+import sys
+
+if len(sys.argv) < 2:
+ print "Usage: %s <node_deg_vectors>" % sys.argv[0]
+ sys.exit(1)
+
+filein=sys.argv[1]
+
+M = -1
+
+with open(filein,"r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+ elems = [int(x) for x in l.strip(" \n").split(" ")]
+ if (M == -1):
+ M = len(elems)
+ sum_elems = 0
+ part = 0
+ for val in elems:
+ sum_elems += val
+ part += val*val
+ if sum_elems > 0:
+ part = M * 1.0 / (M -1) * (1 - part * 1.0 / (sum_elems * sum_elems))
+ else:
+ part = 0.0
+ print sum_elems, part
+
+
diff --git a/structure/metrics/cartography_from_layers.py b/structure/metrics/cartography_from_layers.py
new file mode 100644
index 0000000..b9c4584
--- /dev/null
+++ b/structure/metrics/cartography_from_layers.py
@@ -0,0 +1,54 @@
+####
+##
+## Compute the participation coefficient of each node of a multiplex
+##
+##
+
+import sys
+import networkx as net
+import collections
+from compiler.ast import flatten
+
+
+
+if len(sys.argv) < 3:
+ print "Usage: %s <layer1> <layer2> [<layer3>...]" % sys.argv[0]
+ sys.exit(1)
+
+
+layers = []
+
+for f in sys.argv[1:]:
+ G = net.Graph(net.read_edgelist(f, nodetype=int))
+ layers.append(G)
+
+nodes = flatten([x for j in layers for x in j.nodes()])
+#nodes.sort()
+nodes = set(nodes)
+
+M = len(layers)
+
+#print nodes
+
+for n in nodes:
+ deg_alpha_square = 0
+ deg = 0
+ col = 0
+ print n,
+ for l in layers:
+ val = l.degree([n])
+ if not val:
+ col = 2* col
+ continue
+ col *= 2
+ col += 1
+ val = val[n]
+ deg += val
+ deg_alpha_square += val*val
+ if deg > 0:
+ print deg, 1.0 * M / (M-1) * (1.0 - 1.0 * deg_alpha_square/(deg * deg)) , col
+ else:
+ print 0 , 0, col
+
+
+
diff --git a/structure/metrics/edge_overlap.py b/structure/metrics/edge_overlap.py
new file mode 100644
index 0000000..1b27f55
--- /dev/null
+++ b/structure/metrics/edge_overlap.py
@@ -0,0 +1,43 @@
+####
+##
+## Compute the edge overlap of each edge of the multiplex, i.e. the
+## number of times that each edge appears in the multiplex
+##
+##
+
+import sys
+
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> [<layer2>...]" % sys.argv[0]
+ sys.exit(1)
+
+max_N = -1
+
+all_edges = {}
+
+layer_ID = -1
+
+for layer in sys.argv[1:]:
+ layer_ID += 1
+ with open(layer, "r") as lines:
+ for l in lines:
+ if l[0] == "#":
+ continue
+ s, d = [int(x) for x in l.strip(" \n").split(" ")[:2]]
+ if s > d:
+ tmp = s
+ s = d
+ d = tmp
+ if (s,d) in all_edges:
+ all_edges[(s,d)].append(layer_ID)
+ else:
+ all_edges[(s,d)] = [layer_ID]
+
+K = len(all_edges.keys())
+
+for k in all_edges.keys():
+ val = len(set(all_edges[(k)]))
+ s,d = k
+ print s, d, val
+
diff --git a/structure/metrics/intersect_layers.py b/structure/metrics/intersect_layers.py
new file mode 100644
index 0000000..a9da00c
--- /dev/null
+++ b/structure/metrics/intersect_layers.py
@@ -0,0 +1,47 @@
+####
+##
+## Take a certain number of networks as input and give as output the
+## corresponding intersection graph, i.e. the graph in which an edge
+## exists between i and j if and only if it exists in ALL the graphs
+## given as input
+##
+
+
+import sys
+
+
+if len(sys.argv)< 3:
+ print "Usage: %s <file1> <file2> [<file3>....]" % sys.argv[0]
+ sys.exit(1)
+
+
+graphs = {}
+
+num = 0
+
+for fname in sys.argv[1:]:
+
+ lines = open(fname).readlines()
+ graphs[num] = []
+ for l in lines:
+ s,d = [int(x) for x in l.strip(" \n").split(" ")][:2]
+ if s > d:
+ graphs[num].append((d,s))
+ else:
+ graphs[num].append((d,s))
+ num += 1
+
+#print graphs
+
+
+for edge in graphs[0]:
+ to_print = True
+ for i in range(1, num):
+ if edge not in graphs[i]:
+ to_print = False
+ break
+ if to_print:
+ s,d = edge
+ print s,d
+
+
diff --git a/structure/metrics/overlap_degree.py b/structure/metrics/overlap_degree.py
new file mode 100644
index 0000000..fa69434
--- /dev/null
+++ b/structure/metrics/overlap_degree.py
@@ -0,0 +1,47 @@
+####
+##
+## Compute the overlapping degree for each node and the corresponding
+## z-score
+##
+##
+
+import sys
+import numpy
+
+
+if len(sys.argv) < 2:
+ print "Usage: %s <layer1> <layer2> [<layer3>...]" % sys.argv[0]
+ sys.exit(1)
+
+
+nodes = {}
+
+for f in sys.argv[1:]:
+
+ lines = open(f).readlines()
+
+ for l in lines:
+ if l[0] == "#":
+ continue
+ s,d = [int(x) for x in l.strip(" \n").split(" ")]
+ if nodes.has_key(s):
+ nodes[s] +=1
+ else:
+ nodes[s] = 1
+ if nodes.has_key(d):
+ nodes[d] +=1
+ else:
+ nodes[d] = 1
+
+
+degrees = nodes.values()
+avg_deg = numpy.mean(degrees)
+std_deg = numpy.std(degrees)
+
+#print avg_deg, std_deg
+
+keys = nodes.keys()
+keys.sort()
+
+for n in keys:
+ print n, nodes[n], (nodes[n] - avg_deg)/std_deg
diff --git a/structure/reinforcement/reinforcement.py b/structure/reinforcement/reinforcement.py
new file mode 100644
index 0000000..5939e59
--- /dev/null
+++ b/structure/reinforcement/reinforcement.py
@@ -0,0 +1,61 @@
+import networkx as nx
+import sys
+
+
+if __name__ == "__main__":
+
+
+ if len(sys.argv) < 6:
+ print "Usage: %s <layer1> <layer2> <N bins> <min_value> <max_value>" % sys.argv[0]
+ sys.exit(1)
+
+
+ filename_a = sys.argv[1]
+ filename_b = sys.argv[2]
+
+ intervals=int(sys.argv[3])
+ minvalue=float(sys.argv[4])
+ maxvalue=float(sys.argv[5])
+
+
+ tot_a = []
+ pos_a = []
+ for t in range (intervals):
+ tot_a.append(0)
+ pos_a.append(0)
+
+
+
+
+ fa=open(filename_a, 'r')
+ Ga=nx.read_adjlist(fa)
+
+
+ fb=open(filename_b, 'r')
+ Gb=nx.read_weighted_edgelist(fb)
+
+
+
+
+ for u,v in Gb.edges():
+ Gbw=Gb[u][v]['weight']
+ for i in range (intervals):
+ a=minvalue+float(maxvalue-minvalue)*float(i)/intervals
+ b=minvalue+float(maxvalue-minvalue)*float(i+1)/intervals
+ if (Gbw>a and Gbw<b):
+ tot_a[i]+=1
+ break
+ if (Ga.has_edge(u,v)==True):
+ pos_a[i]+=1
+
+ freq_a=[]
+ for i in range (intervals):
+ if (tot_a[i]>0):
+ freq_a.append(float(pos_a[i])/tot_a[i])
+ else:
+ freq_a.append(0)
+ print "#bin_minvalue bin_maxvalue frequence"
+ for i in range (intervals):
+ a=minvalue+float(maxvalue-minvalue)*float(i)/intervals
+ b=minvalue+float(maxvalue-minvalue)*float(i+1)/intervals
+ print a, b, freq_a[i]