MEXXEXX DM

[seanlaw]

@mexxexx said:

However, if you feel like you need support please feel free to ask, I'd like to help add some new functionality.

[seanlaw]

@mexxexx I’ve really enjoyed working with you and I can certainly use and welcome your help! Let's use this issue to have an ongoing an open dialogue (like a long message thread).

Depending on your availability, I could definitely use some help:

• adding new tutorials
• maintaining NaN/inf code (and adding it to new code) as I still get confused when I don’t touch it in a while
• completing existing issues (or breaking them into smaller ones so you feel more comfortable contributing)
• basically anything that you’d be interested in adding (i.e., new functionality) would be welcome 😊

My goal for STUMPY over the next while will not be to add features. Instead, I would like to focus on making the computation of the matrix profile rock solid and consistent across the code base and then to communicate the many use cases through illustrative (code + visual) examples. Essentially, I want to continue to strengthen the current foundation so that others can build off of it in the future. I would like to keep the scope tight but I trust your judgement.

Your assistance would be greatly appreciated!

[mihailescum]

Hi! I think its a good approach to really have the current functionality working and adding tutorials.

I can definitely assist in maintaining NaN and inf values, also I can look at it for scrimp, once you have restructured the code. Also for the other issues of course, but it is always a bit difficult to completely understand a new topic.

[seanlaw]

Awesome! Btw, I was able to add NaN/inf support to SCRIMP and all unit tests by following what you had done with stumpy.stump. It was rather painless thanks to your earlier work!

In fact, I converted SCRIMP into a generator so that one can quickly and easily iterate and incrementally (and randomly) update the matrix profile. The next step is to add preSCRIMP (and then we can release SCRIMP++ as a new feature).

After that, I’d like to add NaN/inf support to GPU_STUMP and then focus on adding tutorials. Again, you are welcome to take on any of this stuff and I will support you where you need me to! So don’t feel like you are going at it alone.

[seanlaw]

Now that I have SCRIMP (not preSCRIMP and SCRIMP++) completed, I have learned a lot about how to traverse the matrix across the diagonals (instead of row-wise). Indeed stumpy.scrimp is about the same speed (for self-joins) as stumpy.stump. Additionally, I have discovered a lot about how to parallel program using Numba’s prange when you don’t have direct access to each thread and so this will come in handy later.

The diagonal processing is fundamental to improving the speed but the more significant speed up (and added numerical stability) comes from the new way to compute the pairwise distances using correlation/covariance. From my rough estimates, this is about a 5x speedup but it’ll be worth profiling. I think this is one place that you might be able to help and that you may find interesting.

[seanlaw]

Also for the other issues of course, but it is always a bit difficult to completely understand a new topic.

I don’t dispute this but, since I trust you, I also don’t think that it is bad to add individual byte sized functions (or even isolated modules) to the code base that is based on the algorithms/pseudocode from an existing paper. I think SCRIMP in its current form is a good example of this where I started a module and it was a single function with tests. Then I slowly added parallelization and finally NaN/inf support. Eventually, we’ll have preSCRIMP and then SCRIMP++ will be done. But SCRIMP is the foundation. While we obviously want to avoid starting 10 things at the same time, I think it’s okay to have some “work in progress” code as long as there are unit tests to go along with it. I don’t consider a function to be “public” until it is added our docs (and I’d prefer having a tutorial over that). I have been known to start a module with tests and commit it and then decided to redesign it because I figured out a better way. No big deal.

I see something similar with the noisy data issue. Personally, I would add the distance calculation first (along with appropriate tests), assume that the choice is boolean for now, and then eventually add a threshold/float/array. I’ll try not to get into the weeds here with this one but hopefully you get my point that we can iterate together.

[mihailescum]

Awesome! Btw, I was able to add NaN/inf support to SCRIMP and all unit tests by following what you had done with stumpy.stump. It was rather painless thanks to your earlier work!

Great, I'm happy to hear this!

In fact, I converted SCRIMP into a generator so that one can quickly and easily iterate and incrementally (and randomly) update the matrix profile. The next step is to add preSCRIMP (and then we can release SCRIMP++ as a new feature).

