Newbie at Stumpy (motif in a live timeserie) #680

silentrain23 · 2022-09-16T12:03:13Z

silentrain23
Sep 16, 2022

Hi guys,
many thanks for this wonderful library.
I am trying to make it work for a project that has:

live data via MQTT of a flow of 5 samples/second (to be done)
training period when the flow is inited (length of the training time depends on the flow, can be 10 seconds or 2 minutes) (to be done)
motif recognition during the flow (hard to solve!)

this is the motif (a standard one, taken from a previous sampling) I am going to find in the time-serie

this is the time serie

this is the code I am trying

%matplotlib inline

import pandas as pd
import stumpy
from sklearn import preprocessing
import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

sample = np.array([649,642,652,650,642,648,646,641,648,643,644,635,643,634,639,640,649,644,651,651,654,642,639,649,650,644,624,631,648,649,643,647,654,649,634,627,630,632,631,638,632,637,644,649,649,641,1112,1095,1076,1052,632,644,635,654,629,631,634,638,640,639,645,662,651,644,643,647,642,633,631,626,635,649,647,646,646,650,644,640,638,645,647,639,627,635,638,633,651,647,652,643,644,647,641,631,636,642,647,623,639,1057,1065,1061,1074,1078,1076,1082,1082,1079,1076,1080,1069,1061,1063,1075,1061,1067,1059,1070,1058,1069,1057,1050,1045,1062,1056,1063,1064,1079,1087,1083,1078,1067,1073,1053,1089,1063,1072,1074,1071,1065,1032,1027,1047,1040,1033,1017,916,815,804,799,825,1078,1075,1055,1053,1059,1074,1075,1068,1074,1073,1066,1070,1081,1070,1060,1037,926,828,778,778,1059,1089,1081,1087,1093,1096,1085,1086,1075,1079,1062,1038,1040,1053,1050,1059,1065,1070,1073,1072,1053,1050,630,632,631,638,631,627,635,635,628,630,633,648,648,647,643,646,643,647,633,640,630,628,631,1077,1082,1078,757,637,814,920,642,1093,644,638,635,628,639,642,639,829,672,1088,1067,1136,667,650,645,643,642,641,634,643,636,631,632,638,629,629,633,639,642,647])

