Finding multiple motifs with mStump

[SaVoAMP]

Hey Sean,

I am currently working with Matrix Profiles for my bachelor thesis where I am trying to find motifs on multivariate sensor data. Therefore I am interested in finding motifs with your mStump algorithm.
Here is a screenshot of an example data set where the labels are the motifs I am looking for. The data was labeled manually and in the scope of my thesis I am trying to find out if it is possible to automate part of this process with the use of Matrix Profiles.

I have already worked threw your tutorial (they were all very helpful for my understanding because of their "hands-on" approach) and plotted some morif pairs. But as you can see on the screenshot there are multiple motifs that I need to find.
While reading your discussion #392 I was wondering if it suffices to change the number of motifs I want to find in this line of code:
motifs_idx = np.argsort(mps, axis=1)[:, :2].
Since afterwards I have to compute the subspace to find out which k out of the d time series dimensions contains the desired motifs. Thus I am wondering if it suffices to compute the subspace as previously:
S = stumpy.subspace(dimensions, m, motif_idx, nearest_neighbor_idx, k)
or if I would need to increase the number of neighbors (say the second, third, ... nearest neighbor) in the function thus changing the subspace?

Secondly, I was wondering if there is a way to find the motifs with differing time length m, as seen in the screenshot using the multi dimensional matrix profile. I have read every paper of the original authors and they seem to offer different solutions for one dimensional cases only (instead of paper 6 of course, where they introduce mStamp). Is there any hope to find all my motifs using the multi dimensional approach or do I have to map them to a single dimension somehow and work with the algorithms that are available for one dimension? Such as taking the square root of x²+y²+z² respectively and try to merge the remaining dimensions with AB-joins somehow.

Thanks in advance!
Sarah

[seanlaw]

@SaVoAMP Thank you for your question and welcome to the STUMPY community!

So, you're absolutely right that the tooling for multi-dimensional matrix profiles is definitely lacking since it isn't possible/easy to choose k without human intervention. Additionally, multi-dimensional matrix profile stuff has been really HARD (especially for me!). This is still an active area of research and we're figuring things out along the way.

Perhaps, in your case, you are only interested in motifs for k=9 (i.e., including all time series dimensions) Given that you only have 9 dimensions in your time series, are you able to plot the multi-dimensional matrix profiles along with the raw time series (see the first figure in this section of the tutorial)?

To answer some of your questions:

I was wondering if it suffices to change the number of motifs I want to find in this line of code:
motifs_idx = np.argsort(mps, axis=1)[:, :2]

Yes, for motifs of length m, if you increase 2 to, say 10 then it should return the index locations for the top10 (local) motifs along each dimension. And then it would be up to you to decide on what k you want to use in order to retrieve the appropriate motif index location and then cross-reference this in your raw time series:

mps, indices = stumpy.mstump(df, m)
motifs_idx = np.argsort(mps, axis=1)[:, :10]

Thus I am wondering if it suffices to compute the subspace as previously:
S = stumpy.subspace(dimensions, m, motif_idx, nearest_neighbor_idx, k)

Yes, for each column in motifs_idx, you'll need to determine its corresponding k. So, continuing from the example above, you'd have 10 ks (one for each motif index). Then, for each independent motif index, you'll compute the corresponding subspace as you had shown. Note that it is completely possible to end up with different number of ks for the 10 different motif indices!

Now, say, for the 4th motif index (i.e., the 4th column of motifs_idx) we find that k = 1 (i.e., 2 dimensions participate in defining the motif) and, out of the set of time series, only T2 and T3 are important. Armed with this information, we can actually go hunting for this 4th motif in our original time series.

Note: From this point forward, none of this is supported and is subject to change without notice. In other words, this is purely experimental so please use at your own risk!

# Completely un-tested! Do not use this code!
import stumpy
import pandas as pd
import numpy as np
from stumpy.mstump import _query_mstump_profile
from stumpy import config, core

df = pd.read_csv("https://zenodo.org/record/4328047/files/toy.csv?download=1")
m = 30
mps, indices = stumpy.mstump(df, m)
motifs_idx = np.argsort(mps, axis=1)[:, :10]

# Assume that `k = 1` for the 4th motif index (`motifs_idx[:, 3]` - remember that it is zero based)
k = 1

