path: root/src/fitmle/fitmle.c

/**
 *  This program is free software: you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see
 *  <http://www.gnu.org/licenses/>.
 *
 *  (c) Vincenzo Nicosia 2009-2017 -- <v.nicosia@qmul.ac.uk>
 * 
 *  This file is part of NetBunch, a package for complex network
 *  analysis and modelling. For more information please visit:
 *
 *             http://www.complex-networks.net/
 *
 *  If you use this software, please add a reference to 
 *
 *               V. Latora, V. Nicosia, G. Russo             
 *   "Complex Networks: Principles, Methods and Applications"
 *              Cambridge University Press (2017) 
 *                ISBN: 9781107103184
 *
 ***********************************************************************
 *
 *  This program fits a set of values provided as input with a
 *  power-law function using the Maximum-Likelihood Estimator (MLE).
 *
 *  References: 
 * 
 * [1] A. Clauset, C. R. Shalizi, and M. E. J. Newman. "Power-law
 *     distributions in empirical data". SIAM Rev. 51, (2007),
 *     661-703.  
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>

#include "utils.h"

void usage(char *argv[]){

  printf("********************************************************************\n"
         "**                                                                **\n"
         "**                        -*-   fitmle   -*-                      **\n"
         "**                                                                **\n"
         "**    Fit a set of values with a power-law function using the     **\n"
         "**    Maximum-Likelihood Estimator (MLE). The program implements  **\n"
         "**    the MLE for both continuous and discrete values, and        **\n"
         "**    selects the appropriate one automatically.                  **\n"
         "**                                                                **\n"
         "**    The input file 'data_in' contains one value on each row.    **\n"
         "**    If 'data_in' is '-' (dash), read the values from the        **\n"
         "**    standard input (STDIN).                                     **\n"
         "**                                                                **\n"
         "**    The program prints on output a single row, in the format:   **\n"
         "**                                                                **\n"
         "**       gamma x_min ks                                           **\n"
         "**                                                                **\n"
         "**    where 'gamma' is the exponent of the best fit, 'x_min' is   **\n"
         "**    the smallest of the input values after which the power-law  **\n"
         "**    hypothesis hold, and 'ks' is the value of the KS statistics **\n"
         "**    for the best fit.                                           **\n"
         "**                                                                **\n"
         "**    The second (optional) parameter 'tol' sets the tolerance    **\n"
         "**    on the acceptable statistical error of the fit (set to      **\n"
         "**    0.1 if no value is provided).                               **\n"
         "**                                                                **\n"
         "**    If the third parameter is 'TEST', the program uses the      **\n"
         "**    the Kolmogorov-Smirnov test on a series of bootstrapped     **\n"
         "**    power-law sequences to estimate the p-value of the fit, and **\n"
         "**    the output is in the format:                                **\n"
         "**                                                                **\n"
         "**       gamma x_min ks p_value                                   **\n"
         "**                                                                **\n"
         "**    Higher values of 'p_value' correspond to better fits.       **\n"
         "**                                                                **\n"
         "**    If 'num_test' is provided, use that number of bootstrapped  **\n"
         "**    sequences to compute the p-value. Otherwise, use 100        **\n"
         "**    sequences.                                                  **\n"
         "**                                                                **\n"
         "********************************************************************\n"
         " This is Free Software - You can use and distribute it under \n"
         " the terms of the GNU General Public License, version 3 or later\n\n"
         " (c) Vincenzo Nicosia 2010-2017 (v.nicosia@qmul.ac.uk)\n\n"
         "********************************************************************\n\n"
         );
  printf("Usage: %s <data_in> [<tol> [ TEST [<num_test>]]]\n", argv[0]);
}


void load_data(char *fname, double **data, unsigned int *N, char sort){
  
  FILE *fin;
  char error_str[256];
  char buff[256];
  char *ptr;
  double x;
  unsigned int size;
  
  size = 10;

  *data = malloc(size * sizeof(double));

  *N=0;
  if (!strcmp(fname, "-")){
    fin = stdin;
  }
  else{
    fin = fopen(fname, "r");
  }
  if (!fin){
    sprintf(error_str, "Error opening file %s", fname);
    perror(error_str);
    exit(3);
  }
  
  while(fgets(buff, 256, fin)){
    ptr = strtok(buff, " ");
    if (ptr[0] == '#')
      continue;
    x = atof(ptr);
    if (*N == size){
      size += 10;
      *data = realloc(*data, size*sizeof(double));
    }
    (*data)[*N] = x;
    *N += 1;
  }
  *data = realloc(*data, (*N)*sizeof(double));
  if (sort){
    qsort(*data, *N, sizeof(double), compare_double);
  }
  fclose(fin);
  return;
}

/** 
 *  find the first position at which the element x appears in the
 *  "data" array using binary search (the array is sorted in
 *  increasing order of x)
 */