original_readings = np.array([648,651,653,651,641,646,648,640,640,642,635,647,643,635,646,646,1093,1080,1083,1080,1087,1098,1091,1092,1078,1101,1096,1088,1073,1070,1067,1067,1069,924,654,648,646,633,639,640,636,636,643,650,650,646,653,1104,1071,1067,1067,1087,1076,1076,1056,1069,1069,1053,1057,1083,1083,1078,1086,1087,1080,1064,1055,1050,1050,1058,1059,1072,1075,1078,1079,1085,1087,1085,1088,1070,1072,1076,1070,1072,1076,1072,1076,1074,1072,1070,1077,1067,1081,1085,1082,1084,1081,1092,1089,1079,1071,1063,1060,1064,1085,1088,1095,1095,1086,1077,1077,1071,1067,1047,1061,1057,1061,1067,1064,1059,1065,1071,1054,1041,1006,859,825,804,760,1050,1069,1070,1067,1057,1066,1078,1087,1090,1088,1096,1102,1077,1087,1070,1045,994,871,796,785,978,1086,1071,1087,1064,1069,1067,1067,1071,1057,1068,1063,1054,1055,1069,1055,1064,1081,1049,1071,1069,1089,1078,1080,1087,1085,1080,1074,1081,1088,677,637,639,637,643,651,650,654,653,661,652,648,636,643,640,634,640,644,652,644,645,648,652,641,642,647,644,636,634,640,641,633,631,642,639,634,638,638,649,653,649,646,652,654,645,641,639,634,641,647,641,647,641,645,642,649,645,643,653,646,645,647,641,643,641,648,639,643,648,653,649,652,659,658,649,651,639,645,636,638,635,640,641,649,652,646,658,648,643,646,656,646,642,648,642,647,648,650,670,648,641,642,651,643,630,634,647,654,648,652,660,662,654,654,654,651,646,637,637,639,805,1073,1065,1039,663,912,919,651,1104,1117,1053,653,638,649,646,638,648,1025,1066,1074,1097,1097,1093,1083,1091,1069,1063,1054,1068,1063,1080,1074,1086,1094,1098,1093,753,640,640,644,638,630,638,638,637,637,641,645,650,657,659,651,638,639,644,645,638,632,638,648,644,648,657,659,654,652,649,651,921,1086,1098,1089,1082,1059,1065,1059,1063,1075,1087,1081,1088,1098,1083,1086,1087,1055,1070,1076,1081,1080,1071,1076,1080,1065,1074,1060,1065,1063,1062,1065,1088,1098,1093,1092,1093,1093,1086,1077,1061,1048,1069,1063,1063,1071,1081,1079,1092,1105,1091,1094,1096,1092,1091,1088,1066,1063,1084,1071,1066,1068,1058,1055,1049,1047,1066,1054,1060,1072,1077,1086,1087,1088,1081,1068,1069,1056,1055,1068,1069,1054,1064,1066,1069,1067,1075,1073,1080,1067,1076,1074,1068,1061,1074,1081,1076,1080,1087,1094,705,661,1086,1077,1066,1061,1050,1043,1045,1061,1068,1070,1065,1071,1081,1074,1053,1049,970,862,818,814,801,814,810,807,822,805,805,811,820,811,800,800,1059,1052,1067,1052,1059,1053,1069,1079,1072,1082,1084,1083,1082,1087,1088,1083,1072,697,655,645,644,644,651,646,651,658,645,652,644,641,662,650,631,644,654,635,642,650,656,659,648,649,650,647,640,637,637,640,633,640,778,1057,1067,1063,1078,1083,1096,1079,1095,1067,1082,1088,1073,1078,1076,1056,1075,660,661,650,641,640,645,645,633,641,636,649,1070,1046,1061,1074,1083,1079,1067,1081,1082,1067,1068,1078,1068,1067,1078,1083,1088,1085,1072,1064,1072,1064,1067,1056,1055,1075,1072,1081,1080,1064,1058,1065,1073,1066,1059,1058,1053,1071,1069,1073,1089,1069,1055,1071,1067,1064,1028,909,794,789,800,922,1063,1072,1089,1087,1091,1089,1087,1075,1077,1108,1083,1062,1067,1069,1044,1009,902,822,788,777,1073,1069,1094,1079,1080,1065,1073,1055,1053,1042,1059,1058,1067,1074,1083,1081,1088,1088,1076,1077,1070,1055,1056,1059,1054,1064,1066,1071,1076,689,644,646,646,656,646,640,645,638,632,638,645,655,646,649,1066,1087,1087,679,1092,1096,652,700,637,1106,668,1092,1062,1060,1070,1065,1065,1058,1059,1056,1063,1074,1072,1070,1077,1075,1088,1075,1069,1066,651,644,639,641,649,653,650,653,653,661,658,642,641,642,646,637,631,639,640,643,639,648,662,648,646,653,653,657,879,1090,1094,1080,1082,1080,1075,1074,1066,1055,1067,1077,1074,1082,1097,1096,1088,1083,1082,1088,1083,1074,1064,1058,1050,1049,1060,1069,1070,1084,1085,1083,1097,1080,1071,1068,1069,1059,1074,1057,1055,1056,1069,1074,1063,1071,1052,1066,1063,1097,1099,1087,1085,1086,1068,1068,1067,1080,1062,1071,1079,1090,1094,1084,1076,1065,1071,1058,1052,1052,1062,1057,1073,1069,1070,1074,1066,1070,1070,1072,1065,1059,1066,1074,1071,1060,1073,1085,1078,1076,1078,1083,1071,1069,1057,1055,1068,1061,1064,1077,1080,1090,1094,1091,1082,1095,1083,1076,1056,1068,1059,1053,1054,1049,1068,1049,1063,1066,1073,1055,1052,1060,1037,972,859,815,800,795,796,797,801,783,784,796,809,813,811,815,812,1060,1083,1084,1068,1074,1057,1055,1056,1050,1050,1063,1060,1064,1079,1079,1069,1070,937,633,634,636,644,636,635,642,645,639,642,646,656,646,646,646,634,636,636,628,645,636,638,636,642,639,654,646,646,649,653,651,640,647,632,632,647,637,648,645,645,656,655,644,645,648,648,636,628,638,651,646,648,647,651,658,645,648,641,650,643,639,645,650,651,658,640,645,639,635,635,640,632,635,639,645,636,638,648,656,654,646,645,646,646,641,638,631,629,635,631,634,639,641,646,641,650,651,641,655,644,640,640,639,641,645,648,637,652,653,643,639,637,643,639,631,629,641,638,642,642,648,653,645,649,640,641,634,635,631,645,644,642,637,636,640,642,641,648,641,641,636,646,648,911,1078,1079,1069,1066,1061,1068,1067,1051,1065,1068,1085,1068,1079,1085,1083,1090,1078,1080,1068,1065,1048,1066,1055,1070,1050,1070,1061,1067,1076,1064,1059,1079,1073,1078,1062,1065,1061,1041,1048,1061,1060,1068,1071,1085,1077,1079,1081,1087,1069,1063,1052,1049,1060,1064,1062,1071,1082,1082,1073,1067,1072,1060,1044,1053,1079,1062,1071,1060,1059,1064,1063,1062,1077,1067,1060,1073,1059,1062,1040,1027,1029,1044,1033,1018,987,880,837,797,748,1077,1089,1079,1055,1055,1059,1058,1055,1047,1067,1065,1058,1050,1064,1067,1045,1003,858,798,764,955,1066,1051,1055,1049,1069,1068,1074,1083,1082,1079,1068,1060,1055,1044,1065,1042,1047,1058,1074,1073,1069,1077,1071,1071,1080,1078,1084,1071,1075,679,628,643,640,648,643,643,639,639,645,651,1078,1082,1089,750,632,684,885,647,634,875,663,633,632,653,1061,1062,1063,1056,1071,1084,1085,1093,1078,1086,1076,1076,1075,1073,1064,1056,1065,1072,774,647,641,640,649,638,636,638,645,633,632,638,640,636,632,631,636,639,651,667,649,651,654,645,641,637,633,634,631,634,1065,1072,1070,1070,1073,1079,1068,1102,1071,1067,1066,1085,1060,1058,1065,1078,1071,1075,1074,1089,1090,1089,1082,1079,1062,1067,1050,1049,1054,1044,1067,1060,1072,1079,1083,1069,1066,1087,1066,1070,1080,1067,1062,1046,1061,1071,1063,1057,1042,1090,1089,1092,1091,1104,1089,1079,1069,1069,1049,1062,1073,1078,1087,1075,1070,1077,1068,1063,1050,1057,1055,1064,1047,1056,1052,1048,1059,1056,1064,1062,1075,1075,1073,1069,1048,1046,1050,1061,1053,1058,1072,1072,1063,1069,1060,1084,1096,1092,1093,1063,1077,1064,1055,1055,1054,1060,657,641,652,1082,1063,1078,1059,1077,1074,1066,1061,1061,1052,1052,1050,1048,1048,1048,1031,989,863,827,818,823,812,807,804,801,797,795,796,783,791,789,795,1032,1075,1075,1074,1066,1081,1071,1069,1050,1068,1054,1065,1067,1077,1068,1081,1084,883,647,638,633,639,642,634,647,650,655,653,656,648,649,646,662,641,635,644,635,634,643,643,634,637,642,1031,1057,1054,1072,1076,1085,1075,1077,1083,1072,1070,1066,1075,1060,1052,1050,708,656,645,642,655,656,647,644,645,650,638,640,634,640,1086,1059,1055,1048,1060,1050,1050,1065,1049,1063,1062,1052,1060,1061,1061,1073,1065,1075,1063,1073,1067,1056,1058,1061,1052,1048,1046,1051,1055,1059,1073,1059,1076,1080,1081,1074,1075,1086,1068,1063,1064,1074,1061,1048,1060,1049,1066,1064,1067,1086,1078,1082,1069,1079,1057,1070,1078,1081,1063,1045,1049,1071,1072,1074,1074,1083,1072,1079,1075,1084,1077,1062,1053,1053,1042,1032,983,872,804,777,750,1070,1075,1081,1089,1077,1076,1078,1068,1059,1060,1064,1054,1056,1073,1075,1058,995,864,800,762,866,1041,1061,1066,1067,1078,1077,1084,1084,1087,1087,1079,1067,1053,1051,1052,1042,1052,1049,1029,1057,1057,1056,1052,1075,1062,1052,1073,1058,1056,1052,1088,1088,1068,1054,1054,670,640,643,641,651,652,649,650,650,654,652,649,640,649,644,635,635,1068,1068,1028,649,769,674,692,642,941,698,999,1069,1087,1080,1083,1054,635,642,638,630,633,1074,1069,1066,1059,1055,1056,1056,1061,1060,1070,1072,1078,1085,789,643,638,631,640,636,634,625,640,640,634,642,645,639,637,634,642,635,633,630,638,638,627,634,639,644,641,647,1079,1095,1089,1087,1081,1089,1074,1076,1061,1055,1075,1056,1067,1078,1070,1073,1061,1062,1062,1071,1062,1074,1055,1054,1055,1057,1064,1066,1061,1077,1071,1072,1067,1092,1086,1081,1079,1084,1070,1067,1066,1040,1042,1050,1056,1071,1070,1025,1081,1103,1078,1088,1079,1083,1057,1078,1075,1065,1067,1091,1070,1073,1080,1076,1063,1056,1065,1049,1051,1052,1055,1067,1068,1059,1079,1083,1060,1066,1068,1063,1062,1048,1053,1052,1045,1055,1062,1056,1071,1079,1065,1067,1070,1065,1064,1066,1080,1075,1075,1074,1090,1085,1071,1079,866,638,1113,1060,1056,1048,1045,1051,1047,1060,1077,1059,1075,1079,1079,1072,1053,1031,1030,929])

