Can't make a relaxation oscillator

I’m trying to create a relaxation oscillator based upon the following two equations:

When the input I > 0, the equations above create an oscillation between -2 to 2 for x and 0 to 4 (or more depends on the input) for y.
When the input I<0, there is no oscillation.
I can construct the oscillator using MATLAB, but not with Nengo.

The problem with Nengo is that the ensemble does not always achieve the oscillation successfully.
Sometimes, y stucks at the maximum radius, and sometimes the oscillation traps at a random point.

here is my code:

import nengo
import numpy as np

I =.5
epsilon=.2
gamma=6.0
beta=0.1

tau=.6
syn = 0.05

model = nengo.Network(label='Relaxation Oscillator')
with model:
    inp = nengo.Node(I)
    oscillator = nengo.Ensemble(n_neurons= 150,
                                dimensions=2,
                                radius=4,
                                neuron_type= nengo.LIFRate(tau_rc=0.02, tau_ref=0.002))

     # osc to osc connection
    def feedback(x):
        x,y = x
        dx =  3 * x - x**3 + 2 - y
        dy = epsilon * (gamma * (1 + np.tanh(x / beta)) - y)
        return [tau*dx+x,tau*dy+y]
    nengo.Connection(oscillator, oscillator, function=feedback, synapse=syn)
    # inp to osc connection
    nengo.Connection(inp,oscillator[0], transform = tau)

Where did I do wrong?

Some things I would try:

  • tau and syn should be the same
  • Add synapse=syn to the input connection
  • Try a smaller timestep (dt)
  • Change neuron_type to nengo.Direct() to verify that behaviour is ideal (if not, the network is still not implementing Principle 3 correctly)
1 Like

Thanks for your suggestions arvoelke.

I have made tau and syn be the same and added synapse = syn to the input connection, but unfortunately they don’t solve the problem. I know that you are probably correct about tau and syn should be the same, but I think different value of tau and syn would only make the oscillation faster or slower?

I am not sure how to try a smaller timestep, could you tell me a little bit more about it?

After using nengo.Direct() instead of nengo.LIFRate(), I was able to get the oscillation every time I run the code :grin:. I guess it means that I do have the correct code, but if I want to use LIFRate, I will have to play around its properties…?

1 Like

This is the “magic” of Principle 3. The tau that you provide on the transform(s) actually cancel out with the tau on the synapse, and so the speed will not be changed. Setting them to be the same means that the dynamics that you get are the dynamics that you want.

I’m not sure if the nengo GUI has an option to control the simulation time-step. You can run your simulations as a separate script, and use with nengo.Simulator(model, dt=dt) as sim: sim.run(T) to run the simulation for T seconds with a time-step of dt. This is what I normally do and then use matplotlib / pylab / seaborn to visualize plots, but this requires a bit of extra coding. However… since Direct mode works the dt should not be an issue (see below).

There should not be a large difference between Direct and LIFRate unless your ensemble is incapable of properly representing the 2D space or accurately approximating the required function. This might mean you need to increase the number of neurons and/or the radius.

Thanks arvoelke, I just corrected my code as you suggested, but I am still not getting the oscillation.

Here is what I did. I used Direct and gradually increased the radius of my ensemble to find the radius I should be using. I found the maximum of y can reach to 6.5, so I changed my radius to 8 (just to be safe).

Then I switched to LIFRate and gradually increased the number of neurons in my ensemble. I increased to 1000 but the oscillation still can fail sometimes.

Is there anything else I could try? :sob:

(I hope this doesn’t mean the equations are too complicated for nengo…)

Could you please re-post the newest version of your code that fails sometimes? I will take a closer look.

Thanks.

# <editor-fold desc="...imports">
import nengo
import numpy as np
import matplotlib.pyplot as plt
# </editor-fold>

# <editor-fold desc="...constants">
I =.5
epsilon=.2
gamma=6.0
beta=0.1
syn = 9
# </editor-fold>
net = nengo.Network(label='Relaxation Oscillator')
with net:
    inp = nengo.Node(I)
    oscillator = nengo.Ensemble(n_neurons= 500,
                                dimensions=2,
                                radius=8,
                                neuron_type= nengo.LIFRate())

     # osc to osc connection
    def feedback(x):
        x,y = x
        dx =  3 * x - x**3 + 2 - y
        dy = epsilon * (gamma * (1 + np.tanh(x / beta)) - y)
        return [syn*dx+x,syn*dy+y]
    nengo.Connection(oscillator, oscillator, function=feedback, synapse=syn)
    # inp to osc connection
    nengo.Connection(inp,oscillator[0], transform = syn, synapse=syn)
    x_pr = nengo.Probe(oscillator[0], synapse=0.01)
    y_pr = nengo.Probe(oscillator[1], synapse=0.01)
with nengo.Simulator(net) as sim:
    sim.run(40)
# <editor-fold desc="...plot">
t = sim.trange()
# xy activities
fig1 = plt.figure(figsize=(10, 5))
ax1 = fig1.add_subplot(1, 2, 1)
ax1.plot(t, sim.data[x_pr], label="x activity", color='b')
ax1.plot(t, sim.data[y_pr], label="y activity", color='r')
ax1.set_title("Activities of the Oscillator")
ax1.set_xlabel("Time (s)")
ax1.set_ylabel("x&y activities")
ax1.set_ylim(-3, 7)
ax1.legend()

