Improve Snippet function

[NimaSarajpoor]

As I was working with snippets, I noticed the three following issues:

1- Apparently the indices of the snippet start exactly at the multiplies of m, which may not be the case for some data.

2- The function provides only one profile per snippet. However, that snippet repeats throughout the whole time series as higher fraction means it has been repeated more. So, it should be useful to have that information in the output, where I can see all the profiles (the starting index) of each snippet.

3- The fraction output is not sorted. I think it is better to sort it in descending order and accordingly provide the other outputs for the same order.

[seanlaw]

@ninimama If it's okay with you, I'd like to move this over to the discussion section as I think it would be healthy/useful to discuss/debate each of the points above. We should keep in mind that STUMPY doesn't try to be everything for everyone and our goal is to reproduce the published work.

[NimaSarajpoor]

@seanlaw

Sure. I think it is good to move 1 & 3 to discussion. However, if you would like to reproduce the results of the paper, I think the snippet function should produce the indices of a set of subsequences similar to each snippet (e.g: check out Fig 24 of the paper or Fig 19).

However, if you think the ultimate goal of snippets is to merely provide a set of representative sequences rather than how they appear in the time series, then we can move 2 to discussion as well.

[seanlaw]

@ninimama I haven't looked into this part much so maybe you can help me understand how the indices are calculated? It appears that it has do with finding cross-over points.

However, if you think the ultimate goal of snippets is to merely provide a set of representative sequences rather than how they appear in the time series, then we can move 2 to discussion as well.

I think this is what I'm trying to get at with an open discussion. I am open to being convinced :)

[NimaSarajpoor]

This is part of the snippet code. If I understand correctly, mask contains the indices of distance profile corresponds to the snippet (in boolean). And, by summing over that, we calculate the fraction.

So, if my understanding is correct, all we need to do is to return those indices.

Of course there might be some overlaps between one subsequence to the next, but we don't need to worry about that. We just need to return those indices. So, coloring them should give us the Fig. 24 of the paper.

Please feel free to close this one and move the questions to the discussion.

[seanlaw]

1- Apparently the indices of the snippet start exactly at the multiplies of m, which may not be the case for some data.

I did notice this when I initially implemented it according to the paper and vaguely remember thinking "for some reason, this is good enough". I will need some time to think about this. Maybe you can describe what you had in mind?

2- The function provides only one profile per snippet. However, that snippet repeats throughout the whole time series as higher fraction means it has been repeated more. So, it should be useful to have that information in the output, where I can see all the profiles (the starting index) of each snippet.

Thank you for providing more context/information. I came across a solution in the past where you can pass in a mask and have it return the start (inclusive) and stop (exclusive) indices where the mask is True

import numpy as np

def _get_mask_slices(mask):
    m1 = np.r_[0, mask]
    m2 = np.r_[mask, 0]
    starts, = np.where(~m1 & m2)
    ends, = np.where(m1 & ~m2)
    return np.c_[starts, ends]

if __name__ == '__main__':
    mask = np.array([ True,  True, False,  True,  True,  True,  True, False, False, False], dtype=bool)
    slices = _get_mask_slices(mask)
    print(slices)
    # array([[0, 2],
    #        [3, 7]])

This is a two column start and (exclusive) stop indices and, perhaps, you could append a third column with k? Maybe we should call these snippets_regimes in the code. Would you mind confirming if this would work?

Additionally, any changes to snippets.py would also need to happen in aampdist_snippets.py (this is our non-normalized version - i.e., without z-normalization).

3- The fraction output is not sorted. I think it is better to sort it in descending order and accordingly provide the other outputs for the same order.

I don't know if this (below) is correct but were you thinking of something like:

sorted_idx = np.flip(np.argsort(snippets_fractions, kind='mergesort')) 
snippets = snippets[sorted_idx]
snippets_indices = snippets_indices[sorted_idx]
snippets_profiles = snippets_profiles[sorted_idx]
snippets_fractions = snippets_fractions[sorted_idx]
# Do not re-order snippets_areas!!

