Finding the top-k motif (similar subsequence) of a given subsequence in T

[cyuab]

Given a time series T, I am interested in finding the top-K similar subsequence of a given subsequence T' in T.
In addition, I would like to find those before T'.

Suppose |T| =50 and m =5,
for the subsequence T'= T[30, 30+5], I would like to find the subsequences that are similar to T' and before.
I am thinking of using either the left matrix profile or stumpi by ignoring those points after T'.
But I don't know how to obtain the top-2 and top-3 and etc of T'
When I use stump and set k>1, it seems the returning top-k motifs are overlapping. I think it does not consider the exclusion zone requirement.

[seanlaw]

@cyuab Thank you for your question and welcome to the STUMPY community.

When I use stump and set k>1, it seems the returning top-k motifs are overlapping.

Yes, that sounds correct. Given that the matrix profile is a brute force pairwise distance calculation, this means that ALL distances between T' = T[30 : 30 + 5] and T are computed with the exception that the distances are NOT computed in any particular order (assume that the pairwise distances are computed at random) and only the top k are retained. That is, the first k (random) distances are computed and retained and, when a new distance is computed and found to be smaller than one of the existing distances, it replaces that larger distance. Because of this, at run time, it is impossible to know:

1. whether a future, newly identified subsequence will overlap with any existing (previously identified) subsequence that is already in your top-k
2. when there is overlap (between the matches), which subsequence to prefer to keep in the top-k as this would require making assumptions for the user

Thus, we chose speed over committing to making unnecessary overlap assumptions for the user. And, instead, the user can most likely increase k (i.e., double the value of k) and manually eliminate the subsequence according to their desired assumptions of what "overlap" means to them.

I think it does not consider the exclusion zone requirement.

Typically (there are exceptions, see below!), the "exclusion zone" only refers to the zone that is close to the original T = T[30 : 30 + 5]subsequence and does not have anything to do with other possible subsequence matches. All of the subsequence matches that fall outside of the exclusion zone ofT'` should be retained.

I would like to find the subsequences that are similar to T' and before. I am thinking of using either the left matrix profile

