'Backpropagation and the brain' realized in cortical error neuron microcircuits.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § 2. Theory › 2.8. Synaptic plasticity ↔ numpy_model/src/microcircuit.py, lines 1973–2036 · score 0.64 · pre synaptic rate, post synaptic, synapses, computational, microcircuit, plasticity
[2] § 6. Methods › 6.5. Simulations › 6.5.1. Cart pole and DMS task. ↔ numpy_model/cart-pole plot.ipynb, lines 85–111 · score 0.59 · penalty log, cart pole, epoch, reset, dt

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 2,444 lines · 86 KB · no license · 1 match

import numpy as np
import copy
import inspect
from scipy import signal
#import matplotlib.pyplot as plt
import logging
import os.path as path
# define activation functions
def linear(x, slope=1.0, offset=0.0):
x = x-offset
return np.array(x*slope)
def d_linear(x, slope=1.0, offset=0.0):
x = x-offset
return np.ones_like(x)*slope
def relu(x, slope=1.0, offset=0.0):
# offset is a bias
return np.maximum(x-offset, 0, np.array(x))*slope
def d_relu(x, slope=1.0, offset=0.0):
return np.heaviside(x-offset, 0)*slope
def soft_relu(x, slope=1.0, offset=0.0):
x = x-offset
exp = np.exp(-10.0*x)
return np.log((1 + exp)/exp)/10.0*slope
def d_soft_relu(x, slope=1.0, offset=0.0):
x = x-offset
return 1 / (1 + np.exp(-10.0*x))*slope
def logistic(x, slope=1.0, offset=0.0):
return 1/(1 + np.exp(-(x-offset)*slope))
def d_logistic(x, slope=1.0, offset=0.0):
y = logistic(x, slope=slope, offset=offset)
return y * (1.0 - y)
def tanh(x, slope=1.0, offset=0.0):
return np.tanh((x-offset)*slope)
def d_tanh(x, slope=1.0, offset=0.0):
y = tanh(x*slope, offset=offset)
return 1 - y**2
# def tanh_offset(x, slope=1.0):
# # tanh offset by +1 to produce positive output only
# return np.tanh(x*slope) + 1
# def d_tanh_offset(x, slope=1.0):
# y = tanh(x*slope)
# return 1 - y**2
def hard_sigmoid(x, slope=1.0, offset=0.0):
# see defintion at torch.nn.Hardsigmoid
x = x-offset
return np.maximum(0, np.minimum(1, x/6 + 1/2), np.array(x))*slope
def d_hard_sigmoid(x, slope=1.0, offset=0.0):
x = x-offset
return np.heaviside(x/6 + 1/2, 0) * np.heaviside(1 - (x/6 + 1/2), 0)*slope
# define dict between activations and derivatives
dict_d_activation = {
"linear": d_linear,
"relu": d_relu,
"soft_relu": d_soft_relu,
"logistic": d_logistic,
"tanh": d_tanh,
# "tanh_offset": d_tanh_offset,
"hard_sigmoid": d_hard_sigmoid
}
# cosine similarity between tensors
def cos_sim(A, B):
if A.ndim == 1 and B.ndim == 1:
return A.T @ B / np.linalg.norm(A) / np.linalg.norm(B)
else:
return np.trace(A.T @ B) / np.linalg.norm(A) / np.linalg.norm(B)
def dist(A, B):
return np.linalg.norm(A-B)
def deg(cos, precision=1e-6):
if np.abs(1.0 - cos) < precision:
return 0.0
else:
# calculates angle in deg from cosine
return np.arccos(cos) * 180 / np.pi
def deepcopy_array(array):
""" makes a deep copy of an array of np-arrays """
out = [nparray.copy() for nparray in array]
return out.copy()
def MSE(output, target):
return np.linalg.norm(output - target)**2
def moving_average(x, w):
return np.convolve(x, np.ones(w), 'valid') / w
def label_to_onehot(label, classes):
onehot = np.zeros(classes)
onehot[label] = 1
return onehot
def unison_shuffled_copies(a, b, rng):
"""
Shuffles two np arrays in the same way
by mtrw at SO
"""
assert len(a) == len(b)
p = rng.permutation(len(a))
return a[p], b[p]
def convert_eta_to_matrix_for_skip(eta_fw, WPP, layers):
"""
Converts a learning rate array of shape (# of WPP x # of WPP)
to an array of shape (# of neurons # of neurons) == shape of WPP
"""
# pad eta_fw with zeros to align shape with WPP
eta_fw = np.array(eta_fw)
eta_fw = np.pad(eta_fw, ((1,0),(0,1)), 'constant')
# repeat to match number of neurons in each layer
eta_fw = np.repeat(eta_fw, layers, axis=0)
eta_fw = np.repeat(eta_fw, layers, axis=1)
return eta_fw
def init_nonhierarchical_WPP(WPP_hierarchical_range, layers, WPP_skip_connection_range):
"""
Generates a WPP matrix (for class errormc_model) which connects all representation units.
WPP_hierarchical_range: uniform range for init values (hierarchical feedback)
WPP_skip_connection_range: uniform range for init values (skip connections)
"""
WPP_init = np.zeros(shape=(sum(layers), sum(layers)))
# we populate WPP_init the same way as BPP_init (below), then transpose
for idx, _ in enumerate(layers[:-1]):
lower_row = np.sum(layers[:idx], dtype=int)
upper_row = np.sum(layers[:idx+1])
lower_col = np.sum(layers[:idx+1])
upper_col = np.sum(layers[:idx+2])
# generate hierarchical part of WPP
WPP_init[lower_row:upper_row, lower_col:upper_col] = np.random.uniform(WPP_hierarchical_range[0], \
WPP_hierarchical_range[1], size=(layers[idx], layers[idx+1]))
# generate skip connections of WPP
WPP_init[:lower_col, upper_col:] = np.random.uniform(WPP_skip_connection_range[0], \
WPP_skip_connection_range[1], size=(WPP_init[:lower_col, upper_col:].shape))
return [WPP_init.T]
def init_nonhierarchical_BII(BII_hierarchical_range, error_layers, BII_skip_connection_range):
"""
Generates a BII matrix (for class errormc_model) which connects all error units.
BII_hierarchical_range: uniform range for init values (hierarchical feedback)
BII_skip_connection_range: uniform range for init values (skip connections)
"""
# filter out "none" layers (no error units in first layer)
BII_init = np.zeros(shape=(sum(filter(None, error_layers)), sum(filter(None, error_layers))))
for idx, _ in enumerate(error_layers[:-1]):
lower_row = np.sum(error_layers[:idx], dtype=int)
upper_row = np.sum(error_layers[:idx+1])
lower_col = np.sum(error_layers[:idx+1])
upper_col = np.sum(error_layers[:idx+2])
# generate hierarchical part of BII
BII_init[lower_row:upper_row, lower_col:upper_col] = np.random.uniform(BII_hierarchical_range[0], \
BII_hierarchical_range[1], size=(error_layers[idx], error_layers[idx+1]))
# generate skip connections of BII
BII_init[:lower_col, upper_col:] = np.random.uniform(BII_skip_connection_range[0], \
BII_skip_connection_range[1], size=(BII_init[:lower_col, upper_col:].shape))
return [BII_init]
def init_viz_cort_conn(layers, error_layers,
init_WPP_range, init_BII_range):
# load connectivity matrix based on Markov et al., 2014
cwd = path.dirname(path.abspath(__file__))
conn = np.genfromtxt(cwd + '/connectivity.csv', delimiter=',').T
if len(layers[1:]) > len(conn) or len(error_layers[1:]) > len(conn):
raise ValueError(f"Too many layers, connectivity only implemented for {len(conn)} layers.")
# rescale skip connections using hierarchical connections and connectivity matrix
assert len(layers) == len(error_layers)
fw_rescale_mask = conn.copy()
# insert row/column for input
fw_rescale_mask = np.pad(fw_rescale_mask,((1, 0), (1, 0)), mode='constant')
fw_rescale_mask[1,0] = 1.0
WPP_init = np.zeros(shape=(sum(layers), sum(layers)))
for idx, _ in enumerate(layers):
for idy, _ in enumerate(layers):
lower_row = np.sum(layers[:idx], dtype=int)
upper_row = np.sum(layers[:idx+1], dtype=int)
lower_col = np.sum(layers[:idy], dtype=int)
upper_col = np.sum(layers[:idy+1], dtype=int)
weight_lower = init_WPP_range[0] * fw_rescale_mask[idx,idy]
weight_upper = init_WPP_range[1] * fw_rescale_mask[idx,idy]
# generate hierarchical part of WPP
WPP_init[lower_row:upper_row, lower_col:upper_col] = np.random.uniform(weight_lower, \
weight_upper, size=(layers[idx], layers[idy]))
# ignore first entry (no backprojection from V1)
error_layers = error_layers[1:]
# no modifications for input
bw_rescale_mask = conn.copy()
BII_init = np.zeros(shape=(sum(error_layers), sum(error_layers)))
for idx, _ in enumerate(error_layers):
for idy, _ in enumerate(error_layers):
lower_row = np.sum(error_layers[:idx], dtype=int)
upper_row = np.sum(error_layers[:idx+1], dtype=int)
lower_col = np.sum(error_layers[:idy], dtype=int)
upper_col = np.sum(error_layers[:idy+1], dtype=int)
weight_lower = init_BII_range[0] * bw_rescale_mask[idx,idy]
weight_upper = init_BII_range[1] * bw_rescale_mask[idx,idy]
# generate hierarchical part of WPP
BII_init[lower_row:upper_row, lower_col:upper_col] = np.random.uniform(weight_lower, \
weight_upper, size=(error_layers[idx], error_layers[idy]))
WPP_init = np.tril(WPP_init)
BII_init = np.triu(BII_init)
return [WPP_init], [BII_init]
# create matrix with list of matrices on diagonal
# user7328723 @ SO
def diag_mat(rem=[], result=np.empty((0, 0))):
if not rem:
return result
m = rem.pop(0)
result = np.block(
[
[result, np.zeros((result.shape[0], m.shape[1]))],
[np.zeros((m.shape[0], result.shape[1])), m],
]
)
return diag_mat(rem, result)
def convert_learning_rate_for_skip(eta, method):
"""
Takes a list of learning rates, e.g. eta=[1,1,1],
and converts it into a 2d array needed for skip connections
by using three different rules:
- fill_diag: uses eta to fill only the diagonal, all others = 0
- fill_all: all learning rates *to* a given layer are equal
- fill_scaled: uses viz cx conn to scale learning rate
Returns an array of entries:
[[eta10, 0, 0, ...],
[eta20, eta21, 0, ...],
[eta30, eta31, eta32, ...]]
"""
assert method in ["fill_diag", "fill_all", "fill_scaled"], \
f"Unknown conversion method for learning rate (not in 'fill_diag', 'fill_all', 'fill_scaled')"
logging.info(f"Adapting eta_fw according to {method}")
eta = np.array(eta)
# converts the learning
if method == 'fill_diag':
eta = np.diag(eta)
elif method == 'fill_all':
eta = np.repeat(eta[:,np.newaxis], len(eta), axis=1)
elif method == 'fill_scaled':
cwd = path.dirname(path.abspath(__file__))
conn = np.genfromtxt(cwd + '/connectivity.csv', delimiter=',').T
# insert row/column for input
conn = np.pad(conn,((1, 0), (1, 0)), mode='constant')
conn[1,0] = 1.0
conn = np.tril(conn[1:][:len(eta),:len(eta)])
eta = np.diag(eta) @ conn
return np.tril(eta)
class base_model:
""" This class implements a generic microcircuit model """
def __init__(self, bw_connection_mode, dWPP_use_activation, dt, Tpres,
model, activation, layers, uP_init, uI_init, WPP_init,
WIP_init, BPP_init, BPI_init, gl, gden, gbas, gapi, gnI,
gntgt, eta_fw, eta_bw, eta_PI, eta_IP, seed=123, WT_noise=0.0):
self.rng = np.random.RandomState(seed)
# connection_mode: skip or layered
self.bw_connection_mode = bw_connection_mode
# whether to use activation in updates of WPP
self.dWPP_use_activation = dWPP_use_activation
self.model = model # FA, BP or PBP
self.layers = layers
self.uP = deepcopy_array(uP_init)
self.uI = deepcopy_array(uI_init)
# we also set up a buffer of voltages,
# which corresponds to the value at the last time step
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
# if a list of activations has been passed, use it
if isinstance(activation, list):
self.activation = activation
# else, set same activation for all layers
else:
self.activation = [activation for layer in layers[1:]]
self.d_activation = [dict_d_activation[activation.__name__] for activation in self.activation]
# whether the target provided should be a rate or voltage
self.rate_target = False
# define the compartment voltages
self.vbas = [np.zeros_like(uP) for uP in self.uP]
self.vden = [np.zeros_like(uI) for uI in self.uI]
self.vapi = [np.zeros_like(uP) for uP in self.uP[:-1]]
# and make copies
self.vbas_old = deepcopy_array(self.vbas)
self.vden_old = deepcopy_array(self.vden)
self.vapi_old = deepcopy_array(self.vapi)
self.WPP = deepcopy_array(WPP_init)
self.WIP = deepcopy_array(WIP_init)
self.BPP = deepcopy_array(BPP_init)
self.BPI = deepcopy_array(BPI_init)
self.dWPP = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWIP = [np.zeros(shape=WIP.shape) for WIP in self.WIP]
self.dBPP = [np.zeros(shape=BPP.shape) for BPP in self.BPP]
self.dBPI = [np.zeros(shape=BPI.shape) for BPI in self.BPI]
# set perfect transpose weights for BP
if self.model == "BP":
self.set_weights(BPP = [WPP.T for WPP in self.WPP[1:]])
# noise level in setting transpose weights
self.WT_noise = WT_noise
# noise matrix (calculated once)
self.BPP_noise = [self.rng.uniform(-self.WT_noise, self.WT_noise, size=BPP.shape) for BPP in self.BPP]
for i, _ in enumerate(self.BPP):
# add noise
self.BPP[i] += self.BPP_noise[i]
self.gl = gl
self.gden = gden
self.gbas = gbas
self.gapi = gapi
self.gnI = gnI
self.gntgt = gntgt
self.Time = 0 # initialize a model timer
self.dt = dt
self.Tpres = Tpres
self.taueffP, self.taueffP_notgt, self.taueffI = self.calc_taueff()
# learning rates
self.eta_fw = eta_fw
self.eta_bw = eta_bw
self.eta_IP = eta_IP
self.eta_PI = eta_PI
# calculate lookahead
self.uP_breve = [self.prospective_voltage(
self.uP[i],
self.uP_old[i],
self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [self.prospective_voltage(
self.uI[i],
self.uI_old[i],
self.taueffI[i]) for i in range(len(self.uI))]
# calculate rate of lookahead: phi(ubreve)
self.rP_breve = [self.activation[i](self.uP_breve[i])
for i in range(len(self.uP_breve))]
try:
self.rI_breve = [self.activation[i+1](self.uI_breve[i])
for i in range(len(self.uI_breve))]
except (IndexError, ValueError):
logging.info("rI_breve not defined")
self.r0 = np.zeros(self.layers[0])
self.dWPP_r_low_pass = False
self.dWPP_post_low_pass = False
def init_record(self, rec_per_steps=1, rec_MSE=False,
rec_error=False, rec_input=False, rec_target=False,
rec_WPP=False, rec_WIP=False, rec_BPP=False,
rec_BII=False,
rec_BPI=False, rec_dWPP=False, rec_dWIP=False,
rec_dBPP=False, rec_dBPI=False, rec_dBII=False, rec_uP=False,
rec_uP_breve=False, rec_rP_breve=False,
rec_rP_breve_HI=False, rec_uI=False,
rec_uI_breve=False, rec_rI_breve=False, rec_vbas=False,
rec_vapi=False, rec_vapi_noise=False, rec_noise=False,
rec_epsilon=False, rec_epsilon_LO=False,
rec_lat_mismatch=False):
# records the values of the variables given in var_array
# e.g. WPP, BPP, uP_breve
# rec_per_steps sets after how many steps data is recorded
if rec_MSE:
self.MSE_time_series = []
self.MSE_val_time_series = []
self.MSE_test_time_series = []
if rec_error:
self.error_time_series = []
if rec_input:
self.input_time_series = []
if rec_target:
self.target_time_series = []
if rec_WPP:
self.WPP_time_series = []
if rec_WIP:
self.WIP_time_series = []
if rec_BPP:
self.BPP_time_series = []
if rec_BII:
self.BII_time_series = []
if rec_BPI:
self.BPI_time_series = []
if rec_dWPP:
self.dWPP_time_series = []
if rec_dWIP:
self.dWIP_time_series = []
if rec_dBPP:
self.dBPP_time_series = []
if rec_dBPI:
self.dBPI_time_series = []
if rec_dBII:
self.dBII_time_series = []
if rec_uP:
self.uP_time_series = []
if rec_uP_breve:
self.uP_breve_time_series = []
if rec_rP_breve:
self.rP_breve_time_series = []
if rec_rP_breve_HI:
self.rP_breve_HI_time_series = []
if rec_uI:
self.uI_time_series = []
if rec_uI_breve:
self.uI_breve_time_series = []
if rec_rI_breve:
self.rI_breve_time_series = []
if rec_vbas:
self.vbas_time_series = []
if rec_vapi:
self.vapi_time_series = []
if rec_vapi_noise:
self.vapi_noise_time_series = []
if rec_noise:
self.noise_time_series = []
if rec_epsilon:
self.epsilon_time_series = []
if rec_epsilon_LO:
self.epsilon_LO_time_series = []
if rec_lat_mismatch:
self.lat_mismatch_time_series = []
self.rec_per_steps = rec_per_steps
self.rec_counter = 0
def record_step(self, target=None, MSE_only=False, testing=False, validation=False):
if hasattr(self, 'MSE_time_series') and target is not None:
if hasattr(self, "rec_rate_MSE"):
# if MSE should be rate, and target is rate already
if self.rec_rate_MSE and self.rate_target:
mse = MSE(self.rP_breve[-1], target)
# if MSE should be rate, but target is voltage:
# convert target to rate before recording
elif self.rec_rate_MSE and not self.rate_target:
mse = MSE(self.rP_breve[-1], self.activation[-1](target[-1]))
# if MSE should be voltage, and rate is voltage
elif not self.rec_rate_MSE and not self.rate_target:
mse = MSE(self.uP_breve[-1], target)
else:
raise ValueError("MSE cannot be based on voltage if target is rate")
elif self.rate_target:
mse = MSE(self.rP_breve[-1], target)
else:
mse = MSE(self.uP_breve[-1], target)
if testing:
self.MSE_test_time_series.append(mse.copy())
elif validation:
self.MSE_val_time_series.append(mse.copy())
else:
self.MSE_time_series.append(mse.copy())
if not MSE_only:
if hasattr(self, 'error_time_series') and target is not None:
if self.rate_target:
self.error_time_series.append(
self.rP_breve[-1] - target
)
else:
self.error_time_series.append(
self.uP_breve[-1] - target
)
if hasattr(self, 'input_time_series'):
self.input_time_series.append(copy.deepcopy(self.r0))
if hasattr(self, 'target_time_series') and target is not None:
self.target_time_series.append(target.copy())
if hasattr(self, 'WPP_time_series'):
self.WPP_time_series.append(copy.deepcopy(self.WPP))
if hasattr(self, 'WIP_time_series'):
self.WIP_time_series.append(copy.deepcopy(self.WIP))
if hasattr(self, 'BPP_time_series'):
self.BPP_time_series.append(copy.deepcopy(self.BPP))
if hasattr(self, 'BII_time_series'):
self.BII_time_series.append(copy.deepcopy(self.BII))
if hasattr(self, 'BPI_time_series'):
self.BPI_time_series.append(copy.deepcopy(self.BPI))
if hasattr(self, 'dWPP_time_series') and hasattr(self, 'dWPP'):
self.dWPP_time_series.append(copy.deepcopy(self.dWPP))
if hasattr(self, 'dWIP_time_series')and hasattr(self, 'dWIP'):
self.dWIP_time_series.append(copy.deepcopy(self.dWIP))
if hasattr(self, 'dBPP_time_series')and hasattr(self, 'dBPP'):
self.dBPP_time_series.append(copy.deepcopy(self.dBPP))
if hasattr(self, 'dBPI_time_series')and hasattr(self, 'dBPI'):
self.dBPI_time_series.append(copy.deepcopy(self.dBPI))
if hasattr(self, 'dBII_time_series')and hasattr(self, 'dBII'):
self.dBII_time_series.append(copy.deepcopy(self.dBII))
if hasattr(self, 'uP_time_series'):
self.uP_time_series.append(copy.deepcopy(self.uP))
if hasattr(self, 'uP_breve_time_series'):
self.uP_breve_time_series.append(copy.deepcopy(self.uP_breve))
if hasattr(self, 'rP_breve_time_series'):
self.rP_breve_time_series.append(copy.deepcopy(self.rP_breve))
if hasattr(self, 'rP_breve_HI_time_series'):
self.rP_breve_HI_time_series.append(copy.deepcopy(self.rP_breve_HI))
if hasattr(self, 'uI_time_series'):
self.uI_time_series.append(copy.deepcopy(self.uI))
if hasattr(self, 'uI_breve_time_series'):
self.uI_breve_time_series.append(copy.deepcopy(self.uI_breve))
if hasattr(self, 'rI_breve_time_series'):
self.rI_breve_time_series.append(copy.deepcopy(self.rI_breve))
if hasattr(self, 'vbas_time_series'):
self.vbas_time_series.append(copy.deepcopy(self.vbas))
if hasattr(self, 'vapi_time_series'):
self.vapi_time_series.append(copy.deepcopy(self.vapi))
if hasattr(self, 'vapi_noise_time_series'):
self.vapi_noise_time_series.append(copy.deepcopy(self.vapi_noise))
if hasattr(self, 'noise_time_series'):
self.noise_time_series.append(copy.deepcopy(self.noise))
if hasattr(self, 'epsilon_time_series'):
self.epsilon_time_series.append(copy.deepcopy(self.epsilon))
if hasattr(self, 'epsilon_LO_time_series'):
self.epsilon_LO_time_series.append(copy.deepcopy(self.epsilon_LO))
if hasattr(self, 'lat_mismatch_time_series'):
vbashat = [self.gbas / (self.gl + self.gbas + self.gapi) * vbas for vbas in self.vbas]
vbashat[-1] = (self.gl + self.gbas + self.gapi) / (self.gl + self.gbas) * vbashat[-1]
lat_mismatch = [uI_breve - vbashat for uI_breve, vbashat in zip(self.uI_breve, vbashat[1:])]
self.lat_mismatch_time_series.append(copy.deepcopy(lat_mismatch))
def calc_taueff(self):
# calculate tau_eff for pyramidals and interneuron
# taueffP is one value per layer
taueffP = []
for i in self.uP:
taueffP.append(1 / (self.gl + self.gbas + self.gapi))
taueffP[-1] = 1 / (self.gl + self.gbas + self.gntgt)
# tau_eff for output layer in absence of target
taueffP_notgt = [1 / (self.gl + self.gbas)]
taueffI = []
for i in self.uI:
taueffI.append(1 / (self.gl + self.gden + self.gnI))
return taueffP, taueffP_notgt, taueffI
def get_conductances(self):
return self.gl, self.gden, self.gbas, self.gapi, self.gnI, self.gntgt
def get_weights(self):
return self.WPP, self.WIP, self.BPP, self.BPI
def set_weights(self, model=None, WPP=None, WIP=None, BPP=None,
BPI=None):
# if another model is given, copy its weights
if hasattr(model, '__dict__'):
WPP, WIP, BPP, BPI = model.get_weights()
logging.info(f"Copying weights from model {model}")
if WPP is not None: self.WPP = deepcopy_array(WPP)
if WIP is not None: self.WIP = deepcopy_array(WIP)
if BPP is not None: self.BPP = deepcopy_array(BPP)
if BPI is not None: self.BPI = deepcopy_array(BPI)
def set_self_predicting_state(self):
# set WIP and BPI to values corresponding to self-predicting state
for i in range(len(self.BPP)):
self.BPI[i] = - self.BPP[i].copy()
for i in range(len(self.WIP)-1):
self.WIP[i] = self.gbas * (self.gl + self.gden) / (
self.gden * (self.gl + self.gbas +
self.gapi)) * self.WPP[i+1].copy()
if len(self.layers) > 2:
self.WIP[-1] = self.gbas * (self.gl + self.gden) / (self.gden *
(self.gl + self.gbas)) * self.WPP[-1].copy()
def get_voltages(self):
return self.uP, self.uI
def get_old_voltages(self):
return self.uP_old, self.uI_old
def get_breve_voltages(self):
return self.uP_breve, self.uI_breve
def set_voltages(self, model=None, uP=None, uP_old=None, uP_breve=None, uI=None, uI_old=None, uI_breve=None):
# if another model is given, copy its voltages
if hasattr(model, '__dict__'):
uP, uI = model.get_voltages()
uP_old, uI_old = model.get_old_voltages()
uP_breve, uI_breve = model.get_breve_voltages()
logging.info(f"Copying voltages from model {model}")
if uP is not None:
for i in range(len(self.layers)-1):
self.uP[i] = copy.deepcopy(uP[i])
if uP_old is not None:
for i in range(len(self.layers)-1):
self.uP_old[i] = copy.deepcopy(uP_old[i])
if uP_breve is not None:
for i in range(len(self.layers)-1):
self.uP_breve[i] = copy.deepcopy(uP_breve[i])
if uI is not None:
self.uI = copy.deepcopy(uI)
if uI_old is not None:
self.uI_old = copy.deepcopy(uI_old)
if uI_breve is not None:
self.uI_breve = copy.deepcopy(uI_breve)
def calc_vapi(self, rPvec, BPP_mat, rIvec, BPI_mat):
"""
returns apical voltages in pyramidals of a given layer
input: rPvec: vector of rates from pyramidal voltages in
output layer
WPP_mat: matrix connecting pyramidal to pyramidal
rIvec: vector of rates from interneuron voltages in output layer
BPI_mat: matrix connecting interneurons to pyramidals
"""
return BPP_mat @ rPvec + BPI_mat @ rIvec
def calc_vbas(self, rPvec, WPP_mat):
"""
returns basal voltages in pyramidals of a given layer
input: rPvec: vector of rates from pyramidal voltages in
layer below
WPP_mat: matrix connecting pyramidal to pyramidal
"""
return WPP_mat @ rPvec
def calc_vden(self, rPvec, WIP_mat):
"""
returns dendritic voltages in inteneurons
input: rPvec: vector of rates from pyramidal voltages in
layer below
WIP_mat: matrix connecting pyramidal to pyramidal
"""
return WIP_mat @ rPvec
def prospective_voltage(self, uvec, uvec_old, tau, dt=None):
"""
returns an approximation of the lookahead of voltage
vector u at current time
"""
if dt == None:
dt = self.dt
return uvec_old + tau * (uvec - uvec_old) / dt
def evolve_system(self, r0=None, u_tgt=None, learn_weights=True,
learn_lat_weights=True, learn_bw_weights=False,
record=True, testing=False, validation=False, compare_dWPP=False):
"""
evolves the system by one time step:
updates synaptic weights and voltages given input rate r0
"""
# increase timer by dt and round float to nearest dt
self.Time = np.round(self.Time + self.dt,
decimals=int(np.round(-np.log10(self.dt))))
if testing or validation:
self.duP, self.duI = self.evolve_voltages(r0, u_tgt=None) # includes recalc of rP_breve
else:
self.duP, self.duI = self.evolve_voltages(r0, u_tgt) # includes recalc of rP_breve
if learn_weights or learn_bw_weights or learn_lat_weights or compare_dWPP:
self.dWPP, self.dWIP, self.dBPP, self.dBPI = self.evolve_synapses(r0)
# apply evolution
for i in range(len(self.duP)):
self.uP[i] += self.duP[i]
for i in range(len(self.duI)):
self.uI[i] += self.duI[i]
if learn_weights:
for i in range(len(self.dWPP)):
self.WPP[i] += self.dWPP[i]
if learn_lat_weights:
for i in range(len(self.dWIP)):
self.WIP[i] += self.dWIP[i]
for i in range(len(self.dBPI)):
self.BPI[i] += self.dBPI[i]
if learn_bw_weights:
for i in range(len(self.dBPP)):
self.BPP[i] += self.dBPP[i]
# logging.warning("SETTING SPS AFTER dW")
# self.set_self_predicting_state()
# record step
if hasattr(self, 'rec_per_steps') and record:
self.rec_counter += 1
if self.rec_counter % self.rec_per_steps == 0:
self.rec_counter = 0
self.record_step(target=u_tgt)
# during testing and validation, we record MSE for all steps
if (testing or validation) and record:
self.record_step(target=u_tgt, MSE_only=True, testing=testing, validation=validation)
def evolve_voltages(self, r0=None, u_tgt=None):
"""
Evolves the pyramidal and interneuron voltages by one dt
using r0 as input rates
"""
self.duP = [np.zeros(shape=uP.shape) for uP in self.uP]
self.duI = [np.zeros(shape=uI.shape) for uI in self.uI]
# same for dendritic voltages and rates
self.rP_breve_old = deepcopy_array(self.rP_breve)
self.rI_breve_old = deepcopy_array(self.rI_breve)
if self.r0 is not None:
self.r0_old = self.r0.copy()
self.vbas_old = deepcopy_array(self.vbas)
self.vden_old = deepcopy_array(self.vden)
self.vapi_old = deepcopy_array(self.vapi)
# calculate lookahead
self.uP_breve = [self.prospective_voltage(self.uP[i],
self.uP_old[i],
self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [self.prospective_voltage(self.uI[i],
self.uI_old[i],
self.taueffI[i]) for i in range(len(self.uI))]
# calculate rate of lookahead: phi(ubreve)
self.rP_breve = [self.activation[i](self.uP_breve[i])
for i in range(len(self.uP_breve))]
self.rI_breve = [self.activation[i+1](self.uI_breve[i])
for i in range(len(self.uI_breve))]
self.r0 = r0
# before modifying uP and uI, we need to save copies
# for future calculation of u_breve
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
self.vbas, self.vapi, self.vden = self.calc_dendritic_updates(r0, u_tgt)
self.duP, self.duI = self.calc_somatic_updates(u_tgt)
return self.duP, self.duI
def calc_dendritic_updates(self, r0=None, u_tgt=None):
# calculate dendritic voltages from lookahead
if r0 is not None:
self.vbas[0] = self.WPP[0] @ self.r0
for i in range(1, len(self.layers)-1):
self.vbas[i] = self.calc_vbas(self.rP_breve[i-1],
self.WPP[i])
for i in range(len(self.WIP)):
if self.bw_connection_mode == 'skip':
self.vden[0] = self.calc_vden(self.rP_breve[-2],
self.WIP[-1])
elif self.bw_connection_mode == 'layered':
self.vden[i] = self.calc_vden(self.rP_breve[i],
self.WIP[i])
for i in range(len(self.layers)-2):
if self.bw_connection_mode == 'skip':
self.vapi[i] = self.calc_vapi(self.rP_breve[-1],
self.BPP[i],
self.rI_breve[-1],
self.BPI[i])
elif self.bw_connection_mode == 'layered':
self.vapi[i] = self.calc_vapi(self.rP_breve[i+1],
self.BPP[i],
self.rI_breve[i],
self.BPI[i])
return self.vbas, self.vapi, self.vden
def calc_somatic_updates(self, u_tgt=None):
"""
calculates somatic updates from dendritic potentials
"""
# update somatic potentials
for i in range(len(self.uI)):
ueffI = self.taueffI[i] * (self.gden *
self.vden[i] + self.gnI * self.uP_breve[i+1])
delta_uI = (ueffI - self.uI[i]) / self.taueffI[i]
self.duI[i] = self.dt * delta_uI
for i in range(0, len(self.layers)-2):
ueffP = self.taueffP[i] * (self.gbas *
self.vbas[i] + self.gapi * self.vapi[i])
delta_uP = (ueffP - self.uP[i]) / self.taueffP[i]
self.duP[i] = self.dt * delta_uP
if u_tgt is not None:
ueffP = self.taueffP[-1] * (self.gbas *
self.vbas[-1] + self.gntgt * u_tgt[-1])
delta_uP = (ueffP - self.uP[-1]) / self.taueffP[-1]
else:
ueffP = self.taueffP_notgt[-1] * (self.gbas *
self.vbas[-1])
delta_uP = (ueffP - self.uP[-1]) / self.taueffP[-1]
self.duP[-1] = self.dt * delta_uP
return self.duP, self.duI
def evolve_synapses(self, r0, learn_WIP=True, learn_BPI=True):
"""
evolves all synapses by a dt
plasticity of WPP
"""
self.dWPP = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWIP = [np.zeros(shape=WIP.shape) for WIP in self.WIP]
self.dBPP = [np.zeros(shape=BPP.shape) for BPP in self.BPP]
self.dBPI = [np.zeros(shape=BPI.shape) for BPI in self.BPI]
if self.dWPP_use_activation:
if self.dWPP_r_low_pass:
if r0 is not None:
self.r0_LO_old += self.dt / self.tauLO * (self.r0_old - self.r0_LO_old)
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.rP_breve[0] - self.activation[0](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]),
self.r0_LO_old)
for i in range(1, len(self.WPP)-1):
self.r_LO_old[i-1] += self.dt / self.tauLO * (self.rP_breve_old[i-1] - self.r_LO_old[i-1])
# hidden layers
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.rP_breve[i] - self.activation[i](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]),
self.r_LO_old[i-1])
# output layer
self.r_LO_old[-2] += self.dt / self.tauLO * (self.rP_breve_old[-2] - self.r_LO_old[-2])
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.rP_breve[-1] - self.activation[-1](self.gbas / (self.gl + self.gbas) * self.vbas_old[-1]),
self.r_LO_old[-2])
elif self.dWPP_post_low_pass:
if r0 is not None:
self.dWPP_post_LO_old[0] += self.dt / self.tauLO * \
(self.rP_breve[0] - self.activation[0](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]) \
- self.dWPP_post_LO_old[0])
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(self.dWPP_post_LO_old[0], self.r0_old)
for i in range(1, len(self.WPP)-1):
self.dWPP_post_LO_old[i] += self.dt / self.tauLO * \
(self.rP_breve[i] - self.activation[i](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]) \
- self.dWPP_post_LO_old[i])
# hidden layers
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(self.dWPP_post_LO_old[i], self.rP_breve_old[i-1])
# output layer
self.dWPP_post_LO_old[-1] += self.dt / self.tauLO * \
(self.rP_breve[-1] - self.activation[-1](self.gbas / (self.gl + self.gbas) * self.vbas_old[-1]) \
- self.dWPP_post_LO_old[-1])
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(self.dWPP_post_LO_old[-1], self.rP_breve_old[-2])
else:
if len(self.layers) == 2:
# logging.info("updating WPP0")
# input layer
if r0 is not None:
# if the model no hidden layers
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.rP_breve[0] - self.activation[0](self.gbas / (self.gl + self.gbas) * self.vbas_old[0]),
self.r0_old)
else:
# input layer
if r0 is not None:
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.rP_breve[0] - self.activation[0](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]),
self.r0_old)
# hidden layers
for i in range(1, len(self.WPP)-1):
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.rP_breve[i] - self.activation[i](self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]),
self.rP_breve_old[i-1])
# output layer
# logging.info("updating WPP-1")
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.rP_breve[-1] - self.activation[-1](self.gbas / (self.gl + self.gbas) * self.vbas_old[-1]),
self.rP_breve_old[-2])
else:
if self.dWPP_r_low_pass:
if r0 is not None:
self.r0_LO_old += self.dt / self.tauLO * (self.r0_old - self.r0_LO_old)
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uP_breve[0] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]),
self.r0_LO_old)
for i in range(1, len(self.WPP)-1):
self.r_LO_old[i-1] += self.dt / self.tauLO * (self.rP_breve_old[i-1] - self.r_LO_old[i-1])
# hidden layers
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.uP_breve[i] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]),
self.r_LO_old[i-1])
# output layer
self.r_LO_old[-2] += self.dt / self.tauLO * (self.rP_breve_old[-2] - self.r_LO_old[-2])
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.uP_breve[-1] - (self.gbas / (self.gl + self.gbas) * self.vbas_old[-1]),
self.r_LO_old[-2])
elif self.dWPP_post_low_pass:
if r0 is not None:
self.dWPP_post_LO_old[0] += self.dt / self.tauLO * \
(self.uP_breve[0] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]) \
- self.dWPP_post_LO_old[0])
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(self.dWPP_post_LO_old[0], self.r0_old)
for i in range(1, len(self.WPP)-1):
self.dWPP_post_LO_old[i] += self.dt / self.tauLO * \
(self.uP_breve[i] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]) \
- self.dWPP_post_LO_old[i])
# hidden layers
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(self.dWPP_post_LO_old[i], self.r_old[i-1])
# output layer
self.dWPP_post_LO_old[-1] += self.dt / self.tauLO * \
(self.uP_breve[-1] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[-1]) \
- self.dWPP_post_LO_old[-1])
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(self.dWPP_post_LO_old[-1], self.r_old[-2])
else:
if len(self.layers) == 2:
# input layer
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uP_breve[0] - (self.gbas / (self.gl + self.gbas) * self.vbas_old[0]),
self.r0_old)
else:
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uP_breve[0] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]),
self.r0_old)
# hidden layers
for i in range(1, len(self.WPP)-1):
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.uP_breve[i] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[i]),
self.rP_breve_old[i-1])
# output layer
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.uP_breve[-1] - (self.gbas / (self.gl + self.gbas) * self.vbas_old[-1]),
self.rP_breve_old[-2])
"""
plasticity of WIP
"""
if learn_WIP:
if self.bw_connection_mode == 'skip':
self.dWIP[-1] = self.dt * self.eta_IP[-1] * np.outer(
self.rI_breve[-1] - self.activation[-1](self.gden / (self.gl + self.gden) * self.vden_old[-1]),
self.rP_breve_old[-2])
elif self.bw_connection_mode == 'layered':
for i in range(len(self.WIP)):
self.dWIP[i] = self.dt * self.eta_IP[i] * np.outer(
self.rI_breve[i] - self.activation[i+1](self.gden / (self.gl + self.gden) * self.vden_old[i]),
self.rP_breve_old[i])
"""
plasticity of BPI
"""
if learn_BPI:
for i in range(0, len(self.BPI)):
if self.eta_PI[i] != 0.0:
if self.bw_connection_mode == 'skip':
self.dBPI[i] = self.dt * self.eta_PI[i] * np.outer(-self.vapi_old[i], self.rI_breve_old[-1])
elif self.bw_connection_mode == 'layered':
self.dBPI[i] = self.dt * self.eta_PI[i] * np.outer(-self.vapi_old[i], self.rI_breve_old[i])
"""
plasticity of BPP
"""
if self.model == 'FA':
# do nothing
pass
elif self.model == 'BP':
# self.set_weights(BPP = [WPP.T for WPP in self.WPP[1:]])
# add noise to transpose
BPP_noised = [BPP_noise + WPP.T for BPP_noise, WPP in zip(self.BPP_noise, self.WPP[1:])]
self.set_weights(BPP=BPP_noised)
# lateral weight is still perfect copy
self.set_weights(BPI = [-BPP for BPP in self.BPP])
return self.dWPP, self.dWIP, self.dBPP, self.dBPI
class noise_model(base_model):
""" This class inherits all properties from the base model class and adds the function to add noise """
def __init__(self, bw_connection_mode, dWPP_use_activation, dt, dtxi, tauHP, tauLO, Tpres, noise_scale, alpha,
inter_low_pass, pyr_hi_pass, dWPP_low_pass, dWPP_r_low_pass, dWPP_post_low_pass, gate_regularizer,
noise_type, noise_mode,
model, activation, layers,
uP_init, uI_init, WPP_init, WIP_init, BPP_init, BPI_init,
gl, gden, gbas, gapi, gnI, gntgt,
eta_fw, eta_bw, eta_PI, eta_IP, seed=123, **kwargs):
# init base_model with same settings
super().__init__(bw_connection_mode=bw_connection_mode, dWPP_use_activation=dWPP_use_activation, dt=dt, Tpres=Tpres,
model=model, activation=activation, layers=layers,
uP_init=uP_init, uI_init=uI_init,
WPP_init=WPP_init, WIP_init=WIP_init, BPP_init=BPP_init, BPI_init=BPI_init,
gl=gl, gden=gden, gbas=gbas, gapi=gapi, gnI=gnI, gntgt=gntgt,
eta_fw=eta_fw, eta_bw=eta_bw, eta_PI=eta_PI, eta_IP=eta_IP, seed=seed)
self.rng = np.random.RandomState(seed)
self.d_activation = [dict_d_activation[activation.__name__] for activation in self.activation]
# type of noise (OU or white)
self.noise_type = noise_type
# mode of noise injection (order vapi or uP or uP_adative)
self.noise_mode = noise_mode
# for uP_adaptive, we need epsilon: measures angle between BPP, WPP.T
self.epsilon = [np.float64(1.0) for BPP in self.BPP]
if noise_mode == 'uP_adaptive':
self.noise_deg = kwargs.get('noise_deg')
self.tau_eps = kwargs.get('taueps')
if noise_type == 'OU':
self.tauxi = kwargs.get('tauxi')
# low-pass filtered version of epsilon
self.epsilon_LO = deepcopy_array(self.epsilon)
# whether to low-pass filter the interneuron dendritic input
self.inter_low_pass = inter_low_pass
# whether to high-pass filter rPbreve for updates of BPP
self.pyr_hi_pass = pyr_hi_pass
# whether to low-pass filter updates of WPP
self.dWPP_low_pass = dWPP_low_pass
self.dWPP_r_low_pass = dWPP_r_low_pass
self.dWPP_post_low_pass = dWPP_post_low_pass
# whether to gate application of the regularizer
self.gate_regularizer = gate_regularizer
# whether to use phi' B phi' as regularizer
self.varphi_regularizer = kwargs.get('varphi_regularizer', False)
if self.varphi_regularizer:
self.d_rP = [np.diag(np.zeros(shape=uP.shape)) for uP in self.uP]
# noise time scale
self.dtxi = dtxi
# decimals of dt
self.dt_decimals = int(np.round(-np.log10(self.dt)))
# synaptic time constant (sets the low-pass filter of interneuron)
self.tauHP = tauHP
self.tauLO = tauLO
# gaussian noise properties
self.noise_scale = noise_scale
self.noise = [np.zeros(shape=uP.shape) for uP in self.uP]
# self.noise_breve = [np.zeros(shape=uP.shape) for uP in self.uP]
# we need a new variable: vapi after noise has been added
# i.e. vapi = BPP rP + BPI rI (as usual), and vapi_noise = vapi + noise
self.vapi_noise = deepcopy_array(self.vapi)
# init a counter for time steps after which to resample noise
self.noise_counter = 0
self.noise_total_counts = np.round(self.dtxi / self.dt, decimals=self.dt_decimals)
# init a high-pass filtered version of rP_breve
self.rP_breve_HI = deepcopy_array(self.rP_breve)
# init a low-pass filtered version of dWPP
self.dWPP_LO = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWPP_post_LO_old = [np.zeros(shape=rP_breve.shape) for rP_breve in self.rP_breve]
self.r0_LO_old = np.zeros(shape=layers[0])
self.r_LO_old = [np.zeros(shape=rP_breve.shape) for rP_breve in self.rP_breve]
# regularizer for backward weights
self.alpha = alpha
def evolve_system(self, r0=None, u_tgt=None, learn_weights=True, learn_lat_weights=True, learn_bw_weights=True, \
record=True, testing=False, validation=False, compare_dWPP=False):
"""
This overwrites the vanilla system evolution and implements
additional noise
"""
# # update which backwards weights to learn
# if learn_bw_weights and self.Time % self.Tbw == 0:
# # logging.info(f"Current time: {self.Time}s")
# self.active_bw_syn = 0 if self.active_bw_syn == len(self.BPP) - 1 else self.active_bw_syn + 1
# # logging.info(f"Learning backward weights to layer {self.active_bw_syn + 1}")
# self.noise_counter = 0
# calculate voltage evolution, including low pass on interneuron synapses
# see calc_dendritic updates below
if (testing or validation):
self.duP, self.duI = self.evolve_voltages(r0, u_tgt=None, inject_noise=learn_bw_weights) # includes recalc of rP_breve
else:
self.duP, self.duI = self.evolve_voltages(r0, u_tgt, inject_noise=learn_bw_weights) # includes recalc of rP_breve
if learn_weights or learn_lat_weights or compare_dWPP:
self.dWPP, self.dWIP, _, self.dBPI = self.evolve_synapses(r0)
if learn_bw_weights or compare_dWPP:
self.dBPP = self.evolve_bw_synapses()
# apply evolution
for i in range(len(self.duP)):
self.uP[i] += self.duP[i]
for i in range(len(self.duI)):
self.uI [i]+= self.duI[i]
if learn_weights:
if self.dWPP_low_pass:
# calculate lo-passed update of WPP
self.dWPP_LO = self.calc_dWPP_LO()
for i in range(len(self.dWPP_LO)):
self.WPP[i] += self.dWPP_LO[i]
else:
for i in range(len(self.dWPP)):
self.WPP[i] += self.dWPP[i]
if learn_lat_weights:
for i in range(len(self.dWIP)):
self.WIP[i] += self.dWIP[i]
for i in range(len(self.dBPI)):
self.BPI[i] += self.dBPI[i]
if learn_bw_weights:
for i in range(len(self.dBPP)):
self.BPP[i] += self.dBPP[i]
# record step
if hasattr(self, 'rec_per_steps') and record:
self.rec_counter += 1
if self.rec_counter % self.rec_per_steps == 0:
self.rec_counter = 0
self.record_step(target=u_tgt)
# during testing and validation, we record MSE for all steps
if (testing or validation) and record:
self.record_step(target=u_tgt, MSE_only=True, testing=testing, validation=validation)
# increase timer
self.Time = np.round(self.Time + self.dt, decimals=self.dt_decimals)
# calculate hi-passed rP_breve for synapse BPP
self.rP_breve_HI = self.calc_rP_breve_HI()
def evolve_voltages(self, r0=None, u_tgt=None, inject_noise=False):
"""
Overwrites voltage evolution:
Evolves the pyramidal and interneuron voltages by one dt
using r0 as input rates
>> Injects noise into vapi
"""
self.duP = [np.zeros(shape=uP.shape) for uP in self.uP]
self.duI = [np.zeros(shape=uI.shape) for uI in self.uI]
# same for dendritic voltages and rates
self.rP_breve_old = deepcopy_array(self.rP_breve)
self.rI_breve_old = deepcopy_array(self.rI_breve)
if self.r0 is not None:
self.r0_old = self.r0.copy()
self.vbas_old = deepcopy_array(self.vbas)
self.vden_old = deepcopy_array(self.vden)
self.vapi_old = deepcopy_array(self.vapi)
# calculate lookahead
self.uP_breve = [self.prospective_voltage(self.uP[i], self.uP_old[i], self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [self.prospective_voltage(self.uI[i], self.uI_old[i], self.taueffI[i]) for i in range(len(self.uI))]
# calculate rate of lookahead: phi(ubreve)
self.rP_breve = [self.activation[i](self.uP_breve[i]) for i in range(len(self.uP_breve))]
self.rI_breve = [self.activation[i+1](self.uI_breve[i]) for i in range(len(self.uI_breve))]
self.r0 = r0
# before modifying uP and uI, we need to save copies
# for future calculation of u_breve
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
self.vbas, self.vapi, self.vapi_noise, self.vden = self.calc_dendritic_updates(r0, u_tgt)
# inject noise into newly calculated vapi before calculating update du
if inject_noise:
for i in range(len(self.BPP)):
self.generate_noise(layer=i, noise_type=self.noise_type, noise_scale=self.noise_scale)
self.duP, self.duI = self.calc_somatic_updates(u_tgt, inject_noise=inject_noise)
return self.duP, self.duI
def calc_somatic_updates(self, u_tgt=None, inject_noise=False):
"""
this overwrites the somatic update rules
only difference to super: uses vapi_noise instead of vapi
"""
# update somatic potentials
for i in range(len(self.uI)):
ueffI = self.taueffI[i] * (self.gden * self.vden[i] + self.gnI * self.uP_breve[i+1])
delta_uI = (ueffI - self.uI[i]) / self.taueffI[i]
self.duI[i] = self.dt * delta_uI
for i in range(0, len(self.layers)-2):
if inject_noise:
ueffP = self.taueffP[i] * (self.gbas * self.vbas[i] + self.gapi * self.vapi_noise[i])
else:
ueffP = self.taueffP[i] * (self.gbas * self.vbas[i] + self.gapi * self.vapi[i])
delta_uP = (ueffP - self.uP[i]) / self.taueffP[i]
self.duP[i] = self.dt * delta_uP
if u_tgt is not None:
ueffP = self.taueffP[-1] * (self.gbas * self.vbas[-1] + self.gntgt * u_tgt[-1])
delta_uP = (ueffP - self.uP[-1]) / self.taueffP[-1]
else:
ueffP = self.taueffP_notgt[-1] * (self.gbas * self.vbas[-1])
delta_uP = (ueffP - self.uP[-1]) / self.taueffP[-1]
self.duP[-1] = self.dt * delta_uP
return self.duP, self.duI
def generate_noise(self, layer, noise_type, noise_scale):
"""
this function generates noise for a given layer
the noise is added to the apical potential in function calculate_dendritic_updates
there are two noise types:
- hold_white_noise: adds steps of width dtxi and height sampled from normal distribution
- OU: ornstein-uhlenbeck noise, i.e. low-pass filtered white noise updated at every dtxi
"""
# if dtxi timesteps have passed, sample new noise
if np.all(self.noise[layer] == 0) or self.noise_counter % self.noise_total_counts == 0:
if self.noise_mode == 'uP_adaptive':
None # FIX THIS MODE
# # iterate over hidden layers
# for i in range(len(self.layers)-2):
# # first, we reset all noise values
# self.noise[i] = np.zeros(shape=self.uP[i].shape)
# self.vapi_noise[i] = self.vapi[i]
# # 'layer' is the layer which currently has active noise injection/bw learning
# if i == layer:
# # calculate Jacobian alignment factor epsilon
# self.epsilon[layer] = 1/2 * (1 - self.uP_breve[layer] @ self.BPP[layer] @ self.rP_breve[-1] \
# / np.linalg.norm(self.uP_breve[layer]) / np.linalg.norm(self.BPP[layer] @ self.rP_breve[-1]))
# # update low-pass filtered version of epsilon, with filter time-constant Tpres
# self.epsilon_LO[layer] += self.dt/self.tau_eps * (self.epsilon[layer] - self.epsilon_LO[layer])
# # generate noise, rescaled with epsilon
# # if epsilon is below threshold, do not inject noise
# if self.epsilon_LO[layer] > 1/2 * (1 - np.cos(self.noise_deg * np.pi/180)):
# white_noise = noise_scale[layer] * self.epsilon_LO[layer] * np.array([np.random.normal(0, np.abs(x)) for x in self.uP[layer]])
elif self.noise_mode == 'const':
# add noise with magnitude of rescaled uP
white_noise = noise_scale[layer] * self.rng.normal(0, 1, size=self.uP[layer].shape)
elif self.noise_mode == 'uP':
# add noise with magnitude of rescaled uP
white_noise = noise_scale[layer] * np.array([self.rng.normal(0, np.abs(x)) for x in self.uP[layer]])
elif self.noise_mode == 'uP_breve':
# add noise with magnitude of rescaled uP_breve
white_noise = noise_scale[layer] * np.array([self.rng.normal(0, np.abs(x)) for x in self.uP_breve[layer]])
elif self.noise_mode == 'vapi':
# add noise with magnitude of rescaled vapi
white_noise = noise_scale[layer] * np.array([self.rng.normal(0, np.abs(x)) for x in self.vapi[layer]])
# the noise will be added to vapi, depending on the mode
if self.noise_type == 'hold_white_noise':
self.noise[layer] = white_noise
elif self.noise_type == 'OU':
self.noise[layer] += 1 / (self.tauxi) * (- self.dt * self.noise[layer] + np.sqrt(self.dt * self.tauxi) * white_noise)
self.noise_counter = 0
self.noise_counter += 1
def calc_rP_breve_HI(self):
# updates the high-passed instantaneous rate rP_breve_HI which is used to update BPP
# High-pass has the form d v_out = d v_in - dt/tau * v_out
for i in range(len(self.rP_breve)):
self.rP_breve_HI[i] += (self.rP_breve[i] - self.rP_breve_old[i]) - self.dt / self.tauHP * self.rP_breve_HI[i]
# self.rP_breve_HI[-1] += (self.rP_breve[-1] - self.rP_breve_old[-1]) - self.dt / self.tauHP * self.rP_breve_HI[-1]
return self.rP_breve_HI
def calc_dWPP_LO(self):
# updates the low-passed update of WPP
# Low-pass has the form d v_out = dt/tau (v_in - v_out)
for i in range(len(self.dWPP_LO)):
self.dWPP_LO[i] += self.dt / self.tauLO * (self.dWPP[i] - self.dWPP_LO[i])
return self.dWPP_LO
def calc_dendritic_updates(self, r0=None, u_tgt=None):
"""
this overwrites the dendritic updates by adding a low-pass on
the interneuron voltages
and adds vapi_noise somatic voltage after noise injection
"""
# calculate dendritic voltages from lookahead
if r0 is not None:
self.vbas[0] = self.WPP[0] @ self.r0
for i in range(1, len(self.layers)-1):
self.vbas[i] = self.calc_vbas(self.rP_breve[i-1], self.WPP[i])
for i in range(len(self.WIP)):
if self.bw_connection_mode == 'skip':
# add slow response to dendritic compartment of interneurons
if self.inter_low_pass:
self.vden[0] += self.dt / self.tauLO * (self.calc_vden(self.rP_breve[-2], self.WIP[-1]) - self.vden[0])
else:
# else, instant response
self.vden[0] = self.calc_vden(self.rP_breve[-2], self.WIP[-1])
elif self.bw_connection_mode == 'layered':
if self.inter_low_pass:
self.vden[i] += self.dt / self.tauLO * (self.calc_vden(self.rP_breve[i], self.WIP[i]) - self.vden[i])
else:
self.vden[i] = self.calc_vden(self.rP_breve[i], self.WIP[i])
for i in range(0, len(self.layers)-2):
if self.bw_connection_mode == 'skip':
self.vapi[i] = self.calc_vapi(self.rP_breve[-1], self.BPP[i], self.rI_breve[-1], self.BPI[i])
elif self.bw_connection_mode == 'layered':
self.vapi[i] = self.calc_vapi(self.rP_breve[i+1], self.BPP[i], self.rI_breve[i], self.BPI[i])
self.vapi_noise[i] = self.vapi[i] + self.noise[i]
return self.vbas, self.vapi, self.vapi_noise, self.vden
def evolve_bw_synapses(self):
# evolve the synapses of BPPs
self.dBPP = [np.zeros(shape=BPP.shape) for BPP in self.BPP]
if self.varphi_regularizer:
for i in range(len(self.BPP)):
self.d_rP[i] = np.diag(self.d_activation[i](self.uP_breve[i]))
for i in range(len(self.BPP)):
if self.bw_connection_mode == 'skip':
if self.pyr_hi_pass:
r_pre = self.rP_breve_HI[-1]
else:
r_pre = self.rP_breve[-1]
elif self.bw_connection_mode == 'layered':
if self.pyr_hi_pass:
r_pre = self.rP_breve_HI[i+1]
else:
r_pre = self.rP_breve[i+1]
if self.model == "PAL":
self.dBPP[i] = self.dt * self.eta_bw[i] * np.outer(
self.noise[i], r_pre
)
# add regularizer with gate or without
if self.gate_regularizer:
self.dBPP[i] -= self.dt * self.alpha[i] * \
self.eta_bw[i] * self.BPP[i] * d_relu(r_pre)
else:
if self.varphi_regularizer:
self.dBPP[i] -= self.dt * self.alpha[i] * self.eta_bw[i] * self.d_rP[i] @ self.BPP[i] @ self.d_rP[i+1]
else:
# vanilla regularizer
self.dBPP[i] -= self.dt * self.alpha[i] * self.eta_bw[i] * self.BPP[i]
return self.dBPP
class errormc_model(base_model):
""" This class inherits all properties from the base model class and changes it to the error microcircuit """
def __init__(self, fw_connection_mode, bw_connection_mode, dWPP_use_activation,
varphi_transfer, dt, dtxi, tauHP, tauLO, Tpres,
noise_scale, alpha,
pyr_hi_pass, dWPP_low_pass,
noise_type, noise_mode,
model, activation, error_activation, layers,
uP_init, uI_init, WPP_init, WIP_init, BII_init, BPI_init,
gl, gden, gbas, gapi, gntgt,
eta_fw, eta_bw, eta_PI, eta_IP, seed=123, WT_noise=0.0, **kwargs):
# init base_model with same settings
super().__init__(bw_connection_mode=bw_connection_mode, dWPP_use_activation=dWPP_use_activation, dt=dt, Tpres=Tpres,
model=model, activation=activation, layers=layers,
uP_init=uP_init, uI_init=uI_init,
WPP_init=WPP_init, WIP_init=WIP_init, BPP_init=WPP_init, BPI_init=BPI_init,
gl=gl, gden=gden, gbas=gbas, gapi=gapi, gnI=0.0, gntgt=gntgt,
eta_fw=eta_fw, eta_bw=eta_bw, eta_PI=eta_PI, eta_IP=eta_IP, seed=seed, WT_noise=WT_noise)
self.rng = np.random.RandomState(seed)
# forward connection_mode: skip or layered
self.fw_connection_mode = fw_connection_mode
# need to provide learning rate in correct shape
if self.fw_connection_mode == 'skip':
# check that correct shape is passed (list of lists)
assert np.array(eta_fw).shape == (len(layers)-1, len(layers)-1), \
f"Forward connection is 'skip', but eta_fw is not correct shape. " \
f"Provide eta_fw in the shape of (len(layers)-1, len(layers)-1) ({(len(layers)-1, len(layers)-1)}). "\
f"Provided eta_fw: {eta_fw}"
self.eta_fw = np.array(self.eta_fw)
# self.eta_fw = convert_eta_to_matrix_for_skip(eta_fw, WPP_init, layers)
elif self.fw_connection_mode == 'layered':
# check that correct shape is passed (list)
assert ((np.array(eta_fw).shape[0] == len(layers)-1) and (np.array(eta_fw).ndim == 1)), \
f"Forward connection is 'layered', but eta_fw is not correct shape. " \
f"Provide eta_fw in the shape of (len(layers)-1)=={len(layers)-1}. " \
f"Provided eta_fw: {eta_fw}"
# check that same number of vectors has been passed
if len(uP_init) != len(uI_init):
raise ValueError(f"Error and representation init voltages do not have same number of entries")
# assert that skip connections are passed correctly
if len(BII_init) == 1 and len(layers) > 3 and bw_connection_mode == 'layered':
if model != 'BP':
raise ValueError(f"BII_init only has one array, but more are required for {len(layers)} layers \
using bw_connection_mode='layered'. Did you mean bw_connection_mode='skip'?")
else:
self.set_weights(BII=WPP_init[1:])
logging.info("Setting BII = WPP.T")
if self.activation[-1] is not linear:
logging.info("Output layer activation is not linear -- make sure that targets are rates")
# new connectivity
self.BII = deepcopy_array(BII_init)
self.dBII = [np.zeros(shape=BII.shape) for BII in self.BII]
self.d_activation = [dict_d_activation[activation.__name__] for activation in self.activation]
# time constants are different for this model, as errors are given to error units
self.taueffP, self.taueffI_notgt, self.taueffI = self.calc_taueff()
# the output neurons also have a vapi
self.vapi = [np.zeros_like(uP) for uP in self.uP]
self.error_layers = [0] + [len(v) for v in self.uI]
# if a list of error activations has been passed, use it
if isinstance(error_activation, list):
self.error_activation = error_activation
# else, set same activation for all layers
else:
self.error_activation = [error_activation for layer in layers[1:]]
# whether the target provided should be a rate or voltage
self.rate_target = True
# calculate rate of lookahead: phi(ubreve)
logging.info("Defining rI_breve")
self.rI_breve = [self.error_activation[i](self.uI_breve[i])
for i in range(len(self.uI_breve))]
# derivative of rep units
self.d_rP_breve = [np.zeros_like(rP_breve) for rP_breve in self.rP_breve]
# whether to transfer varphi from rep units to error units
self.varphi_transfer = varphi_transfer
if varphi_transfer and dWPP_use_activation:
logging.info("varphi_transfer AND dWPP_use_activation both enabled. Do you want this?")
# type of noise (OU or white)
self.noise_type = noise_type
# mode of noise injection (order vapi or uP or uP_adative)
self.noise_mode = noise_mode
# for uP_adaptive, we need epsilon: measures angle between BPP, WPP.T
if noise_type == 'OU':
self.tauxi = kwargs.get('tauxi')
# disabling of prospectivity for rate calcuation (uP_breve will still be used in dWPP!)
if "rep_lookahead" in kwargs:
self.rep_lookahead = kwargs.get('rep_lookahead')
if "error_lookahead" in kwargs:
self.error_lookahead = kwargs.get('error_lookahead')
# whether to high-pass filter rPbreve for updates of BPP
self.pyr_hi_pass = pyr_hi_pass
# whether to low-pass filter updates of WPP
self.dWPP_low_pass = dWPP_low_pass
# whether to use phi' B phi' as regularizer
self.varphi_regularizer = kwargs.get('varphi_regularizer', False)
if self.varphi_regularizer:
self.d_rP = [np.diag(np.zeros(shape=uP.shape)) for uP in self.uP]
# noise time scale
self.dtxi = dtxi
# decimals of dt
self.dt_decimals = int(np.round(-np.log10(self.dt)))
# synaptic time constant (sets the low-pass filter of interneuron)
self.tauHP = tauHP
self.tauLO = tauLO
# gaussian noise properties
self.noise_scale = noise_scale
self.noise = [np.zeros(shape=uP.shape) for uP in self.uP]
# self.noise_breve = [np.zeros(shape=uP.shape) for uP in self.uP]
# we need a new variable: vapi after noise has been added
# i.e. vapi = BPP rP + BPI rI (as usual), and vapi_noise = vapi + noise
self.vapi_noise = deepcopy_array(self.vapi)
# init a counter for time steps after which to resample noise
self.noise_counter = 0
self.noise_total_counts = np.round(self.dtxi / self.dt, decimals=self.dt_decimals)
# init a high-pass filtered version of rP_breve
self.rP_breve_HI = deepcopy_array(self.rP_breve)
# init a low-pass filtered version of dWPP
self.dWPP_LO = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWPP_post_LO_old = [np.zeros(shape=rP_breve.shape) for rP_breve in self.rP_breve]
self.r0_LO_old = np.zeros(shape=layers[0])
self.r_LO_old = [np.zeros(shape=rP_breve.shape) for rP_breve in self.rP_breve]
# regularizer for backward weights
self.alpha = alpha
# determine whether lateral weights are identity
self.lateral_is_identity = True
if not self.error_layers[1:] == self.layers[1:]:
self.lateral_is_identity = False
for i in range(len(self.WIP)):
if not np.all(self.WIP[i] == np.eye(N=self.WIP[i].shape[0], M=self.WIP[i].shape[1])):
self.lateral_is_identity = False
if not np.all(self.BPI[i] == np.eye(N=self.BPI[i].shape[0], M=self.BPI[i].shape[1])):
self.lateral_is_identity = False
# noise level in setting transpose weights
self.WT_noise = WT_noise
# set transpose weights for BP
if self.model == "BP":
if self.bw_connection_mode == 'layered':
# determine BII = WPP.T, taking unequal dims into account
for i in range(0, len(self.BPI)-1):
WIP_eye = np.eye(N=self.WIP[i+1].shape[0], M=self.WIP[i+1].shape[1])
BPI_eye = np.eye(N=self.BPI[i].shape[0], M=self.BPI[i].shape[1])
# perfect transpose
self.BII[i] = (WIP_eye @ self.WPP[i+1] @ BPI_eye).T
# noise matrix (calculated once)
self.BII_noise = [self.rng.uniform(-self.WT_noise, self.WT_noise, size=BII.shape) for BII in self.BII]
# add noise
self.BII[i] += self.BII_noise[i]
elif self.bw_connection_mode == 'skip':
# determine BII = WPP.T, taking unequal dims into account
WIP_eye = [np.eye(N=WIP.shape[0], M=WIP.shape[1]) for WIP in self.WIP]
BPI_eye = [np.eye(N=BPI.shape[0], M=BPI.shape[1]) for BPI in self.BPI]
WIP_eye = diag_mat(WIP_eye)
BPI_eye = diag_mat(BPI_eye)
# perfect transpose
self.BII[0] = (WIP_eye @ self.WPP[0][self.layers[0]:,self.layers[0]:] @ BPI_eye).T
# noise matrix (calculated once)
self.BII_noise = [self.rng.uniform(-self.WT_noise, self.WT_noise, size=BII.shape) for BII in self.BII]
# add noise
self.BII[0] += self.BII_noise[0]
self.BII[0] = np.triu(self.BII[0])
def calc_taueff(self):
# calculate tau_eff for pyramidals and error units
# taueffP is one value per layer
taueffP = []
for i in self.uP:
taueffP.append(1 / (self.gl + self.gbas + self.gapi))
taueffI = []
for i in self.uI:
taueffI.append(1 / (self.gl + self.gden))
taueffI[-1] = 1 / (self.gl + self.gntgt)
# tau_eff for output layer error units in absence of target
taueffI_notgt = [1 / (self.gl + self.gntgt)]
return taueffP, taueffI_notgt, taueffI
def evolve_system(self, r0=None, u_tgt=None, learn_weights=True, learn_lat_weights=True, learn_bw_weights=False, \
record=True, testing=False, validation=False, compare_dWPP=False):
"""
This overwrites the vanilla system evolution and implements
additional noise
"""
# calculate voltage evolution, including low pass on interneuron synapses
# see calc_dendritic updates below
if testing or validation:
self.duP, self.duI = self.evolve_voltages(r0, u_tgt=None, inject_noise=learn_bw_weights) # includes recalc of rP_breve
else:
self.duP, self.duI = self.evolve_voltages(r0, u_tgt, inject_noise=learn_bw_weights) # includes recalc of rP_breve
if learn_weights or learn_bw_weights or learn_lat_weights or compare_dWPP:
self.dWPP, self.dWIP, _, self.dBPI = self.evolve_synapses(r0)
# apply evolution
for i in range(len(self.duP)):
self.uP[i] += self.duP[i]
for i in range(len(self.duI)):
self.uI [i]+= self.duI[i]
if learn_weights:
if self.dWPP_low_pass:
# calculate lo-passed update of WPP
self.dWPP_LO = self.calc_dWPP_LO()
for i in range(len(self.dWPP_LO)):
self.WPP[i] += self.dWPP_LO[i]
else:
for i in range(len(self.dWPP)):
self.WPP[i] += self.dWPP[i]
if learn_lat_weights:
for i in range(len(self.dWIP)):
self.WIP[i] += self.dWIP[i]
for i in range(len(self.dBPI)):
self.BPI[i] += self.dBPI[i]
# record step
if hasattr(self, 'rec_per_steps') and record:
self.rec_counter += 1
if self.rec_counter % self.rec_per_steps == 0:
self.rec_counter = 0
self.record_step(target=u_tgt)
# during testing and validation, we record MSE for all steps
if (testing or validation) and record:
self.record_step(target=u_tgt, MSE_only=True, testing=testing, validation=validation)
# increase timer
self.Time = np.round(self.Time + self.dt, decimals=self.dt_decimals)
# calculate hi-passed rP_breve for synapse BPP
# self.rP_breve_HI = self.calc_rP_breve_HI()
def evolve_voltages(self, r0=None, u_tgt=None, inject_noise=False):
"""
Overwrites voltage evolution:
Evolves the pyramidal and interneuron voltages by one dt
using r0 as input rates
>> Injects noise into vapi
"""
self.duP = [np.zeros(shape=uP.shape) for uP in self.uP]
self.duI = [np.zeros(shape=uI.shape) for uI in self.uI]
# same for dendritic voltages and rates
self.rP_breve_old = deepcopy_array(self.rP_breve)
self.rI_breve_old = deepcopy_array(self.rI_breve)
if self.r0 is not None:
self.r0_old = self.r0.copy()
self.vbas_old = deepcopy_array(self.vbas)
self.vden_old = deepcopy_array(self.vden)
self.vapi_old = deepcopy_array(self.vapi)
# calculate lookahead
self.uP_breve = [self.prospective_voltage(self.uP[i], self.uP_old[i], self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [self.prospective_voltage(self.uI[i], self.uI_old[i], self.taueffI[i]) for i in range(len(self.uI))]
if u_tgt is None:
self.uI_breve[-1] = self.prospective_voltage(self.uI[-1], self.uI_old[-1], self.taueffI_notgt[-1])
# calculate rate of lookahead: phi(ubreve)
if hasattr(self,"rep_lookahead"):
self.rP_breve = [self.activation[i](self.uP_breve[i]) if self.rep_lookahead[i]
else self.activation[i](self.uP[i])
for i in range(len(self.uP_breve))]
else:
self.rP_breve = [self.activation[i](self.uP_breve[i]) for i in range(len(self.uP_breve))]
self.d_rP_breve = [self.d_activation[i](self.uP_breve[i]) for i in range(len(self.uP_breve))]
if hasattr(self,"error_lookahead"):
self.rI_breve = [self.error_activation[i](self.uI_breve[i]) if self.error_lookahead[i]
else self.error_activation[i](self.uI[i])
for i in range(len(self.uI_breve))]
else:
self.rI_breve = [self.error_activation[i](self.uI_breve[i]) for i in range(len(self.uI_breve))]
self.r0 = r0
# before modifying uP and uI, we need to save copies
# for future calculation of u_breve
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
self.vbas, self.vapi, self.vapi_noise, self.vden = self.calc_dendritic_updates(r0, u_tgt)
# inject noise into newly calculated vapi before calculating update du
if inject_noise:
for i in range(len(self.BII)):
self.generate_noise(layer=i, noise_type=self.noise_type, noise_scale=self.noise_scale)
self.duP, self.duI = self.calc_somatic_updates(u_tgt, inject_noise=inject_noise)
return self.duP, self.duI
def calc_somatic_updates(self, u_tgt=None, inject_noise=False):
"""
this overwrites the somatic update rules
difference to super: uses vapi_noise instead of vapi
and changes circuit to error mc
"""
# update somatic potentials
for i in range(len(self.uI)-1):
ueffI = self.taueffI[i] * (self.gden * self.vden[i])
delta_uI = (ueffI - self.uI[i]) / self.taueffI[i]
self.duI[i] = self.dt * delta_uI
if u_tgt is not None:
ueffI = self.taueffI[-1] * (self.gntgt * self.WIP[-1] @ np.diag(self.d_rP_breve[-1]) @ (u_tgt[-1] - self.rP_breve[-1]))
delta_uI = (ueffI - self.uI[-1]) / self.taueffI[-1]
else:
ueffI = 0.0
delta_uI = (ueffI - self.uI[-1]) / self.taueffI_notgt[-1]
self.duI[-1] = self.dt * delta_uI
for i in range(0, len(self.layers)-1):
if inject_noise:
ueffP = self.taueffP[i] * (self.gbas * self.vbas[i] + self.gapi * self.vapi_noise[i])
else:
ueffP = self.taueffP[i] * (self.gbas * self.vbas[i] + self.gapi * self.vapi[i])
delta_uP = (ueffP - self.uP[i]) / self.taueffP[i]
self.duP[i] = self.dt * delta_uP
return self.duP, self.duI
def calc_vbas(self, WPP_rP_breve_below):
"""
returns basal voltages in pyramidals of a given layer
WPP_rP_breve_below: activity from rep unit below times weight to current layer
"""
return WPP_rP_breve_below
def calc_vapi(self, rIvec, BPI_mat):
"""
returns apical voltages in pyramidals of a given layer for error mc
rIvec: vector of rates from interneuron voltages in output layer
BPI_mat: matrix connecting interneurons to pyramidals
"""
return BPI_mat @ rIvec
def calc_vden(self, deriv_rPvec, WIP_mat, BII_rI_breve_above):
"""
returns dendritic voltages in error units
deriv_rPvec: vector of rates from pyramidal voltages in
layer below
WIP_mat: matrix connecting pyramidal to pyramidal
BII_rI_breve_above: activity from error unit above times weight to current layer
"""
if self.varphi_transfer:
vden = np.diag(WIP_mat @ deriv_rPvec) @ BII_rI_breve_above
else:
vden = BII_rI_breve_above
return vden
def generate_noise(self, layer, noise_type, noise_scale):
"""
this function generates noise for a given layer
the noise is added to the apical potential in function calculate_dendritic_updates
there are two noise types:
- OU: ornstein-uhlenbeck noise, i.e. low-pass filtered white noise updated at every dtxi
"""
# if dtxi timesteps have passed, sample new noise
if np.all(self.noise[layer] == 0) or self.noise_counter % self.noise_total_counts == 0:
if self.noise_type == 'OU':
self.noise[layer] += 1 / (self.tauxi) * (- self.dt * self.noise[layer] + np.sqrt(self.dt * self.tauxi) * white_noise)
self.noise_counter = 0
self.noise_counter += 1
def calc_dWPP_LO(self):
# updates the low-passed update of WPP
# Low-pass has the form d v_out = dt/tau (v_in - v_out)
for i in range(len(self.dWPP_LO)):
self.dWPP_LO[i] += self.dt / self.tauLO * (self.dWPP[i] - self.dWPP_LO[i])
return self.dWPP_LO
def calc_dendritic_updates(self, r0=None, u_tgt=None):
"""
this overwrites the dendritic updates by using the error mc dynamics
and adds vapi_noise somatic voltage after noise injection
"""
if self.fw_connection_mode == 'skip':
# if skip fw connection matrix has been given, first calculate total representation vector
# and multiply with total WPP matrix
r0 = np.zeros(self.layers[0]) if r0 is None else np.array(r0)
# concatenate input with hidden layer rates and multiply with forward weights
total_WPP_rP_breve_below = self.WPP[0] @ np.concatenate([r0] + self.rP_breve)
# the first entries are weight connecting to input, and do not correspond to neurons
total_WPP_rP_breve_below = total_WPP_rP_breve_below[self.layers[0]:]
pos_idx = 0
for i, curr_rP_breve in enumerate(self.rP_breve):
# select correct lines from total forward projected representation units
WPP_rP_breve_below = total_WPP_rP_breve_below[pos_idx:pos_idx+len(curr_rP_breve)]
self.vbas[i] = self.calc_vbas(WPP_rP_breve_below)
pos_idx += len(curr_rP_breve)
elif self.fw_connection_mode == 'layered':
if r0 is not None:
self.vbas[0] = self.WPP[0] @ self.r0
for i in range(1, len(self.layers)-1):
self.vbas[i] = self.calc_vbas(self.WPP[i] @ self.rP_breve[i-1])
# if skip bw connection matrix has been given, first calculate total error vector
# and multiply with total BII matrix
if self.bw_connection_mode == 'skip':
assert np.all(self.BII[0] == np.triu(self.BII[0])), "BII has none-triu entries"
total_BII_rI_breve_above = self.BII[0] @ np.concatenate(self.rI_breve)
pos_idx = 0
for i, curr_rI_breve in enumerate(self.rI_breve[:-1]):
if self.bw_connection_mode == 'skip':
# select correct lines from total backprojected error vector
BII_rI_breve_above = total_BII_rI_breve_above[pos_idx:pos_idx+len(curr_rI_breve)]
self.vden[i] = self.calc_vden(self.d_rP_breve[i],
self.WIP[i], BII_rI_breve_above)
pos_idx += len(curr_rI_breve)
elif self.bw_connection_mode == 'layered':
self.vden[i] = self.calc_vden(self.d_rP_breve[i],
self.WIP[i], self.BII[i] @ self.rI_breve[i+1])
for i in range(0, len(self.layers)-1):
self.vapi[i] = self.calc_vapi(self.rI_breve[i], self.BPI[i])
self.vapi_noise[i] = self.vapi[i] + self.noise[i]
return self.vbas, self.vapi, self.vapi_noise, self.vden
def evolve_synapses(self, r0):
"""
evolves all synapses by a dt
plasticity of WPP
"""
self.dWPP = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWIP = [np.zeros(shape=WIP.shape) for WIP in self.WIP]
self.dBII = [np.zeros(shape=BII.shape) for BII in self.BII]
self.dBPI = [np.zeros(shape=BPI.shape) for BPI in self.BPI]
if self.fw_connection_mode == 'skip':
# if skip fw connection matrix has been given, first calculate total representation vector
# and multiply with total WPP matrix
r0 = np.zeros(self.layers[0]) if self.r0 is None else self.r0
r0_old = np.zeros(self.layers[0]) if self.r0_old is None else self.r0_old
# take all pre-synaptic rates (from previous dt)
rates_old_pre = [r0_old] + self.rP_breve_old
# take all postsynaptic rates (r0 will be skipped, just for ease of computation)
if self.dWPP_use_activation:
rates_post = [r0] + self.rP_breve
else:
rates_post = [r0] + self.uP_breve
# for all post-synaptic rates
post_pos_idx = len(r0)
for i, rate_post in enumerate(rates_post):
if i == 0:
# skip r0
continue
pre_pos_idx = 0
# for all pre-synaptic rates until the layer below the current one
for j, rate_old_pre in enumerate(rates_old_pre[:i]):
if self.dWPP_use_activation:
post_diff = rate_post - self.activation[i-1](self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[i-1])
else:
post_diff = rate_post - (self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[i-1])
tmp_dWPP = self.dt * self.eta_fw[i-1,j] * np.outer(post_diff, rate_old_pre)
# tmp_dWPP = np.outer(post_diff, rate_old_pre)
self.dWPP[0][post_pos_idx:post_pos_idx+len(rate_post),
pre_pos_idx:pre_pos_idx+len(rate_old_pre)] = tmp_dWPP
pre_pos_idx += len(rate_old_pre)
post_pos_idx += len(rate_post)
elif self.fw_connection_mode == 'layered':
if self.dWPP_use_activation:
# input layer
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.rP_breve[0] - \
self.activation[0](self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[0]),
self.r0_old)
# hidden layers
for i in range(1, len(self.WPP)-1):
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.rP_breve[i] - \
self.activation[i](self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[i]),
self.rP_breve_old[i-1])
# output layer
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.rP_breve[-1] - self.activation[-1](self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[-1]),
self.rP_breve_old[-2])
else:
# input layer
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uP_breve[0] - (self.gbas / (self.gl + self.gbas + self.gapi) * self.vbas_old[0]),
self.r0_old)
# hidden layers
for i in range(1, len(self.WPP)-1):
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.uP_breve[i] - (self.gbas / \
(self.gl + self.gbas + self.gapi) * \
self.vbas_old[i]),
self.rP_breve_old[i-1])
# output layer
if len(self.layers) > 2:
# logging.info("updating WPP-1")
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.uP_breve[-1] - (self.gbas / \
(self.gl + self.gbas + self.gapi) \
* self.vbas_old[-1]),
self.rP_breve_old[-2])
"""
plasticity of WIP: weight from rep unit to error unit
"""
for i in range(len(self.WIP)):
if self.eta_IP[i] != 0.0:
self.dWIP[i] = - self.dt * self.eta_IP[i] * (self.WIP[i] - np.eye(N=self.WIP[i].shape[0], M=self.WIP[i].shape[1]))
"""
plasticity of BPI: weight from error unit to rep unit
"""
for i in range(0, len(self.BPI)):
if self.eta_PI[i] != 0.0:
self.dBPI[i] = - self.dt * self.eta_PI[i] * (self.BPI[i] - np.eye(N=self.BPI[i].shape[0], M=self.BPI[i].shape[1]))
"""
plasticity of BII
"""
if self.model == 'FA':
# do nothing
pass
# set transpose weights for BP
if self.model == "BP":
if self.bw_connection_mode == 'layered':
# determine BII = WPP.T, taking unequal dims into account
for i in range(0, len(self.BPI)-1):
WIP_eye = np.eye(N=self.WIP[i+1].shape[0], M=self.WIP[i+1].shape[1])
BPI_eye = np.eye(N=self.BPI[i].shape[0], M=self.BPI[i].shape[1])
# perfect transpose
self.BII[i] = (WIP_eye @ self.WPP[i+1] @ BPI_eye).T
# add noise
self.BII[i] += self.BII_noise[i]
elif self.bw_connection_mode == 'skip':
# determine BII = WPP.T, taking unequal dims into account
WIP_eye = [np.eye(N=WIP.shape[0], M=WIP.shape[1]) for WIP in self.WIP]
BPI_eye = [np.eye(N=BPI.shape[0], M=BPI.shape[1]) for BPI in self.BPI]
WIP_eye = diag_mat(WIP_eye)
BPI_eye = diag_mat(BPI_eye)
# perfect transpose
self.BII[0] = (WIP_eye @ self.WPP[0][self.layers[0]:,self.layers[0]:] @ BPI_eye).T
# add noise
self.BII[0] += self.BII_noise[0]
self.BII[0] = np.triu(self.BII[0])
return self.dWPP, self.dWIP, self.dBII, self.dBPI
def get_weights(self):
return self.WPP, self.WIP, self.BII, self.BPI
def set_weights(self, model=None, WPP=None, WIP=None, BPP=None,
BPI=None, BII=None):
# if another model is given, copy its weights
if hasattr(model, '__dict__'):
WPP, WIP, BII, BPI = model.get_weights()
logging.info(f"Copying weights from model {model}")
if WPP is not None: self.WPP = deepcopy_array(WPP)
if WIP is not None: self.WIP = deepcopy_array(WIP)
if BPP is not None: self.BPP = deepcopy_array(BPP)
if BPI is not None: self.BPI = deepcopy_array(BPI)
if BII is not None: self.BII = deepcopy_array(BII)
class ann_model(base_model):
"""
This class inherits all properties from the base model class and changes it to an ann
Here, uP is equal to uP_breve, and only determined by basal input (gapi = 0)
uI represents the error and is directly backpropagated
"""
def __init__(self, dt, Tpres,
model, activation, layers,
uP_init, uI_init, WPP_init, BPP_init,
gl, gbas,
eta_fw, seed=123, **kwargs):
# init base_model with same settings
super().__init__(bw_connection_mode='layered', dWPP_use_activation=False, dt=dt, Tpres=Tpres,
model=model, activation=activation, layers=layers,
uP_init=uP_init, uI_init=uI_init,
# use WPP_init as placeholder for other weights
WPP_init=WPP_init, WIP_init=WPP_init, BPP_init=BPP_init, BPI_init=BPP_init,
gl=gl, gden=0.0, gbas=gbas, gapi=0.0, gnI=0.0, gntgt=0.0,
eta_fw=eta_fw, eta_bw=None, eta_PI=None, eta_IP=None, seed=seed)
self.rng = np.random.RandomState(seed)
# whether the target provided should be a rate or voltage
self.rate_target = True
if self.activation[-1] is not linear:
logging.info("Output layer activation is not linear -- make sure that targets are rates")
self.d_activation = [dict_d_activation[activation.__name__] for activation in self.activation]
# derivative of rep units
self.d_rP_breve = [np.zeros_like(rP_breve) for rP_breve in self.rP_breve]
def prospective_voltage(self, uvec, uvec_old, tau, dt=None):
"""
for this class, prospective voltage is the same as instantaneous voltage
"""
return uvec
def evolve_voltages(self, r0=None, u_tgt=None):
"""
Evolves the pyramidal and interneuron voltages by one dt
using r0 as input rates
"""
self.duP = [np.zeros(shape=uP.shape) for uP in self.uP]
self.duI = [np.zeros(shape=uI.shape) for uI in self.uI]
# same for dendritic voltages and rates
self.rP_breve_old = deepcopy_array(self.rP_breve)
if self.r0 is not None:
self.r0_old = self.r0.copy()
self.vbas_old = deepcopy_array(self.vbas)
# calculate lookahead
self.uP_breve = [self.prospective_voltage(self.uP[i],
self.uP_old[i],
self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [self.prospective_voltage(self.uI[i],
self.uI_old[i],
self.taueffI[i]) for i in range(len(self.uI))]
# calculate rate of lookahead: phi(ubreve)
self.rP_breve = [self.activation[i](self.uP_breve[i]) for i in range(len(self.uP_breve))]
self.d_rP_breve = [self.d_activation[i](self.uP_breve[i]) for i in range(len(self.uP_breve))]
self.r0 = r0
# before modifying uP and uI, we need to save copies
# for future calculation of u_breve
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
self.duP, self.duI = self.calc_somatic_updates(u_tgt)
return self.duP, self.duI
def calc_somatic_updates(self, u_tgt=None):
"""
calculates somatic updates from dendritic potentials
for ann: this is just bottom-up input for uP
whereas errors are represented by uI
"""
# update errors
for i in range(len(self.uI)-1):
self.duI[i] = np.diag(self.d_rP_breve[i]) @ self.BPP[i] @ self.uI[i+1] - self.uI[i]
if u_tgt is not None:
self.d_rP_breve[-1] = self.d_activation[-1](self.uP_breve[-1])
out_error = np.diag(self.d_rP_breve[-1]) @ (u_tgt[-1] - self.rP_breve[-1])
else:
out_error = 0.0
self.duI[-1] = out_error - self.uI[-1]
# rep units
self.duP[0] = self.gbas / (self.gbas + self.gl) * self.WPP[0] @ self.r0 - self.uP[0]
for i in range(1, len(self.layers)-1):
self.duP[i] = self.gbas / (self.gbas + self.gl) * self.WPP[i] @ self.rP_breve[i-1] - self.uP[i]
return self.duP, self.duI
def evolve_synapses(self, r0):
"""
evolves all synapses by a dt
plasticity of WPP
"""
self.dWPP = [np.zeros(shape=WPP.shape) for WPP in self.WPP]
self.dWIP = [np.zeros(shape=WIP.shape) for WIP in self.WIP]
self.dBPP = [np.zeros(shape=BPP.shape) for BPP in self.BPP]
self.dBPI = [np.zeros(shape=BPI.shape) for BPI in self.BPI]
if len(self.layers) == 2:
# input layer
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uI_breve[0], self.r0_old)
else:
if r0 is not None:
# logging.info("updating WPP0")
self.dWPP[0] = self.dt * self.eta_fw[0] * np.outer(
self.uI_breve[0], self.r0_old)
# hidden layers
for i in range(1, len(self.WPP)-1):
# logging.info(f"updating WPP{i}")
self.dWPP[i] = self.dt * self.eta_fw[i] * np.outer(
self.uI_breve[i], self.rP_breve_old[i-1])
# output layer
# logging.info("updating WPP-1")
if len(self.layers) > 2:
self.dWPP[-1] = self.dt * self.eta_fw[-1] * np.outer(
self.uI_breve[-1], self.rP_breve_old[-2])
"""
plasticity of BPP
"""
if self.model == 'FA':
# do nothing
pass
elif self.model == 'BP':
self.set_weights(BPP = [WPP.T for WPP in self.WPP[1:]])
return self.dWPP, self.dWIP, self.dBPP, self.dBPI
def set_self_predicting_state(self):
pass
class dPC_model(base_model):
""" This class inherits all properties from the base model class and changes it to the dendritic PC microcircuit """
def __init__(self, bw_connection_mode, dWPP_use_activation,
dt, Tpres,
model, activation, layers,
uP_init, uI_init, WPP_init, WIP_init, BPP_init, BPI_init,
gl, gden, gbas, gapi, gntgt,
eta_fw, eta_bw, eta_PI, eta_IP, seed=123, WT_noise=0.0, **kwargs):
# init base_model with same settings
super().__init__(bw_connection_mode=bw_connection_mode, dWPP_use_activation=dWPP_use_activation, dt=dt, Tpres=Tpres,
model=model, activation=activation, layers=layers,
uP_init=uP_init, uI_init=uI_init,
WPP_init=WPP_init, WIP_init=WIP_init, BPP_init=WPP_init, BPI_init=BPI_init,
gl=gl, gden=gden, gbas=gbas, gapi=gapi, gnI=0.0, gntgt=gntgt,
eta_fw=eta_fw, eta_bw=eta_bw, eta_PI=eta_PI, eta_IP=eta_IP, seed=seed, WT_noise=WT_noise)
# overwrite uI_breve to be copy of rP_breve in same layer
self.uI_breve = [rP_breve for rP_breve in self.rP_breve]
# linear activation on interneurons
self.rI_breve = self.uI_breve
# reshape BPI to be square mat of current layer entries
self.BPI = [BPI @ WIP for WIP, BPI in zip(self.WIP, self.BPI)]
def evolve_voltages(self, r0=None, u_tgt=None):
"""
Overwrites voltage evolution:
Only difference is that rI_breve = uI_breve = - rP_breve
"""
self.duP = [np.zeros(shape=uP.shape) for uP in self.uP]
self.duI = [np.zeros(shape=uI.shape) for uI in self.uI]
# same for dendritic voltages and rates
self.rP_breve_old = deepcopy_array(self.rP_breve)
self.rI_breve_old = deepcopy_array(self.rI_breve)
if self.r0 is not None:
self.r0_old = self.r0.copy()
self.vbas_old = deepcopy_array(self.vbas)
self.vden_old = deepcopy_array(self.vden)
self.vapi_old = deepcopy_array(self.vapi)
# calculate lookahead
self.uP_breve = [self.prospective_voltage(self.uP[i],
self.uP_old[i],
self.taueffP[i]) for i in range(len(self.uP))]
self.uI_breve = [rP_breve for rP_breve in self.rP_breve] # modification here
# calculate rate of lookahead: phi(ubreve)
self.rP_breve = [self.activation[i](self.uP_breve[i])
for i in range(len(self.uP_breve))]
self.rI_breve = self.uI_breve # modification here
self.r0 = r0
# before modifying uP and uI, we need to save copies
# for future calculation of u_breve
self.uP_old = deepcopy_array(self.uP)
self.uI_old = deepcopy_array(self.uI)
self.vbas, self.vapi, _ = self.calc_dendritic_updates(r0, u_tgt)
self.duP, _ = self.calc_somatic_updates(u_tgt)
return self.duP, self.duI
def evolve_system(self, r0=None, u_tgt=None, learn_weights=True,
learn_lat_weights=True, learn_bw_weights=False,
record=True, testing=False, validation=False, compare_dWPP=False):
"""
evolves the system by one time step:
updates synaptic weights and voltages given input rate r0
Disables application of dWIP and duI
"""
# increase timer by dt and round float to nearest dt
self.Time = np.round(self.Time + self.dt,
decimals=int(np.round(-np.log10(self.dt))))
if testing or validation:
self.duP, _ = self.evolve_voltages(r0, u_tgt=None) # includes recalc of rP_breve
else:
self.duP, _ = self.evolve_voltages(r0, u_tgt) # includes recalc of rP_breve
if learn_weights or learn_bw_weights or learn_lat_weights or compare_dWPP:
# base_model evolve_synapse sets BPI = - BPP for BP.
# we don't want this here, so we buffer and reset after evolve_synapse
if self.model == 'BP':
BPI_buffer = deepcopy_array(self.BPI)
self.dWPP, _, self.dBPP, self.dBPI = self.evolve_synapses(r0, learn_WIP=False)
if self.model == 'BP':
self.set_weights(BPI = BPI_buffer)
# apply evolution
for i in range(len(self.duP)):
self.uP[i] += self.duP[i]
# for i in range(len(self.duI)):
# self.uI[i] += self.duI[i]
if learn_weights:
for i in range(len(self.dWPP)):
self.WPP[i] += self.dWPP[i]
if learn_lat_weights:
# for i in range(len(self.dWIP)):
# self.WIP[i] += self.dWIP[i]
for i in range(len(self.dBPI)):
self.BPI[i] += self.dBPI[i]
if learn_bw_weights:
for i in range(len(self.dBPP)):
self.BPP[i] += self.dBPP[i]
# logging.warning("SETTING SPS AFTER dW")
# self.set_self_predicting_state()
# record step
if hasattr(self, 'rec_per_steps') and record:
self.rec_counter += 1
if self.rec_counter % self.rec_per_steps == 0:
self.rec_counter = 0
self.record_step(target=u_tgt)
# during testing and validation, we record MSE for all steps
if (testing or validation) and record:
self.record_step(target=u_tgt, MSE_only=True, testing=testing, validation=validation)
def set_self_predicting_state(self):
for i in range(len(self.WIP)-1):
self.BPI[i] = - self.gbas/(self.gl + self.gbas + self.gapi) * self.BPP[i].copy() @ self.WPP[i+1].copy()
if len(self.layers) > 2:
self.BPI[-1] = - self.gbas/(self.gl + self.gbas) * self.BPP[-1].copy() @ self.WPP[-1].copy()

microcircuit.py at commit 7baef7d, no license · at the source

Overview

Authors: Kevin Max^1,2, Ismael Jaras², Arno Granier², Katharina A. Wilmes^2,3, Mihai A. Petrovici²

ORCID iDs: Kevin Max, Mihai A. Petrovici

Institutions: University of Bern (Switzerland); Okinawa Institute of Science and Technology Graduate University (Japan); University of Zurich (Switzerland); ETH Zurich (Switzerland); Institute of Neuroinformatics (Switzerland)

Journal: PLoS computational biology, volume 22, issue 4, article e1014164

Dates: received 8 September 2025; accepted 25 March 2026; published online 17 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1371/journal.pcbi.1014164 · PMID 41996333 · PMCID PMC13089762 · OpenAlex W4412970320

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.07.11.664263

Status: code verified

Categories: computational modeling (no new data) (modality), human (organism)

MeSH: Brain*, Models, Neurological*, Nerve Net*, Neurons*, Action Potentials, Algorithms, Animals, Computational Biology, Computer Simulation, Humans, Learning, Pyramidal Cells, Visual Cortex (* major topic)

Journal subjects: Biology and Life Sciences, Cell Biology, Cellular Types, Animal Cells, Neurons, Neuroscience, Cellular Neuroscience, Neuronal Dendrites, Engineering and Technology, Electrical Engineering, Electrical Circuits, Microcircuits, Computational Biology, Computational Neuroscience, Coding Mechanisms, Cognitive Science, Cognitive Psychology, Learning, Psychology, Social Sciences, Learning and Memory, Computer and Information Sciences, Neural Networks, Anatomy, Brain, Visual Cortex, Medicine and Health Sciences, People and Places, Population Groupings, Professions, Teachers

Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (Swiss National Science Foundation) (225643); Volkswagen Foundation (VolkswagenStiftung); EC | HORIZON EUROPE Framework Programme (Horizon Europe) (101147319); Manfred Stärk-Stiftung; Swiss Federal Commission for Scholarships for Foreign Students (ESKAS)

Citations: cited by 1 paper (Europe PMC); 124 references in the paper

Abstract

Neural responses to mismatches between expected and actual stimuli have been widely reported across different species. How does the brain use such error signals for learning? While global error signals can be useful, their ability to learn complex computation at the scale observed in the brain is lacking. In comparison, more local, neuron-specific error signals enable superior performance, but their computation and propagation remain unclear. Motivated by the breakthrough of deep learning, this has inspired the ‘backpropagation and the brain’ hypothesis, i.e., that the brain implements a form of the error backpropagation algorithm. In this work, we introduce a biologically motivated, multi-area cortical microcircuit model, implementing error backpropagation under consideration of recent physiological evidence. We model populations of cortical pyramidal cells acting as representation and error neurons, with bio-plausible local and inter-area connectivity, guided by experimental observations of connectivity of the primate visual cortex. In our model, all information transfer is biologically motivated, inference and learning occur without phases, and network dynamics demonstrably approximate those of error backpropagation. We show the capabilities of our model on a wide range of benchmarks, and compare to other models, such as dendritic hierarchical predictive coding. In particular, our model addresses shortcomings of other theories in terms of scalability to many cortical areas. Finally, we make concrete predictions, which differentiate it from other theories, and which can be tested experimentally.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

Zenodo 16909564

License: CC-BY-4.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: “Data Availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (14 files), Matplotlib (8 files), SciPy (5 files), PyTorch (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)

29 September 2026: the link answers (HTTP 200)

22 files

numpy_model/DMS_task.py, Python, 480 lines
numpy_model/cart-pole plot.ipynb, Jupyter, 364 lines
numpy_model/cartpole/cartpole.py, Python, 224 lines
numpy_model/cartpole/controller.py, Python, 82 lines
numpy_model/generate figs from runs.ipynb, Jupyter, 1,612 lines
numpy_model/run_Fig3_multilayer_comparison.sh, Shell, 9 lines
numpy_model/run_Fig4_populations.sh, Shell, 5 lines
numpy_model/run_Fig5_noise.sh, Shell, 11 lines
numpy_model/run_FigA1_multilayer_comparison_hierarchical.sh, Shell, 9 lines
numpy_model/run_FigA2_multilayer_comparison_ideal_lat_inh.sh, Shell, 9 lines
numpy_model/run_FigC1_nonlinear_dPC.sh, Shell, 29 lines
numpy_model/runner.py, Python, 257 lines
numpy_model/slurm_script.sh, Shell, 20 lines
numpy_model/src/compare.py, Python, 528 lines
numpy_model/src/init_MC.py, Python, 613 lines
numpy_model/src/init_signals.py, Python, 108 lines
numpy_model/src/microcircuit.py, Python, 2,444 lines
numpy_model/src/plot_exp.py, Python, 720 lines
numpy_model/src/run_exp.py, Python, 258 lines
numpy_model/src/save_exp.py, Python, 17 lines
numpy_model/varphi' transfer.ipynb, Jupyter, 236 lines
README.md, Text, 16 lines

At the source: doi.org/10.5281/zenodo.16909564

kma-code/error-neuron-microcircuits

License: none: the authors keep all their rights

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Commit: 7baef7db7acece086d362ff0b71fe29454e693c9, 2 February 2026

Languages: Python (11), Shell (7), Jupyter (3)

Size: 325 files, 21 scripts

Software Heritage: not archived

Found in: the Zenodo archive record

Holds: README, environment (numpy_model/requirements.txt), 3 notebooks

Not found: license file, CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (14 files), Matplotlib (8 files), SciPy (5 files), PyTorch (1 file)

Availability: 1 check, the latest on 29 September 2026: the link answers

29 September 2026: the link answers

22 files

numpy_model/DMS_task.py, Python, 486 lines
numpy_model/cart-pole plot.ipynb, Jupyter, 364 lines, 1 match
numpy_model/cartpole/cartpole.py, Python, 224 lines
numpy_model/cartpole/controller.py, Python, 82 lines
numpy_model/generate figs from runs.ipynb, Jupyter, 1,612 lines
numpy_model/run_Fig3_multilayer_comparison.sh, Shell, 9 lines
numpy_model/run_Fig4_populations.sh, Shell, 5 lines
numpy_model/run_Fig5_noise.sh, Shell, 11 lines
numpy_model/run_FigA1_multilayer_comparison_hierarchical.sh, Shell, 9 lines
numpy_model/run_FigA2_multilayer_comparison_ideal_lat_inh.sh, Shell, 9 lines
numpy_model/run_FigC1_nonlinear_dPC.sh, Shell, 29 lines
numpy_model/runner.py, Python, 257 lines
numpy_model/slurm_script.sh, Shell, 20 lines
numpy_model/src/compare.py, Python, 528 lines
numpy_model/src/init_MC.py, Python, 613 lines
numpy_model/src/init_signals.py, Python, 108 lines
numpy_model/src/microcircuit.py, Python, 2,444 lines, 1 match
numpy_model/src/plot_exp.py, Python, 720 lines
numpy_model/src/run_exp.py, Python, 258 lines
numpy_model/src/save_exp.py, Python, 17 lines
numpy_model/varphi' transfer.ipynb, Jupyter, 236 lines
README.md, Text, 16 lines

At the source: github.com/kma-code/Error-Neuron-Microcircuits/tree/v1

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
42 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

Code for the network simulations, generating the datasets and producing all figures is accessible under DOI https://doi.org/10.5281/zenodo.16909564.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 13 MeSH terms, 5 funders, 98 references.

Cite

This paper

Max, K., Jaras, I., Granier, A., Wilmes, K. A., & Petrovici, M. A. (2026). 'Backpropagation and the brain' realized in cortical error neuron microcircuits. PLoS computational biology, 22(4), e1014164. https://doi.org/10.1371/journal.pcbi.1014164

BibTeX

@article{max2026backpropagation, author = {Max, Kevin and Jaras, Ismael and Granier, Arno and Wilmes, Katharina A. and Petrovici, Mihai A.}, title = {{'Backpropagation and the brain' realized in cortical error neuron microcircuits}}, journal = {PLoS computational biology}, year = {2026}, month = apr, volume = {22}, number = {4}, pages = {e1014164}, publisher = {PLOS}, issn = {1553-734X}, doi = {10.1371/journal.pcbi.1014164}, url = {https://doi.org/10.1371/journal.pcbi.1014164}, pmid = {41996333}, pmcid = {PMC13089762} }

RIS

TY - JOUR AU - Max, Kevin AU - Jaras, Ismael AU - Granier, Arno AU - Wilmes, Katharina A. AU - Petrovici, Mihai A. TI - 'Backpropagation and the brain' realized in cortical error neuron microcircuits T2 - PLoS computational biology J2 - PLoS Comput Biol PY - 2026 DA - 2026/04/17 VL - 22 IS - 4 SP - e1014164 SN - 1553-734X PB - PLOS DO - 10.1371/journal.pcbi.1014164 UR - https://doi.org/10.1371/journal.pcbi.1014164 LA - en ER -

CSL-JSON

{ "id": "10.1371/journal.pcbi.1014164", "type": "article-journal", "title": "'Backpropagation and the brain' realized in cortical error neuron microcircuits", "container-title": "PLoS computational biology", "author": [ { "family": "Max", "given": "Kevin" }, { "family": "Jaras", "given": "Ismael" }, { "family": "Granier", "given": "Arno" }, { "family": "Wilmes", "given": "Katharina A." }, { "family": "Petrovici", "given": "Mihai A." } ], "container-title-short": "PLoS Comput Biol", "volume": "22", "issue": "4", "page": "e1014164", "DOI": "10.1371/journal.pcbi.1014164", "PMID": "41996333", "PMCID": "PMC13089762", "ISSN": "1553-734X", "publisher": "PLOS", "URL": "https://doi.org/10.1371/journal.pcbi.1014164", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 17 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.