Add MASS V4 Challenger

[seanlaw]

While looking up the original MASS code, I noticed a new MASS V4. It claims to be:

1. Numerically stable (especially for longer time series) and without producing imaginary numbers
2. Exploits Direct Cosine Transform (DCT), which is suppose to make MASS faster

% This code is created by Sheng Zhong and Abdullah Mueen
% The overall time complexity of the code is O(n log n). The code is free to use for research purposes.
% The code does not produce imaginary numbers due to numerical errors 
% k should greater than or equals to floor((3m+1)/2)

function [dist] = MASS_V4(T, Q, k)
    n = length(T);
    m = length(Q);
    Q = zNorm(Q);
    dist = [];
    batch=get_batch_size(k, m);
    
    for j = 1 : batch-m+1 : n-m+1
        right = j + batch - 1;
        if right >= n
            right = n;
        end 
        dot_p = DCT_dot_product(T(j:right), Q);
        sigmaT = movstd(T(j:right),[m-1 0],1);
        d = sqrt(2*(m-(dot_p./sigmaT(m:end)))) ;        
        dist = [dist; d];
    end
    
end

function dct_batch_size = get_batch_size(target_batch_size, m)
    dct_batch_size = floor((2*target_batch_size-2)/3) -1;
    if dct_batch_size < m
        dct_batch_size = m;
    end
    pad_length = dct_batch_size + floor((dct_batch_size-m+1)/2) + floor((m+1)/2);
    while pad_length < target_batch_size
        dct_batch_size = dct_batch_size + 1;
        pad_length = dct_batch_size + floor((dct_batch_size-m+1)/2) + floor((m+1)/2);
    end
    if (pad_length > target_batch_size)
        dct_batch_size = dct_batch_size - 1;
    end
end


function dot_p = DCT_dot_product(x, y)
    n = length(x);
    m = length(y);
    [x_pad, y_pad, si] = DCT_padding(x,y);
    
    N = length(x_pad);
    
    Xc = dct(x_pad, 'Type', 2);
    Yc = dct(y_pad, 'Type', 2);
    
    dct_product = Xc .* Yc;
    dct_product(N + 1) = 0;
    dct_product(1) = dct_product(1) * sqrt(2);
    dot_p = dct(dct_product, 'Type', 1);
    dot_p(1) = dot_p(1) * 2;
    dot_p = sqrt(2*N) * dot_p(si: si+n-m);
end

function [x_pad, y_pad, start_index] = DCT_padding(x, y)
    n = length(x);
    m = length(y);
    p2=floor((n-m+1)/2);
    p1=p2+floor((m+1)/2);
    p4=n-m+p1-p2;
    x_pad = zeros(n+p1, 1);
    x_pad(1+p1:end) = x;
    y_pad = zeros(m+p2+p4, 1);
    y_pad(1+p2:1+p2+m-1) = y;
    start_index = p1-p2+1;
end

[NimaSarajpoor]

@seanlaw
Noticed that when I was looking at their webpage. Should have created an issue for it! Thanks for creating this! 👍

[NimaSarajpoor]

The MASS webpage also point to the authors' MASS paper published in 2024. See section 4.4 to read more about "MASS V4".

[seanlaw]

V4 executes 46% to 48% more FLOPs than V3 due to the extra padding with zeros and most importantly, the DCT always executes more FLOPs than FFT (Johnson and Frigo 2007; Shao and Johnson 2008)

More FLOPs generally means "faster" but, in this case, they mean "Total number of floating point operations to compute the sliding dot product" (so "more FLOPs" is actually "slower" - i.e., a smaller number is faster in the table below)!

Based on this, V3 is the fastest and then followed by V4. V1 and V2 are considered "FFT"/"DFT" based methods and slightly slower. One thing to note is that there is no guarantee that STUMPY has implemented V2 because we haven't directly performed zero-padding on Q so that it matches the length of T. We've simply left everything up to the scipy.signal.convolve function. I wonder if we can achieve the same/similar 2x speed up observed in MASS_V2 by:

1. Manually zero-padding Q to be the same length as T
2. Apply pocketfft rfft/irfft ourselves
3. I wonder if this will match the performance of FFTW?

Note that while MASS_V3 is the fastest, MASS_V4 offers much better stability (no imaginary numbers) and it is still faster than FFT-based MASS methods. So, MASS_V4 is a reasonable compromise (assuming that it works) but we should also test the MASS_V2 with scipy to see if it helps at all

[seanlaw]

I wonder if we can achieve the same/similar 2x speed up observed in MASS_V2 by

Never mind. After some tinkering, I realized that we've already done exactly this in pocketfft_r2c_c2r_sdp

Nonetheless, I think MASS_V4 is still worth pursuing. Even if a SciPy-based MASS_V4 is still slower than FFTW/pyFFTW, it might still be faster than scipy.signal.convolve or pocketfft because it would be using a faster dct method. What do you think @NimaSarajpoor?

[NimaSarajpoor]

Never mind. After some tinkering, I realized that we've already done exactly this in pocketfft_r2c_c2r_sdp

Right! And, our version should be better as we use next_fast_len(n) instead of n

I think MASS_V4 is still worth pursuing. Even if a SciPy-based MASS_V4 is still slower than FFTW/pyFFTW, it might still be faster than scipy.signal.convolve or pocketfft

Right! Will check this out in the upcoming days...