NEF: Dynamics

Chapter 12.1.3

Learn about the third principle of NEF: dynamics

Read Chapter 12.1.3

Python / Nengo demonstration

Imports:

import nengo
import matplotlib.pyplot as plt
from nengo.processes import Piecewise
from nengo.utils.ensemble import tuning_curves
from nengo.utils.ensemble import sorted_neurons

Recurrent computing of f(x) = x + 1

model = nengo.Network()

with model:
    ensA = nengo.Ensemble(100, dimensions=1)
    ensB = nengo.Ensemble(100, dimensions=1)
    ensC = nengo.Ensemble(100, dimensions=1)
    
    def feedback(x):
        return x+1
    
    nengo.Connection(ensA, ensA, function=feedback, synapse = 0.1)
    nengo.Connection(ensB, ensB, function=feedback, synapse = 0.2)
    nengo.Connection(ensC, ensC, function=feedback, synapse = 0.3)

    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(.5)

plt.plot(sim.trange(), sim.data[ensA_p], label="$tau=0.1$")
plt.plot(sim.trange(), sim.data[ensB_p], label='$tau=0.2$')
plt.plot(sim.trange(), sim.data[ensC_p], label='$tau=0.3$')
plt.legend()
plt.ylabel("Output")
plt.xlabel("Time")
plt.ylim(0,1.5);

Result:

Recurrent computing of f(x) = -x

with model:
    
    stim = nengo.Node(Piecewise({0:1, .2:-1, .4:0}))
    
    def feedback(x):
        return -x
    
    ensA = nengo.Ensemble(100, dimensions=1)
    ensB = nengo.Ensemble(100, dimensions=1)
    ensC = nengo.Ensemble(100, dimensions=1)
    
    nengo.Connection(stim, ensA)
    nengo.Connection(stim, ensB)
    nengo.Connection(stim, ensC)
    
    nengo.Connection(ensA, ensA, function=feedback, synapse = 0.1)
    nengo.Connection(ensB, ensB, function=feedback, synapse = 0.2)
    nengo.Connection(ensC, ensC, function=feedback, synapse = 0.3)
    
    ensA_p = nengo.Probe(ensA, synapse=.01)
    ensB_p = nengo.Probe(ensB, synapse=.01)
    ensC_p = nengo.Probe(ensC, synapse=.01)
    stim_p = nengo.Probe(stim, synapse=.01)
    
    
sim = nengo.Simulator(model)
sim.run(.6)

plt.plot(sim.trange(), sim.data[stim_p], 'r', label='stimulus')
plt.plot(sim.trange(), sim.data[ensA_p], label="$tau=0.01$")
plt.plot(sim.trange(), sim.data[ensB_p], label='$tau=0.02$')
plt.plot(sim.trange(), sim.data[ensC_p], label='$tau=0.03$')
plt.legend()
plt.ylabel("Output")
plt.xlabel("Time")

Result:

Recurrent computing of f(x) = x^2

with model:
    
    stim = nengo.Node(Piecewise({.1:.2, .2:.4, 0.6:0}))
    
    def feedback(x):
        return x*x
    
    ensA = nengo.Ensemble(100, dimensions=1)
    ensB = nengo.Ensemble(100, dimensions=1)
    ensC = nengo.Ensemble(100, dimensions=1)
    
    nengo.Connection(stim, ensA)
    nengo.Connection(stim, ensB)
    nengo.Connection(stim, ensC)
    
    nengo.Connection(ensA, ensA, function=feedback, synapse = 0.1)
    nengo.Connection(ensB, ensB, function=feedback, synapse = 0.2)
    nengo.Connection(ensC, ensC, function=feedback, synapse = 0.3)
    
    ensA_p = nengo.Probe(ensA, synapse=.01)
    ensB_p = nengo.Probe(ensB, synapse=.01)
    ensC_p = nengo.Probe(ensC, synapse=.01)
    stim_p = nengo.Probe(stim, synapse=.01)
    
    
sim = nengo.Simulator(model)
sim.run(1)