Again, I have not verified exactly which arrays need to also be sorted so it would be great if you could help me verify. Regardless, I think this change would be acceptable

[seanlaw]

In case it matters, I wanted to share some of my thought process. Usually, I'm looking to see if proposed changes would break backwards compatibility and affect other users who are using the code. Generally speaking, in 2. above, I think the proposed change would break a user's code because an additional array (snippets_regimes) would need to be returned and so the API would be different. However, since snippets is brand new and not yet part of an official release then we should be fine (i.e., there is little to no impact on users)!

[NimaSarajpoor]

1- Apparently the indices of the snippet start exactly at the multiplies of m, which may not be the case for some data.

I did notice this when I initially implemented it according to the paper and vaguely remember thinking "for some reason, this is good enough". I will need some time to think about this. Maybe you can describe what you had in mind?

So, let's say I have a time series, T, of length n, and I am looking for a snippet of size m. And, let's assume everything is good, and we have the following structure in our time series: snippet1 - snippet2 - snippet1 - ....

What will happen if, for some reason, I try to find the snippets on T[(m/2):]?

This is what I got when trying to apply snippet function on T[100:] on the tutorial. Could you please check it out yourself and let me know if that's the same output you get? Please let me know if my example makes sense.

This is a two column start and (exclusive) stop indices and, perhaps, you could append a third column with k? Maybe we should call these snippets_regimes in the code. Would you mind confirming if this would work?

Yes, that would be a nice output. As you said, we can make it a 3-dim array to include the snippets_regimes for all snippets.

If I may ask, what did you mean when you said "Would you mind confirming if this would work?" ? Should I apply it to a time series to see if the output makes sense?

By the way, I mentioned the same problem to another package, matrix-profile-foundation, about a year ago when I was working on my first paper. I was so busy, and unfortunately I couldn't follow up with them. However, I realized they fix it and name it neighbors. However, I think regime might be better since neighbors might be a little misleading.

3- The fraction output is not sorted. I think it is better to sort it in descending order and accordingly provide the other outputs for the same order.

I don't know if this (below) is correct but were you thinking of something like:

I think I made a mistake here. Sorry for that. It seems the top-k motif is already sorted in your output. I will take a look at the paper and the code again.

[NimaSarajpoor]

Usually, I'm looking to see if proposed changes would break backwards compatibility and affect other users who are using the code.

You are right. Otherwise it can become some sort of pain for former users. Will keep that in mind. Thanks for letting me know.

[seanlaw]

Yes, that would be a nice output. As you said, we can make it a 3-dim array to include the snippets_regimes for all snippets.
If I may ask, what did you mean when you said "Would you mind confirming if this would work?" ? Should I apply it to a time series to see if the output makes sense?

Yes, that's what I meant.

Maybe the problem is caused by something else.

Why do you think it's a problem? I guess one should assess how much each snippet covers the time series?

[NimaSarajpoor]

Yes, that's what I meant.

I will check it out.

I guess one should assess how much each snippet covers the time series?

I should have better used the word "meaningful" when I was talking about snippets. I think you are referring to the fraction values when you said assessing how much each snippet covers the time series. Am I right? If yes, then it can help us when we get the result for k=3. However, it doesn't help us when k=2. I expect the algorithm to provide the two distinct subsequences when we only consider k=2 snippet.

However, as you said, STUMPY doesn't try to be everything for everyone, and the output might be good enough. Maybe I am thinking subjectively here and the result does make sense.

[seanlaw]

I should have better used the word "meaningful" when I was talking about snippets. I think you are referring to the fraction values when you said assessing how much each snippet covers the time series. Am I right? If yes, then it can help us when we get the result for k=3. However, it doesn't help us when k=2. I expect the algorithm to provide the two distinct subsequences when we only consider k=2 snippet.

I am curious why that is? For k=2 I expect the algorithm to generate the two snippets (let's call these s1 and s2) that most represent the data (measured by fraction). I understand that it looks like s1 and s2 are quite similar but I would guess that the combination of s1 and s3 OR s2 and s3 would make up a smaller fraction than s1 and s2. Are you able to confirm this? Or do you find that the fraction for either s1ands3ORs2ands3is actually LARGER thans1ands2`? Looks can be deceiving (and we've added random noise and damping to the data, which can affect the fraction) so I would like to let the numbers help us quantify.