original_readings_p = preprocessing.minmax_scale(original_readings)
sample = preprocessing.minmax_scale(sample)

bins = np.arange(-.25, 2, .5)
original_readings_d = (np.digitize(original_readings_p, bins)-1)/2
sample = (np.digitize(sample, bins)-1)/2

original_readings_d = original_readings_d.astype(float)
sample = sample.astype (float)

#start STUMPY process

m = 310    #distance between the first dotted lines
mp = stumpy.stump(original_readings_d.astype(float), m)
window_size = m

    
motifs = stumpy.motifs(
    original_readings_d,
    mp[:, 0],
    max_matches=len(original_readings_d)//window_size
)
    
    
motif_idxs = motifs[1][0].tolist()
nearest_neighbor_idxs = [ mp[motif_idx, 1] for motif_idx in motif_idxs ]

print (nearest_neighbor_idxs)

print(f"The nearest neighbor is located at index {nearest_neighbor_idx}")

print(f"The motif is located at index {motif_idx}")
fig, axs = plt.subplots(2, sharex=True, gridspec_kw={'hspace': 0})
plt.suptitle('MOTIF DISCOVERY', fontsize='30')

axs[0].plot(original_readings_d)
axs[0].set_ylabel('READINGS', fontsize='20')
#rect = Rectangle((motif_idx, 0), m, 1.2, facecolor='lightgrey')
#axs[0].add_patch(rect)
#rect = Rectangle((nearest_neighbor_idx, 0), m, 1.2, facecolor='lightgrey')
#axs[0].add_patch(rect)
# draw the division of the pattern
axs[0].axvline(x=40, linestyle="dashed", color= "red")
axs[0].axvline(x=350, linestyle="dashed", color= "red")
axs[0].axvline(x=740, linestyle="dashed", color= "red")
axs[0].axvline(x=1250, linestyle="dashed", color= "red")
axs[0].axvline(x=1680, linestyle="dashed", color= "red")


