1# -*- coding: utf-8 -*-
2#
3# if_curve.py
4#
5# This file is part of NEST.
6#
7# Copyright (C) 2004 The NEST Initiative
8#
9# NEST is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 2 of the License, or
12# (at your option) any later version.
13#
14# NEST is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with NEST.  If not, see <http://www.gnu.org/licenses/>.
21
22"""
23IF curve example
24----------------
25
26This example illustrates how to measure the I-F curve of a neuron.
27The program creates a small group of neurons and injects a noisy current
28:math:`I(t) = I_mean + I_std*W(t)`
29where :math:`W(t)` is a white noise process.
30The program systematically drives the current through a series of  values in
31the two-dimensional `(I_mean, I_std)` space and measures the firing rate of
32the neurons.
33
34In this example, we measure the I-F curve of the adaptive exponential
35integrate and fire neuron (``aeif_cond_exp``), but any other neuron model that
36accepts current inputs is possible. The model and its parameters are
37supplied when the IF_curve object is created.
38
39"""
40
41import numpy
42import nest
43import shelve
44
45###############################################################################
46# Here we define which model and the neuron parameters to use for measuring
47# the transfer function.
48
49model = 'aeif_cond_exp'
50params = {'a': 4.0,
51          'b': 80.8,
52          'V_th': -50.4,
53          'Delta_T': 2.0,
54          'I_e': 0.0,
55          'C_m': 281.0,
56          'g_L': 30.0,
57          'V_reset': -70.6,
58          'tau_w': 144.0,
59          't_ref': 5.0,
60          'V_peak': -40.0,
61          'E_L': -70.6,
62          'E_ex': 0.,
63          'E_in': -70.}
64
65
66class IF_curve():
67
68    t_inter_trial = 200.  # Interval between two successive measurement trials
69    t_sim = 1000.         # Duration of a measurement trial
70    n_neurons = 100       # Number of neurons
71    n_threads = 4         # Nubmer of threads to run the simulation
72
73    def __init__(self, model, params=None):
74        self.model = model
75        self.params = params
76        self.build()
77        self.connect()
78
79    def build(self):
80        #######################################################################
81        #  We reset NEST to delete information from previous simulations
82        # and adjust the number of threads.
83
84        nest.ResetKernel()
85        nest.local_num_threads = self.n_threads
86
87        #######################################################################
88        # We create neurons and devices with specified parameters.
89
90        self.neuron = nest.Create(self.model, self.n_neurons, self.params)
91        self.noise = nest.Create('noise_generator')
92        self.spike_recorder = nest.Create('spike_recorder')
93
94    def connect(self):
95        #######################################################################
96        # We connect the noisy current to the neurons and the neurons to
97        # the spike recorders.
98
99        nest.Connect(self.noise, self.neuron, 'all_to_all')
100        nest.Connect(self.neuron, self.spike_recorder, 'all_to_all')
101
102    def output_rate(self, mean, std):
103        self.build()
104        self.connect()
105
106        #######################################################################
107        # We adjust the parameters of the noise according to the current
108        # values.
109
110        self.noise.set(mean=mean, std=std, start=0.0, stop=1000., origin=0.)
111
112        # We simulate the network and calculate the rate.
113
114        nest.Simulate(self.t_sim)
115        rate = self.spike_recorder.n_events * 1000. / (1. * self.n_neurons * self.t_sim)
116        return rate
117
118    def compute_transfer(self, i_mean=(400.0, 900.0, 50.0),
119                         i_std=(0.0, 600.0, 50.0)):
120        #######################################################################
121        # We loop through all possible combinations of `(I_mean, I_sigma)`
122        # and measure the output rate of the neuron.
123
124        self.i_range = numpy.arange(*i_mean)
125        self.std_range = numpy.arange(*i_std)
126        self.rate = numpy.zeros((self.i_range.size, self.std_range.size))
127        nest.set_verbosity('M_WARNING')
128        for n, i in enumerate(self.i_range):
129            print('I  =  {0}'.format(i))
130            for m, std in enumerate(self.std_range):
131                self.rate[n, m] = self.output_rate(i, std)
132
133
134transfer = IF_curve(model, params)
135transfer.compute_transfer()
136
137###############################################################################
138# After the simulation is finished, we store the data into a file for
139# later analysis.
140
141with shelve.open(model + '_transfer.dat') as dat:
142    dat['I_mean'] = transfer.i_range
143    dat['I_std'] = transfer.std_range
144    dat['rate'] = transfer.rate
145