Description of problem
I am running a custom brian2 network inside a genetic algorithm. It works when everything is in the same python script, but after refactoring to be across multiple modular scripts that import functions from others as needed, I get the following error, which I cannot find referenced anywhere on the internet:
Error encountered with object named ‘readout_layer_subexpression_update’.
Note that removing the readout layer from the model causes the same error with the reservoir layer, which is declared before the readout layer. I would really appreciate help as I have tried to fix this for many hours and am on a deadline. I reluctantly posted here as a last resort.
Minimal code to reproduce problem
I am not sure this is possible to provide since my simulation consists of 10 separate python scripts with advanced algorithms, functionality, and data structures (not clear how to condense to a minimal reproducible example). If necessary, I am willing to share some more code. However, the following function is likely the culprit, called in the function following it:
def create_neuron_group(layer_size, model, model_type, parameters, name="neurongroup*"):
"""Create a neuron group with excitatory/inhibitory subgroups"""
if "reservoir" in name:
neurone_model=model["reservoir_neurone_model"]
elif "readout" in name:
neurone_model=model["readout_neurone_model"]
else:
neurone_model=model["neurone_model"]
if model_type == "LIF":
G = NeuronGroup(
N=layer_size,
model=neurone_model,
method="euler",
threshold=f"V>{parameters['theta']/mV}*mV",
reset=model["reset_script"],
refractory=parameters["refractory"],
name=name
)
G.R = parameters["resistance"]
G.tau_v = parameters["tau_v"]
G.baseline_rate = parameters["baseline_rate"]
# Split into subgroups
n_inhibitory = round(parameters["proportion_inhibitory"] * layer_size)
G_i = G[:n_inhibitory]
G_e = G[n_inhibitory:]
if "reservoir" in name:
G.Vy_readout_ee_synapses = 0*volt
G.Vy_readout_ei_synapses = 0*volt
G.Vy_readout_ie_synapses = 0*volt
G.Vy_readout_ii_synapses = 0*volt
G.Vy_readout_synapses_efferent = 0*volt
G.Vy_train_synapses_positive = 0*volt
G.Vy_train_synapses_negative = 0*volt
elif "readout" in name:
G.Vy_reservoir_ee_synapses = 0*volt
G.Vy_reservoir_ei_synapses = 0*volt
G.Vy_reservoir_ie_synapses = 0*volt
G.Vy_reservoir_ii_synapses = 0*volt
G.Vy_input_synapses_positive = 0*volt
G.Vy_input_synapses_negative = 0*volt
print(("network_builder line 94 reached", G), flush=True)
return G, G_e, G_i
def build_and_run_network(genome, model_config):
"""Build and run the complete network - main simulation function"""
from models import get_model_equations
from data_loader import prepare_all_data
import warnings
start_scope()
#clear_cache('cython')
prefs.codegen.target = "numpy"
set_device('runtime')
#prefs.codegen.target = "cython" # Enable C++ compilation for speed
#warnings.filterwarnings('ignore')
# Get model equations and prepare data
model = get_model_equations(genome['model_type'])
data = prepare_all_data(model_config=model_config)
parameters = decode_genome_to_parameters(genome, model_config)
# Make data available in Brian2 namespace - not called directly
input_data_positive = data['input_data']['positive']
input_data_negative = data['input_data']['negative']
train_data_positive = data['train_data']['positive']
train_data_negative = data['train_data']['negative']
#globals()['input_data_positive'] = data['input_data']['positive']
#globals()['input_data_negative'] = data['input_data']['negative']
#globals()['train_data_positive'] = data['train_data']['positive']
#globals()['train_data_negative'] = data['train_data']['negative']
namespace = {
# Data arrays
'input_data_positive': input_data_positive,
'input_data_negative': input_data_negative,
'train_data_positive': train_data_positive,
'train_data_negative': train_data_negative,
# Common units and constants
'mV': mV,
'ms': ms,
'Hz': Hz,
'pA': pA,
'volt': volt,
'amp': amp,
'second': second,
'ohm': ohm,
'Mohm': Mohm,
# Parameters that might be referenced in equations
'frequency_factor': parameters['frequency_factor'],
#'baseline_rate': parameters.get('baseline_rate', 0.01),
# Mathematical functions that might be used in equations
'rand': np.random.rand,
'randn': np.random.randn,
'exp': np.exp,
'log': np.log,
'sqrt': np.sqrt,
'sin': np.sin,
'cos': np.cos,
'clip': np.clip,
'poisson': np.random.poisson,
# Additional parameters for equations
'tau_v': parameters.get('tau_v', 5*ms),
'theta': parameters.get('theta', 3*mV),
'refractory': parameters.get('refractory', 0.15*ms),
'resistance': parameters.get('resistance', 1*Mohm),
# For Izhikevich model
'a': parameters.get('a', 0.02),
'b': parameters.get('b', 0.2),
'c': parameters.get('c', -65),
'd': parameters.get('d', 8),
# STDP parameters
'w_max': parameters.get('w_max', 1.0),
'pre_trace_tau': parameters.get('pre_trace_tau', 20*ms),
'post_trace_tau': parameters.get('post_trace_tau', 20*ms),
'pre_trace_step': parameters.get('pre_trace_step', 0.1),
'STDP_asymmetry': parameters.get('STDP_asymmetry', 1.0),
# Synaptic parameters
'tau_nt': parameters.get('tau_nt', 3*ms),
'tau_ntx': parameters.get('tau_ntx', [100*ms, 100*ms, 100*ms, 100*ms]),
'u': parameters.get('u', 0.1),
# Weight signs
'ee_w_sign': parameters.get('ee_w_sign', 1.0),
'ei_w_sign': parameters.get('ei_w_sign', 5.0),
'ie_w_sign': parameters.get('ie_w_sign', -5.0),
'ii_w_sign': parameters.get('ii_w_sign', -3.0),
# Connectivity parameters
'input_connectivity': parameters.get('input_connectivity', 0.3),
'reservoir_connectivity': parameters.get('reservoir_connectivity', 0.1),
'readout_connectivity': parameters.get('readout_connectivity', 0.1),
}
# Reservoir Layer
reservoir_layer, reservoir_layer_excitatory, reservoir_layer_inhibitory = create_neuron_group(parameters['reservoir_size'], model, genome['model_type'], parameters, name="reservoir_layer")
reservoir_ee_synapses = create_synapses(reservoir_layer_excitatory, reservoir_layer_excitatory, model, parameters, parameters['reservoir_connectivity'], "i!=j", parameters['ee_w_sign'], 'reservoir_ee_synapses')
reservoir_ei_synapses = create_synapses(reservoir_layer_excitatory, reservoir_layer_inhibitory, model, parameters, parameters['reservoir_connectivity'], "True", parameters['ei_w_sign'], 'reservoir_ei_synapses')
reservoir_ie_synapses = create_synapses(reservoir_layer_inhibitory, reservoir_layer_excitatory, model, parameters, parameters['reservoir_connectivity'], "True", parameters['ie_w_sign'], 'reservoir_ie_synapses')
reservoir_ii_synapses = create_synapses(reservoir_layer_inhibitory, reservoir_layer_inhibitory, model, parameters, parameters['reservoir_connectivity'], "i!=j", parameters['ii_w_sign'], 'reservoir_ii_synapses')
# Input Layer
input_layer_positive = PoissonGroup(N=parameters['input_layer_size'], rates="input_data_positive(t)")
input_layer_negative = PoissonGroup(N=parameters['input_layer_size'], rates="input_data_negative(t)")
input_synapses_positive = create_synapses(input_layer_positive, reservoir_layer_excitatory, model, parameters, parameters['input_connectivity'], "True", parameters['ee_w_sign'], 'input_synapses_positive')
input_synapses_negative = create_synapses(input_layer_negative, reservoir_layer_excitatory, model, parameters, parameters['input_connectivity'], "True", parameters['ee_w_sign'], 'input_synapses_negative')
# Readout Layer # TODO subexpression_update error occurs here
readout_layer, readout_layer_excitatory, readout_layer_inhibitory = create_neuron_group(2*parameters['readout_layer_size'], model, genome['model_type'], parameters, name="readout_layer")
readout_excitatory_size = len(readout_layer_excitatory)
readout_excitatory_positive = readout_layer_excitatory[:round((1-parameters["proportion_inhibitory"])*parameters["readout_layer_size"])]
readout_excitatory_negative = readout_layer_excitatory[round((1-parameters["proportion_inhibitory"])*parameters["readout_layer_size"]):]
readout_ee_synapses = create_synapses(readout_layer_excitatory, readout_layer_excitatory, model, parameters, parameters["reservoir_connectivity"], "i!=j", parameters["ee_w_sign"], name="readout_ee_synapses")
readout_ei_synapses = create_synapses(readout_layer_excitatory, readout_layer_inhibitory, model, parameters, parameters["reservoir_connectivity"], "True", parameters["ei_w_sign"], name="readout_ei_synapses")
readout_ie_synapses = create_synapses(readout_layer_inhibitory, readout_layer_excitatory, model, parameters, parameters["reservoir_connectivity"], "True", parameters["ie_w_sign"], name="readout_ie_synapses")
readout_ii_synapses = create_synapses(readout_layer_inhibitory, readout_layer_inhibitory, model, parameters, parameters["reservoir_connectivity"], "i!=j", parameters["ii_w_sign"], name="readout_ii_synapses")
readout_synapses_efferent = create_synapses(reservoir_layer, readout_layer, model, parameters, parameters["input_connectivity"], "True", parameters["ee_w_sign"], name="readout_synapses_efferent")
# Training layer
train_layer_positive = PoissonGroup(N=parameters["readout_layer_size"], rates="train_data_positive(t)")
train_layer_negative = PoissonGroup(N=parameters["readout_layer_size"], rates="train_data_negative(t)")
train_synapses_positive = create_synapses(train_layer_positive, readout_excitatory_positive, model, parameters, parameters["input_connectivity"], "True", parameters["ee_w_sign"], name="train_synapses_positive")
train_synapses_negative = create_synapses(train_layer_negative, readout_excitatory_negative, model, parameters, parameters["input_connectivity"], "True", parameters["ee_w_sign"], name="train_synapses_negative")
# Split readout excitatory into positive/negative
# Create monitors
spike_monitor_readout = SpikeMonitor(readout_layer, name="spike_monitor_readout")
# Run simulation
train_runtime = data['train_data']['train_size'] * parameters['dt']
test_runtime = parameters['runtime'] - train_runtime
network = Network([
# Reservoir
reservoir_layer,
reservoir_ee_synapses,
reservoir_ei_synapses,
reservoir_ie_synapses,
reservoir_ii_synapses,
# Input
input_layer_positive,
input_layer_negative,
input_synapses_positive,
input_synapses_negative,
# Readout
readout_layer,
readout_excitatory_positive,
readout_excitatory_negative,
readout_ee_synapses,
readout_ei_synapses,
readout_ie_synapses,
readout_ii_synapses,
readout_synapses_efferent,
spike_monitor_readout,
# Training
train_layer_positive,
train_layer_negative,
train_synapses_positive,
train_synapses_negative,
])
# Training phase
print(("network_builder line 360 reached:"), flush=True)
try:
network.run(train_runtime)#, namespace=namespace)
except Exception as e:
exc_type, exc_value, exc_traceback = sys.exc_info()
print("ERRORFLAG:", e.with_traceback(exc_traceback), flush=True)
print(("network_builder line 362 reached"), flush=True)
# Freeze plasticity for testing
#for synapse in synapses.values():
# synapse.plastic = False
network["reservoir_ee_synapses"].plastic = False
network["reservoir_ei_synapses"].plastic = False
network["reservoir_ie_synapses"].plastic = False
network["reservoir_ii_synapses"].plastic = False
network["readout_ee_synapses"].plastic = False
network["readout_ei_synapses"].plastic = False
network["readout_ie_synapses"].plastic = False
network["readout_ii_synapses"].plastic = False
network["readout_synapses_efferent"].plastic = False
network["train_synapses_positive"].plastic = False
network["train_synapses_negative"].plastic = False
# Testing phase
network.run(test_runtime)#, namespace=namespace)
# Return results
return {
#'readout_positive_spikes': network["spike_monitor_readout_positive"],
#'readout_negative_spikes': network["spike_monitor_readout_negative"],
'spike_monitor_readout': network["spike_monitor_readout"],
'series_length': parameters['runtime'] / parameters['dt'],
'readout_size': readout_excitatory_size // 2,
'proportion_inhibitory': parameters['proportion_inhibitory']
}
The network models are as follows (note the Vy_{}_post):
def get_LIF_model():
"""Get Leaky Integrate-and-Fire model equations"""
return {
# Neurone Models
"reservoir_neurone_model" : '''
dV/dt = (R*I-V
+ Vy_reservoir_ee_synapses
+ Vy_reservoir_ei_synapses
+ Vy_reservoir_ie_synapses
+ Vy_reservoir_ii_synapses
+ Vy_input_synapses_positive
+ Vy_input_synapses_negative)/tau_v : volt (unless refractory)
I = poisson(baseline_rate*(frequency_factor/Hz))*pA: amp (constant over dt)
Vy_reservoir_ee_synapses : volt
Vy_reservoir_ei_synapses : volt
Vy_reservoir_ie_synapses : volt
Vy_reservoir_ii_synapses : volt
Vy_input_synapses_positive : volt
Vy_input_synapses_negative : volt
R : ohm
tau_v : second
baseline_rate : 1 (constant)
''',
"readout_neurone_model" : '''
dV/dt = (R*I-V
+ Vy_readout_ee_synapses
+ Vy_readout_ei_synapses
+ Vy_readout_ie_synapses
+ Vy_readout_ii_synapses
+ Vy_readout_synapses_efferent
+ Vy_train_synapses_positive
+ Vy_train_synapses_negative)/tau_v : volt (unless refractory)
I = poisson(baseline_rate*(frequency_factor/Hz))*pA: amp (constant over dt)
Vy_readout_ee_synapses : volt
Vy_readout_ei_synapses : volt
Vy_readout_ie_synapses : volt
Vy_readout_ii_synapses : volt
Vy_readout_synapses_efferent : volt
Vy_train_synapses_positive : volt
Vy_train_synapses_negative : volt
R : ohm
tau_v : second
baseline_rate : 1 (constant)
''',
"neurone_model" : '''
dV/dt = (R*I-V
+ Vy_reservoir_ee_synapses
+ Vy_reservoir_ei_synapses
+ Vy_reservoir_ie_synapses
+ Vy_reservoir_ii_synapses
+ Vy_input_synapses_positive
+ Vy_input_synapses_negative
+ Vy_readout_ee_synapses
+ Vy_readout_ei_synapses
+ Vy_readout_ie_synapses
+ Vy_readout_ii_synapses
+ Vy_readout_synapses_efferent
+ Vy_train_synapses_positive
+ Vy_train_synapses_negative)/tau_v : volt (unless refractory)
I = poisson(baseline_rate*(frequency_factor/Hz))*pA: amp (constant over dt)
Vy_reservoir_ee_synapses : volt
Vy_reservoir_ei_synapses : volt
Vy_reservoir_ie_synapses : volt
Vy_reservoir_ii_synapses : volt
Vy_input_synapses_positive : volt
Vy_input_synapses_negative : volt
Vy_readout_ee_synapses : volt
Vy_readout_ei_synapses : volt
Vy_readout_ie_synapses : volt
Vy_readout_ii_synapses : volt
Vy_readout_synapses_efferent : volt
Vy_train_synapses_positive : volt
Vy_train_synapses_negative : volt
R : ohm
tau_v : second
baseline_rate : 1 (constant)
''',
"reset_script" : '''
V = 0*mV
''',
"synapse_model" : '''
dpre_trace/dt = -pre_trace/pre_trace_tau : 1 (event-driven)
dpost_trace/dt = -post_trace/post_trace_tau : 1 (event-driven)
w : 1
w_max : 1
pre_trace_tau : second
post_trace_tau : second
pre_trace_step : 1
STDP_asymmetry : 1 (constant)
post_trace_step = -pre_trace_step*pre_trace_tau/post_trace_tau*STDP_asymmetry : 1
dy/dt = -y/tau_nt*plastic : 1 (clock-driven)
dx/dt = z/tau_ntx*plastic : 1 (clock-driven)
dz/dt = y/tau_nt*plastic - z/tau_ntx : 1 (clock-driven)
Vy_{}_post = y*w_sign : volt (summed)
u : 1
tau_nt : second
tau_ntx : second
w_sign : volt (constant)
plastic : 1
''',
# Synaptic Models
"presynaptic_script" : '''
V_post += w*w_sign
y += u*x*plastic
x -= u*x*plastic
pre_trace += pre_trace_step
dw = (w_max-w)*post_trace*plastic
w = clip(w+dw, 0, w_max)
''',
"postsynaptic_script" : '''
post_trace += post_trace_step
dw = (w_max-w)*pre_trace*plastic
w = clip(w+dw, 0, w_max)
''',
"rebalance_script" : '''
V_post -= 0*mv
'''
}
What you have already tried
- explicitly setting the namespace
- sending the namespace to every declared
NeuronGroupand/orSynapseobjects - sending the namespace to
network.run() - using magic function configuration (just
run()instead ofnetwork.run()) - writing input data vector
TimedArraystoglobal() - debugging by printing out model f-string formatting (
.format(name)) - removing the readout layer and various synapses
- transplanting functions that worked in the singular python script
- segregating model equations by readout and reservoir
- initialising unused
Vy_{}_postparameters to 0*mv - using
clear_cache(“cython”)and/or removingstart_scope() - finding similar issues (a Google search with “subexpression_update” leads only to a single Brian documentation page about the source code that is unhelpful)
- consulting an LLM at maximum settings
- some other stuff I can’t recall at present
Expected output (if relevant)
- The network runs without an error
Actual output (if relevant)
- Error encountered with object named ‘readout_layer_subexpression_update’.
Full traceback of error (if relevant)
I suspect the “full” traceback is being cut off, and I tried to find the “debug log” that is mentioned in the error output, and also to find where brian2 would store such a log, but I could not. If there is a location where I can find it, I will be happy to retrieve and post it.
Error in simulation: Error encountered with object named ‘readout_layer_subexpression_update’.
Object was created here (most recent call only, full details in debug log):
File ‘[filepath]\network_builder.py’, line 72, in create_neuron_group
G = NeuronGroup(An error occurred when preparing an object. (See above for original error message and traceback.)