Copyright CRC Press 2023 - present CRC Press
Author: Dr Stephen Lynch National Teaching Fellow FIMA SFHEA
Copyright CRC Press 2023 - present CRC Press
Author: Dr Stephen Lynch National Teaching Fellow FIMA SFHEA
# 16.1: Bifurcation of a Fitzhugh-Nagumo Limit Cycle.
# Hopf (safe) bifurcations.
from matplotlib import pyplot as plt
from matplotlib.animation import ArtistAnimation
import numpy as np
from scipy.integrate import odeint
fig = plt.figure()
myimages = []
theta , epsilon , gamma = 0.14 , 0.01 , 2.54
tmax = 200
def FN_ode(x, t):
return [i + x[0] * (x[0] - theta) * (1 - x[0]) - x[1], \
epsilon * (x[0] - gamma * x[1])]
time = np.arange(0, tmax, 0.1)
x0 = [0, 0.1]
for i in np.arange(0, 0.2, 0.01):
xs = odeint(FN_ode, x0, time)
imgplot = plt.plot(xs[:,0], xs[:, 1], 'r-')
myimages.append(imgplot)
anim = ArtistAnimation(fig, myimages, interval=500, blit=False, repeat_delay=100)
plt.close()
plt.show()
from IPython.display import HTML
HTML(anim.to_jshtml())
# 16.2: Morris-Lecar Model of a neuron, three limit cycles.
# The black limit cycle is unstable, the red and blue are stable.
# Can you produce an animation of a bifurcation from one limit cycle to another?
from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import odeint
fig = plt.figure()
i,gfast,V1,V2,V3,V4=82,20,-1.2,18,-20.5,10
ENa,EK,gslow,gleak,Eleak,phi,C=50,-100,20,2,-70,0.15,2
tmax = 100
def ML_ode(x, t):
return [(i-gfast*0.5*(1+np.tanh((x[0]-V1)/V2))*(x[0]-ENa)-gslow*x[1]*(x[0]-EK)-gleak*(x[0]-Eleak))/C, \
phi*(0.5*(1+np.tanh((x[0]-V3)/V4))-x[1])*np.cosh((x[0]-V3)/(2*V4))]
time = np.linspace(0, tmax, 1000)
x0=[-40,0.335]
xs1 = odeint(ML_ode, x0, time)
imgplot = plt.plot(xs1[:,0], xs1[:, 1], 'r-')
x0=[-40,0.175]
xs2 = odeint(ML_ode, x0, -time)
imgplot = plt.plot(xs2[:,0], xs2[:, 1], 'k-')
x0=[-40,0.03]
xs3 = odeint(ML_ode, x0, time)
imgplot = plt.plot(xs3[:,0], xs3[:, 1], 'b-')
plt.xlabel("V", fontsize=15)
plt.ylabel("N", fontsize=15)
plt.tick_params(labelsize=15)
plt.show()
# 16.3: Fitzhugh-Nagumo full adder
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Input current.
def Input_1(t):
return 1*(t>500)-1*(t>1000)+1*(t>2000)-1*(t>2500)+1*(t>3000)
def Input_2(t):
return 1*(t>1000)-1*(t>1500)+1*(t>2500)
def Input_3(t):
return 1*(t>1500)-1*(t>3000)+1*(t>3500)
# Constants
theta , gamma , epsilon , Tmax =0.1 , 0.1 , 0.1 , 2000
m , c =-100 , 60
w1 , w2 , x1 = 0.8 , 0.45 , 1.5 #weights
t=np.arange(0.0, 4000.0, 0.1)
def FN_ODEs(X,t):
u1,v1,u2,v2,u3,v3,usum,vsum,ucarry,vcarry=X
du1 = -u1*(u1-theta)*(u1-1)-v1+Input_1(t) #I_1
dv1 = epsilon*(u1-gamma*v1)
du2 = -u2*(u2-theta)*(u2-1)-v2+Input_2(t) #I_2
dv2 = epsilon*(u2-gamma*v2)
du3 = -u3*(u3-theta)*(u3-1)-v3+Input_3(t) #I_3
dv3 = epsilon*(u3-gamma*v3)
dusum = -usum*(usum-theta)*(usum-1)-vsum+w1/(1+np.exp(m*v1+c))+ w1/(1+np.exp(m*v2+c))+w1/(1+np.exp(m*v3+c))-x1/(1+np.exp(m*vcarry+c)) #O_1
dvsum = epsilon*(usum-gamma*vsum)
ducarry = -ucarry*(ucarry-theta)*(ucarry-1)-vcarry+w2/(1+np.exp(m*v1+c))+ w2/(1+np.exp(m*v2+c)) + w2/(1+np.exp(m*v3+c)) #O_2
dvcarry = epsilon*(ucarry-gamma*vcarry)
return du1,dv1,du2,dv2,du3,dv3,dusum,dvsum,ducarry,dvcarry
# Solve the coupled ODEs.
X = odeint(FN_ODEs,[0.01,0.01,0.01,0.01,0.01,0.01,0,0,0,0],t,rtol=1e-6)
plt.subplots_adjust(hspace = 1)
plt.figure(1)
plt.rcParams["font.size"] = "14"
plt.subplot(5,1,1)
plt.title("Fitzhugh-Nagumo Full-Adder")
plt.plot(t, X[: , 0], "b")
plt.ylim(-1, 1.5)
plt.ylabel("I$_1$")
plt.subplot(5,1,2)
plt.plot(t, X[:,2], "b")
plt.ylim(-1, 1.5)
plt.ylabel("I$_2$")
plt.subplot(5,1,3)
plt.plot(t, X[:,4], "b")
plt.ylim(0, 1)
plt.ylim(-1, 1.5)
plt.ylabel("I$_3$")
plt.subplot(5,1,4)
plt.plot(t, X[:,6], "g")
plt.ylim(-1, 1.5)
plt.ylabel("O$_1$")
plt.subplot(5,1,5)
plt.plot(t, X[:,8], "g")
plt.ylim(0, 1)
plt.ylim(-1, 1.5)
plt.ylabel("O$_2$")
plt.xlabel("Time")
plt.show()
# 16.4: The Fitzhugh-Nagumo Ballistic SR Flip-Flop with Noise.
# Program takes 1 minute to compile.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint # Input current.
def Input_1(t):
return 1*(t>1000)-1*(t>1010)
def Input_2(t):
return 1*(t>2000)-1*(t>2010)
def Ic1(t):
return 1*(t>500)
def Ic2(t):
return 1
# Constants
theta , gamma , epsilon , Tmax = 0.1 , 0.1 , 0.1 , 2000
m , c = -100 , 60
t=np.arange(0.0, 3000.0, 0.1)
length = len(t)
gauss = 2.5*np.random.normal(0,0.05,length*4) # Gaussian noise array
# Function to add Gaussian noise to system.
def noise(t,n):
interval = round(10*t-1)*n
return gauss[interval]
w1 , x1 = 0.5 , -1
def FN_ODEs(X,t):
u1,v1,u2,v2,u3,v3,u4,v4=X
du1 = -u1*(u1-theta)*(u1-1)-v1+Input_1(t) +noise(t,1) #I_1
dv1 = epsilon*(u1-gamma*v1)
du2 = -u2*(u2-theta)*(u2-1)-v2+Input_2(t) +noise(t,2) #I_2
dv2 = epsilon*(u2-gamma*v2)
du3 = -u3*(u3-theta)*(u3-1)-v3+w1/(1+np.exp(m*v1+c))+x1/(1+np.exp(m*v4+c)) + Ic1(t) + noise(t,3) #O_1
dv3 = epsilon*(u3-gamma*v3)
du4 = -u4*(u4-theta)*(u4-1)-v4+0.45/(1+np.exp(m*v2+c))+x1/(1+np.exp(m*v3+c)) + Ic2(t) + noise(t,4) #O_2
dv4 = epsilon*(u4-gamma*v4)
return du1,dv1,du2,dv2,du3,dv3,du4,dv4
X = odeint(FN_ODEs,[0.01,0.01,0.01,0.01,0,0,0,0],t,rtol=1e-6)
u1 = X[:,0];v1 = X[:,1];u2 = X[:,2];v2 = X[:,3];
u3 = X[:,4];v3 = X[:,5];u4 = X[:,6];v4 = X[:,7]
plt.subplots_adjust(hspace = 1)
plt.figure(1)
plt.rcParams["font.size"] = "14"
plt.subplot(4,1,1)
plt.title("Fitzhugh-Nagumo SR Flip-Flop With Noise")
plt.plot(t, u1, "b")
plt.ylim(-1, 1.5)
plt.ylabel("I$_1$")
plt.subplot(4,1,2)
plt.plot(t, u2, "b")
plt.ylim(-1, 1.5)
plt.ylabel("I$_2$")
plt.subplot(4,1,3)
plt.plot(t, u3, "g")
plt.ylim(0, 1)
plt.ylim(-1, 1.5)
plt.ylabel("O$_1$")
plt.subplot(4,1,4)
plt.plot(t, u4, "g")
plt.ylim(-1, 1.5)
plt.ylabel("O$_2$")
plt.xlabel("Time")
plt.show()
# 17.1: ANN for an XOR Logic Gate.
import numpy as np
w11 , w12 , w21 , w22 , w13 , w23 = 60 , 80 , 60 , 80 , -60 , 60
b1 , b2 , b3 = -90 , -40 , -30
def sigmoid(x):
return 1 / (1 + np.exp(-x))
def XOR(x1 , x2):
h1 = x1 * w11 + x2 * w21 + b1
h2 = x1 * w12 + x2 * w22 + b2
o1 = sigmoid(h1) * w13 + sigmoid(h2) * w23 + b3
return sigmoid(o1)
print("XOR(0,0)= " , XOR(0 ,0))
print("XOR(0,1)= " , XOR(0 ,1))
print("XOR(1,0)= " , XOR(1 ,0))
print("XOR(1,1)= " , XOR(1 ,1))
XOR(0,0)= 9.357622968839299e-14 XOR(0,1)= 0.9999999999999065 XOR(1,0)= 0.9999999999999065 XOR(1,1)= 9.357622968891759e-14
# 17.2: Backpropagation, keep the biases constant in this case.
import numpy as np
w11,w12,w21,w22,w13,w23 = 0.2,0.15,0.25,0.3,0.15,0.1
b1 , b2 , b3 = -1 , -1 , -1
yt , eta = 0 , 0.1
x1 , x2 = 1 , 1
def sigmoid(v):
return 1 / (1 + np.exp(- v))
h1 = x1 * w11 + x2 * w21 + b1
h2 = x1 * w12 + x2 * w22 + b2
o1 = sigmoid(h1) * w13 + sigmoid(h2) * w23 + b3
y = sigmoid(o1)
print("y = ", y)
# Backpropagate.
dErrdw13=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*sigmoid(h1)
w13 = w13 - eta * dErrdw13
print("w13_new = ", w13)
dErrdw23=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*sigmoid(h2)
w23 = w23 - eta * dErrdw23
print("w23_new = ", w23)
dErrdw11=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*w13*sigmoid(h1)*(1-sigmoid(h1))*x1
w11 = w11 - eta * dErrdw11
print("w11_new = ", w11)
dErrdw12=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*w23*sigmoid(h2)*(1-sigmoid(h2))*x1
w12 = w12 - eta * dErrdw12
print("w12_new = ", w12)
dErrdw21=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*w13*sigmoid(h1)*(1-sigmoid(h1))*x2
w21 = w21 - eta * dErrdw21
print("w21_new = ", w21)
dErrdw22=(yt-y)*sigmoid(o1)*(1-sigmoid(o1))*w23*sigmoid(h1)*(1-sigmoid(h1))*x2
w22 = w22 - eta * dErrdw22
print("w22_new = ", w22)
y = 0.28729994077761756 w13_new = 0.1521522763401024 w23_new = 0.10215227634010242 w11_new = 0.20020766275648338 w12_new = 0.15013942100503588 w21_new = 0.2502076627564834 w22_new = 0.3001394210050359
# 17.3: Boston Housing with a Hidden Layer.
# Run in Google Colab.
# Results do not get any better after three hidden neurons.
import numpy as np
import matplotlib.pyplot as plt
# Import keras so that we can access the Boston housing data.
import tensorflow as tf
from tensorflow import keras
## Parameters
num_epochs = 5000
max_num_hidden = 3
eta = 0.05
## Functions
def UniformRandomMatrix (rows, cols):
res = [[np.random.uniform () for c in range (cols)] for r in range (rows)]
return np.matrix (res)
## Load the data
dataset = keras.datasets.boston_housing
(train_X, train_y), (_,_) = dataset.load_data (test_split = 0)
train_X = np.matrix (train_X)
train_y = np.matrix (train_y)
(num_samples, num_inputs) = train_X.shape
# Add bias
bias = np.ones (num_samples)
bias = np.matrix (bias).transpose ()
train_X = np.append (train_X, bias, axis=1)
## Normalize data
for i in range (num_inputs):
col = train_X[:,i]
train_X[:,i] = (col - col.mean()) / col.std()
miny = train_y.min ()
maxy = train_y.max ()
mean = (maxy + miny)/2
std = (maxy - miny)/2
train_y = (train_y - mean)/std
# Adjust for bias column
num_inputs += 1
# Adjust for sample size
eta /= num_samples
## Test various hidden node counts
for num_hidden in range (1, max_num_hidden + 1):
## Initialise weights
np.random.seed (123456)
w_hidden = 0.1*UniformRandomMatrix (num_inputs, num_hidden)
w_output = 0.1*UniformRandomMatrix (num_hidden+1, 1)
## Iterate
mse = []
for _ in range (num_epochs):
# Outputs
phi = np.append (bias, np.tanh (train_X*w_hidden), axis=1)
y = phi * w_output
err = y - train_y.transpose ()
# Gradients
g_output = phi.transpose() * err
phi_range = np.array (phi [:, range (1, num_hidden+1)])
w_output_range = w_output [range (1, num_hidden+1), 0].transpose()
err_term = np.array (err*w_output_range)
g_hidden = train_X.transpose() * np.matrix((1 - phi_range**2)*err_term)
# Update weights
w_output -= eta * g_output
w_hidden -= eta * g_hidden
mse.append (err.var ())
## Plot
plt.plot (range (num_epochs), mse, label="Hidden Neurons = {0}".format (num_hidden))
plt.legend()
plt.ylabel("MSE")
plt.xlabel("Epochs")
plt.show()
Downloading data from https://storage.googleapis.com/tensorflow/tf-keras-datasets/boston_housing.npz 57026/57026 [==============================] - 0s 0us/step
# 17.4: Bifurcation Diagrams for a Two-Neuron Module.
# Hysteresis - system is still in steady state, w12=3.
# See Figure 17.7(b).
from matplotlib import pyplot as plt
import numpy as np
# Parameters for part (i). Change parameters to get bifurcation
# diagrams for parts (ii), (iii), and (iv).
b2, w11, w21, w12, a , b = -1 , 1.5 , 5 , 3 , 1 , 0.1
start, max = -5, 10
half_N = 1999
N = 2 * half_N + 1
N1 = 1 + half_N
xs_up, xs_down = [], []
x0, y0 = -2, -8
ns_up = np.arange(half_N)
ns_down = np.arange(N1, N)
# Ramp b1 up
for n in ns_up:
b1 = start + n*max / half_N
x = b1 + w11 * np.tanh(a * x0) + w12 * np.tanh(b * y0)
y = b2 + w21 * np.tanh(a * x0)
xn = x
x0 , y0 = x , y
xs_up.append([n, xn])
xs_up = np.array(xs_up)
# Ramp b1 down
for n in ns_down:
b1 = start + 2*max - n*max / half_N
x = b1 + w11 * np.tanh(a * x0) + w12 * np.tanh(b * y0)
y = b2 + w21 * np.tanh(a * x0)
xn = x
x0 , y0 = x , y
xs_down.append([N-n, xn])
xs_down = np.array(xs_down)
fig, ax = plt.subplots()
xtick_labels = np.linspace(start, max, 7)
ax.set_xticks([(-start + x) / max * N1 for x in xtick_labels])
ax.set_xticklabels(["{:.1f}".format(xtick) for xtick in xtick_labels])
plt.plot(xs_up[:, 0], xs_up[:, 1], "r.", markersize=0.1)
plt.plot(xs_down[:, 0], xs_down[:,1], "b.", markersize=0.1)
plt.xlabel(r"$b_1$", fontsize=15)
plt.ylabel(r"$x_n$", fontsize=15)
plt.tick_params(labelsize=15)
plt.show()
There is a large bistable region for $-1 < b_1 < 2$, approximately. The solutions are in steady-state. Compare with the stability diagram in Figure 17.7(b).
# 17.4: Bifurcation Diagrams for a Two-Neuron Module.
# Periodic and quasiperiodic behaviour for w12=-5.
# See Figure 17.7(b).
from matplotlib import pyplot as plt
import numpy as np
# Parameters for part (ii). Change parameters to get bifurcation
# diagrams for parts (i), (iii), and (iv).
b2, w11, w21, w12, a , b = -1 , 1.5 , 5 , -5 , 1 , 0.1
start, max = -5, 10
half_N = 1999
N = 2 * half_N + 1
N1 = 1 + half_N
xs_up, xs_down = [], []
x0, y0 = -2, -8
ns_up = np.arange(half_N)
ns_down = np.arange(N1, N)
# Ramp b1 up
for n in ns_up:
b1 = start + n*max / half_N
x = b1 + w11 * np.tanh(a * x0) + w12 * np.tanh(b * y0)
y = b2 + w21 * np.tanh(a * x0)
xn = x
x0 , y0 = x , y
xs_up.append([n, xn])
xs_up = np.array(xs_up)
# Ramp b1 down
for n in ns_down:
b1 = start + 2*max - n*max / half_N
x = b1 + w11 * np.tanh(a * x0) + w12 * np.tanh(b * y0)
y = b2 + w21 * np.tanh(a * x0)
xn = x
x0 , y0 = x , y
xs_down.append([N-n, xn])
xs_down = np.array(xs_down)
fig, ax = plt.subplots()
xtick_labels = np.linspace(start, max, 7)
ax.set_xticks([(-start + x) / max * N1 for x in xtick_labels])
ax.set_xticklabels(["{:.1f}".format(xtick) for xtick in xtick_labels])
plt.plot(xs_up[:, 0], xs_up[:, 1], "r.", markersize=0.5)
plt.plot(xs_down[:, 0], xs_down[:,1], "b.", markersize=0.5)
plt.xlabel(r"$b_1$", fontsize=15)
plt.ylabel(r"$x_n$", fontsize=15)
plt.tick_params(labelsize=15)
plt.show()