123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671 |
- /*
- * 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
|