switched from pyximport to distutils

3 files changed, 35 insertions, 123 deletions

diff --git a/pyspike/cython_distance.pyx b/pyspike/cython_distance.pyx
index 1a6d24a..6edcc01 100644
--- a/pyspike/cython_distance.pyx
+++ b/pyspike/cython_distance.pyx
@@ -3,7 +3,7 @@
 #cython: cdivision=True
 
 """
-cython_distances.py
+cython_distances.pyx
 
 cython implementation of the isi- and spike-distance
 
@@ -159,13 +159,14 @@ def spike_distance_cython(double[:] t1,
             # first calculate the previous interval end value
             dt_p1 = dt_f1 # the previous time now was the following time before
             s1 = dt_p1
-            s2 = (dt_p2*(t2[index2+1]-t1[index1]) + dt_f2*(t1[index1]-t2[index2])) / isi2
-            y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            s2 = (dt_p2*(t2[index2+1]-t1[index1]) + 
+                  dt_f2*(t1[index1]-t2[index2])) / isi2
+            y_ends[index-1] = (s1*isi2 + s2*isi1)/(0.5*(isi1+isi2)*(isi1+isi2))
             # now the next interval start value
             dt_f1 = get_min_dist_cython(t1[index1+1], t2, N2, index2)
             isi1 = t1[index1+1]-t1[index1]
             # s2 is the same as above, thus we can compute y2 immediately
-            y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            y_starts[index] = (s1*isi2 + s2*isi1)/(0.5*(isi1+isi2)*(isi1+isi2))
         elif t1[index1+1] > t2[index2+1]:
             index2 += 1
             if index2+1 >= N2:
@@ -173,15 +174,16 @@ def spike_distance_cython(double[:] t1,
             spike_events[index] = t2[index2]
             # first calculate the previous interval end value
             dt_p2 = dt_f2 # the previous time now was the following time before
-            s1 = (dt_p1*(t1[index1+1]-t2[index2]) + dt_f1*(t2[index2]-t1[index1])) / isi1
+            s1 = (dt_p1*(t1[index1+1]-t2[index2]) + 
+                  dt_f1*(t2[index2]-t1[index1])) / isi1
             s2 = dt_p2
-            y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)*(isi1+isi2))
             # now the next interval start value
             dt_f2 = get_min_dist_cython(t2[index2+1], t1, N1, index1)
             #s2 = dt_f2
             isi2 = t2[index2+1]-t2[index2]
             # s2 is the same as above, thus we can compute y2 immediately
-            y_starts[index] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+            y_starts[index] = (s1*isi2 + s2*isi1)/(0.5*(isi1+isi2)*(isi1+isi2))
         else: # t1[index1+1] == t2[index2+1] - generate only one event
             index1 += 1
             index2 += 1
@@ -204,7 +206,7 @@ def spike_distance_cython(double[:] t1,
     isi2 = max(t2[N2-1]-t2[N2-2], t2[N2-2]-t2[N2-3])
     s1 = dt_p1*(t1[N1-1]-t1[N1-2])/isi1
     s2 = dt_p2*(t2[N2-1]-t2[N2-2])/isi2
-    y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2)
+    y_ends[index-1] = (s1*isi2 + s2*isi1) / (0.5*(isi1+isi2)*(isi1+isi2))
     # use only the data added above 
     # could be less than original length due to equal spike times
     
diff --git a/pyspike/distances.py b/pyspike/distances.py
index 7a85471..79278b4 100644
--- a/pyspike/distances.py
+++ b/pyspike/distances.py
@@ -22,6 +22,7 @@ def add_auxiliary_spikes( spike_train, T_end , T_start=0.0):
     - T_start: start time of the observation interval (default 0.0)
     Returns:
     - spike train with additional spikes at T_start and T_end.
+
     """
     assert spike_train[0] >= T_start, \
            "Spike train has events before the given start time"
@@ -53,9 +54,7 @@ def isi_distance(spikes1, spikes2):
     assert spikes1[-1]==spikes2[-1], \
         "Given spike trains seems not to have auxiliary spikes!"
 
-    # compile and load cython implementation
-    import pyximport
-    pyximport.install(setup_args={'include_dirs': [np.get_include()]})
+    # cython implementation
     from cython_distance import isi_distance_cython
 
     times, values = isi_distance_cython(spikes1, spikes2)
