Converting some manualtests to tests used by installcheck

.travis.yml

@@ -27,7 +27,7 @@ language: minimal
 sudo: required
 dist: trusty
 #https://blog.travis-ci.com/2017-06-19-trusty-updates-2017-Q2
-group: edge
+
 env:
   matrix:
     # We don't have to run a full matrix here, because most of the options are

pynest/examples/cross_check_mip_corrdet.py

@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+#
+# cross_check_mip_corrdet.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/>.
+
+import nest
+from matplotlib.pylab import *
+
+'''
+Auto- and crosscorrelation functions for spike trains.
+
+A time bin of size tbin is centered around the time difference it
+represents. If the correlation function is calculated for tau in
+[-tau_max, tau_max], the pair events contributing to the left-most
+bin are those for which tau in [-tau_max-tbin/2, tau_max+tbin/2) and
+so on.
+
+Correlate two spike trains with each other assumes spike times to be ordered in
+time. tau > 0 means spike2 is later than spike1
+
+tau_max: maximum time lag in ms correlation function
+tbin:    bin size
+spike1:  first spike train [tspike...]
+spike2:  second spike train [tspike...]
+'''
+
+
+def corr_spikes_sorted(spike1, spike2, tbin, tau_max, h):
+    tau_max_i = int(tau_max / h)
+    tbin_i = int(tbin / h)
+
+    cross = zeros(int(2 * tau_max_i / tbin_i + 1), 'd')
+
+    j0 = 0
+
+    for spki in spike1:
+        j = j0
+        while j < len(spike2) and spike2[j] - spki < -tau_max_i - tbin_i / 2.0:
+            j += 1
+        j0 = j
+
+        while j < len(spike2) and spike2[j] - spki < tau_max_i + tbin_i / 2.0:
+            cross[int(
+                (spike2[j] - spki + tau_max_i + 0.5 * tbin_i) / tbin_i)] += 1.0
+            j += 1
+
+    return cross
+
+nest.ResetKernel()
+
+h = 0.1             # Computation step size in ms
+T = 100000.0        # Total duration
+delta_tau = 10.0
+tau_max = 100.0
+pc = 0.5
+nu = 100.0
+
+# grng_seed is 0 because test data was produced for seed = 0
+nest.SetKernelStatus({'local_num_threads': 1, 'resolution': h,
+                      'overwrite_files': True, 'grng_seed': 0})
+
+# Set up network, connect and simulate
+mg = nest.Create('mip_generator')
+nest.SetStatus(mg, {'rate': nu, 'p_copy': pc})
+
+cd = nest.Create('correlation_detector')
+nest.SetStatus(cd, {'tau_max': tau_max, 'delta_tau': delta_tau})
+
+sd = nest.Create('spike_detector')
+nest.SetStatus(sd, {'withtime': True,
+                    'withgid': True, 'time_in_steps': True})
+
+pn1 = nest.Create('parrot_neuron')
+pn2 = nest.Create('parrot_neuron')
+
+nest.Connect(mg, pn1)
+nest.Connect(mg, pn2)
+nest.Connect(pn1, sd)
+nest.Connect(pn2, sd)
+
+nest.SetDefaults('static_synapse', {'weight': 1.0, 'receptor_type': 0})
+nest.Connect(pn1, cd)
+
+nest.SetDefaults('static_synapse', {'weight': 1.0, 'receptor_type': 1})
+nest.Connect(pn2, cd)
+
+nest.Simulate(T)
+
+n_events = nest.GetStatus(cd)[0]['n_events']
+n1 = n_events[0]
+n2 = n_events[1]
+
+lmbd1 = (n1 / (T - tau_max)) * 1000.0
+lmbd2 = (n2 / (T - tau_max)) * 1000.0
+
+h = 0.1
+tau_max = 100.0  # ms correlation window
+t_bin = 10.0  # ms bin size
+
+spikes = nest.GetStatus(sd)[0]['events']['senders']
+
+sp1 = find(spikes[:] == 4)
+sp2 = find(spikes[:] == 5)
+
+# Find crosscorrolation
+cross = corr_spikes_sorted(sp1, sp2, t_bin, tau_max, h)
+
+print("Crosscorrelation:")
+print(cross)
+print("Sum of crosscorrelation:")
+print(sum(cross))

