ne_per_lib.c

/*===========================================================================
 NetEvo Foundation Library
 Copyright (C) 2009, 2010 Thomas E. Gorochowski <tgorochowski@me.com>
 Bristol Centre for Complexity Sciences, University of Bristol, Bristol, UK
 ---------------------------------------------------------------------------- 
 NetEvo is a computing framework designed to allow researchers to investigate 
 evolutionary aspects of dynamical complex networks. By providing tools to 
 easily integrate each of these factors in a coherent way, it is hoped a 
 greater understanding can be gained of key attributes and features displayed 
 by complex systems.
 
 NetEvo is open-source software released under the Open Source Initiative 
 (OSI) approved Non-Profit Open Software License ("Non-Profit OSL") 3.0. 
 Detailed information about this licence can be found in the COPYING file 
 included as part of the source distribution.
 
 This library 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.
 ============================================================================*/


#include "ne_per_lib.h"
#include <math.h>
#include <igraph.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>


ne_per_t * ne_per_top_eigenratio_alloc (void)
{
   /* No parameters */
   ne_per_t *per;
   
   if ((per = (ne_per_t *)malloc(sizeof(ne_per_t))) == NULL) {
      ne_error("ne_per_top_eigenratio_alloc: Could not allocate memory.");
      return NULL;
   }
   
   /* Configure the strucutre */
   per->type = NE_PER_TOP;
   per->fn = &ne_per_top_eigenratio;
   per->params = NULL;
   
   return per;
}


double ne_per_top_eigenratio (ne_sys_t *sys, ne_dyn_vec_t *dyn, ne_real_t *params)
{
   double r = 0.0;
   gsl_matrix* gLocal = NULL;
   gsl_eigen_symm_workspace* w = NULL;
   gsl_vector* v = NULL;
   
   /* Generate a GSL Laplacian Matrix as easier to work with */
   gLocal = ne_sys_laplacian_to_gsl_matrix_alloc(sys);
   
   /* Setup the eigensystem space and calculate the ordered eigenvalues */
   w = gsl_eigen_symm_alloc(gLocal->size1);
   v = gsl_vector_alloc(gLocal->size1);
   gsl_eigen_symm(gLocal, v, w);
   gsl_sort_vector(v);
   
   /* Calculate the eigen ratio lambda_N/lambda_2 */
   r = (double)gsl_vector_get(v, gLocal->size1-1)/(double)gsl_vector_get(v, 1);
   
   gsl_eigen_symm_free(w);
   gsl_vector_free(v);
   gsl_matrix_free(gLocal);
   
   return r;
}


ne_per_t * ne_per_dyn_orderparam_alloc (ne_real_t delta)
{
   /* 0 - delta */
   ne_real_t *newParams;
   ne_per_t *per;
   
   if ((per = (ne_per_t *)malloc(sizeof(ne_per_t))) == NULL) {
      ne_error("ne_per_top_eigenratio_alloc: Could not allocate memory.");
      return NULL;
   }
   
   if ((newParams = (ne_real_t *)malloc(sizeof(ne_real_t))) == NULL) {
      ne_error("ne_per_top_eigenratio_alloc: Could not allocate memory.");
      free(per);
      return NULL;
   }
   newParams[0] = delta;
   
   /* Configure the parameter strucutre */
   per->type = NE_PER_DYN;
   per->fn = &ne_per_dyn_orderparam;
   per->params = newParams;
   
   return per;
}


double ne_per_dyn_orderparam (ne_sys_t *sys, ne_dyn_vec_t *dyn, ne_real_t *params)
{
   /* 0 - delta (acceptable synchronisation error) */
   
   return ne_per_dyn_orderparam_calc_mu(sys, dyn, dyn->curSize-1, params[0]);
}


double ne_per_dyn_orderparam_integral (ne_sys_t *sys, ne_dyn_vec_t *dyn, ne_real_t *params)
{
   /* 0 - delta (acceptable synchronisation error) */
   
   int t;
   double y1, y2, tDiff, curSec, outSum;
   
   outSum = 0;
   
   for ( t=0; t<dyn->curSize; t++ ) {
      
      y1 = ne_per_dyn_orderparam_calc_mu(sys, dyn, t, params[0]);
      y2 = ne_per_dyn_orderparam_calc_mu(sys, dyn, t+1, params[0]);
      
      tDiff = dyn->data[t+1][0] - dyn->data[t][0];
      
      if ( y1 > y2 ) 
         curSec = (y1 * tDiff) - (((y1 - y2) * tDiff) / 2.0);
      else
         curSec = (y2 * tDiff) - (((y2 - y1) * tDiff) / 2.0);
      
      outSum += curSec;
   }
   
   return outSum;
}


