Library available on PyPI
✅ Installable directly with:

pip install snn-reservoir-py==2.0.0

This class implements a Spiking Neural Network (SNN) model based on LIF neurons and a reservoir architecture, with either a small-world topology generated using the Watts-Strogatz model, or a random-uniform topology generated using an Erdős–Rényi model. The typical operational flow is as follows:

1. Instantiate an SNN network.
2. Provide a binary matrix input_spike_times with shape (num_input_neurons, time_steps) where 1 indicates stimulation of an input neuron at a specific time.
3. At each timestep, a current is injected into the activated input neurons.

⚠️ Input neurons correspond to the first input_spike_times.shape[0] indices.
⚠️ Output neurons are randomly selected among hidden neurons (non-input), unless you set them manually via set_output_neurons().

Optionally, the model supports STDP-based synaptic plasticity, synapse pruning, and an auto-tuning membrane threshold mechanism for homeostasis.

STDP is implemented using pre- and post-synaptic traces that decay exponentially:

x_pre(t+1) = x_pre(t) * exp(-Δt / tau_plus)
x_post(t+1) = x_post(t) * exp(-Δt / tau_minus)

When a spike occurs, the corresponding trace is incremented by 1.
Synaptic weights are updated according to either:

• Multiplicative rule
  LTP: w ← w + η * A_plus * x_pre * (1 - w/W_max)
  LTD: w ← w * (1 - (η * A_minus / W_max) * x_post)

• Additive rule
  LTP: w ← w + η * A_plus * x_pre
  LTD: w ← w - η * A_minus * x_post

where tau_plus/tau_minus are decay constants, A_plus/A_minus are learning amplitudes, η is the learning rate, and W_max is the maximum synaptic weight.

A minimal end-to-end example is available here:
➡️ Example repo: https://github.com/RuggeroFreddi/Example-SNNpy

• (self, trace_tau = 10, reset_trace=False)
  Runs the simulation for the full duration and returns a 2D NumPy array (np.ndarray):
  Params: trace_tau (>0) is the exponential decay time constant (in time steps) for the per-neuron trace; reset_trace=True zeros the trace at the start of the run. Returns: a binary matrix of shape [time x output_neurons], where each row is a timestep and each column an output neuron (1 = spike, 0 = no spike).

• get_spike_time_lists_output()
  Returns a list of lists: each sublist contains the timesteps at which each output neuron generated a spike.

• set_input_spike_times(input_spike_times)
  Sets the binary input spike matrix (shape [input_neurons (0/1) x time]).
  Automatically updates the simulation duration (duration) if not already set.

• set_membrane_potentials(membrane_potentials)
  Sets the initial membrane potentials for each neuron.
  The array must be 1D with length equal to the number of neurons (num_neurons).

• set_output_neurons(indices)
  Sets the output neurons by specifying their indices (array of integers).

• reset()
  Resets the internal state of the network to allow for a new simulation.
  Restores initial membrane potentials, clears spike matrix, spike count, and refractory timers, while preserving all other parameters (adjacency_matrix, output_neurons, etc.).

These methods implement synaptic plasticity (STDP) and homeostatic regulation (threshold auto-tuning, weight rescaling, pruning).

• rescale_synaptic_weights_to_mean(target_mean)
  Scales all non-zero synaptic weights so that their mean equals the given target_mean.
  Preserves the relative distribution of weights while adapting their global strength.
  (Homeostatic mechanism: keeps synaptic activity balanced.)

• compute_critical_threshold(use_abs_weights=False) / apply_critical_threshold(use_abs_weights=False)
  Compute / apply a recommended firing threshold based on:
  θ = (mean synaptic weight × mean in-degree) + (2 × input_mean_current × refractory_period).
  (see A Mean-Field Approach to Criticality in Spiking Neural Networks for Reservoir Computing)
  (Homeostatic mechanism: adapts membrane threshold automatically to balance activity.)

