Browse Source

add simulation codes from github.com/gtrensch/RigorousNeuralNetworkSimulations

Robin Gutzen 5 years ago
parent
commit
3df4f1095b
30 changed files with 134139 additions and 0 deletions
  1. 1000 0
      simulation_code/C_model/data/conMatrix.dat
  2. 1000 0
      simulation_code/C_model/data/delayMatrix.dat
  3. 12508 0
      simulation_code/C_model/data/pythonConWeightDelay_after1h.py
  4. 12508 0
      simulation_code/C_model/data/pythonConWeightDelay_after2h.py
  5. 12508 0
      simulation_code/C_model/data/pythonConWeightDelay_after3h.py
  6. 12508 0
      simulation_code/C_model/data/pythonConWeightDelay_after4h.py
  7. 12508 0
      simulation_code/C_model/data/pythonConWeightDelay_after5h.py
  8. 60946 0
      simulation_code/C_model/data/randomNetworkInput.dat
  9. 1000 0
      simulation_code/C_model/data/weightMatrix_after1h.dat
  10. 1000 0
      simulation_code/C_model/data/weightMatrix_after2h.dat
  11. 1000 0
      simulation_code/C_model/data/weightMatrix_after3h.dat
  12. 1000 0
      simulation_code/C_model/data/weightMatrix_after4h.dat
  13. 1000 0
      simulation_code/C_model/data/weightMatrix_after5h.dat
  14. 4 0
      simulation_code/C_model/src/CMakeLists.txt
  15. 126 0
      simulation_code/C_model/src/globals.h
  16. 671 0
      simulation_code/C_model/src/main.cpp
  17. 58 0
      simulation_code/C_model/src/params.h
  18. 60 0
      simulation_code/C_model/src/params_create_5_selected_states.h
  19. 54 0
      simulation_code/C_model/src/params_simulate_selected_state_1.h
  20. 54 0
      simulation_code/C_model/src/params_simulate_selected_state_2.h
  21. 54 0
      simulation_code/C_model/src/params_simulate_selected_state_3.h
  22. 54 0
      simulation_code/C_model/src/params_simulate_selected_state_4.h
  23. 54 0
      simulation_code/C_model/src/params_simulate_selected_state_5.h
  24. 1077 0
      simulation_code/C_model/src/poly_spnet.cpp
  25. 627 0
      simulation_code/C_model/src/utils.cpp
  26. 61 0
      simulation_code/C_model/src/utils.h
  27. 251 0
      simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.c
  28. 30 0
      simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.h
  29. 299 0
      simulation_code/SpiNNaker_model/PyNN/polychronizationNetwork.py
  30. 119 0
      simulation_code/SpiNNaker_model/PyNN/singleNeuronDynamics.py

File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/conMatrix.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/delayMatrix.dat


File diff suppressed because it is too large
+ 12508 - 0
simulation_code/C_model/data/pythonConWeightDelay_after1h.py


File diff suppressed because it is too large
+ 12508 - 0
simulation_code/C_model/data/pythonConWeightDelay_after2h.py


File diff suppressed because it is too large
+ 12508 - 0
simulation_code/C_model/data/pythonConWeightDelay_after3h.py


File diff suppressed because it is too large
+ 12508 - 0
simulation_code/C_model/data/pythonConWeightDelay_after4h.py


File diff suppressed because it is too large
+ 12508 - 0
simulation_code/C_model/data/pythonConWeightDelay_after5h.py


File diff suppressed because it is too large
+ 60946 - 0
simulation_code/C_model/data/randomNetworkInput.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/weightMatrix_after1h.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/weightMatrix_after2h.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/weightMatrix_after3h.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/weightMatrix_after4h.dat


File diff suppressed because it is too large
+ 1000 - 0
simulation_code/C_model/data/weightMatrix_after5h.dat


+ 4 - 0
simulation_code/C_model/src/CMakeLists.txt

@@ -0,0 +1,4 @@
+cmake_minimum_required(VERSION 2.8)
+set(CMAKE_CXX_STANDARD 11)
+project(RigorousNeuralNetworkSimulations)
+add_executable(c_model main.cpp utils.cpp)

+ 126 - 0
simulation_code/C_model/src/globals.h