plt.plot(sim.trange(), sim.data[stim_p], 'r', label='stimulus')
plt.plot(sim.trange(), sim.data[ensA_p], label="$tau=0.01$")
plt.plot(sim.trange(), sim.data[ensB_p], label='$tau=0.02$')
plt.plot(sim.trange(), sim.data[ensC_p], label='$tau=0.03$')
plt.legend()
plt.ylabel("Output")
plt.xlabel("Time")

Result:

Integrator

import matplotlib
font = {'family' : 'normal',
        'weight' : 'bold',
        'size'   : 12}

matplotlib.rc('font', **font)

tau = 0.001
tau2 = 0.01
tau3 = 1

model = nengo.Network('Eye control', seed=8)

with model:
    stim = nengo.Node(Piecewise({.3:1, .6:0 }))
    velocity = nengo.Ensemble(100, dimensions=1)
    position = nengo.Ensemble(200, dimensions=1)
    position2 = nengo.Ensemble(200, dimensions=1)
    position3 = nengo.Ensemble(200, dimensions=1)
    
    def feedback(x):
        return 1*x
    
    conn = nengo.Connection(stim, velocity)
    conn = nengo.Connection(velocity, position, transform=tau, synapse=tau)
    conn = nengo.Connection(position, position, function=feedback, synapse=tau)
    spikes_p = nengo.Probe(velocity.neurons, 'spikes')
    
    conn2 = nengo.Connection(velocity, position2, transform=tau2, synapse=tau2)
    conn2 = nengo.Connection(position2, position2, function=feedback, synapse=tau2)
    
    conn3 = nengo.Connection(velocity, position3, transform=tau3, synapse=tau3)
    conn3 = nengo.Connection(position2, position3, function=feedback, synapse=tau3)

    stim_p = nengo.Probe(stim)
    velocity_p = nengo.Probe(velocity, synapse=.01)
    position_p = nengo.Probe(position, synapse=.01)
    position_p2 = nengo.Probe(position2, synapse=.01)
    position_p3 = nengo.Probe(position3, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(1)

plt.figure()
plt.plot(sim.trange(), sim.data[stim_p], label = "Input", linewidth=4, color='black')
plt.plot(sim.trange(), sim.data[velocity_p], label = "velocity", linewidth=2)
plt.plot(sim.trange(), sim.data[position_p], label = "Position (tau=0.001)", linewidth=2)
plt.plot(sim.trange(), sim.data[position_p2], label = "position (tau=0.1)", linewidth=2)
plt.plot(sim.trange(), sim.data[position_p3], label = "position (tau=1)", linewidth=2)
plt.ylabel("Output")
plt.xlabel("Time")

Result:

Leaky integrator

tau = 0.1
tau_c = 2.0

model = nengo.Network('Eye control', seed=5)

with model:
    stim = nengo.Node(Piecewise({.3:1, .6:0 }))
    velocity = nengo.Ensemble(100, dimensions=1)
    position = nengo.Ensemble(200, dimensions=1)
    
    def feedback(x):
        return (-tau/tau_c + 1)*x
    
    conn = nengo.Connection(stim, velocity)
    conn = nengo.Connection(velocity, position, transform=tau, synapse=tau)
    conn = nengo.Connection(position, position, function=feedback, synapse=tau)

    stim_p = nengo.Probe(stim)
    position_p = nengo.Probe(position, synapse=.01)
    velocity_p = nengo.Probe(velocity, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(3)

plt.plot(sim.trange(), sim.data[stim_p], label = "stim")
plt.plot(sim.trange(), sim.data[position_p], label = "position")
plt.plot(sim.trange(), sim.data[velocity_p], label = "velocity")
plt.ylabel("Output")
plt.xlabel("Time")
plt.legend(loc="best");

Result:

Controlled leaky integrator

tau = 0.1

model = nengo.Network('Controlled integrator', seed=1)

with model:
    vel = nengo.Node(Piecewise({.2:1.5, .5:0 }))
    dec = nengo.Node(Piecewise({.7:.2, .9:0 }))
    
    velocity = nengo.Ensemble(100, dimensions=1)
    decay = nengo.Ensemble(100, dimensions=1)
    position = nengo.Ensemble(400, dimensions=2)
    
    def feedback(x):
        return -x[1]*x[0]+x[0], 0
    
    conn = nengo.Connection(vel, velocity)
    conn = nengo.Connection(dec, decay)
    conn = nengo.Connection(velocity, position[0], transform=tau, synapse=tau)
    conn = nengo.Connection(decay, position[1], synapse=0.01)
    conn = nengo.Connection(position, position, function=feedback, synapse=tau)

    position_p = nengo.Probe(position[0], synapse=.01)
    velocity_p = nengo.Probe(velocity, synapse=.01)
    decay_p = nengo.Probe(decay, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(1)

plt.plot(sim.trange(), sim.data[decay_p])
plt.lineObjects = plt.plot(sim.trange(), sim.data[position_p])
plt.plot(sim.trange(), sim.data[velocity_p])
plt.ylabel("Output")
plt.xlabel("Time")
plt.legend(('decay','position','velocity'),loc="best");

Result:

Oscillator

model = nengo.Network('Oscillator')

freq = -0.25

with model:
    stim = nengo.Node(lambda t: [.5,.5] if t<.02 else [0,0])
    osc = nengo.Ensemble(200, dimensions=2)
    
    def feedback(x):
        return x[0]+freq*x[1], -freq*x[0]+x[1]
    
    nengo.Connection(osc, osc, function=feedback, synapse=.01)
    nengo.Connection(stim, osc)
    
    stim_p = nengo.Probe(stim)
    osc_p = nengo.Probe(osc, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(.5)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(sim.trange(), sim.data[osc_p]);
plt.plot(sim.trange(), sim.data[stim_p], 'r', label = "stim", linewidth=4 )
plt.xlabel('Time (s)')
plt.ylabel('State value')       
plt.subplot(1,2,2)
plt.plot(sim.data[osc_p][:,0],sim.data[osc_p][:,1])
plt.xlabel('$x_0$')
plt.ylabel('$x_1$');

Result:

Cnotrolled oscillator

model = nengo.Network('Oscillator')

freq = -0.25

with model:
    stim = nengo.Node(lambda t: [.5,.5] if t<.02 else [0,0])
    freq_ctrl = nengo.Node(Piecewise({0:-0.1, 4:-.2, 8:-0.3}))
    osc = nengo.Ensemble(200, dimensions=3)
    
    def feedback(x):
        return x[0]+x[2]*x[1], -x[2]*x[0]+x[1], 0
    
    nengo.Connection(osc, osc, function=feedback, synapse=.01)
    nengo.Connection(stim, osc[0:2])
    nengo.Connection(freq_ctrl, osc[2])
    
    stim_p = nengo.Probe(stim)
    osc_p = nengo.Probe(osc, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(12)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(sim.trange(), sim.data[osc_p]);
plt.plot(sim.trange(), sim.data[stim_p], 'r', label = "stim", linewidth=4 )
plt.xlabel('Time (s)')
plt.ylabel('State value')       
plt.subplot(1,2,2)
plt.plot(sim.data[osc_p][:,0],sim.data[osc_p][:,1])
plt.xlabel('$x_0$')
plt.ylabel('$x_1$');

Result:

Point attractor

model = nengo.Network('Oscillator')

freq = -0.25

with model:
    
    stim = nengo.Node(lambda t: [.5,.5] if t<.01 else [0,0])
    osc = nengo.Ensemble(2000, dimensions=2)
    
    def feedback(x):
        p1 = 0.5
        p2 = 0.8
        return x[0]-(x[0]-p1), x[1] - (x[1]-p2)
    
    nengo.Connection(osc, osc, function=feedback, synapse=.01)
    nengo.Connection(stim, osc)
    
    stim_p = nengo.Probe(stim)
    osc_p = nengo.Probe(osc, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(.5)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(sim.trange(), sim.data[osc_p]);
plt.plot(sim.trange(), sim.data[stim_p], 'r', label = "stim", linewidth=4 )
plt.xlabel('Time (s)')
plt.ylabel('State value')       
plt.subplot(1,2,2)
plt.plot(sim.data[osc_p][:,0],sim.data[osc_p][:,1])
plt.xlabel('$x_0$')
plt.ylabel('$x_1$');
plt.xlim(0,1)
plt.ylim(0,1)

Result:

2D plane attractor (2D ensemble)

model = nengo.Network(label='2D Plane Attractor', seed=4)

N = 200 #600 
tau = 0.01

with model:
    stim = nengo.Node(Piecewise({.3:[1, 0], .5:[0, 0], .7:[0, -1], .9:[0,0]}))
    neurons = nengo.Ensemble(N, dimensions=2)
    neurons2 = nengo.Ensemble(N*10, dimensions=2)

    nengo.Connection(stim, neurons, transform=tau, synapse=tau)
    nengo.Connection(neurons, neurons, synapse=tau)
    nengo.Connection(stim, neurons2, transform=tau, synapse=tau)
    nengo.Connection(neurons2, neurons2, synapse=tau)

    stim_p = nengo.Probe(stim)
    neurons_p = nengo.Probe(neurons, synapse=.01)
    neurons_p2 = nengo.Probe(neurons2, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(4)

t=sim.trange()

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(t, sim.data[stim_p], label = "stim")
plt.plot(t, sim.data[neurons_p], label = "position")
plt.ylabel("Output")
plt.xlabel("Time")
plt.legend(loc="best")
plt.title('# neurons = 200')

plt.subplot(1,2,2)
plt.plot(t, sim.data[stim_p], label = "stim")
plt.plot(t, sim.data[neurons_p2], label = "position")
plt.ylabel("Output")
plt.xlabel("Time")
plt.legend(loc="best")
plt.title('# neurons = 2000')

Results:

2D plane attractor (ensemble array)

N = 500 #neurons per sub_ensemble
tau = 0.01

with model:
    stim = nengo.Node(Piecewise({.5:[1, 0], 1:[0, 0], 2:[0, -1], 2.5:[0,0]}))

    neurons = nengo.networks.EnsembleArray(N, n_ensembles=2, seed=6)
    
    nengo.Connection(stim, neurons.input, transform=tau, synapse=tau)
    nengo.Connection(neurons.output, neurons.input, synapse=tau)

    stim_p = nengo.Probe(stim)
    neurons_p = nengo.Probe(neurons.output, synapse=.01)
    
sim = nengo.Simulator(model)
sim.run(4)

t=sim.trange()

plt.plot(t, sim.data[stim_p], label = "stim")
plt.plot(t, sim.data[neurons_p], label = "position");
plt.ylabel("Output")
plt.xlabel("Time")
plt.legend(loc="best")

Result:

Lorenz attractor

model = nengo.Network('Lorenz Attractor', seed=3)

with model:
    
    x = nengo.Ensemble(n_neurons=600, dimensions=3, radius=30)

    synapse = 0.1
    def lorenz(x):
        sigma = 10
        beta = 8.0/3
        rho = 28

        dx0 = -sigma * x[0] + sigma * x[1]
        dx1 = -x[0] * x[2] - x[1]
        dx2 = x[0] * x[1] - beta * (x[2] + rho) - rho

        return [dx0 * synapse + x[0],
                dx1 * synapse + x[1],
                dx2 * synapse + x[2]]

    nengo.Connection(x, x, synapse=synapse, function=lorenz)

    lorenz_p = nengo.Probe(x, synapse=tau)
    
sim = nengo.Simulator(model)
sim.run(14)

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.plot(sim.trange(), sim.data[lorenz_p][:,0], label='$x_0$')
plt.plot(sim.trange(), sim.data[lorenz_p][:,1], label='$x_1$')
plt.plot(sim.trange(), sim.data[lorenz_p][:,2], label='$x_2$')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('State value')
plt.subplot(1,2,2)
plt.plot(sim.data[lorenz_p][:,0],sim.data[lorenz_p][:,1])
plt.xlabel('$x_0$') 
plt.ylabel('$x_1$');

Result: