Skip to content

Commit

Permalink
fixed 115 (NeuralEnsemble#179)
Browse files Browse the repository at this point in the history 
*  fixed NeuralEnsemble#115
  • Loading branch information
@EmanueleLucrezia @alperyeg
EmanueleLucrezia authored and alperyeg committed Oct 24, 2018
1 parent d203397 commit 74ef4a2
Show file tree
Hide file tree
Showing 2 changed files with 70 additions and 33 deletions.
86 changes: 59 additions & 27 deletions elephant/statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ def fanofactor(spiketrains):
return fano


def lv(v):
def lv(v, with_nan=False):
"""
Calculate the measure of local variation LV for
a sequence of time intervals between events.
Expand All @@ -214,6 +214,12 @@ def lv(v):
v : quantity array, numpy array or list
Vector of consecutive time intervals
with_nans : bool, optional
If `True`, lv of an spike train with less than two spikes
results in a `NaN` value and a warning is raised.
If `False`, an value error is raised.
Default: `True`
Returns
-------
Expand All @@ -222,7 +228,7 @@ def lv(v):
Raises
------
AttributeError :
ValueError :
If an empty list is specified, or if the sequence has less
than two entries, an AttributeError will be raised.
ValueError :
Expand All @@ -239,18 +245,30 @@ def lv(v):
"""
# convert to array, cast to float
v = np.asarray(v)
# ensure the input ia a vector
if len(v.shape) > 1:
raise ValueError("Input shape is larger than 1. Please provide "
"a vector as an input.")

# ensure we have enough entries
if v.size < 2:
raise AttributeError("Input size is too small. Please provide "
"an input with more than 1 entry.")
if with_nan:
warnings.warn("Input size is too small. Please provide "
"an input with more than 1 entry. lv returns 'NaN'"
"since the argument `with_nan` is True")
return np.NaN

else:
raise ValueError("Input size is too small. Please provide "
"an input with more than 1 entry. lv returned any"
"value since the argument `with_nan` is False")

# calculate LV and return result
# raise error if input is multi-dimensional
return 3. * np.mean(np.power(np.diff(v) / (v[:-1] + v[1:]), 2))


def cv2(v):
def cv2(v, with_nan=False):
"""
Calculate the measure of CV2 for a sequence of time intervals between
events.
Expand All @@ -274,17 +292,23 @@ def cv2(v):
v : quantity array, numpy array or list
Vector of consecutive time intervals
with_nans : bool, optional
If `True`, cv2 with less than two spikes results in a `NaN` value
and a warning is raised.
If `False`, an attribute error is raised.
Default: `True`
Returns
-------
cv2 : float
The CV2 of the inter-spike interval of the input sequence.
Raises
------
AttributeError :
ValueError :
If an empty list is specified, or if the sequence has less
than two entries, an AttributeError will be raised.
AttributeError :
ValueError :
Only vector inputs are supported. If a matrix is passed to the
function an AttributeError will be raised.
Expand All @@ -299,13 +323,20 @@ def cv2(v):

# ensure the input ia a vector
if len(v.shape) > 1:
raise AttributeError("Input shape is larger than 1. Please provide "
"a vector in input.")
raise ValueError("Input shape is larger than 1. Please provide "
"a vector as an input.")

# ensure we have enough entries
if v.size < 2:
raise AttributeError("Input size is too small. Please provide "
"an input with more than 1 entry.")
if with_nan:
warnings.warn("Input size is too small. Please provide"
"an input with more than 1 entry. cv2 returns `NaN`"
"since the argument `with_nan` is `True`")
return np.NaN
else:
raise ValueError("Input size is too small. Please provide "
"an input with more than 1 entry. cv2 returns any"
"value since the argument `with_nan` is `False`")

# calculate CV2 and return result
return 2. * np.mean(np.absolute(np.diff(v)) / (v[:-1] + v[1:]))
Expand All @@ -314,7 +345,7 @@ def cv2(v):
# sigma2kw and kw2sigma only needed for oldfct_instantaneous_rate!
# to finally be taken out of Elephant

def sigma2kw(form): # pragma: no cover
def sigma2kw(form): # pragma: no cover
warnings.simplefilter('always', DeprecationWarning)
warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
if form.upper() == 'BOX':
Expand All @@ -333,14 +364,14 @@ def sigma2kw(form): # pragma: no cover
return coeff


def kw2sigma(form): # pragma: no cover
def kw2sigma(form): # pragma: no cover
warnings.simplefilter('always', DeprecationWarning)
warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
return 1/sigma2kw(form)


# to finally be taken out of Elephant
def make_kernel(form, sigma, sampling_period, direction=1): # pragma: no cover
def make_kernel(form, sigma, sampling_period, direction=1): # pragma: no cover
"""
Creates kernel functions for convolution.
Expand Down Expand Up @@ -426,7 +457,7 @@ def make_kernel(form, sigma, sampling_period, direction=1): # pragma: no cover
warnings.warn("deprecated", DeprecationWarning, stacklevel=2)
forms_abbreviated = np.array(['BOX', 'TRI', 'GAU', 'EPA', 'EXP', 'ALP'])
forms_verbose = np.array(['boxcar', 'triangle', 'gaussian', 'epanechnikov',
'exponential', 'alpha'])
'exponential', 'alpha'])
if form in forms_verbose:
form = forms_abbreviated[forms_verbose == form][0]