[NimaSarajpoor]

According to the result I got for k=3 snippets, the percentages of s2 and s3 are the largest (please see picture below)

=====================================

but I would guess that the combination of s1 and s3 OR s2 and s3 would make up a smaller fraction than s1 and s2.

Are the s1 and s2 snippets you mentioned in the end of your sentence above are the ones you get when k=2 as you mentioned in the first sentence of your paragraph? If yes, since the summation of their fraction values is 1, it is going to be larger than summation of the fraction values of any other two snippets in k=3.

In fact, the summation of the fraction values of all snippets is 1 for any k snippets. That's why the fraction values of k=2 doesn't help us to detect it.

=====================================

I should note that a low fraction value doesn't necessarily mean that the corresponding pattern is redundant. According to the paper, the authors stated that the changes in the areas can help the user with finding the proper k. I plotted it for different number of snippets of T[100: ] and I got the following figure:

As you can see, it is not so obvious which k is good.

=====================================

NOTE:

There is no random state in the function warp_add_noise. So, sometimes I can get correct snippets with k=2 and sometimes I cannot.
Should the noise have such impact on the result?

[seanlaw]

Should the noise have such impact on the result?

I don't know. Maybe.

As you can see, it is not so obvious which k is good.

Ahhh, that makes sense. Since there is no big drop in the plot of k vs area, this sort of suggests that the input parameters (namely m) might be bad? So, one might want to determine the "best" m first (I am working on the pan matrix profile for this) before computing snippets. Perhaps, that is the real take away message here or maybe I'm not able to understand your point. What do you think?

[NimaSarajpoor]

@seanlaw:

I don't know. Maybe.

I will investigate it (by checking out the MPdist profile values)

this sort of suggests that the input parameters (namely m) might be bad?

Could you please elaborate on this? Finding the optimal window size (m) might help; however, I think this observation is not because of that. Because k=2 already works on T but not on T[100: ].

I tried to remove the first 50 elements of the time series and this is what I got:

Here, it suggests k=2 is good, but again the discovered patterns are not satisfactory in my perspective.

Is it possible that the code considers only multipliers of m as the indices of subsequence, it misses some opportunities in finding the correct pattern?

[seanlaw]

I also received this response from the first author:

You might think of snippets as finding two clusters in the data. If you remove part of the data then there is no guarantee to find the same snippets in your data. A very simple example is if your snippet is part of the data that you are removing then you cannot find it after you remove that part.
Also, when I compute the snippet, I use non-overlapping window size (computation reason) now if you cut part of the data, then if your previous snippet was at position 200 and you remove 50 points of data, then your snippet cannot be found in position 150(150= 200-50).

[NimaSarajpoor]

will be quite expensive.

I didn't know about that. It seems the first author have the same opinion and didn't consider overlapping scenario due to computation reason.

one would need to play around more with the snippets input parameters.

So, I change mpdist_percentage from 0.05 (default) to 0.5, and it can discover the representative patterns.

Do you have any suggestions for how one might address this without significantly increasing the computational time?

At this moment, I have no solid knowledge on mpdist itself. I should put some time to go through the mpdist paper to understand how that parameter affects the final answer in snippets.

If what you said about the input parameters is true (i.e. they can resolve this issue for any other cases/data sets), then one can do grid search and compare some scores(?) or visually see the output. Problem Solved!

