ne_sup_sa.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_sup_sa.h"
#include "ne_evo.h"
#include <time.h>
#include <gsl/gsl_rng.h>
#include <math.h>

/* Shared random number generator */
extern gsl_rng* neRng;


ne_sup_sa_params_t * ne_sup_sa_params_alloc (int initialTrials, double tempReduce, int mainTrails,
                                    int acceptTrials, int acceptRunsNoChange, double minTemp,
                                    int maxIterations, double (*initialTempFn) (double),
                                    ne_per_t *QFn, ne_mut_t *mutationFn, ne_sim_fn_t *simFn,
                                    void *simParams, ne_dyn_vec_t *initCond, 
                                    ne_sup_sa_analysis_fn_t *analysisFn, FILE *evoFile)
{
   ne_sup_sa_params_t *params;
   
   if ((params = (ne_sup_sa_params_t *)malloc(sizeof(ne_sup_sa_params_t))) == NULL) {
      ne_error("ne_sup_sa_alloc: Could not allocate memory for parameters.");
      return NULL;
   }
   
   params->initialTrials = initialTrials;
   params->tempReduce = tempReduce;
   params->mainTrials = mainTrails;
   params->acceptTrials = acceptTrials;
   params->acceptRunsNoChange = acceptRunsNoChange;
   params->minTemp = minTemp;
   params->maxIterations = maxIterations;
   params->initialTempFn = initialTempFn;
   params->QFn = QFn;
   params->mutationFn = mutationFn;
   params->simFn = simFn;
   params->simParams = simParams;
   params->initCond = initCond;
   params->analysisFn = analysisFn;
   params->evoFile = evoFile;
   
   return params;
}


void ne_sup_sa_params_free (ne_sup_sa_params_t *params)
{
   free(params);
}


ne_sys_t * ne_sup_sa_run (ne_sys_t *sys, ne_sup_sa_params_t *params)
{  
   int i, accepts, noChange, iteration;
   double T, dQmax;
   ne_sup_sa_trail_results_t *results = NULL;
   double QChosen;
   
   /* Initialise variables */
   iteration = 1;
   
   /* Create structure used to return the results of each trial */
   results = (ne_sup_sa_trail_results_t *)malloc(sizeof(ne_sup_sa_trail_results_t));
   results->curDyn = NULL;
   results->Q1 = -1;
   
   /* First trials at initial T of infinity (or very large in our case) */
   T = 10000000000000.0;
   dQmax = 0.0;
   
   if (params->analysisFn != NULL) {
      (params->analysisFn)(sys, NULL, iteration, T, -1, TRUE);
   }
   
   if (params->evoFile != NULL) {
      ne_evo_initial(sys, params->evoFile);
   }
   
   /* ===================== FIND INITIAL TEMPERATURE ===================== */
   
   /* Carry out a number of initial trails to estimate starting temperature */
   for ( i=0; i<params->initialTrials; i++ ) {
      
      /* Try a trail at infinite temperature */
      sys = ne_sup_sa_trail(sys, T, results, params);
      
      /* Check to see if accepted and output result */
      if ( results->a == 1 ) {
         if (params->analysisFn != NULL) {
            (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q2, FALSE);
         }
      }
      else {
         if (params->analysisFn != NULL) {
            (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q1, FALSE);
         }
      }
      
      /* Find the worse result for the performance measure to initialise inital temp */
      dQmax = fmax(results->dQ, dQmax);
      
      /* See if we have already found an optimum and stop if so */
      if ((results->a == 1 && results->Q2 == 0.0) ||
          (results->a == 0 && results->Q1 == 0.0)) {
         return sys;
      }
      
      iteration++;
   }

   /* Update tempertaure based on these trials */
   if (dQmax != 0.0) {
      T = (params->initialTempFn)(dQmax);
   }
   else {
      /* No difference could be found so pick initial temp based on performance measure */
      T = results->Q1 * 0.2;
   }
   
   /* ===================== SIMULATED ANNEALING PROCESS ===================== */
   
   /* Flag signifying the number of runs with no change */
   noChange = 0;

   /* Check that the temperature has not already been set to 0 */
   if ( T != 0.0 ) {
      
      /* Continue to lower temperature by tempReducePer% each run until no changes in acceptRunsNoChange runs */
      do {
         
         /* Run for mainTrials more trials or acceptTrials accepting trials */
         accepts = 0;
         for ( i=0; i<params->mainTrials; i++ ) {
         
            /* Check if we have reached the maximum number of iterations */
            if (iteration > params->maxIterations) {
               break;
            }           
            
            /* Run a trail */
            sys = ne_sup_sa_trail(sys, T, results, params);
            
            /* Check to see if accepted and output result */
            if ( results->a == 1 ){
               if (params->analysisFn != NULL) {
                  (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q2, FALSE);
               }
               QChosen = results->Q2;
            }
            else{
               QChosen = results->Q1;
               if (params->analysisFn != NULL) {
                  (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q1, FALSE);
               }
            }
            
            iteration++;

            /* Update accepting counters */
            if ( (int)results->a == 1 ){
               accepts++;
               if ( results->Q2 == 0.0 )
                  break;
            }
            else{
               if ( results->Q1 == 0.0 )
                  break;
            }
            if ( accepts >= params->acceptTrials ) 
               break;
         }
         
         /* Update the counter to check if no changes are made at lower temps */
         if ( accepts == 0 )
            noChange++;
         else
            noChange = 0;
         
         /* Reduce temperature */
         T = T*params->tempReduce;
         
      } while (noChange <= params->acceptRunsNoChange && 
               T > params->minTemp && 
               iteration <= params->maxIterations && 
               QChosen != 0.0);
   }
   
   
   /* Output the final results */
   if ( results->a == 1 ){
      if (params->analysisFn != NULL) {
         (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q2, TRUE);
      }
   }
   else{
      if (params->analysisFn != NULL) {
         (params->analysisFn)(sys, results->curDyn, iteration, T, results->Q1, TRUE);
      }
   }
   
   /* Free used memory */
   ne_dyn_vec_free(results->curDyn);
   free(results);
   
   /* Return the new ne_sys */
   return sys;
}