Yes, I believe that the left matrix profile is what you want. Note that you will need to use the left-index to retrieve the correct subsequence and then to recompute the distance (don't forget to z-normalize where necessary!)

But I don't know how to obtain the top-2 and top-3 and etc of T'

Normally, if you already know what T' is, then you can usually get pretty far by using the built in stumpy.match function (see this tutorial). Since you are looking for subsequences to the left of T' then all you need to do is account for this and simply modify T accordingly:

# Untested example but hopefully you can get the point
import numpy as np
import stumpy
from stumpy import config, core
np.random.seed(0) # Set random seed for reproducibility
T = np.random.rand(1000)
m = 5
start_idx = 300
stop_idx = start_idx + m
T_prime = T[start_idx : stop_idx] # T[300 : 300 + 5]
T_copy = T.copy() # Make a copy of T
# Apply exclusion zone AND clear all values to the right of the start index
excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
core.apply_exclusion_zone(T_copy, start_idx, excl_zone, np.nan)
T_copy[start_idx:] = np.nan
matches = stumpy.match(T_prime, T_copy, max_matches=3) # Note that this makes an internal assumption and applies the exclusion zone that you are seeking
# array([[np.float64(0.4800008207705901), np.int64(44)],
# [np.float64(0.6182357545060001), np.int64(228)],
# [np.float64(0.8576581684735803), np.int64(275)]], dtype=object)

Please feel free to ask any clarifying questions and I'll try my best to help.

[cyuab]

I want to comment on the above code snippet
The above snippet is not consistent with what we have discussed.
For example:

np.random.seed(0)
T = np.random.rand(20)
m = 5
mp = stumpy.stump(T, m)
mp.left_I_
# mparray([-1, -1, -1, 0, 1, 2, 3, 2, 1, 4, 2, 0, 7, 1, 4, 3])
mp.left_I_[3]
# 0

• Since m=5, given excl_zone = int(np.ceil(m / 4)) (from here), excl_zone = 2.
• The left nearest neighbor of T[3:8] is T[0: 5]. (p.s., left inclusive and right exclusive indexing)
  • 0 is not within the 2-neighbors of 3
  • since the starting indices of T[3:8] (i.e., 3) and T[2: 7] (i.e., 2) overlap; (i.e., 2 is within the 2-neighbors of 3)
  • since the starting indices of T[3: 8] (i.e., 3) and T[1: 6] (i.e., 1) overlap. (i.e., 1 is within the 2-neighbors of 3)
• So mp.left_I_[3] = 0

The above code (set start_idx = 3) shows:

start_idx = 3
stop_idx = start_idx + m
T_prime = T[start_idx : stop_idx] T_copy = T.copy() # Make a copy of T
# Apply exclusion zone AND clear all values to the right of the start index
excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
core.apply_exclusion_zone(T_copy, start_idx, excl_zone, np.nan)
T_copy[start_idx:] = np.nan
matches = stumpy.match(T_prime, T_copy, max_matches=3)
matches
# array([], dtype=object)

This second result is not consistent with the first result.

The correct version:

start_idx = 3
stop_idx = start_idx + m
T_prime = T[start_idx : stop_idx] T_copy = T.copy() # Make a copy of T
excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
# core.apply_exclusion_zone(T_copy, start_idx, excl_zone, np.nan)
T_copy[start_idx + excl_zone:] = np.nan
matches = stumpy.match(T_prime, T_copy, max_matches=1) # Note that this makes an internal assumption and applies the exclusion zone that you are seeking
print(matches) # [[np.float64(4.197447479312608) np.int64(0)]]

[seanlaw]

Since m=5, given excl_zone = int(np.ceil(m / 4)) (from here), excl_zone = 2.
The left nearest neighbor of T[3:8] is T[0: 5]. (p.s., left inclusive and right exclusive indexing)
0 is not within the 2-neighbors of 3
since the starting indices of T[3:8] (i.e., 3) and T[2: 7] (i.e., 2) overlap; (i.e., 2 is within the 2-neighbors of 3)
since the starting indices of T[3: 8] (i.e., 3) and T[1: 6] (i.e., 1) overlap. (i.e., 1 is within the 2-neighbors of 3)
So mp.left_I_[3] = 0

You're absolutely right! Thank you for the comment/correction and apologies for the confusion. Truth-be-told, internally, we only ever apply core.apply_exclusion_zone to the results of a matrix profile or distance profile and it is never applied to a raw time series input. I blindly applied it to an input time series and I should've been more careful/thoughtful and I hope it didn't cause you too much frustration! Sorry for my laziness. 😓

T_copy[start_idx + excl_zone:] = np.nan

Ashamedly, I believe that the more precise "masking range" should be T[start_idx - excl_zone - 1 + m:] = np.nan

import numpy as np
from stumpy import config, core
import stumpy
np.random.seed(0)
T = np.random.rand(20)
m = 5
start_idx = 3
stop_idx = start_idx + m
T_prime = T[start_idx : stop_idx]
T_copy = T.copy() # Make a copy of T
# Apply exclusion zone AND clear all values to the right of the start index
excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
#core.apply_exclusion_zone(T_copy, start_idx, excl_zone, np.nan)
T_copy[start_idx - excl_zone - 1 + m:] = np.nan
matches = stumpy.match(T_prime, T_copy)
print(matches)

To break down start_idx - excl_zone - 1 + m a bit more so that we can better understand it:

1. First, let's rearrange some of the terms and re-write this as start_idx - (excl_zone + 1) + m
2. Next, the start_idx anchors our reference index location
3. Then, relative to the start_idx, we subtract (excl_zone + 1) to identify the index location of the LAST VALID subsequence that is "eligible" to be considered as a "match" (i.e., all subsequences between index 0 and start_idx - (excl_zone + 1) would be an allowed subsequence index to find a match for)
4. Finally, now that we've identified the LAST VALID subsequence index location, we add m to ensure that (i) the entire length of this (last valid) subsequence is "allowed" for a distance calculation and (ii) it leaves us at the start idx of the first INVALID subsequence index
5. Lastly, this only applies for left-matrix-profiles (i.e., not full matrix profiles or right matrix profiles!!!)

Please let me know if you disagree with anything

[cyuab]

Yes, I think your masking range is correct. Thanks for the clarification.

[cyuab]

@seanlaw Thank you so much for your time and response.
I have two things that I want to discuss further.

(1)

Typically (there are exceptions, see below!), the "exclusion zone" only refers to the zone that is close to the original T = T[30 : 30 + 5]subsequence and does not have anything to do with other possible subsequence matches.

It seems more natural for me that the "exclusion zone" applies to all the subsequences that are identified as a member of Motif.
According to (p. 15 of slides of https://www.cs.ucr.edu/%7Eeamonn/TimeSeriesMotifs/):

In the above figure, there is a GUI panel.
In the first row, there are two groups of subsequences.
The first group is "Saw motif" (2 subsequences).
The other group is "Sine morif" (3 subsequences).
The panel shows that there is no meaningful motif for the Top-3 motif.
I think they have enforced an "exclusion zone" around all the selected motifs.
Otherwise, the Top-3 motif should consist of (the non-selected "Sine subsequence") and (the selected "slightly-shift" version of the selected "Sine subsequence").

So, I am wondering if we should impose the "exclusion zone" on the k-nearest neighbors of T' so that these neighbors are not just the "slightly-shift" version of the other one?

(2)
For a given time series $T$, I would like to find the top-K nearest neighbor of each subsequence $T'$ with length = $m$ in $T$AND the neighbors are before $T'$.
i.e., the subsequences will be T[0, 0+m], T[1, 1+m], ... T[-m, len(T)]
I am thinking about how to implement it. The following is my sketch idea

Using stumpi.
In each iteration, we treat the last m elements as T' and find the top-K nearest neighbors of it. (Using stumpy.match)
Add a new time point t to T to grow it
Repeat

Do you have any suggestions/ comments about it.
Thanks in advance.

[seanlaw]

So, I am wondering if we should impose the "exclusion zone" on the k-nearest neighbors of T' so that these neighbors are not just the "slightly-shift" version of the other one?

If this is what you want, you can achieve this by simply increasing the number of k nearest-neighbors to return and then systematically pruning the ones that are a slight shift from each other. However, this is not an assumption that we can make for the user as there are situations when a user prefers to identify the slightly shifted matches and it is easier to post-process rather than to throw the information away.

Note that both the stumpy.motifs and stumpy.match functions apply the exclusion zone that you are looking for. However, they are performed on a full matrix profile rather than on a left-matrix-profile. You can examine how these functions were developed and come up with your own requirements. They were created not to do handle all cases but to serve as an example to build from.

Do you have any suggestions/ comments about it.

So, stumpi actually gives you _left_P_ and left_I_ attributes for every existing and incoming subsequence (but only for the top-1) and you should probably choose some sort of distance threshold such that if any of the _left_P_ values are below a given threshold then you can manually compute the (left) distance profile for that new, incoming subsequence and its corresponding matches. Without understanding your actual use case, I would typically argue that one rarely ever want to know the top-k matches for EVERY subsequence because not every subsequence is interesting (and most have a 1-nearest-neighbor that is far away anyway) and that it is best to save your costly computations only to subsequences that have a minimum of one "close" nearest neighbor (again, knowing that ALL other subsequences MUST BE farther away than the 1-nearest neighbor).

[cyuab]

Hi,

I have studied your suggested tutorial (see this tutorial), stumpy.match definitely suits my needs by applying an exclusion zone to each of the newly discovered neighbors.

I have a follow-up idea/question about this.
According to (p. 61 of slides of https://www.cs.ucr.edu/%7Eeamonn/TimeSeriesMotifs/):

We can "enrich" (include more members) the Top-1 Motif pair by adding other members that are similar to the center (defined by $A$ and $B$) with respect to a threshold (= 2 * ED($A$, $B$)) to form a "motif family"
I am thinking we should also add these members ($C$) in our exclusion zone, so that when we calculate the Top-2 motif and its family members, we can find unique occurrences.

Is it possible to achieve this in the current implementation of STUMPY?
Or should I add their exclusion zone manually as your code suggested?

excl_zone = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
core.apply_exclusion_zone(T_copy, start_idx, excl_zone, np.nan)

Where start_idx will be the newly discovered members

[seanlaw]

Or should I add their exclusion zone manually as your code suggested?
Where start_idx will be the newly discovered members

Yes, I would apply the exclusion zone one at a time as you add more matches to the "motif family". Having said that, you're starting to customize your analysis and so i recommend looking at the stumpy.match code to learn more about how things are currently being handled and, hopefully, you can see that there's nothing special going on behind the scenes:

stumpy/stumpy/motifs.py

Lines 532 to 585 in 534488d

	Q=core._preprocess(Q)
	T=core._preprocess(T)
	ifnp.any(np.isnan(Q)) ornp.any(np.isinf(Q)): # pragma: no cover
	raiseValueError("Q contains illegal values (NaN or inf)")
	m=Q.shape[-1]
	excl_zone=int(np.ceil(m/config.STUMPY_EXCL_ZONE_DENOM))
	Q_subseq_isconstant=core.process_isconstant(Q, m, Q_subseq_isconstant)
	T_subseq_isconstant=core.process_isconstant(T, m, T_subseq_isconstant)
	ifQ.ndim==1:
	Q=np.expand_dims(Q, 0)
	ifQ_subseq_isconstant.ndim==1:
	Q_subseq_isconstant=np.expand_dims(Q_subseq_isconstant, 0)
	ifT.ndim==1:
	T=np.expand_dims(T, 0)
	ifT_subseq_isconstant.ndim==1:
	T_subseq_isconstant=np.expand_dims(T_subseq_isconstant, 0)
	T[np.isinf(T)] =np.nan
	ifM_TisNoneorΣ_TisNone:
	M_T, Σ_T=core.compute_mean_std(T, m)
	T[np.isnan(T)] =0
	iflen(M_T.shape) ==1:
	M_T=np.expand_dims(M_T, 0)
	iflen(Σ_T.shape) ==1:
	Σ_T=np.expand_dims(Σ_T, 0)
	d, n=T.shape
	D=np.empty((d, n-m+1))
	foriinrange(d):
	D[i, :] =core.mass(
	Q[i],
	T[i],
	M_T[i],
	Σ_T[i],
	T_subseq_isconstant=T_subseq_isconstant[i],
	Q_subseq_isconstant=Q_subseq_isconstant[i],
	query_idx=query_idx,
	)
	D=np.mean(D, axis=0)
	returncore._find_matches(
	D,
	excl_zone,
	max_distance=max_distance,
	max_matches=max_matches,
	query_idx=query_idx,
	atol=atol,
	)

To reiterate, the stumpy.match function wasn't meant to serve all purposes but should work for the majority of the cases and, hopefully, provide inspiration for you to come up with your own solution. As you've already identified, once you understand the basic matrix profile data structure, there's a lot that you can do with it and it doesn't usually take too much extra work to analyze the results to your liking.

[cyuab]

Thank you so much for your help! Matrix Profile and STUMPY are interesting to play with and do useful stuff!