diff options
Diffstat (limited to 'models')
-rw-r--r-- | models/correlations/Makefile | 15 | ||||
-rw-r--r-- | models/correlations/fit_utils.c | 382 | ||||
-rw-r--r-- | models/correlations/fit_utils.h | 25 | ||||
-rw-r--r-- | models/correlations/rank_utils.c | 217 | ||||
-rw-r--r-- | models/correlations/rank_utils.h | 44 | ||||
-rw-r--r-- | models/correlations/tune_qnn_adaptive.c | 208 | ||||
-rw-r--r-- | models/correlations/tune_rho.c | 142 | ||||
-rw-r--r-- | models/growth/Makefile | 26 | ||||
-rw-r--r-- | models/growth/nibilab_linear_delay.c | 414 | ||||
-rw-r--r-- | models/growth/nibilab_linear_delay_mix.c | 457 | ||||
-rw-r--r-- | models/growth/nibilab_linear_delta.c | 403 | ||||
-rw-r--r-- | models/growth/nibilab_linear_random_times.c | 417 | ||||
-rw-r--r-- | models/growth/nibilab_nonlinear.c | 493 | ||||
-rw-r--r-- | models/growth/node_deg_over_time.py | 82 | ||||
-rw-r--r-- | models/nullmodels/model_MDM.py | 66 | ||||
-rw-r--r-- | models/nullmodels/model_MSM.py | 65 | ||||
-rw-r--r-- | models/nullmodels/model_hypergeometric.py | 56 | ||||
-rw-r--r-- | models/nullmodels/model_layer_growth.py | 121 |
18 files changed, 3633 insertions, 0 deletions
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 + |