int find_pos_x(double x,double *data, unsigned int N){
  int H, L, M;

  L = 0;
  H = N-1;
  
  while(L<=H){
    M = (H + L)/2;
    if (x == data[M]){
      while(M >=0 && x == data[M]) M--;
      if (M==-1){
        return 0;
      }
      return M+1; /* we return the index of the first appearance of element x*/
    }
    else if (x < data[M]){
      H = M-1;
    }
    else if (x > data[M]){
      L = M+1;
    }
  }
  return -1; /* the value is not present */
}

double fit_alpha(double *data, unsigned int N, double xmin, unsigned int idx){
  
  double alpha = 0;
  int i;
  
  for(i = idx; i < N; i++){
    alpha += log(data[i]*1.0/(xmin-0.5));
  }
  alpha = 1 + (N - idx) * 1.0 / alpha;
  return alpha;
}


double fit_alpha_continuous(double *data, unsigned int N, double xmin, unsigned int idx){
  
  double alpha = 0;
  int i;
  
  for(i = idx; i < N; i++){
    alpha += log(data[i]*1.0/(xmin));
  }
  alpha = 1 + (N - idx) * 1.0 / alpha;
  return alpha;
}


/*
 *
 * Kolmogorov-Smirnov test
 *
 */

double KS (double *data, unsigned int N, double alpha, int idx){
  
  double c_data, c_theo, val;
  double xmin, x;
  double dist, max_dist = -1;
  int i;

  c_data = c_theo = 0.0;
  xmin = data[idx];

  for (i=idx; i<N;){
    x = data[i];
    while(i<N && data[i] == x){
      val = xmin * 1.0 / data[i];
      c_theo = 1.0 -  pow(val, alpha-1.0);
      dist = fabs(c_theo - c_data);

      if (dist > max_dist){
        max_dist = dist;
      }
      i++;
    }
    c_data = 1.0 * (i- idx)/(N-idx);
  }
  return max_dist;
}

void dump_data_double(double *v, unsigned int N){
  int i;

  if (N < 1){
    return;
  }
  printf("%g",v[0]);
  for(i=1; i<N; i++){
    printf(" %g",v[i]);
  }
  printf("\n");
}


/* This is the same as best_fit, but for continuous-valued time-series */
void best_fit_continuous(double *data, unsigned int N, double *alpha, double *xmin, double tol){
  
  double min_x, max_x, x;
  double cur_alpha, best_alpha, best_x, ks_test, best_ks;
  int idx,i;
  

 

  qsort(data, N, sizeof(double), compare_double);

  min_x = data[0];
  best_x = data[0];
  max_x = data[N-1];
  best_ks = 10000000;
  cur_alpha = 0.0;
  best_alpha = 0.0;
  idx = 0;
  for(x=min_x, i=0; x<=max_x && (cur_alpha)/sqrt(1.0 * (N-idx))< tol; ){
    idx = find_pos_x(x, data, N);
    if (idx <0){ 
      i ++;
      continue;
    }
    cur_alpha = fit_alpha_continuous(data, N, x, idx);
    ks_test = KS(data, N, cur_alpha, idx);
    if(ks_test < best_ks){
      best_ks = ks_test;
      best_alpha=cur_alpha;
      best_x = data[idx];
    }
    while(data[i] == x && i <N){
      i ++;
    }
    x = data[i];
  }

  *alpha = best_alpha;
  *xmin = best_x;
  return;
}

/* Fit a discrete power-law */

