Spike Order support (#39)0.6.0

* reorganized directionality module

* further refactoring of directionality

* completed python directionality backend

* added SPIKE-Sync based filtering

new function filter_by_spike_sync removes spikes that have a multi-variate
Spike Sync value below some threshold

not yet fully tested, python backend missing.

* spike sync filtering, cython sim ann

Added function for filtering out events based on a threshold for the spike
sync values. Usefull for focusing on synchronous events during directionality
analysis.
Also added cython version of simulated annealing for performance.

* added coincidence single profile to python backend

missing function in python backend added, identified and fixed a bug in the
implementation as well

* updated test case to new spike sync behavior

* python3 fixes

* another python3 fix

* reorganized directionality module

* further refactoring of directionality

* completed python directionality backend

* added SPIKE-Sync based filtering

new function filter_by_spike_sync removes spikes that have a multi-variate
Spike Sync value below some threshold

not yet fully tested, python backend missing.

* spike sync filtering, cython sim ann

Added function for filtering out events based on a threshold for the spike
sync values. Usefull for focusing on synchronous events during directionality
analysis.
Also added cython version of simulated annealing for performance.

* added coincidence single profile to python backend

missing function in python backend added, identified and fixed a bug in the
implementation as well

* updated test case to new spike sync behavior

* python3 fixes

* another python3 fix

* Fix absolute imports in directionality measures

* remove commented code

* Add directionality to docs, bump version

* Clean up directionality module, add doxy.

* Remove debug print from tests

* Fix bug in calling Python backend

* Fix incorrect integrals in PieceWiseConstFunc (#36)

* Add (some currently failing) tests for PieceWiseConstFunc.integral

* Fix implementation of PieceWiseConstFunc.integral

Just by adding a special condition for when we are only taking an
integral "between" two edges of a PieceWiseConstFunc
All tests now pass.
Fixes #33.

* Add PieceWiseConstFunc.integral tests for ValueError

* Add testing bounds of integral
* Raise ValueError in function implementation

* Fix incorrect integrals in PieceWiseLinFunc (#38)

Integrals of piece-wise linear functions were incorrect if the
requested interval lies completely between two support points.
This has been fixed, and a unit test exercising this behavior
was added.

Fixes #38

* Add Spike Order example and Tutorial section

Adds an example computing spike order profile and the optimal
spike train order. Also adds a section on spike train order to the
tutorial.

1 files changed, 82 insertions, 0 deletions

diff --git a/pyspike/cython/cython_simulated_annealing.pyx b/pyspike/cython/cython_simulated_annealing.pyx
new file mode 100644
index 0000000..be9423c
--- /dev/null
+++ b/pyspike/cython/cython_simulated_annealing.pyx
@@ -0,0 +1,82 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: cdivision=True
+
+"""
+cython_simulated_annealing.pyx
+
+cython implementation of a simulated annealing algorithm to find the optimal
+spike train order
+
+Note: using cython memoryviews (e.g. double[:]) instead of ndarray objects
+improves the performance of spike_distance by a factor of 10!
+
+Copyright 2015, Mario Mulansky <mario.mulansky@gmx.net>
+
+Distributed under the BSD License
+
+"""
+
+"""
+To test whether things can be optimized: remove all yellow stuff
+in the html output::
+
+  cython -a cython_simulated_annealing.pyx
+
+which gives:
+
+  cython_simulated_annealing.html
+
+"""
+
+import numpy as np
+cimport numpy as np
+
+from libc.math cimport exp
+from libc.math cimport fmod
+from libc.stdlib cimport rand
+from libc.stdlib cimport RAND_MAX
+
+DTYPE = np.float
+ctypedef np.float_t DTYPE_t
+
+
+def sim_ann_cython(double[:, :] D, double T_start, double T_end, double alpha):
+
+    cdef long N = len(D)
+    cdef double A = np.sum(np.triu(D, 0))
+    cdef long[:] p = np.arange(N)
+    cdef double T = T_start
+    cdef long iterations
+    cdef long succ_iter
+    cdef long total_iter = 0
+    cdef double delta_A
+    cdef long ind1
+    cdef long ind2
+
+    while T > T_end:
+        iterations = 0
+        succ_iter = 0
+        # equilibrate for 100*N steps or 10*N successful steps
+        while iterations < 100*N and succ_iter < 10*N:
+            # exchange two rows and cols
+            # ind1 = np.random.randint(N-1)
+            ind1 = rand() % (N-1)
+            if ind1 < N-1:
+                ind2 = ind1+1
+            else:  # this can never happen!
+                ind2 = 0
+            delta_A = -2*D[p[ind1], p[ind2]]
+            if delta_A > 0.0 or exp(delta_A/T) > ((1.0*rand()) / RAND_MAX):
+                # swap indices
+                p[ind1], p[ind2] = p[ind2], p[ind1]
+                A += delta_A
+                succ_iter += 1
+            iterations += 1
+        total_iter += iterations
+        T *= alpha   # cool down
+        if succ_iter == 0:
+            # no successful step -> we believe we have converged
+            break
+
+    return p, A, total_iter