"""
This module contains functions to create network topologies and initialize weights. Currently, only the
random connectivity topology is implemented.
"""
from os.path import join as osjoin
from os.path import exists
import numpy as np
from scipy import signal
from scipy.sparse import coo_array
from scipy.stats import gmean
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV, RidgeClassifier
from sklearn.svm import SVC, LinearSVC
from sklearn.multiclass import OneVsRestClassifier, OneVsOneClassifier
from sklearn.manifold import TSNE, Isomap, LocallyLinearEmbedding
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler, OrdinalEncoder
from tqdm import tqdm
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['savefig.bbox'] = 'tight'
mpl.colormaps.unregister("hcmap")
hetero_colors = [
"#ae3832",
"#9fbfe2",
"#5d95d0",
"#336090",
]
hcmap = mpl.colors.ListedColormap(hetero_colors, name="hcmap")
mpl.colormaps.register(cmap=hcmap)
def_rng = np.random.default_rng()
## Viz UTILS
[docs]
def polish(fig, axs, lw=1, tight=True):
if type(axs) == type({}):
axs = list(axs.values())
fig.align_xlabels(axs)
fig.align_ylabels(axs)
# fig.align_titles(axs)
axs = np.array(axs)
for ax in axs.ravel():
simpleaxis(ax, lw)
if tight:
plt.tight_layout()
[docs]
def simpleaxis(ax, lw=2):
for axis in ['top','bottom','left','right']:
ax.spines[axis].set_linewidth(lw)
ax.xaxis.set_tick_params(width=lw)
ax.yaxis.set_tick_params(width=lw)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
[docs]
def get_cax(ax, width=3., dist=.05):
"""
Returns a human-friendly colorbar axis from a given ax, to be used with
matplotlib's colorbar method.
"""
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
cax = inset_axes(ax, width=f"{width}%", height="100%",
loc='lower left', bbox_to_anchor=(1 + dist, 0., 1, 1),
bbox_transform= ax.transAxes, borderpad=0)
return cax
[docs]
def create_topology(N, N_syn, allow_autapse=False, topology='random', rng=None):
"""
Create sparse connectivity matrix following specified topology rules.
Parameters
----------
N : int
Number of neurons
N_syn : int
Number of synapses (determines sparsity)
allow_autapse : bool, optional
Whether to allow self-connections (default False)
topology : str, optional
Type of topology to create (default 'random')
Currently supported: 'random'
rng : numpy.random.Generator, optional
Random number generator instance
Returns
-------
scipy.sparse.csr_matrix
Sparse adjacency matrix of shape (N, N)
"""
if rng is None:
rng = def_rng
if topology == 'random':
idxs = np.random.choice(N**2, N_syn, replace=False)
idxs = np.sort(idxs) # it is important to sort out the indices for sparse representation
rows = (idxs//N).astype(int)
cols = idxs % N
if not allow_autapse:
# For each post neuron which is connected to itself, we pick a
# new presynapse from 0 to N-1 which is different from all
# presynapses of that particular neuron
autapses = cols[rows == cols] # index in sparse matrix
new_pres = []
for autapse in autapses:
pres_of_autapse = cols[rows == autapse] # all pres of a post
new_pre = autapse
while new_pre in pres_of_autapse:
new_pre = np.random.randint(N)
# print(autapse, pres_of_autapse, new_pre)
new_pres.append( new_pre )
cols[np.where(rows == cols)] = new_pres
A_adj = coo_array((np.ones_like(rows), (rows, cols)),
shape=(N, N)).tocsr()
assert(N_syn == A_adj.nnz)
else:
raise NotImplementedError(f'Currently we only support random topology. {topology} was requested.')
return A_adj
[docs]
def intialize_weights(N_input, N, N_output,
p, f, mu_e, sigma_0,
allow_autapse=False,
return_adj=True, rng=None):
"""
Creates three weight matrices:
- the feedforward weights ``w_ff`` mapping inputs to neurons,
- the recurrent weights ``w`` between neurons, and
- an empty feedback weight array ``w_fb`` (a placeholder for future use).
The recurrent weights are initialized to ensure excitatory/inhibitory balance, whereas
the feedforward weights are initialized from a centered normal distribution.
Parameters
----------
N_input : int
Number of input dimensions
N : int
Number of neurons
N_output : int
Number of output dimensions
p : float
Connection probability
f : float
Fraction of excitatory neurons
mu_e : float
Mean excitatory weight. The mean inhibitory weight is computed based on f to ensure E/I balance.
sigma_0 : float
Recurrent weight standard deviation (same for both excitatory and inhibitory)
return_adj : bool, optional
Whether to return adjacency matrices. Default is True.
rng : numpy.random.Generator, optional
Random number generator instance
Returns
-------
w_ff : ndarray
Input weights of shape (N, N_input)
w : ndarray
Recurrent weights of shape (N, N)
w_fb : ndarray
Feedback weights of shape (N, N_output)
.. note::
- The recurrent weights are scaled by :math:`1/\\sqrt{Np}` to ensure the expected recurrent input
variance, independent of the network size and connection probability, is equal to ``sigma_0``.
- The feedforward weights are scaled by :math:`1/\\sqrt{dim_{input}}` to ensure that feedforward
input variance on each neuron is equal to ``1`` independent of the input dimension.
- The feedforward and recurrent gains (:math:`J_u` and :math:`J`, respectively), will be applied
later during integrating the dynamical system.
"""
if rng is None:
rng = def_rng
# treat special cases for mu_i to ensure E/I balance
if f==0: # no excitatory
mu_i = -mu_e
elif f==1: # no inhibitory
mu_i = 0 # doesn't matter as there's no inhibitory pop anyways
else:
mu_i = -f/(1-f) * mu_e
std_e = std_i = sigma_0
N_syn = int(round(N**2 * p))
Ne = int(N*f) # number of neurons
# Ni = N - Ne # number of neurons
## FEED-FORWARD WEIGHTS ##
w_in = rng.normal(0,1, size=(N, N_input))
w_in /= np.sqrt(N_input)
# even though we know Ne and Ni, the randomness of A_adj
# precludes us from knowing which first elements of A_adj
# are in the inhibitory/excitatory population. Thus, we
# look at their pre and post indices. (row is post)
A_adj = create_topology(N=N, N_syn=N_syn, allow_autapse=allow_autapse)
posts, pres = A_adj.nonzero()
w_rec = np.zeros(N_syn)
pres_e = pres<Ne
w_rec[pres_e] = rng.normal(mu_e, std_e, size=sum(pres_e))
w_rec[~pres_e] = rng.normal(mu_i, std_i, size=sum(~pres_e))
w_rec /= np.sqrt(N*p) # to make variance independent of N and p
## FEEDBACK WEIGHTS ##
w_out = rng.normal(0,1, size=(N_output, N))
w_out /= np.sqrt(N_output)
if return_adj:
return w_in, w_rec, w_out, A_adj
else:
return w_in, w_rec, w_out
# def extract_base_tau_unstream(
# us, dt, filt=True, savepath='', savename=None,
# pad=False, trial_onset=None, trial_offset=None,
# ):
# # print('WARNING: DOUBLE CHECK THAT YOU ARE PASSING THE TRAIN SET!')
# # temporal axis must be 0 for conveniance
# if us.ndim == 1:
# _us = us[:, np.newaxis]
# elif us.ndim == 2:
# _us = us
# # heuristic to discover the shape and size of
# # features and time
# dim0, dim1 = _us.shape
# if dim0 < dim1:
# _us = _us.T
# if pad == True:
# assert trial_onset is not None and trial_offset is not None
# assert len(trial_offset) == len(trial_onset)
# ntrial = len(trial_offset)
# max_length = (trial_offset - trial_onset).max()
# _us_pad = np.zeros((max_length, ntrial, _us.shape[1])) # T x B x features
# for n, (start, end) in enumerate(zip(trial_onset, trial_offset)):
# _us_pad[:end-start, n, :] = _us[start:end, :]
# _us = _us_pad
# # print(_us.shape)
# _us -= _us.mean(axis=0) # centering the input
# freqs, Pxx = signal.periodogram(_us, 1/dt, axis=0)
# if filt:
# b, a = signal.ellip(4, 0.01, 120, 0.125) # Filter to be applied.
# # print(Pxx.shape)
# Pxx = signal.filtfilt(b, a, Pxx, method="gust", axis=0)
# # print(Pxx.shape)
# # print(Pxx[0].ravel())
# if Pxx.ndim == 3:
# # Pxx = Pxx.mean(axis=1) # average over all trials
# Pxx = Pxx.reshape(len(freqs), -1) # flattening
# fmax = freqs[np.argmax(Pxx, axis=0)]
# base_freq = fmax.mean()
# base_tau = 2*np.pi/base_freq
# if savename is not None:
# plt.figure(figsize=(6,3))
# for k in range(min(Pxx.shape[1], 20)):
# plt.loglog(freqs, Pxx[:, k]/max(Pxx[:, k]), lw=.5, alpha=.6)
# plt.xlabel('Frequency [1/t unit]')
# plt.ylabel('Power density')
# plt.savefig(osjoin(savepath, f'spectra_{savename}.png'))
# plt.close('all')
# plt.figure(figsize=(3,2))
# plt.hist(fmax, bins=100);
# plt.axvline(base_freq, c='k', label='base freq')
# plt.xlabel('Freq [1/time unit]')
# plt.legend()
# plt.savefig(osjoin(savepath, f'freqs_{savename}.png'))
# # plt.close('all')
# return base_tau, fmax
[docs]
def get_spaced_hom_taus(lens, dt=1, num=10, ):
"""
We select taus smaller than the minimal length and taus between
the minimal and maximal length.
"""
# if min(lens) < max(lens):
# qs_between = np.linspace(0, 1, num//2)
# taus_between = np.quantile(lens, qs_between).astype(int)
# # taus smaller
# # we drop the first (0, unphysical) and the last (min(lens),
# # included in the other bit).
# taus_small = np.linspace(0, lens.min(), num//2+2, dtype=int)[1:-1]
# # taus between min and max of lens (quantiled)
# taus = np.concatenate([taus_small, taus_between])
# else:
# taus = np.linspace(0, lens.min(), num+1, dtype=int)[1:]
taus = np.linspace(0, lens.max(), num+1, dtype=int)[1:]
return taus * dt
# def extract_base_tau(us, dt, #modality_dim=1,
# filt=True, savepath='', savename=None):
# # print('WARNING: DOUBLE CHECK THAT YOU ARE PASSING THE TRAIN SET!')
# # temporal axis must be 0 for conveniance
# if us.ndim == 1:
# _us = us[:, np.newaxis]
# elif us.ndim == 2:
# _us = us
# # heuristic to discover the shape and size of
# # features and time
# dim0, dim1 = _us.shape
# if dim0 < dim1:
# _us = _us.T
# elif us.ndim == 3: # (B, T, D) tensor
# _us = us.swapaxes(0,1) # B, T, D --> T, B, D
# # print(_us.shape)
# _us -= _us.mean(axis=0) # centering the input
# freqs, Pxx = signal.periodogram(_us, fs=1/dt, axis=0)
# if filt:
# b, a = signal.ellip(4, 0.01, 120, 0.125) # Filter to be applied.
# # print(Pxx.shape)
# Pxx = signal.filtfilt(b, a, Pxx, method="gust", axis=0)
# # print(Pxx.shape)
# fmax = freqs[np.argmax(Pxx, axis=0)]
# # fmax = np.array(fmax)
# base_freq = fmax.mean()
# base_tau = 2*np.pi/base_freq
# # fmax = []
# # for i, u in enumerate(_us):
# # f, Pxx_den = signal.periodogram(u, 1/dt)
# # if filt:
# # Pxx_den = signal.filtfilt(b, a, Pxx_den, method="gust")
# # fmax.append(f[np.argmax(Pxx_den)])
# # if plot_spectrum:
# # plt.loglog(f, Pxx_den, lw=.2, alpha=.6)
# # plt.savefig(osjoin(savepath, f'spectrum_{i}_{savename}.png'))
# # plt.close('all')
# if savename is not None:
# plt.figure(figsize=(6,3))
# for k in range(_us.shape[1]):
# plt.loglog(freqs, Pxx[:, k]/max(Pxx[:, k]), lw=.5, alpha=.6)
# plt.xlabel('Frequency [1/t unit]')
# plt.ylabel('Power density')
# plt.savefig(osjoin(savepath, f'spectra_{savename}.png'))
# plt.close('all')
# plt.figure(figsize=(3,2))
# plt.hist(fmax, bins=100);
# plt.axvline(base_freq, c='k', label='base freq')
# plt.xlabel('Freq (Hz)')
# plt.legend()
# plt.savefig(osjoin(savepath, f'freqs_{savename}.png'))
# # plt.close('all')
# return base_tau, fmax
[docs]
def sample_taus(size, mu, std, distr='lognormal',
clip_below=False, rng=None, tol_prcntg=-1, verbose=False,
savepath='.', savename=None):
"""
The given mu and std are the desired mean and standard deviation and
are different from the parameters passed to rng. It can also handle
the case of std=0, by returning a constant vector of the desired size.
"""
if rng is None:
rng = np.random.default_rng()
if std == 0:
taus = np.ones(size)*mu
else:
if distr == 'lognormal':
std_ = np.sqrt(np.log((std**2) / (mu**2) + 1))
mu_ = np.log(mu) - (std_**2) / 2
# mu_ = np.log(mu**2/np.sqrt(mu**2 + std**2))
# std_ = np.sqrt(np.log(1+(std/mu)**2))
rng_ = rng.lognormal
kwargs = {'sigma':std_, 'mean': mu_}
elif distr == 'normal':
rng_ = rng.normal
kwargs = {'loc':mu, 'scale':std}
elif distr == 'gamma':
k = (mu / std)**2
theta = mu / k
rng_ = rng.gamma
kwargs = {'shape':k, 'scale':theta}
elif distr == 'uniform':
# The mean and variance can be used to find the min and max
# of the distribution.
# mean = 1/2 (low + high)
# var = 1/12 (high - low)**2
low = mu - 3*(std**2)/mu
high = mu + 3*(std**2)/mu
rng_ = rng.uniform
kwargs = {'low':low, 'high':high}
else:
raise NotImplementedError('Only gamma and lognoraml are supported.')
taus = rng_(size=size, **kwargs)
if tol_prcntg>0:
# Empirical mean and variance, due to sampling can be
# different from the prescribed values. To handle this,
# we simply redraw taus. This only matters at small N
err_mu = abs(taus.mean()-mu)/mu
err_std= abs(taus.std()-std)/(std+1e-12) * (std > 0)
correction = 0
while (err_mu > tol_prcntg) or (err_std>tol_prcntg):
if correction % 10000 == 0:
print(f'{correction} is made.')
taus = rng_(size=size, **kwargs)
err_mu = (taus.mean()-mu)/mu
err_std= (taus.std()-std)/(std+1e-12) * (std > 0)
correction+=1
if bool(clip_below):
resample = taus[taus < clip_below]
resample_size = len(resample)
while resample_size>0:
taus[taus < clip_below] = rng_(size=resample_size, **kwargs)
resample = taus[taus < clip_below]
resample_size = len(resample)
if savename is not None:
plt.figure(figsize=(3.5,2.5))
plt.hist(taus, bins=100, color='C0')
plt.axvline(clip_below, c='r', ls='--', label='lower bound')
plt.axvline(mu, c='k', label='mean')
plt.xlabel(r'$\tau$ [t unit]')
plt.legend()
plt.savefig(osjoin(savepath, f'taus_{savename}.png'))
plt.close('all')
return taus
[docs]
def generate_net(N_input, N=50,
p=.1, f=0.8, mu_e=1, sigma_0=1,
allow_autapse=False,
savepath='.', savename='net',
**kwargs):
# A_adj = create_topology(N, int(np.rint(N**2 * p)))
w_in, w_rec, _, A_adj = intialize_weights(
N = N,
p=p, f=f, mu_e=mu_e, sigma_0=sigma_0,
N_input=N_input,
N_output=1, # dummy
return_adj=True,
allow_autapse=allow_autapse,
**kwargs
)
if savename is not None:
fig, axs = plt.subplots(1, 2, figsize=(8,3))
axs[0].spy(A_adj.toarray(), origin='lower')
axs[0].set_xlabel('pre')
axs[0].set_ylabel('post')
axs[1].hist(w_in.ravel(), bins=100, alpha=0.5, density=True, label='w_in');
axs[1].hist(w_rec.ravel(), bins=100, alpha=0.5, density=True, label='w_rec');
axs[1].set_xlabel('weight value')
axs[1].legend()
plt.tight_layout()
plt.savefig(osjoin(savepath, f'{savename}.png'))
plt.close('all')
return A_adj, w_in, w_rec
[docs]
def activation_func(x, af, x0=1., thr=0, **kwargs):
"""
Apply a specified activation function to input data with optional scaling and thresholding.
The valid activation fucntions are:
- 'lin': linear
- 'relu': rectified linear unit
- 'caplin' : same as 'lin' but capped between 0 and 1
- 'caprelu': same as 'relu' but capped between 0 and 1
- 'capquad': quadratic function capped between 0 and 1
- 'capcube': cubic function capped between 0 and 1
- 'capcube': cubic function capped between 0 and 1
- 'capexp': exponential function capped between 0 and 1
- 'tanh': tangent hyperbolic
- `symsig`: an even sigmoidal function centered at t=0 with range between -1 and 1.
- 'possig': sigmoidal function centered at y=0.5 with range between 0 and 1 .
- 'possig_plus': same as 'possig' but only for values larger than threshold and zero otherwise.
- 'bell' : non-normalized Gaussian activation function
- 'camel': a gaussian multiplied by a pure quadratic function, resulting in two identical bumps
at (±1, 1).
- 'switch': a binary activation. Zero for values smaller than thr and 1 otherwise.
Parameters
----------
x : float or ndarray
Input array or value to be transformed
af: str
Name of the activation function to apply (e.g., 'relu', 'tanh', 'possig', etc.)
x0: float
Scaling factor for the input (default: 1.0)
thr: float
Threshold (bias) value to shift the input (default: 0)
kwargs: dict
Additional keyword arguments for specific activation functions
Returns
-------
float or ndarray
The result of applying the chosen activation function to the input
"""
if af=='lin': # linear
return x/x0
elif af=='relu': # relu
# return np.maximum(0, (x-thr)/x0)
return np.clip((x-thr)/x0, 0, None)
elif af=='caprelu': # capped relu
return np.clip((x-thr)/x0, 0, 1)
elif af=='caplin': # capped linear
return np.clip((x-thr)/x0, 0, 1)
elif af=='capquad': # capped linear
return np.clip(((x-thr)/x0)**2, 0, 1)
elif af=='capcube': # capped linear
return np.clip(((x-thr)/x0)**3, 0, 1)
elif af=='capexp': # capped exponential
return np.exp(np.minimum((x-thr)/x0, 0))
elif af=='tanh': # hyperbolic tangent
return np.tanh((x-thr)/x0)
elif af=='possig': # positive sigmoid
"""Clipping at 11 is to avoid overflow"""
return 1./(1+np.exp(-np.clip((x-thr)/x0,-60, 60) ))
elif af=='possig_plus': # positive sigmoid for positive values only, zero elsewhere
return 1./(1+np.exp(-np.clip((x-thr)/x0,-60, 60))) * (x>thr)
elif af=='symsig': # symmetric sigmoid
"""Clipping at 11 is to avoid overflow"""
return 2./(1+np.exp(-np.clip((x-thr)/x0,-60, 60) )) -1
elif af=='bell': # capped linear
"""Clipping at 11 is to avoid overflow"""
return np.exp(- np.clip((x-thr)/x0,-60, 60)**2)
elif af=='camel':
"""Clipping at 11 is to avoid overflow"""
return np.exp(- np.clip((x-thr)/x0,-60, 60)**2 + 1) * (x/x0)**2
elif af=='switch':
return np.where(x-thr>x0, 1, 0)
else:
NotImplementedError('Activation function {} is not implemented.'.format(af))
[docs]
def activation_func_deriv(x, af, x0=1., thr=0, **kwargs):
"""
Apply a specified activation function to input data with optional scaling and thresholding.
The valid activation fucntions are:
- 'lin': linear
- 'relu': rectified linear unit
- 'caplin' : same as 'lin' but capped between 0 and 1
- 'caprelu': same as 'relu' but capped between 0 and 1
- 'capquad': quadratic function capped between 0 and 1
- 'capcube': cubic function capped between 0 and 1
- 'capcube': cubic function capped between 0 and 1
- 'capexp': exponential function capped between 0 and 1
- 'tanh': tangent hyperbolic
- `symsig`: an even sigmoidal function centered at t=0 with range between -1 and 1.
- 'possig': sigmoidal function centered at y=0.5 with range between 0 and 1 .
- 'possig_plus': same as 'possig' but only for values larger than threshold and zero otherwise.
- 'bell' : non-normalized Gaussian activation function
- 'camel': a gaussian multiplied by a pure quadratic function, resulting in two identical bumps
at (±1, 1).
- 'switch': a binary activation. Zero for values smaller than thr and 1 otherwise.
Parameters
----------
x : float or ndarray
Input array or value to be transformed
af: str
Name of the activation function to apply (e.g., 'relu', 'tanh', 'possig', etc.)
x0: float
Scaling factor for the input (default: 1.0)
thr: float
Threshold (bias) value to shift the input (default: 0)
kwargs: dict
Additional keyword arguments for specific activation functions
Returns
-------
float or ndarray
The result of applying the chosen activation function to the input
"""
if af=='lin': # linear
return 1/x0
elif af=='relu': # relu
# return np.maximum(0, (x-thr)/x0)
return np.clip(1/x0, 0, None)
elif af=='caprelu': # capped relu
return np.choose( ((x-thr)/x0 < 1)*(x>thr), [0, 1/x0])
elif af=='caplin': # capped linear
return np.choose(abs(x-thr)/x0 < 1, [0, 1/x0])
# elif af=='tanh': # hyperbolic tangent
# return np.tanh((x-thr)/x0)
elif af=='possig': # positive sigmoid
"""Clipping at 11 is to avoid overflow"""
return activation_func(x, af)*(1 - activation_func(x, af)) / x0
# elif af=='possig_plus': # positive sigmoid for positive values only, zero elsewhere
# return 1./(1+np.exp(-np.clip((x-thr)/x0,-60, 60))) * (x>thr)
# elif af=='symsig': # symmetric sigmoid
# """Clipping at 11 is to avoid overflow"""
# return 2./(1+np.exp(-np.clip((x-thr)/x0,-60, 60) )) -1
# elif af=='bell': # capped linear
# """Clipping at 11 is to avoid overflow"""
# return np.exp(- np.clip((x-thr)/x0,-60, 60)**2)
# elif af=='camel':
# """Clipping at 11 is to avoid overflow"""
# return np.exp(- np.clip((x-thr)/x0,-60, 60)**2 + 1) * (x/x0)**2
# elif af=='switch':
# return np.where(x-thr>x0, 1, 0)
else:
NotImplementedError('Activation function {} is not implemented.'.format(af))
[docs]
def integrate(u, taus, A_adj, w_in, w_rec, dt,
J=1., Ju=1., af='possig',
reset_dim=0, u_reset=None):
def activation(x):
if af=='possig':
return 1./(1+np.exp(-np.clip(x, -11, 11)))
elif af == 'tanh':
return np.tanh(x)
elif af == 'relu':
return np.maximum(0, x)
else:
return x
dt_over_taus = dt / taus
rec_adj = J * A_adj._with_data(w_rec)
in_adj = Ju * w_in
def step(vt, ut):
return vt + dt_over_taus * (rec_adj @ activation(vt) + in_adj @ ut - vt)
# input value at which an action (decision) is to be made
if u_reset is None:
if reset_dim is None:
u_reset = np.nan
else:
u_reset = u[:,reset_dim].max()
N = A_adj.shape[0]
vs = np.empty((len(u), N))
v = np.zeros(N)
for t, ut in enumerate(tqdm(u)):
if ut[reset_dim] == u_reset:
# reset state upon fixation
v *= 0
else:
# propagate integration
v = step(v, ut)
vs[t] = v
# print(np.linalg.norm(ut), np.linalg.norm(v))
return vs
[docs]
def classifier(
name='logreg',
cv=10, random_state=0, #l1_ratios=(0,),
solver='newton-cg', max_iter=200,
scoring='accuracy', C=1,
#scoring='f1_micro', mode = 'binary',
**kwargs):
if name=='logreg':
if cv > 0:
return LogisticRegressionCV(
cv=cv, #l1_ratios=l1_ratios,
solver=solver,
max_iter=max_iter,
scoring=scoring,
use_legacy_attributes=True, # supresses certain deprecation warning
**kwargs
)
else:
return LogisticRegression(
solver=solver,
max_iter=max_iter,
C=C,
)
elif name=='sv':
return SVC(break_ties=True, **kwargs)
elif name=='ridgeclf':
return RidgeClassifier(**kwargs)
# if mode =='binary':
# elif mode =='ovr':
# return OneVsRestClassifier(
# LogisticRegressionCV(cv=cv, l1_ratios=l1_ratios,
# max_iter=max_iter, scoring='f1_micro',
# **kwargs)#multi_class='multinomial')
# )
# else:
# return OneVsOneClassifier(
# LogisticRegressionCV(cv=cv, l1_ratios=l1_ratios,
# max_iter=max_iter,
# **kwargs)#multi_class='multinomial')
# )
[docs]
def unstream(quant, Tsample, Nsample):
if np.ndim(quant)==3:
# assumes the shape is (Nsample, Tsample, num_features)
return quant
elif np.ndim(quant)==2:
# assumes the shape is (Nsample*Tsample, num_features)
return np.stack([quant[n*Tsample: (1+n)*Tsample] for n in range(Nsample)])
[docs]
def plot_taus(tau_hom, tau_het, savepath, savename):
plt.figure()
plt.hist(tau_hom, histtype='step')
plt.hist(tau_het, histtype='step')
plt.xlabel(r'$\tau$')
plt.ylabel('Count')
plt.yscale('log')
plt.savefig(osjoin(savepath, f'{savename}.png'))
plt.close('all')
[docs]
def plot_dataset_hist(y, savepath='.', savename=''):
plt.hist(y)
plt.xlabel('Class ID')
plt.ylabel('Count')
plt.grid()
plt.tight_layout()
plt.savefig(osjoin(savepath, f'histo_{savename}.png'))
plt.close('all')
[docs]
def plot_u_stream(u, lengths, input_labels=None,
figsize = (8,5), tmin=0, tmax=500,
savepath='.', savename='',
):
fig, ax = plt.subplots(1,1, sharex=True, figsize=figsize)
if lengths is None:
xmarks = []
elif isinstance(lengths, (list, np.ndarray)):
xmarks = np.cumsum(np.concatenate(([0], lengths)))
elif isinstance(lengths, int):
xmarks = np.arange(0, len(u), lengths)
else:
pass
if len(xmarks) > 0:
xmarks = xmarks[(max(0, tmin)<=xmarks) *
(xmarks<min(tmax, len(u)))]
_u = u[max(0, tmin): min(tmax, len(u))]
m = ax.pcolormesh(_u.T, cmap='rainbow')
plt.colorbar(m, ax=ax, label='input')
ax.plot(xmarks, 1.02*np.ones_like(xmarks), # slightly below axis
marker='v', markersize=5, color='k', ls='',
transform=ax.get_xaxis_transform(), clip_on=False)
if input_labels is not None:
ax.set_yticks(np.arange(input_labels)+.5, input_labels)
ax.set_ylabel('Stim channel')
ax.set_xlabel('Time')
# ax.set_xlim(tmin, tmax)
# plt.tight_layout()
# plt.savefig(osjoin(savepath, f'u_{savename}.png'))
# plt.close('all')
[docs]
def plot_uvx_stream(u, v, x,
labels=['input', 'voltage', 'activation'],
cmaps = ['rainbow', 'PRGn', 'Grays'],
input_labels=[],
figsize = (8,8), tmin=0, tmax=500,
lengths = None,
savepath='.', savename='',
):
fig, axs = plt.subplots(3,1, sharex=True, figsize=figsize)
if lengths is None:
xmarks = []
elif isinstance(lengths, (list, np.ndarray)):
xmarks = np.cumsum(np.concatenate(([0], lengths)))
elif isinstance(lengths, int):
xmarks = np.arange(0, len(u), lengths)
else:
pass
if lengths is not None:
if tmax < np.cumsum(lengths[:3])[-1]:
tmax = np.cumsum(lengths[:3])[-1] + 20
if len(xmarks) > 0:
xmarks = xmarks[(max(0, tmin)<=xmarks) *
(xmarks<min(tmax, len(u)))]
# if trial_switch_at is not None:
# if isinstance(trial_switch_at, int):
# # interpreted as the length of each trial
# xmarks = np.arange(1, tmax//trial_switch_at + 1) * trial_switch_at
# elif isinstance(trial_switch_at, (list, np.ndarray)):
# # interpreted as the index of trials onset
# trial_switch_at = np.array(trial_switch_at)
# xmarks = trial_switch_at[trial_switch_at<tmax]
# # end_of_go = np.argwhere(np.diff(u[:,0])<0).squeeze() + 1
# # end_of_go = np.argwhere(np.diff(trial_switch_at)<0).squeeze()+1
for ax, quant, label, cmap in zip(axs, [u,v,x], labels, cmaps):
_q = quant[max(0, tmin): min(tmax, len(u))]
# vmin, vmax = None, None
# if 'volt' in label:
# mx = np.quantile(abs(_q), .95)#.max()
# vmin, vmax = -mx, mx
# if 'input' in label:
# vmin, vmax = np.quantile(_q, [0.01, .99])
mx = np.quantile(abs(_q), .95)#.max()
vmin, vmax = -mx, mx
m = ax.pcolormesh(_q.T, cmap=cmap, vmin=vmin, vmax=vmax)
plt.colorbar(m, ax=ax, label=label)
if len(xmarks) > 0:
ax.plot(xmarks, 1.02*np.ones_like(xmarks), # slightly below axis
marker='v', markersize=5, color='k', ls='',
transform=ax.get_xaxis_transform(), clip_on=False)
if len(input_labels) > 0:
udim = u.shape[1]
axs[0].set_yticks(np.arange(udim)+.5, input_labels)
axs[0].set_ylabel('Stim channel')
axs[1].set_ylabel('Neuron ID')
axs[2].set_ylabel('Neuron ID')
axs[2].set_xlabel('Time')
axs[2].set_xlim(tmin, tmax)
plt.tight_layout()
plt.savefig(osjoin(savepath, f'{savename}_uvx.png'))
plt.close('all')
[docs]
def plot_manifold(u, v, x, y, lengths,
manifold='pca', standarized=True,
figsize = (8,3), savepath='.', savename='',
seed = 2026,
):
"""
Gets stream-version of inputs u, voltages v, and activations x, as well
as the class label of each sample y, breaks the streamed arrays into
a flatten array containing fearues of a single sample, and runs a
manifold algorithm for cluster identificaiton. Points are colored by their
class. Array of classes should be ordinal numbers.
"""
if manifold=='tsne':
man = TSNE(n_components=3, random_state=seed)
elif manifold == 'isomap':
man = Isomap()
elif manifold == 'lle':
man = LocallyLinearEmbedding()
elif manifold == 'pca':
man = PCA(n_components=3, )
else:
raise ValueError(f'manifold {manifold} not recognized.')
Nsample = len(y)
is1d = u.shape[-1] == 1
print(Nsample, y.shape, y)
if not isinstance(y[0], int):
ord_enc = OrdinalEncoder()
y = ord_enc.fit_transform(y[:, np.newaxis])
if isinstance(lengths, int):
# interpreted as the length of each sample
# Tsample = x.shape[0]//len(y)
_u = unstream(u, lengths, Nsample)#.reshape(Nsample, -1)
_v = unstream(v, lengths, Nsample)#.reshape(Nsample, -1)
_x = unstream(x, lengths, Nsample)#.reshape(Nsample, -1)
elif isinstance(lengths, (list, np.ndarray)):
_, _u =get_snapshot(lengths, nsnapshot=10, arr=u)
_, _v =get_snapshot(lengths, nsnapshot=10, arr=v)
_, _x =get_snapshot(lengths, nsnapshot=10, arr=x)
# onsets = np.cumsum(np.concatenate(([0], lengths)))
# # We sample 10 points from each trial to create arrays of
# # equal lengths
# _u = []
# _v = []
# _x = []
# for i in range(Nsample):
# onset = onsets[i]
# offset = onsets[i+1] if i < Nsample-1 else len(u)
# snapshots = np.linspace(onset, offset, 10, endpoint=False, dtype=int)
# _u.append(np.array([u[s] for s in snapshots]))
# _v.append(np.array([v[s] for s in snapshots]))
# _x.append(np.array([x[s] for s in snapshots]))
# # print(len(_u), len(_v), len(_x), len(trial_switch_at))
# # print(_u[0].shape, _v[0].shape, _x[0].shape)
# # shape = (Nsample-1, 10, num-features)
_u = np.array(_u)
_v = np.array(_v)
_x = np.array(_x)
print(_u.shape) # Nsample x Time x intput_dim
print(_v.shape) # Nsample x Time x network_dim
print(_x.shape) # Nsample x Time x network_dim
T = _u.shape[1]
# _u = np.array(_u).reshape(Nsample, -1)
# _v = np.array(_v).reshape(Nsample, -1)
# _x = np.array(_x).reshape(Nsample, -1)
if is1d:
_u = _u.squeeze() # Nsample x time
# combines trials and time, uses dim as features
_u = _u.reshape(-1, _u.shape[-1])
_v = _v.reshape(-1, _v.shape[-1])
_x = _x.reshape(-1, _x.shape[-1])
print('after reshape')
print(_u.shape) # Nsample x Time x intput_dim
print(_v.shape) # Nsample x Time x network_dim
print(_x.shape) # Nsample x Time x network_dim
if standarized:
ss = StandardScaler()
_u = ss.fit_transform(_u)
_v = ss.fit_transform(_v)
_x = ss.fit_transform(_x)
print('identifying manifolds')
_ut = man.fit_transform(_u)
_vt = man.fit_transform(_v)
_xt = man.fit_transform(_x)
if not is1d:
_ut = _ut.reshape(Nsample, T, -1)
_vt = _vt.reshape(Nsample, T, -1)
_xt = _xt.reshape(Nsample, T, -1)
print(y.shape, _ut.shape, _vt.shape, _xt.shape)
c = plt.cm.tab20(y/y.max())[:Nsample]
s = np.linspace(2, 5, T)
fig, axs = plt.subplots(1, 3, figsize=figsize, subplot_kw={'projection': '3d'})
# axs[0].scatter(*_ut.T, c=c, s=2)
# axs[1].scatter(*_vt.T, c=c, s=2)
# axs[2].scatter(*_xt.T, c=c, s=2)
print(y.shape, _ut.shape, _vt.shape, _xt.shape)
for i in range(Nsample):
# axs[0].scatter(*_ut[i].T, c=c[i], s=s)
# axs[1].scatter(*_vt[i].T, c=c[i], s=s)
# axs[2].scatter(*_xt[i].T, c=c[i], s=s)
if not is1d:
axs[0].plot(*_ut[i].T, c=c[i], lw=.1)
axs[1].plot(*_vt[i].T, c=c[i], lw=.1)
axs[2].plot(*_xt[i].T, c=c[i], lw=.1)
if not is1d:
axs[0].scatter(*_ut[-1,0,:].T, c='g', s=1, )
axs[1].scatter(*_vt[-1,0,:].T, c='g', s=1, )
axs[2].scatter(*_xt[-1,0,:].T, c='g', s=1, )
if not is1d:
axs[0].scatter(*_ut[-1,-1,:], c='r', s=1,)
axs[1].scatter(*_vt[-1,-1,:], c='r', s=1,)
axs[2].scatter(*_xt[-1,-1,:], c='r', s=1,)
axs[0].set_title('inputs')
axs[1].set_title('voltages')
axs[2].set_title('activations')
for ax in axs:
ax.set_xlabel(f'{manifold} 1')
ax.set_ylabel(f'{manifold} 2')
ax.set_zlabel(f'{manifold} 3')
plt.tight_layout()
plt.savefig(osjoin(savepath, f'{manifold}_{savename}.png'))
plt.close('all')
[docs]
def plot_score_summary(score_trn_hom,
score_trn_het,
score_tst_hom,
score_tst_het,
savepath, savename,
chance_score=None,
pooling_funcs=None,
Cregs=None
):
if pooling_funcs is None:
pooling_funcs = ['']
score_trn_hom = score_trn_hom[:, np.newaxis]
score_trn_het = score_trn_het[:, np.newaxis]
score_tst_hom = score_tst_hom[:, np.newaxis]
score_tst_het = score_tst_het[:, np.newaxis]
for pool_idx, pooling_func in enumerate(pooling_funcs):
fig = plt.figure(figsize=(4,3))
ax = plt.gca()
ax.hist(score_trn_hom[:, pool_idx], histtype='step', color=hetero_colors[0], ls='--', density=True)
ax.hist(score_trn_het[:, pool_idx], histtype='step', color=hetero_colors[3], ls='--', density=True)
ax.hist(score_tst_hom[:, pool_idx], histtype='step', color=hetero_colors[0], label='homogeneous', density=True, )
ax.hist(score_tst_het[:, pool_idx], histtype='step', color=hetero_colors[3], label='heterogeneous', density=True)
ax.plot([],[], color='k', ls='--', label='train')
ax.plot([],[], color='k', label='test' )
if chance_score is not None:
ax.axvline(chance_score, color='g', label='chance score')
ax.legend(ncols=3)
ax.set_xlabel('Classification score')
ax.set_ylabel('Denstiy')
plt.tight_layout()
plt.savefig(osjoin(savepath, f'score_summary_{savename}_{pooling_func}.png'))
[docs]
def plot_score_summary_of_spaced_hom(score_trn_hom, score_trn_het,
score_tst_hom, score_tst_het,
spaced_taus_hom, pooling_funcs, Cregs,
savepath, savename,
chance_score=None,
):
ntaus = len(spaced_taus_hom)
print(Cregs)
# score shape: (Ntau0, Nnet, Npool, Cregs)
for hom, het, name in zip([score_trn_hom, score_tst_hom],
[score_trn_het, score_tst_het],
['trn', 'tst']):
fig, axs = plt.subplots(len(pooling_funcs), len(Cregs),
sharex=True, sharey=True,
figsize=(16, 5.5))
for cidx, Creg in enumerate(Cregs):
for pidx, pooling_func in enumerate(pooling_funcs):
ax = axs[pidx, cidx]
ax.set_title(f'C={Creg:.2f} | pooling = {pooling_func}')
hom_mean = hom[:, :, pidx, cidx].mean(axis=1)
het_mean = het[:, :, pidx, cidx].mean(axis=1)
hom_std = hom[:, :, pidx, cidx].std(axis=1)
het_std = het[:, :, pidx, cidx].std(axis=1)
x0 = np.arange(ntaus)
ax.bar(x0-.15, hom_mean, width=.3, color=hcmap.colors[0], label='homog')
ax.bar(x0+.15, het_mean, width=.3, color=hcmap.colors[3], label='heterog')
ax.errorbar(x0-.15, hom_mean, yerr=hom_std, color='k', ls='')
ax.errorbar(x0+.15, het_mean, yerr=het_std, color='k', ls='')
if chance_score is not None:
ax.axhline(chance_score, label='chance score', ls='--', color='k')
axs[0,0].legend(loc='upper right')
for ax in axs[-1, :]:
ax.set_xticks(np.arange(ntaus), np.arange(0.1, 1.1, 1/ntaus).round(1))
ax.set_xlabel(r'$E[\tau]/T_{\max}$')
for ax in axs[:, 0]:
ax.set_ylim(0, 1.05)
ax.set_ylabel('score')
polish(fig, axs, )
# plt.tight_layout()
plt.savefig(osjoin(savepath, f'score_summary_spaced_{name}_{savename}.png'),
bbox_inches='tight')
# def plot_score_summary_of_spaced_hom(score_trn_hom, score_trn_het,
# score_tst_hom, score_tst_het,
# spaced_taus_hom,
# pooling_funcs,
# savepath, savename,
# chance_score=None,
# ):
# for pool_idx, pooling_func in enumerate(pooling_funcs):
# fig, axs= plt.subplots(2, 1, sharex=True, )
# ntaus = len(spaced_taus_hom)
# for ax, name, shom, shet in zip(axs, ['train','test'],
# [score_trn_hom[..., pool_idx], score_tst_hom[..., pool_idx]],
# [score_trn_het[..., pool_idx], score_tst_het[..., pool_idx]],
# ):
# hom_mean = shom.mean(axis=1)
# het_mean = shet.mean(axis=1)
# hom_std = shom.std(axis=1)
# het_std = shet.std(axis=1)
# x0 = np.arange(ntaus)
# ax.bar(x0-.1, hom_mean, width=.2, color=hcmap.colors[0], label='homog')
# ax.bar(x0+.1, het_mean, width=.2, color=hcmap.colors[3], label='heterog')
# ax.errorbar(x0-.1, hom_mean, yerr=hom_std, color='k', ls='')
# ax.errorbar(x0+.1, het_mean, yerr=het_std, color='k', ls='')
# ax.set_ylabel(f'{name} score')
# ax.set_ylim(0, 1.05)
# if chance_score is not None:
# ax.axhline(chance_score, label='chance score', ls='--', color='k')
# axs[1].set_xticks(np.arange(10), np.arange(0.1, 1.1, .1).round(1))
# axs[1].set_xlabel(r'$E[\tau]/T_{\max}$')
# fig.suptitle(f'"{pooling_func}" pooling')
# plt.legend(ncols=3, loc='best')
# plt.tight_layout()
# plt.savefig(osjoin(savepath, f'score_summary_spaced_{pooling_func}_{savename}.png'))
[docs]
def get_snapshot(lengths, nsnapshot=10, arr=None, flatten=False):
N_smaples = len(lengths)
onsets = np.cumsum(np.concatenate(([0], lengths)))
idxs = []
for i in range(N_smaples):
onset = onsets[i]
offset = onsets[i+1] if i < N_smaples-1 else len(arr)
# might include repetative index, but that's inevitable for
# very short trials. Also, we omit the onset as it is where
# the voltages are being reset. Instead, we increase the
# num_snapshop by one.
idxs.append(np.linspace(onset, offset, nsnapshot+1,
endpoint=False, dtype=int)[1:])
if arr is not None:
snapped = np.stack([arr[i] for i in idxs])
if flatten:
snapped = snapped.reshape(N_smaples, -1)
# snapped.append(np.array([arr[s] for s in snaptimes]))
else:
snapped = None
# return np.array(snapped).reshape(N_smaples, -1)
return idxs, snapped
# def get_pool(lengths, arr, func='max', n=-1, nsnapshot=10, flatten=False):
# N_smaples = len(lengths)
# onsets = np.cumsum(np.concatenate(([0], lengths))) # with endpoint
# offsets = onsets[1:]
# onsets = onsets[:-1]
# idxs = []
# pooled = []
# for i in range(N_smaples):
# onset = onsets[i]
# offset = offsets[i]
# # offset = onsets[i+1] if i < N_smaples-1 else len(arr)
# sub = arr[onset:offset]
# if isinstance(func, str):
# if func == 'max':
# pooled.append(sub.max(axis=0))
# idxs.append(np.argmax(sub, axis=0))
# elif func == 'min':
# pooled.append(sub.min(axis=0))
# idxs.append(np.argmin(sub, axis=0))
# elif func == 'avg' or func == 'mean':
# pooled.append(sub.mean(axis=0))
# idxs.append([])
# elif func == 'std':
# pooled.append(sub.std(axis=0))
# idxs.append([])
# elif func == 'last':
# pooled.append(sub[-1])
# idxs.append([])
# elif func == 'first': # don't use this!
# pooled.append(sub[0])
# idxs.append([])
# elif func == 'nth':
# pooled.append(sub[n])
# idxs.append([])
# elif func == 'snapshot':
# snap_idxs = np.linspace(onset, offset, nsnapshot+1, endpoint=False, dtype=int)
# snap_idxs = snap_idxs[1:] # dropping the onset snapshot
# snaped = np.stack([sub[i] for i in snap_idxs])
# pooled.append(snaped)
# idxs.append(snap_idxs)
# else:
# raise ValueError(f'func = {func} is not a valid pooling function.')
# else:
# pool_feature = []
# pool_index = []
# for _func in func:
# if _func == 'max':
# pool_feature.append(sub.max(axis=0))
# pool_index.append(np.argmax(sub, axis=0))
# elif _func == 'min':
# pool_feature.append(sub.min(axis=0))
# pool_index.append(np.argmin(sub, axis=0))
# elif _func == 'avg' or func == 'mean':
# pool_feature.append(sub.mean(axis=0))
# pool_index.append([])
# elif _func == 'std':
# pool_feature.append(sub.std(axis=0))
# pool_index.append([])
# elif _func == 'last':
# pool_feature.append(sub[-1])
# pool_index.append([])
# elif _func == 'first': # don't use this!
# pool_feature.append(sub[0])
# pool_index.append([])
# elif _func == 'nth':
# pool_feature.append(sub[n])
# pool_index.append([])
# elif func == 'snapshot':
# snap_idxs = np.linspace(onset, offset, nsnapshot+1, endpoint=False, dtype=int)
# snap_idxs = snap_idxs[1:] # dropping the onset snapshot
# snaped = np.stack([sub[i] for i in snap_idxs])
# pool_feature.append(snaped)
# pool_index.append(snap_idxs)
# else:
# raise ValueError(f'func = {func} is not a valid pooling function.')
# pooled.append(np.concat(pool_feature))
# idxs.append(np.concat(pool_index))
# pooled = np.stack(pooled)
# if flatten:
# pooled = pooled.reshape(N_smaples, -1)
# return idxs, pooled
[docs]
def pool(arr, trial_lengths, func='max', n=-1, nsnapshot=10, golength=10, flatten=True):
onsets = np.cumsum(np.concatenate(([0], trial_lengths))) # with endpoint
offsets = onsets[1:]
onsets = onsets[:-1]
# sanity checks
if func=='go':
assert golength is not None
elif func == 'snapshpt':
assert nsnapshot is not None
elif func == 'nth':
assert n is not None
pooled = []
for onset, offset in zip(onsets, offsets):
sub = arr[onset:offset]
if isinstance(func, str):
if func == 'max':
pooled.append(sub.max(axis=0))
elif func == 'min':
pooled.append(sub.min(axis=0))
elif func == 'avg' or func == 'mean':
pooled.append(sub.mean(axis=0))
elif func == 'std':
pooled.append(sub.std(axis=0))
elif func == 'last':
pooled.append(sub[-1])
elif func == 'first': # don't use this!
pooled.append(sub[0])
elif func == 'nth':
pooled.append(sub[n])
elif func == 'go':
pooled.append(sub[-golength:])
elif func == 'snapshot':
snap_idxs = np.linspace(0, len(sub), nsnapshot+1, endpoint=False, dtype=int)
snap_idxs = snap_idxs[1:] # dropping the onset snapshot
snaped = np.stack([sub[i] for i in snap_idxs])
pooled.append(snaped)
else:
raise ValueError(f'func = {func} is not a valid pooling function.')
else:
pool_feature = []
for _func in func:
if _func == 'max':
pool_feature.append(sub.max(axis=0))
elif _func == 'min':
pool_feature.append(sub.min(axis=0))
elif _func == 'avg' or func == 'mean':
pool_feature.append(sub.mean(axis=0))
elif _func == 'std':
pool_feature.append(sub.std(axis=0))
elif _func == 'last':
pool_feature.append(sub[-1])
elif _func == 'first': # don't use this!
pool_feature.append(sub[0])
elif _func == 'nth':
pool_feature.append(sub[n])
elif func == 'go':
pool_feature.append(sub[-golength:])
elif func == 'snapshot':
snap_idxs = np.linspace(0, len(sub), nsnapshot+1, endpoint=False, dtype=int)
snap_idxs = snap_idxs[1:] # dropping the onset snapshot
snaped = np.stack([sub[i] for i in snap_idxs])
pooled.append(snaped)
else:
raise ValueError(f'func = {func} is not a valid pooling function.')
pooled.append(np.concat(pool_feature))
pooled = np.stack(pooled)
if flatten:
N_smaples = len(trial_lengths)
pooled = pooled.reshape(N_smaples, -1)
return pooled
[docs]
def get_trial_start_end(lens):
onsets = np.cumsum(np.concatenate(([0], lens))) # with endpoint
offsets = onsets[1:]
onsets = onsets[:-1]
return onsets, offsets
[docs]
def get_trial_stack(arr, lens):
onsets = np.cumsum(np.concatenate(([0], lens))) # with endpoint
offsets = onsets[1:]
onsets = onsets[:-1]
return [arr[on:off] for on, off in zip(onsets, offsets)]
[docs]
def get_statistics(us, lens, per_length=True):
trials = get_trial_stack(us, lens)
means = np.array([np.mean(trial, axis=0) for trial in trials])
stds = np.array([np.std(trial, axis=0) for trial in trials])
norms = np.array([np.linalg.norm(trial, axis=0) for trial in trials])
if per_length:
means /= lens[:, np.newaxis]
stds /= lens[:, np.newaxis]
norms /= lens[:, np.newaxis]
return means, stds, norms
# def get_len_quantiles(lens, num_qs=5):
# qs = np.linspace(0, 1, num_qs) # dropping 0
# return np.quantile(lens, qs).astype(int)
[docs]
def set_het_level(value='paper', lens=None,):
if value == 'paper':
h = np.sqrt(10) # similar to the paper
else:
assert lens is not None
if value == 'minlen':
h = lens.min()
elif value == 'maxlen':
h = lens.max()
elif value == 'midlen':
h = (lens.max() - lens.min())/2
return h
[docs]
def plot_taus_summary(het_taus, hom_taus, savepath, savename='tau_summary'):
Ntaus, Nnets, N = het_taus.shape
assert Ntaus == len(hom_taus)
fig, ax = plt.subplots(1, 1, figsize=(8,2.6))
for tau_idx, hom_tau in enumerate(hom_taus):
ax.hist(het_taus[tau_idx].ravel(), bins=100,
histtype='step', color=f'C{tau_idx}', zorder=-tau_idx)
ax.axvline(hom_tau, color=f'C{tau_idx}', ls='--')
# ax.set_xscale('log')
# ax.set_yscale('log')
ax.set_xlabel(r'$\tau$')
ax.set_ylabel('Frequency')
plt.savefig(osjoin(savepath, f'{savename}.png'))
plt.close('all')
[docs]
def collect_taus(savepath, Ntau0, Nnet):
taus = []
for tau0_idx in range(Ntau0):
for net_idx in range(Nnet):
tau_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_taus_het.npy')
taus.append(np.load(tau_path))
return np.array(taus).reshape(Ntau0, Nnet, -1)
[docs]
def collect_scores(savepath, Ntau0, Nnet, Npool, NCregs=1):
trn_homs = []#np.zeros((Ntau0, Nnet, Npool, NCregs))
tst_homs = []#np.zeros((Ntau0, Nnet, Npool, NCregs))
trn_hets = []#np.zeros((Ntau0, Nnet, Npool, NCregs))
tst_hets = []#np.zeros((Ntau0, Nnet, Npool, NCregs))
for tau0_idx in range(Ntau0):
for net_idx in range(Nnet):
trn_homs_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_trn_hom.npy')
tst_homs_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_tst_hom.npy')
trn_hets_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_trn_het.npy')
tst_hets_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_tst_het.npy')
# trn_homs[tau0_idx, net_idx] = np.load(trn_homs_path)
# tst_homs[tau0_idx, net_idx] = np.load(tst_homs_path)
# trn_hets[tau0_idx, net_idx] = np.load(trn_hets_path)
# tst_hets[tau0_idx, net_idx] = np.load(tst_hets_path)
trn_homs.append(np.load(trn_homs_path))
tst_homs.append(np.load(tst_homs_path))
trn_hets.append(np.load(trn_hets_path))
tst_hets.append(np.load(tst_hets_path))
trn_homs = np.array(trn_homs).reshape(Ntau0, Nnet, Npool, -1).squeeze()
tst_homs = np.array(tst_homs).reshape(Ntau0, Nnet, Npool, -1).squeeze()
trn_hets = np.array(trn_hets).reshape(Ntau0, Nnet, Npool, -1).squeeze()
tst_hets = np.array(tst_hets).reshape(Ntau0, Nnet, Npool, -1).squeeze()
return trn_homs, tst_homs, trn_hets, tst_hets
[docs]
def sample_integrate_classify_save(
dt, tau0,
h_value, h_is_multip, #tau_min, tau_max,
N, p, J, Ju,
us_trn, us_tst,
lens_trn, lens_tst,
ys_trn, ys_tst,
pooling_funcs, pooling_kwargs, af,
cv, Creg, clf_name,
prefix, savepath,
plot_uvx=True, plot_man=True,
net_idx=0, tau0_idx=0,
reset_dim=0, u_reset=None):
trn_homs_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_trn_hom.npy')
tst_homs_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_tst_hom.npy')
trn_hets_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_trn_het.npy')
tst_hsts_path = osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_tst_het.npy')
if exists(trn_homs_path) and exists(tst_homs_path) and \
exists(trn_hets_path) and exists(tst_hsts_path):
print(f'tau_{tau0_idx}_net_{net_idx}: Exists. SKIPPED!')
return None
# sample_taus
tau_hom = np.ones(N) * tau0
if h_value == 'uniform':
# tau_min = max(tau_min* lens_trn.max(), 1)
# tau_max = min(tau_max* lens_trn.max(), lens_trn.max())
# tau_het = np.random.uniform(tau_min, tau_max, N) * dt
tau_het = np.random.uniform(1, lens_trn.max(), N) * dt
else:
h = set_het_level(value=h_value, lens=lens_trn)
std= h*tau0 if h_is_multip else h*dt
tau_het = sample_taus(mu=tau0, std=std, size=N,
clip_below=dt, savename=None)
print(tau_hom.mean(), tau_het.mean(), tau_het.std())
# initialize net
net_settings = dict(N=N, p=p, N_input=us_trn.shape[1], savename=None)
A_adj, w_in, w_rec = generate_net(**net_settings)
# w_in /= us_trn.std(axis=0) # scale the network
# integrating
int_settings = dict(A_adj=A_adj, w_in=w_in, w_rec=w_rec,
af=af, J=J, Ju=Ju, dt=dt,
reset_dim=reset_dim, u_reset=u_reset)
trn_homs, tst_homs, trn_hets, tst_hets = [], [], [], []
if isinstance(Creg, (int, float)):
Creg = [Creg]
trn_homs = np.zeros((len(pooling_funcs), len(Creg)))
tst_homs = np.zeros((len(pooling_funcs), len(Creg)))
trn_hets = np.zeros((len(pooling_funcs), len(Creg)))
tst_hets = np.zeros((len(pooling_funcs), len(Creg)))
for taus, name in zip([tau_hom, tau_het], ['hom','het']):
vs_trn = integrate(u=us_trn, taus=taus, **int_settings).astype(np.float32)
vs_tst = integrate(u=us_tst, taus=taus, **int_settings).astype(np.float32)
xs_trn = activation_func(vs_trn, af).astype(np.float32)
xs_tst = activation_func(vs_tst, af).astype(np.float32)
# plotting
if plot_uvx:
if net_idx == 0:
# input_labels = get_input_channel_names(us_trn)
uvx_name = f'{prefix}_tau{tau0_idx}'
plot_uvx_stream(us_trn, vs_trn, xs_trn, lengths=lens_trn,
tmax=1500, #input_labels = input_labels,
savepath=savepath, savename=f'{uvx_name}_{name}')
# plot_uvx_stream(us_tst, vs_tst, xs_tst, lengths=lens_tst,
# tmax=1500, #input_labels = input_labels,
# savepath=savepath, savename=f'{uvx_name}_tst_{name}')
if plot_man:
if (tau0_idx == 0 or tau0_idx == 9) and net_idx == 0:
man_name = f'{prefix}_tau{tau0_idx}'
plot_manifold(us_trn, vs_trn, xs_trn, ys_trn, lens_trn,
savepath=savepath, savename=f'{man_name}_{name}')
# free up space
del vs_trn, vs_tst
for pool_idx, func in enumerate(pooling_funcs):
# pooling states
pooling_config = pooling_kwargs
pooling_config.update(dict(func=func, ))
Xs_trn = pool(xs_trn, trial_lengths=lens_trn, **pooling_config)
Xs_tst = pool(xs_tst, trial_lengths=lens_tst, **pooling_config)
for cidx, C in enumerate(Creg):
clf = classifier(cv=cv, C=C, name=clf_name)
clf.fit(Xs_trn, ys_trn)
score_trn = clf.score(Xs_trn, ys_trn)
score_tst = clf.score(Xs_tst, ys_tst)
# report
print(f'\t {func}- tau{tau0_idx} | net{net_idx} | C={Creg[cidx]} -- {name}: {score_tst}')
eval(f'trn_{name}s')[pool_idx, cidx] = score_trn
eval(f'tst_{name}s')[pool_idx, cidx] = score_tst
del xs_trn, xs_tst, Xs_trn, Xs_tst
np.save(osjoin(savepath, f'tau_{tau0_idx}_net_{net_idx}_taus_het.npy'), tau_het)
np.save(trn_homs_path, trn_homs)
np.save(tst_homs_path, tst_homs)
np.save(trn_hets_path, trn_hets)
np.save(tst_hsts_path, tst_hets)
print(f'tau_{tau0_idx}_net_{net_idx}: DONE!')