ne_sim_gslode.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_sim_gslode.h"
#include "ne_sim.h"
#include <gsl/gsl_rng.h>
#include <string.h>


ne_sim_gslode_params_t * ne_sim_gslode_params_alloc (const gsl_odeiv_step_type *step, 
                                   double eps_abs, double eps_rel, double length, 
                                   ne_bool_t fixedStep, double initStep, double minStep)
{  
   ne_sim_gslode_params_t *params;
   
   if ((params = (ne_sim_gslode_params_t *)malloc(sizeof(ne_sim_gslode_params_t))) == NULL) {
      ne_error("ne_sim_gslode_alloc: Could not allocate memory for parameters");
      return NULL;
   }
   
   params->length    = length;
   params->step      = step;
   params->initStep  = initStep;
   params->minStep   = minStep;
   params->fixedStep = fixedStep;
   params->eps_abs   = eps_abs;
   params->eps_rel   = eps_rel;
   
   return params;
}


void ne_sim_gslode_params_free (ne_sim_gslode_params_t *params) 
{
   /* Free the full structure */
   free(params);
}


ne_dyn_vec_t * ne_sim_gslode (ne_sys_t *sys, void *p, double *iy, FILE *file)
{  
   double     t, t1, h, hMin, *y, *yErr, *yCopy, *dydtIn, *dydtOut;
   int        fixedStep, status, i, stepStatus, dimension;
   ne_dyn_vec_t *vecOut;
   ne_sim_gslode_params_t *params;
   gsl_odeiv_system odeSystem;
   gsl_odeiv_step *step;
   gsl_odeiv_control *control;
   gsl_odeiv_evolve *evolve;
   
   params = (ne_sim_gslode_params_t *)p;
   dimension = (ne_sys_nodes(sys) * ne_sys_node_states(sys)) +  
            (ne_sys_edges(sys) * ne_sys_edge_states(sys));
   
   /* Build the ODE system to run the simulation */
   odeSystem.function  = &ne_sim_gslode_sys_fn;
   odeSystem.jacobian  = NULL;
   odeSystem.dimension = dimension;
   odeSystem.params    = (void *)sys;
   
   /* Define the parameters required by NetEvo when simulating */
   step    = gsl_odeiv_step_alloc(params->step, dimension);
   control = gsl_odeiv_control_y_new(params->eps_abs, params->eps_rel);
   evolve  = gsl_odeiv_evolve_alloc(dimension);

   /* Allocate dynamic vector to hold output */
   vecOut = ne_dyn_vec_alloc(dimension + 1);
   
   /* Setup internal variables and allocate required memory */
   t         = 0.0;
   t1        = params->length;
   h         = params->initStep;
   hMin      = params->minStep;
   fixedStep = params->fixedStep;
   y = (double*)malloc(sizeof(double)*(dimension));
   
   /* Copy accross the initial starting point */
   memcpy(y, iy, sizeof(double)*dimension);
   yCopy = (double*)malloc((sizeof(double)*dimension)+1);
   memcpy(&yCopy[1], y, sizeof(double)*dimension);
   yCopy[0] = t;
   ne_dyn_vec_append_item(vecOut, yCopy);

   if ( fixedStep == TRUE ) {
      
      /* Allocate memory to hold the error */
      yErr = (double*)malloc((sizeof(double)*dimension));
      dydtIn = (double*)malloc((sizeof(double)*dimension));
      dydtOut = (double*)malloc((sizeof(double)*dimension));
      
      /* Calculate the initial dydt value */
      GSL_ODEIV_FN_EVAL(&odeSystem, t, y, dydtIn);
      
      while ( t < t1 ) {
         
         status = gsl_odeiv_step_apply(step, t, h, y, yErr, dydtIn, dydtOut, &odeSystem);
         
         /* Update the current time to the new timepoint */
         t += h;
         
         /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
         /* TODO: Check errors and return error if found to be outside range */
         /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
         
         if (status != GSL_SUCCESS) {
            break;
         }
         
         /* Copy output of previous step to input of next */
         for ( i=0; i<dimension; i++ ) {
            dydtIn[i] = dydtOut[i];
         }
         
         /* Make a copy of the results and append to the output vector */
         yCopy = (double*)malloc((sizeof(double)*dimension)+1);
         memcpy(&yCopy[1], y, sizeof(double)*dimension);
         yCopy[0] = t;
         ne_dyn_vec_append_item(vecOut, yCopy);
      }
   }
   else {
      
      if ( hMin > 0 ) {
         /* Allocate memory to hold the error */
         yErr = (double*)malloc((sizeof(double)*dimension));
         dydtIn = (double*)malloc((sizeof(double)*dimension));
         dydtOut = (double*)malloc((sizeof(double)*dimension));
         
         /* Calculate the initial dydt value */
         GSL_ODEIV_FN_EVAL(&odeSystem, t, y, dydtIn);
         
         /* No error for the first point */
         for ( i=0; i<dimension; i++ ) {
            yErr[i] = 0;
         }
      }

      while ( t < t1 ) {
         
         if ( hMin > 0 ) {
               
            /* Must maintain a minimum h value */
            stepStatus = gsl_odeiv_control_hadjust(control, step, y, yErr, dydtIn, &h);
            
            if ( stepStatus == GSL_ODEIV_HADJ_DEC && h < hMin ) {
               h = hMin;
               ne_error("ne_dyn_ode_solver_evolve - Estimated step size is smaller than given minimum. Error bounds will be exceeded.\n");
            }
            
            status = gsl_odeiv_step_apply(step, t, h, y, yErr, dydtIn, dydtOut, &odeSystem);
            
            /* Update the current time */
            t += h;
            
            if (status != GSL_SUCCESS) {
               break;
            }
            
            for ( i=0; i<dimension; i++ ) {
               dydtIn[i] = dydtOut[i];
            }
            
            yCopy = (double*)malloc((sizeof(double)*dimension)+1);
            memcpy(&yCopy[1], y, sizeof(double)*dimension);
            yCopy[0] = t;
            ne_dyn_vec_append_item(vecOut, yCopy);
            
            /* Ensure we do not exceed t1 */
            if ( t > t1 ) {
               t = t1;
            }           
         }
         else {
            
            /* Allow GSL to evolve the system to maintain bounded errors */
            status = gsl_odeiv_evolve_apply(evolve, control, step, &odeSystem, &t, t1, &h, y);
         }
            
         if (status != GSL_SUCCESS) {
            break;
         }
         
         yCopy = (double*)malloc((sizeof(double)*dimension)+1);
         memcpy(&yCopy[1], y, sizeof(double)*dimension);
         yCopy[0] = t;
         ne_dyn_vec_append_item(vecOut, yCopy);
      }
   }
   
   /* Check that simulation was successful */
   if (status != GSL_SUCCESS) {
      ne_dyn_vec_free(vecOut);
      free(y);
      return NULL;
   }
   
   if ( file != NULL ) {
      ne_sim_fprint_dyn_vec(vecOut, file);
   }  
      
   /* Free used memory and return dynamical output */
   free(y);
   gsl_odeiv_step_free(step);
   gsl_odeiv_control_free(control);
   gsl_odeiv_evolve_free(evolve);
   return vecOut;
}


