SNN.c

#include <ncurses.h>
#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

/* --- Simulation parameters --- */
#define NUM_NEURONS         10      // number of neurons
#define NUM_SYNAPSES        20      // number of synaptic connections
#define SIM_TIME            1000    // total simulation time in ms
#define DT                  1       // simulation time step (ms)
#define VISUAL_DELAY        30      // delay (ms) between visual updates

/* --- Neuron parameters --- */
#define THRESHOLD           1.0     // firing threshold
#define RESET_POTENTIAL     0.0     // potential after spike
#define REFRACTORY_PERIOD   5       // refractory period (ms)
#define LEAK_FACTOR         0.95    // leak (decay) factor per time step

/* --- STDP parameters --- */
#define A_PLUS              0.01    // potentiation factor
#define A_MINUS             -0.012  // depression factor
#define TAU_PLUS            20.0    // time constant for potentiation (ms)
#define TAU_MINUS           20.0    // time constant for depression (ms)
#define STDP_WINDOW         50      // maximum spike-timing difference (ms)

/* --- Neurotransmitter dynamics --- */
#define NT_INITIAL          1.0     // full availability
#define NT_RECOVERY_RATE    0.01    // recovery rate per ms

/* --- Visualization parameters --- */
#define BAR_WIDTH           20      // width of the potential bar

/* --- Conduction delay & Noise parameters --- */
#define MAX_DELAY           5       // maximum synaptic conduction delay (ms)
#define NOISE_CHANCE        5       // 5% chance for noise pulse per neuron per ms
#define NOISE_AMPLITUDE     0.2     // amplitude of noise pulse

/* --- Structures --- */

/* Each neuron tracks its membrane potential, refractory state, and spike count. */
typedef struct {
    double potential;
    int refractory;   // remaining refractory time (ms)
    int spike_count;  // total spikes produced (for summary)
} Neuron;

/* Each synapse connects two neurons and tracks spike times for STDP, and now a delay. */
typedef struct {
    int pre;                // index of pre-synaptic neuron
    int post;               // index of post-synaptic neuron
    double weight;          // synaptic weight (0 to 1)
    double neurotransmitter; // available neurotransmitter (0 to 1)
    int delay;              // conduction delay (ms)
    int last_pre_spike;     // last spike time of pre neuron
    int last_post_spike;    // last spike time of post neuron
} Synapse;

/* --- Helper function: draw a bar representing the neuron potential --- */
void draw_potential_bar(int row, int col, double potential) {
    int i;
    int filled = (int)((potential / THRESHOLD) * BAR_WIDTH);
    if (filled > BAR_WIDTH) filled = BAR_WIDTH;
    mvprintw(row, col, "[");
    for (i = 0; i < BAR_WIDTH; i++) {
        if (i < filled)
            printw("#");
        else
            printw(" ");
    }
    printw("]");
}