Expand Down Expand Up @@ -506,8 +537,8 @@ def make_kernel(form, sigma, sampling_period, direction=1): # pragma: no cover

# to finally be taken out of Elephant
def oldfct_instantaneous_rate(spiketrain, sampling_period, form,
sigma='auto', t_start=None, t_stop=None,
acausal=True, trim=False): # pragma: no cover
sigma='auto', t_start=None, t_stop=None,
acausal=True, trim=False): # pragma: no cover
"""
Estimate instantaneous firing rate by kernel convolution.
Expand Down Expand Up @@ -644,15 +675,14 @@ def oldfct_instantaneous_rate(spiketrain, sampling_period, form,
t_stop = t_stop - ((kernel.size) - m_idx) * spiketrain.units

rate = neo.AnalogSignal(signal=r.reshape(r.size, 1),
sampling_period=sampling_period,
units=pq.Hz, t_start=t_start)
sampling_period=sampling_period,
units=pq.Hz, t_start=t_start)

return rate, sigma


def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
cutoff=5.0, t_start=None, t_stop=None, trim=False):

"""
Estimates instantaneous firing rate by kernel convolution.
Expand Down Expand Up @@ -734,7 +764,8 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
"""
# Merge spike trains if list of spike trains given:
if isinstance(spiketrain, list):
_check_consistency_of_spiketrainlist(spiketrain, t_start=t_start, t_stop=t_stop)
_check_consistency_of_spiketrainlist(
spiketrain, t_start=t_start, t_stop=t_stop)
if t_start is None:
t_start = spiketrain[0].t_start
if t_stop is None:
Expand Down Expand Up @@ -798,7 +829,8 @@ def instantaneous_rate(spiketrain, sampling_period, kernel='auto',
raise TypeError("trim must be bool!")

# main function:
units = pq.CompoundUnit("%s*s" % str(sampling_period.rescale('s').magnitude))
units = pq.CompoundUnit(
"%s*s" % str(sampling_period.rescale('s').magnitude))
spiketrain = spiketrain.rescale(units)
if t_start is None:
t_start = spiketrain.t_start
Expand Down Expand Up @@ -951,8 +983,8 @@ def time_histogram(spiketrains, binsize, t_start=None, t_stop=None,
raise ValueError('Parameter output is not valid.')

return neo.AnalogSignal(signal=bin_hist.reshape(bin_hist.size, 1),
sampling_period=binsize, units=bin_hist.units,
t_start=t_start)
sampling_period=binsize, units=bin_hist.units,
t_start=t_start)


def complexity_pdf(spiketrains, binsize):
Expand Down Expand Up @@ -1092,7 +1124,7 @@ def cost_function(x, N, w, dt):
yh = np.abs(fftkernel(x, w / dt)) # density
# formula for density
C = np.sum(yh ** 2) * dt - 2 * np.sum(yh * x) * \
dt + 2 / np.sqrt(2 * np.pi) / w / N
dt + 2 / np.sqrt(2 * np.pi) / w / N
C = C * N * N
# formula for rate
# C = dt*sum( yh.^2 - 2*yh.*y_hist + 2/sqrt(2*pi)/w*y_hist )
Expand Down Expand Up @@ -1125,7 +1157,7 @@ def sskernel(spiketimes, tin=None, w=None, bootstrap=False):
'C': cost functions of w,
'confb95': (lower bootstrap confidence level, upper bootstrap confidence level),
'yb': bootstrap samples.
If no optimal kernel could be found, all entries of the dictionary are set
to None.
Expand Down Expand Up @@ -1190,7 +1222,7 @@ def sskernel(spiketimes, tin=None, w=None, bootstrap=False):
f2, y2 = cost_function(yhist, N, logexp(c2), dt)
k = 0
while (np.abs(b - a) > (tolerance * (np.abs(c1) + np.abs(c2))))\
and (k < imax):
and (k < imax):
if f1 < f2:
b = c2
c2 = c1
Expand Down
17 changes: 11 additions & 6 deletions elephant/test/test_statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
import elephant.statistics as es
import elephant.kernels as kernels
import warnings
import math


class isi_TestCase(unittest.TestCase):
Expand Down Expand Up @@ -336,10 +337,14 @@ def test_lv_with_list(self):

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, [])
self.assertRaises(ValueError, es.lv, 1)
self.assertRaises(ValueError, es.lv, np.array([seq, seq]))


def test_2short_spike_train(self):
seq = [1]
self.assertTrue(math.isnan(es.lv(seq, with_nan=True)))


class CV2TestCase(unittest.TestCase):
def setUp(self):
Expand Down Expand Up @@ -370,9 +375,9 @@ def test_cv2_with_list(self):

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, es.cv2, [])
self.assertRaises(ValueError, es.cv2, 1)
self.assertRaises(ValueError, es.cv2, np.array([seq, seq]))


class RateEstimationTestCase(unittest.TestCase):
Expand Down

0 comments on commit 74ef4a2

Please sign in to comment.