T, M_T, Σ_T = core.preprocess(df[['T2', 'T3']], m)

excl_zone = int(
        np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM)
    )

# Not supported! Subject to change without notice! Use at your own risk!
P, I = _query_mstump_profile(
    motifs_idx[k, 3], T, T, m, excl_zone, M_T, Σ_T, M_T, Σ_T
)

print(P[k, :])

This should produce the corresponding multi-dimensional matrix profile, P, and multi-dimensional matrix profile indices, I, but this is only to look for the motif within df['T2', 'T3'] that is located starting at motifs_idx[k, 3]. Then, since P[k, :] is 1-dimensional, then you might be able to use stumpy.match to find the other patterns that match. Then, you can repeat this for the other motifs indices in motifs_idx.

This is the first time that I've thought through this and I'm sure that I've made some mistakes or forgotten something so please ask any clarifying questions!

Secondly, I was wondering if there is a way to find the motifs with differing time length m, as seen in the screenshot using the multi dimensional matrix profile.

Unfortunately, no. You'll have to re-run mstump with a different m and then go through the same process above.

I hope this is helpful!

[SaVoAMP]

Just saw that you edited your answer.

It's not entirely clear to me what you mean by motifs_idx[:, "Number of neighbor, starting at 2] or motifs_idx[k][2...last_neighbor]

I think that there might be a conceptual misunderstanding or a terminology mixup here. To make sure that we're on the same page, a motif is loosely defined as a subsequence that has another subsequence that (within some threshold) looks similar (i.e., a "close" nearest neighbor). So, as I had mentioned above, for preciseness, we should probably avoid using "nearest neighbor" when referring to the information that is captured in motifs_idx. Instead, I'd restate it as, "and I would iterate over all my desired motifs (say, 10 of them) and, for the ith motif, I would recalculate k_s[i] and S[i] by using motifs_idx[:, i] and similarly I would change the index location for the nearest neighbor in S[i] to stumpy.subspace(df, m, motifs_idx[k][i], indices[k][motifs_idx[k][i]], k)" (note the use of indices here).

Let me know if any of that may not be clear and we can discuss it further.

Yes, thanks a lot! By reasoning through your last answer I understood and adjusted the code accordingly 😄

