-
Notifications
You must be signed in to change notification settings - Fork 9
/
calibrate_SRM.py
executable file
·130 lines (104 loc) · 4.61 KB
/
calibrate_SRM.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
#!/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 12:07:15 PM 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_rd,
Cm = cm,
I_inj = bias,
g = g_leak,
tau_rec = tau_rec)
print "Creating Population"
#beta arbitrary choice
beta = 2024467142.18
gamma = 8183.40
neuron_group_hidden = NeuronGroup((N),
model=eqs,
threshold = SimpleFunThreshold( exp_prob_beta_gamma(defaultclock.dt, beta, g_leak, gamma), state = 'v'),
refractory = t_ref)
#set injection currents
#---------------------- Connections and Synapses
#Inject Noise for neural sampling
#--------------------------- Monitors
Mh = SpikeMonitor(neuron_group_hidden)
net = Network([neuron_group_hidden, Mh])
print "running"
net.run(t_sim)
mmon_slice = time_slice(Mh, 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()
show()
#reconstruct pdf