@@ -0,0 +1,126 @@
+#ifndef __GLOBALS_H__
+#define __GLOBALS_H__
+/*
+ *  globals.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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 "params.h"
+
+#if( __RUN_REFACTORED_VERSION__ )
+
+#ifdef __IMPORT_GLOBAL_VARIABLES__
+  #define EXTERN extern
+#else
+  #define EXTERN
+#endif
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   N E T W O R K   P A R A M E T E R S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define NUM_EXCITATORY_NEURONS      (int)(800)
+#define NUM_INHIBITORY_NEURONS      (int)(200)
+#define NUM_SYNAPSES_PER_NEURON     (int)(100)
+#define MAX_SYNAPSE_DELAY           (int)(20)
+#define MAX_SYNAPTIC_STRENGTH       (double)(10.0)
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   M A C R O   D E F I N I T I O N S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define GET_RANDOM_INT(max)         (int(floor(max * ( 0.9999999 * double(rand()) / RAND_MAX))))
+
+#define NUM_TOTAL_NEURONS           (int)(NUM_EXCITATORY_NEURONS + NUM_INHIBITORY_NEURONS)
+#define MAX_NUM_FIRINGS             (int)(100 * NUM_TOTAL_NEURONS)
+#define INIT_EXC_SYNAPTIC_WEIGHT    (double)(6.0)
+#define INIT_INH_SYNAPTIC_WEIGHT    (double)(-5.0)
+
+#define RS_A                        (double)(0.02)
+#define RS_D                        (double)(8.0)
+#define RS_V_INIT                   (double)(-65.0)
+#define RS_U_INIT                   (double)(0.2 * RS_V_INIT)
+
+#define FS_A                        (double)(0.1)
+#define FS_D                        (double)(2.0)
+#define FS_V_INIT                   (double)(-65.0)
+#define FS_U_INIT                   (double)(0.2 * FS_V_INIT)
+
+#define RS_FS_B                     (double)(0.2)
+#define RS_FS_C                     (double)(-65.0)
+
+#define RS_FS_THR                   (double)(30.0)
+
+#define I_EXT                       (double)(20.0)
+
+#define RC_NOID                     (int)(-1)
+#define RC_ERROR_EXIT               (int)(-2)
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   G L O B A L   D A T A   S T R U C T U R E S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+EXTERN int      matrixOfPostSynapticNeurons        [ NUM_TOTAL_NEURONS ][ NUM_SYNAPSES_PER_NEURON ];
+EXTERN double   matrixOfSynapticWeights            [ NUM_TOTAL_NEURONS ][ NUM_SYNAPSES_PER_NEURON ];
+EXTERN double   matrixOfSynapticWeights_derivatives[ NUM_TOTAL_NEURONS ][ NUM_SYNAPSES_PER_NEURON ];
+
+EXTERN short    numEntriesPerDelay[NUM_TOTAL_NEURONS][MAX_SYNAPSE_DELAY];
+EXTERN short    listOfSynsByNeuronAndDelay[NUM_TOTAL_NEURONS][MAX_SYNAPSE_DELAY][NUM_SYNAPSES_PER_NEURON];
+
+EXTERN int      numPreSynapticNeuronsOfTarget[ NUM_TOTAL_NEURONS ];
+
+// REFACTOR COMMENT: assumption that it fits
+EXTERN int      listOfPresynapticNeurons[ NUM_TOTAL_NEURONS ][ 3 * NUM_SYNAPSES_PER_NEURON ];
+EXTERN int      listOfPresynapticDelays [ NUM_TOTAL_NEURONS ][ 3 * NUM_SYNAPSES_PER_NEURON ];
+
+// REFACTOR COMMENT: dead code removed
+// double*  listOfPointersToSynapticWeights  [NUM_TOTAL_NEURONS][3 * NUM_SYNAPSES_PER_NEURON];
+
+EXTERN double*  listOfPointersToSynapticWeights_derivatives[NUM_TOTAL_NEURONS][3 * NUM_SYNAPSES_PER_NEURON];
+
+EXTERN double	  LTP[NUM_TOTAL_NEURONS][1001 + MAX_SYNAPSE_DELAY];
+EXTERN double   LTD[NUM_TOTAL_NEURONS];
+
+EXTERN int      numFirings;
+EXTERN int      firings[MAX_NUM_FIRINGS][2];
+#define TIME    (int)(0)
+#define NEURON  (int)(1)
+
+EXTERN double   I_ext[NUM_TOTAL_NEURONS];
+
+EXTERN double   a[NUM_TOTAL_NEURONS];                                                // neuronal dynamics parameters
+EXTERN double   d[NUM_TOTAL_NEURONS];
+EXTERN double   v[NUM_TOTAL_NEURONS];                                                // activity variables
+EXTERN double   u[NUM_TOTAL_NEURONS];
+
+#if( ODE_SOLVER_REFINEMENT )
+EXTERN int      spike[NUM_TOTAL_NEURONS];                                            // spike count in an interval  (precise threshold detection)
+#endif
+
+#if( LOG_MIN_MAX_V_U )
+EXTERN double   v_max;
+EXTERN double   v_min;
+EXTERN double   u_max;
+EXTERN double   u_min;
+#endif
+
+EXTERN FILE*    pFileStimulusOutput;
+EXTERN FILE*    pFileStimulusInput;
+
+#endif   // __RUN_REFACTORED_VERSION__
+#endif   // __GLOBALS_H__

+ 671 - 0
simulation_code/C_model/src/main.cpp

@@ -0,0 +1,671 @@
+/*
+ *  main.cpp
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  This source code is based on the poly_spnet.cpp and spnet.cpp source code available at
+ *  https://www.izhikevich.org/publications/spnet.htm.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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 "params.h"
+
+#if( __RUN_REFACTORED_VERSION__ )
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <cstring>
+
+#include "globals.h"
+#include "utils.h"
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   C H E C K   P A R A M E T E R   S E T T I N G S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#ifndef SIM_TIME
+  #error SIM_TIME not defined!
+#endif
+
+#ifdef GENERATE_NETWORK_FROM_EXTERNAL_DATA
+  #if( GENERATE_NETWORK_FROM_EXTERNAL_DATA )
+    #ifndef INFILE_CONNECTION_MATRIX
+      #error INFILE_CONNECTION_MATRIX not defined!
+    #endif
+    #ifndef INFILE_DELAY_MATRIX
+      #error INFILE_DELAY_MATRIX not defined!
+    #endif
+    #ifndef INFILE_WEIGHT_MATRIX
+      #error INFILE_WEIGHT_MATRIX not defined!
+    #endif
+  #else
+    #ifndef OUTFILE_CONNECTIONS
+      #error OUTFILE_CONNECTIONS not defined!
+    #endif
+    #ifndef OUTFILE_DELAYS
+      #error OUTFILE_DELAYS not defined!
+    #endif
+    #ifndef OUTFILE_STIMULUS
+      #error OUTFILE_STIMULUS not defined!
+    #endif
+  #endif
+#else
+  #error GENERATE_NETWORK_FROM_EXTERNAL_DATA not defined!
+#endif
+
+#ifdef USE_EXTERNAL_STIMULUS
+  #if( USE_EXTERNAL_STIMULUS )
+    #ifndef INFILE_STIMULUS
+      #error INFILE_STIMULUS not defined!
+    #endif
+  #endif
+#else
+  #error USE_EXTERNAL_STIMULUS not defined!
+#endif
+
+#ifndef USE_STDP
+  #error USE_STDP not defined!
+#endif
+
+#ifndef OUTFILE_FIRINGS
+  #error OUTFILE_FIRINGS not defined!
+#endif
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   F O R W A R D   D E C L A R A T I O N S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void InitializeNetwork();
+void InitializeSimulation();
+void FinalizeSimulation();
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   M A I N   E N T R Y
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+int main() {
+
+  printf( "\n\n%s\n\n", PARAM_INFO_STRING );
+
+  printf( "[INFO]  Running simulation for %d seconds\n", SIM_TIME );
+
+#if( ODE_SOLVER_REFINEMENT )
+  for( int i = 0; i < NUM_TOTAL_NEURONS; ++i ) {
+    spike[i] = 0;
+  }
+#endif
+
+  InitializeSimulation();
+  InitializeNetwork();
+
+  for( int simTimeSecond = 0; simTimeSecond < SIM_TIME; ++simTimeSecond ) {          // simulation loop [seconds]
+    for( int t = 0; t < 1000; ++t ) {                                                // simulation loop [milliseconds]
+
+      // clear I_ext[] and select a random neuron for input
+      for( int i = 0; i < NUM_TOTAL_NEURONS; ++i ) {
+        I_ext[i] = 0.0;
+      }
+
+      //  REFACTOR COMMENT: the original implementation does not allow networks sizes < 1000 neurons
+      //                    for( int idx = 0; idx < NUM_TOTAL_NEURONS / 1000; ++idx ) { ... }
+#if( USE_EXTERNAL_STIMULUS )
+      int inputNeuron = GetNextExternalStimulusFromFile(INFILE_STIMULUS, simTimeSecond, t);
+#else
+      int inputNeuron = GET_RANDOM_INT( NUM_TOTAL_NEURONS );
+#endif
+      if( inputNeuron > 0 && inputNeuron != RC_NOID) {
+        I_ext[inputNeuron] = I_EXT;
+      }
+
+#ifdef OUTFILE_STIMULUS
+      RecordRandomStimulusToFile( OUTFILE_STIMULUS, simTimeSecond, t, inputNeuron );
+#endif
+
+#if( ODE_SOLVER_REFINEMENT )
+      for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+        if( spike[n] > 0 ) {                                                         // Does neuron n has fired?
+
+          if( spike[n] > 1 ) printf( "[INFO]  More than 1 spike in simulation interval.\n" );
+          spike[n] = 0;
+
+          LTP[n][t + MAX_SYNAPSE_DELAY] = 0.1;
+          LTD[n] = 0.12;
+
+          for( int idx = 0; idx < numPreSynapticNeuronsOfTarget[n]; ++idx ) {        // LTP: pre-synaptic spike precedes post-synaptic spike
+            int preSynNeuron = listOfPresynapticNeurons[n][idx];
+            int preSynNeuron_correspondingDelay = listOfPresynapticDelays[n][idx];
+
+            *listOfPointersToSynapticWeights_derivatives[n][idx] += LTP[preSynNeuron][t + MAX_SYNAPSE_DELAY - preSynNeuron_correspondingDelay - 1];
+          }
+
+          firings[numFirings][TIME] = t;
+          firings[numFirings][NEURON] = n;
+          numFirings++;
+          if( numFirings == MAX_NUM_FIRINGS) {
+            printf( "[WARNING]  Two many spikes at t = %d (ignoring all)\n", t );
+            numFirings = 1;
+          }
+        }
+      }
+#else
+      for( int n= 0; n < NUM_TOTAL_NEURONS; ++n ) {
+        if( v[n] >= RS_FS_THR ) {                                                    // Does neuron n has fired?
+
+          v[n] = RS_FS_C;                                                            // threshold dynamics
+          u[n] += d[n];
+
+          LTP[n][t + MAX_SYNAPSE_DELAY] = 0.1;
+          LTD[n] = 0.12;
+
+          for ( int idx = 0; idx < numPreSynapticNeuronsOfTarget[n]; ++idx ) {       // LTP: pre-synaptic spike precedes post-synaptic spike
+            int preSynNeuron = listOfPresynapticNeurons[n][idx];
+            int preSynNeuron_correspondingDelay = listOfPresynapticDelays[n][idx];
+
+            *listOfPointersToSynapticWeights_derivatives[n][idx] += LTP[preSynNeuron][t + MAX_SYNAPSE_DELAY - preSynNeuron_correspondingDelay - 1];
+          }
+
+          firings[numFirings][TIME]   = t;
+          firings[numFirings][NEURON] = n;
+          numFirings++;
+          if( numFirings == MAX_NUM_FIRINGS ) {
+            printf( "[INFO]  Two many spikes at t = %d (ignoring all).\n", t );
+            numFirings = 1;
+          }
+        }
+      }
+#endif
+
+      //  Go back through firings as far as MAX_SYNAPSE_DELAY lasts. Calculate I_ext for the next time step.
+      //  Take the delays of the synapses into account. I_ext applies in the next simulation step.
+      //
+      //                          |        |                    |
+      //     |                          |     |  |  |     |  |  |  <-- spikes
+      //  ---+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--> simulation time
+      //    -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9 10 11 12
+      //                                                        t
+      //
+      //                                         |              |
+      //                                         +--------------+  MAX_SYNAPSE_DELAY (e.g., 5)
+
+      int idx = numFirings - 1;                                                      // array index starts with 0
+      while( t - firings[idx][TIME] < MAX_SYNAPSE_DELAY ) {
+
+        int neuronFired = firings[idx][NEURON];
+        int timeSinceNeuronFired = t - firings[idx][TIME];
+        int numEntriesOfDelay = numEntriesPerDelay[neuronFired][timeSinceNeuronFired];
+
+        for( int i = 0; i < numEntriesOfDelay; ++i ) {
+          int s = listOfSynsByNeuronAndDelay[neuronFired][timeSinceNeuronFired][i];
+          int postSynNeuron = matrixOfPostSynapticNeurons[neuronFired][s];
+
+          I_ext[postSynNeuron] += matrixOfSynapticWeights[neuronFired][s];
+
+          if( neuronFired < NUM_EXCITATORY_NEURONS) {                                  // LTD: this spike is before postsynaptic spikes
+            matrixOfSynapticWeights_derivatives[neuronFired][s] -= LTD[postSynNeuron];
+          }
+        }
+        idx--;
+      }
+
+      // = = = = = = = = = = = = = = = = = = = = = =
+      // = UPDATE NETWORK STATE
+      // = = = = = = = = = = = = = = = = = = = = = =
+      for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+
+#if( ODE_SOLVER_REFINEMENT )
+        for( int intCycles = 0; intCycles < ODE_SOLVER_STEPS; ++intCycles ) {
+          v[n] += (1.0 / ODE_SOLVER_STEPS) * ((0.04 * v[n] + 5) * v[n] + 140 - u[n] + I_ext[n]);
+          u[n] += (1.0 / ODE_SOLVER_STEPS) * (a[n] * (RS_FS_B * v[n] - u[n]));
+
+  #if( LOG_MIN_MAX_V_U )
+          // log min and max values of v(t) and u(t)
+          if( v[n] > v_max ) v_max = v[n];
+          if( v[n] < v_min ) v_min = v[n];
+          if( u[n] > u_max ) u_max = u[n];
+          if( u[n] < u_min ) u_min = u[n];
+  #endif
+
+          // threshold detection for exact integration
+          if( v[n] >= 30.0 ) {
+            v[n] = RS_FS_C;
+            u[n] += d[n];
+            spike[n]++;                                                              // remember the spike event, which is aligned to next grid point
+          }
+        }
+#else     // original implementation
+        v[n] += 0.5 * ((0.04 * v[n] + 5) * v[n] + 140 - u[n] + I_ext[n]);            // for numerical stability
+        v[n] += 0.5 * ((0.04 * v[n] + 5) * v[n] + 140 - u[n] + I_ext[n]);            // time step is 0.5 ms
+        u[n] += a[n] * (RS_FS_B * v[n] - u[n]);
+
+        #if( LOG_MIN_MAX_V_U )                                                       // log min and max values of v(t) and u(t)
+          if(v[n] > v_max)  v_max = v[n];
+          if(v[n] < v_min)  v_min = v[n];
+          if(u[n] > u_max)  u_max = u[n];
+          if(u[n] < u_min)  u_min = u[n];
+          #endif
+#endif
+
+        LTP[n][t + MAX_SYNAPSE_DELAY + 1] = 0.95 * LTP[n][t + MAX_SYNAPSE_DELAY];
+        LTD[n] *= 0.95;
+      }
+    }                                                                                // end of simulation loop [milliseconds]
+
+    // every minute: report on the average firing rate of the excitatory and inhibitory population
+    if( simTimeSecond > 0 && (simTimeSecond + 1) % 60 == 0 ) {
+      int spikesTotalExcitatory = 0;
+      int spikesTotalInhibitory = 0;
+      for( int idx = 0; idx < numFirings; ++idx ) {
+        if( firings[idx][NEURON] < NUM_EXCITATORY_NEURONS) {
+          spikesTotalExcitatory++;
+        }
+        else {
+          spikesTotalInhibitory++;
+        }
+      }
+      printf( "[INFO]  Simulation time: %d seconds\n", simTimeSecond + 1 );
+      printf( "[INFO]  ... Firing rate (excitatory) = %f\n",
+              (double) spikesTotalExcitatory / (double) NUM_EXCITATORY_NEURONS);
+      printf( "[INFO]  ... Firing rate (inhibitory) = %f\n",
+              (double) spikesTotalInhibitory / (double) NUM_INHIBITORY_NEURONS);
+
+#if(LOG_MIN_MAX_V_U)
+      printf( "[INFO]  ... v_max = %f\n", v_max );
+      printf( "[INFO]  ... v_min = %f\n", v_min );
+      printf( "[INFO]  ... u_max = %f\n", u_max );
+      printf( "[INFO]  ... u_min = %f\n", u_min );
+#endif
+    }
+
+    // = = = = = = = = = = = = = = = = = = = = = =
+    // = SAVE FIVE SELECTED NETWORK STATES
+    // = = = = = = = = = = = = = = = = = = = = = =
+#ifdef SELECTED_STATE_1_AFTER_N_SECONDS
+  #ifndef OUTFILE_STATE_1_WEIGHT_MATRIX
+    #error OUTFILE_STATE_1_WEIGHT_MATRIX not defined!
+  #endif
+    if( simTimeSecond == SELECTED_STATE_1_AFTER_N_SECONDS - 1 ) {
+      ExportWeightMatrixToFile( OUTFILE_STATE_1_WEIGHT_MATRIX );
+    }
+#endif
+
+#ifdef SELECTED_STATE_2_AFTER_N_SECONDS
+  #ifndef OUTFILE_STATE_2_WEIGHT_MATRIX
+    #error OUTFILE_STATE_2_WEIGHT_MATRIX not defined!
+  #endif
+    if( simTimeSecond == SELECTED_STATE_2_AFTER_N_SECONDS - 1 ) {
+      ExportWeightMatrixToFile( OUTFILE_STATE_2_WEIGHT_MATRIX );
+    }
+#endif
+
+#ifdef SELECTED_STATE_3_AFTER_N_SECONDS
+  #ifndef OUTFILE_STATE_3_WEIGHT_MATRIX
+    #error OUTFILE_STATE_3_WEIGHT_MATRIX not defined!
+  #endif
+    if( simTimeSecond == SELECTED_STATE_3_AFTER_N_SECONDS - 1 ) {
+      ExportWeightMatrixToFile( OUTFILE_STATE_3_WEIGHT_MATRIX );
+    }
+#endif
+
+#ifdef SELECTED_STATE_4_AFTER_N_SECONDS
+  #ifndef OUTFILE_STATE_4_WEIGHT_MATRIX
+    #error OUTFILE_STATE_4_WEIGHT_MATRIX not defined!
+  #endif
+    if( simTimeSecond == SELECTED_STATE_4_AFTER_N_SECONDS - 1 ) {
+      ExportWeightMatrixToFile( OUTFILE_STATE_4_WEIGHT_MATRIX );
+    }
+#endif
+
+#ifdef SELECTED_STATE_5_AFTER_N_SECONDS
+  #ifndef OUTFILE_STATE_5_WEIGHT_MATRIX
+    #error OUTFILE_STATE_5_WEIGHT_MATRIX not defined!
+  #endif
+    if( simTimeSecond == SELECTED_STATE_5_AFTER_N_SECONDS - 1 ) {
+      ExportWeightMatrixToFile( OUTFILE_STATE_5_WEIGHT_MATRIX );
+    }
+#endif
+
+    RecordNetworkActivityToFile( OUTFILE_FIRINGS, simTimeSecond, numFirings );
+
+    // = = = = = = = = = = = = = = = = = = = = = =
+    // = PREPARE NEXT SIMULATION TIME STEP
+    // = = = = = = = = = = = = = = = = = = = = = =
+    for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+      for( int i = 0; i < MAX_SYNAPSE_DELAY + 1; ++i ) {
+        LTP[n][i] = LTP[n][1000 + i];
+      }
+    }
+
+    int k = numFirings - 1;
+    while( 1000 - firings[k][TIME] < MAX_SYNAPSE_DELAY) {
+      k--;
+    }
+    for( int i = 1; i < numFirings - k; ++i ) {
+      firings[i][TIME] = firings[k + i][TIME] - 1000;
+      firings[i][NEURON] = firings[k + i][NEURON];
+    }
+    numFirings = numFirings - k;
+
+#if( USE_STDP )
+    // = = = = = = = = = = = = = = = = = = = = = =
+    // = UPDATE SYNAPTIC WEIGHTS
+    // = = = = = = = = = = = = = = = = = = = = = =
+    // only excitatory connections are modified
+    for( int n = 0; n < NUM_EXCITATORY_NEURONS; ++n ) {
+      for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+
+        // REFACTOR COMMENT:
+        // The following two lines have a different order in spnet.cpp.
+        matrixOfSynapticWeights_derivatives[n][s] *= 0.9;
+        matrixOfSynapticWeights[n][s] += 0.01 + matrixOfSynapticWeights_derivatives[n][s];
+
+        if( matrixOfSynapticWeights[n][s] > MAX_SYNAPTIC_STRENGTH) {
+          matrixOfSynapticWeights[n][s] = MAX_SYNAPTIC_STRENGTH;
+        }
+
+        if( matrixOfSynapticWeights[n][s] < 0 ) {
+          matrixOfSynapticWeights[n][s] = 0.0;
+        }
+      }
+    }
+#endif
+  }                                                                                  // end of simulation loop [seconds]
+
+  FinalizeSimulation();
+
+  printf( "\n\nSimulation terminated normally! \n\n" );
+
+  return (0);
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Initialize the polychronization network
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void InitializeNetwork() {
+
+  // initialize izhikevich neuron parameters
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    if( n < NUM_EXCITATORY_NEURONS ) {                                               // excitatory neurons, RS type
+      a[n] = RS_A;
+      d[n] = RS_D;
+      v[n] = RS_V_INIT;
+      u[n] = RS_U_INIT;
+    }
+    else {                                                                           // inhibitory neurons, FS type
+      a[n] = FS_A;
+      d[n] = FS_D;
+      v[n] = FS_V_INIT;
+      u[n] = FS_U_INIT;
+    }
+  }
+
+
+  // = = = = = = = = = = = = = = = = = = = = = =
+  // = CONNECTION MATRIX
+  // = = = = = = = = = = = = = = = = = = = = = =
+  #if( GENERATE_NETWORK_FROM_EXTERNAL_DATA )
+  ImportConnectionMatrixFromFile( INFILE_CONNECTION_MATRIX );
+#else
+  // fill matrixOfPostSynapticNeurons[n][n] with random target neurons
+  // avoid self-assignments and multiple connections
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+
+      bool duplicateTarget = false;
+      bool selfAssigned = false;
+      int randomTargetNeuron = 0;
+
+      do {
+        duplicateTarget = false;
+        selfAssigned = false;
+
+        if( n < NUM_EXCITATORY_NEURONS ) {                                           // excitatory neurons can connect to all neurons
+          randomTargetNeuron = GET_RANDOM_INT( NUM_TOTAL_NEURONS );
+        }
+        else {                                                                       // inhibitory neurons can connect to excitatory neurons only
+          randomTargetNeuron = GET_RANDOM_INT( NUM_EXCITATORY_NEURONS );
+        }
+        if( randomTargetNeuron == n ) {
+          selfAssigned = true;
+        }
+        for( int i = 0; i < s; ++i ) {
+          if( matrixOfPostSynapticNeurons[n][i] == randomTargetNeuron ) {
+            duplicateTarget = true;
+          }
+        }
+      } while( duplicateTarget || selfAssigned );
+
+      matrixOfPostSynapticNeurons[n][s] = randomTargetNeuron;
+    }
+  }
+#endif
+
+#ifdef OUTFILE_CONNECTIONS
+  ExportConnectionMatrixToFile( OUTFILE_CONNECTIONS );
+#endif
+#if __DEBUG__
+  PrintMatrixOfPostSynapticNeurons();
+#endif
+
+
+  // = = = = = = = = = = = = = = = = = = = = = =
+  // = WEIGHT MATRIX
+  // = = = = = = = = = = = = = = = = = = = = = =
+#if(GENERATE_NETWORK_FROM_EXTERNAL_DATA)
+  ImportWeightMatrixFromFile(INFILE_WEIGHT_MATRIX);
+#else
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+
+      if( n < NUM_EXCITATORY_NEURONS) {
+        matrixOfSynapticWeights[n][s] = INIT_EXC_SYNAPTIC_WEIGHT;
+      }
+      else {
+        matrixOfSynapticWeights[n][s] = INIT_INH_SYNAPTIC_WEIGHT;
+      }
+      matrixOfSynapticWeights_derivatives[n][s] = 0.0;
+    }
+  }
+#endif
+
+#ifdef OUTFILE_WEIGHTS_INITIAL
+  ExportWeightMatrixToFile(OUTFILE_WEIGHTS_INITIAL);
+#endif
+#if __DEBUG__
+  PrintMatrixOfSynapticWeights();
+#endif
+
+
+  // = = = = = = = = = = = = = = = = = = = = = =
+  // = DELAY MATRIX
+  // = = = = = = = = = = = = = = = = = = = = = =
+#if(GENERATE_NETWORK_FROM_EXTERNAL_DATA)
+  ImportDelayMatrixFromFile(INFILE_DELAY_MATRIX);
+#else
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    short synapseCorrespondingToDelay = 0;
+    for( int delayIdx = 0; delayIdx < MAX_SYNAPSE_DELAY; ++delayIdx ) {
+      numEntriesPerDelay[n][delayIdx] = 0;
+    }
+
+    if( n < NUM_EXCITATORY_NEURONS ) {
+      // For the excitatory neurons the delays are drawn from a uniform integer distribution.
+      // e.g.:
+      // delayIdx:     0      1      2      3      4      5       6 ..........      each slot corresponds to a delay value, i.e., 0 -> 1 ms, 1 -> 2 ms  etc.
+      //               0, 1   2, 3   4, 5   6, 7   8, 9   9, 10   11, 12   ...      each slot has space for NUM_SYNAPSES_PER_NEURON (which is not needed) and
+      //                                                                            holds the synapse numbers corresponding to that delay
+      // REFACTOR COMMENT: In the Matlab implementation this is random.
+      //                   This algorithm implicitly creates as many entries in delays as NUM_SYNAPSES_PER_NEURON.
+      //                   This only works if NUM_SYNAPSES_PER_NEURON is a multiple of MAX_SYNAPSE_DELAY!
+      for( int delayIdx = 0; delayIdx < MAX_SYNAPSE_DELAY; ++delayIdx ) {
+        for( int synListIdx = 0; synListIdx < (NUM_SYNAPSES_PER_NEURON / MAX_SYNAPSE_DELAY); ++synListIdx ) {
+          listOfSynsByNeuronAndDelay[n][delayIdx][synListIdx] = synapseCorrespondingToDelay++;
+          numEntriesPerDelay[n][delayIdx]++;
+        }
+      }
+    }
+    else {
+      // the inhibitory synapse delay is always 1 ms, thus, all synapses are added to delayIdx 0.
+      for( int synListIdx = 0; synListIdx < NUM_SYNAPSES_PER_NEURON; ++synListIdx ) {
+        listOfSynsByNeuronAndDelay[n][0][synListIdx] = synapseCorrespondingToDelay++;
+        numEntriesPerDelay[n][0]++;
+      }
+    }
+  }
+#endif
+
+  // for all neurons: find the excitatory presynaptic neurons and delays and store them in lists
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    numPreSynapticNeuronsOfTarget[n] = 0;
+  }
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int excPreSynNeuron = 0; excPreSynNeuron < NUM_EXCITATORY_NEURONS; ++excPreSynNeuron ) {
+      for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+        int targetNeuron = matrixOfPostSynapticNeurons[excPreSynNeuron][s];
+
+        if( targetNeuron == n ) {
+          int idx = numPreSynapticNeuronsOfTarget[n];
+          listOfPresynapticNeurons[n][idx] = excPreSynNeuron;
+
+          // determine the delay of a connection and store it in a list
+          for( int delayIdx = 0; delayIdx < MAX_SYNAPSE_DELAY; ++delayIdx ) {
+            int excPreSynNeuron_numEntriesOfDelay = numEntriesPerDelay[excPreSynNeuron][delayIdx];
+
+            for( int i = 0; i < excPreSynNeuron_numEntriesOfDelay; ++i ) {
+              short synapse = listOfSynsByNeuronAndDelay[excPreSynNeuron][delayIdx][i];
+              int postSynNeuron = matrixOfPostSynapticNeurons[excPreSynNeuron][synapse];
+
+              if( postSynNeuron == n ) {
+                int idx = numPreSynapticNeuronsOfTarget[n];
+                listOfPresynapticDelays[n][idx] = delayIdx;
+                // REFACTOR COMMENT: delayIdx represents the delay value, where 0 represents 1ms, thus,
+                //                   one needs to distinguish between a 1ms delay and an empty entry.
+                //                   This is done by tracking the number of entries in delays_numEntriesPerDelay.
+              }
+            }
+          }
+
+          // maintain list of pointers into the weight matrices; used for weight update
+          {
+            int idx = numPreSynapticNeuronsOfTarget[n];
+            // REFACTOR COMMENT: dead code removed
+            // listOfPointersToSynapticWeights[n][idx] = &matrixOfSynapticWeights[excPreSynNeuron][s];
+            listOfPointersToSynapticWeights_derivatives[n][idx] = &matrixOfSynapticWeights_derivatives[excPreSynNeuron][s];
+            idx++;
+            numPreSynapticNeuronsOfTarget[n] = idx;
+          }
+        }
+      }
+    }
+  }
+
+#if __DEBUG__
+  PrintMatrixOfSynapticDelays();
+#endif
+
+#ifdef OUTFILE_DELAYS
+  ExportDelayMatrixToFile( OUTFILE_DELAYS );
+#endif
+
+  // REFACTOR COMMENT: the array is much larger and remains uninitialized
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int d = 0; d < 1 + MAX_SYNAPSE_DELAY; ++d ) {
+      LTP[n][d] = 0.0;
+    }
+  }
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    LTD[n] = 0.0;
+  }
+
+  // REFACTOR COMMENT: just for the algorithm; does not contribute to the network activity
+  numFirings = 1;                                                                    // dummy spike ...
+  firings[0][TIME] = -MAX_SYNAPSE_DELAY;                                             // ... at -MAX_SYNAPSE_DELAY for ...
+  firings[0][NEURON] = 0;                                                            // ... neuron n = 0
+
+#ifdef OUTFILE_CON_WEIGHT_DELAY_INITIAL
+  ExportConnectionMatrixWeightAndDelay(OUTFILE_CON_WEIGHT_DELAY_INITIAL);
+#endif
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Initialize simulation: set seed, delete previously created files, open files
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void InitializeSimulation() {
+  srand( 0 );                                                                        // set a seed (repeatability)
+
+#if( LOG_MIN_MAX_V_U )
+  v_max = 0.0;
+  v_min = 0.0;
+  u_max = 0.0;
+  u_min = 0.0;
+#endif
+
+#ifdef OUTFILE_STIMULUS
+  DeleteFile( OUTFILE_STIMULUS );
+#endif
+
+#ifdef OUTFILE_CONNECTIONS
+  DeleteFile( OUTFILE_CONNECTIONS );
+#endif
+
+#ifdef OUTFILE_DELAYS
+  DeleteFile( OUTFILE_DELAYS );
+#endif
+
+#ifdef OUTFILE_CON_WEIGHT_DELAY_INITIAL
+  DeleteFile(OUTFILE_CON_WEIGHT_DELAY_INITIAL);
+#endif
+
+#ifdef OUTFILE_WEIGHTS_INITIAL
+  DeleteFile(OUTFILE_WEIGHTS_INITIAL);
+#endif
+
+  DeleteFile( OUTFILE_FIRINGS );
+
+  // open files which are required throughout the entire simulation
+#ifdef OUTFILE_STIMULUS
+  pFileStimulusOutput = fopen( OUTFILE_STIMULUS, "w" );
+  if( pFileStimulusOutput == nullptr ) {
+    printf( "[ERROR] Failed to open stimulus output file. Filename: %s\n", OUTFILE_STIMULUS );
+    exit( RC_ERROR_EXIT );
+  }
+#endif
+
+#if(USE_EXTERNAL_STIMULUS)
+  pFileStimulusInput = fopen(INFILE_STIMULUS,"r" );
+  if(pFileStimulusInput == nullptr) {
+    printf( "[ERROR] Failed to open stimulus input file. Filename: %s\n", INFILE_STIMULUS );
+    exit( RC_ERROR_EXIT );
+  }
+#endif
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Finalize simulation: free resources
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void FinalizeSimulation() {
+  if( pFileStimulusInput != nullptr ) fclose( pFileStimulusInput );
+  if( pFileStimulusOutput != nullptr ) fclose( pFileStimulusOutput );
+}
+
+
+#else                                                                                // original version
+  // original version of the Izhikevich polychronization model available for download at
+  // https://www.izhikevich.org/publications/spnet.htm
+  #include "poly_spnet.cpp"
+#endif

+ 58 - 0
simulation_code/C_model/src/params.h

@@ -0,0 +1,58 @@
+#ifndef __PARAMS_H__
+#define __PARAMS_H__
+/*
+ *  params.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   S I M U L A T I O N   S E T U P
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#include "params_create_5_selected_states.h"
+// #include "params_simulate_selected_state_1.h"
+// #include "params_simulate_selected_state_2.h"
+// #include "params_simulate_selected_state_3.h"
+// #include "params_simulate_selected_state_4.h"
+// #include "params_simulate_selected_state_5.h"
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   S O U R C E   C O D E   V E R S I O N   A N D   D E B U G
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define __DEBUG__                  false        // print connection, weight and delay matrix
+#define __RUN_REFACTORED_VERSION__ true         // run the refactored version of poly_spnet.cpp
+                                                // available for download at
+                                                // https://www.izhikevich.org/publications/spnet.htm
+                                                // set to false to run the original version
+#define LOG_MIN_MAX_V_U            true
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   O D E   S O L V E R   P A R A M E T E R S
+// =
+// =   The parameters do also apply to the poly_spnet.cpp source.
+// =   It has been modified to be able to run with the refined ODE solver to
+// =   test for the number of polychronous groups.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define ODE_SOLVER_REFINEMENT      true
+#define ODE_SOLVER_STEPS           (int)(16)
+
+
+#endif   // __PARAMS_H__

+ 60 - 0
simulation_code/C_model/src/params_create_5_selected_states.h

@@ -0,0 +1,60 @@
+#ifndef __PARAMS_CREATE_5_SELECTED_STATES_H__
+#define __PARAMS_CREATE_5_SELECTED_STATES_H__
+/*
+ *  params_create_5_selected_states.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Select five network states, after 1, 2, 3, 4, and 5 hours.
+// = Run the simulation for 5+ hours with STDP on, and save the five network states, that is, the connection and
+// = delay matrix (same for all states), and the weight matrices: w(t1), ... w(t5).
+// = Additionally, the random input data is recorded as well as the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "CREATE FIVE NETWORK STATES WITH STDP"
+
+#define SIM_TIME                              (int)(60*60*5+120)      // seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   false
+#define USE_EXTERNAL_STIMULUS                 false
+#define USE_STDP                              true
+
+#define SELECTED_STATE_1_AFTER_N_SECONDS      (int)(60*60*1)          // seconds
+#define SELECTED_STATE_2_AFTER_N_SECONDS      (int)(60*60*2)
+#define SELECTED_STATE_3_AFTER_N_SECONDS      (int)(60*60*3)
+#define SELECTED_STATE_4_AFTER_N_SECONDS      (int)(60*60*4)
+#define SELECTED_STATE_5_AFTER_N_SECONDS      (int)(60*60*5)
+
+#define OUTFILE_CONNECTIONS                   "../data/conMatrix.dat"
+#define OUTFILE_DELAYS                        "../data/delayMatrix.dat"
+
+#define OUTFILE_STATE_1_WEIGHT_MATRIX         "../data/weightMatrix_after1h.dat"
+#define OUTFILE_STATE_2_WEIGHT_MATRIX         "../data/weightMatrix_after2h.dat"
+#define OUTFILE_STATE_3_WEIGHT_MATRIX         "../data/weightMatrix_after3h.dat"
+#define OUTFILE_STATE_4_WEIGHT_MATRIX         "../data/weightMatrix_after4h.dat"
+#define OUTFILE_STATE_5_WEIGHT_MATRIX         "../data/weightMatrix_after5h.dat"
+
+#define OUTFILE_STIMULUS                      "../data/randomNetworkInput.dat"
+#define OUTFILE_FIRINGS                       "../data/firings.dat"
+
+#endif   // __PARAMS_CREATE_5_SELECTED_STATES_H__

+ 54 - 0
simulation_code/C_model/src/params_simulate_selected_state_1.h

@@ -0,0 +1,54 @@
+#ifndef __PARAMS_SIMULATE_SELECTED_STATE_1_H__
+#define __PARAMS_SIMULATE_SELECTED_STATE_1_H__
+/*
+ *  params_simulate_selected_state_1.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Reload the 1st selected network state, run the simulation for 60 seconds, and record the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "RELOAD 1st NETWORK STATE AND SIMULATE W/O STDP"
+
+#define SIM_TIME                              (int)(60)      // in seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   true
+#define USE_EXTERNAL_STIMULUS                 true
+#define USE_STDP                              false
+
+#define INFILE_CONNECTION_MATRIX              "../data/conMatrix.dat"
+#define INFILE_DELAY_MATRIX                   "../data/delayMatrix.dat"
+#define INFILE_STIMULUS                       "../data/randomNetworkInput.dat"
+#define INFILE_WEIGHT_MATRIX                  "../data/weightMatrix_after1h.dat"
+
+#define OUTFILE_FIRINGS                       "../data/firings_after1h.dat"
+
+// can be specified for verification purposes
+// #define OUTFILE_CONNECTIONS                "../data/conMatrix_initial.dat"
+// #define OUTFILE_WEIGHTS_INITIAL            "../data/weightMatrix_initial.dat"
+// #define OUTFILE_DELAYS                     "../data/delayMatrix_initial.dat"
+
+// generate connection lists for import into PyNN for SpiNNaker
+#define OUTFILE_CON_WEIGHT_DELAY_INITIAL      "../data/pythonConWeightDelay_after1h.py"
+
+#endif   // __PARAMS_SIMULATE_SELECTED_STATE_1_H__

+ 54 - 0
simulation_code/C_model/src/params_simulate_selected_state_2.h

@@ -0,0 +1,54 @@
+#ifndef __PARAMS_SIMULATE_SELECTED_STATE_2_H__
+#define __PARAMS_SIMULATE_SELECTED_STATE_2_H__
+/*
+ *  params_simulate_selected_state_2.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Reload the 1nd selected network state, run the simulation for 60 seconds, and record the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "RELOAD 2nd NETWORK STATE AND SIMULATE W/O STDP"
+
+#define SIM_TIME                              (int)(60)      // in seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   true
+#define USE_EXTERNAL_STIMULUS                 true
+#define USE_STDP                              false
+
+#define INFILE_CONNECTION_MATRIX              "../data/conMatrix.dat"
+#define INFILE_DELAY_MATRIX                   "../data/delayMatrix.dat"
+#define INFILE_STIMULUS                       "../data/randomNetworkInput.dat"
+#define INFILE_WEIGHT_MATRIX                  "../data/weightMatrix_after2h.dat"
+
+#define OUTFILE_FIRINGS                       "../data/firings_after2h.dat"
+
+// can be specified for verification purposes
+// #define OUTFILE_CONNECTIONS                "../data/conMatrix_initial.dat"
+// #define OUTFILE_WEIGHTS_INITIAL            "../data/weightMatrix_initial.dat"
+// #define OUTFILE_DELAYS                     "../data/delayMatrix_initial.dat"
+
+// generate connection lists for import into PyNN for SpiNNaker
+#define OUTFILE_CON_WEIGHT_DELAY_INITIAL      "../data/pythonConWeightDelay_after2h.py"
+
+#endif   // __PARAMS_SIMULATE_SELECTED_STATE_2_H__

+ 54 - 0
simulation_code/C_model/src/params_simulate_selected_state_3.h

@@ -0,0 +1,54 @@
+#ifndef __PARAMS_SIMULATE_SELECTED_STATE_3_H__
+#define __PARAMS_SIMULATE_SELECTED_STATE_3_H__
+/*
+ *  params_simulate_selected_state_3.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Reload the 3rd selected network state, run the simulation for 60 seconds, and record the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "RELOAD 3rd NETWORK STATE AND SIMULATE W/O STDP"
+
+#define SIM_TIME                              (int)(60)      // in seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   true
+#define USE_EXTERNAL_STIMULUS                 true
+#define USE_STDP                              false
+
+#define INFILE_CONNECTION_MATRIX              "../data/conMatrix.dat"
+#define INFILE_DELAY_MATRIX                   "../data/delayMatrix.dat"
+#define INFILE_STIMULUS                       "../data/randomNetworkInput.dat"
+#define INFILE_WEIGHT_MATRIX                  "../data/weightMatrix_after3h.dat"
+
+#define OUTFILE_FIRINGS                       "../data/firings_after3h.dat"
+
+// can be specified for verification purposes
+// #define OUTFILE_CONNECTIONS                "../data/conMatrix_initial.dat"
+// #define OUTFILE_WEIGHTS_INITIAL            "../data/weightMatrix_initial.dat"
+// #define OUTFILE_DELAYS                     "../data/delayMatrix_initial.dat"
+
+// generate connection lists for import into PyNN for SpiNNaker
+#define OUTFILE_CON_WEIGHT_DELAY_INITIAL      "../data/pythonConWeightDelay_after3h.py"
+
+#endif   // __PARAMS_SIMULATE_SELECTED_STATE_3_H__

+ 54 - 0
simulation_code/C_model/src/params_simulate_selected_state_4.h

@@ -0,0 +1,54 @@
+#ifndef __PARAMS_SIMULATE_SELECTED_STATE_4_H__
+#define __PARAMS_SIMULATE_SELECTED_STATE_4_H__
+/*
+ *  params_simulate_selected_state_4.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Reload the 4th selected network state, run the simulation for 60 seconds, and record the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "RELOAD 4th NETWORK STATE AND SIMULATE W/O STDP"
+
+#define SIM_TIME                              (int)(60)      // in seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   true
+#define USE_EXTERNAL_STIMULUS                 true
+#define USE_STDP                              false
+
+#define INFILE_CONNECTION_MATRIX              "../data/conMatrix.dat"
+#define INFILE_DELAY_MATRIX                   "../data/delayMatrix.dat"
+#define INFILE_STIMULUS                       "../data/randomNetworkInput.dat"
+#define INFILE_WEIGHT_MATRIX                  "../data/weightMatrix_after4h.dat"
+
+#define OUTFILE_FIRINGS                       "../data/firings_after4h.dat"
+
+// can be specified for verification purposes
+// #define OUTFILE_CONNECTIONS                "../data/conMatrix_initial.dat"
+// #define OUTFILE_WEIGHTS_INITIAL            "../data/weightMatrix_initial.dat"
+// #define OUTFILE_DELAYS                     "../data/delayMatrix_initial.dat"
+
+// generate connection lists for import into PyNN for SpiNNaker
+#define OUTFILE_CON_WEIGHT_DELAY_INITIAL      "../data/pythonConWeightDelay_after4h.py"
+
+#endif   // __PARAMS_SIMULATE_SELECTED_STATE_4_H__

+ 54 - 0
simulation_code/C_model/src/params_simulate_selected_state_5.h

@@ -0,0 +1,54 @@
+#ifndef __PARAMS_SIMULATE_SELECTED_STATE_5_H__
+#define __PARAMS_SIMULATE_SELECTED_STATE_5_H__
+/*
+ *  params_simulate_selected_state_5.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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/>.
+ *
+ */
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = P R O G R A M   P A R A M E T E R S
+// =
+// = Reload the 5th selected network state, run the simulation for 60 seconds, and record the network activity data.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define PARAM_INFO_STRING                     "RELOAD 5th NETWORK STATE AND SIMULATE W/O STDP"
+
+#define SIM_TIME                              (int)(60)      // in seconds
+
+#define GENERATE_NETWORK_FROM_EXTERNAL_DATA   true
+#define USE_EXTERNAL_STIMULUS                 true
+#define USE_STDP                              false
+
+#define INFILE_CONNECTION_MATRIX              "../data/conMatrix.dat"
+#define INFILE_DELAY_MATRIX                   "../data/delayMatrix.dat"
+#define INFILE_STIMULUS                       "../data/randomNetworkInput.dat"
+#define INFILE_WEIGHT_MATRIX                  "../data/weightMatrix_after5h.dat"
+
+#define OUTFILE_FIRINGS                       "../data/firings_after5h.dat"
+
+// can be specified for verification purposes
+// #define OUTFILE_CONNECTIONS                "../data/conMatrix_initial.dat"
+// #define OUTFILE_WEIGHTS_INITIAL            "../data/weightMatrix_initial.dat"
+// #define OUTFILE_DELAYS                     "../data/delayMatrix_initial.dat"
+
+// generate connection lists for import into PyNN for SpiNNaker
+#define OUTFILE_CON_WEIGHT_DELAY_INITIAL      "../data/pythonConWeightDelay_after5h.py"
+
+#endif   // __PARAMS_SIMULATE_SELECTED_STATE_5_H__

