First commit of MAMMULT code

1 files changed, 208 insertions, 0 deletions

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);
+  
+}