Fixed #26 Use less-aggressive planning flag as default (#27)

* less intensive fft planning as default

* use multiply param out to save result

* add matlab sdp code

* minor improvement in comment

* minor changes

Author: NimaSarajpoor

matlab_sdp.m

@@ -0,0 +1,27 @@
+% This function computes the sliding dot product between
+% a query Q and a time series T using the FFT method.
+
+% See TABLE I (old) here: https://www.cs.ucr.edu/~eamonn/PID4481997_extend_Matrix%20Profile_I.pdf
+% See MASS_V2 here: https://www.cs.unm.edu/~mueen/FastestSimilaritySearch.html
+% MASS_V2 is used in [DAMP_2_0.m](https://drive.google.com/file/d/1EPDhFXQ2goTJ5m_x1S84-KNpQFsRNVmf/view?usp=sharing)
+
+function [z] = fft_sdp(Q, T)
+    % Input:
+    %   Q - query vector
+    %   T - time series vector
+    % Output:
+    %   z - sliding dot product of Q and T
+
+    m = length(Q);
+    n = length(T);
+
+    Q = Q(end:-1:1);  % Reverse the query
+    Q(m+1:n) = 0;  % Append zeros
+
+    X = fft(T);
+    Y = fft(Q);
+    Z = X.*Y;
+    z = ifft(Z);
+    z = z(m:n);
+
+end

sdp/pyfftw_sdp.py

@@ -21,7 +21,7 @@ def __init__(self, max_n=2**20):
         self.rfft_objects = {}
         self.irfft_objects = {}
 
-    def __call__(self, Q, T, n_threads=1, planning_flag="FFTW_MEASURE"):
+    def __call__(self, Q, T, n_threads=1, planning_flag="FFTW_ESTIMATE"):
         """
         Compute the sliding dot product between `Q` and `T` using FFTW via pyfftw.
 
@@ -98,7 +98,7 @@ def __call__(self, Q, T, n_threads=1, planning_flag="FFTW_MEASURE"):
         # RFFT(Q)
         # Scale by 1/next_fast_n to account for
         # FFTW's unnormalized inverse FFT via execute()
-        real_arr[:m] = Q[::-1] / next_fast_n
+        np.multiply(Q[::-1], 1.0 / next_fast_n, out=real_arr[:m])
         real_arr[m:] = 0.0
         rfft_obj.execute()  # output is in complex_arr
 
@@ -114,10 +114,10 @@ def __call__(self, Q, T, n_threads=1, planning_flag="FFTW_MEASURE"):
 _sliding_dot_product = SLIDING_DOT_PRODUCT()
 
 
-def setup(Q, T, n_threads=1, planning_flag="FFTW_MEASURE"):
+def setup(Q, T, n_threads=1, planning_flag="FFTW_ESTIMATE"):
     _sliding_dot_product(Q, T, n_threads=n_threads, planning_flag=planning_flag)
     return
 
 
-def sliding_dot_product(Q, T, n_threads=1, planning_flag="FFTW_MEASURE"):
+def sliding_dot_product(Q, T, n_threads=1, planning_flag="FFTW_ESTIMATE"):
     return _sliding_dot_product(Q, T, n_threads=n_threads, planning_flag=planning_flag)