Fixed #34 Customizing scipy's oaconvolve

[NimaSarajpoor]

This PR is to address #34.

[gitnotebooks]

Review these changes at https://app.gitnotebooks.com/stumpy-dev/sliding_dot_product/pull/35

[NimaSarajpoor]

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

# in timing.py, I change timeout to 5.0 when `len(T) >= 2^20`

The challenger is the customized version of scipy's oaconvolve.

Observations:

• The challenger outperforms scipy's oaconvolve
• For len(Q) <= 2 ^16 (and len(Q)>= 2^7), challenger outperforms the others for the most part.
• For len(Q) > 2^16, pocketfft outperforms the others for the most part.

For me, the important one is the first bullet point. Of the four optimization opportunities mentioned in this comment, I've addressed 1, 2, and 3 in this PR. The last item, which is about adjusting the number of multiplication for real-valued arrays, can be explored next.

[seanlaw]

As a gentle reminder, even if we can do things faster, we will never (??) remove the public scipy convolution functions from STUMPY because they should be our last resort fallback (in case the alternatives, that may use private functions, raise an error). Does that make sense?

[NimaSarajpoor]

Good reminder. It makes sense!!

[NimaSarajpoor]

This line needs to be modified if, at a later time, we decide to move the proposal to a new file (module).

[NimaSarajpoor]

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

# in timing.py, I change timeout to 5.0 when `len(T) >= 2^20`

@seanlaw
"Challenger" seems to be the winner for most cases, and I think it is worth it to include it. What do you think? Also, can you please review the script? I've made major changes. The private objects are now only r2c and c2r. IMO, the script looks cleaner now.

[seanlaw]

@NimaSarajpoor I've left some comments but would still like another pass after you've cleaned things up further

I do agree that, for the most part, things look clean. I think it still lacks clarity as to what is happening or why the logic is coded in this way

[NimaSarajpoor]

why?

[NimaSarajpoor]

As I mentioned before, this is purely based on scipy's logic:

https://github.com/scipy/scipy/blob/8c75ae75176236f233824e9a0483c26a69e6dfec/scipy/signal/_signaltools.py#L748-L750

However, the reason is not clear to me.

[NimaSarajpoor]

The overlap-add method adds the last m-1 elements of a block to the first m-1 elements of the next block. Therefore, if conv_block_size is at least 2 * (m-1), then the head and tail of each block will not have common elements and the current implementation works. If conv_block_size is less than 2 * (m-1), the current implementation fails. So: conv_block_size >= 2(m-1)

Also, note that conv_block_size should be <n. Otherwise, there is no point in splitting T into blocks.

Therefore: 2(m-1) <= conv_block_size < n , which gives: m < n/2 + 1

---

Note:
The current implementation shows conv_block_size = max(conv_block_size, m). However, as mentioned above, the sdp function fails if conv_block_size is < 2 * (m-1). So, we need to:

(1) use conv_block_size = max(conv_block_size, 2 * (m-1)) instead of conv_block_size = max(conv_block_size, m)

(2) Or, we need to revise the implementation so that it can handle cases where conv_block_size is < 2(m-1).

[NimaSarajpoor]

Also: I submitted a question in stack overflow regarding the condition if m >= n / 2:

[NimaSarajpoor]

The overlap-add method adds the last m-1 elements of a block to the first m-1 elements of the next block. Therefore, if conv_block_size is at least 2 * (m-1), then the head and tail of each block will not have common elements and the current implementation works.

If the condition conv_block_size < 2 * (m-1) is not met, the overlap-add approach should still work. See example below.

Example:
Let's compute the sliding dot product between T=[1, 2, 3, 4] and Q=[10, 20, 30]. The sliding dot product is: [140, 200].

The sliding dot product between T and Q is equivalent to the valid convolution between T and Qr, which is the reverse of Q. Let's use overlap-add method and set conv_block_size to m=len(Qr) == 3.

# Each block has zero-padding of length `m-1`

block1: 1, 0, 0 
block2: 2, 0, 0
block3: 3, 0, 0
block4: 4, 0, 0

# reverse of Q
Qr = [30, 20, 10]

In overlap-add method, we compute the circular convolution between each block and Qr. We can use the flip-and-slide method. This gives:

out_block1: 30, 20, 10
out_block2: 60, 40, 20
out_block3: 90, 60, 30
out_block4: 120, 80, 40

We now add m-1 == 2 elements of each block to the first m-1 == 2 elements of next block.

out: out_block1 & out_blokc2 --> 30, 80, 50, 20
out: out & out_block3 --> 30, 80, 140, 80, 30
out: out & out_block4 --> 30, 80, 140, 200, 110, 40

Get the slice (m-1, n), which is (2, 4):
[140, 200]

So, the overlap-add method works if conv_block_size is set to m. The current implementation uses vectorized operation to perform the overlap-add between blocks. However, currently, the vectorized operation used in our implementation fails when conv_block_size is set to a value that is less than 2 * (m - 1).