double ne_per_dyn_orderparam_calc_mu (ne_sys_t *sys, ne_dyn_vec_t *dyn, int t, double delta)
{
   int i, j, e, N;
   double mu, d;
   
   mu = 0;
   N  = igraph_vcount(sys->graph);
   
   for ( i=0; i<N; i++ ) {
      for ( j=0; j<N; j++ ) {
         
         if ( i != j ) {
            
            /* Calculate the euclidean distance between the two nodes */
            d = 0;
            for ( e=0; e<sys->nodeStates; e++ ) {
               
               if ( isnan(dyn->data[t][(i*sys->nodeStates)+e+1]) || 
                   isnan(dyn->data[t][(j*sys->nodeStates)+e+1]) ) {
                  return 1.0;
               }
               
               d += pow(dyn->data[t][(i*sys->nodeStates)+e+1] - dyn->data[t][(j*sys->nodeStates)+e+1], 2.0);
            }
            d = sqrt(d);
            
            /* Calculate the heavyside function with delta error range */
            if ( d - delta >= 0 )
               mu += 100;
            else
               mu += 0;
         }
      }
   }
   
   /* Normalise */
   mu = mu * (1.0 / (N * (N - 1)));
   
   return mu;
}


ne_per_t * ne_per_dyn_bounded_alloc (ne_real_t tEnd)
{
   /* 0 - end period length to calculate bound using */
   ne_real_t *newParams;
   ne_per_t *per;
   
   if ((per = (ne_per_t *)malloc(sizeof(ne_per_t))) == NULL) {
      ne_error("ne_per_top_eigenratio_alloc: Could not allocate memory.");
      return NULL;
   }
   
   if ((newParams = (ne_real_t *)malloc(sizeof(ne_real_t))) == NULL) {
      ne_error("ne_per_top_eigenratio_alloc: Could not allocate memory.");
      return NULL;
   }
   newParams[0] = tEnd;
   
   /* Configure the parameter strucutre */
   per->type = NE_PER_DYN;
   per->fn = &ne_per_dyn_bounded;
   per->params = newParams;
   
   return per;
}


double ne_per_dyn_bounded (ne_sys_t *sys, ne_dyn_vec_t *dyn, ne_real_t *params)
{
   long int t, i, j, e;
   double d = 0, curBound = -1, tStart, indStart;
   
   curBound = 0.0;
   
   /* Find the start time in the dynamics output */
   tStart = dyn->data[dyn->curSize-1][0] - params[0];
   
   /* Search for the index at which this time corrisponds */
   indStart = -1;
   for (i=0; i<dyn->curSize; i++) {
      
      if (dyn->data[i][0] > tStart) {
         indStart = i-1;
         break;
      }
   }
   
   /* Check that an index was found */
   if (indStart == -1) {
      
      /* Error occured */
      ne_error("ne_per_dyn_bounded - Could not find start index.\n");
      return NE_HUGE_NUMBER;
   }
   
   /* For each time point calculate the largest bound */
   for (t=indStart; t<dyn->curSize; t++) {
      
      
      /* For each pair of vertices (i, j) where (i != j) */
      for (i=0; i<igraph_vcount(sys->graph); i++) {
         for (j=0; j<igraph_vcount(sys->graph); j++) {
            
            /* Check that a distance will exist */
            if (i != j) {
               
               /* Calculate the euclidean distance between the two nodes */
               d = 0;
               for (e=0; e<sys->nodeStates; e++) {
                  
                  if (isnan(dyn->data[t][(i*sys->nodeStates)+e+1]) || 
                      isnan(dyn->data[t][(j*sys->nodeStates)+e+1]) ) {
                     return NE_HUGE_NUMBER;
                  }
                  
                  d += pow(dyn->data[t][(i*sys->nodeStates)+e+1] - dyn->data[t][(j*sys->nodeStates)+e+1], 2.0);
               }
               d = sqrt(d);
               
               if (d > curBound) {
                  curBound = d;
               }
            }
         }
      }
   }
   
   return curBound;
}