This PR addresses issue #592. In this PR, we want to extend the function stump and the related ones so that it returns Top-K Nearest Neighbors Matrix Profile (i.e. the k smallest values for each distance profile and their corresponding indices).

What I have done so far:

• Wrote a naïve implementation of stump_topk
• Provided a new unit test for function stump

Notes: (1) np.searchsort() is used in the naive.stump_topk. Another naive approach can be traversing distance matrix row-by-row, and using np.argpartition() followed by `np.argsort()'.

(2) I think considering the first k columns in the output of stump_topk is cleaner, and also easier for the user to get access to those top-k values.

I can/will think about a clean way to change the structure of output just before returning it so that the first four columns becomes the same for all k. However, I think that makes it hard to use the output later in other modules.

Add To-Do List

• [x] Add parameter k to stump
• [x] Add parameter k to stumped
• [x] Add parameter k to gpu_stump
• [x] Add parameter k to scrump
• [x] Add parameter k to stumpi
• [ ] Run published tutorials to see if they still work with no problem

Side notes:

• Add parameter k to non-normalized versions.
• Check docstrings for functions that require matrix profile P as input. For instance, in stumpy.motifs, we may say something like... "P is (top-1) 1D array matrix profile"

Discovering motifs under uniform scaling

From #85, I add Discovering motifs under uniform scaling tutorial from 3.1 of The Matrix Profile Top Ten paper. Add docs/Tutorial_Matrix_Profile_Top_Ten.ipynb

It still needs modification to plot, markdown, coding...., but please have a look through and let me know if you find problem.

Pull Request Checklist

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

• [x] Fork, clone, and checkout the newest version of the code
• [x] Create a new branch
• [x] Make necessary code changes
• [ ] Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• [ ] Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• [ ] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• [ ] Run black . in the root stumpy directory
• [ ] Run flake8 . in the root stumpy directory
• [ ] Run ./setup.sh && ./test.sh in the root stumpy directory
• [ ] Reference a Github issue (and create one if one doesn't already exist)

In this PR, we would like to implement MERLIN algorithm that discovers discords of arbitrary length. Although the MATLAB implementation of the faster version of MERLIN is available on MERLIN: SUPPORT, we, for now, implement the original version as proposed in the MERLIN.

What I have done so far:

• Add Introduction
• Explain the advantage / disadvantage of MatrixProfile for discovering discords
• Explain two core ideas of MERLIN algorithm
• Implement the first phase of DRAG (DRAG is the first part of MERLIN)

NOTE: (1) I already implemented MERLIN algorithm and enhanced it a little bit. Both MERLIN and STUMPY: MatrixProfile gives EXACTLY the same output regarding the discords indices and their distances to their NN. The discord NN indices are the same in almost all cases.

(2) In terms of performance, MERLIN outperforms STUMPY for long time series.

Pull Request Checklist

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

• [x] Fork, clone, and checkout the newest version of the code
• [x] Create a new branch
• [x] Make necessary code changes
• [x] Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• [x] Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• [x] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• [x] Run black . in the root stumpy directory
• [x] Run flake8 . in the root stumpy directory
• [x] Run ./setup.sh && ./test.sh in the root stumpy directory
• [x] Reference a Github issue (and create one if one doesn't already exist)

Pull Request Checklist

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

• [x] Fork, clone, and checkout the newest version of the code
• [x] Create a new branch
• [x] Make necessary code changes
• [x] Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• [x] Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• [x] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• [x] Run black . in the root stumpy directory
• [x] Run flake8 . in the root stumpy directory
• [x] Run ./setup.sh && ./test.sh in the root stumpy directory
• [x] Reference a Github issue (and create one if one doesn't already exist) #552

Added first draft of the Annotation Vectors Tutorial, requested in issue #177, for feedback. The tutorial includes all 4 examples explored in the Matrix Profile V paper (https://www.cs.ucr.edu/~hdau001/guided_motif_search/). Planning to remove 2/3 out of the 4 examples. Let me know which ones I should keep/remove.

Pull Request Checklist

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

• [x] Fork, clone, and checkout the newest version of the code
• [x] Create a new branch
• [x] Make necessary code changes
• [x] Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• [x] Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• [x] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• [x] Run black . in the root stumpy directory
• [x] Run flake8 . in the root stumpy directory
• [x] Run ./setup.sh && ./test.sh in the root stumpy directory
• [x] Reference a Github issue (and create one if one doesn't already exist)

I gave the motif discovery module a go (issue #184), what do you think? (Ignore all changes outside motifs.py and test_motifs.py, as we discuss this in another PR)

There are still more tests to be done, but I'm not sure how to approach it yet. Ideally I'd like a naive_k_motifs method.

This is a work-in-progress pull request addressing issue #248, i.e. implementation of a method that reduces cache misses while computing the Matrix Profile using stump.

How it works right now

Let us first briefly explore how stump works right now and what's wrong with it.

Consider a time series of length 20 and a window of 6. Then we know the distance matrix will be a 15 x 15 matrix (20 - 6 + 1 = 15). Furthermore, the exclusion zone will be ±2 (np.ceil(6 / 4) = 2). Thus the distance matrix will look something like this:

The grey area is the exclusion zone.

Currently, if this time series and window are passed into stump, each of the diagonals in the matrix will be processed one by one:

The traversal will be along the yellow dividers.

The benefit with this method is the benefit we get from computing by diagonals. Computing each distance is a lot simpler given the distance from the previous index in the diagonal. For eg, once we compute the distance for point (3, 0), we can use it to calculate the distance for point (4, 1).

The issue with this is, the cache is more likely to store points that are closer to each other, and so traversing an entire diagonal is likely to cause too many cache misses.

How we want it to work

While traversing diagonally is certainly the right approach, traversing entire diagonals will negatively affect performance. Instead, we want to split the matrix into smaller tiles, each of which will be divides into groups of diagonals.

Suppose the tile size is set to 4 and the number of diagonals to be iterated together is set to 2. Then, the distance matrix will be divided like this:

The matrix is split into 4 x 4 tiles (except near the edges). Tiles with blue borders indicate that there are distances to traverse within those tiles, tiles with red borders indicate that we should skip those tiles as there are no distances to traverse within them. Further more, each of the diagonals within the blue tiles are split into groups of 2 (or less). Each diagonal groups are traversed together.

For example, if we were to traverse the tile colored light pink, we would do so in the following order:

(8, 4)
(8, 5)
(9, 5)
(9, 6)
(10, 6)
(10, 7)
(11, 7)
(9, 4)
(10, 4)
(10, 5)
(11, 5)
(11, 6)
(11, 4)

What's done so far

The changes in the alvii147/cache_optimization branch is the first step to accomplishing this.

Two new functions were added to core.py, _get_tiles and _get_tiles_ranges, and a new module was added named stump_tiles.py.

core._get_tiles

_get_tiles(m, n_A, n_B, diags, tile_size)

Parameters
----------
m : int
    Window size

n_A : ndarray
    The length of time series `T_A`

n_B : ndarray
    The length of time series `T_B`

diags : ndarray
    The diagonal indices of interest

tile_size : int
    Maximum length of each tile

_get_tiles can be used to split the distance matrix into tiles of a given size. It returns tiles, an array of 2-column arrays, with each element representing each tile. The contents of each 2-column array are described below.

2-column array contents:

[ y_offset                x_offset            ]
[tile_height              tile_width          ]
[lower_diagonal_index     upper_diagonal_index]
[tile_weight              tile_validity       ]

x_offset, y_offset: position of the tile with respect to the distance matrix
tile_height, tile_width: dimensions of the tile
lower_diagonal_index, upper_diagonal_index: range of diagonals to be traversed within the tile
tile_weight: number of distances to be computed within the tile
tile_validity: 0 if no distances within tile are to be computed, 1 otherwise

The first 4 tiles of the distance matrix considered above should return the following result:

[[ 0  0]
 [ 4  4]
 [ 3  4]
 [ 1  1]]

[[ 0  4]
 [ 4  4]
 [-1  4]
 [13  1]]

[[ 0  8]
 [ 4  4]
 [-3  4]
 [16  1]]

[[ 0 12]
 [ 4  3]
 [-3  2]
 [12  1]]

...

core._get_tiles_ranges

_get_tiles_ranges(tiles, n_threads)

Split given `tiles` into `n_threads`.

Parameters
----------
tiles : ndarray
    An array of two column arrays representing each tile

n_threads : int
    Number of threads to split the tiles into

Returns
-------
tiles_ranges : ndarray
    A two-dimensional array representing tile indices for each thread

_get_tiles_ranges can be used to split the tiles that we obtain from _get_tiles into n_threads given threads. This is done by sorting the tiles in descending order of tile_weight, then assigning them one by one to each thread. Once we reach the end of the thread, we keep assigning the tiles to the threads in reverse order.

For example, suppose we have tiles of weights in descending order [7, 6, 5, 3, 2, 1] and we have 3 threads to divide these into. We start by assigning 7, 6, and 5 to the first three threads respectively. Then the third thread gets 3, second thread gets 2, and first thread gets 1. So the division between the threads would look something like this:

First thread: [7 1]
Second thread: [6 2]
Third thread: [5 3]

This way, the workload is most equally divided.

If we were to run _get_tiles on the previously described example, and then divide those tiles into 3 threads using _get_tiles_ranges, we would obtain the following array, where each row are the indices of the tiles assigned to each thread:

[[ 2, 11, 10, 13, 12, -1]
 [ 6,  3,  5, 14,  9, -1]
 [ 1,  7,  0, 15,  8,  4]]

NOTE: -1 is just a blank slot.

stump_tiles

The stump_tiles module is a modified version of stump. While stump calls _stump which calls _compute_diagonal, stump_tiles calls _stump_tiles, which calls _compute_tiles.

_stump_tiles will call core._get_tiles and core._get_tiles_ranges to get tiles data and pass each group of tiles to _compute_tiles in a single thread. _compute_tiles will process the given tiles and divide each tile into "chunks" of tiles that are to be processed together.

The default tile size, STUMPY_TILE_SIZE and the default number of diagonal to traverse together, STUMPY_DIAGONALS_CHUNK_SIZE are current defined in config.py.

...
STUMPY_TILE_SIZE = 1000
STUMPY_DIAGONALS_CHUNK_SIZE = 100

Testing

The only testing module added so far is test_stump_tiles.py, which is copied over and modified from test_stump.py. Use the following command to run them:

export NUMBA_DISABLE_JIT=1
export NUMBA_ENABLE_CUDASIM=1
py.test -x -W ignore::RuntimeWarning -W ignore::DeprecationWarning tests/test_stump_tiles.py

unset NUMBA_DISABLE_JIT
unset NUMBA_ENABLE_CUDASIM
py.test -x -W ignore::RuntimeWarning -W ignore::DeprecationWarning tests/test_stump_tiles.py

All test cases in test_stump_tiles are currently PASSING.

Performance

Currently stump_tiles performs very poorly on relatively smaller time series, and marginally faster on relatively larger ones. Outlined below are the performances comparisons stump and stump_tiles for different time series lengths.

Time Series Length | Window Size m | stump Performance (seconds) | stump_tiles Performance (seconds) | Percentage Improvement --- | --- | --- | --- | --- 1000 | 10 | 0.016 | 0.016 | 0.0% 2500 | 10 | 0.016 | 0.031 | -100.0% 5000 | 10 | 0.06 | 0.078 | -30.77% 10000 | 10 | 0.269 | 0.219 | 18.61% 25000 | 10 | 1.533 | 1.366 | 10.91% 50000 | 10 | 6.262 | 5.206 | 16.87% 100000 | 10 | 25.189 | 20.589 | 18.27% 200000 | 10 | 104.818 | 84.124 | 19.74% 250000 | 10 | 164.197 | 130.789 | 20.35%

You may use the following script for testing performance:

import stumpy
import numpy as np
import numpy.testing as npt
import time

def performance_test(m, ts_len):
    ts = np.random.rand(ts_len)

    print('\nTime series length:', ts_len)
    print('Window length:', m)

    time_elapsed = time.time()
    stump_mp = stumpy.stump(ts, m=m)
    time_elapsed = time.time() - time_elapsed
    print('stump performance (seconds):', time_elapsed)

    time_elapsed = time.time()
    stump_tiles_mp = stumpy.stump_tiles(ts, m=m)
    time_elapsed = time.time() - time_elapsed
    print('stump_tiles performance (seconds):', time_elapsed)

    npt.assert_almost_equal(stump_tiles_mp, stump_mp)

if __name__ == "__main__":
    stumpy.stump(np.random.rand(1000), m=200)
    stumpy.stump_tiles(np.random.rand(1000), m=200)

    parameters = (
        (10, 1000),
        (10, 2500),
        (10, 5000),
        (100, 10000),
        (100, 100000),
        (100, 200000),
        (100, 250000),
    )
    for param in parameters:
        performance_test(*param)

Note that the performance may vary with the values of STUMPY_TILE_SIZE and STUMPY_DIAGONALS_CHUNK_SIZE, as well as with different systems.

Okay, but what now?

There's still a lot to do, here are a few things that come to mind.

Time different components

Use the time module to time different functions that we call, namely _get_tiles, _get_tiles_ranges and _compute_tiles, find out where exactly we lose time. If any of the components take time that is exponentially increasing with the number of distances, that's a good place for possible optimization.

Optimize _get_tiles_ranges

I am almost certain _get_tiles_ranges can be simplified. Currently this function sorts the tiles by weight, which is very expensive and may not be worth it. One thing to remember, the results of _get_tiles_ranges need not be perfect, the tiles don't need to be perfectly equally split among the threads. There's no point in dividing it perfectly if we spend a lot of time doing it. A roughly equal division is also good enough if it saves time.

Optimize _compute_tiles

I will also not be surprised if _compute_tiles isn't perfectly optimized. There may be room for improvement with the way the iterations are done.

Find good values for STUMPY_TILE_SIZE and STUMPY_DIAGONALS_CHUNK_SIZE

This is a tough one, since it's system dependent. Section 4.2.2 of the Time Series Analysis with Matrix Profile on HPC Systems Paper may help with this.

Make the code NumPy and Numba friendly

Some places may use native Python functions, such as range. These should be replaced with NumPy equivalents, such as np.arange.

This uses the mtDNA example. Ostinato transcribed from table 2 in paper. Performance is really slow.

resolves #277, resolves #278

Pull Request Checklist

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

• [x] Fork, clone, and checkout the newest version of the code
• [x] Create a new branch
• [x] Make necessary code changes
• [x] Install black (i.e., python -m pip install black or conda install black)
• [x] Install flake8 (i.e., python -m pip install flake8 or conda install flake8)
• [x] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install pytest-cov)
• [x] Run black . in the root stumpy directory
• [x] Run flake8 . in the root stumpy directory
• [x] Run ./setup.sh && ./test.sh in the root stumpy directory
• [x] Reference a Github issue (and create one if one doesn't already exist)

Given that there is interest for 1-D discord/motif search, it may be useful to have M-dimensional motif/discord search capabilities.

According to the last paragraph in Section IV E of the mSTAMP paper:

In the case where multiple semantically meaningful k-dimensional motifs are presented in the multidimensional time series (e.g., Fig. 5), we can just interactively apply the MDL-based method to discover the motif. There are two steps in each iteration: 1) apply the MDL-based method to find the k-dimensional motif with the minimum bit size and 2) remove the found k-dimensional motif by replacing the matrix profile values of the found motif (and its trivial match) to infinity. If we apply the two steps above to the time series shown in Fig. 5, the 3-dimensional motif would be discovered in the first iteration, and the 2-dimensional motif would be discovered in the second iteration. In terms of the terminal condition for the iterative method, it can be either be an input for the user or a more advanced technique could be applied. Due to space limitations, we will have to leave the discussion on termination condition to future works. An example of applying such iterative algorithm on real-world physical activity monitoring time series is shown in Section V.E.

So, this means that it is possible to find different motifs at the same index location but with a different subset of dimensions!

Added my first attempts at optimizing stump using unsigned integer indexing. I added a temporarity stumpy.stump_uint module, which (if successfully optimized) may replace stump. But I haven't had too much luck, maybe someone could help me out here. Here's what I've done so far and its results:

Time series length | window size | stump mean execution time | stump_uint mean execution time | percentage improvement --- | --- | --- | --- | --- 1000 | 10 | 0.00966 | 0.01066 | -10.35% 2500 | 10 | 0.02866 | 0.02933 | -2.323% 5000 | 10 | 0.079 | 0.095 | -20.252% 10000 | 100 | 0.27566 | 0.30233 | -9.674% 100000 | 100 | 26.511 | 27.163 | -2.458% 100000 | 5000 | 25.706 | 27.825 | -8.245% 200000 | 100 | 112.324 | 108.871 | 3.074% 200000 | 5000 | 105.614 | 104.584 | 0.976% 250000 | 100 | 169.058 | 169.176 | -0.07%

Doesn't look like there's much of a difference just yet, but I'm nearly certain this is just me optimizing the wrong parts. Let me know if you have feedback, feel free to give this a shot yourself.

See this discussion

Currently, STUMPY is poised to accept/handle a ray distributed cluster to compute a matrix profile. This would be purely experimental but it should be possible with a little work.

• [ ] stumped/aamped
• [ ] mstumped/maamped
• [ ] stimped/aamp_stimped
• [ ] mpdisted/aampdisted
• [ ] Add ray example to docstrings

It seems there are a few inconsistencies in the naming of variables and the way they are being passed to functions.

(i) In stump, we have this: https://github.com/TDAmeritrade/stumpy/blob/f078424350ab97d02a517da9e29b254309b27574/stumpy/stump.py#L634-L650

So, T_A corresponds to Q-related variables, and T_B corresponds to T-related variables. However, we see something different in gpu_stump: https://github.com/TDAmeritrade/stumpy/blob/f078424350ab97d02a517da9e29b254309b27574/stumpy/gpu_stump.py#L577-L578

(ii) Also, the module stump itself shows some sort of inconsistency: https://github.com/TDAmeritrade/stumpy/blob/f078424350ab97d02a517da9e29b254309b27574/stumpy/stump.py#L678-L695

Since we pass the arguments as _stump(T_A, T_B, ...), we should expect to see the same order in the rest of the arguments. However, this is not true here. As you may have noticed already, M_T (which corresponds to T_B) is located before μ_Q (which corresponds to T_A).

In my opinion, these inconsistencies occur because we, unintentionally, confuse ourselves by using the T_A, T_B naming and Q and T naming. I think we should usually consider T_A as the sequence that contains queries (hence Q), and T_B is T.

So, if we do:

Q, μ_Q, σ_Q = core.preprocess(T_A, m)
T, M_T, Σ_T = core.preprocess(T_B, m) 

We will be consistent and we will not be confused when we try to be consistent regarding the order of arguments in a function. (Btw, I just provided two examples. I haven't checked all modules yet)

This PR is an attempt to handle constant subsequences in different modules by using the outcome of core.rolling_isconstant(T, m).

• [ ] check docstrings
• [ ] check order of arguments (in private / public APIs)
• [ ] make sure GPU modules pass all tests when CUDA is enabled
• [ ] make sure all notebooks are not broken

It seems that this

https://github.com/TDAmeritrade/stumpy/blob/8da71ecec45e947cf448eb3dba8a75c98f5dbfe2/stumpy/scrump.py#L209-L211

is never executed, and even if it is executed, there will be no change!

The reason is that: (1) all elements of P_squared are initially set to np.inf and their corresponding indices in I are initially set to -1. (2) P_squared is sorted ascendingly. Thus, after finding at least one finite distance, the value P_squared[thread_idx, i, 0] becomes finite. This is because the smallest distance will be at index 0.

Therefore, after finding the first finite distance between subseq i and one of its neighbor, the value at P_squared[thread_idx, i, 0] becomes finite. Note that if we never find a finite distance for subseq i, then all values in P_squared[thread_idx, i, :] remain np.inf and their corresponding indices remain -1. Hence, if-block results in no change.

We have a similar situation in scraamp. see Lines 194-197

The only drawback is that P_squared is a variable that is passed to the function _prescrump. Hence, we may need to improve the docstring to let user know that the default value in P_squared is np.inf, and its corresponding index in I is -1.

SWAMP/ DTW MP discussion

This is a pull-request for information exchange for #466. About SWAMP, DTW MP, LB MP and related algorithm to reduce calculation cost.

Pull Request Checklist

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

• [ ] Fork, clone, and checkout the newest version of the code
• [ ] Create a new branch
• [ ] Make necessary code changes
• [ ] Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• [ ] Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• [ ] Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• [ ] Run black . in the root stumpy directory
• [ ] Run flake8 . in the root stumpy directory
• [ ] Run ./setup.sh && ./test.sh in the root stumpy directory
• [ ] Reference a Github issue (and create one if one doesn't already exist)