Python: sliding window of variable width

Question

I'm writing a program in Python that's processing some data generated during experiments, and it needs to estimate the slope of the data. I've written a piece of code that does this quite nicely, but it's horribly slow (and I'm not very patient). Let me explain how this code works:

1) It grabs a small piece of data of size dx (starting with 3 datapoints)

2) It evaluates whether the difference (i.e. |y(x+dx)-y(x-dx)| ) is larger than a certain minimum value (40x std. dev. of noise)

3) If the difference is large enough, it will calculate the slope using OLS regression. If the difference is too small, it will increase dx and redo the loop with this new dx

4) This continues for all the datapoints

[See updated code further down]

For a datasize of about 100k measurements, this takes about 40 minutes, whereas the rest of the program (it does more processing than just this bit) takes about 10 seconds. I am certain there is a much more efficient way of doing these operations, could you guys please help me out?

Thanks

EDIT:

Ok, so I've got the problem solved by using only binary searches, limiting the number of allowed steps by 200. I thank everyone for their input and I selected the answer that helped me most.

FINAL UPDATED CODE:

def slope(self, data, time):
    (wave1, wave2) = wt.dwt(data, "db3")
    std = 2*np.std(wave2)
    e = std/0.05
    de = 5*std
    N = len(data)
    slopes = np.ones(shape=(N,))
    data2 = np.concatenate((-data[::-1]+2*data[0], data, -data[::-1]+2*data[N-1]))
    time2 = np.concatenate((-time[::-1]+2*time[0], time, -time[::-1]+2*time[N-1]))
    for n in xrange(N+1, 2*N):     
        left = N+1
        right = 2*N
        for i in xrange(200):
            mid = int(0.5*(left+right))
            diff = np.abs(data2[n-mid+N]-data2[n+mid-N])
            if diff >= e:
                if diff < e + de:  
                    break
                right = mid - 1
                continue
            left = mid + 1
        leftlim = n - mid + N
        rightlim = n + mid - N
        y = data2[leftlim:rightlim:int(0.05*(rightlim-leftlim)+1)]
        x = time2[leftlim:rightlim:int(0.05*(rightlim-leftlim)+1)]
        xavg = np.average(x)
        yavg = np.average(y)
        xlen = len(x)
        slopes[n-N] = (np.dot(x,y)-xavg*yavg*xlen)/(np.dot(x,x)-xavg*xavg*xlen)
    return np.array(slopes) 

Answer 1

Your comments suggest that you need to find a better method to estimate i_k+1 given i_k. No knowledge of values in data would yield to the naive algorithm:

At each iteration for n, leave i at previous value, and see if the abs(data[start]-data[end]) value is less than e. If it is, leave i at its previous value, and find your new one by incrementing it by 1 as you do now. If it is greater, or equal, do a binary search on i to find the appropriate value. You can possibly do a binary search forwards, but finding a good candidate upper limit without knowledge of data can prove to be difficult. This algorithm won't perform worse than your current estimation method.

If you know that data is kind of smooth (no sudden jumps, and hence a smooth plot for all i values) and monotonically increasing, you can replace the binary search with a search backwards by decrementing its value by 1 instead.

Answer 2

How to optimize this will depend on some properties of your data, but here are some ideas:

1. Have you tried profiling the code? Using one of the Python profilers can give you some useful information about what's taking the most time. Often, a piece of code you've just written will have one biggest bottleneck, and it's not always obvious which piece it is; profiling lets you figure that out and attack the main bottleneck first.

2. Do you know what typical values of i are? If you have some idea, you can speed things up by starting with i greater than 0 (as @vhallac noted), or by increasing i by larger amounts — if you often see big values for i, increase i by 2 or 3 at a time; if the distribution of is has a long tail, try doubling it each time; etc.

3. Do you need all the data when doing the least squares regression? If that function call is the bottleneck, you may be able to speed it up by using only some of the data in the range. Suppose, for instance, that at a particular point, you need i to be 200 to see a large enough (above-noise) change in the data. But you may not need all 400 points to get a good estimate of the slope — just using 10 or 20 points, evenly spaced in the start:end range, may be sufficient, and might speed up the code a lot.

Answer 3

I work with Python for similar analyses, and have a few suggestions to make. I didn't look at the details of your code, just to your problem statement:

1) It grabs a small piece of data of size dx (starting with 3 datapoints)

2) It evaluates whether the difference (i.e. |y(x+dx)-y(x-dx)| ) is larger than a certain minimum value (40x std. dev. of noise)

3) If the difference is large enough, it will calculate the slope using OLS regression. If the difference is too small, it will increase dx and redo the loop with this new dx

4) This continues for all the datapoints

I think the more obvious reason for slow execution is the LOOPING nature of your code, when perhaps you could use the VECTORIZED (array-based operations) nature of Numpy.

For step 1, instead of taking pairs of points, you can perform directly `data[3:] - data[-3:] and get all the differences in a single array operation;

For step 2, you can use the result from array-based tests like numpy.argwhere(data > threshold) instead of testing every element inside some loop;

Step 3 sounds conceptually wrong to me. You say that if the difference is too small, it will increase dx. But if the difference is small, the resulting slope would be small because it IS actually small. Then, getting a small value is the right result, and artificially increasing dx to get a "better" result might not be what you want. Well, it might actually be what you want, but you should consider this. I would suggest that you calculate the slope for a fixed dx across the whole data, and then take the resulting array of slopes to select your regions of interest (for example, using data_slope[numpy.argwhere(data_slope > minimum_slope)].

Hope this helps!