""" python_backend.py Collection of python functions that can be used instead of the cython implementation. Copyright 2014, Mario Mulansky Distributed under the BSD License """ import numpy as np from pyspike import PieceWiseConstFunc, PieceWiseLinFunc ############################################################ # isi_distance_python ############################################################ def isi_distance_python(s1, s2): """ Plain Python implementation of the isi distance. """ # compute the interspike interval nu1 = s1[1:] - s1[:-1] nu2 = s2[1:] - s2[:-1] # compute the isi-distance spike_events = np.empty(len(nu1) + len(nu2)) spike_events[0] = s1[0] # the values have one entry less - the number of intervals between events isi_values = np.empty(len(spike_events) - 1) # add the distance of the first events # isi_values[0] = nu1[0]/nu2[0] - 1.0 if nu1[0] <= nu2[0] \ # else 1.0 - nu2[0]/nu1[0] isi_values[0] = (nu1[0] - nu2[0]) / max(nu1[0], nu2[0]) index1 = 0 index2 = 0 index = 1 while True: # check which spike is next - from s1 or s2 if s1[index1+1] < s2[index2+1]: index1 += 1 # break condition relies on existence of spikes at T_end if index1 >= len(nu1): break spike_events[index] = s1[index1] elif s1[index1+1] > s2[index2+1]: index2 += 1 if index2 >= len(nu2): break spike_events[index] = s2[index2] else: # s1[index1 + 1] == s2[index2 + 1] index1 += 1 index2 += 1 if (index1 >= len(nu1)) or (index2 >= len(nu2)): break spike_events[index] = s1[index1] # compute the corresponding isi-distance isi_values[index] = (nu1[index1] - nu2[index2]) / \ max(nu1[index1], nu2[index2]) index += 1 # the last event is the interval end spike_events[index] = s1[-1] # use only the data added above # could be less than original length due to equal spike times return PieceWiseConstFunc(spike_events[:index + 1], isi_values[:index]) ############################################################ # get_min_dist ############################################################ def get_min_dist(spike_time, spike_train, start_index=0): """ Returns the minimal distance |spike_time - spike_train[i]| with i>=start_index. """ d = abs(spike_time - spike_train[start_index]) start_index += 1 while start_index < len(spike_train): d_temp = abs(spike_time - spike_train[start_index]) if d_temp > d: break else: d = d_temp start_index += 1 return d ############################################################ # spike_distance_python ############################################################ def spike_distance_python(spikes1, spikes2): """ Computes the instantaneous spike-distance S_spike (t) of the two given spike trains. The spike trains are expected to have auxiliary spikes at the beginning and end of the interval. Use the function add_auxiliary_spikes to add those spikes to the spike train. Args: - spikes1, spikes2: ordered arrays of spike times with auxiliary spikes. Returns: - PieceWiseLinFunc describing the spike-distance. """ # check for auxiliary spikes - first and last spikes should be identical assert spikes1[0] == spikes2[0], \ "Given spike trains seems not to have auxiliary spikes!" assert spikes1[-1] == spikes2[-1], \ "Given spike trains seems not to have auxiliary spikes!" # shorter variables t1 = spikes1 t2 = spikes2 spike_events = np.empty(len(t1) + len(t2) - 2) spike_events[0] = t1[0] y_starts = np.empty(len(spike_events) - 1) y_ends = np.empty(len(spike_events) - 1) index1 = 0 index2 = 0 index = 1 dt_p1 = 0.0 dt_f1 = get_min_dist(t1[1], t2, 0) dt_p2 = 0.0 dt_f2 = get_min_dist(t2[1], t1, 0) isi1 = max(t1[1]-t1[0], t1[2]-t1[1]) isi2 = max(t2[1]-t2[0], t2[2]-t2[1]) s1 = dt_f1*(t1[1]-t1[0])/isi1 s2 = dt_f2*(t2[1]-t2[0])/isi2 y_starts[0] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) while True: # print(index, index1, index2) if t1[index1+1] < t2[index2+1]: index1 += 1 # break condition relies on existence of spikes at T_end if index1+1 >= len(t1): break spike_events[index] = t1[index1] # first calculate the previous interval end value dt_p1 = dt_f1 # the previous time 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) # now the next interval start value dt_f1 = get_min_dist(t1[index1+1], t2, 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) elif t1[index1+1] > t2[index2+1]: index2 += 1 if index2+1 >= len(t2): break spike_events[index] = t2[index2] # first calculate the previous interval end value dt_p2 = dt_f2 # the previous time was the following time before 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) # now the next interval start value dt_f2 = get_min_dist(t2[index2+1], t1, 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) else: # t1[index1+1] == t2[index2+1] - generate only one event index1 += 1 index2 += 1 if (index1+1 >= len(t1)) or (index2+1 >= len(t2)): break assert dt_f2 == 0.0 assert dt_f1 == 0.0 spike_events[index] = t1[index1] y_ends[index-1] = 0.0 y_starts[index] = 0.0 dt_p1 = 0.0 dt_p2 = 0.0 dt_f1 = get_min_dist(t1[index1+1], t2, index2) dt_f2 = get_min_dist(t2[index2+1], t1, index1) isi1 = t1[index1+1]-t1[index1] isi2 = t2[index2+1]-t2[index2] index += 1 # the last event is the interval end spike_events[index] = t1[-1] # the ending value of the last interval isi1 = max(t1[-1]-t1[-2], t1[-2]-t1[-3]) isi2 = max(t2[-1]-t2[-2], t2[-2]-t2[-3]) s1 = dt_p1*(t1[-1]-t1[-2])/isi1 s2 = dt_p2*(t2[-1]-t2[-2])/isi2 y_ends[index-1] = (s1*isi2 + s2*isi1) / ((isi1+isi2)**2/2) # use only the data added above # could be less than original length due to equal spike times return PieceWiseLinFunc(spike_events[:index+1], y_starts[:index], y_ends[:index]) ############################################################ # add_piece_wise_const_python ############################################################ def add_piece_wise_const_python(x1, y1, x2, y2): x_new = np.empty(len(x1) + len(x2)) y_new = np.empty(len(x_new)-1) x_new[0] = x1[0] y_new[0] = y1[0] + y2[0] index1 = 0 index2 = 0 index = 0 while (index1+1 < len(y1)) and (index2+1 < len(y2)): index += 1 # print(index1+1, x1[index1+1], y1[index1+1], x_new[index]) if x1[index1+1] < x2[index2+1]: index1 += 1 x_new[index] = x1[index1] elif x1[index1+1] > x2[index2+1]: index2 += 1 x_new[index] = x2[index2] else: # x1[index1+1] == x2[index2+1]: index1 += 1 index2 += 1 x_new[index] = x1[index1] y_new[index] = y1[index1] + y2[index2] # one array reached the end -> copy the contents of the other to the end if index1+1 < len(y1): x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] y_new[index+1:index+1+len(y1)-index1-1] = y1[index1+1:] + y2[-1] index += len(x1)-index1-2 elif index2+1 < len(y2): x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] y_new[index+1:index+1+len(y2)-index2-1] = y2[index2+1:] + y1[-1] index += len(x2)-index2-2 else: # both arrays reached the end simultaneously # only the last x-value missing x_new[index+1] = x1[-1] # the last value is again the end of the interval # x_new[index+1] = x1[-1] # only use the data that was actually filled return x_new[:index+2], y_new[:index+1] ############################################################ # add_piece_lin_const_python ############################################################ def add_piece_wise_lin_python(x1, y11, y12, x2, y21, y22): x_new = np.empty(len(x1) + len(x2)) y1_new = np.empty(len(x_new)-1) y2_new = np.empty_like(y1_new) x_new[0] = x1[0] y1_new[0] = y11[0] + y21[0] index1 = 0 # index for self index2 = 0 # index for f index = 0 # index for new while (index1+1 < len(y11)) and (index2+1 < len(y21)): # print(index1+1, x1[index1+1], self.y[index1+1], x_new[index]) if x1[index1+1] < x2[index2+1]: # first compute the end value of the previous interval # linear interpolation of the interval y = y21[index2] + (y22[index2]-y21[index2]) * \ (x1[index1+1]-x2[index2]) / (x2[index2+1]-x2[index2]) y2_new[index] = y12[index1] + y index1 += 1 index += 1 x_new[index] = x1[index1] # and the starting value for the next interval y1_new[index] = y11[index1] + y elif x1[index1+1] > x2[index2+1]: # first compute the end value of the previous interval # linear interpolation of the interval y = y11[index1] + (y12[index1]-y11[index1]) * \ (x2[index2+1]-x1[index1]) / \ (x1[index1+1]-x1[index1]) y2_new[index] = y22[index2] + y index2 += 1 index += 1 x_new[index] = x2[index2] # and the starting value for the next interval y1_new[index] = y21[index2] + y else: # x1[index1+1] == x2[index2+1]: y2_new[index] = y12[index1] + y22[index2] index1 += 1 index2 += 1 index += 1 x_new[index] = x1[index1] y1_new[index] = y11[index1] + y21[index2] # one array reached the end -> copy the contents of the other to the end if index1+1 < len(y11): # compute the linear interpolations values y = y21[index2] + (y22[index2]-y21[index2]) * \ (x1[index1+1:-1]-x2[index2]) / (x2[index2+1]-x2[index2]) x_new[index+1:index+1+len(x1)-index1-1] = x1[index1+1:] y1_new[index+1:index+1+len(y11)-index1-1] = y11[index1+1:]+y y2_new[index:index+len(y12)-index1-1] = y12[index1:-1] + y index += len(x1)-index1-2 elif index2+1 < len(y21): # compute the linear interpolations values y = y11[index1] + (y12[index1]-y11[index1]) * \ (x2[index2+1:-1]-x1[index1]) / \ (x1[index1+1]-x1[index1]) x_new[index+1:index+1+len(x2)-index2-1] = x2[index2+1:] y1_new[index+1:index+1+len(y21)-index2-1] = y21[index2+1:] + y y2_new[index:index+len(y22)-index2-1] = y22[index2:-1] + y index += len(x2)-index2-2 else: # both arrays reached the end simultaneously # only the last x-value missing x_new[index+1] = x1[-1] # finally, the end value for the last interval y2_new[index] = y12[-1]+y22[-1] # only use the data that was actually filled return x_new[:index+2], y1_new[:index+1], y2_new[:index+1]