@@ -81,9 +80,7 @@ def spike_distance(spikes1, spikes2):
     assert spikes1[-1]==spikes2[-1], \
         "Given spike trains seems not to have auxiliary spikes!"
 
-    # compile and load cython implementation
-    import pyximport
-    pyximport.install(setup_args={'include_dirs': [np.get_include()]})
+    # cython implementation
     from cython_distance import spike_distance_cython
 
     times, y_starts, y_ends = spike_distance_cython(spikes1, spikes2)
@@ -95,8 +92,8 @@ def spike_distance(spikes1, spikes2):
 # multi_distance
 ############################################################
 def multi_distance(spike_trains, pair_distance_func, indices=None):
-    """ Internal implementation detail, use isi_distance_multi or 
-    spike_distance_multi.
+    """ Internal implementation detail, don't call this function directly,
+    use isi_distance_multi or spike_distance_multi instead.
 
     Computes the multi-variate distance for a set of spike-trains using the
     pair_dist_func to compute pair-wise distances. That is it computes the 
diff --git a/pyspike/function.py b/pyspike/function.py
index 26ca4b2..5444c36 100644
--- a/pyspike/function.py
+++ b/pyspike/function.py
@@ -82,46 +82,15 @@ class PieceWiseConstFunc:
         """
         assert self.x[0] == f.x[0], "The functions have different intervals"
         assert self.x[-1] == f.x[-1], "The functions have different intervals"
-        x_new = np.empty(len(self.x) + len(f.x))
-        y_new = np.empty(len(x_new)-1)
-        x_new[0] = self.x[0]
-        y_new[0] = self.y[0] + f.y[0]
-        index1 = 0
-        index2 = 0
-        index = 0
-        while (index1+1 < len(self.y)) and (index2+1 < len(f.y)):
-            index += 1
-            # print(index1+1, self.x[index1+1], self.y[index1+1], x_new[index])
-            if self.x[index1+1] < f.x[index2+1]:
-                index1 += 1
-                x_new[index] = self.x[index1]
-            elif self.x[index1+1] > f.x[index2+1]:
-                index2 += 1
-                x_new[index] = f.x[index2]
-            else: # self.x[index1+1] == f.x[index2+1]:
-                index1 += 1
-                index2 += 1
-                x_new[index] = self.x[index1]
-            y_new[index] = self.y[index1] + f.y[index2]
-        # one array reached the end -> copy the contents of the other to the end
-        if index1+1 < len(self.y):
-            x_new[index+1:index+1+len(self.x)-index1-1] = self.x[index1+1:]
-            y_new[index+1:index+1+len(self.y)-index1-1] = self.y[index1+1:] + \
-                                                          f.y[-1]
-            index += len(self.x)-index1-2
-        elif index2+1 < len(f.y):
-            x_new[index+1:index+1+len(f.x)-index2-1] = f.x[index2+1:]
-            y_new[index+1:index+1+len(f.y)-index2-1] = f.y[index2+1:] + \
-                                                       self.y[-1]
-            index += len(f.x)-index2-2
-        else: # both arrays reached the end simultaneously
-            # only the last x-value missing
-            x_new[index+1] = self.x[-1]
-        # the last value is again the end of the interval
-        # x_new[index+1] = self.x[-1]
-        # only use the data that was actually filled
-        self.x = x_new[:index+2]
-        self.y = y_new[:index+1]
+        
+        # python implementation
+        # from python_backend import add_piece_wise_const_python
+        # self.x, self.y = add_piece_wise_const_python(self.x, self.y, f.x, f.y)
+
+        # cython version
+        from cython_add import add_piece_wise_const_cython
+        self.x, self.y = add_piece_wise_const_cython(self.x, self.y, f.x, f.y)
+
 
     def mul_scalar(self, fac):
         """ Multiplies the function with a scalar value
