# -*- coding: utf-8 -*-
#
# test_growth_curves.py
#
# This file is part of NEST.
#
# Copyright (C) 2004 The NEST Initiative
#
# NEST is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#
# NEST is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NEST. If not, see <http://www.gnu.org/licenses/>.
from scipy.integrate import quad
import math
import numpy
from numpy import testing
import unittest
import nest
import time
HAVE_OPENMP = nest.sli_func("is_threaded")
[docs]class SynapticElementIntegrator(object):
"""
Generic class which describes how to compute the number of
Synaptic Element based on Ca value
Each derived class should overwrite the get_se(self, t) method
"""
def __init__(self, tau_ca=10000.0, beta_ca=0.001):
"""
Constructor
:param tau_ca (float): time constant of Ca decay
:param beta_ca (float): each spike increase Ca value by this value
"""
self.tau_ca = tau_ca
self.beta_ca = beta_ca
self.t_minus = 0
self.ca_minus = 0
self.se_minus = 0
[docs] def reset(self):
self.t_minus = 0
self.ca_minus = 0
self.se_minus = 0
[docs] def handle_spike(self, t):
"""
Add beta_ca to the value of Ca at t = spike time
Also update the number of synaptic element
:param t (float): spike time
"""
assert t >= self.t_minus
# Update the number of synaptic element
self.se_minus = self.get_se(t)
# update Ca value
self.ca_minus = self.get_ca(t) + self.beta_ca
self.t_minus = t
[docs] def get_ca(self, t):
"""
:param t (float): current time
:return: Ca value
"""
assert t >= self.t_minus
ca = self.ca_minus * math.exp((self.t_minus - t) / self.tau_ca)
if ca > 0:
return ca
else:
return 0
[docs] def get_se(self, t):
"""
:param t (float): current time
:return: Number of synaptic element
Should be overwritten
"""
return 0.0
[docs]class LinearExactSEI(SynapticElementIntegrator):
"""
Compute the number of synaptic element corresponding to a
linear growth curve
dse/dCa = nu * (1 - Ca/eps)
Use the exact solution
"""
def __init__(self, eps=0.7, growth_rate=1.0, *args, **kwargs):
"""
Constructor
:param eps: fix point
:param growth_rate: scaling of the growth curve
.. seealso:: SynapticElementIntegrator()
"""
super(LinearExactSEI, self).__init__(*args, **kwargs)
self.eps = eps
self.growth_rate = growth_rate
[docs] def get_se(self, t):
"""
:param t (float): current time
:return: Number of synaptic element
"""
assert t >= self.t_minus
se = 1 / self.eps * (
self.growth_rate * self.tau_ca * (
self.get_ca(t) - self.ca_minus
) + self.growth_rate * self.eps * (t - self.t_minus)
) + self.se_minus
if se > 0:
return se
else:
return 0
[docs]class LinearNumericSEI(SynapticElementIntegrator):
"""
Compute the number of synaptic element corresponding to a
linear growth curve
dse/dCa = nu * (1 - Ca/eps)
Use numerical integration (see scipy.integrate.quad)
"""
def __init__(self, eps=0.7, growth_rate=1.0, *args, **kwargs):
"""
Constructor
:param eps: fix point
:param growth_rate: scaling of the growth curve
.. seealso:: SynapticElementIntegrator()
"""
super(LinearNumericSEI, self).__init__(*args, **kwargs)
self.eps = eps
self.growth_rate = growth_rate
[docs] def get_se(self, t):
"""
:param t (float): current time
:return: Number of synaptic element
"""
assert t >= self.t_minus
se = self.se_minus + quad(self.growth_curve, self.t_minus, t)[0]
if se > 0:
return se
else:
return 0
[docs] def growth_curve(self, t):
return self.growth_rate * (1.0 - (self.get_ca(t) / self.eps))
[docs]class GaussianNumericSEI(SynapticElementIntegrator):
"""
Compute the number of synaptic element corresponding to a
linear growth curve
dse/dCa = nu * (2 * exp( ((Ca - xi)/zeta)^2 ) - 1)
with:
xi = (eta + eps) / 2.0
zeta = (eta - eps) / (2.0 * sqrt(ln(2.0)))
Use numerical integration (see scipy.integrate.quad)
"""
def __init__(self, eta=0.1, eps=0.7, growth_rate=1.0, *args, **kwargs):
"""
Constructor
:param eps: low fix point
:param eta: high fix point
:param growth_rate: scaling of the growth curve
.. seealso:: SynapticElementIntegrator()
"""
super(GaussianNumericSEI, self).__init__(*args, **kwargs)
self.zeta = (eta - eps) / (2.0 * math.sqrt(math.log(2.0)))
self.xi = (eta + eps) / 2.0
self.growth_rate = growth_rate
[docs] def get_se(self, t):
"""
:param t (float): current time
:return: Number of synaptic element
"""
assert t >= self.t_minus
se = self.se_minus + quad(self.growth_curve, self.t_minus, t)[0]
if se > 0:
return se
else:
return 0
[docs] def growth_curve(self, t):
return self.growth_rate * (
2 * math.exp(
- math.pow((self.get_ca(t) - self.xi) / self.zeta, 2)
) - 1
)
[docs]class SigmoidNumericSEI(SynapticElementIntegrator):
"""
Compute the number of synaptic element corresponding to a
sigmoid growth curve
dse/dCa = nu * ((2.0 / exp( (Ca - eps)/psi)) - 1.0)
Use numerical integration (see scipy.integrate.quad)
"""
def __init__(self, eps=0.7, growth_rate=1.0, psi=0.1, *args, **kwargs):
"""
Constructor
:param eps: set point
:param psi: controls width of growth curve
:param growth_rate: scaling of the growth curve
.. seealso:: SynapticElementIntegrator()
"""
super(SigmoidNumericSEI, self).__init__(*args, **kwargs)
self.eps = eps
self.psi = psi
self.growth_rate = growth_rate
[docs] def get_se(self, t):
"""
:param t (float): current time
:return: Number of synaptic element
"""
assert t >= self.t_minus
se = self.se_minus + quad(self.growth_curve, self.t_minus, t)[0]
if se > 0:
return se
else:
return 0
[docs] def growth_curve(self, t):
return self.growth_rate * (
(2.0 / (1.0 + math.exp(
(self.get_ca(t) - self.eps) / self.psi
))) - 1.0
)
[docs]@unittest.skipIf(not HAVE_OPENMP, 'NEST was compiled without multi-threading')
class TestGrowthCurve(unittest.TestCase):
"""
Unittest class to test the GrowthCurve used with nest
"""
[docs] def setUp(self):
nest.ResetKernel()
nest.SetKernelStatus({"total_num_virtual_procs": 4})
nest.ResetNetwork()
nest.set_verbosity('M_DEBUG')
self.sim_time = 10000
self.sim_step = 100
nest.SetKernelStatus(
{'structural_plasticity_update_interval': self.sim_time + 1})
self.se_integrator = []
self.sim_steps = None
self.ca_nest = None
self.ca_python = None
self.se_nest = None
self.se_python = None
# build
self.pop = nest.Create('iaf_psc_alpha', 10)
local_nodes = nest.GetNodes([0], {'model': 'iaf_psc_alpha'}, True)
self.local_nodes = local_nodes[0]
self.spike_detector = nest.Create('spike_detector')
nest.Connect(self.pop, self.spike_detector, 'all_to_all')
noise = nest.Create('poisson_generator')
nest.SetStatus(noise, {"rate": 800000.0})
nest.Connect(noise, self.pop, 'all_to_all')
[docs] def simulate(self):
self.sim_steps = numpy.arange(0, self.sim_time, self.sim_step)
self.ca_nest = numpy.zeros(
(len(self.local_nodes), len(self.sim_steps)))
self.ca_python = numpy.zeros(
(len(self.se_integrator), len(self.sim_steps)))
self.se_nest = numpy.zeros(
(len(self.local_nodes), len(self.sim_steps)))
self.se_python = numpy.zeros(
(len(self.se_integrator), len(self.sim_steps)))
start = time.clock()
for t_i, t in enumerate(self.sim_steps):
for n_i, n in enumerate(self.local_nodes):
self.ca_nest[n_i][t_i], synaptic_elements = nest.GetStatus(
[n], ('Ca', 'synaptic_elements'))[0]
self.se_nest[n_i][t_i] = synaptic_elements['se']['z']
nest.Simulate(self.sim_step)
start = time.clock()
tmp = nest.GetStatus(self.spike_detector, 'events')[0]
spikes_all = tmp['times']
senders_all = tmp['senders']
for n_i, n in enumerate(self.local_nodes):
spikes = spikes_all[senders_all == n]
[sei.reset() for sei in self.se_integrator]
spike_i = 0
for t_i, t in enumerate(self.sim_steps):
while spike_i < len(spikes) and spikes[spike_i] <= t:
[sei.handle_spike(spikes[spike_i])
for sei in self.se_integrator]
spike_i += 1
for sei_i, sei in enumerate(self.se_integrator):
self.ca_python[sei_i, t_i] = sei.get_ca(t)
self.se_python[sei_i, t_i] = sei.get_se(t)
for sei_i, sei in enumerate(self.se_integrator):
testing.assert_almost_equal(
self.ca_nest[n_i], self.ca_python[sei_i], decimal=5)
testing.assert_almost_equal(
self.se_nest[n_i], self.se_python[sei_i], decimal=5)
[docs] def test_linear_growth_curve(self):
beta_ca = 0.0001
tau_ca = 10000.0
growth_rate = 0.0001
eps = 0.10
nest.SetStatus(
self.local_nodes,
{
'beta_Ca': beta_ca,
'tau_Ca': tau_ca,
'synaptic_elements': {
'se': {
'growth_curve': 'linear',
'growth_rate': growth_rate,
'eps': eps,
'z': 0.0
}
}
}
)
self.se_integrator.append(LinearExactSEI(
tau_ca=tau_ca, beta_ca=beta_ca, eps=eps, growth_rate=growth_rate))
self.se_integrator.append(LinearNumericSEI(
tau_ca=tau_ca, beta_ca=beta_ca, eps=eps, growth_rate=growth_rate))
self.simulate()
# check that we got the same values from one run to another
# expected = self.se_nest[:, 10]
# print(self.se_nest[:, 10].__repr__())
expected = numpy.array([
0.08376263, 0.08374046, 0.08376031, 0.08376756, 0.08375428,
0.08378699, 0.08376784, 0.08369779, 0.08374215, 0.08370484
])
for n in self.pop:
if n in self.local_nodes:
testing.assert_almost_equal(
self.se_nest[self.local_nodes.index(n), 10], expected[
self.pop.index(n)],
decimal=8)
[docs] def test_gaussian_growth_curve(self):
beta_ca = 0.0001
tau_ca = 10000.0
growth_rate = 0.0001
eta = 0.05
eps = 0.10
nest.SetStatus(
self.local_nodes,
{
'beta_Ca': beta_ca,
'tau_Ca': tau_ca,
'synaptic_elements': {
'se': {
'growth_curve': 'gaussian',
'growth_rate': growth_rate,
'eta': eta, 'eps': eps, 'z': 0.0
}
}
}
)
self.se_integrator.append(
GaussianNumericSEI(tau_ca=tau_ca, beta_ca=beta_ca,
eta=eta, eps=eps, growth_rate=growth_rate))
self.simulate()
# check that we got the same values from one run to another
# expected = self.se_nest[:, 30]
# print(self.se_nest[:, 30].__repr__())
expected = numpy.array([
0.10044035, 0.10062526, 0.1003149, 0.10046311, 0.1005713,
0.10031755, 0.10032216, 0.10040191, 0.10058179, 0.10068598
])
for n in self.pop:
if n in self.local_nodes:
testing.assert_almost_equal(
self.se_nest[self.local_nodes.index(n), 30], expected[
self.pop.index(n)],
decimal=5)
[docs] def test_sigmoid_growth_curve(self):
beta_ca = 0.0001
tau_ca = 10000.0
growth_rate = 0.0001
eps = 0.10
psi = 0.10
nest.SetStatus(
self.local_nodes,
{
'beta_Ca': beta_ca,
'tau_Ca': tau_ca,
'synaptic_elements': {
'se': {
'growth_curve': 'sigmoid',
'growth_rate': growth_rate,
'eps': eps, 'psi': 0.1, 'z': 0.0
}
}
}
)
self.se_integrator.append(
SigmoidNumericSEI(tau_ca=tau_ca, beta_ca=beta_ca,
eps=eps, psi=psi, growth_rate=growth_rate))
self.simulate()
# check that we got the same values from one run to another
# expected = self.se_nest[:, 30]
# print(self.se_nest[:, 30].__repr__())
expected = numpy.array([
0.07801164, 0.07796841, 0.07807825, 0.07797382, 0.07802574,
0.07805961, 0.07808139, 0.07794451, 0.07799474, 0.07794458
])
for n in self.pop:
if n in self.local_nodes:
testing.assert_almost_equal(
self.se_nest[self.local_nodes.index(n), 30], expected[
self.pop.index(n)],
decimal=5)
[docs]def suite():
test_suite = unittest.makeSuite(TestGrowthCurve, 'test')
return test_suite
if __name__ == '__main__':
unittest.main()