for i in range(len(nearest_neighbor_idxs)):
    rect = Rectangle((nearest_neighbor_idxs[i], 0), 250, 40, facecolor='lightgrey')
    axs[0].add_patch(rect)
    axs[1].axvline(x=nearest_neighbor_idxs[i], linestyle="dashed")
    rect = Rectangle((nearest_neighbor_idxs[i], 0), 250, 40, facecolor='lightgrey')
    axs[0].add_patch(rect)

axs[1].set_xlabel('Time', fontsize ='20')
axs[1].set_ylabel('Matrix Profile', fontsize='20')
axs[1].axvline(x=motif_idx, linestyle="dashed")

#axs[1].axvline(x=nearest_neighbor_idx, linestyle="dashed")
axs[1].plot(mp[:, 0])
plt.show()

At the end, I expect to match 4 motifs (the dotted line is the slice of them in the timeserie) but I am getting very confused about.

Can you support me @seanlaw and team?

Thank u in advance

seanlaw · 2022-09-16T17:51:11Z

seanlaw
Sep 16, 2022
Maintainer

@silentrain23 I'm having a bit of trouble understanding what it is that you are trying to do. However, based on your post, it appears that you already have a known pattern that you are looking for. Perhaps, you can simply search for that known pattern directly?

Please see this tutorial