/* --- Main simulation with TUI --- */
int main(void) {
    int t, i;
    
    /* Initialize ncurses */
    initscr();
    cbreak();
    noecho();
    curs_set(0);
    if (has_colors()) {
        start_color();
        init_pair(1, COLOR_WHITE, COLOR_BLACK);  // Normal text
        init_pair(2, COLOR_GREEN, COLOR_BLACK);  // Firing neurons
        init_pair(3, COLOR_RED, COLOR_BLACK);    // Refractory neurons
    }

    /* Seed the random number generator */
    srand((unsigned) time(NULL));

    /* Initialize neurons */
    Neuron neurons[NUM_NEURONS];
    for (i = 0; i < NUM_NEURONS; i++) {
        neurons[i].potential = 0.0;
        neurons[i].refractory = 0;
        neurons[i].spike_count = 0;
    }
    
    /* Define external pulse intervals (ms) for each neuron */
    int pulse_intervals[NUM_NEURONS] = { 100, 90, 110, 80, 120, 70, 130, 60, 140, 50 };

    /* Initialize synapses with random connectivity and delays */
    Synapse synapses[NUM_SYNAPSES];
    for (i = 0; i < NUM_SYNAPSES; i++) {
        synapses[i].pre = rand() % NUM_NEURONS;
        synapses[i].post = rand() % NUM_NEURONS;
        while (synapses[i].post == synapses[i].pre) {
            synapses[i].post = rand() % NUM_NEURONS;
        }
        synapses[i].weight = ((double)rand() / RAND_MAX) * 0.5; // initial weight between 0 and 0.5
        synapses[i].neurotransmitter = NT_INITIAL;
        synapses[i].delay = (rand() % MAX_DELAY) + 1;  // delay between 1 and MAX_DELAY ms
        synapses[i].last_pre_spike = -1000;
        synapses[i].last_post_spike = -1000;
    }
    
    /* Delay buffer for synaptic currents (circular buffer) */
    double delay_buffer[MAX_DELAY][NUM_NEURONS];
    for (int d = 0; d < MAX_DELAY; d++) {
        for (i = 0; i < NUM_NEURONS; i++) {
            delay_buffer[d][i] = 0.0;
        }
    }
    int current_delay_index = 0;
    
    /* Array to record whether a neuron fired in the current time step */
    int fired[NUM_NEURONS] = {0};
    
    /* --- Main simulation loop --- */
    for (t = 0; t < SIM_TIME; t += DT) {
        clear();

        /* Step 1: Add delayed synaptic currents from the delay buffer */
        for (i = 0; i < NUM_NEURONS; i++) {
            if (neurons[i].refractory <= 0) {
                neurons[i].potential += delay_buffer[current_delay_index][i];
            }
            delay_buffer[current_delay_index][i] = 0.0;  // clear after applying
        }
        
        /* Step 2: External input via scheduled pulses */
        for (i = 0; i < NUM_NEURONS; i++) {
            if (t % pulse_intervals[i] == 0 && neurons[i].refractory <= 0) {
                neurons[i].potential += 1.0;
            }
        }
        
        /* Step 3: External noise input for background activity */
        for (i = 0; i < NUM_NEURONS; i++) {
            if (rand() % 100 < NOISE_CHANCE && neurons[i].refractory <= 0) {
                neurons[i].potential += NOISE_AMPLITUDE;
            }
        }
        
        /* Step 4: Update neuron states and check for spikes */
        for (i = 0; i < NUM_NEURONS; i++) {
            if (neurons[i].refractory > 0)
                neurons[i].refractory--;
            
            if (neurons[i].potential >= THRESHOLD && neurons[i].refractory <= 0) {
                fired[i] = 1;
                neurons[i].spike_count++;
                neurons[i].potential = RESET_POTENTIAL;
                neurons[i].refractory = REFRACTORY_PERIOD;
            } else {
                fired[i] = 0;
                neurons[i].potential *= LEAK_FACTOR;
            }
        }
        
        /* Step 5: Process synapses: deliver spikes with delay and apply STDP */
        for (i = 0; i < NUM_SYNAPSES; i++) {
            int pre = synapses[i].pre;
            int post = synapses[i].post;
            
            if (fired[pre]) {
                double effective_signal = synapses[i].weight * synapses[i].neurotransmitter;
                int target_index = (current_delay_index + synapses[i].delay) % MAX_DELAY;
                delay_buffer[target_index][post] += effective_signal;
                synapses[i].neurotransmitter = 0.0;
                synapses[i].last_pre_spike = t;
            }
            if (fired[post]) {
                synapses[i].last_post_spike = t;
            }
            
            /* Recover neurotransmitter */
            if (synapses[i].neurotransmitter < NT_INITIAL) {
                synapses[i].neurotransmitter += NT_RECOVERY_RATE;
                if (synapses[i].neurotransmitter > NT_INITIAL)
                    synapses[i].neurotransmitter = NT_INITIAL;
            }
            
            /* Apply STDP rules based on spike timing differences */
            int dt_pre_post = synapses[i].last_post_spike - synapses[i].last_pre_spike;
            if (dt_pre_post > 0 && dt_pre_post < STDP_WINDOW) {
                double dw = A_PLUS * exp(-((double)dt_pre_post) / TAU_PLUS);
                synapses[i].weight += dw;
            }
            int dt_post_pre = synapses[i].last_pre_spike - synapses[i].last_post_spike;
            if (dt_post_pre > 0 && dt_post_pre < STDP_WINDOW) {
                double dw = A_MINUS * exp(-((double)dt_post_pre) / TAU_MINUS);
                synapses[i].weight += dw;
            }
            if (synapses[i].weight < 0)
                synapses[i].weight = 0;
            if (synapses[i].weight > 1)
                synapses[i].weight = 1;
        }
        
        /* Step 6: Draw the simulation state using ncurses */
        attron(COLOR_PAIR(1));
        mvprintw(0, 0, "Spiking Neural Network Simulation (Time: %d ms)", t);
        mvprintw(1, 0, "---------------------------------------------------");
        
        for (i = 0; i < NUM_NEURONS; i++) {
            int row = 3 + i;
            mvprintw(row, 0, "Neuron %2d: ", i);
            if (fired[i]) {
                attron(COLOR_PAIR(2));
                printw("FIRING! ");
                attroff(COLOR_PAIR(2));
            } else if (neurons[i].refractory > 0) {
                attron(COLOR_PAIR(3));
                printw("Refractory ");
                attroff(COLOR_PAIR(3));
            } else {
                printw("         ");
            }
            draw_potential_bar(row, 20, neurons[i].potential);
            mvprintw(row, 42, " V=%.2f, Spikes=%d", neurons[i].potential, neurons[i].spike_count);
        }
        
        mvprintw(3 + NUM_NEURONS + 1, 0, "Press 'q' to quit early.");
        refresh();
        
        timeout(VISUAL_DELAY);
        int ch = getch();
        if (ch == 'q' || ch == 'Q')
            break;
        
        /* Advance the circular delay buffer index */
        current_delay_index = (current_delay_index + 1) % MAX_DELAY;
    }
    
    /* --- End-of-Simulation Summary --- */
    clear();
    mvprintw(0, 0, "Simulation Complete!\n\n");
    mvprintw(1, 0, "Summary Statistics:\n");
    for (i = 0; i < NUM_NEURONS; i++) {
        mvprintw(3 + i, 0, "Neuron %2d fired %d times, External pulse interval: %d ms",
                 i, neurons[i].spike_count, pulse_intervals[i]);
    }
    
    mvprintw(3 + NUM_NEURONS + 2, 0, "Explanation:");
    mvprintw(3 + NUM_NEURONS + 3, 0,
             "This simulation demonstrates a spiking neural network where each neuron");
    mvprintw(3 + NUM_NEURONS + 4, 0,
             "receives scheduled external pulses, random noise input, and delayed synaptic");
    mvprintw(3 + NUM_NEURONS + 5, 0,
             "inputs. Neurons integrate these inputs until they reach a threshold, fire,");
    mvprintw(3 + NUM_NEURONS + 6, 0,
             "reset, and enter a refractory period. STDP and neurotransmitter dynamics");
    mvprintw(3 + NUM_NEURONS + 7, 0,
             "further modulate the network behavior.");
    mvprintw(3 + NUM_NEURONS + 9, 0, "Press any key to exit.");
    refresh();
    timeout(-1);
    getch();
    
    endwin();
    return 0;
}