pynest/examples/hh_phaseplane.py

@@ -0,0 +1,164 @@
+# -*- coding: utf-8 -*-
+#
+# hh_phaseplane.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/>.
+
+'''
+  hh_phaseplane makes a numerical phase-plane analysis of the Hodgkin-Huxley
+  neuron (iaf_psc_alpha). Dynamics is investigated in the V-n space (see remark
+  below). A constant DC can be specified  and its influence on the nullclines
+  can be studied.
+
+  REMARK
+  To make the two-dimensional analysis possible, the (four-dimensional)
+  Hodgkin-Huxley formalism needs to be artificially reduced to two dimensions,
+  in this case by 'clamping' the two other variables, m an h, to
+  constant values (m_eq and h_eq).
+'''
+
+import nest
+from matplotlib import pyplot as plt
+
+
+amplitude = 100.  # Set externally applied current amplitude in pA
+dt = 0.1  # simulation step length [ms]
+
+nest.ResetKernel()
+nest.set_verbosity('M_ERROR')
+
+nest.SetKernelStatus({'resolution': dt})
+neuron = nest.Create('hh_psc_alpha')
+
+# Numerically obtain equilibrium state
+nest.Simulate(1000)
+
+m_eq = nest.GetStatus(neuron)[0]['Act_m']
+h_eq = nest.GetStatus(neuron)[0]['Act_h']
+
+nest.SetStatus(neuron, {'I_e': amplitude})  # Apply external current
+
+# Scan state space
+print('Scanning phase space')
+
+V_new_vec = []
+n_new_vec = []
+# x will contain the phase-plane data as a vector field
+x = []
+count = 0
+for V in range(-100, 42, 2):
+    n_V = []
+    n_n = []
+    for n in range(10, 81):
+        # Set V_m and n
+        nest.SetStatus(neuron, {'V_m': V*1.0, 'Inact_n': n/100.0,
+                                'Act_m': m_eq, 'Act_h': h_eq})
+        # Find state
+        V_m = nest.GetStatus(neuron)[0]['V_m']
+        Inact_n = nest.GetStatus(neuron)[0]['Inact_n']
+
+        # Simulate a short while
+        nest.Simulate(dt)
+
+        # Find difference between new state and old state
+        V_m_new = nest.GetStatus(neuron)[0]['V_m'] - V*1.0
+        Inact_n_new = nest.GetStatus(neuron)[0]['Inact_n'] - n/100.0
+
+        # Store in vector for later analysis
+        n_V.append(abs(V_m_new))
+        n_n.append(abs(Inact_n_new))
+        x.append([V_m, Inact_n, V_m_new, Inact_n_new])
+
+        if count % 10 == 0:
+            # Write updated state next to old state
+            print('')
+            print('Vm:  \t', V_m)
+            print('new Vm:\t', V_m_new)
+            print('Inact_n:', Inact_n)
+            print('new Inact_n:', Inact_n_new)
+
+        count += 1
+    # Store in vector for later analysis
+    V_new_vec.append(n_V)
+    n_new_vec.append(n_n)
+
+# Set state for AP generation
+nest.SetStatus(neuron, {'V_m': -34., 'Inact_n': 0.2,
+                        'Act_m': m_eq, 'Act_h': h_eq})
+
+print('')
+print('AP-trajectory')
+# ap will contain the trace of a single action potential as one possible
+# numerical solution in the vector field
+ap = []
+for i in range(1, 1001):
+    # Find state
+    V_m = nest.GetStatus(neuron)[0]['V_m']
+    Inact_n = nest.GetStatus(neuron)[0]['Inact_n']
+
+    if i % 10 == 0:
+        # Write new state next to old state
+        print('Vm: \t', V_m)
+        print('Inact_n:', Inact_n)
+    ap.append([V_m, Inact_n])
+
+    # Simulate again
+    nest.SetStatus(neuron, {'Act_m': m_eq, 'Act_h': h_eq})
+    nest.Simulate(dt)
+
+# Make analysis
+print('')
+print('Plot analysis')
+
+V_matrix = [list(x) for x in zip(*V_new_vec)]
+n_matrix = [list(x) for x in zip(*n_new_vec)]
+n_vec = [x/100. for x in range(10, 81)]
+V_vec = [x*1. for x in range(-100, 42, 2)]
+
+nullcline_V = []
+nullcline_n = []
+
+print('Searching nullclines')
+for i in range(0, len(V_vec)):
+    index = V_matrix[:][i].index(min(V_matrix[:][i]))
+    if index != 0 and index != len(n_vec):
+        nullcline_V.append([V_vec[i], n_vec[index]])
+
+    index = n_matrix[:][i].index(min(n_matrix[:][i]))
+    if index != 0 and index != len(n_vec):
+        nullcline_n.append([V_vec[i], n_vec[index]])
+
+print('Plotting vector field')
+factor = 0.1
+for i in range(0, count, 3):
+    plt.plot([x[i][0], x[i][0] + factor*x[i][2]],
+             [x[i][1], x[i][1] + factor*x[i][3]], color=[0.6, 0.6, 0.6])
+
+plt.plot(nullcline_V[:][0], nullcline_V[:][1], linewidth=2.0)
+plt.plot(nullcline_n[:][0], nullcline_n[:][1], linewidth=2.0)
+
+plt.xlim([V_vec[0], V_vec[-1]])
+plt.ylim([n_vec[0], n_vec[-1]])
+
+plt.plot(ap[:][0], ap[:][1], color='black', linewidth=1.0)
+
+plt.xlabel('Membrane potential V [mV]')
+plt.ylabel('Inactivation variable n')
+plt.title('Phase space of the Hodgkin-Huxley Neuron')
+
+plt.show()