int ne_sim_gslode_sys_fn (double t, const double y[], double dydt[], void *params)
{  
   int nodes, edges, nStates, eStates, i, j;
   ne_dyn_node_t *nodeDyn;
   ne_dyn_edge_t *edgeDyn;
   ne_sys_t *sys;
   
   sys = (ne_sys_t *)params;
   
   nodes = ne_sys_nodes(sys);
   edges = ne_sys_edges(sys);
   nStates = ne_sys_node_states(sys);
   eStates = ne_sys_edge_states(sys);
   
   /* Calculate the dynamics for each node */
   for (i=0; i<nodes; i++) {
      nodeDyn = ne_sys_dyn_node(sys, i);
      (nodeDyn->fn)(i, (void *)sys, t, y, dydt, ne_sys_dyn_node_params(sys,i));
   }
   
   /* Calculate the dynamics for each edge */
   if (igraph_is_directed(sys->graph)) {
      for (j=0; j<edges; j++) {
         edgeDyn = ne_sys_dyn_eid(sys, j);
         (edgeDyn->fn)(j, (int)IGRAPH_FROM(sys->graph, (igraph_integer_t)j), (int)IGRAPH_TO(sys->graph, (igraph_integer_t)j), 
                    (void *)sys, t, y, dydt, ne_sys_dyn_eid_params(sys,j));
      }
   }
   else {
      for (j=0; j<edges; j++) {
         edgeDyn = ne_sys_dyn_eid(sys, j);
         (edgeDyn->fn)(j, (int)IGRAPH_FROM(sys->graph, (igraph_integer_t)j), (int)IGRAPH_TO(sys->graph, (igraph_integer_t)j), 
                    (void *)sys, t, y, dydt, ne_sys_dyn_eid_params(sys,j));
         (edgeDyn->fn)(j, (int)IGRAPH_TO(sys->graph, (igraph_integer_t)j), (int)IGRAPH_FROM(sys->graph, (igraph_integer_t)j), 
                    (void *)sys, t, y, dydt, ne_sys_dyn_eid_params(sys,j));
      }
   }
   
   return GSL_SUCCESS;
}