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.shcheck_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.

[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).

[NimaSarajpoor]

@seanlaw
Ok... I think this should be ready now for another round of review.

[seanlaw]

@NimaSarajpoor I think everything looks good. Just a few comments for you to consider

[NimaSarajpoor]

For our future reference:

(R)FFT and I(R)FFT are implemented in a way that:

• IFFT(FFT(x)) = x
• IRFFT(RFFT(x)) = x

By convention, (R)FFT is scaled by 1 (i.e. no scaling), and I(R)FFT is scaled by 1 / len(x)

Example 1

import numpy as np
n_T = 4
T = np.random.rand(n_T)
R = np.fft.rfft(T)
IR = np.fft.irfft(R)
np.testing.assert_allclose(T, IR, rtol=1e-5, atol=1e-8)

Example 2

Let's use pyfftw instead.

import numpy as np
import pyfftw
n_T = 4
real_view = pyfftw.empty_aligned(n_T, dtype="float64")
complex_view = pyfftw.empty_aligned(n_T//2 + 1, dtype="complex128")
planning_flag = "FFTW_ESTIMATE"
n_threads = 1
rfft_obj = pyfftw.FFTW(
input_array=real_view,
output_array=complex_view,
direction="FFTW_FORWARD",
flags=(planning_flag,),
threads=n_threads,
)
irfft_obj = pyfftw.FFTW(
input_array=complex_view,
output_array=real_view,
direction="FFTW_BACKWARD",
flags=(planning_flag,),
threads=n_threads,
)
T = np.random.rand(n_T)
real_view[:] = T
rfft_obj()
# rfft_obj.output_array contains the RFFT of T
complex_view[:] = rfft_obj.output_array.copy()
irfft_obj()
# irfft_obj.output_array contains the IRFFT of the RFFT of T
np.testing.assert_allclose(T, irfft_obj.output_array, rtol=1e-5, atol=1e-8)

Example 3

Let's use pyfftw instead. However, this time, we will use the thin wrapper "execute" that does Fourier transform only without any scaling.

import numpy as np
import pyfftw
n_T = 4
real_view = pyfftw.empty_aligned(n_T, dtype="float64")
complex_view = pyfftw.empty_aligned(n_T//2 + 1, dtype="complex128")
planning_flag = "FFTW_ESTIMATE"
n_threads = 1
rfft_obj = pyfftw.FFTW(
input_array=real_view,
output_array=complex_view,
direction="FFTW_FORWARD",
flags=(planning_flag,),
threads=n_threads,
)
irfft_obj = pyfftw.FFTW(
input_array=complex_view,
output_array=real_view,
direction="FFTW_BACKWARD",
flags=(planning_flag,),
threads=n_threads,
)
T = np.random.rand(n_T)
real_view[:] = T
rfft_obj.execute()
# rfft_obj.output_array contains the RFFT of T
complex_view[:] = rfft_obj.output_array.copy()
irfft_obj.execute()
# irfft_obj.output_array contains the IRFFT of the RFFT of T
scaling_factor = T / irfft_obj.output_array
print("Scaling factor:", scaling_factor)

which gives: Scaling factor: [0.25 0.25 0.25 0.25]. As observed, scaling factor is 1/len(T).

[seanlaw]

For our future reference:

Should this information be recorded in https://github.com/stumpy-dev/sliding_dot_product/blob/main/docs/Tutorial_FFT_Based_SDP.ipynb instead??

[NimaSarajpoor]

For our future reference:

Should this information be recorded in https://github.com/stumpy-dev/sliding_dot_product/blob/main/docs/Tutorial_FFT_Based_SDP.ipynb instead??

YES, it should! Will add the info to the notebook (will update this comment)

---

Update: Created PR stumpy-dev/sliding_dot_product#39 to address this.

[NimaSarajpoor]

@seanlaw
I made minor changes in 3cd46d7. How do you feel about this PR now? Please let me know if anything else needs to be addressed in this PR.

[seanlaw]

How do you feel about this PR now? Please let me know if anything else needs to be addressed in this PR.

I think it is good and appears to complete:

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

If so, then I think it is ready?

[NimaSarajpoor]

Perfect! I will check one last time just in case and then merge.

[NimaSarajpoor]

@seanlaw
Thanks for all the inputs! Going to merge this now.