some refactoring and bug fix for issue #1

.gitignore

@@ -58,3 +58,6 @@ docs/_build/
 .swp
 .swo
 *~
+
+# Ide configuration
+.idea

ledapy/analyse.py

@@ -90,24 +90,23 @@ def succnz(data, crit, fac, sr):
     """
     Original location: analyze/deconvolution/succnz.m
     """
-    nzL = npa([])
-    cntr = 0
     n = len(data)
-    for i in range(n):
-        if abs(data[i]) > crit:
-            cntr = cntr + 1
-        else:
-            if cntr > 0:
-                nzL = np.hstack((nzL, cntr))
-                cntr = 0
-    if cntr > 0:
-        nzL = np.hstack((nzL, cntr))
-    if not len(nzL) == 0:
-        snz = np.sum(pow((nzL / sr), fac)) / (n / sr)
-    else:
-        snz = 0
-    return snz
 
+    abovecrit = np.abs(data) > crit
+    nzidx = np.flatnonzero(np.diff(abovecrit)) + 1
+
+    if len(nzidx) == 0:
+        return 0
+
+    if abovecrit[0] == 1:
+        nzidx = np.hstack((1, nzidx))
+
+    if abovecrit[-1] == 1:
+        nzidx = np.hstack((nzidx, n+1))
+
+    nzL = nzidx[np.arange(1, len(nzidx), 2)] - nzidx[np.arange(0, len(nzidx), 2)]
+
+    return np.sum(pow(nzL / sr, fac) / (n / sr))
 
 def cgd_linesearch(x, error0, direction, error_fcn, h):
     """
@@ -348,7 +347,7 @@ def sdeco_interimpulsefit(driver, kernel, minL, maxL):
             t_idx = np.logical_and(iif_t > groundtime[i] - tonicGridSize, iif_t < t[-1] - 1)
             grid_idx = np.logical_and(t > groundtime[i] - tonicGridSize, t < t[-1] - 1)
         else:
-            t_idx = np.logical_and(iif_t > groundtime[i] - tonicGridSize, iif_t <= groundtime[i] + tonicGridSize / 2)
+            t_idx = np.logical_and(iif_t > groundtime[i] - tonicGridSize / 2, iif_t <= groundtime[i] + tonicGridSize / 2)
             grid_idx = np.logical_and(t > groundtime[i] - tonicGridSize / 2, t <= groundtime[i] + tonicGridSize / 2)
         # Estimate groundlevel at groundtime
         if len(np.flatnonzero(t_idx)) > 2:

ledapy/deconvolution.py

@@ -62,17 +62,16 @@ def sdeconv_analysis(x, estim_tonic=1):
     """
     Original location: deconvolution/sdeconv_analysis.m
     """
-    newx = npa(x)
     # Check Limits
-    newx[0] = utils.withinlimits(newx[0], leda2.settings.tauMin, 10)
-    newx[1] = utils.withinlimits(newx[1], leda2.settings.tauMin, 20)
+    x[0] = utils.withinlimits(x[0], leda2.settings.tauMin, 10)
+    x[1] = utils.withinlimits(x[1], leda2.settings.tauMin, 20)
 
-    if newx[1] < newx[0]:  # tau1 < tau2
-        newx = newx[::-1]
-    if abs(newx[0] - newx[1]) < leda2.settings.tauMinDiff:
-        newx[1] = newx[1] + leda2.settings.tauMinDiff
+    if x[1] < x[0]:  # tau1 < tau2
+        x = x[::-1]
+    if abs(x[0] - x[1]) < leda2.settings.tauMinDiff:
+        x[1] = x[1] + leda2.settings.tauMinDiff
 
-    tau = npa(newx[0:2])
+    tau = npa(x[0:2])
 
     data = npa(leda2.analysis0.target.d)
     t = npa(leda2.analysis0.target.t)
@@ -91,7 +90,7 @@ def sdeconv_analysis(x, estim_tonic=1):
     bg = analyse.bateman_gauss(tb, 5, 1, 2, 40, .4)
     idx = np.argmax(bg)
 
-    prefix = bg[:idx + 2] / bg[idx + 2] * d[0]  # +10
+    prefix = bg[:idx + 1] / bg[idx + 1] * d[0]  # +10
     prefix = npa(prefix)
     prefix_idx = np.flatnonzero(prefix)
     prefix = prefix[prefix_idx]
@@ -144,7 +143,7 @@ def sdeconv_analysis(x, estim_tonic=1):
     err_MSE = analyse.fiterror(data, tonicData + phasicData, 0)
     err_RMSE = np.sqrt(err_MSE)
     err_chi2 = err_RMSE / leda2.data.conductance_error
-    err1d = deverror(phasicDriver, [0, .2])
+    err1d = deverror(phasicDriver, .2)
     err1s = analyse.succnz(phasicDriver, np.maximum(.01, np.max(phasicDriver) / 20), 2, sr)
     phasicDriverNeg = npa(phasicDriver)
     phasicDriverNeg[np.flatnonzero(phasicDriverNeg > 0)] = 0
@@ -174,18 +173,16 @@ def sdeconv_analysis(x, estim_tonic=1):
     leda2.analysis0.error.discreteness = np.hstack((err_discreteness, 0))
     leda2.analysis0.error.negativity = err_negativity
     leda2.analysis0.error.compound = err
-    return (err, newx)
+    return (err, x)
 
 
 def deverror(v, elim):
     if not isinstance(v, np.ndarray):
         raise ValueError('v is not a numpy array')
 
-    err = 0.0
-    err += np.sum(v[np.logical_and(v > elim[0], v < 0)] / elim[0])
-    err += np.sum(v[np.logical_and(v < elim[1], v > 0)] / elim[1])
-    err += np.prod(np.size(np.flatnonzero(np.logical_or(v <= elim[0], v >= elim[1]))))
-    return err / np.prod(np.size(v))
+    idx = np.logical_and(v > 0, v < elim)
+    err = 1 + ( np.sum(v[idx]) / elim - np.sum(idx) ) / len(v)
+    return err
 
 
 def deconv_optimize(x0, nr_iv):

ledapy/utils.py

@@ -158,8 +158,8 @@ def smooth(data, winwidth_in):
     if not len(data.shape) == 1:
         raise ValueError('data is not a vector. Shape: ' + str(data.shape))
     paddata = npa(np.hstack((data[0], data, data[-1])))  # pad to remove border errors
-    winwidth = int(np.floor(winwidth_in // 2) * 2)  # force even winsize for odd window
-    window = norm.pdf(np.arange(0, winwidth + 1), winwidth // 2, winwidth // 8)
+    winwidth = int(np.floor(winwidth_in / 2) * 2)  # force even winsize for odd window
+    window = norm.pdf(np.arange(0, winwidth + 1), winwidth / 2 + 1, winwidth / 8)
     window = window / np.sum(window)  # normalize window
 
     data_ext = np.hstack(((np.zeros(winwidth // 2) + 1) * paddata[0],

setup.py

@@ -2,7 +2,7 @@
 
 setuptools.setup(
     name="ledapy",
-    version="1.0",
+    version="1.1",
     author="Marco Filetti",
     description="Ledapy is a minimal port of Ledalab (www.ledalab.de) for Python",
     url="https://github.com/HIIT/ledapy",