void best_fit(double *data, unsigned int N, double *alpha, double *xmin, double tol){
  
  double min_x, max_x, x;
  double cur_alpha, best_alpha, best_x, ks_test, best_ks;
  int idx, i;
  

  
  qsort(data, N, sizeof(double), compare_double);
  min_x = data[0];
  best_x = data[0];
  max_x = data[N-1];
  best_ks = 10000000;
  cur_alpha = 0.0;
  best_alpha = 0.0;
  idx = 0;
  for(x=min_x, i = 0; x<=max_x && (cur_alpha)/sqrt(1.0 * (N-idx))< tol;){
    idx = find_pos_x(x, data, N);

    if (idx<0) {
      i ++;
      continue;
    }
    cur_alpha = fit_alpha(data, N, x, idx);

    ks_test = KS(data, N, cur_alpha, idx);
    if(ks_test < best_ks){
      best_ks = ks_test;
      best_alpha=cur_alpha;
      best_x = data[idx];
    }
    while(data[i] == x && i<N){
      i ++;
    }
    x = data[i];
  }

  *alpha = best_alpha;
  *xmin = best_x;
  return;
}

void sample_powerlaw(double *data, unsigned int N, double alpha, double xmin,double *v){
  int i, last_idx, r;
  double val;
  
  i=0;
  last_idx = 0;
  while(data[i] < xmin){
    last_idx ++;
    i++;
  }
  
  i=0;
  while(i <last_idx){
    r = rand() % (last_idx + 1);
    v[i] = data[r];
    i++;
  }
  
  for(; i<N; i++){
    val = 1.0 * rand()/RAND_MAX;
    v[i] = floor(xmin * pow(1.0-val, -1.0/(alpha-1)));

  }
  
}


double test_powerlaw(double *data, unsigned int N, double alpha, double xmin, 
                     double ks, int num_test, 
                     void (*f)(double *, unsigned int , double *, double *, double)){
  
  int num = 0, i, idx;
  double *v, cur_alpha, cur_xmin, cur_ks;


  v = malloc(N * sizeof(double));

  for(i=0; i<num_test; i++){
    sample_powerlaw(data, N, alpha, xmin, v);
    qsort(v, N, sizeof(double), compare_double);
    f(v, N, &cur_alpha, &cur_xmin, 0.1); 
    idx = find_pos_x(cur_xmin, v, N);
    cur_ks = KS(v, N, cur_alpha, idx);
    if (cur_ks > ks){
      num += 1;
    }
  }
  free(v);
  return 1.0 * num / num_test;
}

/* We assume that the data set has already been sorted in ascending
   order */
int is_continuous(double *data, int N){
  int i;

  for (i=0; i<N; i++){
    if (data[i] - (double)(int)data[i] != 0.0)
      return 1;
  }
  /* return 1 only if we need a real-valued fit....*/
  return 0;
}

/* we assume that the data set has already been sorted in ascending
   order */
double renormalise(double *data, unsigned int N){
  
  int i;
  double scaling = 1.0;
  
  if (data[0] < 1.0 ){
    scaling = 1.0/data[0];
    for (i=0; i<N; i++){
      data[i] *= scaling;
    }
  }
  return scaling;
}


int main(int argc, char *argv[]){
  
  double *data=NULL;
  unsigned int N, i;
  double alpha, xmin, ks, tol, p_value, scaling_fact;
  char test = 0;
  int num_test;
  
  void  (*fit_func)(double *, unsigned int , double *, double *, double);
  

  if (argc < 2){
    usage(argv);
    exit(1);
  }
  
  if (argc > 2)
    tol = atof(argv[2]);
  else
    tol = 0.1;

  if (argc > 3 && !my_strcasecmp(argv[3], "TEST")){
    test = 1;
  }
  
  if(argc > 4){
    num_test = atoi(argv[4]);
  }
  else
    num_test = 100;
  
  srand(time(NULL));
  

  
  load_data(argv[1], &data, &N, 1);

  if(is_continuous(data, N)){
    fprintf(stderr, "Using continuous fit\n");
    fit_func = best_fit_continuous;
  }
  else{
    fprintf(stderr, "Using discrete fit\n");
    fit_func = best_fit;
  }

  scaling_fact = 1.0;

  fit_func(data, N, &alpha, &xmin, tol);
  i = find_pos_x(xmin, data, N);
  ks = KS(data, N, alpha, i);
  if (test){
    p_value = test_powerlaw(data, N, alpha, xmin, ks, num_test, fit_func);
    printf("%g %g %g %g\n", alpha,  xmin/scaling_fact, ks, p_value);
  }
  else{
    printf("%g %g %g\n", alpha, xmin / scaling_fact, ks);
  }
  free(data);
}