I now also understand what this part of code is needed for:
P, I = _query_mstump_profile( motifs_idx[k, 3], T, T, m, excl_zone, M_T, Σ_T, M_T, Σ_T )
Unfortunately it doesnt work to call P[k, :]` since P is one dimensional.

Here is my full code for reference:

recordings = df.groupby(['recording_id'])
m = 100  # UNCLEAR HOW TO CHOOSE M (vor 1 second -> m=100)
reformat = 0  # Variable to reformat the time axis for each recording
desired_motifs = 10  # Number of desired motifs

# Iterate over each recording
for i, recording in enumerate(recordings):
    recording = recording[1]

    # # Visualize the data
    # visualize_data()

    # Get dimensions (x, y and z of acceleration, gyroscope and magnetic field --> 9 dims)
    dimensions = recording[['acc_x', 'acc_y', 'acc_z', 'gyr_x', 'gyr_y', 'gyr_z', 'mag_x', 'mag_y', 'mag_z']].copy()

    # Compute multi-dim. matrix profile
    # (Each row of the array corresponds to each MP for a given dimension --> first row: 1D MP, second row: 2D MP)
    mps, indices = stumpy.mstump(dimensions, m)

    # Get index of the k-dim. motif pair
    motifs_idx = np.argsort(mps, axis=1)[:, :desired_motifs]  # motifs must have the same mp value for the first indices

    for motif_number in range(desired_motifs):
  
        # Choose dimension k by automatically detecting turning point position
        k_s = KneeLocator(range(mps.shape[0]), mps[range(mps.shape[0]), motifs_idx[:, motif_number]], curve="convex",
                          direction="increasing")
        k = k_s.knee 
        #k_s.plot_knee()

        # Compute subspace to find out which k out of the d time series dimensions contains the desired motif pair
        motif_idx = motifs_idx[k, motif_number]
        nearest_neighbor_idx = indices[k, motif_idx] 
        S = stumpy.subspace(dimensions, m, motif_idx, nearest_neighbor_idx, k)
        print(f"For k = {k}, the {k + 1}-dimensional subspace includes subsequences from {dimensions.columns[S].values}")

        # query = dimensions[['mag_x', 'mag_y', 'mag_z']].copy()
        # T, M_T, Σ_T = core.preprocess(query, m)
        #
        # excl_zone = int(
        #     np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM)
        # )
        #
        # P, I = _query_mstump_profile(
        #     motif_idx, T, T, m, excl_zone, M_T, Σ_T, M_T, Σ_T
        # )
        #
        # print(P[k, :])

        # Plot the k-dim. motif pairs
        f, ax = plt.subplots(k + 1, 1, figsize=(20, 20), sharex=True)
        plt.suptitle(f"{motif_number + 1}. Motif")

        for j, s in enumerate(S):
            tmp = dimensions[dimensions.columns[s]]
            ax[j].plot(tmp)
            ax[j].set_title(f"{dimensions.columns.values[s]}")  # CHANGE NAME
            ax[j].plot(tmp[motif_idx: motif_idx + m], c='orange', linewidth=4)
            ax[j].plot(tmp[nearest_neighbor_idx: nearest_neighbor_idx + m], c='r',
                       linewidth=4)
            ax[j].axvline(x=motif_idx+reformat, linestyle='dashed', c='k', linewidth=2)
            ax[j].axvline(x=nearest_neighbor_idx+reformat, linestyle='dashed', c='k', linewidth=2)
        plt.show()
    reformat = reformat + dimensions.shape[0]
    plt.show()

As you can see I iterate over each recording (belongs to different data frames), calculate mps and indices with mstump and get motifs_idx for the first 10 columns (belonging to desired_motifs). Afterwards I iterate through each number in desired_motifs and calculate k and then S with S = stumpy.subspace(dimensions, m, motif_idx, nearest_neighbor_idx, k) where motif_idx = motifs_idx[k, motif_number] and nearest_neighbor_idx = indices[k, motif_idx].

I commented out the middle part because of the reasons mentioned in my previous reply. 'Reformat' just adjusts the axis of my data.

Thanks for all of your suggestions!

[seanlaw]

Okay, I think the next missing piece is that we can't simply do motifs_idx[:, desired_motifs] since, as you had correctly, pointed out, would give you motifs that are very close in sequence/time/index value to each other. So, what you'd need to do is to "disqualify" subsequences that are "near" a motif that has already been identified by:

1. Identify the first motif, curr_motif_idx = np.argsort(mps, axis=1)[:, 0]
2. Use this motif as your query to compute a multi-dimensional distance profile (not matrix profile) - see the _get_multi_distance_profile function above
3. Identify the "best matches" using the distance profile
4. Set the exclusion zone surrounding curr_motif_idx to np.inf for mps via core.apply_exclusion_zone(mps, curr_motif_idx, excl_zone)
5. Now repeat step 1 but note that since we set mps values to np.inf, that means that the next motif will not be discovered anywhere surrounding the first motif index. And, as you search for more and more motifs, more and more of mps gets "disqualified" since exclusion zones surrounding previous motifs will have been set to np.inf

We are essentially doing all of this in stumpy.motifs but it currently only works for 1-D time series since multi-dimensional time series require much more care as you can see.

[SaVoAMP]

It´s working, thank you very much! I think we are getting closer 😃

I exchanged this part P, I = _query_mstump_profile( motif_idx, T, T, m, excl_zone, M_T, Σ_T, M_T, Σ_T ) with the new code you suggested and plotted the first 10 nearest neighbors in addition to the already existing motif pair (some appear to be pretty close to each other).

[seanlaw]

Great! It's not going to be perfect but I think you'll be mostly there once you start "disqualifying" regions in mps that are near a previously identified motif.

Note that with the distance profile, you'll also have to apply a "disqualification" trick where you set column values in D that are similar to your motif (i.e., a "match") to np.inf before you search for the next similar "match". I strongly recommend looking through this motifs.py file to see how 1) a motif is identified and then parts of the matrix profile is disqualified and 2) how to find matches in the distance profile, D, while also disqualifying locations that have already been identified as a "match".

After this discussion is completed, I will try to make some of this a little easier for the next person. You are delving into uncharted territory!

[SaVoAMP]

After this discussion is completed, I will try to make some of this a little easier for the next person. You are delving into uncharted territory!

Sounds good! I will answer as soon as I have tried your latest tips (hopefully tomorrow)!

[SaVoAMP]

Hey,

Now I am able to find the 10 nearest neighbors for each motif! Thank you very much!
Here is my code:

S = stumpy.subspace(dimensions, m, motif_idx, nearest_neighbor_idx, k)
print(f"For k = {k}, the {k + 1}-dimensional subspace includes subsequences from {dimensions.columns[S].values}")

subspace = [dimensions.columns.values[s] for s in S]

query = dimensions[subspace].copy()
T, M_T, Σ_T = core.preprocess(query, m)

excl_zone = int(
        np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM)
    )

# Compute multi-dimensional distance profile
D = _get_multi_distance_profile(
    motif_idx, T, T, m, excl_zone, M_T, Σ_T, M_T, Σ_T
)

matches = []  # stores the distance and idx of the max_matches nearest neighbors
max_matches = 10  # 10 nearest neighbors
candidate_idx = np.argmin(D[k, :])
while len(matches) < max_matches:
    matches.append([D[k, :][candidate_idx], candidate_idx])
    core.apply_exclusion_zone(D[k, :], candidate_idx, excl_zone)
    candidate_idx = np.argmin(D[k, :])

However I still have problems to disqualify subsequences that are "near" a motif that has already been identified since core.apply_exclusion_zone(mps, candidate_idx, excl_zone) is not working.
Here the corresponding code:

motif_indices = []
motif_distances = []
candidate_idx = np.argsort(mps, axis=1)[:, 0]  # Identify the first motif as candidate motif
while len(motif_indices) < desired_motifs:
    profile_value = mps[:, candidate_idx]
    motif_indices.append(candidate_idx)
    motif_distances.append(profile_value)
    core.apply_exclusion_zone(mps, candidate_idx, excl_zone)
    candidate_idx = np.argsort(mps, axis=1)[:, 0]

What am I doing wrong? (probably it's just a really stupid mistake)

[seanlaw]

Since k can be selected automatically now, shouldn't we generalize the functions stumpy.motifs and stumpy.matches to the multidimensional case for our motif discovery tutorial?

Let's discuss in #187

[seanlaw]

In case it matters, I just updated the tutorial and added a new section on MDL:

https://stumpy.readthedocs.io/en/latest/Tutorial_Multidimensional_Motif_Discovery.html#Finding-k---A-Minimum-Description-Length-Approach

Note that you'll need to clone the latest development version of STUMPY to use it. Any and all feedback is welcome before we cut a new release version!

[SaVoAMP]

I haven't tested it yet, but I like the tutorial part!
I want to read up on MDL a bit in the days, then I can certainly give better feedback! At the latest in the tutorial I will work with MDL anyway.

The only comments I could give now are rather small typos:

Unfortunately, the naive elbow-finding method above isn’t robust and cannot be easily automated so the authors of the Matrix Profile VI paper recommends computing somethng called the “minimum description length” (MDL)

In the first sentence you forgot the i in something and I would write recommend instead of recommends since I think this refers to the authors.
Otherwise, I find it only slightly confusing that the subspace is partly capitalized (S) and partly lowercase (s).

So in the code blocks there is a capital S, but in the explanations there is a small (pink) s.
But these are all just total trifles. I think the tutorial part is really good!

[seanlaw]

In the first sentence you forgot the i in something and I would write recommend instead of recommends since I think this refers to the authors.

Perfect! I've corrected that now

So in the code blocks there is a capital S, but in the explanations there is a small (pink) s.

That is odd and I can see where that could lead to confusion. In my notebook, the pink S is capitalized in my browser. When I copy and paste "and S is" from that sentence into another text document, I see a clear difference in the height. Of course, there is no doubt in what what you pasted above in that the pink S certainly appears shorter than it should be. Nonetheless, I've incorporated your feedback and swapped S for subspace to hopefully avoid any problems for future readers.

You may also notice that I simplified a lot of the motif_idx code along the way so it only grabs the top-1 motif and then I use that to retrieve the nearest neighbor from indices. This should be done consistently throughout the notebook.

But these are all just total trifles. I think the tutorial part is really good!

No trifling at all! Your eye for detail is 100% appreciated and other users will too.

[SaVoAMP]

Perfect!
There is only one s left above code block [14].