Specialized version of the cross correlation function

I am writing a specialized version of the cross correlation function as used in neuroscience. The function below is supposed to take a time series data and ask how many of its values fall in specified bins. My function xcorr works but is horrifically slow even. A test data set with 1000 points (and so 0.5*(1000*999) intervals) distributed over 400 bins takes almost ten minutes.

Bottleneck

The bottleneck is in the line counts = array([sum .... I assume it is there because each iteration of the foreach loop searches through the entire vector diffs, which is of length len(first)**2.

def xcorr(first,second,dt=0.001,window=0.2,savename=None):
length = min(len(first),len(second))
diffs = array([first[i]-second[j] for i in xrange(length) for j in xrange(length)])
bins = arange(-(int(window/dt)),int(window/dt))
counts = array[sum((diffs>((i-.5)*dt)) & (diffs<((i+.5)*dt))) for i in bins]
counts -= len(first)**2*dt/float(max(first[-1],second[-1])) #Normalization
if savename:
cPickle.dump((bins*dt,counts),open(savename,'wb'))
return (bins,counts)


2 Answers

Here's a crack at my own question that uses built-in functions from NumPy

Update

The function chain bincount(digitize(data,bins)) is equivalent to histogram, which allows for a more succinct expression.

 def xcorr(first,second,dt=0.001,window=0.2,savename=None):
length = min(len(first),len(second))
diffs = array([first[i]-second[j] for i in xrange(length) for j in xrange(length)])
counts,edges = histogram(diffs,bins=arange(-window,window,dt))
counts -= length**2*dt/float(first[length])
if savename:
cPickle.dump((counts,edges[:-1]),open(savename,'wb'))
return (counts,edges[:-1])


The indexing in the return statement comes because I want to ignore the edge bins.

Going to try to help you do this in a faster fashion rather than write the code for you entirely:

import pandas as pd
data = pd.DataFrame(data)


Now you can easily do your diffs without those slow loops - note my comments should help but basically you are taking advantage of Pandas very fast and simple calculation on the time series; .shift(1) if 1st row is the OLDEST date (example is for the NEWEST date in 1st row .shift(-1)). Likewise you can shift(2) if the calculation is 2 rows away. Modify to use your algorithm:

diffs = (data[0:] - data[0:].shift(-1))/data[0:]


Okay back to your bottleneck issue. Can you tell in your question after your NumPy modification which line is slow? Just counting bins?