Fixed #1111 Add branching logic for sdp

[NimaSarajpoor]

See #1111

Pull Request Checklist

Below is a simple checklist but please do not hesitate to ask for assistance!

• Read our Contributing Guide
• Referenced a Github issue (or create one if one doesn't already exist)
• Read and reviewed all of the comments in the Github issue that you've referenced (along with cross-referenced issues/pull requests) to ensure that the issue still requires a pull request
• Checked that the issue has not already been assigned to anybody else or is already being addressed in another pull request
• Left a meaningful comment on the original Github issue to discuss the detailed approach for your contribution and received confirmation from the maintainers before proceeding with this pull request
• Forked, cloned, and checked out the newest version of the code
• Created a new branch
• Made necessary code changes
• Installed black (i.e., python -m pip install black or conda install -c conda-forge black)
• Installed flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• Installed pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• Ran black --exclude=".*\.ipynb" --extend-exclude=".venv" --diff ./ in the root stumpy directory
• Ran flake8 --extend-exclude=.venv ./ in the root stumpy directory
• Ran ./setup.sh dev && ./test.sh in the root stumpy directory and ensured that all tests are passing locally
• Check this box if AI code assistance was used to generate
  • less than 25% of the code in this pull request
  • 25-50% of the code in this pull request
  • more than 50% of the code in this pull request

Please do not commit any code to avoid/circumvent a failing test and, instead, engage in a discussion (below) to determine the best course of action.

Only request a review after the checklist above is fully completed!

---

This is to address PR 4 as described in #1118 (comment). I've copied the corresponding notes below:

Notice that we haven't utilized ANY of the _X_sliding_dot_product anywhere in the code except for _njit_sliding_dot_product but that hasn't altered any of the code in a meaningful way as it was purely a refactor. Only at the very end in this PR do we start working on allowing the other _X_sliding_dot_product to be used. This is where we focus on the branching logic!

Note that instead of core.sliding_dot_product = sdp._sliding_dot_product, we can redefine core.sliding_dot_product so that it performs checks (e.g., check whether pyfftw is installed BEFORE callings sdp._sliding_dot_product). This way, sdp._sliding_dot_product can make assumptions that the user knows what they're doing while core.sliding_dot_product is slower but is super safe to use (hence it is public)!

[gitnotebooks]

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

[NimaSarajpoor]

As discussed previously in this comment, I am planning to replace core.sliding_dot_product with something like this:

def sliding_dot_product(Q, T, boundaries):      
      m = len(Q)
      n = len(T)
      sdp = None
      for boundary in boundaries:
            (LB_m, UB_m), (LB_n, UB_n), sdp_func = boundary
             if LB_m <= m < UB_m and LB_n <= n < UB_n:
                  return sdp_func(Q, T)
           
      return sdp_func_default(Q, T)

[seanlaw]

def sliding_dot_product(Q, T, boundaries):      
      m = len(Q)
      n = len(T)
      sdp = None
      for boundary in boundaries:
            (LB_m, UB_m), (LB_n, UB_n), sdp_func = boundary
             if LB_m <= m < UB_m and LB_n <= n < UB_n:
                  return sdp_func(Q, T)
           
      return sdp_func_default(Q, T)

Doesn't this mean that every time somebody passes in the boundaries then the sets of bounds get loaded every time you call this function? This logic seems suspicious but I'll wait until you implement it to assess further.

[NimaSarajpoor]

@seanlaw

Doesn't this mean that every time somebody passes in the boundaries then the sets of bounds get loaded every time you call this function?

Right! A user often uses a certain system. So, users probably would like to set the boundaries once and then just compute the sliding dot product by passing Q and T. I thought about three choices:

(1) Config variable: But... we are not actually using boundaries anywhere else in the code.
(2) Class, with __call__ method: But... similar to pyfftw in sdp.py, attributes will be for internal use.
(3) Closure: the approach we took for pyfftw in sdp.py

I decided the last approach, i.e. closure, and I've pushed the code. What do you think?

[seanlaw]

func not fucn

[seanlaw]

"chooseS"

[seanlaw]

There's just something about this that doesn't sit right with me. Here's what I am seeing:

1. I get that when we do sliding_dot_product = make_sliding_dot_product() (let's refer to this left side as sdp1), you are setting sliding_dot_product to the sliding_dot_product function that is being returned from inside of the closure (let's refer to this inside-of-the-closure-function as sdp2)
2. However, the contents of the sliding_dot_product function inside of the closure (sdp2) are never executed until you call sdp1.
3. This means that if you have many Q and T pairs, it is performing boundary comparisons every time rather than performing an O(1) lookup
4. Similarly, if you have many boundaries to check, then this would be slow

So, while the logic is sound, this point:

Doesn't this mean that every time somebody passes in the boundaries then the sets of bounds get loaded every time you call this function?

is still relevant. Yes, the sets of boundaries are supplied/embedded once into sdp2 BUT each call to sdp2 triggers a slow iterative search. Is there a cheap time AND space efficient way to create an O(1) function lookup given m and n? This for loop isn't the right approach. If we have an O(1) function lookup then it's really lightweight and I wouldn't mind going the config route because 99.9% of users will use our default.

Having said that, are we going to have a way for users to run a function and come up with boundaries and functions that are best suited for their hardware?

[NimaSarajpoor]

So, while the logic is sound, this point:

Doesn't this mean that every time somebody passes in the boundaries then the sets of bounds get loaded every time you call this function?

is still relevant. Yes, the sets of boundaries are supplied/embedded once into sdp2 BUT each call to sdp2 triggers a slow iterative search.

Right, I misunderstood your point before.

Is there a cheap time AND space efficient way to create an O(1) function lookup given m and n? This for loop isn't the right approach. If we have an O(1) function lookup then it's really lightweight

Let's replace that for-loop with choose_sdp_func, a function that gets m and n as inputs and return a sdp function as output, and also let's get rid of closure. So, we can have a regular function like this:

def sliding_dot_product(Q, T, choose_sdp_func=None):
    m = len(Q)
    n = len(T)
    if m == n:
        return np.array([np.dot(Q, T)])
 
    if choose_sdp_func is None: 
        # STUMPY DEFAULT
    
    sdp_func = choose_sdp_func(m, n)  # branching logic
    return sdp_func(Q, T)

Is it a bad design? This allows [advanced] users to use different ways to implement branching logic for their choose_sdp_func function:

• A dictionary, keyed by (m,n), and the values are sdp functions. It is O(1)
• Or, if-else logic
• Or, boundaries.
• etc

We can always add a caching mechanism on top of the function choose_sdp_func to make sure it returns output in O(1) for a repeated input (m, n). Regarding STUMPY DEFAULT for choose_sdp_func, we can go with a simple if-else logic.

As a side, regarding the following question you raised before in another comment:

the question remains, how do we allow the flexibility to swap out or modify the if/else logic in that case as well? We should think about this..

I think the proposal should address that

---

Having said that, are we going to have a way for users to run a function and come up with boundaries and functions that are best suited for their hardware?

Finding "boundaries" is tricky. I have a narrower scope (and, tbh, simpler solution 😅) in mind... a function that gets (m, n) and a list of sdp functions as input, and returns the best one for the provided (m, n).

[NimaSarajpoor]

Finding "boundaries" is tricky

Okay, I think it might be tricky for most users too! So, maybe they appreciate if we just come up with a function that helps them identify boundaries approximately. For instance, we just explore cases where m and n are exact power of two.

[seanlaw]

I am mostly interested in what the STUMPY_DEFAULT choose_sdp_func would look like and how efficient it would be?

[NimaSarajpoor]

I wanted to say a simple if-else... but I think it can get ugly later. I think we should do what you proposed before in stumpy-dev/sliding_dot_product#18 (comment), a 2D lookup table. And, I think there is no longer a need for choose_sdp_func !!

def make_sliding_dot_product(
    lengths=None, sdp_lookup=None, sdp_functions=None, sdp_default=None
):
    if sdp_default is None:
        sdp_default = sdp._convolve_sliding_dot_product

    if sdp_functions is None:
        sdp_functions = [
            sdp._njit_sliding_dot_product,
            sdp._pocketfft_sliding_dot_product,
            sdp._convolve_sliding_dot_product,
        ]

        # append if pyfftw is available.
        # len(sdp_functions) can change here but it will be 
        # taken care of later in the code
        if sdp.PYFFTW_IS_AVAILABLE:  # pragma: no cover
            sdp_functions.append(sdp._pyfftw_sliding_dot_product)
      

    if lengths  is None or sdp_lookup is None:
        lengths = set()
        for i in range(3, 2**28 + 1):
            lengths.add(next_fast_len(i, real=True))
        lengths = np.sort(list(lengths), dtype=np.uint32) 
        sdp_lookup = np.empty((len(lengths), len(lengths)), dtype=np.uint8)

        # Filling sdp_lookup
        idx_m = np.searchsorted(lengths, config.STUMPY_NJIT_SDP_Q_LENGTH, side='right')
        sdp_lookup[:idx_m, :] = 0
        sdp_lookup[idx_m:, :] = 4  
    
    # when sdp_lookup contains values that are not in the range of indices of `sdp_functions`
    mask = sdp_lookup >= len(sdp_functions)
    if np.any(mask):
        # raise warning
        sdp_lookup[mask] = -1  # fallback indicator

The # Filling sdp_lookup part has similar issue to what we raised for if-else. I mean... we need to revise it every time as we find better functions for different (len(Q), len(T)) pairs.

So, I feel I need to think about...

are we going to have a way for users to run a function and come up with boundaries and functions that are best suited for their hardware?

In this PR.

[seanlaw]

I think we should do what you proposed before in stumpy-dev/sliding_dot_product#18 (comment), a 2D lookup table

Even using np.searchsorted is O(nlogn), which is undesirable. It would be best to find something that is O(1) if possible.

Filling sdp_lookup

I wouldn't do this inside of make_sliding_dot_product. Instead, I would move the filling code outside of make_sliding_dot_product and this filling code would rarely be used (i.e., a human would manually call it every so often to generate the lookup). Then, the generated lookup table would/should be hardcoded somewhere in a .py file (maybe as a function?). In make_sliding_dot_product, it would simply import the lookup-table-function (but not the "generation" function).