File diff suppressed because it is too large
+ 1077 - 0
simulation_code/C_model/src/poly_spnet.cpp


+ 627 - 0
simulation_code/C_model/src/utils.cpp

@@ -0,0 +1,627 @@
+/*
+ *  utils.cpp
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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 "params.h"
+
+#if( __RUN_REFACTORED_VERSION__ )
+
+#define __IMPORT_GLOBAL_VARIABLES__ true
+
+#include <math.h>
+#include <stdio.h>
+#include <errno.h>
+#include <stdlib.h>
+#include <cstring>
+
+#include "globals.h"
+#include "utils.h"
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   E X P O R T   N E T W O R K   S T A T E S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Export connection matrix to file in printable ASCII format.
+// =
+// =                      001 002                     100         <- synapse
+// = Source neuron 0001   nnn nnn nnn nnn nnn nnn ... nnn
+// = Source neuron 0002   nnn nnn nnn nnn nnn nnn ... nnn
+// = ...
+// = Source neuron 1000   nnn nnn nnn nnn nnn nnn ... nnn
+// =
+// = nnn ... target neuron id
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ExportConnectionMatrixToFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "w" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno);
+    exit(RC_ERROR_EXIT);
+  }
+
+  printf( "[INFO]  Export connection matrix to file: %s\n", pFileName );
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      fprintf( pFile, "%03d ", matrixOfPostSynapticNeurons[n][s] );
+    }
+    fprintf( pFile, "\n" );
+  }
+
+  fclose( pFile );
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Export weight matrix to file in printable ASCII format.
+// =
+// =                      001     002                         100        <- synapse
+// = Source neuron 0001   ww.wwww ww.wwww ww.wwww ww.wwww ... ww.wwww
+// = Source neuron 0002   ww.wwww ww.wwww ww.wwww ww.wwww ... ww.wwww
+// = ...
+// = Source neuron 1000   ww.wwww ww.wwww ww.wwww ww.wwww ... ww.wwww
+// =
+// = ww.wwww ... synaptic strength
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ExportWeightMatrixToFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "w" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno);
+    exit(RC_ERROR_EXIT);
+  }
+
+  printf( "[INFO]  Export weight matrix to file: %s\n", pFileName );
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      fprintf( pFile, "%f ", matrixOfSynapticWeights[n][s] );
+    }
+    fprintf( pFile, "\n" );
+  }
+
+  fclose( pFile );
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Export delay matrix to file in printable ASCII format.
+// =
+// =                      001 002                     100   <- synapse
+// = Source neuron 0001   ddd ddd ddd ddd ddd ddd ... ddd
+// = Source neuron 0002   ddd ddd ddd ddd ddd ddd ... ddd
+// = ...
+// = Source neuron 1000   ddd ddd ddd ddd ddd ddd ... ddd
+// =
+// = ddd ... conduction delay
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ExportDelayMatrixToFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "w" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno);
+    exit(RC_ERROR_EXIT);
+  }
+
+  printf( "[INFO]  Export delay matrix to file: %s\n", pFileName );
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      int delayValue = GetDelayOfConnection( n, s );
+      fprintf( pFile, "%03d ", delayValue );
+    }
+    fprintf( pFile, "\n" );
+  }
+
+  fclose( pFile );
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Export connection, weight, and delay matrix to file for import into PyNN.
+// =
+// =
+// = exc_exc_connections = [
+// =   (   source neuron, target neuron, weight, delay ), ...,
+// =   ... , (   source neuron, target neuron, weight, delay )
+// =                      ]
+// =
+// = exc_inh_connections = [
+// =   (   source neuron, target neuron, weight, delay ), ...,
+// =   ... , (   source neuron, target neuron, weight, delay )
+// =                      ]
+// = inh_exc_connections = [
+// =   (   source neuron, target neuron, weight, delay ), ...,
+// =   ... , (   source neuron, target neuron, weight, delay )
+// =                      ]
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ExportConnectionMatrixWeightAndDelay( const char *pFileName ) {
+  const int CONST_ENTRIES_PER_LINE = 8;
+
+  FILE *pFile = fopen( pFileName, "w" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno);
+    exit(RC_ERROR_EXIT);
+  }
+
+  printf( "[INFO]  Export connection, weight and delay matrix for import into PyNN to file: %s\n", pFileName );
+
+  // excitatory to excitatory connections
+  fprintf( pFile, "exc_exc_connections = [\n" );
+  int entriesCount = 0;
+  for( int excPreSynNeuron = 0; excPreSynNeuron < NUM_EXCITATORY_NEURONS; ++excPreSynNeuron ) {
+    for( int synapse = 0; synapse < NUM_SYNAPSES_PER_NEURON; ++synapse ) {
+      if( matrixOfPostSynapticNeurons[excPreSynNeuron][synapse] < NUM_EXCITATORY_NEURONS ) {
+
+        int excPostSynNeuron = matrixOfPostSynapticNeurons[excPreSynNeuron][synapse];
+        double weight = fabs( matrixOfSynapticWeights[excPreSynNeuron][synapse] );   // PyNN expects positive weight values for inhibitory synapses
+        double delay = GetDelayOfConnection( excPreSynNeuron, synapse );
+
+        fprintf( pFile, "( %3d, %3d, %f, %f ), ", excPreSynNeuron, excPostSynNeuron, weight, delay );
+
+        entriesCount++;
+        if( entriesCount == CONST_ENTRIES_PER_LINE ) {
+          fprintf( pFile, "\n" );
+          entriesCount = 0;
+        }
+      }
+    }
+  }
+  fprintf( pFile, "]\n\n" );
+
+  // excitatory to inhibitory connections
+  fprintf( pFile, "exc_inh_connections = [\n" );
+  entriesCount = 0;
+  for( int excPreSynNeuron = 0; excPreSynNeuron < NUM_EXCITATORY_NEURONS; ++excPreSynNeuron ) {
+    for( int synapse = 0; synapse < NUM_SYNAPSES_PER_NEURON; ++synapse ) {
+      if( matrixOfPostSynapticNeurons[excPreSynNeuron][synapse] >= NUM_EXCITATORY_NEURONS) {
+
+        int inhPostSynNeuron = matrixOfPostSynapticNeurons[excPreSynNeuron][synapse];
+        double weight = fabs( matrixOfSynapticWeights[excPreSynNeuron][synapse] );   // PyNN expects positive weight values for inhibitory synapses
+        double delay = GetDelayOfConnection( excPreSynNeuron, synapse );
+
+        fprintf( pFile, "( %3d, %3d, %f, %f ), "
+               , excPreSynNeuron, inhPostSynNeuron - NUM_EXCITATORY_NEURONS, weight, delay );
+
+        entriesCount++;
+        if( entriesCount == CONST_ENTRIES_PER_LINE ) {
+          fprintf( pFile, "\n" );
+          entriesCount = 0;
+        }
+      }
+    }
+  }
+  fprintf( pFile, "]\n\n" );
+
+  // inhibitory to excitatory connections
+  entriesCount = 0;
+  fprintf( pFile, "inh_exc_connections = [\n" );
+  for( int inhPreSynNeuron = NUM_EXCITATORY_NEURONS; inhPreSynNeuron < NUM_TOTAL_NEURONS; ++inhPreSynNeuron ) {
+    for( int synapse = 0; synapse < NUM_SYNAPSES_PER_NEURON; ++synapse ) {
+      if( matrixOfPostSynapticNeurons[inhPreSynNeuron][synapse] < NUM_EXCITATORY_NEURONS) {
+
+        int excPostSynNeuron = matrixOfPostSynapticNeurons[inhPreSynNeuron][synapse];
+        double weight = fabs( matrixOfSynapticWeights[inhPreSynNeuron][synapse] );   // PyNN expects positive weight values for inhibitory synapses
+        double delay = GetDelayOfConnection( inhPreSynNeuron, synapse );
+
+        fprintf( pFile, "( %3d, %3d, %6.3f, %6.3f ), "
+               , inhPreSynNeuron - NUM_EXCITATORY_NEURONS, excPostSynNeuron, weight, delay );
+
+        entriesCount++;
+        if( entriesCount == CONST_ENTRIES_PER_LINE ) {
+          fprintf( pFile, "\n" );
+          entriesCount = 0;
+        }
+      }
+    }
+  }
+  fprintf( pFile, "]\n\n" );
+
+  fclose( pFile );
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   I M P O R T   N E T W O R K   S T A T E S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Import the connection matrix from file generated with
+// = ExportConnectionMatrixToFile().
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ImportConnectionMatrixFromFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "r" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno);
+    exit(RC_ERROR_EXIT);
+  }
+
+  printf( "[INFO]  Import connection matrix from file: %s\n", pFileName );
+
+  char line[NUM_SYNAPSES_PER_NEURON * SIZE_OF_ENTRY_NEURON] = {};
+  int targetNeuron = 0;
+  char *pToken = nullptr;
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    fgets( line, sizeof( line ), pFile );
+
+    // first token, i.e. synapse 0 target neuron
+    pToken = strtok( line, " " );
+    targetNeuron = atoi( pToken );
+    matrixOfPostSynapticNeurons[n][0] = targetNeuron;
+
+    // subsequent tokens, i.e. synapses 1 .. s
+    for( int s = 1; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      pToken = strtok( NULL, " " );
+      targetNeuron = atoi( pToken );
+      matrixOfPostSynapticNeurons[n][s] = targetNeuron;
+    }
+  }
+
+  fclose( pFile );
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Import the weight matrix from file generated with
+// = ExportWeightMatrixToFile().
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ImportWeightMatrixFromFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "r" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno );
+    exit( RC_ERROR_EXIT );
+  }
+
+  printf( "[INFO]  Import weight matrix from file: %s\n", pFileName );
+
+  char line[NUM_SYNAPSES_PER_NEURON * SIZE_OF_ENTRY_WEIGHT] = {};
+  double weight = 0.0;
+  char *pToken = nullptr;
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    fgets( line, sizeof( line ), pFile );
+
+    // first token, i.e. weight of connection synapse 0 of target neuron n
+    pToken = strtok( line, " " );
+    weight = atof( pToken );
+    matrixOfSynapticWeights[n][0] = weight;
+
+    // subsequent tokens, i.e. weights of connections synapses s of target neuron n
+    for( int s = 1; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      pToken = strtok( NULL, " " );
+      weight = atof( pToken );
+      matrixOfSynapticWeights[n][s] = weight;
+    }
+  }
+
+  fclose( pFile );
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Import the delay matrix from file generated with
+// = ExportDelayMatrixToFile().
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ImportDelayMatrixFromFile( const char *pFileName ) {
+  FILE *pFile = fopen( pFileName, "r" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno );
+    exit( RC_ERROR_EXIT );
+  }
+
+  printf( "[INFO]  Import delay matrix from file: %s\n", pFileName );
+
+  char line[NUM_SYNAPSES_PER_NEURON * SIZE_OF_ENTRY_DELAY] = {};
+  char *pToken = nullptr;
+
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    fgets( line, sizeof( line ), pFile );
+
+    // first token, i.e. delay of synapse 0
+    pToken = strtok( line, " " );
+    int delayIdx = atoi( pToken ) - 1;
+
+    int idx = numEntriesPerDelay[n][delayIdx];
+    listOfSynsByNeuronAndDelay[n][delayIdx][idx] = 0;                        // synapse 0
+    numEntriesPerDelay[n][delayIdx]++;
+
+    // subsequent tokens, i.e. delays of synapses s
+    for( int s = 1; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      pToken = strtok( NULL, " " );
+      delayIdx = atoi( pToken ) - 1;
+
+      idx = numEntriesPerDelay[n][delayIdx];
+      listOfSynsByNeuronAndDelay[n][delayIdx][idx] = s;
+      numEntriesPerDelay[n][delayIdx]++;
+    }
+  }
+
+  fclose( pFile );
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   R E C O R D   E X T E R N A L   S T I M U L U S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Record the random input data the network is stimulated with.
+// =
+// = second  millisecond  id of the neuron that receives the input
+// = ssssss  mmm          nnn
+// = ...
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void RecordRandomStimulusToFile( const char *pFileName, int simTimeSecond, int simTimeMillisecond, int inputNeuron ) {
+  if( pFileStimulusOutput == nullptr ) {
+    printf( "[ERROR] File not open: %s   (errno: %d)\n", pFileName, errno );
+    exit( RC_ERROR_EXIT );
+  }
+
+  if( inputNeuron > 0 ) {
+    fprintf( pFileStimulusOutput, "%06d %03d %03d\n", simTimeSecond, simTimeMillisecond, inputNeuron );
+  }
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   R E A D   E X T E R N A L   S T I M U L U S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Read next record from external stimulus file and return the neuron id.
+// =
+// = second  millisecond  id of the neuron that receives the input
+// = ssssss  mmm          nnn
+// = ...
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+int GetNextExternalStimulusFromFile( const char *pFileName, int simTimeSecond, int t ) {
+  char line[32] = {};
+  char *pToken = nullptr;
+
+  if( pFileStimulusInput == nullptr ) {
+    printf( "[ERROR] File not open: %s   (errno: %d)\n", pFileName, errno );
+    exit( RC_ERROR_EXIT );
+  }
+
+  if( !fgets( line, sizeof( line ), pFileStimulusInput )) {
+    printf( "[ERROR] Unexpected end of file while reading stimulus data. Simulation time: %ds %dms  Filename: %s\n"
+          , simTimeSecond, t, pFileName );
+    exit( RC_ERROR_EXIT );
+  }
+
+  pToken = strtok( line, " " );
+  if( pToken == nullptr ) {
+    return RC_NOID;
+  }
+
+  pToken = strtok( nullptr, " " );
+  if( pToken == nullptr ) {
+    return RC_NOID;
+  }
+
+  pToken = strtok( nullptr, " " );
+  if( pToken == nullptr ) {
+    return RC_NOID;
+  }
+  int neuronId = atoi( pToken );
+
+  return neuronId;
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   R E C O R D   N E T W O R K   A C T I V I T Y   D A T A
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Write the spike times to file in a printable ASCII format.
+// =
+// = second  millisecond  id of the neuron that has fired
+// = ssssss  mmm          nnn
+// = ...
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void RecordNetworkActivityToFile( const char *pFileName, int simulationSecond, int numFirings ) {
+  FILE *pFile = fopen( pFileName, "a+" );
+  if( pFile == nullptr ) {
+    printf( "[ERROR] Failed to open file: %s   (errno: %d)\n", pFileName, errno );
+    exit( RC_ERROR_EXIT );
+  }
+
+  // skip negative times
+  int idx = 0;
+  while( firings[idx][TIME] < 0 ) {
+    idx++;
+  }
+
+  for( ; idx < numFirings; ++idx ) {
+    fprintf( pFile, "%06d %03d %03d\n", simulationSecond, firings[idx][TIME], firings[idx][NEURON] );
+  }
+
+  fclose( pFile );
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   H E L P E R   F U N C T I O N S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Return the delay value of a connection.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+int GetDelayOfConnection( int preSynNeuron, int synapse ) {
+  int delay = 0;
+
+  for( int delayIdx = 0; delayIdx < MAX_SYNAPSE_DELAY; ++delayIdx ) {
+    int preSynNeuron_numEntriesOfDelay = numEntriesPerDelay[preSynNeuron][delayIdx];
+    for( int i = 0; i < preSynNeuron_numEntriesOfDelay; ++i ) {
+      if( synapse == listOfSynsByNeuronAndDelay[preSynNeuron][delayIdx][i] ) {
+        delay = delayIdx + 1;                                                  // index 0 corresponds to 1ms delay
+      }
+    }
+  }
+
+  return delay;
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Delete a file from disk.
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void DeleteFile( const char *fileName ) {
+  if( remove( fileName ) == 0 ) {
+    printf( "[INFO]  File deleted: %s\n", fileName );
+  }
+  else {
+    printf( "[INFO]  File could not be deleted: %s\n", fileName );
+  }
+}
+
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   D E B U G   F U N C T I O N S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Debugprint: connection matrix
+// =
+// = Matrix of Post Synaptic Neurons
+// =
+// =          0    1    2    3    4    5    6    7    8    9  ...  <- synapse
+// = -------------------------------------------------------- ...
+// =  0 :   840  394  783  798  911  197  335  768  277  553  ...
+// =  1 :   950  920  147  881  641  431  619  281  786  307  ...  <- target neuron id
+// =  2 :    76  649  248  629  229  700  316  328  231   74  ...
+// =  3 :   291  180  684  727  139  603  492  838  724  178  ...
+// =  4 :    98  923  169  481  225  826  290  357  878  344  ...
+// =  5 :   324  874  589  637  759  775  794  262  604  470  ...
+// =  6 :   887  933  173  447  487  795  639  965  155  292  ...
+// =  7 :   347  205  522  400  307  679  645  443  269  703  ...
+// =  8 :   452  160  308  433    5  649  126  461   84  780  ...
+// =  9 :   805  749  398  366  394  272  599   68  901  432  ...
+// = 10 :    20   53  897  899   39  419  183  219  778  622  ...
+// = ....
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void PrintMatrixOfPostSynapticNeurons() {
+  printf( "     Matrix of Post Synaptic Neurons\n\n" );
+  printf( "        " );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( " %3d ", s );
+  }
+  printf( "\n" );
+  printf( "-------" );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( "-----" );
+  }
+  printf( " \n" );
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    printf( "%3d :   ", n );
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      printf( " %3d ", matrixOfPostSynapticNeurons[n][s] );
+    }
+    printf( "\n" );
+  }
+  printf( "Press enter ...\n" );
+  getchar();
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Debugprint: weight matrix
+// =
+// = Matrix of Synaptic Weights
+// =
+// =                0         1         2         3         4  ...  <- synapse
+// = --------------------------------------------------------- ...
+// =  0 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  1 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...  <- connection strength
+// =  2 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  3 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  4 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  5 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  6 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  7 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  8 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// =  9 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// = 10 :    6.000000  6.000000  6.000000  6.000000  6.000000  ...
+// = ....
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void PrintMatrixOfSynapticWeights() {
+  printf( "     Matrix of Synaptic Weights\n\n" );
+  printf( "        " );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( " %8d ", s );
+  }
+  printf( "\n" );
+  printf( "-------" );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( "----------" );
+  }
+  printf( " \n" );
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    printf( "%3d :   ", n );
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      printf( " %f ", matrixOfSynapticWeights[n][s] );
+    }
+    printf( "\n" );
+  }
+  printf( "Press enter ...\n" );
+  getchar();
+}
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Debugprint: delay matrix
+// =
+// = Matrix of Synaptic Delays
+// =
+// =           0    1    2    3    4    5    6    7    8    9   10  ...  <- synapse
+// = -------------------------------------------------------------- ...
+// =  0 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  1 :    001  001  001  001  001  002  002  002  002  002  003  ...  <- conduction delay
+// =  2 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  3 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  4 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  5 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  6 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  7 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  8 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// =  9 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// = 10 :    001  001  001  001  001  002  002  002  002  002  003  ...
+// = ....
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void PrintMatrixOfSynapticDelays() {
+  printf( "     Matrix of Synaptic Delays\n\n" );
+  printf( "        " );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( " %8d ", s );
+  }
+  printf( "\n" );
+  printf( "-------" );
+  for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+    printf( "----------" );
+  }
+  printf( " \n" );
+  for( int n = 0; n < NUM_TOTAL_NEURONS; ++n ) {
+    printf( "%3d :   ", n );
+    for( int s = 0; s < NUM_SYNAPSES_PER_NEURON; ++s ) {
+      printf( " %03d ", GetDelayOfConnection( n, s ));
+    }
+    printf( "\n" );
+  }
+  printf( "Press enter ...\n" );
+  getchar();
+}
+
+#endif   // __RUN_REFACTORED_VERSION__

+ 61 - 0
simulation_code/C_model/src/utils.h

@@ -0,0 +1,61 @@
+#ifndef __UTILS_H__
+#define __UTILS_H__
+/*
+ *  utils.h
+ *
+ *  This file is part of the refactored Izhikevich polychronization model application.
+ *
+ *  Copyright (C) 2018, Author: G. Trensch
+ *
+ *  The refactored Izhikevich polychronization model 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 "params.h"
+
+#if( __RUN_REFACTORED_VERSION__ )
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   F O R W A R D   D E C L A R A T I O N S
+/// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+void ExportConnectionMatrixToFile( const char *pFileName );
+void ExportWeightMatrixToFile( const char *pFileName );
+void ExportDelayMatrixToFile( const char *pFileName );
+void ExportConnectionMatrixWeightAndDelay( const char *pFileName );
+
+void ImportConnectionMatrixFromFile( const char *pFileName );
+void ImportWeightMatrixFromFile( const char *pFileName );
+void ImportDelayMatrixFromFile( const char *pFileName );
+
+void RecordNetworkActivityToFile( const char *pFileName, int simulationSecond, int numFirings );
+void RecordRandomStimulusToFile( const char *pFileName, int simTimeSecond, int simTimeMillisecond, int inputNeuron );
+int  GetNextExternalStimulusFromFile( const char *pFileName, int simTimeSecond, int t );
+
+int  GetDelayOfConnection( int preSynNeuron, int synapse );
+void DeleteFile( const char *pFileName );
+
+void PrintMatrixOfPostSynapticNeurons();
+void PrintMatrixOfSynapticWeights();
+void PrintMatrixOfSynapticDelays();
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// =   M A C R O   D E F I N I T I O N S
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+#define SIZE_OF_ENTRY_NEURON       (int)(5)
+#define SIZE_OF_ENTRY_WEIGHT       (int)(11)
+#define SIZE_OF_ENTRY_DELAY        (int)(11)
+
+#endif   // __RUN_REFACTORED_VERSION__
+#endif   // __UTILS_H__

+ 251 - 0
simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.c

@@ -0,0 +1,251 @@
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Set up desired ODE solver implementation
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+#define IMPL01  false    // original SpiNNaker ESR implementation
+#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
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+
+
+#include "neuron_model_izh_impl.h"
+#include <debug.h>
+
+static global_neuron_params_pointer_t global_params;
+
+/*! \brief For linear membrane voltages, 1.5 is the correct value. However
+ * with actual membrane voltage behaviour and tested over an wide range of
+ * use cases 1.85 gives slightly better spike timings.
+ */
+static const REAL SIMPLE_TQ_OFFSET = REAL_CONST(1.85);
+
+#if( IMPL01 )
+static inline void _rk2_kernel_midpoint(REAL h, neuron_pointer_t neuron,
+                                        REAL input_this_timestep) {
+
+    // to match Mathematica names
+    REAL lastV1 = neuron->V;
+    REAL lastU1 = neuron->U;
+    REAL a = neuron->A;
+    REAL b = neuron->B;
+
+    REAL pre_alph = REAL_CONST(140.0) + input_this_timestep - lastU1;
+    REAL alpha = pre_alph
+                 + ( REAL_CONST(5.0) + REAL_CONST(0.0400) * lastV1) * lastV1;
+    REAL eta = lastV1 + REAL_HALF(h * alpha);
+
+    // could be represented as a long fract?
+    REAL beta = REAL_HALF(h * (b * lastV1 - lastU1) * a);
+
+    neuron->V += h * (pre_alph - beta
+                      + ( REAL_CONST(5.0) + REAL_CONST(0.0400) * eta) * eta);
+
+    neuron->U += a * h * (-lastU1 - beta + b * eta);
+}
+#endif
+
+#if( IMPL02 )
+// h = 1.0, i.e., simulation timestep must be set to 1.0 !
+static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron,
+                                        REAL input_this_timestep) {
+   REAL v = neuron->V;
+   REAL u = neuron->U;
+   REAL param_a = neuron->A;
+   REAL param_b = neuron->B;
+
+   v += REAL_HALF(( 0.04k * v + 5.0k) * v + 140.0k - u + input_this_timestep );   // for numerical stability
+   v += REAL_HALF(( 0.04k * v + 5.0k) * v + 140.0k - u + input_this_timestep );   // time step is 0.5 ms
+   u += param_a * ( param_b * v - u );
+
+   neuron->V = v;
+   neuron->U = u;
+}
+#endif
+
+#if( IMPL03 )
+// h = 1.0, i.e., simulation timestep 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
+                                            // Hiding and encoding information in the data to control
+                                            // the program flow is not good practice! 
+   REAL v = neuron->V;
+   REAL u = neuron->U;
+   REAL param_a = neuron->A;
+   REAL param_b = neuron->B;
+   REAL a;
+   REAL b;
+   bool spikeEvent = false;
+   int steps = 16;
+ 
+   for( int cycles = 0; cycles; ++cycles ) {
+     a = 0.04k * v;
+     a *= v;
+
+     b = 5.0k * v + 140.0k - u + input_this_timestep;
+     v += 0.0625k * (a + b);                           // 0.0625 = 1/16
+  
+     u += 0.0625k * (param_a * ( param_b * v - u ));   // 0.0625 = 1/16
+ 
+     if( v >= 30.0k ) {
+       spikeEvent = true;
+       v = neuron->C;
+       u += neuron->D;
+     }
+   }
+ 
+   if( spikeEvent ) {
+     neuron->V = v + THRESHOLD_DETECTED;
+   }
+   else { 
+     neuron->V = v;
+   }
+   neuron->U = u;    
+}
+#endif
+
+#if( IMPL04 )
+// h = 1.0, i.e., simulation timestep 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
+
+   REAL v = neuron->V;
+   REAL u = neuron->U;
+   REAL param_a = neuron->A;
+   REAL param_b = neuron->B;
+
+   REAL param_a_RS = 5.12k;                // 0.02 RS param_a s8.23
+   REAL param_a_FS = 25.6k;                // 0.1  FS param_a s8.23
+   REAL param_b_FSRS = 51.2k;              // 0.2  FS param_a s8.23
+   REAL param_a_s8_23;
+   REAL a;
+   REAL b;
+   REAL c;
+   REAL d;
+   bool spikeEvent = false;
+
+   int steps = 16;
+ 
+   for( int cycles = 0; cycles < (steps + 2); ++cycles ) {  // (steps + 2) to supress unrolling the loop ..
+     if(cycles < steps ) {                                  // .. to prevent ITCM overflow (2 empty cyles)
+                                                            // (Corresponding compiler option did not work.)
+
+       a = 10.24k * v;                                      // byte left-shift
+       a *= 0.00390625k;                                    // byte right-shift
+       a *= v;
+
+       b = 5.0k * v + 140.0k - u + input_this_timestep;
+       v += 0.0625k * (a + b);                              // 0.0625 = 1/16
+  
+       if( param_a > 0.05k ) {                              // identify neuron type
+         param_a_s8_23 = param_a_FS;
+       }
+       else {
+         param_a_s8_23 = param_a_RS;
+       }
+
+       d = param_b_FSRS * v;
+       d *= 0.00390625k;                                   // byte right-shift 
+       d -= u;
+
+       c = 0.0625k * (param_a_s8_23 * ( d ));              // 0.0625 = 1/16
+ 
+       c *= 0.00390625k;                                   // byte right-shift 
+       u += c;
+
+       if( v >= 30.0k ) {
+         spikeEvent = true;
+         v = neuron->C;
+         u += neuron->D;
+       }
+     }
+   }
+ 
+   if( spikeEvent ) {
+     neuron->V = v + THRESHOLD_DETECTED;
+   }
+   else { 
+     neuron->V = v;
+   }
+   neuron->U = u;
+}
+#endif
+
+void neuron_model_set_global_neuron_params(
+        global_neuron_params_pointer_t params) {
+    global_params = params;
+}
+
+
+state_t neuron_model_state_update( input_t exc_input, input_t inh_input, input_t external_bias,
+                                   neuron_pointer_t neuron) {
+
+    input_t input_this_timestep = exc_input - inh_input + external_bias + neuron->I_offset;
+
+#if( IMPL01 )
+    // the best AR update so far
+    _rk2_kernel_midpoint(neuron->this_h, neuron, input_this_timestep);
+#endif
+
+#if( IMPL02 )
+    _originalIzhikevich(neuron->this_h, neuron, input_this_timestep);
+#endif
+
+#if( IMPL03 )
+    _fixedStepSizeEuler(neuron->this_h, neuron, input_this_timestep);
+#endif
+
+#if( IMPL04 )
+    _fixedStepSizeEulerWithFPConv(neuron->this_h, neuron, input_this_timestep);
+#endif
+
+    neuron->this_h = global_params->machine_timestep_ms;
+
+    return neuron->V;
+}
+
+
+void neuron_model_has_spiked(neuron_pointer_t neuron) {
+#if( IMPL01 )
+    // reset membrane voltage
+    neuron->V = neuron->C;
+
+    // offset 2nd state variable
+    neuron->U += neuron->D;
+
+    // simple threshold correction - next timestep (only) gets a bump
+    neuron->this_h = global_params->machine_timestep_ms * SIMPLE_TQ_OFFSET;
+#endif
+
+#if( IMPL02 )
+    neuron->V = neuron->C;
+    neuron->U += neuron->D;    
+#endif
+
+#if( IMPL03 || IMPL04 )
+    neuron->V -= THRESHOLD_DETECTED;
+#endif
+}
+
+state_t neuron_model_get_membrane_voltage(neuron_pointer_t neuron) {
+    return neuron->V;
+}
+
+void neuron_model_print_state_variables(restrict neuron_pointer_t neuron) {
+    log_debug("V = %11.4k ", neuron->V);
+    log_debug("U = %11.4k ", neuron->U);
+}
+
+void neuron_model_print_parameters(restrict neuron_pointer_t neuron) {
+    log_debug("A = %11.4k ", neuron->A);
+    log_debug("B = %11.4k ", neuron->B);
+    log_debug("C = %11.4k ", neuron->C);
+    log_debug("D = %11.4k ", neuron->D);
+
+    log_debug("I = %11.4k \n", neuron->I_offset);
+}
+

