updating simulations codes from github.com/gtrensch/

simulation_code/C_model/src/params.h

@@ -51,8 +51,9 @@
 // =   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)
+// #define ODE_SOLVER_REFINEMENT   false        // corresponds to Iteration I, Trensch et al., Gutzen et al.        
 
+#define ODE_SOLVER_REFINEMENT      true         // corresponds to Iteration II and III, Trensch et al., Gutzen et al.        
+#define ODE_SOLVER_STEPS           (int)(16)
 
 #endif   // __PARAMS_H__

simulation_code/README.md

@@ -1,3 +1,5 @@
+***Copied from [github.com/gtrensch/RigorousNeuralNetworkSimulations](https://github.com/gtrensch/RigorousNeuralNetworkSimulations)***
+
 These are the C/C++ source files and PyNN scripts developed as part of a study to reproduce selected network states of the Izhikevich polychronization model on the SpiNNaker neuromorphic system.
 
 **STATUS**: Several releases for this source code have been published on Zenodo: 

simulation_code/SpiNNaker_model/C/neuron_model_izh_impl.c

@@ -1,14 +1,21 @@
-// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
-// = Set up desired ODE solver implementation
-// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
-
-#define IMPL01  false    // original SpiNNaker ESR implementation (corresponds to Iteration I)
-#define IMPL02  false    // original Izhikevich implementation
-#define IMPL03  false    // fixed-step size (h=1/16) symplectic forward Euler, precise threshold detection
-#define IMPL04  true     // fixed-step size (h=1/16) symplectic forward Euler, precise threshold detection, FP conversion (corresponds to Iteration III)
-
-// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
-
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+// = Set up the desired ODE solver implementation coresponding to the
+// = iterations described in Trensch et al. and Gutzen et al.
+// =
+// = For all implementations the simulation time step must be set to 1.0, i.e,
+// = the value of SIM_TIME_STEP in polychronizationNetwork.py
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
+                                                                                    
+#define IMPL01  false    // SpiNNaker ESR implementation                             corresponds to Iteration I, (i)   Gutzen et al.
+#define IMPL02  false    // original Izhikevich implementation                                                   (ii)  Gutzen et al.
+#define IMPL03  false    // original Izhikevich implementation with FP conversion                                (iii) Gutzen et al.
+
+#define IMPL04  false    // fixed-step size (h=1/16) symplectic forward Euler and    corresponds to Iteration II plus III
+                         // precise threshold detection
+#define IMPL05  true     // fixed-step size (h=1/16) symplectic forward Euler, 
+                         // precise threshold detection and fixed-point conversion
+
+// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 
 
 #include "neuron_model_izh_impl.h"
@@ -16,15 +23,15 @@
 
 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.
- */
+#if( IMPL01 )  
+// Original SpiNNaker implementation
+// Solver: Explicit Solver Reduction (ESR) implementation
+
+// 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) {
+static inline void _rk2_kernel_midpoint(REAL h, neuron_pointer_t neuron, REAL input_this_timestep) {
 
     // to match Mathematica names
     REAL lastV1 = neuron->V;
@@ -33,31 +40,29 @@ static inline void _rk2_kernel_midpoint(REAL h, neuron_pointer_t neuron,
     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 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->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 time-step must be set to 1.0 !
-static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron,
-                                        REAL input_this_timestep) {
+#if( IMPL02 )  
+// Original Izhikevich implementation 
+// Solver: semi-implicit fixed-step size Forward Euler
+
+static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron, REAL input_this_timestep) {
+   // h is set to a fixed value: h = 1.0, i.e., simulation time-step must be set to 1.0!
    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
+   v += REAL_HALF(( 0.04k * v + 5.0k) * v + 140.0k - u + input_this_timestep );   // 2 time steps of 0.5 ms
    u += param_a * ( param_b * v - u );
 
    neuron->V = v;
@@ -66,12 +71,43 @@ static inline void _originalIzhikevich( REAL h, neuron_pointer_t neuron,
 #endif
 
 #if( IMPL03 )
-// h = 1.0, i.e., simulation time-step must be set to 1.0 !
-static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
-                                        REAL input_this_timestep) {
-   #define THRESHOLD_DETECTED  1024.0k      // some high value above threshold to trigger spike event
-                                            // Hiding and encoding information in the data to control
-                                            // the program flow is not good practice! 
+// Original Izhikevich implementation with fixed-point conversion of constant 0.04 (i.e., s16.15 -> s8.23)
+// Solver: semi-implicit fixed-step size Forward Euler
+
+static inline void _originalIzhikevichWithFPConv( 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;
+   REAL a;
+   REAL b;
+ 
+   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 += REAL_HALF(a + b);
+ 
+   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 += REAL_HALF(a + b);
+ 
+   u += param_a * ( param_b * v - u );
+ 
+   neuron->V = v;
+   neuron->U = u;
+}
+#endif
+
+#if( IMPL04 )
+// Refined ODE solver implementation including a precise threshold detection
+// Solver: Fixed-step size (1/16) semi-implicit symplectic Forward Euler
+
+static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron, REAL input_this_timestep) {
+   #define THRESHOLD_DETECTED  1024.0k      // Some high value above threshold, used to trigger a spike event
+                                            // in the SpiNNaker kernel.
    REAL v = neuron->V;
    REAL u = neuron->U;
    REAL param_a = neuron->A;
@@ -81,12 +117,12 @@ static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
    bool spikeEvent = false;
    int steps = 16;
  
-   for( int cycles = 0; cycles; ++cycles ) {
+   for( int cycles = 0; cycles < steps; ++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
+     v += 0.0625k * (a + b);                           // 0.0625 = 1/16, avoid costly division
   
      u += 0.0625k * (param_a * ( param_b * v - u ));   // 0.0625 = 1/16
  
@@ -98,7 +134,11 @@ static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
    }
  
    if( spikeEvent ) {
-     neuron->V = v + THRESHOLD_DETECTED;
+     neuron->V = v + THRESHOLD_DETECTED;    // Overlay the data with the spike event information.
+                                            // Note: Hiding and encoding information in the data to 
+                                            //       control the program flow is not good practice!
+                                            //       This has been done here for the sake of simplicity
+                                            //       and to avoid changes in the SpiNNaker kernel.
    }
    else { 
      neuron->V = v;
@@ -107,20 +147,20 @@ static inline void _fixedStepSizeEuler( REAL h, neuron_pointer_t neuron,
 }
 #endif
 
-#if( IMPL04 )
-// h = 1.0, i.e., simulation time-step must be set to 1.0 !
-static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuron,
-                                                  REAL input_this_timestep) {
-   #define THRESHOLD_DETECTED  1024.0k      // some high value above threshold to trigger spike event
+#if( IMPL05 )
+// Refined ODE solver implementation including a precise threshold detection and fixed-point conversion
+// Solver: Fixed-step size (1/16) semi-implicit symplectic Forward Euler
 
+static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuron, REAL input_this_timestep) {
+   #define THRESHOLD_DETECTED  1024.0k      // Some high value above threshold, used to trigger a spike event
+                                            // in the SpiNNaker kernel.
    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_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 RS param_b s8.23
    REAL param_a_s8_23;
    REAL a;
    REAL b;
@@ -149,12 +189,12 @@ static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuro
        }
 
        d = param_b_FSRS * v;
-       d *= 0.00390625k;                                   // byte right-shift 
+       d *= 0.00390625k;                                    // byte right-shift 
        d -= u;
 
-       c = 0.0625k * (param_a_s8_23 * ( d ));              // 0.0625 = 1/16
+       c = 0.0625k * (param_a_s8_23 * ( d ));               // 0.0625 = 1/16
  
-       c *= 0.00390625k;                                   // byte right-shift 
+       c *= 0.00390625k;                                    // byte right-shift 
        u += c;
 
        if( v >= 30.0k ) {
@@ -166,7 +206,11 @@ static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuro
    }
  
    if( spikeEvent ) {
-     neuron->V = v + THRESHOLD_DETECTED;
+     neuron->V = v + THRESHOLD_DETECTED;    // Overlay the data with the spike event information.
+                                            // Note: Hiding and encoding information in the data to 
+                                            //       control the program flow is not good practice!
+                                            //       This has been done here for the sake of simplicity
+                                            //       and to avoid changes in the SpiNNaker kernel.
    }
    else { 
      neuron->V = v;
@@ -175,8 +219,7 @@ static inline void _fixedStepSizeEulerWithFPConv( REAL h, neuron_pointer_t neuro
 }
 #endif
 
-void neuron_model_set_global_neuron_params(
-        global_neuron_params_pointer_t params) {
+void neuron_model_set_global_neuron_params( global_neuron_params_pointer_t params ) {
     global_params = params;
 }
 
@@ -187,7 +230,6 @@ state_t neuron_model_state_update( input_t exc_input, input_t inh_input, input_t
     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
 
@@ -196,37 +238,36 @@ state_t neuron_model_state_update( input_t exc_input, input_t inh_input, input_t
 #endif
 
 #if( IMPL03 )
-    _fixedStepSizeEuler(neuron->this_h, neuron, input_this_timestep);
+    _originalIzhikevichWithFPConv(neuron->this_h, neuron, input_this_timestep);
 #endif
 
 #if( IMPL04 )
+    _fixedStepSizeEuler(neuron->this_h, neuron, input_this_timestep);
+#endif
+
+#if( IMPL05 )
     _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
+#if( IMPL01 || IMPL02 || IMPL03 )
     neuron->V = neuron->C;
-
-    // offset 2nd state variable
     neuron->U += neuron->D;
+#endif 
 
+#if( IMPL01 )
     // 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 )
+#if( IMPL04 || IMPL05 )
+    // remove metadata from the membran voltage after the spike processing
+    // has been triggered in the SpiNNaker kernel
     neuron->V -= THRESHOLD_DETECTED;
 #endif
 }
@@ -248,4 +289,3 @@ void neuron_model_print_parameters(restrict neuron_pointer_t neuron) {
 
     log_debug("I = %11.4k \n", neuron->I_offset);
 }
-