Add pyfftw sdp

[NimaSarajpoor]

This is to address PR 3 described in #1118 (comment). Have copied the corresponding notes below:

1. Add _PYFFTW_SLIDING_DOT_PRODUCT to sdp.py
2. Add unit tests that demonstrate that _PYFFTW_SLIDING_DOT_PRODUCT matches the results of naive_rolling_window_dot_product
3. Handle the test.sh check_fftw_pyfftw stuff here too

[gitnotebooks]

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

[seanlaw]

@NimaSarajpoor I've left some comments for you to address. I think we can afford to be clearer since all of this pyfftw stuff will be hard to maintain. We should probably be as verbose (and add more comments cross referencing their docs as possible). I'll do another few passes after you've responded and made modifications. I think this pyfftw stuff needs to be crystal clear

[NimaSarajpoor]

@seanlaw
This PR should be ready to be merged if there are no concerns.

[seanlaw]

@NimaSarajpoor This might be a distraction but a (very, very small) part of me was wondering if we should be using a closure instead of a class:

https://realpython.com/python-closure/

I don't know if it makes things easier/harder to read. Given all of the array allocations, maybe a class is the right thing to use. I just wanted to bring it to your attention since I'm not sure we're getting the full benefits of a class since:

1. None of the class attributes are actually being accessed externally
2. We don't really have an external method (except __call__) but we do want to make it a function

Just something to consider. No pressure though.

[NimaSarajpoor]

We don't really have an external method (except call) but we do want to make it a function

This is something that bothered me too but thought the reason for being bothered is just my lack of experience and, also, I didn't know/think about an alternative option.

but a (very, very small) part of me was wondering if we should be using a closure instead of a class:
https://realpython.com/python-closure/

Thanks for sharing this! I will go through this and try to modify the script. Let's see how it will look like.

[NimaSarajpoor]

@seanlaw
The code is modified to use closure. I pushed the changes so that you can see it too. It looks readable to me. Also, this way the user cannot change the attributes that are used internally. I think that should be a good thing. Right?

[seanlaw]

The code is modified to use closure. I pushed the changes so that you can see it too. It looks readable to me.

Yeah, it feels pretty natural! There is one thing that still bothers me and that is "how do you change/reset max_n after the function has been created?"

The natural thing to do is to overwrite the existing sdp._pyfftw_sliding_dot_product with a new call to _make_pyfftw_sliding_dot_proudct along with your desired max_n but feels a bit cumbersome (or hard to remember). Did you have any thoughts about this?

It feels like we're very close.

Also, this way the user cannot change the attributes that are used internally. I think that should be a good thing. Right?

Indeed. This is partially why I didn't like using the class approach!

[NimaSarajpoor]

@seanlaw

There is one thing that still bothers me and that is "how do you change/reset max_n after the function has been created?"

Did you have any thoughts about this?

I did not think about this scenario. I am now thinking in what scenario(s) a user might want to change max_n. Maybe they need to process much larger arrays at some point and they want to increase max_n to a certain value once to cover those lengths (?) or maybe they already processed some large arrays and now there is no need to have large preallocated arrays.

We can define a new function inside _make_pyfftw_sliding_dot_product and then return that object too. Then users can leverage that function to change values. Not sure if that is a good approach though as now users need to track two different objects returned by _make_pyfftw_sliding_dot_product. The option class might be better if we are going to end up with having more than one object returned by _make_pyfftw_sliding_dot_product.

[NimaSarajpoor]

@seanlaw

There is one thing that still bothers me and that is "how do you change/reset max_n after the function has been created?"

Did you have any thoughts about this?

I did not think about this scenario. I am now thinking in what scenario(s) a user might want to change max_n. Maybe they need to process much larger arrays at some point and they want to increase max_n to a certain value once to cover those lengths (?) or maybe they already processed some large arrays and now there is no need to have large preallocated arrays.

We can define a new function inside _make_pyfftw_sliding_dot_product and then return that object too. Then users can leverage that function to change values. Not sure if that is a good approach though as now users need to track two different objects returned by _make_pyfftw_sliding_dot_product. The option class might be better if we are going to end up with having more than one object returned by _make_pyfftw_sliding_dot_product.

Okay... on second thought, I think it should be okay to allow the returned callable object has a "callable feature" for resetting max_n. I think most users do not need to reset the pre-allocated arrays, and this design should be okay. See the changes in the last commit.

[seanlaw]

I think it should be okay to allow the returned callable object has a "callable feature" for resetting max_n

Can you please show me a complete example of how this works (i.e., from instantiating a new object to changing the max_n)?

[NimaSarajpoor]

@seanlaw

I think it should be okay to allow the returned callable object has a "callable feature" for resetting max_n

Can you please show me a complete example of how this works (i.e., from instantiating a new object to changing the max_n)?