+ 30 - 0
simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.h

@@ -0,0 +1,30 @@
+#ifndef _NEURON_MODEL_IZH_CURR_IMPL_H_
+#define _NEURON_MODEL_IZH_CURR_IMPL_H_
+
+#include "neuron_model.h"
+
+typedef struct neuron_t {
+
+    // nominally 'fixed' parameters
+    REAL A;
+    REAL B;
+    REAL C;
+    REAL D;
+
+    // Variable-state parameters
+    REAL V;
+    REAL U;
+
+    // offset current [nA]
+    REAL I_offset;
+
+    // current timestep - simple correction for threshold
+    REAL this_h;
+
+} neuron_t;
+
+typedef struct global_neuron_params_t {
+    REAL machine_timestep_ms;
+} global_neuron_params_t;
+
+#endif   // _NEURON_MODEL_IZH_CURR_IMPL_H_

+ 299 - 0
simulation_code/SpiNNaker_model/PyNN/polychronizationNetwork.py

@@ -0,0 +1,299 @@
+#  polychronizationNetwork.py
+#
+#  This file is part of the SpiNNaker Izhikevich polychronization model implementation.
+#
+#  Copyright (C) 2018, Author: G. Trensch
+#
+#  The SpiNNaker Izhikevich polychronization model implementation 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/>.
+#
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   SpiNNaker Izhikevich polychronization network implementation w/o STDP.
+# =   Simulation of selected network states created with the C model implementation.
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+import spynnaker8 as sim
+from pyNN.utility.plotting import Figure, Panel
+import pylab
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+import seaborn
+
+
+SPIKEPLOT = True
+
+start_time = time.time()
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Import the network state that was generated with the C model
+# =   Parameterize according to the simulated network and state
+# =
+# =                            source  target  weight    delay
+# =   exc_exc_connections = [ (nnn,    nnn,    ww.wwww,  d ), () ... ]
+# =   exc_inh_connections = [ (), () ... ]
+# =   inh_exc_connections = [ (), () ... ]
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+EXTERNAL_CURRENT_FILE_NAME = 'randomNetworkInput.dat'         # random input
+
+
+from pythonConWeightDelay_after1h import *                    # 1st state: connections, weights, and delays
+FIRINGS_FILE_NAME  = 'firings_after1h.dat'                    #            network acivity data output file
+
+# from pythonConWeightDelay_after2h import *                  # 2nd state
+# FIRINGS_FILE_NAME  = 'firings_after2h.dat'                  #
+
+# from pythonConWeightDelay_after3h import *                  # 3rd state
+# FIRINGS_FILE_NAME  = 'firings_after3h.dat'                  #
+
+# from pythonConWeightDelay_after4h import *                  # 4th state
+# FIRINGS_FILE_NAME  = 'firings_after4h.dat'                  #
+
+# from pythonConWeightDelay_after5h import *                  # 5th state
+# FIRINGS_FILE_NAME  = 'firings_after5h.dat'                  #
+
+
+# population sizes have to correspond to the imported network state
+NUM_EXCITATORY_NEURONS  = 800
+NUM_INHIBITORY_NEURONS  = 200
+NUM_NEURONS             = NUM_EXCITATORY_NEURONS + NUM_INHIBITORY_NEURONS
+SYNAPSES_PER_NEURON     = 100
+CURRENT_INJECTION_VALUE = 20.0
+
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Simulation parameters
+# =   The simulation is carried out in 60 cycles with 1000 milliseconds simulation time each
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+SIM_TIME_TOTAL  = 1000 * 60                                  # [ms]
+SIM_TIME_SINGLE = 1000                                       # [ms]
+SIM_CYCLES      = SIM_TIME_TOTAL / SIM_TIME_SINGLE
+
+SIM_TIME_STEP = 1.0                                          # [ms] (1ms = real time)
+SIM_MIN_DELAY = 1.0
+SIM_MAX_DELAY = 144.0
+
+SIM_NEURONS_PER_CORE = 70
+
+sim.setup( timestep = SIM_TIME_STEP, min_delay= SIM_MIN_DELAY, max_delay = SIM_MAX_DELAY )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Neuron model and Izhikevich neuron parameters
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+NEURON_MODEL = sim.Izhikevich
+
+
+# Regular spiking type Izhikevich neuron
+NEURON_PARAMS_EXCITATORY = { 'a':          0.02
+                           , 'b':          0.2
+                           , 'c':        -65.0
+                           , 'd':          8.0
+                           , 'v_init':   -65.0
+                           , 'u_init':   -65.0 * 0.2         # according to the polychronization model
+                           , 'i_offset':   0.0
+                           }
+
+# Fast spiking type Izhikevich neuron
+NEURON_PARAMS_INHIBITORY = { 'a':          0.1
+                           , 'b':          0.2
+                           , 'c':        -65.0
+                           , 'd':          2.0
+                           , 'v_init':   -65.0
+                           , 'u_init':   -65.0 * 0.2        # according to the polychronization model
+                           , 'i_offset':   0.0
+                           }
+
+sim.set_number_of_neurons_per_core( NEURON_MODEL, SIM_NEURONS_PER_CORE )
+
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Populations: create the two populations of the polychronization  model
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+pop_exc = sim.Population( NUM_EXCITATORY_NEURONS, NEURON_MODEL(**NEURON_PARAMS_EXCITATORY) )
+pop_inh = sim.Population( NUM_INHIBITORY_NEURONS, NEURON_MODEL(**NEURON_PARAMS_INHIBITORY) )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Populations: create two spike source arrays for emulating external current and two empty lists,
+# =                filled at begin of every cycle with the spike sequence, that is, the random input data
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+spikes_inp_exc = list()
+for i in range( NUM_EXCITATORY_NEURONS ):
+    spikes_inp_exc.append( [] )
+
+pop_inp_exc = sim.Population( NUM_EXCITATORY_NEURONS, sim.SpikeSourceArray( spike_times = spikes_inp_exc ) )
+
+spikes_inp_inh = list()
+for i in range( NUM_INHIBITORY_NEURONS ):
+    spikes_inp_inh.append( [] )
+
+pop_inp_inh = sim.Population( NUM_INHIBITORY_NEURONS, sim.SpikeSourceArray( spike_times = spikes_inp_inh ) )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Projections: set up the poylchronization network
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+proj_exc_exc = sim.Projection( pop_exc, pop_exc
+                             , sim.FromListConnector( exc_exc_connections )
+                             , sim.StaticSynapse()
+                             , receptor_type = "excitatory" )
+
+proj_exc_inh = sim.Projection( pop_exc, pop_inh
+                             , sim.FromListConnector( exc_inh_connections )
+                             , sim.StaticSynapse()
+                             , receptor_type = "excitatory" )
+
+proj_inh_exc = sim.Projection( pop_inh, pop_exc
+                             , sim.FromListConnector( inh_exc_connections )
+                             , sim.StaticSynapse()
+                             , receptor_type = "inhibitory" )
+
+proj_inp_exc = sim.Projection( pop_inp_exc, pop_exc
+                             , sim.OneToOneConnector()
+                             , sim.StaticSynapse(weight = CURRENT_INJECTION_VALUE)
+                             , receptor_type = "excitatory" )
+
+proj_inp_inh = sim.Projection( pop_inp_inh, pop_inh
+                             , sim.OneToOneConnector()
+                             , sim.StaticSynapse(weight = CURRENT_INJECTION_VALUE)
+                             , receptor_type = "excitatory" )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Set up recording
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+pop_exc.record('spikes')
+pop_inh.record('spikes')
+pop_inp_exc.record('spikes')
+pop_inp_inh.record('spikes')
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Open files
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+inputCurrentFile = open( EXTERNAL_CURRENT_FILE_NAME, 'r' )
+
+if os.path.isfile( FIRINGS_FILE_NAME ):
+    os.remove( FIRINGS_FILE_NAME )
+
+outputSpikesFile = open( FIRINGS_FILE_NAME, 'a' )
+
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Set spike plot properties
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+if SPIKEPLOT:
+    seaborn.set()
+    seaborn.despine()
+    seaborn.set_style("white")
+    seaborn.axes_style("darkgrid")
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Run simulation
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+for cycle in range( SIM_CYCLES ):
+
+    simStartTime = SIM_TIME_SINGLE * cycle
+
+    spikes_inp_exc = list()
+    for i in range(NUM_EXCITATORY_NEURONS):
+        spikes_inp_exc.append([])
+
+    spikes_inp_inh = list()
+    for i in range(NUM_INHIBITORY_NEURONS):
+        spikes_inp_inh.append([])
+
+    # load random input date from file for this cycle
+    for tickCount in range( SIM_TIME_SINGLE ):
+        inputDataForThisTick = inputCurrentFile.readline().split()
+        inputNeuronForThisTick = int(inputDataForThisTick[2])
+
+        if( inputNeuronForThisTick < NUM_EXCITATORY_NEURONS ):
+            spikes_inp_exc[inputNeuronForThisTick].append( tickCount + simStartTime )
+
+        else:
+            inputNeuronForThisTick = inputNeuronForThisTick - NUM_EXCITATORY_NEURONS
+            spikes_inp_inh[inputNeuronForThisTick].append( tickCount + simStartTime )
+
+    # set spike source array with random input data fo this cycle
+    pop_inp_exc.set( spike_times = spikes_inp_exc )
+    pop_inp_inh.set( spike_times = spikes_inp_inh )
+
+
+    sim.run( SIM_TIME_SINGLE )
+
+    # retrieve input spikes from populations, i.e., from spike source arrays
+    spikesInpToExcPop = pop_inp_exc.get_data()
+    spikesInpToInhPop = pop_inp_inh.get_data()
+
+    # retrieve output spikes from populations and clear segment
+    spikesExcPop = pop_exc.get_data('spikes', clear = True)
+    spikesInhPop = pop_inh.get_data('spikes', clear = True)
+
+    simEndTime = sim.get_current_time()
+    print( simStartTime, simEndTime )
+
+    if SPIKEPLOT:
+        Figure( Panel(spikesInpToInhPop.segments[0].spiketrains
+              , xticks = True, yticks = True, markersize = 2.0, xlim = (simStartTime, simEndTime))
+
+              , Panel(spikesInpToExcPop.segments[0].spiketrains
+              , xticks = True, yticks = True, markersize = 2.0, xlim = (simStartTime, simEndTime))
+
+              , Panel(spikesInhPop.segments[0].spiketrains
+              , xticks = True, yticks = True, markersize = 2.0, xlim = (simStartTime, simEndTime))
+
+              , Panel(spikesExcPop.segments[0].spiketrains
+              , xticks = True, yticks = True, markersize = 2.0, xlim = (simStartTime, simEndTime), xlabel = "Time [ms]")
+
+              , title="Plochronization Network", annotations="Simulated with {}".format(sim.name())
+              )
+
+        plt.pause(0.1)
+
+    # write firings to file
+    firings = list()
+    for i in range( SIM_TIME_SINGLE ):
+        firings.append([])
+
+    for neuron in range( NUM_EXCITATORY_NEURONS ):
+        for spikeTimeIdx in range( len(spikesExcPop.segments[0].spiketrains[neuron]) ):
+            firingTime = int(spikesExcPop.segments[0].spiketrains[neuron][spikeTimeIdx]) - simStartTime
+            firings[firingTime].append( neuron )
+
+    for neuron in range( NUM_INHIBITORY_NEURONS ):
+        for spikeTimeIdx in range( len(spikesInhPop.segments[0].spiketrains[neuron]) ):
+            firingTime = int(spikesInhPop.segments[0].spiketrains[neuron][spikeTimeIdx]) - simStartTime
+            firings[firingTime].append( neuron + NUM_EXCITATORY_NEURONS )
+
+    for i in range( SIM_TIME_SINGLE ):
+        simCycle = str( '{:06d}'.format(cycle) )
+        tickInCycle = str('{:06d}'.format(i))
+
+        for n in firings[i]:
+            neuronId = str('{:03d}'.format(n))
+            line = ' '.join([simCycle, tickInCycle, neuronId, "\n"])
+            outputSpikesFile.writelines( line )
+
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   End simulation, close files, show spike plot
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+sim.end()
+
+inputCurrentFile.close()
+outputSpikesFile.close()
+
+print("Simulation time: %s seconds" % (time.time() - start_time))
+
+if SPIKEPLOT:
+    pylab.show()

