Define a sliding inclusion zone to exclude matches "too far in time"

[javierdvalle]

I am new to Stumpy. I have read in the tutorial about the exclusion zone relative to the diagonal to exclude trivial matches. I am not sure if this option exists because I've not seen it in the API docs.

But anyway, I would like to know if there exists the opposite option: defining a sliding inclusion zone to only include matches inside the inclusion zone, and exclude matches outside (too far in time).

For example, with a daily time serie I would like to exclude matches beyond two months.

[seanlaw]

I am new to Stumpy. I have read in the tutorial about the exclusion zone relative to the diagonal to exclude trivial matches. I am not sure if this option exists because I've not seen it in the API docs.

But anyway, I would like to know if there exists the opposite option: defining a sliding inclusion zone to only include matches inside the inclusion zone, and exclude matches outside (too far in time).

For example, with a daily time serie I would like to exclude matches beyond two months.

@javierdvalle Welcome to the STUMPY community and thank you for your question. I can see how this feature could be very useful in certain use cases. Unfortunately, to be consistent with the originally published work, we have chosen a default exclusion zone and there is currently no way to do what you are asking.

However, if you can provide a little more detail about your use case, the time series size/length, and the frequency of your data then I ma be able to offer an alternative solution.

[javierdvalle]

Hi @seanlaw thank you for your response and congratulations for the library!

My use case is the following: We have a one-year mobility dataset with the daily number of people traveling from point A to point B. Actually it comes as a daily origin-destination matrix but to keep it simple we want to start with just one unidimensional time serie.

One time serie looks like this:

There are several patterns there:

• First, the lockdown effect during March-June is super clear and was easy to detect using FLUSS Semantic Segmentation.
• Second, for the typical week the peak number of trips is on Friday, but there are some weeks with bridging days where the peak appears on Thursday, or Wednesday, and we would like to detect those (btw, is there a way to force stumpy to compare full weeks i.e. monday with monday, tuesday with tuesday and so on?).
• Last, the peak with the red dot at 2020-12-03 is not only anomalous because it appears on Thursday, but also because the number of trips was almost double than all the peaks in a two months window. This one is what motivated my question, that is why I would like each week to be compared just with close-in-time weeks.

[seanlaw]

@javierdvalle Thank you for the context.

First, the lockdown effect during March-June is super clear and was easy to detect using FLUSS Semantic Segmentation.

Wonderful!

Second, for the typical week the peak number of trips is on Friday, but there are some weeks with bridging days where the peak appears on Thursday, or Wednesday, and we would like to detect those (btw, is there a way to force stumpy to compare full weeks i.e. monday with monday, tuesday with tuesday and so on?).

Technically, no. I think that what matrix profiles are really great at is when you don't know if there are any conserved behaviors in your data and to find the nearest neighbor for each subsequence in an unconstrained way. Thus, a matrix profile will compare monday with monday, monday with tuesday, monday with wednesday, and so on and return the absolute best match. This is an exhaustive search. In your case, it sounds like you have an idea of where things may be conserved/repeating. There are a few non-traditional ideas that come to mind (and please forgive me as they are not well thought out or are poorly formulated in my head):

Let's say you want to compare full weeks (7 days, starting on a monday and ending on a sunday) and align the days of the week. What you can do is insert np.nan between every sunday/monday and then feed this into STUMPY. These np.nan basically tells STUMPY to return a np.inf distance value if the 7 day sliding subsequence in question contains a np.nan. In other words, this makes it so that only subsequences that start on monday and end on sunday are "valid" and ensures that only monday-sunday subsequences are compared with monday-sunday subsequences. The only challenge here is that you need to account for the indices at these np.nan values so your matrix profile indices will be "off" from your original (unmodified, i.e., no np.nan) data.

Let me know if that makes sense?

Last, the peak with the red dot at 2020-12-03 is not only anomalous because it appears on Thursday, but also because the number of trips was almost double than all the peaks in a two months window. This one is what motivated my question, that is why I would like each week to be compared just with close-in-time weeks.

I think it's important to point out that matrix profiles are really good at identifying repeating shapes and not great at finding single point peaks or anomalies (due to z-normalization of the subsequences). I find that a 7 day window is really pushing it and might even prefer a 14 day window as 7 data points isn't that much data to define a "shape".

On the flip side, since you have daily data (which is "small"), you might just be able to run stumpy.stump multiple times on shorter segments of your data. This is not usually possible if your data was of higher frequency so luck is on your side here. Here's what I propose:

Let's say you have 365 days worth of data and a month is considered 30 days. In your original post, you stated that you wanted to exclude any matches beyond, say, two months (i.e., 60 days), then you can do something like:

import stumpy

T = np.random.rand(365)  # This will represent 365 days of your data. You may want to insert np.nan as we described above
m = 7  # 
for i in range(T.shape[0] - m + 1):
    Q = T[i : i + m]
    mp = stumpy.stump(Q, m, T[i + m : i + m + 60])

This isn't super pretty (and admittedly hacky) but you are basically extracting out one week query chunks, Q, and computing the matrix profile for Q and a 60 day window after Q, T[i + m : i + m + 60]. The matrix profile that is returned will then correspond only to subsequences that have occurred within 60 days AFTER Q. Of course, as a reminder, the matrix profile indices will be "off" and not representative of your full length T so you'll have to perform some bookkeeping on your own.

