adding more description to result generating notebook and the simulation codes

simulation_code/Evaluate_ODE_solver_implementations/README.md

@@ -0,0 +1 @@
+Console application to evaluate different ODE solver strategies for the Izhikevich neuron model.

simulation_code/Evaluate_ODE_solver_implementations/src/CMakeLists.txt

@@ -0,0 +1,17 @@
+cmake_minimum_required (VERSION 2.8)
+
+# + + + set a project name, it is referenced in other cmake function calls
+project( izhikevich )
+
+# + + + set cache variables for a version number
+set( VERSION_MAJOR 1 )
+set( VERSION_MINOR 0 )
+
+# + + + add executable target to the project:   add_executable( <target_name> source1 [source 2] ...)  
+add_executable( izhikevich izhikevich izhikevich.cpp )
+
+# + + + specify  include directories
+include_directories( izhikevich /usr/include/plplot )
+
+# + + + ... and libraries
+target_link_libraries( izhikevich plplotcxxd plplotd gsl gslcblas )

simulation_code/Evaluate_ODE_solver_implementations/src/izhikevich.cpp

@@ -0,0 +1,878 @@
+/*
+ *  izhikevich.cpp
+ *
+ *  This file is part of the Izhikevich console application.
+ *
+ *  Copyright (C) 2016, Author: Guido Trensch
+ *
+ *  The Izhikevich console application 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 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  It 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 application. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   IZHIKEVICH MODEL
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =
+// =   compile with:
+// =   g++ main.cpp -o izhikevich -I/usr/include/plplot  -lplplotcxxd  -lplplotd
+// =
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+#include <plstream.h>
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <cstring>
+#include "izhikevich.h"
+#include "izhikevich_model.h"
+#include "ode_solvers.h"
+
+#define PRINT_DEBUG  printf
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   SIMULATION PARAMETERS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define SIMULATION_TIME                    (float)(600.0)    // in milliseconds
+#define INTEGRATION_STEP_SIZE              (float)(1.000)    // in milliseconds
+#define HIGHRES_INTEGRATION_STEP_SIZE      (float)(0.001)    // in milliseconds
+#define PLOT_SHIFT                         (float)(0.000)    // shift plot for debug, e.g. 5.0
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   FORWARD DECLARATIONS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void drawGraph( CNTXT* pCntxt, plstream* pPlStream, int numDataPoints, PLFLT* pTimesArray, PLFLT* pVoltagesArray, int plot_color, int plot_line_style, char* pLegendText );
+void exitProgram( CNTXT* pCntxt );
+void finalizeContext( CNTXT* pCntxt );
+void finalizePlot( CNTXT* pCntxt, plstream* pPlStream );
+void finalizeSimulation( CNTXT* pCntxt, PLFLT* pTimesArray, PLFLT* pVoltagesArray, float* pCurrentArray );
+void initializeContext( CNTXT** ppCntxt );
+void initializePlot( CNTXT* pCntxt, plstream** pPlStream );
+void initializeSimulation( CNTXT* pCntxt, float simulationTime, float integrationStepSize, PLFLT** ppTimesArray, PLFLT** ppVoltagesArray, PLFLT** ppFrequencyArray, float** ppCurrentArray );
+int  runSimulation( CNTXT* pCntxt, ENUM_ODESOLVER_METHOD method, float simulationTime, float integrationStepSize, PLFLT* pTimesArray, PLFLT* pVoltagesArray, PLFLT* pFrequencyArray, float* pCurrentArray );
+void setCurrentShape( CNTXT* pCntxt, float* pCurrentArray, int numDataPoints, float integrationStepSize );
+
+void processMainArguments( CNTXT* pCntxt, int argc, char** argv );
+void processOPT_C( CNTXT* pCntxt, int argc, char** argv );
+void processOPT_F( CNTXT* pCntxt, int argc, char** argv );
+void processOPT_H( CNTXT* pCntxt, int argc, char** argv );
+void processOPT_M( CNTXT* pCntxt, int argc, char** argv );
+void processOPT_N( CNTXT* pCntxt, int argc, char** argv );
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   MAIN ENTRY
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+int main (int argc, char** argv)
+{
+  CNTXT*     pCntxt          = NULL;
+  PLFLT*     pVoltagesArray  = NULL;
+  PLFLT*     pTimesArray     = NULL;
+  PLFLT*     pFrequencyArray = NULL;
+  float*     pCurrentArray   = NULL;
+  int        numDataPoints   = 0;
+  plstream*  pPlStream       = NULL;
+
+  initializeContext( &pCntxt );
+  processMainArguments( pCntxt, argc, argv );
+  initializePlot( pCntxt, &pPlStream );
+  
+  // PERFORM STANDARD EULER
+  if( pCntxt->args.fOdeSolverStandardEuler ) { 
+    PRINT_DEBUG( "    Debug: + + + perform standard Euler + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, STANDARD_EULER, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    // ***GTR print frequency
+    // drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pFrequencyArray, PLOT_COLOR_RED, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_STANDARD_EULER );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_RED, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_STANDARD_EULER );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // PERFORM SYMPLECTIC EULER
+  if( pCntxt->args.fOdeSolverSymplecticEuler ) {
+    PRINT_DEBUG( "    Debug: + + + symplectic Euler + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, SYMPLECTIC_EULER, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_GREEN, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_SYMPLECTIC_EULER );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }    
+
+  // PERFORM EXPLICIT EULER SPINNAKER
+  if( pCntxt->args.fOdeSolverEulerSpiNNaker ) {
+    PRINT_DEBUG( "    Debug: + + + perform explicit Euler - SpiNNaker implementation + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, EXPLICIT_EULER_SPINNAKER, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_MAGENTA, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_EULER_SPINNAKER );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // EXPLICIT EULER NEST
+  if( pCntxt->args.fOdeSolverEulerNEST ) {
+    PRINT_DEBUG( "    Debug: + + + explicit Euler - NEST implementation + + +\n" );  
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, EXPLICIT_EULER_NEST, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_YELLOW, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_EULER_NEST );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // EXPLICIT IZHIKEVICH NUMERICS NEST
+  if( pCntxt->args.fOdeSolverIzhikevichNEST ) {
+    PRINT_DEBUG( "    Debug: + + + explicit Izhikevich numerics - NEST implementation + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, EXPLICIT_IZHIKEVICH_NUMERICS_NEST, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_CYAN, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_IZHIKEVICH_NEST );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+ 
+  // GSL library
+  if( pCntxt->args.fOdeSolverGSL ) {
+    PRINT_DEBUG( "    Debug: + + + GSL library + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, GSL_LIBRARY, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_GREY, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_GSL );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // Heun
+  if( pCntxt->args.fOdeSolverHeun ) {
+    PRINT_DEBUG( "    Debug: + + + Heun's method + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, HEUN_METHOD, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_BROWN, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_HEUN );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // Euler with adaptive stepsize
+  if( pCntxt->args.fOdeSolverAdaptiveEuler ) { 
+    PRINT_DEBUG( "    Debug: + + + perform adaptive Euler + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, ADAPTIVE_EULER, SIMULATION_TIME, INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_BROWN, PLOT_LINE_FULL, STR_ODE_SOLVER_METHOD_ADAPTIVE_EULER );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  // HIGH RESOLUTION REFERENCE USING STANDARD EULER
+  if( pCntxt->args.fOdeSolverHighResolutionEuler ) {
+    PRINT_DEBUG( "    Debug: + + + reference high resolution standard Euler + + +\n" );
+    initializeSimulation( pCntxt, SIMULATION_TIME, HIGHRES_INTEGRATION_STEP_SIZE, &pTimesArray, &pVoltagesArray, &pFrequencyArray, &pCurrentArray );
+    numDataPoints = runSimulation( pCntxt, STANDARD_EULER, SIMULATION_TIME, HIGHRES_INTEGRATION_STEP_SIZE, pTimesArray, pVoltagesArray, pFrequencyArray, pCurrentArray );
+    drawGraph( pCntxt, pPlStream, numDataPoints, pTimesArray, pVoltagesArray, PLOT_COLOR_BLUE, PLOT_LINE_DOTTED, STR_ODE_SOLVER_METHOD_HIGH_RESOLUTION_EULER );
+    finalizeSimulation( pCntxt, pTimesArray, pVoltagesArray, pCurrentArray );
+  }
+
+  finalizePlot( pCntxt, pPlStream ); 
+  finalizeContext( pCntxt ); 
+  printf( "    Terminated normally\n" );
+  return 0;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+//   HELPER FUNCTIONS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+//   Plot graph, i.e. write to svg file
+//   Multiple graphs can be plotted and overlayed
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void drawGraph( CNTXT* pCntxt, plstream* pPlStream, int numDataPoints, PLFLT* pTimesArray, PLFLT* pVoltagesArray, int plot_color, int plot_line_style, char* pLegendText )
+{
+  PRINT_DEBUG( "    Debug: plot graph, i.e. write to svg file: %s\n", pCntxt->args.outputSvgFilename );
+  pllsty( plot_line_style );
+  plcol0( plot_color );
+  pPlStream->line( numDataPoints, pTimesArray, pVoltagesArray );
+  plptex( 80.0, pCntxt->legendYOffset, 0.0, 0.0, 0.5, pLegendText );
+  pCntxt->legendYOffset += 5;
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Free program context area and exit the program
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void exitProgram( CNTXT* pCntxt )
+{
+  finalizeContext( pCntxt );
+  exit( 0 );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Free program context area
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void finalizeContext( CNTXT* pCntxt )
+{
+  if( pCntxt != NULL ) {
+    free( pCntxt );
+  }
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Delete Plplot stream
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void finalizePlot( CNTXT* pCntxt, plstream* pPlStream )
+{
+  delete pPlStream;
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Free data arrays
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void finalizeSimulation( CNTXT* pCntxt, PLFLT* pTimesArray, PLFLT* pVoltagesArray, float* pCurrentArray )
+{
+  PRINT_DEBUG( "    Debug: finalize\n" );
+  free( pTimesArray );
+  free( pVoltagesArray );
+  free( pCurrentArray );
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Allocate and initialize the program's context area.
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void initializeContext( CNTXT** ppCntxt )
+{
+  PRINT_DEBUG( "    Debug: initialize context\n" );
+  *ppCntxt = (CNTXT*)malloc( S_CNTXT );
+  if( *ppCntxt == NULL ) {
+    printf( "ERRROR: Failed to allocate program context area\n");
+    exitProgram( NULL );
+  }
+  memset( *ppCntxt, NULLCHR, S_CNTXT );
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Set up Plplot
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void initializePlot( CNTXT* pCntxt, plstream** pPlStream )
+{
+  *pPlStream = new plstream();
+  plsdev("svg");
+  plsfnam( pCntxt->args.outputSvgFilename );
+  plscolbg( 255, 255, 255 );
+  plscmap0n(16);
+  (*pPlStream)->init();
+
+  (*pPlStream)->env( 0, SIMULATION_TIME, -80, 50, 0, 0 );
+  //                 |  |                 |   |   |  |
+  //                 |  |                 |   |   |  +---  axis - 0 draw axis box and numeric labels
+  //                 |  |                 |   |   +------  just - 0 axis scales independently
+  //                 |  |                 |   +----------  ymax
+  //                 |  |                 +--------------  ymin
+  //                 |  +--------------------------------  xmax
+  //                 +-----------------------------------  xmin
+
+  (*pPlStream)->lab( "Time [ms]", "Voltage [mv]", pCntxt->neuronTypeName );
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Allocate and set up data structures
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+void initializeSimulation( CNTXT* pCntxt, float simulationTime, float integrationStepSize, PLFLT** ppTimesArray, PLFLT** ppVoltagesArray, PLFLT** ppFrequencyArray, float** ppCurrentArray )
+{
+
+  PRINT_DEBUG( "    Debug: initialize simulation and prepare data arrays\n" );
+
+  PLFLT*  pVoltagesArray  = NULL;
+  PLFLT*  pTimesArray     = NULL;
+  PLFLT*  pFrequencyArray = NULL;
+  float*  pCurrentArray   = NULL;
+  int     numDataPoints   = floor(simulationTime / integrationStepSize) + 10;   // due to rounding errors there might be a few more datapoints returned from the simulation run; just add 10 to be on the save side
+  
+  PRINT_DEBUG( "    Debug: number of data points calculated: %d\n", numDataPoints );
+
+  pVoltagesArray = (PLFLT*)malloc( numDataPoints * sizeof(PLFLT) );
+  pTimesArray = (PLFLT*)malloc( numDataPoints * sizeof(PLFLT) );
+  pFrequencyArray = (PLFLT*)malloc( numDataPoints * sizeof(PLFLT) );
+  pCurrentArray = (float*)malloc( numDataPoints * sizeof(float) );
+
+  if( NULL == pVoltagesArray || NULL == pTimesArray || NULL == pFrequencyArray || NULL == pCurrentArray ) {
+    printf( "ERROR: failed to allocate memory\n" );
+    exitProgram( pCntxt );
+  }
+  memset( pVoltagesArray, 0x00, numDataPoints * sizeof(PLFLT) );
+  memset( pTimesArray, 0x00, numDataPoints * sizeof(PLFLT) );
+  memset( pFrequencyArray, 0x00, numDataPoints * sizeof(PLFLT) );
+  memset( pCurrentArray, 0x00, numDataPoints * sizeof(float) );
+
+  setCurrentShape( pCntxt, pCurrentArray, numDataPoints, integrationStepSize );
+
+  *ppTimesArray = pTimesArray;
+  *ppVoltagesArray = pVoltagesArray;
+  *ppFrequencyArray = pFrequencyArray;
+  *ppCurrentArray = pCurrentArray;
+
+  // set Izhikevich model constants
+  IZHIKEVICH_A = pCntxt->neuronParam_A;
+  IZHIKEVICH_B = pCntxt->neuronParam_B;
+  IZHIKEVICH_C = pCntxt->neuronParam_C;
+  IZHIKEVICH_D = pCntxt->neuronParam_D;
+  IZHIKEVICH_THR = pCntxt->neuronParam_THR;        
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Depending on parameter "method", call the corresponding ODE solver and
+//   perform the simulation 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+int runSimulation( CNTXT* pCntxt, ENUM_ODESOLVER_METHOD method, float simulationTime, float integrationStepSize
+                 , PLFLT* pTimesArray, PLFLT* pVoltagesArray, PLFLT* pFrequencyArray, float* pCurrentArray )
+{
+
+  PFUNCD pODE1dxdt = (PFUNCD)dvdt;      // Izhikevich ODE1 
+  PFUNCD pODE2dydt = (PFUNCD)dudt;      // Izhikevich ODE2 
+  float  v_t = pCntxt->neuronParam_V_INIT;
+  float  u_t = pCntxt->neuronParam_U_INIT;
+  float  v_tplus1 = 0;
+  float  u_tplus1 = 0;
+  float  t = 0;
+  int    dataPointIndex = 0;
+
+  PRINT_DEBUG( "    Debug: run simulation\n" );
+
+  switch( method ) {
+    case STANDARD_EULER:
+      t = PLOT_SHIFT * 1;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        ode2dSolver_StandardEuler( pODE1dxdt, pODE2dydt, integrationStepSize, v_t, u_t, &v_tplus1, &u_tplus1, &pCurrentArray[dataPointIndex] ); 
+        if( threshold( &v_tplus1, &u_tplus1 )) {
+          pFrequencyArray[dataPointIndex] = 30;
+        };
+        v_t = v_tplus1;
+        u_t = u_tplus1;
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case SYMPLECTIC_EULER:
+      t = PLOT_SHIFT * 2;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        ode2dSolver_SymplecticEuler( pODE1dxdt, pODE2dydt, integrationStepSize, v_t, u_t, &v_tplus1, &u_tplus1, &pCurrentArray[dataPointIndex] );
+        threshold( &v_tplus1, &u_tplus1 );
+        v_t = v_tplus1;
+        u_t = u_tplus1;
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case EXPLICIT_EULER_SPINNAKER:
+      t = PLOT_SHIFT * 3;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        explizitSpiNNaker( integrationStepSize, &v_t, &u_t, pCurrentArray[dataPointIndex] );
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case EXPLICIT_EULER_NEST:
+      t = PLOT_SHIFT * 4;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        explizitNEST1( integrationStepSize, &v_t, &u_t, pCurrentArray[dataPointIndex] );
+        t += integrationStepSize;
+        dataPointIndex++;
+      }    
+      break;
+
+    case EXPLICIT_IZHIKEVICH_NUMERICS_NEST:
+      t = PLOT_SHIFT * 5;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        explizitNEST2( integrationStepSize, &v_t, &u_t, pCurrentArray[dataPointIndex] );
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case GSL_LIBRARY:
+      t = PLOT_SHIFT * 6;
+      gslODESolver_Init();
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        gslODESolver_Integrate( integrationStepSize, &v_t, &u_t, pCurrentArray[dataPointIndex] );
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case HEUN_METHOD:
+      t = PLOT_SHIFT * 7;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        ode2dSolver_Heun( pODE1dxdt, pODE2dydt, integrationStepSize, v_t, u_t, &v_tplus1, &u_tplus1, &pCurrentArray[dataPointIndex] ); 
+        threshold( &v_tplus1, &u_tplus1 );
+        v_t = v_tplus1;
+        u_t = u_tplus1;
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+    case ADAPTIVE_EULER:
+      t = PLOT_SHIFT * 8;
+      while( t <= simulationTime ) {
+        pTimesArray[dataPointIndex] = t;     
+        pVoltagesArray[dataPointIndex] = v_t;
+        ode2dSolver_AdaptiveEuler( pODE1dxdt, pODE2dydt, integrationStepSize, v_t, u_t, &v_tplus1, &u_tplus1, &pCurrentArray[dataPointIndex] ); 
+        threshold( &v_tplus1, &u_tplus1 );
+        v_t = v_tplus1;
+        u_t = u_tplus1;
+        t += integrationStepSize;
+        dataPointIndex++;
+      }
+      break;
+
+
+    default:
+      printf( "ERROR: ODE solver method not defined\n");
+      exitProgram( pCntxt );
+  }
+  
+  PRINT_DEBUG( "    Debug: number of data points: %d\n", dataPointIndex );
+  return( dataPointIndex );  // return number of data points, hence the dataPointIndex
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Fill the current-data-array with the desired current-shape.
+//   The size of the array, i.e. the number od datapoints depend on the 
+//   resolution, hence the integration step size.
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void setCurrentShape( CNTXT* pCntxt, float* pCurrentArray, int numDataPoints, float integrationStepSize )
+{
+
+  PRINT_DEBUG( "    Debug: initialize current-shape array\n" );
+
+  switch( pCntxt->args.currentShape ) {
+
+    case CURRENT_STEP_150: 
+      for( int i = 0; i < numDataPoints; ++i ) {
+        pCurrentArray[i] = pCntxt->neuronParam_I_EXT;
+      }
+      for( int i = 0; i < floor( 150 / integrationStepSize ); ++i ) {
+        pCurrentArray[i] = 0;
+      }
+      // printf( "             i     +------------------------      \n" );
+      // printf( "             0 ----+                              \n" );
+      // printf( "                  150                    t [ms]   \n" );  
+      break;
+
+    case CURRENT_PULSE: 
+      for( int i = 0; i < numDataPoints; ++i ) {
+        pCurrentArray[i] = pCntxt->neuronParam_I_EXT;
+      }
+      for( int i = 0; i < floor( 150 / integrationStepSize ); ++i ) {
+        pCurrentArray[i] = 0;
+      }
+      for( int i = floor( 200 / integrationStepSize ) ; i < floor( (200 + 2) / integrationStepSize ); ++i ) {
+        pCurrentArray[i] = pCntxt->neuronParam_I_EXT + 5 ;
+      }
+      // printf( "                              +--+                \n" );
+      // printf( "                              |  | i + 5 pA       \n" );
+      // printf( "             i     +----------+  +----------      \n" );
+      // printf( "             0 ----+                              \n" );
+      // printf( "                  150        200 202     t [ms]   \n" );  
+      break;
+
+    case CURRENT_NEGATIVE_STEP: 
+      for( int i = 0; i < numDataPoints; ++i ) {
+        pCurrentArray[i] = - pCntxt->neuronParam_I_EXT;
+      }
+      for( int i = floor( 150 / integrationStepSize ); i < numDataPoints; ++i ) {
+        pCurrentArray[i] = 0;
+      }
+      // printf( "             0     +------------------------      \n" );
+      // printf( "            -i ----+                              \n" );
+      // printf( "                  150                    t [ms]   \n" );  
+      break;
+
+    default:
+      printf( "ERROR: Failed to set up current shape\n" );
+      exitProgram( pCntxt ) ;
+  }  
+
+  return;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+//   HELPER FUNCTIONS TO PROCESS PROGRAM MAIN ARGUMENTS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Check program options
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+void processMainArguments( CNTXT* pCntxt, int argc, char** argv )
+{
+  while( --argc ) {
+
+    if( argv[argc][0] != '-' ) {
+      printf( "ERROR: Parameter error\n");
+      processOPT_H( pCntxt, argc, argv );
+      exitProgram( pCntxt );
+    }
+
+    switch( argv[argc][1] ) {
+
+      case OPT_C:
+        processOPT_C( pCntxt, argc, argv );
+        break;
+
+      case OPT_F:
+        processOPT_F( pCntxt, argc, argv );
+        break;
+
+      case OPT_H:
+        processOPT_H( pCntxt, argc, argv );
+        exitProgram( pCntxt );
+        break;
+
+      case OPT_M:
+        processOPT_M( pCntxt, argc, argv );
+        break;
+
+      case OPT_N:
+        processOPT_N( pCntxt, argc, argv );
+        break;
+
+      default:
+        printf( "ERROR: Parameter error\n");
+        processOPT_H( pCntxt, argc, argv );
+        exitProgram( pCntxt );
+    }
+  }
+
+  // some parameter checks
+
+  if( pCntxt->args.currentShape == 0 ) {
+    printf( "ERROR: option -c not set\n" );
+    processOPT_H( pCntxt, argc, argv );
+    exitProgram( pCntxt );
+  }
+  if( pCntxt->args.outputSvgFilename[0] == NULLCHR ) {
+    printf( "ERROR: option -f not set\n" );
+    processOPT_H( pCntxt, argc, argv );
+    exitProgram( pCntxt );
+  }
+  if( pCntxt->neuronParam_A == 0 ) {
+    printf( "ERROR: option -n not set\n" );
+    processOPT_H( pCntxt, argc, argv );
+    exitProgram( pCntxt );
+  }
+  if(( pCntxt->args.fOdeSolverStandardEuler || pCntxt->args.fOdeSolverSymplecticEuler ||  pCntxt->args.fOdeSolverEulerSpiNNaker || 
+       pCntxt->args.fOdeSolverEulerNEST || pCntxt->args.fOdeSolverIzhikevichNEST ||  pCntxt->args.fOdeSolverHighResolutionEuler ||
+       pCntxt->args.fOdeSolverGSL || pCntxt->args.fOdeSolverHeun || pCntxt->args.fOdeSolverAdaptiveEuler ) == 0 ) {
+    printf( "ERROR: option -m not set\n" );
+    processOPT_H( pCntxt, argc, argv );
+    exitProgram( pCntxt );
+  }
+
+  printf( "+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n" );
+  printf( "+           INZHIKEVICH MODEL SIMULATION                      +\n" );
+  printf( "+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n" );
+  printf( "Following parameters are in use for this run:\n" );
+  printf( "    SVG Output File Name: %s\n", pCntxt->args.outputSvgFilename );
+  printf( "    Neuron Type         : %s\n", pCntxt->neuronTypeName );
+  
+  switch( pCntxt->args.currentShape ) {
+
+    case CURRENT_STEP_150: 
+      printf( "    Current Shape       : step current at 150ms\n" );
+      printf( "                          i     +------------------------      \n" );
+      printf( "                          0 ----+                              \n" );
+      printf( "                               150                    t [ms]   \n" );  
+      break;
+
+    case CURRENT_PULSE: 
+      printf( "    Current Shape       : step current at 150ms with pulse at 200ms for 2ms\n" );
+      printf( "                                           +--+                \n" );
+      printf( "                                           |  | i + 5 pA       \n" );
+      printf( "                          i     +----------+  +----------      \n" );
+      printf( "                          0 ----+                              \n" );
+      printf( "                               150        200 202     t [ms]   \n" );  
+      break;
+
+    case CURRENT_NEGATIVE_STEP: 
+      printf( "    Current Shape       : negative step current at 150ms\n" );
+      printf( "                           0     +------------------------      \n" );
+      printf( "                          -i ----+                              \n" );
+      printf( "                                150                    t [ms]   \n" );  
+      break;
+  }
+
+  printf( "    ODE solver selected : \n" );
+  if( pCntxt->args.fOdeSolverStandardEuler ) {
+      printf( "                          Standard Euler\n" );
+  }
+  if( pCntxt->args.fOdeSolverSymplecticEuler ) {
+      printf( "                          Symplextic Euler\n" );
+  }
+  if( pCntxt->args.fOdeSolverEulerSpiNNaker ) {
+      printf( "                          Explicit Euler SpiNNaker implementation\n" );
+  }
+  if( pCntxt->args.fOdeSolverEulerNEST ) {
+      printf( "                          Explicit Euler NEST implementation\n" );
+  }
+  if( pCntxt->args.fOdeSolverIzhikevichNEST ) {
+      printf( "                          Izhikevich numerics NEST implementation\n" );
+  }
+  if( pCntxt->args.fOdeSolverGSL ) {
+      printf( "                          GSL library\n" );
+  }
+  if( pCntxt->args.fOdeSolverHeun ) {
+      printf( "                          Heun's method\n" );
+  }
+  if( pCntxt->args.fOdeSolverHighResolutionEuler ) {
+      printf( "                          Standard Euler with a high integration step resolution, i.e. 1 micro second\n" );
+  }
+  if( pCntxt->args.fOdeSolverAdaptiveEuler ) {
+      printf( "                          Standard Euler with an adaptive integration step size\n" );
+  }
+
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Set the current shape
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+void processOPT_C( CNTXT* pCntxt, int argc, char** argv ) 
+{
+  if( strncmp( &argv[argc][2], STR_CURRENT_STEP_150, sizeof( STR_CURRENT_STEP_150 )) == 0 ) {
+    pCntxt->args.currentShape = CURRENT_STEP_150;
+  }
+  else if ( strncmp( &argv[argc][2], STR_CURRENT_PULSE, sizeof( STR_CURRENT_PULSE )) == 0 ) {
+    pCntxt->args.currentShape = CURRENT_PULSE;
+  }
+  else if ( strncmp( &argv[argc][2], STR_CURRENT_NEGATIVE_STEP, sizeof( STR_CURRENT_NEGATIVE_STEP )) == 0 ) {
+    pCntxt->args.currentShape = CURRENT_NEGATIVE_STEP;
+  }
+  else {
+    printf( "ERROR: Parameter error, failed to set current shape\n" );
+    exitProgram( pCntxt );
+  }
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Set the svg output file name
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+void processOPT_F( CNTXT* pCntxt, int argc, char** argv )
+{
+  strncpy( pCntxt->args.outputSvgFilename, &argv[argc][2], S_NAME - 1 );
+  pCntxt->args.outputSvgFilename[S_NAME - 1] = NULLCHR;
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Show help text
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+void processOPT_H( CNTXT* pCntxt, int argc, char** argv )
+{
+  printf( "USAGE  : izhikevich -n<neuron type> -c<current shape> -f<svg output path/file.svg> -m<ODE solver method> ... -m<ODE solver method>\n" );
+  printf( "         -f     filename         svg output file name, e.g. svg\\RS_step150.svg\n" );
+  printf( "         -n     RS               regular spiking\n" );
+  printf( "                IB               intrinsically bursting\n" );
+  printf( "                CH               chattering\n" );
+  printf( "                FS               fast spiking\n" );  
+  printf( "                TC               thalamo cortical\n" );
+  printf( "                RZ               resonance\n" );
+  printf( "                LTS              low threshold spiking\n" );
+  printf( "\n" );      
+  printf( "         -c     step150          step current at 150ms\n" );
+  printf( "                                 i     +------------------------      \n" );
+  printf( "                                 0 ----+                              \n" );
+  printf( "                                      150                    t [ms]   \n" );  
+  printf( "                pulse            step current at 150ms with pulse at 200ms for 2ms\n" );
+  printf( "                                                  +--+                \n" );
+  printf( "                                                  |  | i + 5 pA       \n" );
+  printf( "                                 i     +----------+  +----------      \n" );
+  printf( "                                 0 ----+                              \n" );
+  printf( "                                      150        200 202     t [ms]   \n" );    
+  printf( "                nstep            negative step current at 150ms\n" );
+  printf( "                                  0     +------------------------      \n" );
+  printf( "                                 -i ----+                              \n" );
+  printf( "                                      150                    t [ms]   \n" );  
+  printf( "\n" );
+  printf( "         -m     stdEuler         standard Euler\n" );
+  printf( "                sympEuler        symplectic Euler\n" );
+  printf( "                SpiNNaker        explicit Euler SpiNNaker implementation\n" );
+  printf( "                EulerNEST        explicit Euler NEST implementation\n" );  
+  printf( "                IzhikevichNEST   Izhikevich numerics implementation in NEST\n" );  
+  printf( "                GSL              use GSL library\n" );  
+  printf( "                Heun             Heun method\n" );    
+  printf( "                highResEuler     standard Euler with a high resolution integration step size, i.e. 1 micro second\n" );  
+  printf( "                adaptiveEuler    standard Euler with an adaptive integration step size\n" );  
+  printf( "\n" );
+  printf( "EXAMPLE: compare different ODE solver methods for a regular spiking neuron with a positive external step current at 150 milliseconds\n" );
+  printf( "         izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mstdEuler -msympEuler -mSpiNNaker -mEulerNEST -mIzhikevichNEST -mhighResEuler\n" );
+  printf( "\n" );
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Set ODE solver method(s)
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+void processOPT_M( CNTXT* pCntxt, int argc, char** argv )
+{
+  if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_STANDARD_EULER, sizeof( STR_ODE_SOLVER_METHOD_STANDARD_EULER )) == 0 ) {
+    pCntxt->args.fOdeSolverStandardEuler = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_SYMPLECTIC_EULER, sizeof( STR_ODE_SOLVER_METHOD_SYMPLECTIC_EULER )) == 0 ) {
+    pCntxt->args.fOdeSolverSymplecticEuler = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_EULER_SPINNAKER, sizeof( STR_ODE_SOLVER_METHOD_EULER_SPINNAKER )) == 0 ) {
+    pCntxt->args.fOdeSolverEulerSpiNNaker = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_EULER_NEST, sizeof( STR_ODE_SOLVER_METHOD_EULER_NEST )) == 0 ) {
+    pCntxt->args.fOdeSolverEulerNEST = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_IZHIKEVICH_NEST, sizeof( STR_ODE_SOLVER_METHOD_IZHIKEVICH_NEST )) == 0 ) {
+    pCntxt->args.fOdeSolverIzhikevichNEST = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_GSL, sizeof( STR_ODE_SOLVER_METHOD_GSL )) == 0 ) {
+    pCntxt->args.fOdeSolverGSL = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_HEUN, sizeof( STR_ODE_SOLVER_METHOD_HEUN )) == 0 ) {
+    pCntxt->args.fOdeSolverHeun = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_HIGH_RESOLUTION_EULER, sizeof( STR_ODE_SOLVER_METHOD_HIGH_RESOLUTION_EULER )) == 0 ) {
+    pCntxt->args.fOdeSolverHighResolutionEuler = TRUE;
+  }
+  else if( strncmp( &argv[argc][2], STR_ODE_SOLVER_METHOD_ADAPTIVE_EULER, sizeof( STR_ODE_SOLVER_METHOD_ADAPTIVE_EULER )) == 0 ) {
+    pCntxt->args.fOdeSolverAdaptiveEuler = TRUE;
+  }
+  else {
+    printf( "ERROR: Parameter error, invalid ODE solver method\n" );
+    exitProgram( pCntxt );
+  }    
+
+  return;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
+//   Set model parameters depending on the neuron type
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
+void processOPT_N( CNTXT* pCntxt, int argc, char** argv )
+{
+  if( strncmp( &argv[argc][2], STR_NEURON_TYPE_RS, sizeof( STR_NEURON_TYPE_RS )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_RS;
+    strncpy( pCntxt->neuronTypeName, RS_TYPENAME, sizeof( RS_TYPENAME ));
+    pCntxt->neuronParam_A = RS_A;
+    pCntxt->neuronParam_B = RS_B;
+    pCntxt->neuronParam_C = RS_C;
+    pCntxt->neuronParam_D = RS_D;
+    pCntxt->neuronParam_THR = RS_THR;
+    pCntxt->neuronParam_I_EXT = RS_I_EXT;
+    pCntxt->neuronParam_V_INIT = RS_V_INIT;
+    pCntxt->neuronParam_U_INIT = RS_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_IB, sizeof( STR_NEURON_TYPE_IB )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_IB;
+    strncpy( pCntxt->neuronTypeName, IB_TYPENAME, sizeof( IB_TYPENAME ));
+    pCntxt->neuronParam_A = IB_A;
+    pCntxt->neuronParam_B = IB_B;
+    pCntxt->neuronParam_C = IB_C;
+    pCntxt->neuronParam_D = IB_D;
+    pCntxt->neuronParam_THR = IB_THR;
+    pCntxt->neuronParam_I_EXT = IB_I_EXT;
+    pCntxt->neuronParam_V_INIT = IB_V_INIT;
+    pCntxt->neuronParam_U_INIT = IB_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_CH, sizeof( STR_NEURON_TYPE_CH )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_CH;
+    strncpy( pCntxt->neuronTypeName, CH_TYPENAME, sizeof( CH_TYPENAME ));
+    pCntxt->neuronParam_A = CH_A;
+    pCntxt->neuronParam_B = CH_B;
+    pCntxt->neuronParam_C = CH_C;
+    pCntxt->neuronParam_D = CH_D;
+    pCntxt->neuronParam_THR = CH_THR;
+    pCntxt->neuronParam_I_EXT = CH_I_EXT;
+    pCntxt->neuronParam_V_INIT = CH_V_INIT;
+    pCntxt->neuronParam_U_INIT = CH_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_FS, sizeof( STR_NEURON_TYPE_FS )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_FS;
+    strncpy( pCntxt->neuronTypeName, FS_TYPENAME, sizeof( FS_TYPENAME ));
+    pCntxt->neuronParam_A = FS_A;
+    pCntxt->neuronParam_B = FS_B;
+    pCntxt->neuronParam_C = FS_C;
+    pCntxt->neuronParam_D = FS_D;
+    pCntxt->neuronParam_THR = FS_THR;
+    pCntxt->neuronParam_I_EXT = FS_I_EXT;
+    pCntxt->neuronParam_V_INIT = FS_V_INIT;
+    pCntxt->neuronParam_U_INIT = FS_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_TC, sizeof( STR_NEURON_TYPE_TC )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_TC;
+    strncpy( pCntxt->neuronTypeName, TC_TYPENAME, sizeof( TC_TYPENAME ));
+    pCntxt->neuronParam_A = TC_A;
+    pCntxt->neuronParam_B = TC_B;
+    pCntxt->neuronParam_C = TC_C;
+    pCntxt->neuronParam_D = TC_D;
+    pCntxt->neuronParam_THR = TC_THR;
+    pCntxt->neuronParam_I_EXT = TC_I_EXT;
+    pCntxt->neuronParam_V_INIT = TC_V_INIT;
+    pCntxt->neuronParam_U_INIT = TC_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_RZ, sizeof( STR_NEURON_TYPE_RZ )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_RZ;
+    strncpy( pCntxt->neuronTypeName, RZ_TYPENAME, sizeof( RZ_TYPENAME ));
+    pCntxt->neuronParam_A = RZ_A;
+    pCntxt->neuronParam_B = RZ_B;
+    pCntxt->neuronParam_C = RZ_C;
+    pCntxt->neuronParam_D = RZ_D;
+    pCntxt->neuronParam_THR = RZ_THR;
+    pCntxt->neuronParam_I_EXT = RZ_I_EXT;
+    pCntxt->neuronParam_V_INIT = RZ_V_INIT;
+    pCntxt->neuronParam_U_INIT = RZ_U_INIT;
+  }
+  else if( strncmp( &argv[argc][2], STR_NEURON_TYPE_LTS, sizeof( STR_NEURON_TYPE_LTS )) == 0 ) {
+    pCntxt->args.neuronType = NEURON_TYPE_LTS;
+    strncpy( pCntxt->neuronTypeName, LTS_TYPENAME, sizeof( LTS_TYPENAME ));
+    pCntxt->neuronParam_A = LTS_A;
+    pCntxt->neuronParam_B = LTS_B;
+    pCntxt->neuronParam_C = LTS_C;
+    pCntxt->neuronParam_D = LTS_D;
+    pCntxt->neuronParam_THR = LTS_THR;
+    pCntxt->neuronParam_I_EXT = LTS_I_EXT;
+    pCntxt->neuronParam_V_INIT = LTS_V_INIT;
+    pCntxt->neuronParam_U_INIT = LTS_U_INIT;
+  }  
+  else {
+    printf( "ERROR: Failed to set neuron type and neuron type properties\n" );
+    exitProgram( pCntxt );
+  }
+
+  return;
+}

simulation_code/Evaluate_ODE_solver_implementations/src/izhikevich.h

@@ -0,0 +1,152 @@
+/*
+ *  izhikevich.h
+ *
+ *  This file is part of the Izhikevich console application.
+ *
+ *  Copyright (C) 2016, Author: Guido Trensch
+ *
+ *  The Izhikevich console application 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 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  It 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 application. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   TYPE DEFINITIONS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+typedef struct  _CNTXT      CNTXT;
+typedef struct  _MAIN_ARGS  MAIN_ARGS;
+
+typedef float (*PFUNCD)( ... );
+
+typedef enum {
+  STANDARD_EULER                    = 1,
+  SYMPLECTIC_EULER                  = 2,
+  EXPLICIT_EULER_SPINNAKER          = 3,
+  EXPLICIT_EULER_NEST               = 4,
+  EXPLICIT_IZHIKEVICH_NUMERICS_NEST = 5,
+  GSL_LIBRARY                       = 6,
+  HEUN_METHOD                       = 7,
+  ADAPTIVE_EULER                    = 8
+} ENUM_ODESOLVER_METHOD;
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   MACRO DEFINITIONS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define TRUE               1
+#define FALSE              0
+#define S_CNTXT            sizeof( CNTXT )
+#define NULLCHR            '\0'
+#define S_NAME             256
+
+#define PLOT_COLOR_RED     1
+#define PLOT_COLOR_YELLOW  2
+#define PLOT_COLOR_GREEN   3
+#define PLOT_COLOR_GREY    7
+#define PLOT_COLOR_BROWN   8
+#define PLOT_COLOR_BLUE    9  
+#define PLOT_COLOR_CYAN    11
+#define PLOT_COLOR_MAGENTA 13
+
+#define PLOT_LINE_FULL     1
+#define PLOT_LINE_DOTTED   2
+#define PLOT_LINE_DASHED1  3
+#define PLOT_LINE_DASHED2  5
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   COMMAND LINE PARAMETER DEFINITIONS AND DATA STRUCTURES
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// Option -c ... current shape
+#define OPT_C  'c'
+
+typedef enum {
+  CURRENT_STEP_150      = 1,
+  CURRENT_PULSE         = 2,
+  CURRENT_NEGATIVE_STEP = 3
+} ENUM_CURRENT_SHAPE;
+
+#define STR_CURRENT_STEP_150       "step150"
+#define STR_CURRENT_PULSE          "pulse"
+#define STR_CURRENT_NEGATIVE_STEP  "nstep"
+
+// Option -n ... neuron type
+#define OPT_N  'n'
+
+typedef enum {
+  NEURON_TYPE_RS  = 1,       // regular spiking
+  NEURON_TYPE_IB  = 2,       // intrinsically spiking
+  NEURON_TYPE_CH  = 3,       // chattering
+  NEURON_TYPE_FS  = 4,       // fast spiking
+  NEURON_TYPE_TC  = 5,       // thalamo cortical 
+  NEURON_TYPE_RZ  = 6,       // resonance
+  NEURON_TYPE_LTS = 7        // low threshold spiking
+} ENUM_NEURON_TYPE;
+
+#define STR_NEURON_TYPE_RS   "RS"
+#define STR_NEURON_TYPE_IB   "IB"
+#define STR_NEURON_TYPE_CH   "CH"
+#define STR_NEURON_TYPE_FS   "FS"
+#define STR_NEURON_TYPE_TC   "TC"
+#define STR_NEURON_TYPE_RZ   "RZ"
+#define STR_NEURON_TYPE_LTS  "LTS"
+
+// Option -f ... svg output file name
+#define OPT_F  'f'
+
+// Option -m ... ODE solver method
+#define OPT_M  'm'
+
+#define STR_ODE_SOLVER_METHOD_STANDARD_EULER         (char*)"stdEuler"
+#define STR_ODE_SOLVER_METHOD_SYMPLECTIC_EULER       (char*)"sympEuler"
+#define STR_ODE_SOLVER_METHOD_EULER_SPINNAKER        (char*)"SpiNNaker"  
+#define STR_ODE_SOLVER_METHOD_EULER_NEST             (char*)"EulerNEST"
+#define STR_ODE_SOLVER_METHOD_IZHIKEVICH_NEST        (char*)"IzhikevichNEST"
+#define STR_ODE_SOLVER_METHOD_GSL                    (char*)"GSL"
+#define STR_ODE_SOLVER_METHOD_HEUN                   (char*)"Heun"
+#define STR_ODE_SOLVER_METHOD_HIGH_RESOLUTION_EULER  (char*)"highResEuler"
+#define STR_ODE_SOLVER_METHOD_ADAPTIVE_EULER         (char*)"adaptiveEuler"
+
+// Option -h ... help
+#define OPT_H  'h'
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   DATA STRUCTURES
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+struct _MAIN_ARGS {
+  ENUM_CURRENT_SHAPE  currentShape;
+  ENUM_NEURON_TYPE    neuronType;
+  char                outputSvgFilename[S_NAME];
+  bool                fOdeSolverStandardEuler;
+  bool                fOdeSolverSymplecticEuler;  
+  bool                fOdeSolverEulerSpiNNaker;
+  bool                fOdeSolverEulerNEST;
+  bool                fOdeSolverIzhikevichNEST;
+  bool                fOdeSolverGSL;
+  bool                fOdeSolverHeun;
+  bool                fOdeSolverHighResolutionEuler;
+  bool                fOdeSolverAdaptiveEuler;
+};
+
+struct _CNTXT {
+  MAIN_ARGS  args;
+  char      neuronTypeName[S_NAME];
+  float     neuronParam_A;
+  float     neuronParam_B;
+  float     neuronParam_C;
+  float     neuronParam_D;
+  float     neuronParam_THR;
+  float     neuronParam_I_EXT;
+  float     neuronParam_V_INIT;
+  float     neuronParam_U_INIT;
+  float     legendYOffset;
+};
+

simulation_code/Evaluate_ODE_solver_implementations/src/izhikevich_model.h

@@ -0,0 +1,151 @@
+/*
+ *  izhikevich_model.h
+ *
+ *  This file is part of the Izhikevich console application.
+ *
+ *  Copyright (C) 2016, Author: Guido Trensch
+ *
+ *  The Izhikevich console application 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 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  It 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 application. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   IZHIKEVICH MODEL DEFINITION
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// IZHIKEVICH MODEL CONSTANTS SET UP AS GLOBAL VARIABLES
+float IZHIKEVICH_A   = 0;
+float IZHIKEVICH_B   = 0;
+float IZHIKEVICH_C   = 0;
+float IZHIKEVICH_D   = 0;
+float IZHIKEVICH_THR = 0;
+
+#define TRUE  1
+#define FALSE 0
+
+// ODE 1
+float dvdt( float  v         // IN
+          , float  u         // IN
+          , void*  pI )      // IN  external current (optional for the ODE solver)
+{
+  if( pI == NULL ) {
+    printf( "ERROR: external current expected\n" );
+    exit( -1 );
+  } 
+  // SpiNNaker FP error simulation 
+  // return( 0.03997802 * v * v + 5.0 * v + 140.0 - u + *(float*)pI );
+  return( 0.04 * v * v + 5.0 * v + 140.0 - u + *(float*)pI );
+}
+
+// ODE 2
+float dudt( float  v         // IN
+          , float  u )       // IN
+{
+  return( IZHIKEVICH_A * ( IZHIKEVICH_B * v - u ));
+}
+
+// Threshold
+int threshold( float*  v     // IN / OUT
+              , float*  u )   // IN / OUT             
+{
+  if( *v >= IZHIKEVICH_THR ) {
+    *v = IZHIKEVICH_C;
+    *u += IZHIKEVICH_D;
+    return( TRUE );
+  }
+  else {
+    return( FALSE );    
+  }
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   IZHIKEVICH MODEL PARAMETERS
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// parameter set for regular spiking
+#define RS_A        (float)(0.02)   
+#define RS_B        (float)(0.2)
+#define RS_C        (float)(-65.0)
+#define RS_D        (float)(8.0)
+#define RS_THR      (float)(30.0)
+#define RS_I_EXT    (float)(5.0)
+#define RS_V_INIT   (float)(-75.0)   
+#define RS_U_INIT   (float)(0.0)
+#define RS_TYPENAME "IZHIKEVICH RS NEURON TYPE"
+
+// parameter set for intrinsically bursting
+#define IB_A        (float)(0.02)   
+#define IB_B        (float)(0.2)
+#define IB_C        (float)(-55.0)
+#define IB_D        (float)(4.0)
+#define IB_THR      (float)(30.0)
+#define IB_I_EXT    (float)(5.0)
+#define IB_V_INIT   (float)(-75.0)   
+#define IB_U_INIT   (float)(0.0)   
+#define IB_TYPENAME "IZHIKEVICH IB NEURON TYPE"
+
+// parameter set for chattering
+#define CH_A        (float)(0.02)   
+#define CH_B        (float)(0.2)
+#define CH_C        (float)(-50.0)
+#define CH_D        (float)(2.0)
+#define CH_THR      (float)(30.0)
+#define CH_I_EXT    (float)(5.0)
+#define CH_V_INIT   (float)(-75.0)   
+#define CH_U_INIT   (float)(0.0)
+#define CH_TYPENAME "IZHIKEVICH CH NEURON TYPE"
+
+// parameter set for fast spiking
+#define FS_A        (float)(0.1)   
+#define FS_B        (float)(0.2)
+#define FS_C        (float)(-65.0)
+#define FS_D        (float)(2.0)
+#define FS_THR      (float)(30.0)
+#define FS_I_EXT    (float)(5.0)
+#define FS_V_INIT   (float)(-75.0)   
+#define FS_U_INIT   (float)(0.0)
+#define FS_TYPENAME "IZHIKEVICH FS NEURON TYPE"
+
+// parameter set for thalamo cortical
+#define TC_A        (float)(0.02)   
+#define TC_B        (float)(0.25)
+#define TC_C        (float)(-65.0)
+#define TC_D        (float)(0.05)
+#define TC_THR      (float)(30.0)
+#define TC_I_EXT    (float)(5.0)
+#define TC_V_INIT   (float)(-75.0)   
+#define TC_U_INIT   (float)(0.0)
+#define TC_TYPENAME "IZHIKEVICH TC NEURON TYPE"
+
+// parameter set for resonator
+// requires extra current step at e.g. 150 ms
+#define RZ_A        (float)(0.1)   
+#define RZ_B        (float)(0.25)
+#define RZ_C        (float)(-65.0)
+#define RZ_D        (float)(2.0)
+#define RZ_THR      (float)(30.0)
+#define RZ_I_EXT    (float)(0.7)
+#define RZ_V_INIT   (float)(-75.0)   
+#define RZ_U_INIT   (float)(0.0)
+#define RZ_TYPENAME "IZHIKEVICH RZ NEURON TYPE"
+
+// parameter set for low-threshold spiking
+#define LTS_A        (float)(0.02)   
+#define LTS_B        (float)(0.25)
+#define LTS_C        (float)(-65.0)
+#define LTS_D        (float)(2.0)
+#define LTS_THR      (float)(30.0)
+#define LTS_I_EXT    (float)(5.0)
+#define LTS_V_INIT   (float)(-75.0)   
+#define LTS_U_INIT   (float)(0.0)
+#define LTS_TYPENAME "IZHIKEVICH LTS NEURON TYPE"

simulation_code/Evaluate_ODE_solver_implementations/src/ode_solvers.h

@@ -0,0 +1,384 @@
+/*
+ *  ode_solver.h
+ *
+ *  This file is part of the Izhikevich console application.
+ *
+ *  Copyright (C) 2016, Author: Guido Trensch
+ *
+ *  The Izhikevich console application 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 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  It 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 application. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_odeiv.h>
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   2D ODE SOLVER: STANDARD EULER
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ode2dSolver_StandardEuler( PFUNCD  pODE1dxdt          // IN        function pointers to the 2d coupled ODE system
+                              , PFUNCD  pODE2dydt          // IN 
+                              , float   h                  // IN        step width
+                              , float   x_t                // IN        x(t)
+                              , float   y_t                // IN        y(t)
+                              , float*  x_tplus1           // OUT       x(t+1)
+                              , float*  y_tplus1           // OUT       y(t+1)
+                              , void*   optParm = NULL )   // IN / OUT  Optional parameters
+{
+  // NOTE: There is a problem with the function pointer return values.
+  //       It surprisingly works if all data types are double instead of float !!!
+  //       To circumvent the problem the functions are called directly.
+
+  // *x_tplus1 = x_t + h * pODE1dxdt( x_t, y_t, optParm );
+  *x_tplus1 = x_t + h * dvdt( x_t, y_t, optParm );
+  // *y_tplus1 = y_t + h * pODE2dydt( x_t, y_t, optParm );
+  *y_tplus1 = y_t + h * dudt( x_t, y_t );
+
+  return;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   2D ODE SOLVER: SYMPLECTIC EULER
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ode2dSolver_SymplecticEuler( PFUNCD  pODE1dxdt          // IN        function pointers to the 2d coupled ODE system
+                                , PFUNCD  pODE2dydt          // IN 
+                                , float   h                  // IN        step width
+                                , float   x_t                // IN        x(t)
+                                , float   y_t                // IN        y(t)
+                                , float*  x_tplus1           // OUT       x(t+1)
+                                , float*  y_tplus1           // OUT       y(t+1)
+                                , void*   optParm = NULL )   // IN / OUT  Optional parameters
+{
+  // NOTE: There is a problem with the function pointer return values.
+  //       It surprisingly works if all data types are double instead of float !!!
+  //       To circumvent the problem the functions are called directly.
+
+  // *x_tplus1 = x_t + h * pODE1dxdt( x_t, y_t, optParm );
+  *x_tplus1 = x_t + h * dvdt( x_t, y_t, optParm );
+  // *y_tplus1 = y_t + h * pODE2dydt( *x_tplus1, *y_tplus1, optParm );    // symplectic
+  *y_tplus1 = y_t + h * dudt( *x_tplus1, *y_tplus1 );    // symplectic
+
+  return;
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   EXPLIZIT SOLVER FOR THE IZHIKEVICH ODEs
+// =   SpiNNaker implementation
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void explizitSpiNNaker( float   h
+                      , float*  v
+                      , float*  u
+                      , float   i )
+{  
+  float lastV1 = *v;
+  float lastU1 = *u;
+  float a = IZHIKEVICH_A;
+  float b = IZHIKEVICH_B;
+  float input_this_timestep = i;
+
+  float pre_alph = (140.0) + input_this_timestep - lastU1;
+  float alpha    = pre_alph + ( (5.0) + (0.0400) * lastV1) * lastV1;
+  float eta = lastV1 + 0.5*(h * alpha);
+
+  float beta = 0.5*(h * (b * lastV1 - lastU1) * a); 
+
+  *v += h * (pre_alph - beta + ( (5.0) + (0.0400) * eta) * eta);
+
+  *u += a * h * (-lastU1 - beta + b * eta);         
+
+  threshold( v, u );  
+  return;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   EXPLIZIT SOLVER FOR THE IZHIKEVICH ODEs
+// =   explizit euler ( NEST implementation 1 ) 
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void explizitNEST1( float   h
+                  , float*  v
+                  , float*  u
+                  , float   i )
+{  
+  float v_old = *v;
+  float u_old = *u;
+    
+  float a = IZHIKEVICH_A;
+  float b = IZHIKEVICH_B;
+
+  *v += h * ( 0.04 * v_old * v_old + 5.0 * v_old + 140.0 - u_old + i );  
+  *u += h * a * ( b * v_old - u_old );
+
+  threshold( v, u );  
+  return;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   EXPLIZIT SOLVER FOR THE IZHIKEVICH ODEs
+// =   Izhikevich numerics from 2003 paper ( NEST implementation 2 ) 
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void explizitNEST2( float   h
+                  , float*  v
+                  , float*  u
+                  , float   i )
+{
+    float v_old = *v;
+    float u_old = *u;
+    
+    float a = IZHIKEVICH_A;
+    float b = IZHIKEVICH_B;
+
+    *v += h * 0.5 * ( 0.04 * v_old * v_old + 5.0 * v_old + 140.0 - u_old + i );  // for numetrical stablility
+    *v += h * 0.5 * ( 0.04 * v_old * v_old + 5.0 * v_old + 140.0 - u_old + i );  // see Izhikevich 2003
+    *u += h * a * ( b * v_old - u_old );
+
+    threshold( v, u );
+    return;
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   GNU Scientific Library
+// =   Based on the code generated by NESTML
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+int gslODESolver_EquationDynamics( double t, const double y[], double f[], void* params );
+
+const int V_m_INDEX = 0;
+const int U_m_INDEX = 1;
+
+float param_I_e = 0;
+
+static int countGSL = 0;
+
+// gsl_odeiv.h
+// typedef struct  {
+//   int (* function) (float t, const float y[], float dydt[], void * params);
+//   int (* jacobian) (float t, const float y[], float * dfdy, float dfdt[], void * params);
+//   size_t dimension;
+//   void * params;
+// }
+// gsl_odeiv_system;
+gsl_odeiv_system    systemOfODEs = { 0 };
+gsl_odeiv_step*     pSteppingFunction  = NULL;
+gsl_odeiv_control*  pControlObject     = NULL;
+gsl_odeiv_evolve*   pEvolutionFunction = NULL;
+
+
+void gslODESolver_Integrate( float   h
+                           , float*  v
+                           , float*  u
+                           , float   i )
+{
+  int    gslRc = GSL_SUCCESS;
+  double t = 0;
+  double t1 = h;
+  double stateVector[2];
+  double stepSize = h;          // stepSize will be adjusted by the evolve-function
+
+  param_I_e = i;
+
+  stateVector[V_m_INDEX] = *v;
+  stateVector[U_m_INDEX] = *u;
+  
+
+
+  while( t < t1 ) {                                        // There are 3 loops !!!
+                                                           // - The outer loop of the simulation.
+                                                           // - This loop for t < stepSize < t1. Where stepSize is adjusted each iteration.
+                                                           // - The loop within evolve calling the ODEs "systemOfODEs.function = &gslODESolver_EquationDynamics"
+    gslRc = gsl_odeiv_evolve_apply( pEvolutionFunction
+                                  , pControlObject         // accuracy
+                                  , pSteppingFunction      // the algorithm
+                                  , &systemOfODEs          // equation
+                                  , &t
+                                  , t1
+                                  , &stepSize
+                                  , stateVector );         // needs to be updated in the equations, it is reset in between
+                                                           // It won't work to pass parameters, like synapse currents, to the equations
+                                                           // using the 
+
+    if( gslRc != GSL_SUCCESS ) {
+      printf( "ERROR: gsl error\n" );
+      exit( -1 );
+    }
+
+    // Threshold
+    if( stateVector[V_m_INDEX] >= IZHIKEVICH_THR ) {
+      stateVector[V_m_INDEX] = IZHIKEVICH_C;
+      stateVector[U_m_INDEX] += IZHIKEVICH_D;
+    }
+
+  }
+
+  *v = stateVector[V_m_INDEX];
+  *u = stateVector[U_m_INDEX];
+
+  return;
+}
+
+void gslODESolver_Init()
+{
+  pSteppingFunction = gsl_odeiv_step_alloc( gsl_odeiv_step_rkf45, 2 );      // instance of a stepping function for a two dimensinal system
+  pControlObject = gsl_odeiv_control_y_new( 1e-6, 0.0 );                    // create a new control object which will keep the local error on each step 
+                                                                            // within an absolute error 1e-6 and relative error of 0.0 with respect to 
+                                                                            // the solution y_i(t)
+  pEvolutionFunction = gsl_odeiv_evolve_alloc ( 2 );                        // instance of an evolution function for a two dimensinal system
+  
+  systemOfODEs.function = &gslODESolver_EquationDynamics;
+  systemOfODEs.jacobian = NULL;
+  systemOfODEs.dimension = 2;
+  systemOfODEs.params = &param_I_e;
+ 
+  return;
+}
+
+int gslODESolver_EquationDynamics( double t, const double y[], double f[], void* params )
+{
+  float I_e = *(float*)params;
+
+  f[ V_m_INDEX ] = 0.04 * y[V_m_INDEX] * y[V_m_INDEX] + 5.0 * y[V_m_INDEX] + 140 - y[U_m_INDEX] + I_e;
+  f[ U_m_INDEX ] = IZHIKEVICH_A * ( IZHIKEVICH_B * y[V_m_INDEX] - y[U_m_INDEX] );
+
+  countGSL++;
+  printf( "countGSL: %d\n", countGSL );
+
+  return GSL_SUCCESS;
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   Heun
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ode2dSolver_Heun( PFUNCD  pODE1dxdt          // IN        function pointers to the 2d coupled ODE system
+                     , PFUNCD  pODE2dydt          // IN 
+                     , float   h                  // IN        step width
+                     , float   x_t                // IN        x(t)
+                     , float   y_t                // IN        y(t)
+                     , float*  x_tplus1           // OUT       x(t+1)
+                     , float*  y_tplus1           // OUT       y(t+1)
+                     , void*   optParm = NULL )   // IN / OUT  Optional parameters
+{
+  // NOTE: There is a problem with the function pointer return values.
+  //       It surprisingly works if all data types are double instead of float !!!
+  //       To circumvent the problem the functions are called directly.
+
+  // float x_Euler = x_t + h * pODE1dxdt( x_t, y_t, optParm );
+  float x_Euler = x_t + h * dvdt( x_t, y_t, optParm );
+  // float y_Euler = y_t + h * pODE2dydt( x_t, y_t, optParm );
+  float y_Euler = y_t + h * dudt( x_t, y_t );
+  // *x_tplus1 = x_t + h * pODE1dxdt( 0.5 * (x_t + x_Euler), 0.5 * (y_t + y_Euler), optParm );
+  *x_tplus1 = x_t + h * dvdt( 0.5 * (x_t + x_Euler), 0.5 * (y_t + y_Euler), optParm );
+  // *y_tplus1 = y_t + h * pODE2dydt( 0.5 * (x_t + x_Euler), 0.5 * (y_t + y_Euler), optParm );
+  *y_tplus1 = y_t + h * dudt( 0.5 * (x_t + x_Euler), 0.5 * (y_t + y_Euler) );
+
+  return;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   Standard Euler with adaptive stepsize
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ode2dSolver_AdaptiveEuler( PFUNCD  pODE1dxdt          // IN        function pointers to the 2d coupled ODE system
+                              , PFUNCD  pODE2dydt          // IN 
+                              , float   h                  // IN        step width
+                              , float   x_t                // IN        x(t)
+                              , float   y_t                // IN        y(t)
+                              , float*  x_tplus1           // OUT       x(t+1)
+                              , float*  y_tplus1           // OUT       y(t+1)
+                              , void*   optParm = NULL )   // IN / OUT  Optional parameters
+{
+  float t        = 0;
+  float stepSize = h;   // initial stepsize, will be adjusted in every step
+  float x_interm = x_t;
+  float y_interm = y_t;
+
+  // determine integration step size
+  float error_accepted   = 0.2;
+  float error_calculated = 0;
+  int   k = 1;
+  float x_hi;
+
+  static float simulationTime_ms = 0;
+  int count1 = 0;
+  int count2 = 0;
+  static int count3 = 0;
+  static int count4 = 0;
+
+  float tempError[32] = { 0.0 };
+  int index = 0;
+  int i = 0;
+  int fSpike = FALSE;
+  
+  if( count4 == 0 ) {
+    printf( "[INITIAL VALUES]      h: %11.6f    x_t: %11.6f    y_t: %11.6f   \n", h, x_t, y_t );
+  }
+  
+  do {
+    stepSize = h/k;
+   
+    x_hi = x_t + stepSize * dvdt( x_t, y_t, optParm ); 
+
+    error_calculated =  fabs(fabs(x_hi) - fabs(x_t));
+
+    if( index == 31 ) {
+      tempError[ index ] = 99.99;
+    }
+    else {
+      tempError[ index ] = error_calculated;
+      index++;
+    }
+
+
+    k *= 2;
+    count1++;
+  }  while(!(( error_calculated < error_accepted ) || ( stepSize <= 0.001 )));
+
+  // if( stepSize < 0.001 ) stepSize = 0.001;
+
+  while( t < h ) {
+
+    x_interm = x_interm + stepSize * dvdt( x_interm, y_interm, optParm );
+    y_interm = y_interm + stepSize * dudt( x_interm, y_interm );
+    fSpike = threshold( &x_interm, &y_interm );
+
+    t += stepSize;
+    count2++;
+    count3++;
+
+    if( fSpike ) {
+      printf( "[%8.3f]   > > > S P I K E < < <   +  %11.6f ms \n", simulationTime_ms, t );
+      // break;
+    }
+
+//    t = roundf( t * 1000 ) / 1000;
+  }; 
+
+  #if( 0 )
+    printf( "[%8.3f]   V: %11.6f   U: %11.6f   stepSize: %11.6f   CNT(find stepSize): %6d   CNT(calc): %6d   CNT(total): %6d "
+          , simulationTime_ms, x_interm, y_interm, stepSize, count1, count2, count3 );
+    printf( "      ERRORs: " );
+    for( i = 0; i < index; ++i ) {
+      printf( "%11.6f, ", tempError[ i ] );
+    }
+    printf( "\n" );
+  #endif
+
+  *x_tplus1 = x_interm;
+  *y_tplus1 = y_interm;
+
+  count4++;
+  simulationTime_ms += h;
+  return;
+}

simulation_code/Evaluate_ODE_solver_implementations/src/run.sh

@@ -0,0 +1,30 @@
+#!/bin/sh
+
+cmake CMakeLists.txt 
+
+make clean
+make
+
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mstdEuler -msympEuler -mSpiNNaker -mEulerNEST -mIzhikevichNEST -mGSL -mHeun -mhighResEuler -madaptiveEuler
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -mstdEuler
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -msympEuler
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -mSpiNNaker
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -mEulerNEST
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -mIzhikevichNEST
+./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mSpiNNaker -mGSL
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mhighResEuler -madaptiveEuler
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mSpiNNaker
+#./izhikevich -nCH -cstep150 -fsvg/RS_step150.svg -mhighResEuler -madaptiveEuler -mGSL
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mGSL -madaptiveEuler -mhighResEuler
+
+#./izhikevich -nRS -cstep150 -fsvg/RS_GSL.svg -mGSL
+#./izhikevich -nRS -cstep150 -fsvg/RS_SpiNNaker.svg -mSpiNNaker
+#./izhikevich -nCH -cstep150 -fsvg/CH_Type_IntStepSize0.5ms.svg -mstdEuler -mGSL -madaptiveEuler
+#./izhikevich -nRS -cstep150 -fsvg/RS_step150.svg -mstdEuler
+
+#./izhikevich -nIB -cstep150 -fsvg/IB_step150.svg
+#./izhikevich -nCH -cstep150 -fsvg/CH_step150.svg
+#./izhikevich -nFS -cstep150 -fsvg/FS_step150.svg
+#./izhikevich -nTC -cstep150 -fsvg/TC_step150.svg
+#./izhikevich -nRZ -cstep150 -fsvg/RZ_step150.svg
+#./izhikevich -nLTS -cstep150 -fsvg/LTS_step150.svg

simulation_code/LICENSE

@@ -0,0 +1,339 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Lesser General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+                            NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    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 2 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, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) year name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.

simulation_code/README.md

@@ -0,0 +1,7 @@
+These are the C/C++ source files and PyNN scripts developed as part of a study to reproduce selected network states of the Izhikevich polychronization model on the SpiNNaker neuromorphic system.
+
+**STATUS**: Several releases for this source code have been published on Zenodo: 
+
+Ver. 1.0: [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1435026.svg)](https://doi.org/10.5281/zenodo.1435026)
+
+Ver. 2.0  [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1435831.svg)](https://doi.org/10.5281/zenodo.1435831)

simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.c

@@ -2,10 +2,10 @@
 // = Set up desired ODE solver implementation
 // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 
-#define IMPL01  false    // original SpiNNaker ESR implementation
+#define IMPL01  false    // original SpiNNaker ESR implementation (corresponds to Iteration I)
 #define IMPL02  false    // original Izhikevich implementation
-#define IMPL03  false    // fixed step size (h=1/16) forward Euler, precise threshold detection
-#define IMPL04  true     // fixed step size (h=1/16) forward Euler, precise threshold detection, FP conversion
+#define IMPL03  false    // fixed-step size (h=1/16) symplectic forward Euler, precise threshold detection
+#define IMPL04  true     // fixed-step size (h=1/16) symplectic forward Euler, precise threshold detection, FP conversion (corresponds to Iteration III)
 
 // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 
@@ -48,7 +48,7 @@ static inline void _rk2_kernel_midpoint(REAL h, neuron_pointer_t neuron,
 #endif
 
 #if( IMPL02 )
-// h = 1.0, i.e., simulation timestep must be set to 1.0 !
+// h = 1.0, i.e., simulation time-step must be set to 1.0 !
 static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron,
                                         REAL input_this_timestep) {
    REAL v = neuron->V;
@@ -66,7 +66,7 @@ static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron,
 #endif
 
 #if( IMPL03 )
-// h = 1.0, i.e., simulation timestep must be set to 1.0 !
+// h = 1.0, i.e., simulation time-step must be set to 1.0 !
 static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
                                         REAL input_this_timestep) {
    #define THRESHOLD_DETECTED  1024.0k      // some high value above threshold to trigger spike event
@@ -108,7 +108,7 @@ static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
 #endif
 
 #if( IMPL04 )
-// h = 1.0, i.e., simulation timestep must be set to 1.0 !
+// h = 1.0, i.e., simulation time-step must be set to 1.0 !
 static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuron,
                                                   REAL input_this_timestep) {
    #define THRESHOLD_DETECTED  1024.0k      // some high value above threshold to trigger spike event