diff options
Diffstat (limited to 'models/correlations')
| -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 |
7 files changed, 1033 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); + +} |