123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251 |
- // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
- // = 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);
- }
|