• disable_stdp()
  Disables STDP learning and clears all internal traces and buffers.
  Useful to run the network in static mode after plasticity has been tested.

• prune(fraction)
  Removes the weakest fraction of synapses (by absolute weight).
  For example, prune(0.2) deletes the 20% smallest weights.
  After pruning, statistics (mean, variance, in-degree) are updated automatically.

• reset_synaptic_weights(mean, std=0.1)
  Resets existing synaptic weights by re-sampling them from a normal distribution with the specified mean and standard deviation (std * mean). Can be combined with STDP to restart learning from a fresh initialization.

Each feature extraction method returns a one-dimensional NumPy array (np.ndarray) containing one feature per output neuron, except extract_features_from_spikes(), which returns a dict of feature arrays keyed by feature name (each value is a 1D array over output neurons).

• extract_features_from_spikes() Extracts all the main features from the output neurons (spike count, entropy, ISI, etc.)

• get_spike_counts()
  Total number of spikes per neuron.

• get_mean_spike_times(), get_first_spike_times(), get_last_spike_times()
  Temporal statistics on spike timing.

• get_mean_isi_per_neuron(), get_isi_variance_per_neuron()
  Statistics on inter-spike intervals (ISI).

• get_spike_entropy_per_neuron()
  Entropy of the temporal distribution of spikes.

• get_spike_rates()
  Average spike rate (spikes per time unit).

• get_autocorrelation_first_lag()
  Autocorrelation at lag 1 (indicator of regularity).

• get_burstiness()
  Burstiness index (ISI variability).

• get_spike_symmetry()
  Spike symmetry between the first and second half of the timeline.

• get_spike_histogram_moments()
  Mean, skewness, and kurtosis of the spike time histogram.

• get_burst_counts()
  Number of bursts per neuron. A burst is defined as a contiguous sequence of spikes (value 1) in the neuron's spike train. Useful for assessing the frequency of temporally clustered activity.

• calculate_mean_isi()
  Returns the average ISI across all neurons.

• reset_synaptic_weights(mean, std)
  Regenerates all existing synaptic weights from a normal distribution.

• save_topology(), load_topology()
  Save/load the synaptic weight matrix (adjacency_matrix) to/from disk in .npz format (SciPy sparse matrix).

• save_membrane_potentials(), load_membrane_potentials()
  Save/load membrane potentials in .npy format.

• save_output_neurons(), load_output_neurons()
  Save/load the indices of the output neurons in .npy format.

• set_topology(sparse_matrix), get_topology()
  Set or retrieve the synaptic weight matrix.

• set_membrane_potentials(array), get_membrane_potentials()
  Set or retrieve the membrane potentials.

• set_output_neurons(indices), get_output_neurons()
  Set or retrieve the indices of the output neurons.

ℹ️ Note: The methods load_topology(), load_output_neurons(), and load_membrane_potentials()
are also available as standalone functions, useful for loading data before initializing the SNN object.

from snnpy import SNN, SimulationParams
import numpy as np

# Input spike train (input neurons x duration)
my_input = np.random.randint(0, 2, size=(50, 500), dtype=np.uint8)

# Parameter configuration
params = SimulationParams(
    num_neurons=2000,
    num_output_neurons=35,
    input_spike_times=my_input,
    leak_coefficient=1 / 10000,
    refractory_period=2,
    membrane_threshold=2.0,
    is_random_uniform=False,
    small_world_graph_p=0.2,
    small_world_graph_k=int(0.10 * 2000 * 2),
    mean_weight=0.00745167232 * 1.05
)

# Create and simulate the network
snn = SNN(params)
output = snn.simulate()

# Extract temporal features from output neurons
features = snn.extract_features_from_spikes()

print("Output shape:", output.shape)

print("Available feature keys:", list(features.keys()))
print("Spike counts shape:", features["spike_counts"].shape)