0 replies

silentrain23 · 2022-09-17T10:50:00Z

silentrain23
Sep 17, 2022
Author

Hi @seanlaw thank you for your kind reply.

With match or mass, i had these issues

I took the sample and did some modifications: enlarged by factor 2; added noise; divided by 2 the values

The algo can match the pattern just when it has the same length. Maybe there are some DTW to do on the algo?

Another point: I cannot understand the values of the idxs
k = 6
idxs = np.argpartition(distance_profile, k)[:k]
idxs = idxs[np.argsort(distance_profile[idxs])]

which is [1188 909 1187 910 908 1186]: shouldnt be just two values?

Anyway I tried to match my original_readings and sample arrays with poor results

Thank you for supporting
Emanuele

7 replies

seanlaw Sep 17, 2022
Maintainer

I am having good results now by applying a DTW on my original serie and sample, so I remove all the time-dependent variations on the same state (eg: 3 second at state 1)

That makes sense. Matrix profiles do not take into account any stretching/compression of the time series data. You may be able to perform uniform scaling OR DTW to stretch and compress your data but this is beyond the scope of matrix profiles.

Are there any parameters to do DTW with Stompy?

No, matrix profile has no concept of DTW.

A good thing might be to put a parameter in mass to avoid proximity indexes (some cutoff?)

No, stumpy.mass is a function to help you compute all pairwise distances (i.e., a "distance profile") quickly. What you choose to do with those distances is up to you. It doesn't make any sense to add a parameter here. Instead, you should use stumpy.match instead and adjust the config.STUMPY_EXCL_ZONE_DENOM instead which controls how the "exclusion zone" (i.e., the proximity indices) according to:

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

Then, all indices in the range T[i - excl, i + excl + 1] (surround index i) are ignored. Internally, stumpy.match calls stumpy.mass and handles the "proximity indices" for you. We recommend that you start with that first.

silentrain23 Sep 17, 2022
Author