pynest/examples/hh_psc_alpha.py

@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+#
+# hh_psc_alpha.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/>.
+
+"""
+Example using hh_psc_alpha
+--------------------------
+
+This example produces a rate-response (FI) curve of the Hogkin-Huxley
+neuron  in response to a range of different current (DC) stimulations.
+The result is plotted using matplotlib.
+
+Since a DC input affetcs only the neuron's channel dynamics, this routine
+does not yet check correctness of synaptic response.
+"""
+
+import nest
+import numpy as np
+import matplotlib.pyplot as plt
+
+nest.set_verbosity('M_WARNING')
+nest.ResetKernel()
+
+simtime = 1000
+
+# Amplitude range, in pA
+dcfrom = 0
+dcstep = 20
+dcto = 2000
+
+h = 0.1  # simulation step size in mS
+
+neuron = nest.Create('hh_psc_alpha')
+sd = nest.Create('spike_detector')
+
+nest.SetStatus(sd, {'to_memory': False})
+
+nest.Connect(neuron, sd, syn_spec={'weight': 1.0, 'delay': h})
+
+
+# Simulation loop
+n_data = dcto / float(dcstep)
+amplitudes = np.zeros(n_data)
+event_freqs = np.zeros(n_data)
+for i, amp in enumerate(range(dcfrom, dcto, dcstep)):
+    nest.SetStatus(neuron, {'I_e': float(amp)})
+    print("Simulating with current I={} pA".format(amp))
+    nest.Simulate(1000)  # one second warm-up time for equilibrium state
+    nest.SetStatus(sd, {'n_events': 0})  # then reset spike counts
+    nest.Simulate(simtime)  # another simulation call to record firing rate
+
+    n_events = nest.GetStatus(sd, keys={'n_events'})[0][0]
+    amplitudes[i] = amp
+    event_freqs[i] = n_events / (simtime / 1000.)
+
+plt.plot(amplitudes, event_freqs)
+plt.show()