# phase plane
xmin, xmax, ymin, ymax = -2.5, 2.5, -2, 8
ax2 = fig1.add_subplot(1, 2, 2)
ax2.plot(sim.data[x_pr], sim.data[y_pr], label="Neuron Activity", color='#ffa500', linewidth=.5, marker='x')
X = np.linspace(xmin, xmax)
ax2.plot(X, 3. * X - X ** 3 + 2 + I, label='x-nullcline', color='b')
ax2.plot(X, gamma * (1 + np.tanh(X / beta)), label='y-nullcline', color='r')
ax2.set_title("Phase Plane")
ax2.set_xlabel("x activity")
ax2.set_ylabel("y activity")
ax2.set_ylim(ymin, ymax)
ax2.set_xlim(xmin, xmax)
ax2.legend()
plt.show()
# </editor-fold>

Hmm so this is quite the nice challenge!

The hard part of the problem comes from the fact that the feedback needs to compute the term np.tanh(x / beta). If we plot this, we see that it has a fairly sharp transition from -1 to +1 around 0

By default, Nengo distributes the tuning curves uniformly about the space, and so they are not optimized to accurately compute such functions. With the defaults, you would likely need tens of thousands of neurons to approach the level of accuracy that you need to get the dynamics that you want.

Fortunately, there is a way. The most straight-forward approach is to place all of the tuning curves for the x-coordinate so that they “turn on” at x=0. This point is known as the intercept of each neuron.

from nengo.utils.ensemble import tuning_curves

with nengo.Network() as model:
    x = nengo.Ensemble(
        n_neurons=50, dimensions=1, intercepts=np.zeros(50))

with nengo.Simulator(model) as sim:
    plt.plot(*tuning_curves(x, sim))

plt.ylabel("Firing rate (Hz)")
plt.xlabel("Input scalar, x")
plt.show()

As seen above, this provides tuning curves that more closely match the overall shape of the desired function. Consequently, this ensemble may compute the function with much greater accuracy with fewer numbers of neurons.

To make this work, we now need two separate ensembles, one for the x-coordinate and one for the y-coordinate. We can separate them out because there are no nonlinear interactions between coordinates. We now need 4 separate connections, such that their addition computes the same feedback function as before.

Here is my revised code:

import nengo
import numpy as np
import matplotlib.pyplot as plt

I =.5
epsilon=.2
gamma=6.0
beta=0.1
syn=0.1

n_neurons = 500
neuron_type = nengo.LIFRate()
syn_probe = 0.01

net = nengo.Network(label='Relaxation Oscillator')
with net:
    inp = nengo.Node(I)

    oscillator_x = nengo.Ensemble(
        n_neurons=n_neurons, dimensions=1, neuron_type=neuron_type,
        radius=2, intercepts=np.zeros(n_neurons))

    oscillator_y = nengo.Ensemble(
        n_neurons=n_neurons, dimensions=1, neuron_type=neuron_type,
        radius=7)

    '''
    # osc to osc connection
    def feedback(x):
        x,y = x
        dx =  3 * x - x**3 + 2 - y
        dy = epsilon * (gamma * (1 + np.tanh(x / beta)) - y)
        return [syn*dx+x,syn*dy+y]
    '''

    nengo.Connection(oscillator_x, oscillator_x, synapse=syn,
                     function=lambda x: syn*(3*x - x**3 + 2) + x)
    nengo.Connection(oscillator_x, oscillator_y, synapse=syn,
                     function=lambda x: syn*(epsilon * (gamma * (1 + np.tanh(x / beta)))))
    nengo.Connection(oscillator_y, oscillator_x, synapse=syn,
                     function=lambda y: syn*(-y))
    nengo.Connection(oscillator_y, oscillator_y, synapse=syn,
                     function=lambda y: syn*(epsilon * (-y)) + y)

    # inp to osc connection
    nengo.Connection(inp, oscillator_x, transform=syn, synapse=syn)
    x_pr = nengo.Probe(oscillator_x, synapse=syn_probe)
    y_pr = nengo.Probe(oscillator_y, synapse=syn_probe)

with nengo.Simulator(net) as sim:
    sim.run(40)

t = sim.trange()
# xy activities
fig1 = plt.figure(figsize=(10, 5))
ax1 = fig1.add_subplot(1, 2, 1)
ax1.plot(t, sim.data[x_pr], label="x activity", color='b')
ax1.plot(t, sim.data[y_pr], label="y activity", color='r')
ax1.set_title("Activities of the Oscillator")
ax1.set_xlabel("Time (s)")
ax1.set_ylabel("x&y activities")
ax1.set_ylim(-3, 7)
ax1.legend()

# phase plane
xmin, xmax, ymin, ymax = -2.5, 2.5, -2, 8
ax2 = fig1.add_subplot(1, 2, 2)
ax2.plot(sim.data[x_pr], sim.data[y_pr], label="Neuron Activity", color='#ffa500', linewidth=.5, marker='x')
X = np.linspace(xmin, xmax)
ax2.plot(X, 3. * X - X ** 3 + 2 + I, label='x-nullcline', color='b')
ax2.plot(X, gamma * (1 + np.tanh(X / beta)), label='y-nullcline', color='r')
ax2.set_title("Phase Plane")
ax2.set_xlabel("x activity")
ax2.set_ylabel("y activity")
ax2.set_ylim(ymin, ymax)
ax2.set_xlim(xmin, xmax)
ax2.legend()
plt.show()

Some further optimizations should be possible. For instance, note that we only need the positive portion of the y space, and so half the space is currently wasted. This could be solved by shifting the space to be centered around 0 and cutting the radius in half.

There may also be a better choice or mix of tuning curves for the x ensemble that is better optimized to accurately compute both np.tanh(x / beta))) and 3*x - x**3 + 2 simultaneously. Solving this sort of optimization problem in general is very difficult (it is a nonlinear optimization problem with discontinuous gradients), and so Nengo currently leaves this as a problem for the user.

1 Like