works: I found the 7 patterns in my TS with mass (idxs: [ 19 796 348 631 186 915 472], when k=7; idxs: [ 19 796 348 915 186 631 472 349 632 187 916 473 797 20] when k=14). And this is fine.

I amt trying the match as you suggested but I cannot solve the issue with the STUMPY_EXCL_ZONE_DENOM just with trying values
'm = sample_d.size
print (f"sample size: {m}")

To set the exclusion zone to zero, set the denominator to np.`inf`

stumpy.config.STUMPY_EXCL_ZONE_DENOM = 0.4

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

print (f"EXCL ZONE: {excl_zone}")

#stumpy.config.STUMPY_EXCL_ZONE_DENOM = np.inf

matches_best = stumpy.match(
sample_d,
original_readings_d,
max_distance=np.inf, # set the threshold to infinity include all subsequences
max_matches=28, # find the top n matches

)
print (f"Found {matches_best.size/2} items:\n{matches_best[:,1]}")

stumpy.config.STUMPY_EXCL_ZONE_DENOM = 4 #'

sample size: 50
EXCL ZONE: 125
Found 7.0 items:
[19 796 348 186 631 474 998] which are not exactly the same as 'stumpy.mass' values. I got the same values if I leave all untouched as the sample code

seanlaw Sep 17, 2022
Maintainer

[19 796 348 186 631 474 998] which are not exactly the same as 'stumpy.mass' values. I got the same values if I leave all untouched as the sample code

They are not the same as stumpy.mass because you've set stumpy.config.STUMPY_EXCL_ZONE_DENOM = 0.4. If you set stumpy.config.STUMPY_EXCL_ZONE_DENOM = np.inf right before calling stumpy.match:

stumpy.config.STUMPY_EXCL_ZONE_DENOM = np.inf

matches_best = stumpy.match(
sample_d,
original_readings_d,
max_distance=np.inf, # set the threshold to infinity include all subsequences
max_matches=7, # find the top n matches
)

Then this should produce the same results as when using stumpy.mass(sample_d, original_readings_d)

Notice that the indices [19 796 348 186 631 474 998] does not contain any overlapping regions.

silentrain23 Sep 17, 2022
Author

Right!
Everything is running fine, the question is, as I have to put some fixed value, how to be sure 0.4 is fine as STUMPY_EXCL_ZONE_DENOM?

seanlaw Sep 18, 2022
Maintainer

At the end of the day, you need to make the decision on what the exclusion zone needs to be relative to your matching index:

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

The conventional wisdom is to choose T[i - m/2 : i + m/2 + 1] (i.e., half the length of your pattern before/after your matching index). This is usually reasonable. However, there is no "right" answer as this is data specific, so this is something you'll have to play around with on your own. I'd start with the excl_zone equation and choose what you'd want to exclude and then rearrange the equation to determine what value of config.STUMPY_EXCL_ZONE_DENOM you need to select. Again, this choice must be handled by the user. You have to decide whether you want overlapping (nearby) matches

silentrain23 · 2022-09-19T06:36:44Z

silentrain23
Sep 19, 2022
Author

Yes, @seanlaw this is the only manual setting I have to do. Also the "k" factor is to be defined, but, by now is fine for me.
My goal was to look for a way to do some ML on the data and then decide the "k" and excl_zone without any manual setting. But now is ok.

My effort are now on the live data analysis. Is STUMPI the right place to start?

1 reply

seanlaw Sep 19, 2022
Maintainer

Yes, for streaming/online data, stumpy.stumpi is the right place to start. Note that in that case, you compute a matrix profile and not a distance profile. If you only care about a specific pattern then you can simply proceed with what you've done above. However, if you care about how the matrix profile changes/evolves then you can use stumpy.stumpi.

silentrain23 · 2022-09-26T13:57:08Z

silentrain23
Sep 26, 2022
Author