@@ -204,73 +173,17 @@ class PieceWiseLinFunc:
         """
         assert self.x[0] == f.x[0], "The functions have different intervals"
         assert self.x[-1] == f.x[-1], "The functions have different intervals"
-        x_new = np.empty(len(self.x) + len(f.x))
-        y1_new = np.empty(len(x_new)-1)
-        y2_new = np.empty_like(y1_new)
-        x_new[0] = self.x[0]
-        y1_new[0] = self.y1[0] + f.y1[0]
-        index1 = 0 # index for self
-        index2 = 0 # index for f
-        index = 0  # index for new
-        while (index1+1 < len(self.y1)) and (index2+1 < len(f.y1)):
-            # print(index1+1, self.x[index1+1], self.y[index1+1], x_new[index])
-            if self.x[index1+1] < f.x[index2+1]:
-                # first compute the end value of the previous interval
-                # linear interpolation of the interval
-                y = f.y1[index2] + (f.y2[index2]-f.y1[index2]) * \
-                    (self.x[index1+1]-f.x[index2]) / (f.x[index2+1]-f.x[index2])
-                y2_new[index] = self.y2[index1] + y
-                index1 += 1
-                index += 1
-                x_new[index] = self.x[index1]
-                # and the starting value for the next interval
-                y1_new[index] = self.y1[index1] + y
-            elif self.x[index1+1] > f.x[index2+1]:
-                # first compute the end value of the previous interval
-                # linear interpolation of the interval
-                y = self.y1[index1] + (self.y2[index1]-self.y1[index1]) * \
-                    (f.x[index2+1]-self.x[index1]) / \
-                    (self.x[index1+1]-self.x[index1])
-                y2_new[index] = f.y2[index2] + y
-                index2 += 1
-                index += 1
-                x_new[index] = f.x[index2]
-                # and the starting value for the next interval
-                y1_new[index] = f.y1[index2] + y
-            else: # self.x[index1+1] == f.x[index2+1]:
-                y2_new[index] = self.y2[index1] + f.y2[index2]
-                index1 += 1
-                index2 += 1
-                index += 1
-                x_new[index] = self.x[index1]
-                y1_new[index] = self.y1[index1] + f.y1[index2]
-        # one array reached the end -> copy the contents of the other to the end
-        if index1+1 < len(self.y1):
-            # compute the linear interpolations values
-            y = f.y1[index2] + (f.y2[index2]-f.y1[index2]) * \
-                (self.x[index1+1:-1]-f.x[index2]) / (f.x[index2+1]-f.x[index2])
-            x_new[index+1:index+1+len(self.x)-index1-1] = self.x[index1+1:]
-            y1_new[index+1:index+1+len(self.y1)-index1-1] = self.y1[index1+1:]+y
-            y2_new[index:index+len(self.y2)-index1-1] = self.y2[index1:-1] + y
-            index += len(self.x)-index1-2
-        elif index2+1 < len(f.y1):
-            # compute the linear interpolations values
-            y = self.y1[index1] + (self.y2[index1]-self.y1[index1]) * \
-                (f.x[index2+1:-1]-self.x[index1]) / \
-                (self.x[index1+1]-self.x[index1])
-            x_new[index+1:index+1+len(f.x)-index2-1] = f.x[index2+1:]
-            y1_new[index+1:index+1+len(f.y1)-index2-1] = f.y1[index2+1:] + y
-            y2_new[index:index+len(f.y2)-index2-1] = f.y2[index2:-1] + y
-            index += len(f.x)-index2-2
-        else: # both arrays reached the end simultaneously
-            # only the last x-value missing
-            x_new[index+1] = self.x[-1]
-        # finally, the end value for the last interval
-        y2_new[index] = self.y2[-1]+f.y2[-1]
-        # only use the data that was actually filled
-        self.x = x_new[:index+2]
-        self.y1 = y1_new[:index+1]
-        self.y2 = y2_new[:index+1]
+
+        # python implementation
+        # from python_backend import add_piece_wise_lin_python
+        # self.x, self.y1, self.y2 = add_piece_wise_lin_python( 
+        #     self.x, self.y1, self.y2, f.x, f.y1, f.y2)
+
+        # cython version
+        from cython_add import add_piece_wise_lin_cython
+        self.x, self.y1, self.y2 = add_piece_wise_lin_cython(
+            self.x, self.y1, self.y2, f.x, f.y1, f.y2)
+
 
     def mul_scalar(self, fac):
         """ Multiplies the function with a scalar value