I hope that gets you going in the right direction but let me know if I can help clarify anything!

[seanlaw]

For the last part, just to check I got it: from the mp inside the for loop I can get the maximum on each iteration and that would be my "matrix profile with an inclusion zone", right?

Yes, that is correct. Although, getting the "maximum on each iteration" would grab the local discord (potential anomaly) and not the local motif (most conserved pattern). Essentially, you are extracting out Q = T[i : i + m] and then comparing it to T[i + m : i + m + 60] and computing the matrix profile. Note that since this is an "AB-join" (and not a self-join), there is no concept of an "exclusion zone" (i.e., EVERYTHING is considered "fair game" for nearest neighbor considerations).

In case others stumble across this question, I want to take the opportunity to re-iterate that the matrix profile indices (i.e., mp[:, 1]) will be "off" as they will be zero-indexed based on the length of T[i + m : i + m + 60] (i.e., the matrix profile indices are not relative to the full length of T). So, one will need to take extra care here.

[jgill018]

Thank you so much for asking this question and for the detailed answer. I have a similarly motivated need for an inclusion zone, but a different enough use case that I don't think inserting NANs into the sequence will solve my problem. I am working with a large 1D sequence (approx. 1e9 pts) and need a quick result (<2 seconds on a standard laptop). Again, I am more interested in comparing a subsequence to it's more immediate neighbors (using both an exclusion zone to prevent trivial matches and an inclusion zone to only search a region around the sliding window). I can write a for-loop that generates the sliding window query and search data (AB-join), where the search data is the region around the sliding window (making sure to ignore the trivial self-matches). The matrix profile from inside that for loop would be further reduced, by taking the min distance and populating the final matrix profile output (sorry this is confusing...). Not sure if there is a better way - like the slick NaN solution you described above? Or perhaps the streaming implementation could be used somehow?

[seanlaw]

@jgill018 Welcome to the STUMPY community!

Unfortunately, with 1e9 data points (which is HUGE), you definitely don't want to use standard methods as it will either take many days and/or a lot of resources to compute the full matrix profile. Considering that you are only interested in small regions near each query subsequence, I think it is best to try what you had proposed as it would serve to prune away unnecessary calculations. Though, given the data size, I'm not sure that you'll be able to get it down to <2 seconds and you'll likely need some distributed computing capability.

Your idea of using stumpi is interesting but the matrix profile is only updated if the newly ingressed data point produces a lower distance. I think your best bet is to go the for-loop route but maybe you can share what you were considering and we can brainstorm together?

[seanlaw]

Also, what you mention about the z-normalization is important, because in our use case the absolute number of trips matters a lot, I tried using AAMP and I think got better results.

@javierdvalle Just a heads up that starting in STUMPY v1.8.0, you no longer have to call aamp and instead, you can/should use the normalize=False parameter:

stumpy.stump(T, m, normalize=False)

This is equivalent to:

stumpy.aamp(T, m)

In the background, when normalize=False, stumpy.stump just calls stumpy.aamp for you and hands over the relevant data.

[seanlaw]

@jgill018 I put some more thought into your question and I believe that this should work:

#!/usr/bin/env python

import numpy as np
from stumpy.stump import _stump
from stumpy import core

def jgill_func(T_A, m, incl_zone=1):
    """
    """
    T_B = T_A
    ignore_trivial = True

    (
        T_A,
        μ_Q,
        σ_Q_inverse,
        μ_Q_m_1,
        T_A_subseq_isfinite,
        T_A_subseq_isconstant,
    ) = core.preprocess_diagonal(T_A, m)

    (
        T_B,
        M_T,
        Σ_T_inverse,
        M_T_m_1,
        T_B_subseq_isfinite,
        T_B_subseq_isconstant,
    ) = core.preprocess_diagonal(T_B, m)

    core.check_window_size(m, max_size=min(T_A.shape[0], T_B.shape[0]))

    n_A = T_A.shape[0]
    n_B = T_B.shape[0]
    l = n_A - m + 1

    excl_zone = int(np.ceil(m / 4))
    out = np.empty((l, 4), dtype=object)

    diags = np.arange(excl_zone + 1, max(excl_zone + incl_zone, n_A - m + 1))

    P, I = _stump(
        T_A,
        T_B,
        m,
        M_T,
        μ_Q,
        Σ_T_inverse,
        σ_Q_inverse,
        M_T_m_1,
        μ_Q_m_1,
        T_A_subseq_isfinite,
        T_B_subseq_isfinite,
        T_A_subseq_isconstant,
        T_B_subseq_isconstant,
        diags,
        ignore_trivial,
    )

    out[:, 0] = P[:, 0]
    out[:, 1:] = I

    return out


if __name__ == "__main__":
    T = np.random.rand(10000)
    m = 50

    mp = jgill_func(T, m, incl_zone=5)

Here, the incl_zone is "how many sliding window subsequences adjacent to the exclusion zone" do you want to compute. So, incl_zone=1 would only compute a single subsequence adjacent to the exclusion zone while incl_zone=5 would compute the five windows adjacent to the exclusion zone.

Unfortunately, with your large (1e9) data size, it is impossible to get a result in < 2 seconds 😢