# import
import numpy as np

from stumpy import sdp

# naive function to assert output
def naive_sdp(Q, T):
    out = np.zeros(len(T) - len(Q) + 1, dtype=np.float64)
    for i in range(len(out)):
        out[i] = np.dot(Q, T[i : i + len(Q)])
    return out


# pyfftw (default max_n: 2**20)
pyfftw_sdp = sdp._pyfftw_sliding_dot_product 

# case 1
Q = np.random.rand(2**2)
T = np.random.rand(2**15)
ref = naive_sdp(Q, T)
comp = pyfftw_sdp(Q, T)
np.testing.assert_allclose(comp, ref)

# set pre-allocated arrays' size to 2^10
pyfftw_sdp.reset_arr(2**10)

# case 2 
# len(T) < length of preallocated real-valued array
Q = np.random.rand(2**2)
T = np.random.rand(2**5)
ref = naive_sdp(Q, T)
comp = sdp._pyfftw_sliding_dot_product(Q, T)
np.testing.assert_allclose(comp, ref)

# case 3
# len(T) > length of preallocated real-valued array
Q = np.random.rand(2**2)
T = np.random.rand(2**11)
ref = naive_sdp(Q, T)
comp = sdp._pyfftw_sliding_dot_product(Q, T)
np.testing.assert_allclose(comp, ref)

Maybe I should've used the name set_arr instead of reset_arr. Because, as shown in this example, I "set" it to smaller size rather than just resetting the arrays to their original sizes. (Btw, note that user cannot check if the array is resized. So, I am allowing user to do something that they cannot verify. Is that okay??)

[seanlaw]

Btw, note that user cannot check if the array is resized. So, I am allowing user to do something that they cannot verify. Is that okay??

I don't think I understand what you mean by "user cannot check if the array is resized". Can you articulate what could possibly go wrong? What would somebody need/want to verify?

Do the arrays automatically get resized if, later, somebody uses an array that is larger than max_n?

[NimaSarajpoor]

Do the arrays automatically get resized if, later, somebody uses an array that is larger than max_n?

Yes. Based on the current design, if n=len(T) > max_n, the preallocated array will be resized automatically, and its length will be set to n.

I don't think I understand what you mean by "user cannot check if the array is resized". Can you articulate what could possibly go wrong? What would somebody need/want to verify?

Please ignore. That was stupid!

---

Based on the current design, if n=len(T) > max_n, the preallocated array will be resized automatically

What do you think about this? Now that we are planning to add set_max_n, I was wondering if we should avoid allowing the length of preallocated arrays to be increased automatically. So, if the length of input is greater than max_n, we can create temporary preallocated arrays. This way we do not touch the original preallocated array with length max_n.

[seanlaw]

This way we do not touch the original preallocated array with length max_n.

I guess I would need to consider what the side effect would be if the preallocated were to be automatically resized (i.e., what would be the problem with this?). Currently, I don't see any issues with automatically resizing the preallocated array

[NimaSarajpoor]

@seanlaw
I just left one comment. I think it should be ready to merge after discussing/addressing that comment.

[seanlaw]

@NimaSarajpoor I've left some comments for you to consider. Once completed, please let me know as I'd like to go over the code in more detail before you move on (as it's been a while since I've gone over things).

[seanlaw]

The more I think about this, the more I feel like max_n is the wrong parameter name because "max" implies that the array length cannot exceed this number but we might actually violate that. Perhaps, what we're really talking about is the "initial size/length" and so maybe something like init_size or init_len.

[seanlaw]

Additionally, I notice that we are overloading n in n_threads and next_fast_n so the n is starting to either losing meaning or we are being inconsistent when we add max_n.

I can handle two different meanings for n within the same function but probably not more. I almost feel like I might be okay with nfl = pyfftw.next_fast_len(n) as long as it is defined clearly (i.e., this assignment isn't hidden/buried and I can quickly lookup what nfl is)

[NimaSarajpoor]

Currently, n appears in the following cases:

1. The len of T --> n
2. next_fast_len(n) --> next_fast_n
3. The (initial) len of (real-valued) preallocated array --> max_n
4. number of threads --> n_threads

---

1. Good
2. The value returned by next_fast_len(n) is strongly tied to n. So, I think it should be okay to keep next_fast_n. This is used in several lines. So using a name that reflects its description should help with the readability of the code.
3. As you pointed out, it is better to avoid using max_n. I think init_len is a good choice.
4. Can we use num_threads instead? I can see n_threads is used in other modules though as well.

[seanlaw]

2. Okay, keep it as next_fast_n
3. Cool. Let's do init_len then
4. I think n_threads is fine. We're already use to seeing it elsewhere and I never get it confused with other code.

I think simply changing max_n to init_len should be sufficient