diff options
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] |