-
Notifications
You must be signed in to change notification settings - Fork 0
/
nips_demo.py
449 lines (389 loc) · 20.8 KB
/
nips_demo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
#!/usr/bin/env python
"""
Identifying Dendritic Processing in a [Filter]-[Ideal IAF] neural circuit
This demo illustrates identification of the [Filter] in the
[Filter]-[Ideal IAF] circuit using band-limited input signals, i.e.,
signals that belong to the Paley-Wiener space.
The code below corresponds to Corollary 3 in [1] and was used to generate
Figs. 4-6 in [1]. The employed filter was taken from [2].
Author: Yevgeniy B. Slutskiy <[email protected]>
Revision Author: Chung-Heng Yeh <[email protected]>
Bionet Group, Columbia University
Copyright 2010-2012 Yevgeniy B. Slutskiy and Chung-Heng Yeh
"""
import cim
import matplotlib as mpl
#mpl.use('TKAGG')
import matplotlib.pylab as p
from scipy.signal import fftconvolve
import numpy as np
from math import sqrt, log10, factorial
from time import time
def run_time(msg,timer):
"""
Display the time between the function call and TIMER in "min' sec'" format.
"""
t = time()-timer
print msg + " %d\' %.2f\"" % (int(t/60.), t-60.*int(t/60.) )
def low_pass_filter(dt,t_sig,sig,Omega):
"""
Pass input through a low-pass-filter of bandwidth Omega.
"""
t_len = len(t_sig)
t_sinc = dt*np.arange(-t_len+1,t_len+1)
g = np.sinc(t_sinc*Omega/np.pi)*Omega/np.pi
return dt*fftconvolve(g,sig)[t_len-1:2*t_len-1]
def create_random_band_signal(dt, t_sig, Omega, seed=1):
"""
Create a random band-limited signal.
Parameters
----------
dt : float
Sampling resolution of input signal.
t_sig : array_like of floats
time course of Signal.
Omega : float
bandwidth of the signal.
seed : integer
seed for random number generator.
Returns
-------
v : ndarray of floats
returns the synthesized signal.
Notes
-----
The idea is to first randomly generate samples of signal, and then pass the
samples through a low pass filter (LPF). For readers who are not familiar
with the above concept, please refer to Shannon Sampling Theorem.
"""
np.random.seed(seed) # set the seed of random number generator
v = np.zeros_like(t_sig) # initialize the signal
# generate the random samples
gap = int( (np.pi/Omega)/dt )
v[ np.arange(0,len(v),gap) ] = np.random.randn( len(v)/gap )
# compute the low-pass-filter; the length of LPF is twice the length of
# the signal in order to cover the entire signal during convolution.
t_len = len(t_sig)
t_sinc = dt*np.arange(-t_len+1,t_len+1) # set time course of LPF
g = np.sinc(t_sinc*Omega/np.pi)*Omega/np.pi # compute the LPF
# convolve samples with LPF and extract the proper portion
v = fftconvolve(g,v)[t_len-1:2*t_len-1]
v /= np.max(np.abs(v)) # normalize the signal
return v
def plot_cim_result(dt,t_sig,sig,filt_sig,windows,spk,t_h,h,Ph,h_hat,freq,
tau,fig_name='',fig_title='',fig_size=(4,3)):
"""
Plot the CIM result.
"""
fig = p.figure(num=fig_name,figsize=fig_size,dpi=300, facecolor='w')
fig.suptitle(fig_title,fontsize=24)
# plot the input stimulus u
ax = fig.add_subplot(3,2,1,xlim=(0,t_sig[-1]-t_sig[0]),ylim=(-1,1),
ylabel='Amplitude', xlabel='Time, [s]',
title='$(a)\qquad$Input signal u(t)');
ax.plot(t_sig - t_sig[0], sig,
label='$\Omega = 2\pi\cdot '+ repr(int(freq)) + '$rad/s$\qquad$')
p.legend(loc='lower right')
#plot the periodogram power spectrum estimate of u
ax = fig.add_subplot(3,2,2,
xlim=(-150,150),ylim=(-140,0),
title='$(d)\qquad$Periodogram Power Spectrum Estimate of u(t)');
ax.psd(sig,pad_to=len(t_sig),NFFT=len(t_sig),Fs=1/dt,sides='twosided',
label='supp$(\mathcal{F}u)=[-\Omega,\Omega]\qquad$')
p.legend(loc='lower right')
# plot h, Ph and h_hat (the filter identified by the algorithm)
ax = fig.add_subplot(3,2,3,
xlim=(0,t_sig[-1]-t_sig[0]),ylim=(0,1.2), xlabel='Time, [s]',
yticks=[],yticklabels=[],
title='$(b)\qquad$Output of the [Filter]-[Ideal IAF] neural circuit');
colorset = mpl.cm.hsv(np.arange(len(windows))/float(len(windows)),1)
ax.stem(spk-t_sig[0], np.ones_like(spk),linefmt='k-',
markerfmt='k^',label='$D = 40$Hz')
for color,win in zip(colorset,windows):
ax.axvspan(win[0]-t_sig[0],win[1]-t_sig[0],facecolor=color,
edgecolor='none',label='_nolegend_')
p.legend(loc='lower right')
# plot the periodogram power spectrum estimate of h
ax = fig.add_subplot(3,2,4,
xlim=(-150,150),ylim=(-140,0),#ylabel='Power, [dB]',
title='$(e)\qquad$Periodogram Power Spectrum Estimate of h(t)');
ax.psd(h,pad_to=len(t_sig),NFFT=len(t_sig),Fs=1/dt,sides='twosided',
label='supp$(\mathcal{F}h)\supset[-\Omega,\Omega]\qquad$')
p.legend(loc='lower right')
# plot the periodogram power spectrum estimate of v=u*h
ax = fig.add_subplot(3,2,5,xlim=(t_h[0],t_h[-1]), ylim=(-50,100),
ylabel='Amplitude', xlabel='Time, [s]',
title='$(c)\qquad$Original filter vs. the identified filter')
idx = np.logical_and( t_h >= tau[0], t_h <= tau[1] )
# Normalized RMSE between h and h_hat
h_hhat_err = np.abs(h-h_hat)/max(abs(h)); # compute the absolute error
h_hhat_RMSE = sqrt(dt*np.trapz(h_hhat_err[idx]**2)/(tau[1]-tau[0])) # compute the RMSE
# Normalized RMSE between Ph and h_hat
Ph_hhat_err = np.abs(Ph-h_hat)/max(abs(Ph)); # compute the absolute error
Ph_hhat_RMSE = sqrt(dt*np.trapz(Ph_hhat_err[idx]**2)/(tau[1]-tau[0])) # compute the RMSE
ax.plot(t_h, h,'--k',label='$h,\,$RMSE$(\hat{h},h)$ = %3.2e $\qquad$' % h_hhat_RMSE )
ax.plot(t_h, Ph,'-b',label='$\mathcal{P}h,\,$RMSE$(\hat{h},\mathcal{P}h)$ = %3.2e $\qquad$' % Ph_hhat_RMSE)
ax.plot(t_h, h_hat,'-r',label='$\hat{h}$')
p.legend(loc='upper right')
# plot the periodogram power spectrum estmate of v=u*h
ax = fig.add_subplot(3,2,6,
xlim=(-150,150),ylim=(-140,0),#ylabel='Power, [dB]',
title='$(f)\qquad$Periodogram Power Spectrum Estimate of v(t)');
ax.psd(filt_sig,pad_to=len(t_sig),NFFT=len(t_sig),Fs=1/dt,sides='twosided',
label='supp$(\mathcal{F}v)=[-\Omega,\Omega]\qquad$')
p.legend(loc='lower right')
fig.text(.1,0.05,'[1] $Identifying\ Dendritic\ Processing$, ' +
'Advances in Neural Information Processing Systems 23, pp. 1261-1269, 2010')
if fig_name:
p.savefig(fig_name+'.png',dpi=300)
p.show()
if __name__=='__main__':
# Initialize the demo
tic_demo = time() # initialize the timer for the entire demo
tic_init = time() # initialize the timer for the initialization
dt = 5e-6 # set the time step, [s]
# Specify the filter h to be used
# -------------------------------
# Generate a filter according to Adelson and Bergen in [2]. h has a
# temporal support on the interval [T_1, T_2]
T_1, T_2 = 0., 0.1+dt # set the interval of the impulse response, [s]
t_filt = np.arange(T_1,T_2,dt) # set the length of the impulse response, [s]
a = 200; # set the filter parameter
h = 3*a*np.exp(-a*t_filt)*((a*t_filt)**3/factorial(3)\
-(a*t_filt)**5/factorial(5))
# set zero-padded version of h, for comparison purpose
h_long = np.concatenate( (np.zeros(len(h)//2),h,np.zeros(len(h)//2)) )
t_long = np.concatenate(( dt*np.arange(-(len(h)//2),0),t_filt,
dt*np.arange(len(h),len(h)+len(h)//2)))
# Plot the filter
p.figure(0);
p.axes(xlim=(t_filt[0],t_filt[-1]),xlabel='Time, [s]',ylabel='Amplitude',
title='Impulse response $h(t)$ of the filter');
p.plot(t_filt, h);
p.show()
# Create a band-limited stimulus. The bandwidth W = 2*pi*25 rad/s
# --------------------------------------------------------------
f = 25. # set the input signal bandwidth, [Hz]
W = 2.*np.pi*f # calculate the bandwidth in [rad]
t = np.arange(0.,1.12+dt,dt) # set the time course of the input signal
# fix the state of random number generator for reproducible results
u = create_random_band_signal(dt, t, W, seed=20130101)
# Plot the stimulus
p.figure(1)
p.axes(xlim=(t[0],t[-1]),xlabel='Time, [s]',ylabel='Amplitude',
title='Input Stimulus $u(t)$');
p.plot(t,u)
p.show()
# Compute the filter projection Ph
# --------------------------------
# Ph is the projection of h onto the input signal space. It is the best
# approximation of h that can be recovered.
t_Ph = t_long # set the time course of the projection Ph
Ph = low_pass_filter(dt,t_Ph,h_long,W) # find the projection Ph
# Plot the filter projection Ph
p.figure(2)
p.axes(xlim=(t_Ph[0],t_Ph[-1]),xlabel='Time, [s]',ylabel='Amplitude',
title='Filter $h(t)$ and Filter projection $\mathcal{P}h(t)$');
p.plot(t_Ph, h_long,'--k',label='$h\qquad$' )
p.plot(t_Ph, Ph,'-b',label='$\mathcal{P}h\qquad$')
p.legend(loc='upper right')
p.show()
# Filter the stimulus u
# -------------------------
# Since all signals are finite, the filter output v=u*h is not calculated
# properly on the boundaries. v_proper is that part of v, for which the
# convolution of u and h is computed correctly.
v_proper = dt*fftconvolve(u,h,'valid') # convolve u with h
u_proper = u[len(h)-1:] # get the proper part of v
t_proper = t[len(h)-1:] # get the corresponding time vector
# Plot the filter output
p.figure(3)
p.axes(xlim=(0,t_proper[-1]-t_proper[0]),xlabel='Time, [s]',ylabel='Amplitude',
title='Filter output $v(t)$')
p.plot(t_proper-t_proper[0],v_proper)
p.show()
# Encode the filter output v=u*h with an IAF neuron
# -------------------------------------------------
# Specify parameters of the Ideal IAF neuron
delta = 0.007 # set the threshold
bias = 0.28 # set the bias
kappa = 1. # set capacitance
# # Encode the filter output
spk_train, vol_trace = cim.iaf_encode(dt, t_proper, v_proper, b=bias, d=delta, C=kappa)
run_time('Initialization time:',tic_init) # display the initialization time
# Plot the voltage trace and the associated spike train
p.figure(4)
ax1 = p.subplot2grid((7, 1), (0, 0),rowspan=5,xlim=(0., 1.0),
ylabel='Amplitude', ylim=(min(vol_trace), delta*1.1),
title='Output of the [Filter]-[Ideal IAF] neural circuit')
p.setp(ax1.get_xticklabels(), visible=False)
ax2 = p.subplot2grid((7, 1), (5, 0),rowspan=2,sharex=ax1,
xlabel='Time, [s]', ylim=(0., 1.1),xlim=(0., 1.0),
yticks = (),yticklabels = ())
ax1.plot(t_proper-t_proper[0],vol_trace,'b-',
label='Membrane voltage $v\qquad$')
ax1.plot((0, t_proper[-1]-t_proper[0]),(delta,delta),'--r',
label='Threshold $\delta=' + repr(delta) + '$')
ax1.plot(spk_train-t_proper[0],delta*np.ones_like(spk_train),'ro',
label= '$v(t)=\delta$')
ax1.legend(loc='lower right')
ax2.stem(spk_train-t_proper[0],np.ones_like(spk_train),
linefmt='k-', markerfmt='k^', label='$(t_k)_{k\in Z}$')
ax2.legend(loc='lower right')
p.show()
# Identify the filter projection Ph
# ---------------------------------
# Since the temporal support of the filter h is not known, we identify the
# projection Ph in a window tau. Temporal windows W^i of spikes can be
# chosen arbitrarily. Here we pick W^i so that the
# Nyquist-type condition on the density of spikes is achieved quickly (see
# also Theorem 1 and Remark 1 in [1]).
T_filt_rec = 0.12; # specify the hypothesized length of the impulse response
tau = [-(T_filt_rec - t_filt[-1])/2, 0] # tau is centered around the actual temporal support T_2-T_1
tau[1] = t_filt[-1]-tau[0];
# set the maximum number of windows to be used (could be smaller depending on the simulation)
N_max = 10
# start the algorithm timer
tic = time()
# execute the CIM algorithm
[h_hat, windows] = cim.cim_ideal_iaf(dt, t_Ph, t_proper, u_proper,
W, bias, delta, kappa, spk_train, tau, N_max)
# display execution time
run_time('CIM running time for Fig.4:',tic)
# Generate Fig. 4 of [1]
# ----------------------
plot_cim_result(dt,t_proper,u_proper,v_proper,windows,spk_train,
t_Ph,h_long,Ph,h_hat,f,tau,fig_size = (20,15),fig_name = 'nips_fig_4',
fig_title='NIPS 2010 Figure 4\nA.A. Lazar and Y.B. Slutskiy')
# Generate Fig. 5 of [1]
# ----------------------
# The following procedures are same as above except that the input stimulus
# is now band-limited to 100Hz.
f = 100. # set the input signal bandwidth, [Hz]
W = 2.*np.pi*f # calculate the bandwidth in [rad]
t = np.arange(0.,1.52+dt,dt) # set the time course of the input signal
# fix the state of random number generator for reproducible results
u = create_random_band_signal(dt, t, W, seed=20130101)
# Compute the filter projection Ph
Ph = low_pass_filter(dt,t_Ph,h_long,W)
# Filter the input signal u
v_proper = dt*fftconvolve(u,h,'valid') # convolve u with h
u_proper = u[len(h)-1:] # get the proper part of v
t_proper = t[len(h)-1:] # get the corresponding time vector
# Specify parameters of the Ideal IAF neuron
delta = 0.007 # set the threshold
bias = 0.28 # set the bias
kappa = 1. # set capacitance
# Encode the filter output
spk_train, vol_trace = cim.iaf_encode(dt, t_proper, v_proper, b=bias, d=delta, C=kappa)
N_max = 15 # set the maximum number of windows to be used
tic = time() # start the algorithm timer
[h_hat, windows] = cim.cim_ideal_iaf(dt, t_Ph, t_proper, u_proper,
W, bias, delta, kappa, spk_train,tau, N_max)
# display execution time
run_time('CIM running time for Fig.5:',tic)
# Plot the results
plot_cim_result(dt,t_proper,u_proper,v_proper,windows,spk_train,
t_Ph,h_long,Ph,h_hat,f,tau,fig_size = (20,15),fig_name = 'nips_fig_5',
fig_title='NIPS 2010 Figure 5\nA.A. Lazar and Y.B. Slutskiy')
# Generate Fig. 6 of [1]
# ----------------------
# In Fig. 6a we plot the mean square error (MSE) between the filter
# projection Ph and the identified filter h_hat as a function of the number
# of temporal windows N.
#
# In Fig. 6b we plot the mean square error (MSE) between the original
# filter h and the identified filter h_hat as a function of the input
# signal bandwidth
tic_fig6 = time() # initialize the timer for Fig.6
tic_fig6a = time() # initialize the timer for Fig.6a
f = 100. # set the input signal bandwidth, [Hz]
W = 2.*np.pi*f # calculate the bandwidth in [rad]
t = np.arange(0.,8.5+dt,dt) # set the time course of the input signal
# fix the state of random number generator for reproducible results
u = create_random_band_signal(dt, t, W, seed=20130101)
# Compute the filter projection Ph
Ph = low_pass_filter(dt,t_Ph,h_long,W)
# Filter the input signal u
v_proper = dt*fftconvolve(u,h,'valid') # convolve u with h
u_proper = u[len(h)-1:] # get the proper part of v
t_proper = t[len(h)-1:] # get the corresponding time vector
# Encode the filter output v=u*h with an IAF neuron
delta = 0.007 # set the threshold
kappa = 1. # set capacitance
D = [20., 40., 60.] # set the average spiking density
N_max = 30 # set the maximum number of windows to be used
MSE_N = np.zeros((len(D),N_max)) # initialize the MSE array
idx = np.logical_and( t_Ph >= tau[0], t_Ph <= tau[1] ) # find indices of t for the given frame
for i,d in enumerate(D):
bias = d*delta/kappa; # set the bias
# Encode the filter output v=u*h with an IAF neuron
spk_train, vol = cim.iaf_encode(dt, t_proper, v_proper, b=bias, d=delta, C=kappa)
# Identify the filter projection Ph
[h_hat, windows, h_hat_N] = cim.cim_ideal_iaf(dt, t_Ph, t_proper, u_proper,
W, bias, delta, kappa, spk_train, tau, N_max,
Calc_MSE_N=True)
# Compute the normalized MSE
for j in xrange(len(windows)):
# Normalized RMSE error for Ph - h_hat
Ph_hhat_err = abs(Ph-h_hat_N[j])/max(abs(Ph))
Ph_hhat_RMSE = sqrt(dt*np.trapz(Ph_hhat_err[idx]**2)/(tau[1]-tau[0]))
MSE_N[i,j] = 10*log10( Ph_hhat_RMSE**2 )
run_time('Fig. 6a time: ', tic_fig6a) # display the running time for Fig.6a
# Get data for Fig. 6b
# --------------------
# In the following, for fixed number of windows N, we compute the MSE
# between the original filter h and the identified filter h_hat as a
# function of the input signal bandwidth
tic_fig6b = time() # initialize the timer for Fig.6b
F = np.arange(10.,160.,10.) # specify the bandwidth vector
t = np.arange(0,1.8+dt,dt) # specify the bandwidth vector
N_max = 15 # set the number of windows
# specify Ideal IAF neuron parameters
delta = 0.007 # set the threshold
bias = 0.42 # set the bias
kappa = 1. # set the capacitance
MSE_BW = np.zeros_like(F)
for i,f in enumerate(F):
W = 2*np.pi*f # calculate the bandwidth in [rad]
u = create_random_band_signal(dt, t, W, seed=20130101)
# Filter the input signal u
v_proper = dt*fftconvolve(u,h,'valid') # convolve u with h
u_proper = u[len(h)-1:] # get the proper part of v
t_proper = t[len(h)-1:] # get the corresponding time vector
# Encode the filter output v=u*h with an IAF neuron
spk_train, vol = cim.iaf_encode(dt, t_proper, v_proper, b=bias, d=delta, C=kappa)
# Identify the filter projection Ph
h_hat, windows = cim.cim_ideal_iaf(dt, t_Ph, t_proper, u_proper, W,
bias, delta, kappa, spk_train,tau, N_max)
# Normalized RMSE between h and h_hat
h_hhat_err = np.abs(h_long-h_hat)/max(abs(h_long)); # compute the absolute error
h_hhat_RMSE = sqrt(dt*np.trapz(h_hhat_err[idx]**2)/(tau[1]-tau[0]))
# Compute the RMSE in dB
MSE_BW[i] = 10*log10( h_hhat_RMSE**2 )
run_time('Fig. 6b time: ',tic_fig6b) # display the running time for Fig.6b
# Plot FIg. 6
fig = p.figure( num='NIPS 2010 Fig. 6',figsize=(12,9), dpi=300, facecolor='w')
# Plot Fig, 6-a
ax = fig.add_subplot(2,1,1,xlim=(1,MSE_N.shape[1]),
title = 'MSE$(\hat{h},\mathcal{P}h)$ vs. the Number of Temporal Windows',
xlabel = 'Number of Windows $N$',ylabel = 'MSE$(\hat{h},\mathcal{P}h)$, [dB]')
lineSpec = ['-o','-s','-d'];linecolor = ['k','b','r'];Marker_Size = [5, 5, 7]
for d,mse,spec,color,marker_size in zip(D,MSE_N,lineSpec,linecolor,Marker_Size):
ax.plot(xrange(1,len(mse)+1),mse,spec,color=color,markersize=marker_size,
label='$D = %d$ Hz' % int(d))
ax.axvline(x=2*np.pi*100/(np.pi*d), ymin=-100, ymax=40,color=color,
linestyle='--',label='$\Omega/(%d \pi)$' % d)
p.legend()
# Plot Fig, 6-b
ax = fig.add_subplot(2,1,2,xlim=(10,150),ylim=(-70,0),xticks = F,
title = 'MSE$(\hat{h},h)$ vs. input signal bandwidth',
xlabel = 'Input signal bandwidth $\Omega/(2\pi)$, [Hz]',
ylabel = 'MSE$(\hat{h},h)$, [dB]')
ax.plot( F, MSE_BW, '-bo', markersize = 5,
label='$D = 60\;$Hz, $N = %d \quad$' % N_max)
p.savefig('nips_fig_6.png',dpi=300)
# Finalize the demo
# -------------------
run_time('Total Fig. 6 time: ',tic_fig6); # display the running time for Fig.6
run_time('Demo time: ', tic_demo); # display the demo time