Feature request: speeding up zeta ts test by not always interpolating

[jonahpearl]

Hi there — thanks so much for creating this statistical test and for implementing a python package for it! I'm excited to use it. That said, I'm finding it a bit slow to run, so I wanted to make a request for the future if possible.

I'm looking specifically at zetatstest. My data is a long interpolated timeseries (30 minutes at 500 Hz), with about 100 sets of 100 types of events. It takes 5-10 seconds to run one pass of zetatstest, so that means it would take a few minutes overall. The slowest parts of the code on my data, based on output from a profiler, are getTimeseriesOffsetOne, and within that where it calls getInterpolatedTimeSeries (it spends ~90% of its time there). It looks like zetapy is automatically choosing to re-interpolate the data every run, regardless of whether the data is already smoothly sampled (which it is, in my case) — it makes sense that this would be a bit slow.

To avoid this, I would suggest allowing the user to provide a optional sampling rate parameter, which would indicate that the data are already interpolated and don't need it again. I'd be happy to try implementing this myself, though I don't understand a lot of the other bits of the code (eg here, what exactly is vecTime) and I'm wary of making changes that accidentally mess up the downstream stats.

If you can confirm that the interpolation is purely superficial (i.e. it's just because some data might not be interpolated, rather than, say, because some pre-processing step is un-interpolating the data and it needs to be re-interpolated), and can help me see what vecTime is for, I'm more than happy to have a go at this change.

(There are some other small things that could also provide speed ups, especially on long timeseries data like mine. For example, the findfirst() function is essentially 2*O(n) right now, since it does a boolean comparison on each value and then calls np.where on the result. Since you've already sorted vecTimestamps, you could use a binary search to go faster.)

Thanks!

[JorritMontijn]

Hi Jonah, thank you for your comment! I agree that there are still several optimizations to be made. The following is an explanation plus my own stream of consciousness as to what the optimal solution could be, so apologies if it's somewhat incoherent

First let me explain why a solution isn't simple:

When the sample times are aligned with your event times, you can use direct sample indexing without interpolation, but this only works for the unjittered deviation vector. Once you start jittering the event times (which can be any real number), that alignment between event times and samples will be gone. The way I solve that now is by first defining a reference time vector (this is the vecTime you referred to), and then interpolating the jittered values from the time-series data that have been shifted by some arbitrary amount by the random-jitter procedure. So simply skipping the interpolation, even if you start with data where sample times and event times are aligned, is not possible.

Why work with interpolation in the first place:

In many cases, the event times will not be aligned to sample times, which will require interpolation regardless of what else you might want to do. More importantly, using a reference time vector will ensure that you always have the same number of samples per event-window. Since ZETA is based on the maximum of the deviation, it is dependent on the total number of samples: more samples will lead to higher maxima (on average, statistically speaking). Therefore, you want to keep the number of samples identical across jitter iterations, and between the real data and the jittered data. This has the added benefit of being able to work with heterogeneous sample times (since the interpolation will get rid of any inconsistencies), and even sampling at a super-resolution level by interpolating to, for example 1/100th of the original sample duration, allowing for a higher temporal resolution. So to conclude, the interpolation makes the whole procedure more robust and/or independent of the specific sample acquisition times.

A bad solution:

For the matlab version I was also annoyed by the slow running time, so after the initial release I added an extra route that bypasses the interpolation. This unfortunately turned out to be more involved than I had hoped, and is still imperfect. First, I discretize the jitters to sample-time steps if the number of samples per window is larger than the number of requested resamplings. This check is required because if it's not, then it won't be possible to create enough unique jitters and you would end up with duplicate resamplings. Secondly, I check whether a super-resolution procedure was requested, because if it was, then we'll have to interpolate no matter what. If it's not, then I pass the reference time vector to getInterpolatedTimeSeries, which in turn checks whether the first and last sample times are matches to the first and last reference time vector entries (if they are, I'll assume the rest are as well). Unfortunately, you cannot use exact matching here, since the sample times are stored as floating point numbers. So when you take sample(117385) minus event(116) where both are (e.g.) 10^4 versus sample(100) minus event(1) where both are (e.g.) 10^0, you end up with different precisions, so any isequal will say no. So here you also have to you use some critical value to determine when samples are similar-enough. Anyway, if they match, I use direct indexing; otherwise I won't. This speeds things up considerably and has the advantage of still using the same set of functions to perform the procedure, but with the amount of checks that are required, it's still not as fast as it could be.

A potentially better solution:

In hindsight I think it might be better to just go with a completely new parallel route with separate functions that you can use with an input switch to zetatstest, where you use integer indexing of samples rather than the convoluted floating-point-similar-enough-checking-nonsense I described above. The floating-point entries not matching between global and local times still remains a problem, however, but perhaps we'll just have to accept that as an inaccuracy in favour of computation speed. So then we'd have:

0. check if discrete processing is requested and possible(!) by checking if the number of samples per window is larger than the number of requested resamplings
1. if 0 checks out, then find sample entries corresponding to each event time => create integer lookup-table
2. run discretized deviation vector calculation on real data
3. create randomized integer array of jitters for each random iteration, and shift the event times by x_j samples, then run the discretized deviation vector calculation on each random iteration
4. calculate the p-value and return the results

Do you have thoughts on this?

In addition, the current interpolating procedure could indeed still be sped up. If I understand correctly, I could simply replace
intStartSample = findfirst(vecSpikeTimes >= dblEventT)
with
intStartSample = bisect.bisect_left(vecSpikeTimes , dblEventT)

and
intEndSample = findfirst(vecSpikeTimes > (dblEventT+dblWindowDur))
with
intEndSample = bisect.bisect_right(vecSpikeTimes,dblEventT+dblWindowDur)
?

[jonahpearl]

Hi Jorrit — thanks for the thorough and quick response! A few points from the perspective of timeseries data. Let's assume that the data are densely and evenly sampled. By that I mean, they are sampled fast enough that you are going to capture, with high resolution, the expected changes in the signal; and that they are already interpolated to (or were recorded with) a constant sampling rate.

When the sample times are aligned with your event times, you can use direct sample indexing without interpolation, but this only works for the unjittered deviation vector. Once you start jittering the event times (which can be any real number), that alignment between event times and samples will be gone.

I agree with you about the problem, but I think there is an easier way to solve it. For example, if a neuron is changing its firing rate every few msec, then sampling at 30 kHz (the standard for ephys these days) is going to be more than fast enough to capture those changes. So a very tiny amount of imprecision in your data alignment isn't going to change anything. So if your jittered event time doesn't fall on a real timestamp from the recording, that's ok — just take the nearest timestamp. Since the data are changing slowly wrt each timestep, you don't gain any "jittered-ness" from being overly precise.

In many cases, the event times will not be aligned to sample times, which will require interpolation regardless of what else you might want to do. More importantly, using a reference time vector will ensure that you always have the same number of samples per event-window.

Same comment as above — most of the time, I'd bet it's ok to have the event times be slightly mis-aligned to the sample times, as long as the sample timestamps are dense enough that the imprecision is irrelevant. For example, if you had smoothed neural firing rates sampled at 1 kHz (every msec), and your event time was some fraction of a msec, the firing rate would have to be changing quite fast in order to get an issue with nudging the event time over by 1 minus that fraction.

Since ZETA is based on the maximum of the deviation, it is dependent on the total number of samples: more samples will lead to higher maxima (on average, statistically speaking). Therefore, you want to keep the number of samples identical across jitter iterations, and between the real data and the jittered data.

Again, if the data are densely and evenly sampled, this shouldn't be an issue. (I haven't thought about how all of this translates to spike trains.)

Ok this comment is long enough, will make a new one for discussing implementation details.

[jonahpearl]

0. check if discrete processing is requested and possible(!) by checking if the number of samples per window is larger than the number of requested resamplings

I don't see why it's important to have the resamplings be unique — indeed, statistical bootstrapping (not really what this is, but it's related) relies on resampling data with replacement. I assume that increasing the number of resamplings here has the same effect as bootstrapping for more iterations, i.e. it increases the statistical power of your test. In that case I don't see a reason not to do it with replacement (in the discrete case).

1. if 0 checks out, then find sample entries corresponding to each event time => create integer lookup-table
2. run discretized deviation vector calculation on real data
3. create randomized integer array of jitters for each random iteration, and shift the event times by x_j samples, then run the discretized deviation vector calculation on each random iteration
4. calculate the p-value and return the results

Yes, and I have implemented a fast version that basically does this. It relies on this library (docs) of functions that I've put together for solving this exact sort of problem efficiently. You can see the core of the refactored zetapy implementation here. The main idea is that if you know the sampling rate of the data, you can generate a matrix of aligned traces without doing any calculation / interpolation. That is implemented in the _collate_events function. (This also lets you skip pseudovector creation.) Then for jittering, I sort of hacked my way back to your original implementation where you stitch together the relevant parts of the data, and take jittered windows from that stitched vector. (I also took the liberty of refactoring a few other functions just to make it easier for me to read.)

I made a notebook to test my version and compare it to the original ("og") version, using some fake data. The results are nearly identical, though my p-values are slightly higher (less signf), probably because I had to filter any jittered events that were at the boundaries of the stitched data, and so I just have fewer shuffles overall. I think if we thought carefully about how to deal with boundaries, the implementations would match.

In terms of timing, if you pass a sampling rate, it is ~65x faster (1.33 sec vs 19.7 msec). If you don't, it is still ~20x faster (1.33 sec vs 61.4 msec).

If you think this is potentially useful, let me know and I can open a PR with this branch. I'd be happy to help you re-factor the rest of the python codebase, it seems like it wouldn't be that different, though again, I haven't thought about this wrt spike trains!

[JorritMontijn]

Hi Jorrit — thanks for the thorough and quick response! A few points from the perspective of timeseries data. Let's assume that the data are densely and evenly sampled. By that I mean, they are sampled fast enough that you are going to capture, with high resolution, the expected changes in the signal; and that they are already interpolated to (or were recorded with) a constant sampling rate.

[snip]

Again, if the data are densely and evenly sampled, this shouldn't be an issue. (I haven't thought about how all of this translates to spike trains.)

Ok this comment is long enough, will make a new one for discussing implementation details.

Hi Jonah, I see your point in the case of ephys, but I don't think this will be true for most calcium imaging experiments. That's why I think the use-discrete-samples-switch would be useful. Then you can use the interpolating version for low acquisition rate data (e.g., two-photon GCaMP), and the discrete non-interpolating version for high rate data (e.g., MUA, LFP, EEG).

0. check if discrete processing is requested and possible(!) by checking if the number of samples per window is larger than the number of requested resamplings

I don't see why it's important to have the resamplings be unique — indeed, statistical bootstrapping (not really what this is, but it's related) relies on resampling data with replacement. I assume that increasing the number of resamplings here has the same effect as bootstrapping for more iterations, i.e. it increases the statistical power of your test. In that case I don't see a reason not to do it with replacement (in the discrete case).

Even in the case of resampling with replacement, you still need a sufficient number of samples. Bootstrapping with 2 data points is still going to give you a 3-value 1/4-1/2-1/4 distribution, no matter how many iterations you use. However, I agree that #samples >= #resamplings is not a hard boundary if you allow duplicates. Maybe it would indeed be better to leave this to the judgment of the user instead of imposing an arbitrary threshold.

1. if 0 checks out, then find sample entries corresponding to each event time => create integer lookup-table
2. run discretized deviation vector calculation on real data
3. create randomized integer array of jitters for each random iteration, and shift the event times by x_j samples, then run the discretized deviation vector calculation on each random iteration
4. calculate the p-value and return the results

Yes, and I have implemented a fast version that basically does this. It relies on this library (docs) of functions that I've put together for solving this exact sort of problem efficiently. You can see the core of the refactored zetapy implementation here. The main idea is that if you know the sampling rate of the data, you can generate a matrix of aligned traces without doing any calculation / interpolation. That is implemented in the _collate_events function. (This also lets you skip pseudovector creation.) Then for jittering, I sort of hacked my way back to your original implementation where you stitch together the relevant parts of the data, and take jittered windows from that stitched vector. (I also took the liberty of refactoring a few other functions just to make it easier for me to read.)

I made a notebook to test my version and compare it to the original ("og") version, using some fake data. The results are nearly identical, though my p-values are slightly higher (less signf), probably because I had to filter any jittered events that were at the boundaries of the stitched data, and so I just have fewer shuffles overall. I think if we thought carefully about how to deal with boundaries, the implementations would match.

In terms of timing, if you pass a sampling rate, it is ~65x faster (1.33 sec vs 19.7 msec). If you don't, it is still ~20x faster (1.33 sec vs 61.4 msec).

If you think this is potentially useful, let me know and I can open a PR with this branch. I'd be happy to help you re-factor the rest of the python codebase, it seems like it wouldn't be that different, though again, I haven't thought about this wrt spike trains!

That does sound pretty useful. I've got my hands full at the moment, but I'll take a look after the christmas holidays (early january). I'm a bit worried about the higher p-values though, so if you can fix that, I'll integrate it into the zetapy repository. Do you think you could make a pull request to the dev branch of zetapy by early january, so I can take a look at it?

By the way, thanks a lot for your input, I really appreciate it!

Best wishes,
Jorrit

[jonahpearl]

Sounds good, I can see the concern for Ca2+ imaging! I will have a go at working out the boundary problems to make the pvals closer. Do you have other (real) test data I could run it on as well? I can make a set of automated tests that verify everything is working as expected.