import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_absolute_error
from brian2 import *
# Set random seed for reproducibility
np.random.seed(0)
# Define network parameters
n_input = x_train.shape[1]
n_hidden1 = 200
n_hidden2 = 100
n_hidden3 = 50
n_output = 1
timestep = 5 * ms
duration = 100 * ms
Cm = 1 * uF # Membrane capacitance
gbar_Na = 120 * msiemens # Sodium conductance
gbar_K = 36 * msiemens # Potassium conductance
gbar_L = 0.3 * msiemens # Leak conductance
ENa = 50 * mV # Sodium reversal potential
EK = -77 * mV # Potassium reversal potential
EL = -54.4 * mV # Leak reversal potential
v_rest = -65 * mV # Resting membrane potential
# Create a Network object
net = Network()
# Define neuron model equations (Hodgkin-Huxley)
neuron_eqs = '''
dv/dt = (1/Cm) * (I_Na + I_K + I_L - I_input) : volt
I_Na = gbar_Na * m**3 * h * (ENa - v) : amp
I_K = gbar_K * n**4 * (EK - v) : amp
I_L = gbar_L * (EL - v) : amp
I_input : amp
dm/dt = (alpha_m * (1 - m) - beta_m * m) * Hz : 1
dn/dt = (alpha_n * (1 - n) - beta_n * n) * Hz : 1
dh/dt = (alpha_h * (1 - h) - beta_h * h) * Hz : 1
alpha_m = (0.1/mV) * (10*mV) / exprel((v + 25*mV) / (10*mV)) : Hz
beta_m = 4 * exp((v + 50*mV) / (18*mV)) : Hz
alpha_n = (0.01/mV) * (10*mV) / exprel((v + 10*mV) / (10*mV)) : Hz
beta_n = 0.125 * exp(v / (80*mV)) : Hz
alpha_h = 0.07 * exp(v / (20*mV)) : Hz
beta_h = 1 / (exp((v + 30*mV) / (10*mV)) + 1) : Hz
'''
# Define synapse model equations
synapse_eqs = '''
dw/dt = 0 : 1
'''
# Create neuron groups
input_neurons = NeuronGroup(n_input, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
hidden_neurons1 = NeuronGroup(n_hidden1, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
hidden_neurons2 = NeuronGroup(n_hidden2, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
hidden_neurons3 = NeuronGroup(n_hidden3, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
output_neurons = NeuronGroup(n_output, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
# Create synapses
synapse_input_hidden1 = Synapses(input_neurons, hidden_neurons1, on_pre=' v_post += 0.2')
synapse_hidden1_hidden2 = Synapses(hidden_neurons1, hidden_neurons2, on_pre=' v_post += 0.2')
synapse_hidden2_hidden3 = Synapses(hidden_neurons2, hidden_neurons3, on_pre=' v_post += 0.2')
synapse_hidden3_output = Synapses(hidden_neurons3, output_neurons, on_pre=' v_post += 0.2')
# Connect synapses
synapse_input_hidden1.connect()
synapse_hidden1_hidden2.connect()
synapse_hidden2_hidden3.connect()
synapse_hidden3_output.connect()
# Add neuron groups and synapses to the network
net.add(input_neurons, hidden_neurons1, hidden_neurons2, hidden_neurons3, output_neurons)
net.add(synapse_input_hidden1, synapse_hidden1_hidden2, synapse_hidden2_hidden3, synapse_hidden3_output)
# Define input and output variables
inputs = TimedArray(x_train.T, dt=timestep)
outputs = np.zeros(n_output)
# Define monitors
output_monitor = SpikeMonitor(output_neurons)
# Run simulation
net.run(duration)
# Get spikes and decode output
output_spikes = output_monitor.spike_trains()
for i in range(n_output):
output_train = output_spikes[i]
if len(output_train) > 0:
outputs[i] = output_train[-1]
# Print output values
print('Output:', outputs)
# Decode output values
predicted_train = outputs
# Calculate MAE
mae_train = mean_absolute_error(y_train, predicted_train)
print('Train MAE:', mae_train)
# Define input variables for testing
inputs_test = TimedArray(x_test.T, dt=timestep)
outputs_test = np.zeros(n_output)
# Run simulation for testing
run(duration)
# Get spikes and decode output for testing
output_spikes_test = output_monitor.spike_trains()
for i in range(n_output):
output_test = output_spikes_test[i]
if len(output_test) > 0:
outputs_test[i] = output_test[-1]
# Decode output values for testing
predicted_test = outputs_test
# Calculate MAE for testing
mae_test = mean_absolute_error(y_test, predicted_test)
print('Test MAE:', mae_test)
Hi @sakattack85. I can guess your question from the title of your post, but please in the future, ask a complete question and try to reduce your example code to the minimal code necessary.
For example, the following code would lead to the same error message:
from brian2 import *
Cm = 1 * uF # Membrane capacitance
gbar_Na = 120 * msiemens # Sodium conductance
gbar_K = 36 * msiemens # Potassium conductance
gbar_L = 0.3 * msiemens # Leak conductance
ENa = 50 * mV # Sodium reversal potential
EK = -77 * mV # Potassium reversal potential
EL = -54.4 * mV # Leak reversal potential
v_rest = -65 * mV # Resting membrane potential
# Define neuron model equations (Hodgkin-Huxley)
neuron_eqs = '''
dv/dt = (1/Cm) * (I_Na + I_K + I_L - I_input) : volt
I_Na = gbar_Na * m**3 * h * (ENa - v) : amp
I_K = gbar_K * n**4 * (EK - v) : amp
I_L = gbar_L * (EL - v) : amp
I_input : amp
dm/dt = (alpha_m * (1 - m) - beta_m * m) * Hz : 1
dn/dt = (alpha_n * (1 - n) - beta_n * n) * Hz : 1
dh/dt = (alpha_h * (1 - h) - beta_h * h) * Hz : 1
alpha_m = (0.1/mV) * (10*mV) / exprel((v + 25*mV) / (10*mV)) : Hz
beta_m = 4 * exp((v + 50*mV) / (18*mV)) : Hz
alpha_n = (0.01/mV) * (10*mV) / exprel((v + 10*mV) / (10*mV)) : Hz
beta_n = 0.125 * exp(v / (80*mV)) : Hz
alpha_h = 0.07 * exp(v / (20*mV)) : Hz
beta_h = 1 / (exp((v + 30*mV) / (10*mV)) + 1) : Hz
'''
# Create neuron groups
input_neurons = NeuronGroup(1, neuron_eqs, threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
run(0*ms)
The error you are getting complains about the equation alpha_m = (0.1/mV) * (10*mV) / exprel((v + 25*mV) / (10*mV)) : Hz which states to be in Hz (which would be correct for a rate variable), but the stated equation gives a dimensionless value. In principle, you could correct the dimensions of this equation by multiplying it by Hz, but I think in this case the values are rather stated “per ms” (as often the case in papers that do not use units in the equations), i.e. you should divide the right-hand side by ms (also for beta_m and the other rate variables). Conversely, you’ll need to remove the * Hz in the definitions of m, n, and h for consistent dimensions:
neuron_eqs = '''
dv/dt = (1/Cm) * (I_Na + I_K + I_L - I_input) : volt
I_Na = gbar_Na * m**3 * h * (ENa - v) : amp
I_K = gbar_K * n**4 * (EK - v) : amp
I_L = gbar_L * (EL - v) : amp
I_input : amp
dm/dt = (alpha_m * (1 - m) - beta_m * m) : 1
dn/dt = (alpha_n * (1 - n) - beta_n * n) : 1
dh/dt = (alpha_h * (1 - h) - beta_h * h) : 1
alpha_m = (0.1/mV) * (10*mV) / exprel((v + 25*mV) / (10*mV))/ms : Hz
beta_m = 4 * exp((v + 50*mV) / (18*mV))/ms : Hz
alpha_n = (0.01/mV) * (10*mV) / exprel((v + 10*mV) / (10*mV))/ms : Hz
beta_n = 0.125 * exp(v / (80*mV))/ms : Hz
alpha_h = 0.07 * exp(v / (20*mV))/ms : Hz
beta_h = 1 / (exp((v + 30*mV) / (10*mV)) + 1)/ms : Hz
'''
Hope that helps!
PS: When you paste code, please include it in triple backticks like this:
```
# Here comes the code
from brian2 import *
eqs = '''something'''
# ...
```
Thanks for your reply! Now the problem is that I get only one output and my y_test has [200,1] values. It seems that it only simulates the network once, for the first input . Can you checck the code above? What i am doing is that I have x_train (800, 44) etc. and i try to do regression an calculate val_mae
e.g.
‘’‘’
Generate random input data
input_data = np.random.random((1000, 44))
output_data = np.random.random((1000, 1))
Split the data into training and validation sets
x_train, y_train, x_train, y_test = train_test_split(
input_data, output_data, test_size=0.2, random_state=1
)
‘’’
Hi @sakattack85 I am afraid I don’t yet understand the details of your model. You are defining:
inputs = TimedArray(x_train.T, dt=timestep)
But inputs is never used in your model, and I_input is never set to anything, so your network never gets any input.
In general, I’d strongly suggest to start simple (e.g. with just the neuron model and its input) and verify that the basics work before going to training/testing, multiple layers, etc.
here it is
# Define input variables for testing
inputs_test = TimedArray(x_test, dt=timestep)
outputs_test = np.zeros(800, n_output)) # Initialize array to store output spikes
# Run simulation for testing
for i in range(len(inputs_test)):
input_neurons.I = inputs_test(x_test[i, :], dt=1*ms) * amp # Set input current for each sample
run(duration)
# Get spikes and decode output for each output neuron
output_spikes_test = output_monitor.spike_trains()
for j in range(n_output):
output_test = output_spikes_test[j]
if len(output_test) > 0:
outputs_test[i, j] = output_test[-1]
# Decode output values for testing
predicted_test = outputs_test
Hi again. I am afraid this is not how to use a TimedArray, please have a look at the documentation and the tutorial.
Does your stimulus actually vary over time (for each run)? From your code above, it rather looks as if you have a constant stimulus (a different one for each neuron) for each of the samples. In that case you could either do without TimedArray by assigning input_neurons.I = x_test[i, :] for each run, or you could do a single long run (more efficient) where you use a TimedArray to change the input across samples. For the latter, you’d write something like this:
inputs = TimedArray(x_test.T, dt=duration)
neuron_eqs = '''
...
I_input = inputs(t, i)*nA : amp
'''