Cool! If you want, we can work on this together, e.g. I can doublecheck anything you propose. I found the paper rather confusing, so I'm not sure if I could implement prescrimp myself.

After that, I’d like to add NaN/inf support to GPU_STUMP and then focus on adding tutorials. Again, you are welcome to take on any of this stuff and I will support you where you need me to! So don’t feel like you are going at it alone.

I can write a tutorial on how to use scrimp if you want. I dont have a GPU so I cannot help with this part.

[mihailescum]

Now that I have SCRIMP (not preSCRIMP and SCRIMP++) completed, I have learned a lot about how to traverse the matrix across the diagonals (instead of row-wise). Indeed stumpy.scrimp is about the same speed (for self-joins) as stumpy.stump. Additionally, I have discovered a lot about how to parallel program using Numba’s prange when you don’t have direct access to each thread and so this will come in handy later.

That is good, I also think that implementing algorithms is a great way to learn a language/specifics.

The diagonal processing is fundamental to improving the speed but the more significant speed up (and added numerical stability) comes from the new way to compute the pairwise distances using correlation/covariance. From my rough estimates, this is about a 5x speedup but it’ll be worth profiling. I think this is one place that you might be able to help and that you may find interesting.

Definitely, you just have to clarify what you mean exacltly, maybe in a separate issue.

[seanlaw]

@mexxexx Here are some links:

Google Colab GPUs
Suffix Array Library
Suffix Array YouTube Playlist

[mihailescum]

So about definition 13 in the mstomp profile: for the multidimensional profile you sort the distance profiles column wise, and by this lose the information in which order the distances where combined (i.e. for the two dimensional distance profile the two lowest distances are averaged, but if we have more than two dimensions in our time series, we lose the information about which dimensions played a role here). So the idea of Def 13 is to keep track of it with something called the matrix profile subspace.

Say that you have a three dimensional series and compute the following distances

[1, 2, 3]
[2, 3, 1]
[1, 3, 2]

Then, after sorting you get

[10, 20, 10]
[10, 30, 20]
[20, 30, 30]

And averaging (D_prime_prime)

[10, 20, 10]
[10, 2.5, 15]
[16.6, 26.6, 20]

So for this query the profile would be [10, 10, 16.6]. Now, if you want to find the closest subsequence in two dimensions, how do you know in which dimensions you have to look? You know that the closest distance is 10, but is it in dimensions (1,2), (2,3) or (1,3)?

For our example the MP subspace would be (the column wise sort order9

[1, 1, 2]
[3, 2, 3]
[2, 3, 1]

Hence, as the MP value of 10 os at index zero, we know that for two dimensions we have to look at dimensions (1,3) combined.

[mihailescum]

This isn't the greatest explanation, but maybe it gives a hint of why we need this information.

[seanlaw]

Okay, thank you. This rings a bell. So, it sounds like we'd need to replace the D.sort(axis=0) with sorted_idx = np.argsort(D, axis=0, kind='mergesort') (we need to choose mergesort as the default quicksort is considered unstable/nondeterministic), then store the indices that are associated with each number of dimensions (i.e., one index for D1, two indices for D2, and so on and this applies to each column), and then do D = D[sorted_idx]. I vaguely remember this when I first implemented mstump and I avoided the subspace stuff because I couldn't come up with a good way to store this information inside of a NumPy array. If I'm doing my math correctly, for a D-dimensional time series, the output would be D! * (T.shape[0] - m + 1? Hmm, upon reading it further, it looks like the user would have to choose the k dimensions ahead of time and then only return the subspace, which would be of size k * [T.shape[0] - m + 1. But the million dollar question is "what should k be"?

I think that once we can figure the output then the rest shouldn't be too hard (as you pointed out).

Thanks for explaining it! Btw, should the originally computed distances be?

[10, 20, 30]
[20, 30, 10]
[10, 30, 20]

[seanlaw]

Okay, it's not the best option but apparently you can do:

import numpy as np
subspace = np.ones((5, 10), dtype=np.ndarray) * -1
subspace[2, 0] = np.array([2, 4, 1])

which gives:

array([[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
       [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
       [array([2, 4, 1]), -1, -1, -1, -1, -1, -1, -1, -1, -1],
       [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
       [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1]], dtype=object)