Add PyFFTW as a (Soft) Dependency

[seanlaw]

With pyffttw sliding dot product in a stable place, we should consider adding it as a soft dependency to STUMPY.

Things to consider:

1. How to fallback to Python/SciPy if PyFFTW is unavailable
2. Additional logic for when to switch between different algorithms (including MASS_V4)

[NimaSarajpoor]

can we start simple? and then improve the solution? First, see the performances below:

./timing.py -timeout 1.0 -pmin 7 -pmax 24 pyfftw njit pocketfft_r2c_c2r scipy_oaconvolve > timing.csv

So, based on the provided figure, I think we can go with njit & fftw, and if fftw is not available, we can do: njit & pocketfft:

# imports (for hard dependencies)
 
try:
    import pfftw
    FFTW_IS_AVAILABLE = True
except ImportError:
    FFTW_IS_AVAILABLE = False


def choose_sliding_dot_function(Q, T):
    if len(Q) == len(T):
        spd_module = njit_dot  # special case
    elif len(Q) <= 2 ** 7   # config
        sdp_module = njit_sdp
    else: 
        if FFTW_IS_AVAILABLE:
            sdp_module = pyfftw_sdp
        else:
            sdp_module = pocket_fft
    
     return sdp_module

Does this sound like a reasonable approach? Later... we can change it to use Lookup to choose module/function, and the fallback can also depend on (len(Q), len(T)).

[seanlaw]

Why not:

def choose_sliding_dot_function(Q, T):
    if len(Q) == len(T):
        spd_module = njit_dot  # special case
    elif len(Q) <= 2 ** 7   # config
        sdp_module = njit_sdp
    elif FFTW_IS_AVAILABLE:
        sdp_module = pyfftw_sdp
    else:
        sdp_module = pocket_fft
    
     return sdp_module

the only question that I have is how do we allow users to completely override this if/else logic? For example, let's say they somehow discover that it's faster on their hardware to do elif len(Q) <= 2 ** 9 instead of len(Q) <= 2 ** 7. As I think about it more, the above implementation implies that we check Q and T every time in order to decide which module to return but this doesn't make sense. It should just do the calculation and return the requested value:

def liding_dot_product(Q, T):
    if len(Q) == len(T):
        return njit_dot.sliding_dot_product(Q, T)  # special case
    elif len(Q) <= 2 ** 7   # config
        return njit_sdp.sliding_dot_product(Q, T)
    elif FFTW_IS_AVAILABLE:
        return pyfftw_sdp.sliding_dot_product(Q, T)
    else:
        return pocket_fft.sliding_dot_product(Q, T)

Again, the question remains, how do we allow the flexibility to swap out or modify the if/else logic in that case as well? We should think about this...

[NimaSarajpoor]

Why not:

Was hard for my mind to switch between two different KINDs of condition. We started with checking len(Q). So, it was a bit strange for my mind to suddenly see elif FFTW_IS_AVAILABLE:. So, I thought maybe I should consider a separate if-else for the sake of readability (??)

[NimaSarajpoor]

Again, the question remains, how do we allow the flexibility to swap out or modify the if/else logic in that case as well?

The purpose of the sliding_dot_product function is to compute SDP for ANY (len(Q), len(T)) (using the fastest algo). The domain of this function is $\{ n_{Q}, n_{T} | 2 \le n_{Q} \le n_{T}; n_Q, n_T \in \mathbb{N} \}$.

So, the main question is: how do we want to divide this domain space into different regions? For instance, len(Q) <= 2 ** 7 basically divides the 2D space len(Q), len(T) into two regions with one line. What if we allow users to define regions, and their corresponding algo?

QT_regions = [
     [(LB1_Q, UB1_Q), (LB1_T, UB1_T), sdp_module_1],  
     [(LB2_Q, UB2_Q), (LB2_T,  UB2_T), sdp_module_2],
     [(LB3_Q, UB3_Q), (LB3_T,  UB3_T), sdp_module_3],
      ...
]

def sliding_dot_product(Q, T, QT_regions):
      n_Q = len(Q)
      n_T = len(T)
      sdp = None
      for region in QT_regions:
            (LB_Q, UB_Q), (LB_T, UB_T), sdp_module = region
             if LB_Q <= n_Q < UB_Q and LB_T <= n_T < UB_T:
                  return sdp_module(Q, T)
           
      # when the union of regions is not comprehensive  
      return sdp_module_default(Q, T)

[seanlaw]

1. Maybe just call it boundaries instead of QT_regions since you are using LB and UB
2. To be consistent, I think we should use m = len(Q) and n = len(T)
3. They should be SDP "functions", not SDP "modules"