Tell me more ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.
up vote 2 down vote favorite

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)
share|improve this question

1 Answer

up vote 1 down vote

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.