Create Lookup

[seanlaw]

Continuing this conversation, it would be good to come up with a rough strategy for which algorithm to use depending on the len(T).

To get started, I explored the amount of storage that we might need to keep track of everything in memory:

#!/usr/bin/env python

from scipy.fft import next_fast_len
import numpy as np
from array import array
import sys

lengths = set()
for i in range(2**28):
    lengths.add(next_fast_len(i, real=True))

fast_lengths = np.empty(len(lengths), dtype=np.uint32)
fast_methods = np.empty(len(lengths), dtype=np.uint8)
for i, l in enumerate(sorted(list(lengths))):
    fast_lengths[i] = l
    fast_methods[i] = 1

print(f"{fast_lengths.nbytes / 1024} kb, {sys.getsizeof(fast_lengths.tolist()) / 1024} kb")
print(f"{fast_methods.nbytes / 1024} kb, {sys.getsizeof(fast_methods.tolist()) / 1024} kb")

This returns:

4.9921875 kb, 10.0390625 kb
1.248046875 kb, 10.0390625 kb

So, fast_lengths stores all of the possible len(T) generated by next_fast_len from 0 to 2**28 and fast_methods stores the algorithm index number (for now, I've only set the all of the algorithms to 1). The printed numbers at the bottom corresponds to the amount of space needed for storing things in a numpy array versus a python list (both in-memory).

In conclusion, storing all of this information is very cheap (as long as it is 1D and not 2D). I should probably simplify and store everything in a simple Dict and see how much space that will take up.

[seanlaw]

Okay, storing everything as a Dict only costs 36.0859375 kb (6x more than numpy total) so we should probably do that and avoid unnecessary data structures?

[NimaSarajpoor]

@seanlaw

In #16 (comment):

For len(T) >= len(Q), we have shape >= len(T) >= len(Q), and we pad them to get to the length shape. So, can we say that we are okay to fix len(T), and check performance for one Q?

I see what you mean! If I understand correctly, you're saying that len(T) is always greater than or equal to len(Q), which then implies that next_fast_len(len(T)) always dictates the shape for 2D rfft step. Thus, it never matters what the len(Q) is because it always gets padded to next_fast_len(len(T))!

I noticed that the logic I proposed before probably does NOT work for scipy_oaconvolve_sdp when it internally triggers the overlap-add approach. When len(T) is about 8 or 16 times larger than len(Q), the ovarlap-add approach is triggered. In that case, FFT is applied to a 2D array that contains some chunks of the time series T, and the length of the chunk depends on len(Q).

[seanlaw]

If needed, we can create a 2D lookup table if it makes sense to? Perhaps, we should consider an efficient data structure to hold the information that we need and that allows us to quickly query and find out what algo should be used?

So, there are a few separate things:

1. Given len(Q) and len(T), what algo should we use?
2. How should we store this info?
3. How do we give the user the freedom to modify this info and/or supply their own algo (i.e., pyFFTW)?

Let's propose a reasonably simple solution first that is easy to maintain rather than seeking a perfect solution upfront and then we can refine it from there.

[seanlaw]

Here is a suggestion. I think we can still store a 2D array quite efficiently with:

from scipy.fft import next_fast_len
import numpy as np


lengths = set()
for i in range(2**28):
    lengths.add(next_fast_len(i, real=True))
lengths = np.array(list(lengths), dtype=np.uint32)

# lengths contains 1278 elements
sdp_lookup = np.empty((len(lengths), len(lengths)), dtype=np.uint8)  # This is only 1633284 bytes or ~1.63 mb

sdp_algos = ["njit", "r2r/c2c", "oaconvolve"]  # These should really be direct references to those functions

Then, for each len(Q) and len(T) pair, you can fill sdp_lookup with the index of the algorithm to use (e.g., 0 is "njit", 1 is r2r/c2c, and 2 is "oaconvolve").

Later, when you have a new len(Q) and len(T), one can simply use Q_idx = np.searchsorted(lengths, len(Q)) and T_idx = np.searchsorted(lengths, len(T)) to look up which algo to use in algo_idx = sdp_lookup[Q_idx, T_idx]