Another potential solution is to define a step, gama, which just gives the flexibility to the user to consider overlapping subsequences but with gama step. So, after the subsequence T[i:i+m], the subsequece T[i+gamma:i+gamma+m] will be considered. So, the user can try gamma=100, then 25, then, 20 while acknowledging that it increases the computation time. (Of course, it doesn't guarantee to provide the true snippet as it depends on the gamma!)

I can keep our discussion in mind and put it into my to-do list 😃 and think about it when I get some time in the upcoming days.

[seanlaw]

Instead of gamma, I would just call it step as it is a more intuitive/obvious parameter name. gamma sounds like a statistical value. Of course, where possible, I try to be consistent with the original publication parameter names but this would be a new parameter. In my opinion, this is likely a low priority item.

[NimaSarajpoor]

@seanlaw

I agree. The user can probably resolve this by trying different settings.

Just to confirm, we are probably not going to add the parameter due to its low priority and better to keep the algorithm provided by the author. Right?

So, I can add a small section in the tutorial to inform users of such a potential issue(?). What do you think?

=======================================================

For now, I can go and check out that snippet_regime indices and will update you shortly.

[seanlaw]

Frankly, I don't think it is necessary to draw attention to it. One difference between tutorials and an academic publication is that we don't need to spend time pointing out all of the flaws. Reviewers care about pros/cons. Instead, we are highlighting all of the useful things and users just want to know how to use the tool. I usually like to allow future user feedback to drive modifications to the basic tutorials.

Regarding snippets_regimes, sounds good and looking forward to it! No rush

[NimaSarajpoor]

@seanlaw

Thank you for providing more context/information. I came across a solution in the past where you can pass in a mask and have it return the start (inclusive) and stop (exclusive) indices where the mask is True

import numpy as np

def _get_mask_slices(mask):
    m1 = np.r_[0, mask]
    m2 = np.r_[mask, 0]
    starts, = np.where(~m1 & m2)
    ends, = np.where(m1 & ~m2)
    return np.c_[starts, ends]

if __name__ == '__main__':
    mask = np.array([ True,  True, False,  True,  True,  True,  True, False, False, False], dtype=bool)
    slices = _get_mask_slices(mask)
    print(slices)
    # array([[0, 2],
    #        [3, 7]])

This is a two column start and (exclusive) stop indices and, perhaps, you could append a third column with k? Maybe we should call these snippets_regimes in the code. Would you mind confirming if this would work?

Additionally, any changes to snippets.py would also need to happen in aampdist_snippets.py (this is our non-normalized version - i.e., without z-normalization).

I tried to check it out and this is what I got for the snippets of the toy data in the notebook:

Please note that a subsequence of length m = 200 and starting index: 910 is a common subsequence between two snippets in the figure above.

The regimes is the additional output of the snippets module. It is a list where the i-th item contains the slices of indices of the i-th snippet.

=========================================================================
I was wondering if you could help me with the following process:

I think one approach to check out the functionality of _get_mask_slices(mask) on the snippets is to copy the whole snippet.py module into one of the Jupyter notebook cells and then modify the script and get the results. However, I realized I had to change some lines of the script where you did relative imports. Am I right? And, I should confess that it was confusing, and I didn't put enough time into understand what I should do to call them from the notebook (located in docs) properly!

So, instead, I created the Snippets_Regime branch, which is the sub_branch of Snippets_Tutorial branch; and then, modify snippets.py and aampdist_snippets.py to include the new function and return the snippets_regime in addition to the previous ones. Then, I resolve the issues raised by flake8 (e.g. removing the blank line after the docstring of a function). Next, I run the code from Jupyter, as I thought it would call the module locally from the directory I am working with. But, I got the error regarding the expected number of outputs. So, I did ./setup.sh and it works. However, I skipped ./test.sh because I did it once and got error since the unit test is written for the previous version of the snippets' modules.

I would truly appreciate if you could clarify a few things for me:

1. Is that a proper approach to create sub_branch of the development branch? I couldn't create a new branch under main because the updated version of notebook is in the Snippets_Tutorial branch and not in the main. If it is not the proper approach, should I have first merged my development branch with main and then created another branch (for new PR) under main to modify snippets module? If my approach was okay, then is that okay to merge the changes to my development branch and then push it to the remote repo? (In this case, you will see the changes in the modules and the additional plot of regimes in the notebook. Right?)

2. Although we install stumpy from local directory, the python will call it from a library. That's why I had to do ./setup.sh. Right? So, is that okay to skip the unit test for now? (and, add later if the result makes sense)

3. Can I copy the module into a jupyter cell and modify the module there (and maybe change the name of the snippet function) to instantly see the result? (Is there any way to resolve the issue of relative imports written in the beginning of the modules?)

Sorry for my long questions. I am reading some articles currently to get a better idea of what's going on. I would truly appreciate if you could help me with the aforementioned process/questions.

[seanlaw]

This is a good question. Please allow me some time to put together a proper response

[seanlaw]

The regimes is the additional output of the snippets module. It is a list where the i-th item contains the slices of indices of the i-th snippet.

So, I like to avoid outputting a list whenever possible. Instead, I think that the output should be a three column array where column 1 is the (zero based) ith snippet, the second column is the (inclusive) start index, and the third column is the (exclusive) stop index. The first column might have repeating values if the ith snippet can be found in multiple regions of the data:

[
[0,   0, 200],
[0, 400, 600],
[1, 200, 400],
]

In this example, the first snippet (zero) is found twice in different places in the time series.

I was wondering if you could help me with the following process:
I think one approach to check out the functionality of _get_mask_slices(mask) on the snippets is to copy the whole snippet.py module into one of the Jupyter notebook cells and then modify the script and get the results. However, I realized I had to change some lines of the script where you did relative imports. Am I right? And, I should confess that it was confusing, and I didn't put enough time into understand what I should do to call them from the notebook (located in docs) properly!
So, instead, I created the Snippets_Regime branch, which is the sub_branch of Snippets_Tutorial branch; and then, modify snippets.py and aampdist_snippets.py to include the new function and return the snippets_regime in addition to the previous ones. Then, I resolve the issues raised by flake8 (e.g. removing the blank line after the docstring of a function). Next, I run the code from Jupyter, as I thought it would call the module locally from the directory I am working with. But, I got the error regarding the expected number of outputs. So, I did ./setup.sh and it works. However, I skipped ./test.sh because I did it once and got error since the unit test is written for the previous version of the snippets' modules.

So, you did the right thing by branching. Consider branching as a way to test out ideas and it is normal to branch off of a branch as you've done. However, when one of your ideas works out then you'll want to merge that code back to the parent branch (and eventually up to the main branch). Tracking branches can become complicated so try to only work on 1-2 branches at a time as a best practice before you become more comfortable. Otherwise, you may end up making changes to the wrong branch :(

Is that a proper approach to create sub_branch of the development branch? I couldn't create a new branch under main because the updated version of notebook is in the Snippets_Tutorial branch and not in the main. If it is not the proper approach, should I have first merged my development branch with main and then created another branch (for new PR) under main to modify snippets module?

You have the correct approach! Well done.

If my approach was okay, then is that okay to merge the changes to my development branch and then push it to the remote repo? (In this case, you will see the changes in the modules and the additional plot of regimes in the notebook. Right?)

Yep! That's exactly it. So, you'll likely create a sub-branch when the child branch depends on a parent branch. However, let's say you decide to work on a completely new_feature unrelated to snippets (but you haven't completed the snippets work yet), then you won't sub-branch your snippets_branch to add the new feature. Instead, you'll branch from your main branch in your local fork. In other words the new_feature branch will be at the same "tree depth" as your snippets_branch.

Although we install stumpy from local directory, the python will call it from a library. That's why I had to do ./setup.sh. Right? So, is that okay to skip the unit test for now? (and, add later if the result makes sense)

Yes, that's right. And, yes, for now you can skip the unit tests.

Can I copy the module into a jupyter cell and modify the module there (and maybe change the name of the snippet function) to instantly see the result? (Is there any way to resolve the issue of relative imports written in the beginning of the modules?)

This is not preferred since, as you pointed out, would require relative imports and things just get uglier. I will confess that I will do it this way sometimes when I'm lazy but it is definitely not the correct approach. What you've described above with sub-branching is preferred.

[NimaSarajpoor]

@seanlaw
Thanks a lot for your thorough response. Truly appreciate it.

So, I like to avoid outputting a list whenever possible. Instead, I think that the output should be a three column array where column 1 is the (zero based) ith snippet, the second column is the (inclusive) start index, and the third column is the (exclusive) stop index. The first column might have repeating values if the ith snippet can be found in multiple regions of the data:

I tried to use array, but I couldn't figure it out how I should do so, as the number of slices for each snippet might be different. Now, I can see your approach. I will change the code accordingly.

So, you'll likely create a sub-branch when the child branch depends on a parent branch.

In fact, the child branch (Snippets_Regime) doesn't depend on the parent branch (Snippet_Tutorial). I mean, I could have modified the snippets' modules in a completely new branch. Right? However, then I should have shown the result in the previous version of the Tutorial notebook. Can I merge my Snippets_Tutorial branch to the main branch on my local repo on PC? Then, I can create the branch Snippets_Regime under the main. Will that lead to conflict later, since my main branch is not the same as the origin?

[seanlaw]

In fact, the child branch (Snippets_Regime) doesn't depend on the parent branch (Snippet_Tutorial). I mean, I could have modified the snippets' modules in a completely new branch. Right?

Yes, if Snippets_Regime does not at all depend on Snippet_Tutorial then you could have modified the snippets module in a completely new branch (branched off of your own local main). This is a pretty routine thing to do.

However, then I should have shown the result in the previous version of the Tutorial notebook. Can I merge my Snippets_Tutorial branch to the main branch on my local repo on PC? Then, I can create the branch Snippets_Regime under the main. Will that lead to conflict later, since my main branch is not the same as the origin?

I think this is why it is somewhat important to git fetch (or rebase) to your local main and then continue to pull those changes down to your branches where necessary. The only time there may be major conflicts is if somebody else is modifying the same files as you are OR if somebody changes a function that you depends on. If not, then there will rarely be a conflict. This is also where unit tests really help to ensure that things aren't broken.

I hope I understand correctly but I think you can also make a copy of your modified snippets.py (from Snippets_Regime under Snippet_Tutorial) and put it in your desktop (outside of Git control), then create a new Snippets_Regime branch under main, and then move the copy from your desktop to this new branch (under mai) and commit this new file to this new branch. Then you can go and delete the branch of Snippets_Regimethat underSnippet_Tutorial`.

Does that help? If not, please feel free to ask!

[NimaSarajpoor]

@seanlaw:

A potential bug(?)

If you look at the regimes of snippets I provided in the figure in my previous comment, you will see that the last index of the regime of the second snippet is 1801. But, the length of the whole time series is 2000 (and the last index is 1999). So, if 1801 is the beginning of the snippet, then the last subsequence is Time Series [1801:(1801+m)] (where m=200); however, that will go beyond the index range. So, I had to use m-1 (or we can do min(2000, 1801+m).

But, my point is this:

shouldn't the last index be 1800 ?

[seanlaw]

Nope. In NumPy, the indexing uses an "inclusive" start index (i.e., it includes that index position) but an "exclusive" stop index (i.e., it excludes that index position). So, take this example:

import numpy as np

x = np.arange(10)
print(x[0:5])  # This will print [0, 1, 2, 3, 4] and NOT [0, 1, 2, 3, 4, 5]

Additionally, assuming you use a positive index value, NumPy will automatically stop you from going OVER the end of an array. So, following our example above:

print(x[0:10])  # This will print [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
print(x[0:11])  # This will print [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
print(x[0:100])  # This will print [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
print(x[0:1000])  # This will print [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
print(x[0:100000000000])  # This will print [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

So, your question about T[1801 : 1801 + m] will always stop at the last element in the time series (in this case, 1999). Does that make sense? Even if you accidentally did T[1801 : 1801 + m + 10000000], this would end with 1999 (inclusive).

[NimaSarajpoor]

if name == 'main':
mask = np.array([ True, True, False, True, True, True, True, False, False, False], dtype=bool)
slices = _get_mask_slices(mask)
print(slices)
# array([[0, 2],
# [3, 7]])

Sorry, My bad!!! I knew about the inclusive/exclusive indexing. However, I got confused about the output of the function _get_mask_slices(). I thought the second number in each slice provided as the output is a starting position as well.

My apologies!

[seanlaw]

Ahh, got it.