Hallo, again @seanlaw
I am having trouble with another timeserie to validate the algo.
I know the (approximate) sample (6 values array), I know that in the timeserie (original_readings_d) I have 60 occurrencies, but stumpy.match gets more (0, 0, 0, 0, 0, 0) than sequences with higher values

This is the code

%matplotlib notebook
from ipywidgets import *

import stumpy
import numpy as np
import matplotlib.pyplot as plt

sample_d = np.array ([0. , 1. , 0. , 0.5, 0.5, 0. ])
original_readings_d = np.array([0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. ,
 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. ,
 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5, 0.5, 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0.5,
 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 1. , 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5,
 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5, 0.5,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. ,
 0. , 0.5, 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5,
 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. ,
 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. ,
 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0.5, 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 1. , 0. , 0. ,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 1. , 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. ,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. ,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. ,
 0. , 0. , 0.5, 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5,
 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0. , 0.5, 0. , 0. , 0. , 0. ,
 0. , 0. , 0.5, 0.5, 0.5, 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. ,
 0.5, 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. ,
 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0. , 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0.5, 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. ,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5, 0.5, 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5,
 0.5, 0. , 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0. , 0.5,
 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 1. , 0. , 0.5, 0.5, 0. , 0. , 0. ,
 0. , 0. , 0. , 0.5, 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ,
 0. , 0. ])

# AMPLIFY
np.set_printoptions(threshold=1000000)
original_readings_d = np.select([original_readings_d == 0.5], [1], original_readings_d)
sample_d = np.select([sample_d == 0.5], [1], sample_d)
print (sample_d)

m = sample_d.size
print (f"sample size: {m}")

# To set the exclusion zone to zero, set the denominator to np.`inf`
stumpy.config.STUMPY_EXCL_ZONE_DENOM = 0.5 #era 0.4 con azzurra#3
#stumpy.config.STUMPY_EXCL_ZONE_DENOM = 2 #

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

print (f"EXCL ZONE: {excl_zone}")

#stumpy.config.STUMPY_EXCL_ZONE_DENOM = np.inf

best_matches = stumpy.match(
    sample_d,
    original_readings_d,
    max_distance=np.inf, # set the threshold to infinity include all subsequences
    max_matches=70,      # find the top n matches
)

print (f"Found {best_matches.size/2} items:\n{best_matches[:,1]}")

plt.suptitle('stumpy.match', fontsize='30')
plt.xlabel('Time', fontsize ='20')
plt.ylabel('WT', fontsize='20')
plt.plot(original_readings_d)

for match in best_matches[:,1]:
    plt.plot(range(match, match+len(sample_d)), original_readings_d[match:match+len(sample_d)], lw=4)


stumpy.config.STUMPY_EXCL_ZONE_DENOM = 4 #back to defautl


print (sample_d)

The excl_zone is fine enough, but I cannot understand how to avoid the "6-zeros' matches.
Thank you

3 replies

seanlaw Sep 26, 2022
Maintainer

Your problem is not with STUMPY. The main issue is that your sample_d is closer in (z-normalized) Euclidean distance to [0., 0., 0., 0., 0., 0.] than to [0., 1., 1., 1., 0., 1.].

You can see this by computing the z-norm distances:

import numpy as np
import stumpy

# Compute distance between sample_d and 6-zeros
print(np.linalg.norm(stumpy.core.z_norm(sample_d) - stumpy.core.z_norm([0., 0., 0., 0., 0., 0.])))

# Compute distance sample_d and [0., 1., 1., 1., 0., 1.]
print(np.linalg.norm(stumpy.core.z_norm(sample_d) - stumpy.core.z_norm([0., 1., 1., 1., 0., 1.])))

Which gives:

2.449489742783178
3.4641016151377544

You have a few options. The first thing you can try to do is set any constant subsequence (of length len(sample_d)) to np.nan values:

threshold = 0.0
stddevs = np.nanstd(stumpy.core.rolling_window(original_readings_d, len(sample_d)), axis=1)
for idx, stddev in enumerate(stddevs):
    if stddev <= threshold:
        original_readings_d[idx:idx + len(sample_d)] = np.nan

This is done by computing the standard deviation of each subsequence (rolling window) and any subsequence that has a standard deviation of zero must be constant. And then we set those specific constant subsequence to np.nan and then we can use STUMPY as you've already done above:

# To set the exclusion zone to zero, set the denominator to np.`inf`
stumpy.config.STUMPY_EXCL_ZONE_DENOM = 0.5 #era 0.4 con azzurra#3
#stumpy.config.STUMPY_EXCL_ZONE_DENOM = 2 #

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

print (f"EXCL ZONE: {excl_zone}")

#stumpy.config.STUMPY_EXCL_ZONE_DENOM = np.inf

best_matches = stumpy.match(
    sample_d,
    original_readings_d,
    max_distance=np.inf, # set the threshold to infinity include all subsequences
    max_matches=70,      # find the top n matches
)

print (f"Found {best_matches.size/2} items:\n{best_matches[:,1]}")

plt.suptitle('stumpy.match', fontsize='30')
plt.xlabel('Time', fontsize ='20')
plt.ylabel('WT', fontsize='20')
plt.plot(original_readings_d)

for match in best_matches[:,1]:
    plt.plot(range(match, match+len(sample_d)), original_readings_d[match:match+len(sample_d)], lw=4)


stumpy.config.STUMPY_EXCL_ZONE_DENOM = 4 #back to defautl

silentrain23 Sep 27, 2022
Author

Thank you @seanlaw but to make it work, I had to simplify a lot sample_d from
[0., 1. , 0. , 1. ,1. ,0.]
to
[1., 0., 1., 1.]

now I can grab 65 readings instead of the 12 I had with a longer sample_d

Do you think it's correct?

seanlaw Sep 27, 2022
Maintainer

Unfortunately, I have no idea what you are trying to do. The choice of the window size is domain specific and it is up to you to decide whether it is "right". Anecdotally, a shorter window size usually makes little sense and a single digit window sizes tend to cause problems because it means that your pattern is too short to be unique.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Newbie at Stumpy (motif in a live timeserie) #680

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 4 comments 11 replies

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

Select a reply

Newbie at Stumpy (motif in a live timeserie) #680

silentrain23 Sep 16, 2022

Replies: 4 comments · 11 replies

seanlaw Sep 16, 2022 Maintainer

silentrain23 Sep 17, 2022 Author

seanlaw Sep 17, 2022 Maintainer

silentrain23 Sep 17, 2022 Author

To set the exclusion zone to zero, set the denominator to np.inf

seanlaw Sep 17, 2022 Maintainer

silentrain23 Sep 17, 2022 Author

seanlaw Sep 18, 2022 Maintainer

silentrain23 Sep 19, 2022 Author

seanlaw Sep 19, 2022 Maintainer

silentrain23 Sep 26, 2022 Author

seanlaw Sep 26, 2022 Maintainer

silentrain23 Sep 27, 2022 Author

seanlaw Sep 27, 2022 Maintainer

silentrain23
Sep 16, 2022

Replies: 4 comments 11 replies

seanlaw
Sep 16, 2022
Maintainer

silentrain23
Sep 17, 2022
Author

seanlaw Sep 17, 2022
Maintainer

silentrain23 Sep 17, 2022
Author

To set the exclusion zone to zero, set the denominator to np.`inf`

seanlaw Sep 17, 2022
Maintainer

silentrain23 Sep 17, 2022
Author

seanlaw Sep 18, 2022
Maintainer

silentrain23
Sep 19, 2022
Author

seanlaw Sep 19, 2022
Maintainer

silentrain23
Sep 26, 2022
Author

seanlaw Sep 26, 2022
Maintainer

silentrain23 Sep 27, 2022
Author

seanlaw Sep 27, 2022
Maintainer