Accompanying Readings: Chapter 6
The story so far:
So far we've just talked about neural activity in a single population
In [3]:
%pylab inline
In [4]:
import numpy as np
import nengo
from nengo.dists import Uniform
from nengo.processes import WhiteSignal
from nengo.solvers import LstsqL2
T = 1.0
max_freq = 10
model = nengo.Network('Communication Channel', seed=3)
with model:
stim = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.5))
ensA = nengo.Ensemble(20, dimensions=1, neuron_type=nengo.LIFRate())
ensB = nengo.Ensemble(19, dimensions=1, neuron_type=nengo.LIFRate())
temp = nengo.Ensemble(1, dimensions=1, neuron_type=nengo.LIFRate())
nengo.Connection(stim, ensA)
stim_B = nengo.Connection(stim, ensB)
connectionA = nengo.Connection(ensA, temp) #This is just to generate the decoders
connectionB = nengo.Connection(ensB, temp) #This is just to generate the decoders
stim_p = nengo.Probe(stim)
a_rates = nengo.Probe(ensA.neurons, 'rates')
b_rates = nengo.Probe(ensB.neurons, 'rates')
sim = nengo.Simulator(model, seed=3)
sim.run(T)
x = sim.data[stim_p]
d_i = sim.data[connectionA].weights.T
A_i = sim.data[a_rates]
d_j = sim.data[connectionB].weights.T
A_j = sim.data[b_rates]
#Add noise
A_i = A_i + np.random.normal(scale=0.2*np.max(A_i), size=A_i.shape)
A_j = A_j + np.random.normal(scale=0.2*np.max(A_j), size=A_j.shape)
xhat_i = np.dot(A_i, d_i)
yhat_j = np.dot(A_j, d_j)
t = sim.trange()
figure(figsize=(8,4))
subplot(1,2,1)
plot(t, xhat_i, 'g', label='$\hat{x}$')
plot(t, x, 'b', label='$x$')
legend()
xlabel('time (seconds)')
ylabel('value')
title('first population')
subplot(1,2,2)
plot(t, yhat_j, 'g', label='$\hat{y}$')
plot(t, x, 'b', label='$y$')
legend()
xlabel('time (seconds)')
ylabel('value')
title('second population');
In [22]:
#Have to run previous cells first
model.connections.remove(stim_B)
del stim_B
def xhat_fcn(t):
idx = int(t/sim.dt)
if idx>=1000: idx=999
return xhat_i[idx]
with model:
xhat = nengo.Node(xhat_fcn)
nengo.Connection(xhat, ensB)
xhat_p = nengo.Probe(xhat)
sim = nengo.Simulator(model, seed=3)
sim.run(T)
d_j = sim.data[connectionB].weights.T
A_j = sim.data[b_rates]
A_j = A_j + numpy.random.normal(scale=0.2*numpy.max(A_j), size=A_j.shape)
yhat_j = numpy.dot(A_j, d_j)
t = sim.trange()
figure(figsize=(8,4))
subplot(1,2,1)
plot(t, xhat_i, 'g', label='$\hat{x}$')
plot(t, x, 'b', label='$x$')
legend()
xlabel('time (seconds)')
ylabel('value')
title('$\hat{x}$')
ylim(-1,1)
subplot(1,2,2)
plot(t, yhat_j, 'g', label='$\hat{y}$')
plot(t, x, 'b', label='$x$')
legend()
xlabel('time (seconds)')
ylabel('value')
title('$\hat{y}$ (second population)')
ylim(-1,1);
Let's write down what we've done:
Now let's just do the substitutions:
In other words, we can get the entire weight matrix just by multiplying the decoders from the first population with the encoders from the second population
In [23]:
#Have to run previous cells first
n = nengo.neurons.LIFRate()
alpha_j = sim.data[ensB].gain
bias_j = sim.data[ensB].bias
encoders_j = sim.data[ensB].encoders.T
connection_weights = np.outer(alpha_j*encoders_j, d_i)
J_j = np.dot(connection_weights, sim.data[a_rates].T).T + bias_j
A_j = n.rates(J_j, gain=1, bias=0) #Gain and bias already in the previous line
A_j = A_j + numpy.random.normal(scale=0.2*numpy.max(A_j), size=A_j.shape)
xhat_j = numpy.dot(A_j, d_j)
figure(figsize=(8,4))
subplot(1,2,1)
plot(t, xhat_i, 'g', label='$\hat{x}$')
plot(t, x, 'b', label='$x$')
legend()
xlabel('Time (s)')
ylabel('Value')
title('Decode from A')
ylim(-1,1)
subplot(1,2,2)
plot(t, xhat_j, 'g', label='$\hat{y}$')
plot(t, x, 'b', label='$y$')
legend()
xlabel('Time (s)')
title('Decode from B');
ylim(-1,1);
In [7]:
J_j = numpy.outer(numpy.dot(A_i, d_i), alpha_j*encoders_j)+bias_j
In [24]:
import nengo
from nengo.processes import WhiteNoise
from nengo.utils.matplotlib import rasterplot
T = 1.0
max_freq = 5
model = nengo.Network()
with model:
stim = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3))
ensA = nengo.Ensemble(25, dimensions=1)
ensB = nengo.Ensemble(23, dimensions=1)
nengo.Connection(stim, ensA)
nengo.Connection(ensA, ensB, transform=2) #function=lambda x: 2*x)
stim_p = nengo.Probe(stim)
ensA_p = nengo.Probe(ensA, synapse=.01)
ensB_p = nengo.Probe(ensB, synapse=.01)
ensA_spikes_p = nengo.Probe(ensA.neurons, 'spikes')
ensB_spikes_p = nengo.Probe(ensB.neurons, 'spikes')
sim = nengo.Simulator(model, seed=4)
sim.run(T)
t = sim.trange()
figure(figsize=(8, 6))
subplot(2,1,1)
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, sim.data[ensA_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[ensA_spikes_p], ax=ax.twinx(), color='k', use_eventplot=True)
#axis('tight')
ylabel("Neuron")
subplot(2,1,2)
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, sim.data[ensB_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[ensB_spikes_p], ax=ax.twinx(), color='k', use_eventplot=True)
#axis('tight')
ylabel("Neuron");
In [25]:
import nengo
from nengo.processes import WhiteNoise
from nengo.utils.matplotlib import rasterplot
T = 1.0
max_freq = 5
model = nengo.Network()
with model:
stim = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3))
ensA = nengo.Ensemble(25, dimensions=1)
ensB = nengo.Ensemble(23, dimensions=1)
nengo.Connection(stim, ensA)
nengo.Connection(ensA, ensB, function=lambda x: x**2)
stim_p = nengo.Probe(stim)
ensA_p = nengo.Probe(ensA, synapse=.01)
ensB_p = nengo.Probe(ensB, synapse=.01)
ensA_spikes_p = nengo.Probe(ensA.neurons, 'spikes')
ensB_spikes_p = nengo.Probe(ensB.neurons, 'spikes')
sim = nengo.Simulator(model, seed=4)
sim.run(T)
t = sim.trange()
figure(figsize=(8, 6))
subplot(2,1,1)
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, sim.data[ensA_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[ensA_spikes_p], ax=ax.twinx(), color='k', use_eventplot=True)
#axis('tight')
ylabel("Neuron")
subplot(2,1,2)
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, sim.data[ensB_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[ensB_spikes_p], ax=ax.twinx(), color='k', use_eventplot=True)
#axis('tight')
ylabel("Neuron");
In equations:
In code:
In [ ]:
f_x = my_function(x)
gamma=np.dot(A.T,A)
upsilon_f=np.dot(A.T,f_x)
d_f = np.dot(np.linalg.pinv(gamma),upsilon)
xhat = np.dot(A, d_f)
In [8]:
import nengo
from nengo.processes import WhiteNoise
from nengo.utils.matplotlib import rasterplot
T = 1.0
max_freq = 5
model = nengo.Network()
with model:
stimA = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3, seed=3))
stimB = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3, seed=5))
ensA = nengo.Ensemble(25, dimensions=1)
ensB = nengo.Ensemble(23, dimensions=1)
ensC = nengo.Ensemble(24, dimensions=1)
nengo.Connection(stimA, ensA)
nengo.Connection(stimB, ensB)
nengo.Connection(ensA, ensC)
nengo.Connection(ensB, ensC)
stimA_p = nengo.Probe(stimA)
stimB_p = nengo.Probe(stimB)
ensA_p = nengo.Probe(ensA, synapse=.01)
ensB_p = nengo.Probe(ensB, synapse=.01)
ensC_p = nengo.Probe(ensC, synapse=.01)
sim = nengo.Simulator(model)
sim.run(T)
figure(figsize=(8,6))
plot(t, sim.data[stimA_p],'g', label="$x$")
plot(t, sim.data[ensA_p],'b', label="$\hat{x}$")
plot(t, sim.data[stimB_p],'c', label="$y$")
plot(t, sim.data[ensB_p],'m--', label="$\hat{y}$")
plot(t, sim.data[stimB_p]+sim.data[stimA_p],'r', label="$x+y$")
plot(t, sim.data[ensC_p],'k--', label="$\hat{z}$")
legend(loc='best')
ylabel("Output")
xlabel("Time");
In [11]:
import nengo
from nengo.processes import WhiteNoise
from nengo.utils.matplotlib import rasterplot
T = 1.0
max_freq = 5
model = nengo.Network()
with model:
stimA = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3, seed=3), size_out=2)
stimB = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.3, seed=5), size_out=2)
#stimA = nengo.Node([.3,.5])
#stimB = nengo.Node([.3,-.5])
ensA = nengo.Ensemble(55, dimensions=2)
ensB = nengo.Ensemble(53, dimensions=2)
ensC = nengo.Ensemble(54, dimensions=2)
nengo.Connection(stimA, ensA)
nengo.Connection(stimB, ensB)
nengo.Connection(ensA, ensC)
nengo.Connection(ensB, ensC)
stimA_p = nengo.Probe(stimA)
stimB_p = nengo.Probe(stimB)
ensA_p = nengo.Probe(ensA, synapse=.02)
ensB_p = nengo.Probe(ensB, synapse=.02)
ensC_p = nengo.Probe(ensC, synapse=.02)
sim = nengo.Simulator(model)
sim.run(T)
figure()
plot(sim.data[ensA_p][:,0], sim.data[ensA_p][:,1], 'g', label="$\hat{x}$")
plot(sim.data[ensB_p][:,0], sim.data[ensB_p][:,1], 'm', label="$\hat{y}$")
plot(sim.data[ensC_p][:,0], sim.data[ensC_p][:,1], 'k', label="$\hat{z}$")
legend(loc='best')
figure()
plot(t, sim.data[stimA_p],'g', label="$x$")
plot(t, sim.data[ensA_p],'b', label="$\hat{x}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p],'c', label="$y$")
plot(t, sim.data[ensB_p],'m--', label="$\hat{y}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p]+sim.data[stimA_p],'r', label="$x+y$")
plot(t, sim.data[ensC_p],'k--', label="$\hat{z}$")
legend(loc='best')
ylabel("Output")
xlabel("Time");
If these assumptions don't hold, you have to do some other form of optimization
If you already have a decoder for $x$, you can quickly find a decoder for any linear function of $x$
Volunteer for:
$z = Rx$
\cos \theta & \sin \theta \\
-\sin \theta & \cos \theta
\end{array} \right]$$
$z = x \times y$
In [12]:
import nengo
from nengo.processes import WhiteNoise
from nengo.utils.matplotlib import rasterplot
T = 1.0
max_freq = 5
model = nengo.Network()
with model:
stimA = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.5, seed=3))
stimB = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.5, seed=5))
ensA = nengo.Ensemble(55, dimensions=1)
ensB = nengo.Ensemble(53, dimensions=1)
ensC = nengo.Ensemble(200, dimensions=2)
ensD = nengo.Ensemble(54, dimensions=1)
nengo.Connection(stimA, ensA)
nengo.Connection(stimB, ensB)
nengo.Connection(ensA, ensC, transform=[[1],[0]])
nengo.Connection(ensB, ensC, transform=[[0],[1]])
nengo.Connection(ensC, ensD, function=lambda x: x[0]*x[1])
stimA_p = nengo.Probe(stimA)
stimB_p = nengo.Probe(stimB)
ensA_p = nengo.Probe(ensA, synapse=.01)
ensB_p = nengo.Probe(ensB, synapse=.01)
ensC_p = nengo.Probe(ensC, synapse=.01)
ensD_p = nengo.Probe(ensD, synapse=.01)
sim = nengo.Simulator(model)
sim.run(T)
figure()
plot(t, sim.data[stimA_p],'g', label="$x$")
plot(t, sim.data[ensA_p],'b', label="$\hat{x}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p],'c', label="$y$")
plot(t, sim.data[ensB_p],'m--', label="$\hat{y}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p]*sim.data[stimA_p],'r', label="$x * y$")
plot(t, sim.data[ensD_p],'k--', label="$\hat{z}$")
legend(loc='best')
ylabel("Output")
xlabel("Time");
Multiplication is quite powerful, and has lots of uses
Here's a simple gating example using the same network
In [13]:
with model:
stimB.output = lambda t: 0 if (t<.5) else .5
sim = nengo.Simulator(model)
sim.run(T)
figure()
plot(t, sim.data[stimA_p],'g', label="$x$")
plot(t, sim.data[ensA_p],'b', label="$\hat{x}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p],'c', label="$y$")
plot(t, sim.data[ensB_p],'m--', label="$\hat{y}$")
legend(loc='best')
ylabel("Output")
xlabel("Time")
figure()
plot(t, sim.data[stimB_p]*sim.data[stimA_p],'r', label="$x+y$")
plot(t, sim.data[ensD_p],'k--', label="$\hat{z}$")
legend(loc='best')
ylabel("Output")
xlabel("Time");
In [ ]: