-
Notifications
You must be signed in to change notification settings - Fork 9
/
calibrate.py
executable file
·134 lines (105 loc) · 4.63 KB
/
calibrate.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
#!/bin/python
#-----------------------------------------------------------------------------
# File Name : calibrate.py
# Purpose: Calibrates the parameters of the sigmoid function by fitting the parameters gamma and beta to the transfer function of the IF neuron
#
# Author: Emre Neftci
#
# Creation Date : 24-04-2013
# Last Modified : Thu 22 Jan 2015 11:07:24 AM PST
#
# Copyright : (c) UCSD, Emre Neftci, Srinjoy Das, Bruno Pedroni, Kenneth Kreutz-Delgado, Gert Cauwenberghs
# Licence : GPLv2
#-----------------------------------------------------------------------------
import meta_parameters
meta_parameters.parameters_script = 'parameters_calibrate'
from common import *
N = 20
t_sim = 5000*ms
def wrap_run(bias, t_ref):
defaultclock.reinit()
#----------------------------------------- RBM parameters
#------------------------------------------ Neuron Groups
i_inj = 0*amp #We are trying to calibrate gamma and beta, so i-inj cannot be calculated
eqs = Equations(eqs_str_lif_wnrd, #neuron equations are described in neusa.experimentLib
Cm = cm,
I_inj = bias,
g = g_leak,
sigma = wnsigma,
tau_rec = tau_rec)
neuron_group = NeuronGroup(\
N,
model = eqs,
threshold = 'v>theta*volt',
refractory = t_ref,
reset = 0*volt
)
#set injection currents
#---------------------- Connections and Synapses
#Inject Noise for neural sampling
#--------------------------- Monitors
M = SpikeMonitor(neuron_group)
net = Network([neuron_group, M])
print "running"
net.run(t_sim)
mmon_slice = time_slice(M, t_start=.5*second, t_stop=t_sim)
return float(len(mmon_slice.spikes)) / N / (t_sim-.5*second)
def wrap_run_notref(bias):
return wrap_run(bias, t_ref =0.)
def wrap_run_tref(bias):
return wrap_run(bias, t_ref = t_ref)
if __name__ == '__main__':
#Number of points of the transfer curve
N_run = 12
import multiprocessing
pool = multiprocessing.Pool(8)
d = et.mksavedir()
rates = np.linspace(-6000e-12,00e-12, N_run)
pool_out = pool.map(wrap_run_notref, rates)
rates_out = np.array(pool_out)
rates_sigm = np.linspace(-6000e-12,2000e-12, N_run)
pool_out_sigm = pool.map(wrap_run_tref, rates_sigm)
rates_out_sigm = np.array(pool_out_sigm)
idx = (rates_out>30) * (rates_out<.7/t_ref)
P = -polyfit(rates[idx], log(rates_out[idx]), 1)
idx = (rates_out_sigm>40) * (rates_out_sigm<1./t_ref-40)
P = polyfit(rates_sigm[idx], log(rates_out_sigm[idx]**-1-t_ref), 1)
#Generate plot
from plot_options import *
matplotlib.rcParams['font.size']=15.0
matplotlib.rcParams['figure.subplot.left'] = .2
matplotlib.rcParams['figure.subplot.bottom'] = .27
matplotlib.rcParams['figure.subplot.right'] = .92
figure(figsize = (5.0,3.0))
title('$(\\tau_r = 0\\mathrm{ms})$')
dx_plot = range(np.nonzero(rates_out>.5)[0][0],np.nonzero(rates_out>75)[0][0])
plot(rates, np.exp(-P[1]-P[0]*rates), linewidth=2, color='k', label='$\\gamma\,\\exp(\\beta)$')
plot(rates, rates_out, 'o', markersize=3, color='b',markeredgecolor='b', label='$\\rho(I)$')
xlabel('$I_{inj}\mathrm{[nA]}$',fontsize=17)
ylabel('$\\nu\mathrm{[Hz]}$',fontsize=17)
xticks( xticks()[0][::2],
xticks()[0][::2]*1e9)
yticks([0,500])
ylim([-15,500])
legend(frameon=False, loc=2,prop={'size':15},numpoints=1,borderaxespad=0.1,handletextpad=0.2)
et.savefig('exp.png', format='png', dpi=300)
figure(figsize = (5.0,3.0))
title('$(\\tau_r = 4\\mathrm{ms})$')
plot(rates_sigm, (1./t_ref)/(1+np.exp(P[0]*rates_sigm+P[1]-np.log(t_ref))), linewidth=2, color='k', label='$\\left(\\tau+\\gamma^{-1}\\exp(-\\beta)\\right)^{-1}$')
plot(rates_sigm, rates_out_sigm, 'o', markersize=3, mfc='b',markeredgecolor='b',label='$\\rho(I)$')
xlabel('$I_{inj}\mathrm{[nA]}$',fontsize=17)
ylabel('$\\nu\mathrm{[Hz]}$',fontsize=17)
ylim([-10,60])
yticks([0,int(1./t_ref)/2,int(1./t_ref)])
xticks( np.concatenate([rates_sigm[[0,-1]],xticks()[0]])[::2],
np.concatenate([rates_sigm[[0,-1]],xticks()[0]])[::2]*1e9)
xlim([rates_sigm[0],rates_sigm[-1]])
legend(frameon=False, loc=2,prop={'size':15},numpoints=1,borderaxespad=0.1,handletextpad=0.2)
et.savefig('sigmoid.png', format='png', dpi=300)
et.globaldata.pool_out = pool_out
et.globaldata.pool_out_sigm = pool_out_sigm
et.globaldata.P = P
et.save()
print 'Fitted values are gamma: {0} beta: {1}'.format(np.exp(-P[1]),-P[0])
show()
#reconstruct pdf