In [1]:
# call code to simulate action potentials via FHN biophysical model
import simFHN as fhn
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, 100, .5)
a = 0.7
b = 0.8
[V, w2] = fhn.simFN(a,b,t,True,1)
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 apf_fhn import *
n_particles = 500
In [5]:
import numpy as np
Sigma = .15*np.asarray([[1, .15],[.15, 1]])
Gamma = .12*np.asarray([[1, .15], [.15, 1]])
B = np.diag([1,3])
T = len(t)
x_0 = [0,0]#[0,0]
Obs = np.asarray([obs1]).T
I_ext = 1
In [6]:
# run particle filter
import timeit
start_time = timeit.default_timer()
[w, x, k] = apf(Obs, T, n_particles, 10, B, Sigma, Gamma, x_0, I_ext)
elapsed = timeit.default_timer() - start_time
In [7]:
print "time elapsed: ", elapsed, "seconds or", (elapsed/60.0), "minutes", "\ntime per iteration: ", elapsed/T
In [8]:
# 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 [9]:
# 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 [10]:
# 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 [11]:
x.shape
Out[11]:
In [12]:
# check raw signal before applying smoothing or shifting
time = np.arange(T)
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'))
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'))
plt.subplots_adjust(right=1.5, hspace=.75)
In [20]:
# shift and scale the signal
# TO DO: update code...
predsignal3 = predsignal2 + I_ext
w3 = w2[20:800]
predsignal4 = predsignal3[0:780]
print len(w2), len(predsignal4)
In [23]:
# define a moving average to smooth the signals
def moving_average(a, n=7) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
In [24]:
# Smoothed Signal
plt.subplot(121)
plt.title('apf recovering V')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal1))
plt.plot(V)
plt.subplot(122)
plt.title('apf recovering w')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal2))
plt.plot(w2)
plt.subplots_adjust(right=1.5, hspace=.85)
In [25]:
# Shifted and Scaled
plt.subplot(121)
plt.title('apf recovering V')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal1))
plt.plot(V)
plt.subplot(122)
plt.title('apf recovering w')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(moving_average(predsignal4))
plt.plot(w3)
plt.subplots_adjust(right=1.5, hspace=.85)