In [1]:
# call code to simulate action potentials via FHN biophysical model
import simFHN
import scipy as sp
In [2]:
# pass in parameters to generate and plot the simulated data and phase portrait of the system
%matplotlib inline
t = sp.arange(0.0, 200, .25)
a = 0.7
b = 0.8
[V, w2] = simFHN.simFN(a,b,t,True,.4)
In [3]:
# generate noisy data and plot observations from two neurons over true intracellular membrane potential
import numpy as np
import matplotlib.pyplot as plt
obs1 = V + np.random.normal(0,.1,len(t))
obs2 = V + np.random.normal(0,.15,len(t))
plt.subplot(121)
time = np.arange((len(t)))
lo = plt.plot(time, obs1, 'purple', time, V, 'red')
plt.xlabel('time')
plt.ylabel('signal')
plt.title('noisy measurements vs intracellular membrane potential')
plt.legend(lo, ('measurement 1','intracellular voltage'), loc='lower left')
plt.subplot(122)
lo = plt.plot(time,obs2, 'green', time, V, 'red')
plt.xlabel('time')
plt.ylabel('signal')
plt.title('noisy measurements vs intracellular membrane potential')
plt.legend(lo, ('measurement 2','intracellular voltage'), loc='lower left')
plt.subplots_adjust(right=2.5, hspace=.95)
In [4]:
# import auxiliary particle filter code
from apf20 import *
n_particles = 1000
In [5]:
# proposal_cov_mat = 0.075*np.eye
# g_mean = np.dot(B,X[i,t-1,:])
# fhn order switched, minus sign switched to plus in discretization
# Fixed Euler Discretization (was an extra 1 in I)
# delta_t = 0.25
# I_ext = 0.4
In [6]:
import numpy as np
#A = np.diag([.9, .8]) # A is not used by the algorithm
Sigma = .15*np.asarray([[1, .15],[.15, 1]])
Gamma = .12*np.asarray([[1, .15], [.15, 1]])
B = np.diag([1,3])#np.diag([2,2])
T = len(t)
x_0 = [0,0]#[0,0]
#obs2 = np.ones(len(obs1))
Obs = np.asarray([obs2]).T
I_ext = 0.4
In [7]:
# run particle filter
import timeit
start_time = timeit.default_timer()
[w, x, k] = apf(Obs, T, n_particles, 15, B, Sigma, Gamma, x_0, I_ext)
elapsed = timeit.default_timer() - start_time
In [8]:
print "time elapsed: ", elapsed, "seconds or", (elapsed/60.0), "minutes", "\ntime per iteration: ", elapsed/T
In [9]:
# visualize parameters
import matplotlib.pyplot as plt
%matplotlib inline
#parts = np.array([np.array(xi) for xi in w])
plt.subplot(141)
plt.imshow(w)
plt.xlabel('time')
plt.ylabel('particle weights')
plt.title('weight matrix')
plt.subplot(142)
plt.imshow(x[:,:,0])
plt.xlabel('time')
plt.ylabel('particles')
plt.title('path matrix')
plt.subplot(143)
plt.imshow(x[:,:,1])
plt.xlabel('time')
plt.ylabel('particles')
plt.title('path matrix')
plt.subplot(144)
plt.imshow(k)
plt.xlabel('time')
plt.ylabel('p(y_n | x_{n-1})')
plt.title('posterior')
plt.subplots_adjust(right=2.5, hspace=.75)
In [10]:
# examine particle trajectories over time
plt.subplot(141)
plt.plot(np.transpose(x[:,:,0]), alpha=.01, linewidth=1.5)
plt.xlabel('time')
plt.ylabel('displacement')
plt.title('particle path trajectories over time (dim 1)')
plt.subplot(142)
plt.plot(np.transpose(x[:,:,1]), alpha=.01, linewidth=1.5)
plt.xlabel('time')
plt.ylabel('displacement')
plt.title('particle path trajectories over time (dim 2)')
plt.subplot(143)
plt.plot(x[:,:,0])
plt.xlabel('particle')
plt.ylabel('time')
plt.title('particle variance (dim 1)')
plt.subplot(144)
plt.plot(x[:,:,1])
plt.xlabel('particle')
plt.ylabel('time')
plt.title('particle variance (dim 2)')
plt.subplots_adjust(right=2.5, hspace=.85)
In [11]:
# average over particle trajectories to obtain predicted state means for APF output
predsignal1 = np.mean(x[:,:,0], axis=0)
predsignal2 = np.mean(x[:,:,1], axis=0)
In [12]:
x.shape
Out[12]:
In [21]:
# check true values against standard kalman filter output
time = np.arange(T)
#plt.subplot(121)
#lo = plt.plot(time, hidden, 'r', time, filtered_state_means, 'b')
#plt.xlabel('time')
#plt.ylabel('signal')
#plt.title('kalman filter')
#plt.legend(lo, ('true value','prediction'))
plt.subplot(121)
plt.title('apf recovering V')
lo = plt.plot(time, V, 'r', time, predsignal1, 'b')
plt.xlabel('time')
plt.ylabel('signal')
plt.legend(lo, ('true value','prediction'), loc='lower left')
plt.subplot(122)
plt.title('apf recovering w')
lo = plt.plot(time, w2, 'r', time, predsignal2, 'b')
plt.xlabel('time')
plt.ylabel('signal')
plt.legend(lo, ('true value','prediction'), loc='lower left')
plt.subplots_adjust(right=1.5, hspace=.75)
In [14]:
# shift and scale the signal
# TO DO: update code...
predsignal3 = predsignal2 + 0.4
w3 = w2[25:800]
predsignal4 = predsignal3[0:775]
print len(w3), len(predsignal4)
In [15]:
def moving_average(a, n=10) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
In [22]:
# Recovered signal with smoothing
plt.subplot(121)
plt.title('apf recovering V')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal1))
plt.plot(V)
plt.legend(lo, ('true value','prediction'), loc='lower left')
plt.subplot(122)
plt.title('apf recovering w')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal4))
plt.plot(w3)
plt.legend(lo, ('true value','prediction'), loc='lower left')
plt.subplots_adjust(right=1.5, hspace=.75)
In [20]:
# three spike dynamics
plt.subplot(121)
plt.title('apf recovering V')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal1[0:400]))
plt.plot(V[0:400])
plt.subplot(122)
plt.title('apf recovering w')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal4[0:400]))
plt.plot(w3[0:400])
plt.subplots_adjust(right=1.5, hspace=.75)