|
@@ -2,22 +2,26 @@
|
|
|
"""
|
|
|
Unit tests for the statistics module.
|
|
|
|
|
|
-:copyright: Copyright 2014-2016 by the Elephant team, see AUTHORS.txt.
|
|
|
+:copyright: Copyright 2014-2016 by the Elephant team, see `doc/authors.rst`.
|
|
|
:license: Modified BSD, see LICENSE.txt for details.
|
|
|
"""
|
|
|
from __future__ import division
|
|
|
|
|
|
+import math
|
|
|
+import sys
|
|
|
import unittest
|
|
|
|
|
|
import neo
|
|
|
import numpy as np
|
|
|
-from numpy.testing.utils import assert_array_almost_equal, assert_array_equal
|
|
|
import quantities as pq
|
|
|
import scipy.integrate as spint
|
|
|
+from numpy.testing.utils import assert_array_almost_equal, assert_array_equal
|
|
|
|
|
|
-import elephant.statistics as es
|
|
|
import elephant.kernels as kernels
|
|
|
-import warnings
|
|
|
+from elephant import statistics
|
|
|
+from elephant.spike_train_generation import homogeneous_poisson_process
|
|
|
+
|
|
|
+python_version_major = sys.version_info.major
|
|
|
|
|
|
|
|
|
class isi_TestCase(unittest.TestCase):
|
|
@@ -25,7 +29,7 @@ class isi_TestCase(unittest.TestCase):
|
|
|
self.test_array_2d = np.array([[0.3, 0.56, 0.87, 1.23],
|
|
|
[0.02, 0.71, 1.82, 8.46],
|
|
|
[0.03, 0.14, 0.15, 0.92]])
|
|
|
- self.targ_array_2d_0 = np.array([[-0.28, 0.15, 0.95, 7.23],
|
|
|
+ self.targ_array_2d_0 = np.array([[-0.28, 0.15, 0.95, 7.23],
|
|
|
[0.01, -0.57, -1.67, -7.54]])
|
|
|
self.targ_array_2d_1 = np.array([[0.26, 0.31, 0.36],
|
|
|
[0.69, 1.11, 6.64],
|
|
@@ -39,40 +43,40 @@ class isi_TestCase(unittest.TestCase):
|
|
|
st = neo.SpikeTrain(
|
|
|
self.test_array_1d, units='ms', t_stop=10.0, t_start=0.29)
|
|
|
target = pq.Quantity(self.targ_array_1d, 'ms')
|
|
|
- res = es.isi(st)
|
|
|
+ res = statistics.isi(st)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_isi_with_quantities_1d(self):
|
|
|
st = pq.Quantity(self.test_array_1d, units='ms')
|
|
|
target = pq.Quantity(self.targ_array_1d, 'ms')
|
|
|
- res = es.isi(st)
|
|
|
+ res = statistics.isi(st)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_isi_with_plain_array_1d(self):
|
|
|
st = self.test_array_1d
|
|
|
target = self.targ_array_1d
|
|
|
- res = es.isi(st)
|
|
|
+ res = statistics.isi(st)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_isi_with_plain_array_2d_default(self):
|
|
|
st = self.test_array_2d
|
|
|
target = self.targ_array_2d_default
|
|
|
- res = es.isi(st)
|
|
|
+ res = statistics.isi(st)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_isi_with_plain_array_2d_0(self):
|
|
|
st = self.test_array_2d
|
|
|
target = self.targ_array_2d_0
|
|
|
- res = es.isi(st, axis=0)
|
|
|
+ res = statistics.isi(st, axis=0)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_isi_with_plain_array_2d_1(self):
|
|
|
st = self.test_array_2d
|
|
|
target = self.targ_array_2d_1
|
|
|
- res = es.isi(st, axis=1)
|
|
|
+ res = statistics.isi(st, axis=1)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
@@ -82,15 +86,15 @@ class isi_cv_TestCase(unittest.TestCase):
|
|
|
self.test_array_regular = np.arange(1, 6)
|
|
|
|
|
|
def test_cv_isi_regular_spiketrain_is_zero(self):
|
|
|
- st = neo.SpikeTrain(self.test_array_regular, units='ms', t_stop=10.0)
|
|
|
+ st = neo.SpikeTrain(self.test_array_regular, units='ms', t_stop=10.0)
|
|
|
targ = 0.0
|
|
|
- res = es.cv(es.isi(st))
|
|
|
+ res = statistics.cv(statistics.isi(st))
|
|
|
self.assertEqual(res, targ)
|
|
|
|
|
|
def test_cv_isi_regular_array_is_zero(self):
|
|
|
st = self.test_array_regular
|
|
|
targ = 0.0
|
|
|
- res = es.cv(es.isi(st))
|
|
|
+ res = statistics.cv(statistics.isi(st))
|
|
|
self.assertEqual(res, targ)
|
|
|
|
|
|
|
|
@@ -115,120 +119,144 @@ class mean_firing_rate_TestCase(unittest.TestCase):
|
|
|
self.targ_array_1d = self.targ_array_2d_1[0]
|
|
|
self.max_array_1d = self.max_array_2d_1[0]
|
|
|
|
|
|
+ def test_invalid_input_spiketrain(self):
|
|
|
+ # empty spiketrain
|
|
|
+ self.assertRaises(ValueError, statistics.mean_firing_rate, [])
|
|
|
+ for st_invalid in (None, 0.1):
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate,
|
|
|
+ st_invalid)
|
|
|
+
|
|
|
def test_mean_firing_rate_with_spiketrain(self):
|
|
|
st = neo.SpikeTrain(self.test_array_1d, units='ms', t_stop=10.0)
|
|
|
- target = pq.Quantity(self.targ_array_1d/10., '1/ms')
|
|
|
- res = es.mean_firing_rate(st)
|
|
|
+ target = pq.Quantity(self.targ_array_1d / 10., '1/ms')
|
|
|
+ res = statistics.mean_firing_rate(st)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
+ def test_mean_firing_rate_typical_use_case(self):
|
|
|
+ np.random.seed(92)
|
|
|
+ st = homogeneous_poisson_process(rate=100 * pq.Hz, t_stop=100 * pq.s)
|
|
|
+ rate1 = statistics.mean_firing_rate(st)
|
|
|
+ rate2 = statistics.mean_firing_rate(st, t_start=st.t_start,
|
|
|
+ t_stop=st.t_stop)
|
|
|
+ self.assertEqual(rate1.units, rate2.units)
|
|
|
+ self.assertAlmostEqual(rate1.item(), rate2.item())
|
|
|
+
|
|
|
def test_mean_firing_rate_with_spiketrain_set_ends(self):
|
|
|
st = neo.SpikeTrain(self.test_array_1d, units='ms', t_stop=10.0)
|
|
|
- target = pq.Quantity(2/0.5, '1/ms')
|
|
|
- res = es.mean_firing_rate(st, t_start=0.4, t_stop=0.9)
|
|
|
+ target = pq.Quantity(2 / 0.5, '1/ms')
|
|
|
+ res = statistics.mean_firing_rate(st, t_start=0.4 * pq.ms,
|
|
|
+ t_stop=0.9 * pq.ms)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_quantities_1d(self):
|
|
|
st = pq.Quantity(self.test_array_1d, units='ms')
|
|
|
- target = pq.Quantity(self.targ_array_1d/self.max_array_1d, '1/ms')
|
|
|
- res = es.mean_firing_rate(st)
|
|
|
+ target = pq.Quantity(self.targ_array_1d / self.max_array_1d, '1/ms')
|
|
|
+ res = statistics.mean_firing_rate(st)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_quantities_1d_set_ends(self):
|
|
|
st = pq.Quantity(self.test_array_1d, units='ms')
|
|
|
- target = pq.Quantity(2/0.6, '1/ms')
|
|
|
- res = es.mean_firing_rate(st, t_start=400*pq.us, t_stop=1.)
|
|
|
- assert_array_almost_equal(res, target, decimal=9)
|
|
|
+
|
|
|
+ # t_stop is not a Quantity
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
+ t_start=400 * pq.us, t_stop=1.)
|
|
|
+
|
|
|
+ # t_start is not a Quantity
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
+ t_start=0.4, t_stop=1. * pq.ms)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_1d(self):
|
|
|
st = self.test_array_1d
|
|
|
- target = self.targ_array_1d/self.max_array_1d
|
|
|
- res = es.mean_firing_rate(st)
|
|
|
+ target = self.targ_array_1d / self.max_array_1d
|
|
|
+ res = statistics.mean_firing_rate(st)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_1d_set_ends(self):
|
|
|
st = self.test_array_1d
|
|
|
- target = self.targ_array_1d/(1.23-0.3)
|
|
|
- res = es.mean_firing_rate(st, t_start=0.3, t_stop=1.23)
|
|
|
+ target = self.targ_array_1d / (1.23 - 0.3)
|
|
|
+ res = statistics.mean_firing_rate(st, t_start=0.3, t_stop=1.23)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_2d_default(self):
|
|
|
st = self.test_array_2d
|
|
|
- target = self.targ_array_2d_default/self.max_array_2d_default
|
|
|
- res = es.mean_firing_rate(st)
|
|
|
+ target = self.targ_array_2d_default / self.max_array_2d_default
|
|
|
+ res = statistics.mean_firing_rate(st)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_2d_0(self):
|
|
|
st = self.test_array_2d
|
|
|
- target = self.targ_array_2d_0/self.max_array_2d_0
|
|
|
- res = es.mean_firing_rate(st, axis=0)
|
|
|
+ target = self.targ_array_2d_0 / self.max_array_2d_0
|
|
|
+ res = statistics.mean_firing_rate(st, axis=0)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_2d_1(self):
|
|
|
st = self.test_array_2d
|
|
|
- target = self.targ_array_2d_1/self.max_array_2d_1
|
|
|
- res = es.mean_firing_rate(st, axis=1)
|
|
|
+ target = self.targ_array_2d_1 / self.max_array_2d_1
|
|
|
+ res = statistics.mean_firing_rate(st, axis=1)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_3d_None(self):
|
|
|
st = self.test_array_3d
|
|
|
- target = np.sum(self.test_array_3d, None)/5.
|
|
|
- res = es.mean_firing_rate(st, axis=None, t_stop=5.)
|
|
|
+ target = np.sum(self.test_array_3d, None) / 5.
|
|
|
+ res = statistics.mean_firing_rate(st, axis=None, t_stop=5.)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_3d_0(self):
|
|
|
st = self.test_array_3d
|
|
|
- target = np.sum(self.test_array_3d, 0)/5.
|
|
|
- res = es.mean_firing_rate(st, axis=0, t_stop=5.)
|
|
|
+ target = np.sum(self.test_array_3d, 0) / 5.
|
|
|
+ res = statistics.mean_firing_rate(st, axis=0, t_stop=5.)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_3d_1(self):
|
|
|
st = self.test_array_3d
|
|
|
- target = np.sum(self.test_array_3d, 1)/5.
|
|
|
- res = es.mean_firing_rate(st, axis=1, t_stop=5.)
|
|
|
+ target = np.sum(self.test_array_3d, 1) / 5.
|
|
|
+ res = statistics.mean_firing_rate(st, axis=1, t_stop=5.)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_3d_2(self):
|
|
|
st = self.test_array_3d
|
|
|
- target = np.sum(self.test_array_3d, 2)/5.
|
|
|
- res = es.mean_firing_rate(st, axis=2, t_stop=5.)
|
|
|
+ target = np.sum(self.test_array_3d, 2) / 5.
|
|
|
+ res = statistics.mean_firing_rate(st, axis=2, t_stop=5.)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_2d_1_set_ends(self):
|
|
|
st = self.test_array_2d
|
|
|
- target = np.array([4, 1, 3])/(1.23-0.14)
|
|
|
- res = es.mean_firing_rate(st, axis=1, t_start=0.14, t_stop=1.23)
|
|
|
+ target = np.array([4, 1, 3]) / (1.23 - 0.14)
|
|
|
+ res = statistics.mean_firing_rate(st, axis=1, t_start=0.14,
|
|
|
+ t_stop=1.23)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
def test_mean_firing_rate_with_plain_array_2d_None(self):
|
|
|
st = self.test_array_2d
|
|
|
- target = self.targ_array_2d_None/self.max_array_2d_None
|
|
|
- res = es.mean_firing_rate(st, axis=None)
|
|
|
+ target = self.targ_array_2d_None / self.max_array_2d_None
|
|
|
+ res = statistics.mean_firing_rate(st, axis=None)
|
|
|
assert not isinstance(res, pq.Quantity)
|
|
|
assert_array_almost_equal(res, target, decimal=9)
|
|
|
|
|
|
- def test_mean_firing_rate_with_plain_array_and_units_start_stop_typeerror(self):
|
|
|
+ def test_mean_firing_rate_with_plain_array_and_units_start_stop_typeerror(
|
|
|
+ self):
|
|
|
st = self.test_array_2d
|
|
|
- self.assertRaises(TypeError, es.mean_firing_rate, st,
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
t_start=pq.Quantity(0, 'ms'))
|
|
|
- self.assertRaises(TypeError, es.mean_firing_rate, st,
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
t_stop=pq.Quantity(10, 'ms'))
|
|
|
- self.assertRaises(TypeError, es.mean_firing_rate, st,
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
t_start=pq.Quantity(0, 'ms'),
|
|
|
t_stop=pq.Quantity(10, 'ms'))
|
|
|
- self.assertRaises(TypeError, es.mean_firing_rate, st,
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
t_start=pq.Quantity(0, 'ms'),
|
|
|
t_stop=10.)
|
|
|
- self.assertRaises(TypeError, es.mean_firing_rate, st,
|
|
|
+ self.assertRaises(TypeError, statistics.mean_firing_rate, st,
|
|
|
t_start=0.,
|
|
|
t_stop=pq.Quantity(10, 'ms'))
|
|
|
|
|
@@ -258,121 +286,196 @@ class FanoFactorTestCase(unittest.TestCase):
|
|
|
# Test with list of spiketrains
|
|
|
self.assertEqual(
|
|
|
np.var(self.sp_counts) / np.mean(self.sp_counts),
|
|
|
- es.fanofactor(self.test_spiketrains))
|
|
|
+ statistics.fanofactor(self.test_spiketrains))
|
|
|
|
|
|
# One spiketrain in list
|
|
|
st = self.test_spiketrains[0]
|
|
|
- self.assertEqual(es.fanofactor([st]), 0.0)
|
|
|
+ self.assertEqual(statistics.fanofactor([st]), 0.0)
|
|
|
|
|
|
def test_fanofactor_empty(self):
|
|
|
# Test with empty list
|
|
|
- self.assertTrue(np.isnan(es.fanofactor([])))
|
|
|
- self.assertTrue(np.isnan(es.fanofactor([[]])))
|
|
|
+ self.assertTrue(np.isnan(statistics.fanofactor([])))
|
|
|
+ self.assertTrue(np.isnan(statistics.fanofactor([[]])))
|
|
|
|
|
|
# Test with empty quantity
|
|
|
- self.assertTrue(np.isnan(es.fanofactor([] * pq.ms)))
|
|
|
+ self.assertTrue(np.isnan(statistics.fanofactor([] * pq.ms)))
|
|
|
|
|
|
# Empty spiketrain
|
|
|
st = neo.core.SpikeTrain([] * pq.ms, t_start=0 * pq.ms,
|
|
|
t_stop=1.5 * pq.ms)
|
|
|
- self.assertTrue(np.isnan(es.fanofactor(st)))
|
|
|
+ self.assertTrue(np.isnan(statistics.fanofactor(st)))
|
|
|
|
|
|
def test_fanofactor_spiketrains_same(self):
|
|
|
# Test with same spiketrains in list
|
|
|
sts = [self.test_spiketrains[0]] * 3
|
|
|
- self.assertEqual(es.fanofactor(sts), 0.0)
|
|
|
+ self.assertEqual(statistics.fanofactor(sts), 0.0)
|
|
|
|
|
|
def test_fanofactor_array(self):
|
|
|
- self.assertEqual(es.fanofactor(self.test_array),
|
|
|
+ self.assertEqual(statistics.fanofactor(self.test_array),
|
|
|
np.var(self.sp_counts) / np.mean(self.sp_counts))
|
|
|
|
|
|
def test_fanofactor_array_same(self):
|
|
|
lst = [self.test_array[0]] * 3
|
|
|
- self.assertEqual(es.fanofactor(lst), 0.0)
|
|
|
+ self.assertEqual(statistics.fanofactor(lst), 0.0)
|
|
|
|
|
|
def test_fanofactor_quantity(self):
|
|
|
- self.assertEqual(es.fanofactor(self.test_quantity),
|
|
|
+ self.assertEqual(statistics.fanofactor(self.test_quantity),
|
|
|
np.var(self.sp_counts) / np.mean(self.sp_counts))
|
|
|
|
|
|
def test_fanofactor_quantity_same(self):
|
|
|
lst = [self.test_quantity[0]] * 3
|
|
|
- self.assertEqual(es.fanofactor(lst), 0.0)
|
|
|
+ self.assertEqual(statistics.fanofactor(lst), 0.0)
|
|
|
|
|
|
def test_fanofactor_list(self):
|
|
|
- self.assertEqual(es.fanofactor(self.test_list),
|
|
|
+ self.assertEqual(statistics.fanofactor(self.test_list),
|
|
|
np.var(self.sp_counts) / np.mean(self.sp_counts))
|
|
|
|
|
|
def test_fanofactor_list_same(self):
|
|
|
lst = [self.test_list[0]] * 3
|
|
|
- self.assertEqual(es.fanofactor(lst), 0.0)
|
|
|
+ self.assertEqual(statistics.fanofactor(lst), 0.0)
|
|
|
+
|
|
|
+ @unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
|
|
|
+ def test_fanofactor_different_durations(self):
|
|
|
+ st1 = neo.SpikeTrain([1, 2, 3] * pq.s, t_stop=4 * pq.s)
|
|
|
+ st2 = neo.SpikeTrain([1, 2, 3] * pq.s, t_stop=4.5 * pq.s)
|
|
|
+ self.assertWarns(UserWarning, statistics.fanofactor, (st1, st2))
|
|
|
+
|
|
|
+ def test_fanofactor_wrong_type(self):
|
|
|
+ # warn_tolerance is not a quantity
|
|
|
+ st1 = neo.SpikeTrain([1, 2, 3] * pq.s, t_stop=4 * pq.s)
|
|
|
+ self.assertRaises(TypeError, statistics.fanofactor, [st1],
|
|
|
+ warn_tolerance=1e-4)
|
|
|
|
|
|
|
|
|
class LVTestCase(unittest.TestCase):
|
|
|
def setUp(self):
|
|
|
- self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
|
|
|
- 4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
|
|
|
- 4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
|
|
|
- 5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
|
|
|
- 2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
|
|
|
- 11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
|
|
|
- 1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
|
|
|
- 6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
|
|
|
- 6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
|
|
|
- 5, 7, 5, 6, 8, 11, 33, 10, 7, 4]
|
|
|
+ self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
|
|
|
+ 4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
|
|
|
+ 4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
|
|
|
+ 5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
|
|
|
+ 2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
|
|
|
+ 11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
|
|
|
+ 1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
|
|
|
+ 6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
|
|
|
+ 6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
|
|
|
+ 5, 7, 5, 6, 8, 11, 33, 10, 7, 4]
|
|
|
|
|
|
self.target = 0.971826029994
|
|
|
|
|
|
def test_lv_with_quantities(self):
|
|
|
seq = pq.Quantity(self.test_seq, units='ms')
|
|
|
- assert_array_almost_equal(es.lv(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.lv(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_lv_with_plain_array(self):
|
|
|
seq = np.array(self.test_seq)
|
|
|
- assert_array_almost_equal(es.lv(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.lv(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_lv_with_list(self):
|
|
|
seq = self.test_seq
|
|
|
- assert_array_almost_equal(es.lv(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.lv(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_lv_raise_error(self):
|
|
|
seq = self.test_seq
|
|
|
- self.assertRaises(AttributeError, es.lv, [])
|
|
|
- self.assertRaises(AttributeError, es.lv, 1)
|
|
|
- self.assertRaises(ValueError, es.lv, np.array([seq, seq]))
|
|
|
+ self.assertRaises(ValueError, statistics.lv, [])
|
|
|
+ self.assertRaises(ValueError, statistics.lv, 1)
|
|
|
+ self.assertRaises(ValueError, statistics.lv, np.array([seq, seq]))
|
|
|
+
|
|
|
+ @unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
|
|
|
+ def test_2short_spike_train(self):
|
|
|
+ seq = [1]
|
|
|
+ with self.assertWarns(UserWarning):
|
|
|
+ """
|
|
|
+ Catches UserWarning: Input size is too small. Please provide
|
|
|
+ an input with more than 1 entry.
|
|
|
+ """
|
|
|
+ self.assertTrue(math.isnan(statistics.lv(seq, with_nan=True)))
|
|
|
+
|
|
|
+
|
|
|
+class LVRTestCase(unittest.TestCase):
|
|
|
+ def setUp(self):
|
|
|
+ self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
|
|
|
+ 4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
|
|
|
+ 4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
|
|
|
+ 5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
|
|
|
+ 2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
|
|
|
+ 11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
|
|
|
+ 1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
|
|
|
+ 6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
|
|
|
+ 6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
|
|
|
+ 5, 7, 5, 6, 8, 11, 33, 10, 7, 4]
|
|
|
+
|
|
|
+ self.target = 2.1845363464753134
|
|
|
+
|
|
|
+ def test_lvr_with_quantities(self):
|
|
|
+ seq = pq.Quantity(self.test_seq, units='ms')
|
|
|
+ assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)
|
|
|
+
|
|
|
+ def test_lvr_with_plain_array(self):
|
|
|
+ seq = np.array(self.test_seq)
|
|
|
+ assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)
|
|
|
+
|
|
|
+ def test_lvr_with_list(self):
|
|
|
+ seq = self.test_seq
|
|
|
+ assert_array_almost_equal(statistics.lvr(seq), self.target, decimal=9)
|
|
|
+
|
|
|
+ def test_lvr_raise_error(self):
|
|
|
+ seq = self.test_seq
|
|
|
+ self.assertRaises(ValueError, statistics.lvr, [])
|
|
|
+ self.assertRaises(ValueError, statistics.lvr, 1)
|
|
|
+ self.assertRaises(ValueError, statistics.lvr, np.array([seq, seq]))
|
|
|
+ self.assertRaises(ValueError, statistics.lvr, seq, -1)
|
|
|
+
|
|
|
+ @unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
|
|
|
+ def test_lvr_refractoriness_kwarg(self):
|
|
|
+ seq = np.array(self.test_seq)
|
|
|
+ with self.assertWarns(UserWarning):
|
|
|
+ assert_array_almost_equal(statistics.lvr(seq, R=5),
|
|
|
+ self.target, decimal=9)
|
|
|
+
|
|
|
+ @unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
|
|
|
+ def test_2short_spike_train(self):
|
|
|
+ seq = [1]
|
|
|
+ with self.assertWarns(UserWarning):
|
|
|
+ """
|
|
|
+ Catches UserWarning: Input size is too small. Please provide
|
|
|
+ an input with more than 1 entry.
|
|
|
+ """
|
|
|
+ self.assertTrue(math.isnan(statistics.lvr(seq, with_nan=True)))
|
|
|
+
|
|
|
|
|
|
|
|
|
class CV2TestCase(unittest.TestCase):
|
|
|
def setUp(self):
|
|
|
- self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
|
|
|
- 4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
|
|
|
- 4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
|
|
|
- 5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
|
|
|
- 2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
|
|
|
- 11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
|
|
|
- 1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
|
|
|
- 6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
|
|
|
- 6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
|
|
|
- 5, 7, 5, 6, 8, 11, 33, 10, 7, 4]
|
|
|
+ self.test_seq = [1, 28, 4, 47, 5, 16, 2, 5, 21, 12,
|
|
|
+ 4, 12, 59, 2, 4, 18, 33, 25, 2, 34,
|
|
|
+ 4, 1, 1, 14, 8, 1, 10, 1, 8, 20,
|
|
|
+ 5, 1, 6, 5, 12, 2, 8, 8, 2, 8,
|
|
|
+ 2, 10, 2, 1, 1, 2, 15, 3, 20, 6,
|
|
|
+ 11, 6, 18, 2, 5, 17, 4, 3, 13, 6,
|
|
|
+ 1, 18, 1, 16, 12, 2, 52, 2, 5, 7,
|
|
|
+ 6, 25, 6, 5, 3, 15, 4, 3, 16, 3,
|
|
|
+ 6, 5, 24, 21, 3, 3, 4, 8, 4, 11,
|
|
|
+ 5, 7, 5, 6, 8, 11, 33, 10, 7, 4]
|
|
|
|
|
|
self.target = 1.0022235296529176
|
|
|
|
|
|
def test_cv2_with_quantities(self):
|
|
|
seq = pq.Quantity(self.test_seq, units='ms')
|
|
|
- assert_array_almost_equal(es.cv2(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.cv2(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_cv2_with_plain_array(self):
|
|
|
seq = np.array(self.test_seq)
|
|
|
- assert_array_almost_equal(es.cv2(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.cv2(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_cv2_with_list(self):
|
|
|
seq = self.test_seq
|
|
|
- assert_array_almost_equal(es.cv2(seq), self.target, decimal=9)
|
|
|
+ assert_array_almost_equal(statistics.cv2(seq), self.target, decimal=9)
|
|
|
|
|
|
def test_cv2_raise_error(self):
|
|
|
seq = self.test_seq
|
|
|
- self.assertRaises(AttributeError, es.cv2, [])
|
|
|
- self.assertRaises(AttributeError, es.cv2, 1)
|
|
|
- self.assertRaises(AttributeError, es.cv2, np.array([seq, seq]))
|
|
|
+ self.assertRaises(ValueError, statistics.cv2, [])
|
|
|
+ self.assertRaises(ValueError, statistics.cv2, 1)
|
|
|
+ self.assertRaises(ValueError, statistics.cv2, np.array([seq, seq]))
|
|
|
|
|
|
|
|
|
class RateEstimationTestCase(unittest.TestCase):
|
|
@@ -384,34 +487,31 @@ class RateEstimationTestCase(unittest.TestCase):
|
|
|
self.st_margin = 5.0 # seconds
|
|
|
self.st_rate = 10.0 # Hertz
|
|
|
|
|
|
- st_num_spikes = np.random.poisson(
|
|
|
- self.st_rate*(self.st_dur-2*self.st_margin))
|
|
|
- spike_train = np.random.rand(
|
|
|
- st_num_spikes) * (self.st_dur-2*self.st_margin) + self.st_margin
|
|
|
- spike_train.sort()
|
|
|
+ np.random.seed(19)
|
|
|
+ duration_effective = self.st_dur - 2 * self.st_margin
|
|
|
+ st_num_spikes = np.random.poisson(self.st_rate * duration_effective)
|
|
|
+ spike_train = sorted(
|
|
|
+ np.random.rand(st_num_spikes) *
|
|
|
+ duration_effective +
|
|
|
+ self.st_margin)
|
|
|
|
|
|
# convert spike train into neo objects
|
|
|
- self.spike_train = neo.SpikeTrain(spike_train*pq.s,
|
|
|
- t_start=self.st_tr[0]*pq.s,
|
|
|
- t_stop=self.st_tr[1]*pq.s)
|
|
|
+ self.spike_train = neo.SpikeTrain(spike_train * pq.s,
|
|
|
+ t_start=self.st_tr[0] * pq.s,
|
|
|
+ t_stop=self.st_tr[1] * pq.s)
|
|
|
|
|
|
# generation of a multiply used specific kernel
|
|
|
- self.kernel = kernels.TriangularKernel(sigma=0.03*pq.s)
|
|
|
+ self.kernel = kernels.TriangularKernel(sigma=0.03 * pq.s)
|
|
|
|
|
|
+ @unittest.skipUnless(python_version_major == 3, "assertWarns requires 3.2")
|
|
|
def test_instantaneous_rate_and_warnings(self):
|
|
|
st = self.spike_train
|
|
|
- sampling_period = 0.01*pq.s
|
|
|
- with warnings.catch_warnings(record=True) as w:
|
|
|
- inst_rate = es.instantaneous_rate(
|
|
|
+ sampling_period = 0.01 * pq.s
|
|
|
+ with self.assertWarns(UserWarning):
|
|
|
+ # Catches warning: The width of the kernel was adjusted to a
|
|
|
+ # minimally allowed width.
|
|
|
+ inst_rate = statistics.instantaneous_rate(
|
|
|
st, sampling_period, self.kernel, cutoff=0)
|
|
|
- message1 = "The width of the kernel was adjusted to a minimally " \
|
|
|
- "allowed width."
|
|
|
- message2 = "Instantaneous firing rate approximation contains " \
|
|
|
- "negative values, possibly caused due to machine " \
|
|
|
- "precision errors."
|
|
|
- warning_message = [str(m.message) for m in w]
|
|
|
- self.assertTrue(message1 in warning_message)
|
|
|
- self.assertTrue(message2 in warning_message)
|
|
|
self.assertIsInstance(inst_rate, neo.core.AnalogSignal)
|
|
|
self.assertEqual(
|
|
|
inst_rate.sampling_period.simplified, sampling_period.simplified)
|
|
@@ -421,108 +521,191 @@ class RateEstimationTestCase(unittest.TestCase):
|
|
|
|
|
|
def test_error_instantaneous_rate(self):
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=[1, 2, 3]*pq.s,
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel)
|
|
|
+ TypeError, statistics.instantaneous_rate,
|
|
|
+ spiketrain=[1, 2, 3] * pq.s,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=[1, 2, 3],
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel)
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=[1, 2, 3],
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel)
|
|
|
st = self.spike_train
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
sampling_period=0.01, kernel=self.kernel)
|
|
|
self.assertRaises(
|
|
|
- ValueError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=-0.01*pq.ms, kernel=self.kernel)
|
|
|
+ ValueError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=-0.01 * pq.ms, kernel=self.kernel)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=0.01*pq.ms, kernel='NONE')
|
|
|
- self.assertRaises(TypeError, es.instantaneous_rate, self.spike_train,
|
|
|
- sampling_period=0.01*pq.s, kernel='wrong_string',
|
|
|
- t_start=self.st_tr[0]*pq.s, t_stop=self.st_tr[1]*pq.s,
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel='NONE')
|
|
|
+ self.assertRaises(TypeError, statistics.instantaneous_rate,
|
|
|
+ self.spike_train,
|
|
|
+ sampling_period=0.01 * pq.s, kernel='wrong_string',
|
|
|
+ t_start=self.st_tr[0] * pq.s,
|
|
|
+ t_stop=self.st_tr[1] * pq.s,
|
|
|
trim=False)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel, cutoff=20*pq.ms)
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel,
|
|
|
+ cutoff=20 * pq.ms)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel, t_start=2)
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel, t_start=2)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel, t_stop=20*pq.mV)
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel,
|
|
|
+ t_stop=20 * pq.mV)
|
|
|
self.assertRaises(
|
|
|
- TypeError, es.instantaneous_rate, spiketrain=st,
|
|
|
- sampling_period=0.01*pq.ms, kernel=self.kernel, trim=1)
|
|
|
+ TypeError, statistics.instantaneous_rate, spiketrain=st,
|
|
|
+ sampling_period=0.01 * pq.ms, kernel=self.kernel, trim=1)
|
|
|
|
|
|
def test_rate_estimation_consistency(self):
|
|
|
"""
|
|
|
Test, whether the integral of the rate estimation curve is (almost)
|
|
|
equal to the number of spikes of the spike train.
|
|
|
"""
|
|
|
- kernel_types = [obj for obj in kernels.__dict__.values()
|
|
|
- if isinstance(obj, type) and
|
|
|
- issubclass(obj, kernels.Kernel) and
|
|
|
- hasattr(obj, "_evaluate") and
|
|
|
- obj is not kernels.Kernel and
|
|
|
- obj is not kernels.SymmetricKernel]
|
|
|
- kernel_list = [kernel_type(sigma=0.5*pq.s, invert=False)
|
|
|
- for kernel_type in kernel_types]
|
|
|
- kernel_resolution = 0.01*pq.s
|
|
|
- for kernel in kernel_list:
|
|
|
- rate_estimate_a0 = es.instantaneous_rate(self.spike_train,
|
|
|
- sampling_period=kernel_resolution,
|
|
|
- kernel='auto',
|
|
|
- t_start=self.st_tr[0]*pq.s,
|
|
|
- t_stop=self.st_tr[1]*pq.s,
|
|
|
- trim=False)
|
|
|
-
|
|
|
- rate_estimate0 = es.instantaneous_rate(self.spike_train,
|
|
|
- sampling_period=kernel_resolution,
|
|
|
- kernel=kernel)
|
|
|
-
|
|
|
- rate_estimate1 = es.instantaneous_rate(self.spike_train,
|
|
|
- sampling_period=kernel_resolution,
|
|
|
- kernel=kernel,
|
|
|
- t_start=self.st_tr[0]*pq.s,
|
|
|
- t_stop=self.st_tr[1]*pq.s,
|
|
|
- trim=False)
|
|
|
-
|
|
|
- rate_estimate2 = es.instantaneous_rate(self.spike_train,
|
|
|
- sampling_period=kernel_resolution,
|
|
|
- kernel=kernel,
|
|
|
- t_start=self.st_tr[0]*pq.s,
|
|
|
- t_stop=self.st_tr[1]*pq.s,
|
|
|
- trim=True)
|
|
|
- # test consistency
|
|
|
- rate_estimate_list = [rate_estimate0, rate_estimate1,
|
|
|
- rate_estimate2, rate_estimate_a0]
|
|
|
-
|
|
|
- for rate_estimate in rate_estimate_list:
|
|
|
+ kernel_types = tuple(
|
|
|
+ kern_cls for kern_cls in kernels.__dict__.values()
|
|
|
+ if isinstance(kern_cls, type) and
|
|
|
+ issubclass(kern_cls, kernels.Kernel) and
|
|
|
+ kern_cls is not kernels.Kernel and
|
|
|
+ kern_cls is not kernels.SymmetricKernel)
|
|
|
+ kernels_available = [kern_cls(sigma=0.5 * pq.s, invert=False)
|
|
|
+ for kern_cls in kernel_types]
|
|
|
+ kernels_available.append('auto')
|
|
|
+ kernel_resolution = 0.01 * pq.s
|
|
|
+ for kernel in kernels_available:
|
|
|
+ for center_kernel in (False, True):
|
|
|
+ rate_estimate = statistics.instantaneous_rate(
|
|
|
+ self.spike_train,
|
|
|
+ sampling_period=kernel_resolution,
|
|
|
+ kernel=kernel,
|
|
|
+ t_start=self.st_tr[0] * pq.s,
|
|
|
+ t_stop=self.st_tr[1] * pq.s,
|
|
|
+ trim=False,
|
|
|
+ center_kernel=center_kernel)
|
|
|
num_spikes = len(self.spike_train)
|
|
|
- auc = spint.cumtrapz(y=rate_estimate.magnitude[:, 0],
|
|
|
- x=rate_estimate.times.rescale('s').magnitude)[-1]
|
|
|
- self.assertAlmostEqual(num_spikes, auc, delta=0.05*num_spikes)
|
|
|
+ auc = spint.cumtrapz(
|
|
|
+ y=rate_estimate.magnitude.squeeze(),
|
|
|
+ x=rate_estimate.times.simplified.magnitude)[-1]
|
|
|
+ self.assertAlmostEqual(num_spikes, auc,
|
|
|
+ delta=0.01 * num_spikes)
|
|
|
+
|
|
|
+ def test_not_center_kernel(self):
|
|
|
+ # issue 107
|
|
|
+ t_spike = 1 * pq.s
|
|
|
+ st = neo.SpikeTrain([t_spike], t_start=0 * pq.s, t_stop=2 * pq.s,
|
|
|
+ units=pq.s)
|
|
|
+ kernel = kernels.AlphaKernel(200 * pq.ms)
|
|
|
+ fs = 0.1 * pq.ms
|
|
|
+ rate = statistics.instantaneous_rate(st,
|
|
|
+ sampling_period=fs,
|
|
|
+ kernel=kernel,
|
|
|
+ center_kernel=False)
|
|
|
+ rate_nonzero_index = np.nonzero(rate > 1e-6)[0]
|
|
|
+ # where the mass is concentrated
|
|
|
+ rate_mass = rate.times.rescale(t_spike.units)[rate_nonzero_index]
|
|
|
+ all_after_response_onset = (rate_mass >= t_spike).all()
|
|
|
+ self.assertTrue(all_after_response_onset)
|
|
|
+
|
|
|
+ def test_regression_288(self):
|
|
|
+ np.random.seed(9)
|
|
|
+ sampling_period = 200 * pq.ms
|
|
|
+ spiketrain = homogeneous_poisson_process(10 * pq.Hz,
|
|
|
+ t_start=0 * pq.s,
|
|
|
+ t_stop=10 * pq.s)
|
|
|
+ kernel = kernels.AlphaKernel(sigma=5 * pq.ms, invert=True)
|
|
|
+ rate = statistics.instantaneous_rate(spiketrain,
|
|
|
+ sampling_period=sampling_period,
|
|
|
+ kernel=kernel)
|
|
|
+ self.assertEqual(
|
|
|
+ len(rate), (spiketrain.t_stop / sampling_period).simplified.item())
|
|
|
+
|
|
|
+ # 3 Hz is not a target - it's meant to test the non-negativity of the
|
|
|
+ # result rate; ideally, for smaller sampling rates, the integral
|
|
|
+ # should match the num. of spikes in the spiketrain
|
|
|
+ self.assertGreater(rate.mean(), 3 * pq.Hz)
|
|
|
+
|
|
|
+ def test_spikes_on_edges(self):
|
|
|
+ # this test demonstrates that the trimming (convolve valid mode)
|
|
|
+ # removes the edge spikes, underestimating the true firing rate and
|
|
|
+ # thus is not able to reconstruct the number of spikes in a
|
|
|
+ # spiketrain (see test_rate_estimation_consistency)
|
|
|
+ cutoff = 5
|
|
|
+ sampling_period = 0.01 * pq.s
|
|
|
+ t_spikes = np.array([-cutoff, cutoff]) * pq.s
|
|
|
+ spiketrain = neo.SpikeTrain(t_spikes, t_start=t_spikes[0],
|
|
|
+ t_stop=t_spikes[-1])
|
|
|
+ kernel_types = tuple(
|
|
|
+ kern_cls for kern_cls in kernels.__dict__.values()
|
|
|
+ if isinstance(kern_cls, type) and
|
|
|
+ issubclass(kern_cls, kernels.Kernel) and
|
|
|
+ kern_cls is not kernels.Kernel and
|
|
|
+ kern_cls is not kernels.SymmetricKernel)
|
|
|
+ kernels_available = [kern_cls(sigma=1 * pq.s, invert=False)
|
|
|
+ for kern_cls in kernel_types]
|
|
|
+ for kernel in kernels_available:
|
|
|
+ for center_kernel in (False, True):
|
|
|
+ rate = statistics.instantaneous_rate(
|
|
|
+ spiketrain,
|
|
|
+ sampling_period=sampling_period,
|
|
|
+ kernel=kernel,
|
|
|
+ cutoff=cutoff, trim=True,
|
|
|
+ center_kernel=center_kernel)
|
|
|
+ assert_array_almost_equal(rate.magnitude, 0, decimal=3)
|
|
|
+
|
|
|
+ def test_trim_as_convolve_mode(self):
|
|
|
+ cutoff = 5
|
|
|
+ sampling_period = 0.01 * pq.s
|
|
|
+ t_spikes = np.linspace(-cutoff, cutoff, num=(2 * cutoff + 1)) * pq.s
|
|
|
+ spiketrain = neo.SpikeTrain(t_spikes, t_start=t_spikes[0],
|
|
|
+ t_stop=t_spikes[-1])
|
|
|
+ kernel = kernels.RectangularKernel(sigma=1 * pq.s)
|
|
|
+ assert cutoff > kernel.min_cutoff, "Choose larger cutoff"
|
|
|
+ kernel_types = tuple(
|
|
|
+ kern_cls for kern_cls in kernels.__dict__.values()
|
|
|
+ if isinstance(kern_cls, type) and
|
|
|
+ issubclass(kern_cls, kernels.SymmetricKernel) and
|
|
|
+ kern_cls is not kernels.SymmetricKernel)
|
|
|
+ kernels_symmetric = [kern_cls(sigma=1 * pq.s, invert=False)
|
|
|
+ for kern_cls in kernel_types]
|
|
|
+ for kernel in kernels_symmetric:
|
|
|
+ for trim in (False, True):
|
|
|
+ rate_centered = statistics.instantaneous_rate(
|
|
|
+ spiketrain, sampling_period=sampling_period,
|
|
|
+ kernel=kernel, cutoff=cutoff, trim=trim)
|
|
|
+
|
|
|
+ rate_convolve = statistics.instantaneous_rate(
|
|
|
+ spiketrain, sampling_period=sampling_period,
|
|
|
+ kernel=kernel, cutoff=cutoff, trim=trim,
|
|
|
+ center_kernel=False)
|
|
|
+ assert_array_almost_equal(rate_centered, rate_convolve)
|
|
|
|
|
|
def test_instantaneous_rate_spiketrainlist(self):
|
|
|
- st_num_spikes = np.random.poisson(
|
|
|
- self.st_rate*(self.st_dur-2*self.st_margin))
|
|
|
- spike_train2 = np.random.rand(
|
|
|
- st_num_spikes) * (self.st_dur - 2 * self.st_margin) + self.st_margin
|
|
|
- spike_train2.sort()
|
|
|
+ np.random.seed(19)
|
|
|
+ duration_effective = self.st_dur - 2 * self.st_margin
|
|
|
+ st_num_spikes = np.random.poisson(self.st_rate * duration_effective)
|
|
|
+ spike_train2 = sorted(
|
|
|
+ np.random.rand(st_num_spikes) *
|
|
|
+ duration_effective +
|
|
|
+ self.st_margin)
|
|
|
spike_train2 = neo.SpikeTrain(spike_train2 * pq.s,
|
|
|
t_start=self.st_tr[0] * pq.s,
|
|
|
t_stop=self.st_tr[1] * pq.s)
|
|
|
- st_rate_1 = es.instantaneous_rate(self.spike_train,
|
|
|
- sampling_period=0.01*pq.s,
|
|
|
- kernel=self.kernel)
|
|
|
- st_rate_2 = es.instantaneous_rate(spike_train2,
|
|
|
- sampling_period=0.01*pq.s,
|
|
|
- kernel=self.kernel)
|
|
|
- combined_rate = es.instantaneous_rate([self.spike_train, spike_train2],
|
|
|
- sampling_period=0.01*pq.s,
|
|
|
- kernel=self.kernel)
|
|
|
+ st_rate_1 = statistics.instantaneous_rate(self.spike_train,
|
|
|
+ sampling_period=0.01 * pq.s,
|
|
|
+ kernel=self.kernel)
|
|
|
+ st_rate_2 = statistics.instantaneous_rate(spike_train2,
|
|
|
+ sampling_period=0.01 * pq.s,
|
|
|
+ kernel=self.kernel)
|
|
|
+ combined_rate = statistics.instantaneous_rate(
|
|
|
+ [self.spike_train, spike_train2],
|
|
|
+ sampling_period=0.01 * pq.s,
|
|
|
+ kernel=self.kernel)
|
|
|
summed_rate = st_rate_1 + st_rate_2 # equivalent for identical kernels
|
|
|
- for a, b in zip(combined_rate.magnitude, summed_rate.magnitude):
|
|
|
- self.assertAlmostEqual(a, b, delta=0.0001)
|
|
|
+ # 'time_vector.dtype' in instantaneous_rate() is changed from float64
|
|
|
+ # to float32 which results in 3e-6 abs difference
|
|
|
+ assert_array_almost_equal(combined_rate.magnitude,
|
|
|
+ summed_rate.magnitude, decimal=5)
|
|
|
|
|
|
# Regression test for #144
|
|
|
def test_instantaneous_rate_regression_144(self):
|
|
@@ -530,7 +713,37 @@ class RateEstimationTestCase(unittest.TestCase):
|
|
|
# other, that the optimal kernel cannot be detected. Therefore, the
|
|
|
# function should react with a ValueError.
|
|
|
st = neo.SpikeTrain([2.12, 2.13, 2.15] * pq.s, t_stop=10 * pq.s)
|
|
|
- self.assertRaises(ValueError, es.instantaneous_rate, st, 1 * pq.ms)
|
|
|
+ self.assertRaises(ValueError, statistics.instantaneous_rate, st,
|
|
|
+ 1 * pq.ms)
|
|
|
+
|
|
|
+ # Regression test for #245
|
|
|
+ def test_instantaneous_rate_regression_245(self):
|
|
|
+ # This test makes sure that the correct kernel width is chosen when
|
|
|
+ # selecting 'auto' as kernel
|
|
|
+ spiketrain = neo.SpikeTrain(
|
|
|
+ range(1, 30) * pq.ms, t_start=0 * pq.ms, t_stop=30 * pq.ms)
|
|
|
+
|
|
|
+ # This is the correct procedure to attain the kernel: first, the result
|
|
|
+ # of sskernel retrieves the kernel bandwidth of an optimal Gaussian
|
|
|
+ # kernel in terms of its standard deviation sigma, then uses this value
|
|
|
+ # directly in the function for creating the Gaussian kernel
|
|
|
+ kernel_width_sigma = statistics.optimal_kernel_bandwidth(
|
|
|
+ spiketrain.magnitude, times=None, bootstrap=False)['optw']
|
|
|
+ kernel = kernels.GaussianKernel(kernel_width_sigma * spiketrain.units)
|
|
|
+ result_target = statistics.instantaneous_rate(
|
|
|
+ spiketrain, 10 * pq.ms, kernel=kernel)
|
|
|
+
|
|
|
+ # Here, we check if the 'auto' argument leads to the same operation. In
|
|
|
+ # the regression, it was incorrectly assumed that the sskernel()
|
|
|
+ # function returns the actual bandwidth of the kernel, which is defined
|
|
|
+ # as approximately bandwidth = sigma * 5.5 = sigma * (2 * 2.75).
|
|
|
+ # factor 2.0 connects kernel width with its half width,
|
|
|
+ # factor 2.7 connects half width of Gaussian distribution with
|
|
|
+ # 99% probability mass with its standard deviation.
|
|
|
+ result_automatic = statistics.instantaneous_rate(
|
|
|
+ spiketrain, 10 * pq.ms, kernel='auto')
|
|
|
+
|
|
|
+ assert_array_almost_equal(result_target, result_automatic)
|
|
|
|
|
|
|
|
|
class TimeHistogramTestCase(unittest.TestCase):
|
|
@@ -549,49 +762,53 @@ class TimeHistogramTestCase(unittest.TestCase):
|
|
|
|
|
|
def test_time_histogram(self):
|
|
|
targ = np.array([4, 2, 1, 1, 2, 2, 1, 0, 1, 0])
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=pq.s)
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains, bin_size=pq.s)
|
|
|
assert_array_equal(targ, histogram.magnitude[:, 0])
|
|
|
|
|
|
def test_time_histogram_binary(self):
|
|
|
targ = np.array([2, 2, 1, 1, 2, 2, 1, 0, 1, 0])
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=pq.s,
|
|
|
- binary=True)
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains, bin_size=pq.s,
|
|
|
+ binary=True)
|
|
|
assert_array_equal(targ, histogram.magnitude[:, 0])
|
|
|
|
|
|
def test_time_histogram_tstart_tstop(self):
|
|
|
# Start, stop short range
|
|
|
targ = np.array([2, 1])
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=pq.s,
|
|
|
- t_start=5 * pq.s, t_stop=7 * pq.s)
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains, bin_size=pq.s,
|
|
|
+ t_start=5 * pq.s,
|
|
|
+ t_stop=7 * pq.s)
|
|
|
assert_array_equal(targ, histogram.magnitude[:, 0])
|
|
|
|
|
|
# Test without t_stop
|
|
|
targ = np.array([4, 2, 1, 1, 2, 2, 1, 0, 1, 0])
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=1 * pq.s,
|
|
|
- t_start=0 * pq.s)
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains,
|
|
|
+ bin_size=1 * pq.s,
|
|
|
+ t_start=0 * pq.s)
|
|
|
assert_array_equal(targ, histogram.magnitude[:, 0])
|
|
|
|
|
|
# Test without t_start
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=1 * pq.s,
|
|
|
- t_stop=10 * pq.s)
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains,
|
|
|
+ bin_size=1 * pq.s,
|
|
|
+ t_stop=10 * pq.s)
|
|
|
assert_array_equal(targ, histogram.magnitude[:, 0])
|
|
|
|
|
|
def test_time_histogram_output(self):
|
|
|
# Normalization mean
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=pq.s,
|
|
|
- output='mean')
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains, bin_size=pq.s,
|
|
|
+ output='mean')
|
|
|
targ = np.array([4, 2, 1, 1, 2, 2, 1, 0, 1, 0], dtype=float) / 2
|
|
|
assert_array_equal(targ.reshape(targ.size, 1), histogram.magnitude)
|
|
|
|
|
|
# Normalization rate
|
|
|
- histogram = es.time_histogram(self.spiketrains, binsize=pq.s,
|
|
|
- output='rate')
|
|
|
+ histogram = statistics.time_histogram(self.spiketrains, bin_size=pq.s,
|
|
|
+ output='rate')
|
|
|
assert_array_equal(histogram.view(pq.Quantity),
|
|
|
targ.reshape(targ.size, 1) * 1 / pq.s)
|
|
|
|
|
|
# Normalization unspecified, raises error
|
|
|
- self.assertRaises(ValueError, es.time_histogram, self.spiketrains,
|
|
|
- binsize=pq.s, output=' ')
|
|
|
+ self.assertRaises(ValueError, statistics.time_histogram,
|
|
|
+ self.spiketrains,
|
|
|
+ bin_size=pq.s, output=' ')
|
|
|
|
|
|
|
|
|
class ComplexityPdfTestCase(unittest.TestCase):
|
|
@@ -613,12 +830,13 @@ class ComplexityPdfTestCase(unittest.TestCase):
|
|
|
|
|
|
def test_complexity_pdf(self):
|
|
|
targ = np.array([0.92, 0.01, 0.01, 0.06])
|
|
|
- complexity = es.complexity_pdf(self.spiketrains, binsize=0.1*pq.s)
|
|
|
+ complexity = statistics.complexity_pdf(self.spiketrains,
|
|
|
+ bin_size=0.1 * pq.s)
|
|
|
assert_array_equal(targ, complexity.magnitude[:, 0])
|
|
|
self.assertEqual(1, complexity.magnitude[:, 0].sum())
|
|
|
- self.assertEqual(len(self.spiketrains)+1, len(complexity))
|
|
|
+ self.assertEqual(len(self.spiketrains) + 1, len(complexity))
|
|
|
self.assertIsInstance(complexity, neo.AnalogSignal)
|
|
|
- self.assertEqual(complexity.units, 1*pq.dimensionless)
|
|
|
+ self.assertEqual(complexity.units, 1 * pq.dimensionless)
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|