ne_sys.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_sys.h"
#include <string.h>


ne_sys_t * ne_sys_alloc (ne_bool_t directed, ne_int_t nodes, ne_int_t nodeStates, 
                   ne_int_t edgeStates, ne_bool_t dynNodeVary, ne_bool_t dynEdgeVary, 
                   ne_bool_t dynNodeParamsVary, ne_bool_t dynEdgeParamsVary)
{  
   ne_sys_t *sys;
   igraph_t G;
   igraph_empty(&G, (igraph_integer_t)nodes, (igraph_bool_t)directed);
   sys = ne_sys_alloc_from_igraph(&G, nodeStates, edgeStates, dynNodeVary, dynEdgeVary, 
                           dynNodeParamsVary, dynEdgeParamsVary);
   igraph_destroy(&G);
   return sys;
}


ne_sys_t * ne_sys_alloc_from_igraph (igraph_t *G, ne_int_t nodeStates, ne_int_t edgeStates, 
                            ne_bool_t dynNodeVary, ne_bool_t dynEdgeVary, 
                            ne_bool_t dynNodeParamsVary, ne_bool_t dynEdgeParamsVary)
{  
   ne_sys_t *sys = NULL;
   
   if ((sys = (ne_sys_t *)malloc(sizeof(ne_sys_t))) == NULL) {
      ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.");
      return NULL;
   }

   sys->graph = (igraph_t *)malloc(sizeof(igraph_t));
   if (sys->graph == NULL) {
      ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.graph.");
      free(sys);
      return NULL;
   }
   igraph_copy(sys->graph, G);
   
   sys->nodeStates        = nodeStates;
   sys->edgeStates        = edgeStates;
   sys->nodeDynVary       = dynNodeVary;
   sys->edgeDynVary       = dynEdgeVary;
   sys->nodeDynParamsVary = dynNodeParamsVary;
   sys->edgeDynParamsVary = dynEdgeParamsVary;
   
   if ((sys->nodeParams = (ne_param_t *)calloc(NE_NODE_PARAMS * igraph_vcount(G), sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.nodeParams.");
      igraph_destroy(sys->graph);
      free(sys);
      return NULL;
   }
   
   if ((sys->edgeParams = (ne_param_t *)calloc(NE_EDGE_PARAMS * igraph_ecount(G), 
                                               sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.edgeParams.");
      igraph_destroy(sys->graph);
      free(sys->nodeParams);
      free(sys);
      return NULL;
   }
   
   if (dynNodeVary == TRUE) {
      if ((sys->nodeDyn = (ne_dyn_node_t **)malloc(sizeof(ne_dyn_node_t *) * igraph_vcount(G))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.nodeDyn.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys);
         return NULL;
      }
   }
   else {
      if ((sys->nodeDyn = (ne_dyn_node_t **)malloc(sizeof(ne_dyn_node_t *))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.nodeDyn.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys);
         return NULL;
      }
   }
   
   if (dynEdgeVary == TRUE) {
      if ((sys->edgeDyn = (ne_dyn_edge_t **)malloc(sizeof(ne_dyn_edge_t *) * 
                                                            igraph_ecount(G))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.edgeDyn.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys);
         return NULL;
      }
   }
   else {
      if ((sys->edgeDyn = (ne_dyn_edge_t **)malloc(sizeof(ne_dyn_edge_t *))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.edgeDyn.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys);
         return NULL;
      }
   }
   
   if (dynNodeParamsVary == TRUE) {
      if ((sys->nodeDynParams = (ne_dyn_param_t *)calloc(igraph_vcount(G) * NE_DYN_PARAMS, 
                                                         sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.nodeDynParams.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys->edgeDyn);
         free(sys);
         return NULL;
      }     
   }
   else {
      if ((sys->nodeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS, 
                                                         sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.nodeDynParams.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys->edgeDyn);
         free(sys);
         return NULL;
      }
   }
   
   if (dynEdgeParamsVary == TRUE) {
      if ((sys->edgeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS * igraph_ecount(G), 
                                                         sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.edgeDynParams.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys->edgeDyn);
         free(sys->nodeDynParams);
         free(sys);
         return NULL;
      }     
   }
   else {
      if ((sys->edgeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS, 
                                                         sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_alloc_from_igraph: Could not allocate memory for ne_sys.edgeDynParams.");
         igraph_destroy(sys->graph);
         free(sys->nodeParams);
         free(sys->edgeParams);
         free(sys->nodeDyn);
         free(sys->edgeDyn);
         free(sys->nodeDynParams);
         free(sys);
         return NULL;
      }
   }
   
   return sys;
}


ne_sys_t * ne_sys_clone (ne_sys_t *from) 
{  
   ne_sys_t *to;
   
   /* Initialise the copy with details of the graph */
   to = ne_sys_alloc_from_igraph(from->graph, from->nodeStates, from->edgeStates, 
                           from->nodeDynVary, from->edgeDynVary, from->nodeDynParamsVary, 
                          from->edgeDynParamsVary);
   
   if (to == NULL) {
      ne_error("ne_sys_copy: Could not copy system.");
      return NULL;
   }

   /* Bulk copy the arrays of dynamic data */
   memcpy(to->nodeParams, from->nodeParams, sizeof(ne_param_t) * igraph_vcount(from->graph) * NE_NODE_PARAMS);
   memcpy(to->edgeParams, from->edgeParams, sizeof(ne_param_t) * igraph_ecount(from->graph) * NE_EDGE_PARAMS);
   
   if (from->nodeDynVary == TRUE) {
      memcpy(to->nodeDyn, from->nodeDyn, sizeof(ne_dyn_node_t) * igraph_vcount(from->graph));
   }
   else {
      to->nodeDyn[0] = from->nodeDyn[0];
   }
   
   if (from->edgeDynVary == TRUE) {
      memcpy(to->edgeDyn, from->edgeDyn, sizeof(ne_dyn_edge_t) * igraph_ecount(from->graph));
   }
   else {
      to->edgeDyn[0] = from->edgeDyn[0];
   }
   
   if (from->nodeDynParamsVary == TRUE) {
      memcpy(to->nodeDynParams, from->nodeDynParams, sizeof(ne_dyn_param_t) * 
             igraph_vcount(from->graph) * NE_DYN_PARAMS);
   }
   else {
      memcpy(to->nodeDynParams, from->nodeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   
   if (from->edgeDynParamsVary == TRUE) {
      memcpy(to->edgeDynParams, from->edgeDynParams, sizeof(ne_dyn_param_t) * 
             igraph_ecount(from->graph) * NE_DYN_PARAMS);
   }
   else {
      memcpy(to->edgeDynParams, from->edgeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   
   return to;
}


void ne_sys_free (ne_sys_t *sys)
{
   igraph_destroy(sys->graph);
   free(sys->nodeParams);
   free(sys->edgeParams);
   free(sys->nodeDyn);
   free(sys->edgeDyn);
   free(sys->nodeDynParams);
   free(sys->edgeDynParams);
   free(sys);
}


igraph_t * ne_sys_igraph (ne_sys_t *sys)
{
   return sys->graph;
}


ne_bool_t ne_sys_directed (ne_sys_t *sys) 
{
   return (ne_bool_t)igraph_is_directed(sys->graph);
}


ne_int_t ne_sys_node_states (ne_sys_t *sys)
{
   return sys->nodeStates;
}

ne_int_t ne_sys_edge_states (ne_sys_t *sys)
{
   return sys->edgeStates;
}


ne_int_t ne_sys_nodes (ne_sys_t *sys)
{
   return (ne_int_t)igraph_vcount(sys->graph);
}


ne_int_t ne_sys_edges (ne_sys_t *sys)
{
   return (ne_int_t)igraph_ecount(sys->graph);
}



ne_bool_t ne_sys_dyn_node_vary (ne_sys_t *sys)
{
   return sys->nodeDynVary;
}


ne_bool_t ne_sys_dyn_edge_vary (ne_sys_t *sys)
{
   return sys->edgeDynVary;
}


ne_bool_t ne_sys_dyn_node_param_vary (ne_sys_t *sys)
{
   return sys->nodeDynParamsVary;
}


ne_bool_t ne_sys_dyn_edge_param_vary (ne_sys_t *sys)
{
   return sys->edgeDynParamsVary;
}


ne_bool_t ne_sys_edge_exists (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode)
{
   igraph_integer_t i;
   igraph_vector_t adjV;
   igraph_vector_init(&adjV, 0);
   
   /* Decide which edges to extract from graph */
   if (!igraph_is_directed(sys->graph)) {
      /* Only consider inward and outward edges */
      igraph_neighbors(sys->graph, &adjV, (igraph_integer_t)startNode, IGRAPH_ALL);
   }
   else {   
      /* Only consider outward edges */
      igraph_neighbors(sys->graph, &adjV, (igraph_integer_t)startNode, IGRAPH_OUT);
   }
   
   /* Check each adjacent edge */
   for (i=0; i<igraph_vector_size(&adjV); i++) {
      /* Check to see if correct edge found */
      if (VECTOR(adjV)[(int)i] == (igraph_integer_t)endNode) {
         /* Found the edge, stop searching and return success */
         return TRUE;
      }
   }
   
   igraph_vector_destroy(&adjV);
   
   /* Could not find an edge */
   return FALSE;
}


ne_err_code_t ne_sys_add_node (ne_sys_t *sys, ne_dyn_node_t *dyn, ne_dyn_param_t *params)
{
   ne_param_t *newNodeParams;
   ne_dyn_param_t *newNodeDynParams;
   ne_dyn_node_t **newNodeDyn;
   ne_int_t nodes;
   nodes = (ne_int_t)igraph_vcount(sys->graph);
   
   /* Allocate memory for new parameter arrays */
   if ((newNodeParams = (ne_param_t *)calloc(NE_NODE_PARAMS * (nodes + 1), 
                                             sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_add_node: Could not allocate memory for ne_sys.nodeParams.");
      return NE_MEM_ALLOC_ERROR;
   }
   if (sys->nodeDynParamsVary == TRUE) {
      if ((newNodeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS * (nodes + 1), 
                                                       sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_add_node: Could not allocate memory for ne_sys.nodeDynParams.");
         free(newNodeParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newNodeDynParams, sys->nodeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * nodes);
      free(sys->nodeDynParams);
      sys->nodeDynParams = newNodeDynParams;
   }
   if (sys->nodeDynVary == TRUE) {
      if ((newNodeDyn = (ne_dyn_node_t **)calloc(nodes+1, sizeof(ne_dyn_node_t *))) == NULL) {
         ne_error("ne_sys_add_node: Could not allocate memory for ne_sys.nodeDyn.");
         free(newNodeParams);
         if (sys->nodeDynParamsVary == TRUE)
            free(newNodeDynParams);
         return NE_MEM_ALLOC_ERROR;
      }
      /* Bulk copy the node dynamics */
      memcpy(newNodeDyn, sys->nodeDyn, sizeof(ne_dyn_node_t *) * nodes);
      free(sys->nodeDyn);
      sys->nodeDyn = newNodeDyn;
   }
   
   /* Bulk copy the parameters */
   memcpy(newNodeParams, sys->nodeParams, sizeof(ne_param_t) * NE_NODE_PARAMS * nodes);
   free(sys->nodeParams);
   sys->nodeParams = newNodeParams;

   /* Update the graph */
   igraph_add_vertices(sys->graph, (igraph_integer_t)1, 0);
   
   ne_sys_set_dyn_node(sys, nodes, dyn);
   ne_sys_set_dyn_node_params(sys, nodes, params);
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_add_edge (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                               ne_dyn_edge_t *dyn, ne_dyn_param_t *params)
{
   ne_param_t *newEdgeParams;
   ne_dyn_param_t *newEdgeDynParams;
   ne_dyn_edge_t **newEdgeDyn;
   ne_int_t edges;
   edges = (ne_int_t)igraph_ecount(sys->graph);
   
   if (ne_sys_edge_exists(sys, startNode, endNode)) {
      ne_error("ne_sys_add_edge: Edge (%i, %i) already exists.", startNode, endNode);
      return NE_FAILURE;
   }
   
   /* Allocate memory for new parameter arrays */
   if ((newEdgeParams = (ne_param_t *)calloc(NE_EDGE_PARAMS * (edges + 1), 
                                             sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_add_edge: Could not allocate memory for ne_sys.edgeParams.");
      return NE_MEM_ALLOC_ERROR;
   }
   if (sys->edgeDynParamsVary == TRUE) {
      if ((newEdgeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS * (edges + 1), 
                                                       sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_add_edge: Could not allocate memory for ne_sys.edgeDynParams.");
         free(newEdgeParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newEdgeDynParams, sys->edgeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * edges);
      free(sys->edgeDynParams);
      sys->edgeDynParams = newEdgeDynParams;
   }
   if (sys->edgeDynVary == TRUE) {
      if ((newEdgeDyn = (ne_dyn_edge_t **)calloc(edges+1, sizeof(ne_dyn_edge_t *))) == NULL) {
         ne_error("ne_sys_add_edge: Could not allocate memory for ne_sys.edgeDyn.");
         free(newEdgeParams);
         if (sys->edgeDynParamsVary == TRUE)
            free(newEdgeDynParams);
         return NE_MEM_ALLOC_ERROR;
      }
      /* Bulk copy the edge dynamics */
      memcpy(newEdgeDyn, sys->edgeDyn, sizeof(ne_dyn_edge_t *) * edges);
      free(sys->edgeDyn);
      sys->edgeDyn = newEdgeDyn;
   }
   
   /* Bulk copy the parameters */
   memcpy(newEdgeParams, sys->edgeParams, sizeof(ne_param_t) * NE_EDGE_PARAMS * edges);
   free(sys->edgeParams);
   sys->edgeParams = newEdgeParams;
   
   /* Update the graph */
   igraph_add_edge(sys->graph, (igraph_integer_t)startNode, (igraph_integer_t)endNode);
   
   ne_sys_set_dyn_edge(sys, startNode, endNode, dyn);
   ne_sys_set_dyn_edge_params(sys, startNode, endNode, params);
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_del_node (ne_sys_t *sys, ne_int_t node)
{
   ne_param_t *newNodeParams;
   ne_dyn_param_t *newNodeDynParams;
   ne_dyn_node_t **newNodeDyn;
   igraph_es_t esDel;
   igraph_eit_t eitDel;
   ne_int_t nodes;
   nodes = (ne_int_t)igraph_vcount(sys->graph);
   
   if (nodes == 1) {
      ne_error("ne_sys_del_node: Cannot delete node when only 1 left.");
      return NE_FAILURE;
   }
   
   /* Cycle through each connected edge and remove */
   /* This will be slow as it requires reallocation of memory at each step */
   igraph_es_adj(&esDel, (igraph_integer_t)node, IGRAPH_ALL);
   igraph_eit_create(sys->graph, esDel, &eitDel);
   while (!IGRAPH_EIT_END(eitDel)) {
      ne_sys_del_eid(sys, (ne_int_t)IGRAPH_EIT_GET(eitDel));
      IGRAPH_EIT_NEXT(eitDel);
   }
   
   /* Allocate memory for new parameter arrays */
   if ((newNodeParams = (ne_param_t *)calloc(NE_NODE_PARAMS * (nodes - 1), 
                                   sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_del_node: Could not allocate memory for ne_sys.nodeParams.");
      return NE_MEM_ALLOC_ERROR;
   }
   
   if (sys->nodeDynParamsVary == TRUE) {
      if ((newNodeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS * (nodes - 1), 
                                           sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_del_node: Could not allocate memory for ne_sys.nodeDynParams.");
         free(newNodeParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newNodeDynParams, sys->nodeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * node);
      memcpy(&newNodeDynParams[node], &sys->nodeDynParams[node+1], sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * (nodes-1-node));
      free(sys->nodeDynParams);
      sys->nodeDynParams = newNodeDynParams;
   }
   if (sys->nodeDynVary == TRUE) {
      if ((newNodeDyn = (ne_dyn_node_t **)calloc(nodes-1, sizeof(ne_dyn_node_t *))) == NULL) {
         ne_error("ne_sys_del_node: Could not allocate memory for ne_sys.nodeDyn.");
         free(newNodeParams);
         if (sys->nodeDynParamsVary == TRUE) 
            free(newNodeDynParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newNodeDyn, sys->nodeDyn, sizeof(ne_dyn_node_t *) * node);
      memcpy(&newNodeDyn[node], &sys->nodeDyn[node+1], sizeof(ne_dyn_node_t *) * (nodes-1-node));
      free(sys->nodeDyn);
      sys->nodeDyn = newNodeDyn;
   }
   
   /* Bulk copy removing the element that is removed */
   memcpy(newNodeParams, sys->nodeParams, sizeof(ne_param_t) * NE_NODE_PARAMS * node);
   memcpy(&newNodeParams[node], &sys->nodeParams[node+1], sizeof(ne_param_t) * NE_NODE_PARAMS * (nodes-1-node));
   
   /* Swap arrays and free old ones */
   free(sys->nodeParams);
   sys->nodeParams = newNodeParams;

   /* Remove node from graph */
   igraph_delete_vertices(sys->graph, igraph_vss_1((igraph_integer_t)node));
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_del_edge (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode)
{
   igraph_integer_t eidDel;
   
   /* Edge ID required to calculate indexes that need to be copied */
   igraph_get_eid(sys->graph, &eidDel, (igraph_integer_t)startNode, (igraph_integer_t)endNode, 
                  igraph_is_directed(sys->graph));
   
   return ne_sys_del_eid(sys, (ne_int_t)eidDel);   
}


ne_err_code_t ne_sys_del_eid (ne_sys_t *sys, ne_int_t eid)
{
   ne_param_t *newEdgeParams;
   ne_dyn_param_t *newEdgeDynParams;
   ne_dyn_edge_t **newEdgeDyn;
   ne_int_t edges;
   edges = (ne_int_t)igraph_ecount(sys->graph);
   
   /* Allocate memory for new parameter arrays */
   if ((newEdgeParams = (ne_param_t *)calloc(NE_EDGE_PARAMS * (edges - 1), 
                                             sizeof(ne_param_t))) == NULL) {
      ne_error("ne_sys_del_edge: Could not allocate memory for ne_sys.edgeParams.");
      return NE_MEM_ALLOC_ERROR;
   }
   if (sys->edgeDynParamsVary == TRUE) {
      if ((newEdgeDynParams = (ne_dyn_param_t *)calloc(NE_DYN_PARAMS * (edges - 1), 
                                                       sizeof(ne_dyn_param_t))) == NULL) {
         ne_error("ne_sys_del_edge: Could not allocate memory for ne_sys.edgeDynParams.");
         free(newEdgeParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newEdgeDynParams, sys->edgeDynParams, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * eid);
      memcpy(&newEdgeDynParams[eid], &sys->edgeDynParams[eid+1], sizeof(ne_dyn_param_t) * NE_DYN_PARAMS * (edges-1-eid));
      free(sys->edgeDynParams);
      sys->edgeDynParams = newEdgeDynParams;
      
   }
   if (sys->edgeDynVary == TRUE) {
      if ((newEdgeDyn = (ne_dyn_edge_t **)calloc(edges-1, sizeof(ne_dyn_edge_t *))) == NULL) {
         ne_error("ne_sys_del_edge: Could not allocate memory for ne_sys.edgeDyn.");
         free(newEdgeParams);
         if (sys->edgeDynParamsVary == TRUE)
            free(newEdgeDynParams);
         return NE_MEM_ALLOC_ERROR;
      }
      memcpy(newEdgeDyn, sys->edgeDyn, sizeof(ne_dyn_edge_t *) * eid);
      memcpy(&newEdgeDyn[eid], &sys->edgeDyn[eid+1], sizeof(ne_dyn_edge_t *) * (edges-1-eid));
      free(sys->edgeDyn);
      sys->edgeDyn = newEdgeDyn;
   }
   
   /* Bulk copy removing the element that is removed */
   memcpy(newEdgeParams, sys->edgeParams, sizeof(ne_param_t) * NE_EDGE_PARAMS * eid);
   memcpy(&newEdgeParams[eid], &sys->edgeParams[eid+1], sizeof(ne_param_t) * NE_EDGE_PARAMS * (edges-1-eid));
   
   /* Swap arrays and free old ones */
   free(sys->edgeParams);
   sys->edgeParams = newEdgeParams;
   
   /* Remove edge from graph */
   igraph_delete_edges(sys->graph, igraph_ess_1((igraph_integer_t)eid));
   
   return NE_SUCCESS;   
}


ne_param_t ne_sys_node_param (ne_sys_t *sys, ne_int_t node, ne_int_t param)
{
   return sys->nodeParams[(node * NE_NODE_PARAMS) + param];
}


ne_param_t ne_sys_edge_param (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, ne_int_t param)
{
   igraph_integer_t eid;
   /* Find the edge ID used for the index in the parameters array */
   igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
   return sys->edgeParams[((ne_int_t)eid * NE_EDGE_PARAMS) + param];
}


ne_param_t ne_sys_eid_param (ne_sys_t *sys, ne_int_t eid, ne_int_t param)
{
   return sys->edgeParams[(eid * NE_EDGE_PARAMS) + param];
}


ne_err_code_t ne_sys_set_node_param (ne_sys_t *sys, ne_int_t node, ne_int_t param, ne_param_t value)
{
   sys->nodeParams[(node * NE_NODE_PARAMS) + param] = value;
   return NE_SUCCESS;

}


ne_err_code_t ne_sys_set_edge_param (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                                     ne_int_t param, ne_param_t value)
{
   igraph_integer_t eid;
   /* Find the edge ID used for the index in the parameters array */
   igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
   sys->edgeParams[((ne_int_t)eid * NE_EDGE_PARAMS) + param] = value;
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_eid_param (ne_sys_t *sys, ne_int_t eid, ne_int_t param, 
                                    ne_param_t value)
{
   sys->edgeParams[(eid * NE_EDGE_PARAMS) + param] = value;
   return NE_SUCCESS;
}


ne_dyn_node_t * ne_sys_dyn_node (ne_sys_t *sys, ne_int_t node)
{
   if (sys->nodeDynVary == TRUE) {
      return sys->nodeDyn[node];
   }
   else {
      return sys->nodeDyn[0];
   }
}


ne_dyn_edge_t * ne_sys_dyn_edge (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode)
{
   igraph_integer_t eid;
   
   if (sys->edgeDynVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
      return sys->edgeDyn[(ne_int_t)eid];
   }
   else {
      return sys->edgeDyn[0];
   }
}


ne_dyn_edge_t * ne_sys_dyn_eid (ne_sys_t *sys, ne_int_t eid)
{
   if (sys->edgeDynVary == TRUE) {
      return sys->edgeDyn[eid];
   }
   else {
      return sys->edgeDyn[0];
   }
}


ne_err_code_t ne_sys_set_dyn_node (ne_sys_t *sys, ne_int_t node, ne_dyn_node_t *nodeDyn)
{
   if (sys->nodeDynVary == TRUE) {
      sys->nodeDyn[node] = nodeDyn;
   }
   else {
      sys->nodeDyn[0] = nodeDyn;
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_edge (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                                   ne_dyn_edge_t *edgeDyn)
{
   igraph_integer_t eid;
   
   if (sys->edgeDynVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
      sys->edgeDyn[(ne_int_t)eid] = edgeDyn;
   }
   else {
      sys->edgeDyn[0] = edgeDyn;
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_eid (ne_sys_t *sys, ne_int_t eid, ne_dyn_edge_t *edgeDyn)
{
   if (sys->edgeDynVary == TRUE) {
      sys->edgeDyn[eid] = edgeDyn;
   }
   else {
      sys->edgeDyn[0] = edgeDyn;
   }
   
   return NE_SUCCESS;
}


ne_dyn_param_t ne_sys_dyn_node_param (ne_sys_t *sys, ne_int_t node, ne_int_t param)
{
   if (sys->nodeDynParamsVary == TRUE) {
      return sys->nodeDynParams[(node * NE_DYN_PARAMS) + param];
   }
   else {
      return sys->nodeDynParams[param];
   }
}


ne_dyn_param_t * ne_sys_dyn_node_params (ne_sys_t *sys, ne_int_t node)
{
   if (sys->nodeDynParamsVary == TRUE) {
      return &sys->nodeDynParams[node * NE_DYN_PARAMS];
   }
   else {
      return sys->nodeDynParams;
   }
}


ne_dyn_param_t ne_sys_dyn_edge_param (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                                      ne_int_t param)
{
   igraph_integer_t eid;
   
   if (sys->edgeDynParamsVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
      return sys->edgeDynParams[((ne_int_t)eid * NE_DYN_PARAMS) + param];
   }
   else {
      return sys->edgeDynParams[param];
   }
}


ne_dyn_param_t ne_sys_dyn_eid_param (ne_sys_t *sys, ne_int_t eid, ne_int_t param)
{
   if (sys->edgeDynParamsVary == TRUE) {
      return sys->edgeDynParams[(eid * NE_DYN_PARAMS) + param];
   }
   else {
      return sys->edgeDynParams[param];
   }
}


ne_dyn_param_t * ne_sys_dyn_edge_params (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode)
{  
   igraph_integer_t eid;
   
   if (sys->edgeDynParamsVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
      return &sys->edgeDynParams[(ne_int_t)eid * NE_DYN_PARAMS];
   }
   else {
      return sys->edgeDynParams;
   }
}


ne_dyn_param_t * ne_sys_dyn_eid_params  (ne_sys_t *sys, ne_int_t eid)
{
   if (sys->edgeDynParamsVary == TRUE) {
      return &sys->edgeDynParams[eid * NE_DYN_PARAMS];
   }
   else {
      return sys->edgeDynParams;
   }
}


ne_err_code_t ne_sys_set_dyn_node_param (ne_sys_t *sys, ne_int_t node, ne_int_t param, 
                                         ne_dyn_param_t value)
{  
   if (sys->nodeDynParamsVary == TRUE) {
      sys->nodeDynParams[(node * NE_DYN_PARAMS) + param] = value;
   }
   else {
      sys->nodeDynParams[param] = value;
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_edge_param (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                               ne_int_t param, ne_dyn_param_t value)
{  
   igraph_integer_t eid;
   
   if (sys->edgeDynParamsVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, igraph_is_directed(sys->graph));
      sys->edgeDynParams[((ne_int_t)eid * NE_DYN_PARAMS) + param] = value;
   }
   else {
      sys->edgeDynParams[param] = value;
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_eid_param (ne_sys_t *sys, ne_int_t eid, ne_int_t param, 
                              ne_dyn_param_t value)
{  
   if (sys->edgeDynParamsVary == TRUE) {
      sys->edgeDynParams[(eid * NE_DYN_PARAMS) + param] = value;
   }
   else {
      sys->edgeDynParams[param] = value;
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_node_params (ne_sys_t *sys, ne_int_t node, ne_dyn_param_t *params)
{
   if (sys->nodeDynParamsVary == TRUE) {
      memcpy(&sys->nodeDynParams[node*NE_DYN_PARAMS], params, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   else {
      memcpy(sys->nodeDynParams, params, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_edge_params (ne_sys_t *sys, ne_int_t startNode, ne_int_t endNode, 
                                          ne_dyn_param_t *params)
{
   igraph_integer_t eid;
   
   if (sys->edgeDynParamsVary == TRUE) {
      /* Find the edge ID used for the index in the parameters array */
      igraph_get_eid(sys->graph, &eid, (igraph_integer_t)startNode, (igraph_integer_t)endNode, 
                     igraph_is_directed(sys->graph));
      return ne_sys_set_dyn_eid_params(sys, (ne_int_t)eid, params);
   }
   else {
      memcpy(sys->edgeDynParams, params, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_set_dyn_eid_params (ne_sys_t *sys, ne_int_t eid, ne_dyn_param_t *params)
{
   if (sys->edgeDynParamsVary == TRUE) {
      memcpy(&sys->edgeDynParams[eid*NE_DYN_PARAMS], params, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   else {
      memcpy(sys->edgeDynParams, params, sizeof(ne_dyn_param_t) * NE_DYN_PARAMS);
   }
   return NE_SUCCESS;
}


ne_err_code_t ne_sys_to_GML (ne_sys_t *sys, FILE *file)
{
   int i, j;
   char paramBuf[100];
   igraph_t *G;
   
   G = sys->graph;
   
   /* Remove all current attributes (none should exist) */
   DELALL(G);
   
   /* Set system level attributes */
   SETVAN(G, "neNodeStates", 0, sys->nodeStates);
   SETVAN(G, "neEdgeStates", 0, sys->edgeStates);
   SETVAN(G, "neNodeDynVary", 0, sys->nodeDynVary);
   SETVAN(G, "neEdgeDynVary", 0, sys->edgeDynVary);
   SETVAN(G, "neNodeDynParamsVary", 0, sys->nodeDynParamsVary);
   SETVAN(G, "neEdgeDynParamsVary", 0, sys->edgeDynParamsVary);
   
   /* Set all node parameters */
   for (i=0; i<ne_sys_nodes(sys); i++) {
      for (j=0; j<NE_NODE_PARAMS; j++) {
         sprintf(paramBuf, "neNodeParam%i", j);
         SETVAN(G, paramBuf, i, ne_sys_node_param(sys, i, j));
      }
   }
   
   /* Set node dynamics parameters */
   if (ne_sys_dyn_node_param_vary(sys)) {
      for (i=0; i<ne_sys_nodes(sys); i++) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neNodeDynParam%i", j);
            SETVAN(G, paramBuf, i, ne_sys_dyn_node_param(sys, i, j));
         }
      }
   }
   else {
      for (j=0; j<NE_DYN_PARAMS; j++) {
         sprintf(paramBuf, "neNodeDynParam%i", j);
         SETVAN(G, paramBuf, 0, ne_sys_dyn_node_param(sys, 0, j));
      }
   }
   
   /* Set node dynamics */
   if (ne_sys_dyn_node_vary(sys)) {
      for (i=0; i<ne_sys_nodes(sys); i++) {
         SETVAS(G, "neNodeDyn", i, ne_sys_dyn_node(sys, i)->name);
      }
   }
   else {
      SETVAS(G, "neNodeDyn", 0, ne_sys_dyn_node(sys, 0)->name);
   }
   
   /* Set all edge parameters */
   for (i=0; i<ne_sys_edges(sys); i++) {
      for (j=0; j<NE_EDGE_PARAMS; j++) {
         sprintf(paramBuf, "neEdgeParam%i", j);
         SETEAN(G, paramBuf, i, ne_sys_eid_param(sys, i, j));
      }
   }
   
   /* Set edge dynamics parameters */
   if (ne_sys_dyn_edge_param_vary(sys)) {
      for (i=0; i<ne_sys_edges(sys); i++) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neEdgeDynParam%i", j);
            SETEAN(G, paramBuf, i, ne_sys_dyn_eid_param(sys, i, j));
         }
      }
   }
   else {
      if (ne_sys_edges(sys) > 0) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neEdgeDynParam%i", j);
            SETEAN(G, paramBuf, 0, ne_sys_dyn_eid_param(sys, 0, j));
         }
      }
   }
   
   /* Set edge dynamics */
   if (ne_sys_dyn_edge_vary(sys)) {
      for (i=0; i<ne_sys_edges(sys); i++) {
         SETEAS(G, "neEdgeDyn", i, ne_sys_dyn_eid(sys, i)->name);
      }
   }
   else {
      if (ne_sys_edges(sys) > 0) {
         SETEAS(G, "neEdgeDyn", 0, ne_sys_dyn_eid(sys, 0)->name);
      }
   }
   
   /* Save to file */
   igraph_write_graph_gml(G, file, NULL, "NetEvo");   
   
   /* Remove all attributes */
   DELALL(G);
   
   return NE_SUCCESS;
}


ne_sys_t * ne_sys_from_GML (FILE *file)
{
   int i, j;
   char paramBuf[100];
   igraph_t G;
   ne_sys_t *sys;
   
   /* Read the graph */
   igraph_read_graph_gml(&G, file);
   
   /* Generate a new system from this attributed file */
   sys = ne_sys_alloc_from_igraph(&G, 
                      (ne_int_t)VAN(&G, "neNodeStates", 0), 
                     (ne_int_t)VAN(&G, "neEdgeStates", 0), 
                     (ne_bool_t)VAN(&G, "neNodeDynVary", 0), 
                     (ne_bool_t)VAN(&G, "neEdgeDynVary", 0), 
                     (ne_bool_t)VAN(&G, "neNodeDynParamsVary", 0), 
                     (ne_bool_t)VAN(&G, "neEdgeDynParamsVary", 0));
   
   /* Load node parameters */
   for (i=0; i<ne_sys_nodes(sys); i++) {
      for (j=0; j<NE_NODE_PARAMS; j++) {
         sprintf(paramBuf, "neNodeParam%i", j);
         ne_sys_set_node_param(sys, i, j, (ne_param_t)VAN(&G, paramBuf, i));
      }
   }
   
   /* Set node dynamics parameters */
   if (ne_sys_dyn_node_param_vary(sys)) {
      for (i=0; i<ne_sys_nodes(sys); i++) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neNodeDynParam%i", j);
            ne_sys_set_dyn_node_param(sys, i, j, (ne_dyn_param_t)VAN(&G, paramBuf, i));
         }
      }
   }
   else {
      for (j=0; j<NE_DYN_PARAMS; j++) {
         sprintf(paramBuf, "neNodeDynParam%i", j);
         ne_sys_set_dyn_node_param(sys, 0, j, (ne_dyn_param_t)VAN(&G, paramBuf, 0));
      }
   }
   
   /* Set node dynamics */
   if (ne_sys_dyn_node_vary(sys)) {
      for (i=0; i<ne_sys_nodes(sys); i++) {
         ne_sys_set_dyn_node(sys, i, ne_dyn_find_node_dyn(VAS(&G, "neNodeDyn", i)));
      }
   }
   else {
      ne_sys_set_dyn_node(sys, 0, ne_dyn_find_node_dyn(VAS(&G, "neNodeDyn", 0)));
   }
   
   /* Load edge parameters */
   for (i=0; i<ne_sys_edges(sys); i++) {
      for (j=0; j<NE_EDGE_PARAMS; j++) {
         sprintf(paramBuf, "neEdgeParam%i", j);
         ne_sys_set_eid_param(sys, i, j, (ne_param_t)EAN(&G, paramBuf, i));
      }
   }
   
   /* Load edge dynamics parameters */
   if (ne_sys_dyn_edge_param_vary(sys)) {
      for (i=0; i<ne_sys_edges(sys); i++) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neEdgeDynParam%i", j);
            ne_sys_set_dyn_eid_param(sys, i, j, (ne_dyn_param_t)EAN(&G, paramBuf, i));
         }
      }
   }
   else {
      if (ne_sys_edges(sys) > 0) {
         for (j=0; j<NE_DYN_PARAMS; j++) {
            sprintf(paramBuf, "neEdgeDynParam%i", j);
            ne_sys_set_dyn_eid_param(sys, 0, j, (ne_dyn_param_t)EAN(&G, paramBuf, 0));
         }
      }
   }
   
   /* Load edge dynamics */
   if (ne_sys_dyn_edge_vary(sys)) {
      for (i=0; i<ne_sys_edges(sys); i++) {
         ne_sys_set_dyn_eid(sys, i, ne_dyn_find_edge_dyn(EAS(&G, "neEdgeDyn", i)));
      }
   }
   else {
      if (ne_sys_edges(sys) > 0) {
         ne_sys_set_dyn_eid(sys, 0, ne_dyn_find_edge_dyn(EAS(&G, "neEdgeDyn", 0)));
      }
   }
   
   /* Free memory and free any current attributes */
   igraph_destroy(&G);
   DELALL(sys->graph);
   
   return sys;
}


void ne_sys_print(ne_sys_t *sys)
{
   /* TODO: ne_sys_print */
   ne_log("ne_sys_print: Not yet implemented.");
}


gsl_matrix * ne_sys_adj_to_gsl_matrix_alloc (ne_sys_t *sys)
{
   gsl_matrix *M;
   ne_int_t i, j;
   
   M = gsl_matrix_calloc((size_t)ne_sys_nodes(sys), (size_t)ne_sys_nodes(sys));
   
   for (i=0; i<ne_sys_nodes(sys); i++) {
      for (j=0; j<ne_sys_nodes(sys); j++) {
         if (ne_sys_edge_exists(sys, i, j) == TRUE)
            gsl_matrix_set(M, (size_t)i, (size_t)j, 1.0);
      }
   }

   return M;
}


gsl_matrix * ne_sys_laplacian_to_gsl_matrix_alloc (ne_sys_t *sys)
{
   gsl_matrix *L;
   igraph_matrix_t igraphL;
   ne_int_t i, j;
   
   L = gsl_matrix_calloc((size_t)ne_sys_nodes(sys), (size_t)ne_sys_nodes(sys));
   
   igraph_matrix_init(&igraphL, 0, 0);
   
   /* Find the degree of all nodes */
   igraph_laplacian(sys->graph, &igraphL, FALSE);
   
   /* Might be quicker to query the igraph structure directly for each element */
   for (i=0; i<igraph_vcount(sys->graph); i++) {
      for (j=0; j<igraph_vcount(sys->graph); j++) {
         gsl_matrix_set(L, i, j, MATRIX(igraphL,i,j));
      }  
   }
   
   igraph_matrix_destroy(&igraphL);
   
   return L;
}