+ 119 - 0
simulation_code/SpiNNaker_model/PyNN/singleNeuronDynamics.py

@@ -0,0 +1,119 @@
+#
+#  singleNeuronDynamics.py
+#
+#  This file is part of the SpiNNaker Izhikevich polychronization model implementation.
+#
+#  Copyright (C) 2018, Author: G. Trensch
+#
+#  The SpiNNaker Izhikevich polychronization model implementation 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/>.
+#
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Plot the individual neuron dynamics, i.e., v(t), for a regular spiking (RS type) and a fast spiking
+# =   (FS type) Izhikevich neuron, that is stimulated with an external current of 5.0pA.
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+
+import spynnaker8 as sim
+from pyNN.utility.plotting import Figure, Panel
+import pylab
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Simulation time and resolution
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+SIM_TIME = 1000    # [ms]
+sim.setup( timestep = 1.0 )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Neuron model and Izhikevich neuron parameters
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+NEURON_MODEL = sim.Izhikevich
+
+# Regular spiking type Izhikevich neuron
+NEURON_PARAMS_RS = { 'a':          0.02
+                   , 'b':          0.2
+                   , 'c':        -65.0
+                   , 'd':          8.0
+                   , 'v_init':   -75.0
+                   , 'u_init':     0.0
+                   , 'i_offset':   5.0
+                   }
+
+# Fast spiking type Izhikevich neuron
+NEURON_PARAMS_FS = { 'a':          0.1
+                   , 'b':          0.2
+                   , 'c':        -65.0
+                   , 'd':          2.0
+                   , 'v_init':   -75.0
+                   , 'u_init':     0.0
+                   , 'i_offset':   5.0
+                   }
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Func: record membrane voltages to file in a printable ASCII format
+# =
+# =   t [ms]   v
+# =   000     -75.0 mV
+# =   001     -77.9880371094 mV
+# =   002     -78.662109375 mV
+# =   003     -78.6662902832 mV
+# =   004     -78.4951477051 mV
+# =   ...
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+def writeVtoDisk( population, fileName ):
+    outFile = open(fileName, 'w')
+    x = population.get_data('v')
+    v = x.segments[0].filter(name='v')[0]
+
+    for i in range(len(v)):
+        time = str('{:03d}'.format(i))
+        voltage = str(v[i][0])
+        line = ' '.join([time, voltage, "\n"])
+        outFile.writelines(line)
+
+    outFile.close()
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Create an RS-type and FS-type Izhikevich neuron
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+neuron_RS = sim.Population( 1, NEURON_MODEL(**NEURON_PARAMS_RS) )
+neuron_FS = sim.Population( 1, NEURON_MODEL(**NEURON_PARAMS_FS) )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Set up recording of the membrane voltages for both neurons
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+neuron_RS.record('v')
+neuron_FS.record('v')
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Run simulation
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+sim.run( SIM_TIME )
+
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+# =   Retrieve and plot data, and write data to disk
+# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+v_neuron_RS = neuron_RS.get_data('v').segments[0].filter(name='v')[0]
+v_neuron_FS = neuron_FS.get_data('v').segments[0].filter(name='v')[0]
+
+writeVtoDisk( neuron_RS, "v(t)_RS_type.dat")
+writeVtoDisk( neuron_FS, "v(t)_FS_type.dat")
+
+Figure( Panel(v_neuron_RS, xticks = True, yticks = True, markersize = 2.0, xlabel = "regular spiking")
+      , Panel(v_neuron_FS, xticks = True, yticks = True, markersize = 2.0, xlabel = "fast spiking")
+      , annotations="Simulated with {}".format(sim.name())
+      )
+
+sim.end()
+pylab.show()