ne_sys_t * ne_sup_sa_trail (ne_sys_t *sys, double T, ne_sup_sa_trail_results_t *results, 
                 ne_sup_sa_params_t *params)
{
   int a;
   double dQ;
   igraph_bool_t res;
   ne_dyn_vec_t *dynOut2 = NULL;
   ne_sys_t *S2 = NULL;
   char strBuffer[1024];
   
   strBuffer[0] = '\0';
   
   /* Set inital dQ value incase none is found */
   dQ = 0.0;

   /* Make a copy of the system */
   S2 = ne_sys_clone(sys);

   /* Call mutational function to generate new topology */
   (params->mutationFn->fn)(S2, params->mutationFn->params, strBuffer);

   /* Do not accept changes by default */
   a = 0;
   
   /* Find the Q1 */
   /* Depending on the type of performance measure simulation of dynamics may be required */
   switch (params->QFn->type) {
         
      case NE_PER_TOP:
         if ( results->Q1 == -1 ) {
            results->Q1 = (params->QFn->fn)(sys, NULL, params->QFn->params);
         }
         else {
            if (results->a == 1) {
               results->Q1 = results->Q2;
            }
         }
         break;
         
      case NE_PER_DYN:
      case NE_PER_TOP_AND_DYN:
         if ( results->Q1 == -1 ) {
            results->Q1 = ne_sup_sa_simulate(sys, params);
         }
         else{
            if ( results->a == 1 ) {
               results->Q1 = results->Q2;
            }
         }
         break;
         
      default:
         results->Q1 = -1;
         break;
   }
   
   /* Check that connected and if not reject trial */
   igraph_is_connected(S2->graph, &res, IGRAPH_WEAK);
   if (res == TRUE) {
      
      /* Depending on the type of performance measure simulation of dynamics may be required */
      switch (params->QFn->type) {
            
         case NE_PER_TOP:
            results->Q2 = (params->QFn->fn)(S2, NULL, params->QFn->params);
            dQ = results->Q2 - results->Q1;
            break;
            
         case NE_PER_DYN:
         case NE_PER_TOP_AND_DYN:
            results->Q2 = ne_sup_sa_simulate(S2, params);
            dQ = results->Q2 - results->Q1;     
            break;
            
         default:
            dQ = 0;
            break;
      }
      
      if (dQ < 0.0) {
         
         /* Accept change */
         a = 1;
      }
      else {
         
         /* Check for divide by zero */
         if (T != 0.0 && T > 0.0) {
            
            /* Calculate probability of still accepting */
            /* There are several possible accepting forms we use standard */
            if (gsl_rng_uniform(neRng) <= exp(-dQ/T)) {
               /* Make change even though worse cost */
               a = 1;
            }
         }
      }
   }

   /* Set results to be used by caller */
   results->dQ = dQ;
   results->a = a;
   
   /* Make updates based on outcome of change */
   if ( a == 1 ) {

      /* Update the current dynamics to the new data */
      ne_dyn_vec_free(results->curDyn);
      results->curDyn = dynOut2;
      
      /* Output the evolutionary steps */
      if (params->evoFile != NULL){
         fprintf(params->evoFile, "%s%s", strBuffer, ne_evo_step(NE_STEP_EVO));
         fflush(params->evoFile);
      }
      
      /* Update the master network as change was accepted */
      ne_sys_free(sys);
      return S2;
   }
   else{
      
      /* Free memory used for altered dynamics */
      ne_dyn_vec_free(dynOut2);
      
      /* Change was not successful */
      ne_sys_free(S2);
      return sys;
   }
}


double ne_sup_sa_simulate (ne_sys_t *sys, ne_sup_sa_params_t *params) {
   
   int i;
   ne_dyn_vec_t **dynOuts;
   double *perfVals;
   double perfSum = 0;
   
//#ifdef BENCHMARK
   clock_t ctime_1, ctime_2;
   ctime_1 = clock();
//#endif
   
   perfSum = 0;
   dynOuts = (ne_dyn_vec_t **)malloc(sizeof(ne_dyn_vec_t *)*params->initCond->curSize);
   perfVals = (double *)malloc(sizeof(double)*params->initCond->curSize);
   
   /* Run each simulation using OpenMP */
   #pragma omp parallel for
   for ( i=0; i<params->initCond->curSize; i++ ) {
      
      /* Simulate the system */
      dynOuts[i] = (params->simFn)(sys, params->simParams, params->initCond->data[i], NULL);
      
      /* Calculate the performance measure */
      perfVals[i] += (params->QFn->fn)(sys, dynOuts[i], params->QFn->params);
      
      /* Free simulation memory */
      ne_dyn_vec_free(dynOuts[i]);
   }
   
   for (i=0; i<params->initCond->curSize; i++ ) {
      perfSum += perfVals[i];
   }
   
   free(dynOuts);
   free(perfVals);
   
//#ifdef BENCHMARK
   ctime_2 = clock();
   ne_log("BENCHMARK: ne_sup_sa_simulate - CPU time to simulate dynamics: %f seconds\n", 
         (double)(ctime_2 - ctime_1) / (double)CLOCKS_PER_SEC);
//#endif
   
   /* Return the average performance measure */
   return perfSum / params->initCond->curSize;
}


double ne_sup_sa_int_temp_basic (double qMax) 
{
   return 4*qMax;
}