In [116]:
from jedi import jedi
from jedi.utils import plot, seedutil, func_generator, init_tools
from functools import partial
import random
import types
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from __future__ import division
from scipy.integrate import odeint, ode
from numpy import zeros,ones,eye,tanh,dot,outer,sqrt,linspace, \
cos,pi,hstack,zeros_like,abs,repeat
from numpy.random import uniform,normal,choice
import numpy.linalg as la
from scipy.optimize import minimize
from ipywidgets import interact, fixed
from sklearn.decomposition import PCA
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
In [2]:
reload(jedi)
reload(seedutil)
Out[2]:
In [3]:
# Setting Seeds
seeds = random.sample(range(100000), 1)
In [4]:
t = linspace(-10,10,1000)
plt.subplot(2,1,1)
plt.plot(t, jedi.step_decode(t),'k', lw=4);
plt.ylim(-0.5,1.5)
delta = 2
plt.subplot(2,1,2)
plt.plot(t, jedi.sigmoid(delta, t), 'k', lw=4);
plt.ylim(-.5,1.5);
In [5]:
reload(jedi)
reload(seedutil)
Out[5]:
In [6]:
# sine-wave target
target = lambda t0: cos(2 * pi * t0/.5)
In [8]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 6 # simulation length
tstart = 2 # learning start time
tstop = 4 # learning stop time
rho = 1.5 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
In [23]:
errors = []
for seedling in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
tanh_x = params['tanh_x']
z = params['z']
return (-x + dot(J, tanh_x) + Wz*z)/dt
x, t, z, _, wu,_ = jedi.force(target, model, lr, dt, tmax, tstart, tstop, x0, w, 0)
error = np.abs(z-target(t))
errors.append(error)
errors = np.array(errors)
In [ ]:
In [24]:
# Visualizing activities of first 20 neurons
T = 600
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(10):
ax.plot(t[:T], x[:T, i])
In [25]:
plt.vlines?
In [26]:
errs = errors
burn_in = 1
title = ""
plt.figure(figsize=(12,3))
plt.xlim(0,6)
mean = True
if mean:
errs = np.mean(errs, axis=0)
ymax = 2*np.max(errs[burn_in:])
plt.plot(t[burn_in:], errs[burn_in:], 'k', label="Signal/Output Error")
plt.vlines(tstart,0, ymax, colors='g', label='Train Start')
plt.vlines(tstop,0, ymax, colors='r', label='Train Stop')
plt.ylim(0,ymax)
#plt.xlabel('time', fontweight='bold', fontsize=16)
plt.ylabel('error', fontsize=30)
plt.title(title, fontweight='bold', fontsize=20)
plt.legend();
In [27]:
plt.figure(figsize=(12,4))
plot.target_vs_output_plus_error(t, z, wu, target, offset=1, log=False)
In [28]:
offset = 1; log = False; ylim=None; xlim=None
plt.figure(figsize=(12,5))
plt.subplot(2, 1, 1)
if isinstance(target, types.FunctionType):
plt.plot(t, target(t), '-k', lw=2)
else:
plt.plot(t[offset:], target, '-k', lw=2)
plt.plot(t, z, '-b')
plt.legend(('target', 'output'))
if ylim is None:
ylim = [-1.1,3]
plt.ylim(ylim)
if xlim is None:
plt.xlim([min(t), max(t)])
plt.xticks([])
plt.subplot(2, 1, 2)
plt.plot(t, wu, '-r')
if log:
plt.yscale('log')
if xlim is None:
plt.xlim([min(t), max(t)])
plt.ylabel('$|\Delta w|$', fontsize=30)
plt.xlabel('time', fontsize=30)
Out[28]:
In [29]:
derrors = []
rho = 1.5
for seed in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, p=sparsity, rho=rho)
def model(t0, x, params):
tanh_x = params['tanh_x']
z = params['z']
return (-x + dot(J, tanh_x) + Wz*z)/dt
x, t, z, _, wu,_ = jedi.dforce(jedi.step_decode, target, model, lr, dt, tmax, tstart, tstop, x0, w, 0, pE=pE)
derror = np.abs(z-target(t))
derrors.append(derror)
derrors = np.array(derrors)
In [30]:
# Visualizing activities of first 20 neurons
T = 600
plt.figure(figsize=(12,3))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for i in range(10):
ax.plot(t[:T], x[:T, i])
In [31]:
errs = derrors
burn_in = 1
title = ""
plt.figure(figsize=(12,3))
plt.xlim(0,6)
mean = True
if mean:
errs = np.mean(errs, axis=0)
ymax = 2*np.max(errs[burn_in:])
plt.plot(t[burn_in:], errs[burn_in:], 'k', label="Signal/Output Error")
plt.vlines(tstart,0, ymax, colors='g', label='Train Start')
plt.vlines(tstop,0, ymax, colors='r', label='Train Stop')
plt.ylim(0,ymax)
#plt.xlabel('time', fontweight='bold', fontsize=16)
plt.ylabel('error', fontsize=30)
plt.title(title, fontweight='bold', fontsize=30)
plt.legend();
In [32]:
offset = 1; log = False; ylim=None; xlim=None
plt.figure(figsize=(12,5))
plt.subplot(2, 1, 1)
if isinstance(target, types.FunctionType):
plt.plot(t, target(t), '-k', lw=2)
else:
plt.plot(t[offset:], target, '-k', lw=2)
plt.plot(t, z, '-b')
plt.legend(('target', 'output'))
if ylim is None:
ylim = [-1.1,3]
plt.ylim(ylim)
if xlim is None:
plt.xlim([min(t), max(t)])
plt.xticks([])
plt.subplot(2, 1, 2)
plt.plot(t, wu, '-r')
if log:
plt.yscale('log')
if xlim is None:
plt.xlim([min(t), max(t)])
plt.ylabel('$|\Delta w|$', fontsize=30)
plt.xlabel('time', fontsize=30)
Out[32]:
In [33]:
plt.figure(figsize=(12,5))
plot.signal_error(derrors, t, tstart, tstop, title="DFORCE (Sin Wave)", burn_in=5)
plot.target_vs_output_plus_error(t, z, wu, target, log=False)
In [34]:
plt.figure(figsize=(12,3))
plt.xlim(1,6)
plot.cross_signal_error(errors, derrors, t, tstart, tstop,
title="FORCE vs SFORCE (Sin Wave)", burn_in=100)
In [81]:
reload(func_generator)
Out[81]:
In [101]:
# Plotting inputs and targets
tmax = 10
dt = .01
n = 1
inputs, targets = func_generator.flip_flop_generator(n=n, spikes=[[5,5] for _ in range(n)], t=tmax, dt=dt)
scale = 3
plt.figure(figsize=(12, scale*n))
for i in range(n):
plt.subplot(n, 1, i+1)
plt.plot(range(int(tmax/dt)+2), targets[i], label="Target");
plt.plot(range(int(tmax/dt)+2), inputs[i], label="Input");
plt.ylim(-2,2);
plt.legend();
In [102]:
targets = np.load("../data/stability/flipflop/targets_tmax10sec.npy")
inputs = np.load("../data/stability/flipflop/inputs_tmax10sec.npy")
In [103]:
import seaborn
In [104]:
t = np.linspace(0,10,1002)
In [105]:
seaborn.set_style("white")
plt.figure(figsize=(12,5))
plt.ylim(-1.5,3)
plt.xlim(0,10)
plt.xlabel('time', fontsize=35)
ax = plt.subplot(111)
ax.plot(t,targets,'b', label='target')
ax.plot(t, inputs, 'm', label='input')
ax.fill_between(t[501:650], -1.5, 3, alpha=.3, color='r')
ax.fill_between(t[:501], -1.5, 3, alpha=.1, color='g')
ax.vlines(5,-1.5, ymax)
ax.vlines(6.5,-1.5, ymax)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.legend(prop={'size':30})
Out[105]:
In [194]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstart = 0
tstop = 5 # learning stop time
rho = 1.02 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.05,1,1) # sparsity
I = 1 # input size
In [195]:
errors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, Wi, x0, u, w = init_tools.set_simulation_parameters(seedling, N, I, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
inp = params['inputs'][index]
return (-x + dot(J, tanh_x) + dot(Wi, inp) + Wz*z)/dt
x,t,z,w_learn,wu,_ = jedi.force(targets, model, lr, dt, tmax, tstart, tstop, x0, w, inputs=inputs)
zs.append(z)
wus.append(wu)
error = np.abs(z-np.array(targets))
errors.append(error)
errors = np.array(errors)
In [196]:
plt.figure(figsize=(12,5))
plot.signal_error(errors, t, tstart, tstop, title= "FORCE (Flip-Flop)", burn_in=5)
In [197]:
ind = 0
In [198]:
plt.figure(figsize=(12,5))
if ind < len(seeds):
print("Seed: %d" % ind)
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=0, log=False)
ind+=1
In [199]:
plt.figure(figsize=(12,5))
plt.legend()
plt.subplot(2,1,2)
for i in range(20):
plt.plot(t[:], x[:, i]);
plt.subplot(2,1,1)
plt.plot(t, np.array(z), lw=2, label="output");
In [200]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("FORCE (Flip Flop)");
In [201]:
pca = PCA(n_components=2)
pca.fit(x)
pca_x = pca.transform(x).T
In [202]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [203]:
def q(x):
F = lambda x: -x + dot(J, np.tanh(x)) + Wz*dot(w,np.tanh(x))/dt
return .5*(la.norm(F(x))**2)
In [204]:
ics = x[np.random.randint(0,1000,10)]
In [205]:
minima = []
for ic in ics:
minima.append(minimize(q, ic).x)
In [206]:
pca_minima = pca.transform(np.array(minima)).T
pca_ics = pca.transform(np.array(ics)).T
In [207]:
plt.plot(pca_x[0], pca_x[1], 'ko', alpha=.1);
plt.plot(pca_minima[0], pca_minima[1], 'ro', alpha=.8);
plt.plot(pca_ics[0], pca_ics[1], 'bo', alpha=.2);
In [ ]:
In [152]:
res = minimize(q, ics[0])
In [53]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 5 # learning stop time
rho = 1.02 # spectral radius of J
N = 300 # size of stochastic pool
lr = 1.0 # learning rate
pE = .8 # percent excitatory
sparsity = (.1,1,1) # sparsity
I = 1 # input size
In [54]:
derrors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, Wi, x0, u, w = init_tools.set_simulation_parameters(seedling, N, I, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
inp = params['inputs'][index]
return (-x + dot(J, tanh_x) + dot(Wi, inp) + Wz*z)/dt
x,t,z,w_learn,wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w,
pE=.8, inputs=inputs)
zs.append(z)
wus.append(wu)
derror = np.abs(z-np.array(targets))
derrors.append(derror)
derrors = np.array(derrors)
In [56]:
plt.figure(figsize=(12,5))
plot.signal_error(derrors, t, tstart, tstop, title="SFORCE (Flip-Flop)", burn_in=5)
In [57]:
# Setting seed index
ind = 0
In [58]:
plt.figure(figsize=(12,5))
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=0, log=False)
In [59]:
plt.figure(figsize=(12,6))
plt.legend()
plt.subplot(2,1,2)
for i in range(20):
plt.plot(t[:], x[:, i]);
plt.subplot(2,1,1)
plt.plot(t, np.array(z), lw=2, label="output");
In [60]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("DFORCE (Flip Flop)");
In [61]:
pca = PCA(n_components=2)
pca.fit(x)
pca_x = pca.transform(x).T
In [62]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_2dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
In [ ]:
In [63]:
# cross mean signal error
plt.figure(figsize=(12,5))
plot.cross_signal_error(errors, derrors, t, tstart, tstop,
title="FORCE vs SFORCE (Flip-Flop)", burn_in=5)
In [69]:
# Setting seed index
ind = 0
In [70]:
plt.figure(figsize=(12,5))
print("Seed: %d" % ind)
plot.cross_signal_error(errors[ind], derrors[ind], t, tstart, tstop,
title="FORCE vs SFORCE, (Flip-Flop)",
burn_in=5, mean=False)
ind+=1
In [71]:
# Parameters specified by Abbott 2009.
def lorentz((x, y, z), t0, sigma=10., beta=8./3, rho=28.0):
return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
In [72]:
break_in = 500
T = 1501 # period
x0 = np.random.randn(3) # starting vector
t_= np.linspace(0, 60, T)
lorenz = odeint(lorentz, x0, t_)/20
targets = lorenz[break_in:,0]
t_in = t_[break_in:]
In [73]:
# Visualizing Lorentz attractor
plt.figure(figsize=(12,3))
plt.plot(t_in, targets);
plt.xlim([min(t_in), max(t_in)])
plt.ylim([-1.5,1.5])
plt.xlabel('time', fontsize=14);
plt.ylabel('y', fontsize=14);
In [74]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 8 # learning stop time
rho = 1.2 # spectral radius of J
N = 1000 # size of stochastic pool
lr = 1.0 # learning rate
pE = None # percent excitatory
sparsity = (.1,1,1) # sparsity
In [76]:
errors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, Wi, x0, u, w = init_tools.set_simulation_parameters(seedling, N, I, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, params):
index = params['index']
tanh_x = params['tanh_x']
z = params['z']
inp = params['inputs'][index]
return (-x + dot(J, tanh_x) + dot(Wi, inp) + Wz*z)/dt
x,t,z,w_learn,wu,_ = jedi.force(targets, model, lr, dt, tmax, tstart, tstop, x0, w, inputs=inputs)
zs.append(z)
wus.append(wu)
error = np.abs(z-np.array(targets))
errors.append(error)
errors = np.array(errors)
In [ ]:
plt.figure(figsize=(12,7))
ind = 0
tstart = 1
plot.target_vs_output_plus_error(t[tstart:], zs[ind][tstart:],
wus[ind][tstart:], targets[tstart:], offset=1, ylim=[-1,2])
In [ ]:
# Visualizing activities of first 20 neurons
T = 500
K = 20
plt.figure(figsize=(12,4))
plt.subplot(211)
plt.title("Neuron Dynamics");
for i in range(K):
plt.plot(t[:T], x[:T, i]);
plt.subplot(212)
for i in range(K):
plt.plot(t[-T:], x[-T:, i]);
plt.xlim(t[-T], t[-1])
In [ ]:
plt.figure(figsize=(12,5))
plot.signal_error(errors, t[1:], tstop, title= "FORCE", burn_in=5)
In [ ]:
ind = 0
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("FORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=3)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_3dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
In [ ]:
# simulation parameters for FORCE
dt = .01 # time step
tmax = 10 # simulation length
tstop = 8 # learning stop time
rho = 1.02 # spectral radius of J
N = 1000 # size of stochastic pool
lr = 1.0 # learning rate
pE = None # percent excitatory
sparsity = (.1,1,1) # sparsity
In [ ]:
derrors = []
wus = []
zs = []
for seedling in seeds:
J, Wz, _, x0, u, w = init_tools.set_simulation_parameters(seedling, N, 1, pE=pE, p=sparsity, rho=rho)
# inp & z are dummy variables
def model(t0, x, tanh_x, inp, z):
return (-x + dot(J, tanh_x) + Wz*z)/dt
x, t, z, _, wu,_ = jedi.dforce(jedi.step_decode, targets, model, lr, dt, tmax, tstart, tstop, x0, w, 0, pE=.8)
zs.append(z)
wus.append(wu)
derror = np.abs(z[1:]-np.array(targets))
derrors.append(derror)
derrors = np.array(derrors)
In [ ]:
plt.figure(figsize=(12,5))
ind = 0
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=1, log=False)
In [ ]:
plt.figure(figsize=(12,4))
plot.signal_error(derrors, t[1:], tstop, title= "SFORCE", burn_in=5)
In [ ]:
plt.figure(figsize=(12,5))
print("Seed: %d" % ind)
plot.signal_error(derrors[ind], t[1:], tstop,
title= "SFORCE", burn_in=5, mean=False)
plot.target_vs_output_plus_error(t, zs[ind], wus[ind], targets, offset=1, log=False)
ind+=1
In [ ]:
pca = PCA(n_components=20)
X_pca = pca.fit_transform(x)
plt.figure(figsize=(8,5))
plt.plot(100*pca.explained_variance_/sum(pca.explained_variance_), 'bo');
plt.xlabel("PCs");
plt.ylabel("Percent of Variance Explained");
plt.title("SFORCE (Flip Flop)");
In [ ]:
pca = PCA(n_components=3)
pca.fit(x)
pca_x = pca.transform(x).T
In [ ]:
tmax = t[-1]
tmin = t[0]
tv = t[:-1]
interact(plot.visualize_3dim_state, time=(tmin, tmax, .1), pca_x=fixed(pca_x), tv=fixed(tv));
In [ ]:
In [ ]:
# cross mean signal error
plot.cross_signal_error(errors, derrors, t[1:], tstop,
title="FORCE vs SFORCE (Lorentz)", burn_in=5)
In [ ]: