123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291 |
- // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
- // = 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"
- #include <debug.h>
- static global_neuron_params_pointer_t global_params;
- #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);
- 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 )
- // 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 ); // 2 time steps of 0.5 ms
- u += param_a * ( param_b * v - u );
- neuron->V = v;
- neuron->U = u;
- }
- #endif
- #if( IMPL03 )
- // 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;
- REAL param_b = neuron->B;
- REAL a;
- REAL b;
- bool spikeEvent = false;
- int steps = 16;
-
- 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, avoid costly division
-
- 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; // 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;
- }
- neuron->U = u;
- }
- #endif
- #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_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;
- 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; // 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;
- }
- 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 )
- _rk2_kernel_midpoint(neuron->this_h, neuron, input_this_timestep);
- #endif
- #if( IMPL02 )
- _originalIzhikevich(neuron->this_h, neuron, input_this_timestep);
- #endif
- #if( IMPL03 )
- _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 || IMPL02 || IMPL03 )
- neuron->V = neuron->C;
- 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( IMPL04 || IMPL05 )
- // remove metadata from the membran voltage after the spike processing
- // has been triggered in the SpiNNaker kernel
- 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);
- }
|