Stumpi pre-calculated mp

[tomweb153]

Hi,

I have a question regarding the use of a pre-calculated mp for the stumpi function. For me, it isn't clear how to use the "mp" parameter.

My use case:

• Calculate mp with stumpy.stump for data (1 week of data) -> initially
• Save all information from mp in database table
• Collect a days worth of data
• Load existing mp from database -> every day
• Update existing (1 week long) mp with stumpi -> egress=True
• Search if anomalies have occured in new data
• Save updated mp (still 1 week long) in database table

My restrictions:

• Tools for fully online mp calculation are not available
• Have to rely on database (as data only arrives there)
• Data can only be retrieved day by day (not as it arrives)

Is it possible to update an existing mp that contains information about 1 week of data with 1 day of new data (trying to prevent having to calculate the whole mp every time new data arrives)?
Your documentation says that the pre-calculated mp has to be the same length as the data one wants to calculate anything for. So what does the mp parameter exactly do, when passing a pre-calculated mp from day 1 and new data from day 2 (seems like it is not possible to update something of different lengths).

Here is a very short implementation demonstrating what I'm trying to do. The vertical line in the first plot marks where updatedata starts. I was expecting the mp to have changed after the stumpi usage.

Thank you for your help.
Best Regards
Tom

[seanlaw]

@tomweb153 Thank you for your question and welcome to the STUMPY community. When you say "Collect a days worth of data", can you please clarify if that is a single data point (representing the whole day) or, say, 24 data points (one data point for each hour of that day)?

[tomweb153]

Hello @seanlaw
Thanks for the warm welcome.
I have a dynamic sensor that collects about 2880 readings per day (2 readings per minute, if there is uninterrupted communication). So when I talk about a day of data I mean 2880 readings.

To clarify my use case further:

• 1 week of data is 20160 data points
• 1 day of data is 2880 data Points
• I want to incrementally update my mp
• I only want to keep 1 week of calculated mps

[seanlaw]

Thank you for the extra information @tomweb153

Following what you showed above, let's say you already have one weeks worth of data, you pre-compute the matrix profile for it:

import numpy as np
import stumpy

past_week_data = np.random.rand(20160)
m = 30
mp = stumpy.stump(past_week_data, m, normalize=False)

Hopefully, there is nothing unusual here but please feel free to ask any questions. Now, let's set up our matrix profile stream:

stream = stumpy.stumpi(past_week_data, m, normalize=False, egress=True)

Our stream has been initialized and now we can pass single/individual datapoints to it to update the matrix profile. Let's say we have some new data:

day_one_data = np.random.rand(2880)  # or retrieve from database

The stream has no concept of 2880 data points at this time. It only understands a SINGLE data point so let's update our matrix profile one data point a time:

for i in range(len(day_one_data):
    stream.update(day_one_data[i])

And (I think) that's it! For the second day, you'd do:

day_two_data = np.random.rand(2880)  # or retrieve from database
for i in range(len(day_two_data):
    stream.update(day_two_data[i])

Certainly, you don't need to call the new data day_one_data and day_two_data since, again, the stream only cares about the next single data point that you feed into it in order to update the matrix profile. Please let me know if that makes sense or if you have any further questions.

In case you missed it , I also recommend looking at this tutorial.

[tomweb153]

@seanlaw
Thank you for your response. I have been following your tutorials on the stumpy webpage for the past few months. They have been incredibly helpful, and I want to express my gratitude once again!

I have already considered a solution similar to what you suggested.
However, I believe I would need to initialize a streaming object with data from one week prior to updating my mp, right? If I understand this correctly, it would involve repetitive calculations because I cannot initiate the streaming object, update it, and then wait for new data, all while keeping the streaming object ready.

The challenge here is that I receive daily updates at a specific time from the database where my sensor readings are stored. Additionally, I have multiple sensors for which I need to perform anomaly detection (although I didn't mention this earlier, I'm dealing with approximately 10,000 sensors feeding data into the database). So I am currently searching for a effective solution to do anomaly detection on my sensor readings.
I was hoping that the option to provide the stumpi function with a mp, would solve my problem.

My ideas are mainly based on the following paper, published in 2023 by Eamonn Keogh and his colleagues: "Matrix Profile XXVIII: Discovering Multi-Dimensional Time Series Anomalies with K of N Anomaly Detection". In this context 5..7 sensors can be seen as a single system.

I had hoped to find an efficient way to calculate the mp for each sensor (or system) over the course of one weeks data, save these calculated mp in my database, and simply update them as new data becomes available. It turns out that one week provides the optimal amount of data required for secure anomaly detection within the given systems.

If you are available over the next few days, would you be open to have a video call to further discuss this? Please let me know if you are available, and I can send you an invitation via Teams or Zoom. I'll ensure the timing aligns with my location in Leipzig, Germany and your location, to minimize any inconveniences.

[seanlaw]

Thank you for your response. I have been following your tutorials on the stumpy webpage for the past few months. They have been incredibly helpful, and I want to express my gratitude once again!

Thank you for your kind words. It means a lot!

However, I believe I would need to initialize a streaming object with data from one week prior to updating my mp, right? If I understand this correctly, it would involve repetitive calculations because I cannot initiate the streaming object, update it, and then wait for new data, all while keeping the streaming object ready.

@tomweb153 I'm not sure I fully understand what problem you are facing so please elaborate. There is nothing stopping you from initiating the stream with np.nan data if the data doesn't already exist:

import numpy as np
import stumpy

past_week_data = np.full(20160, np.nan)
m = 30
mp = stumpy.stump(past_week_data, m, normalize=False)

Now, some folks have mentioned using cloudpickle to pickle this stream object and then storing it either locally (on disk) or into some database (I have not tried this myself and so I cannot comment). And then, when the new data arrives, retrieve the last stream object, load it into Python, retrieve the NEW data, and then update the stream accordingly. Note that since the initial "fake" data is np.nan, those will never produce a valid distance for consideration into the matrix profile.

So, Once the first m=30 points have been used to update the stream, then you will start having actual distances in your matrix profile after the 31st data point. Certainly, you'll need to wait for 20160 data points before you are guaranteed to have had a full weeks of data in your matrix profile.

Is this what you mean? If not, please clarify.

Note that at some point, we are moving beyond STUMPY and matrix profiles and going into architecture and production scalability of your problem, which is (sadly) beyond the scope of the package. I imagine that there may be something you can pull together by putting your 10K streaming objects into some sort of cloud storage (e.g. S3) and then create a cloud function (e.g., Lambda) that accepts a new day's data and location of the corresponding streaming object and updates the stream accordingly. I have little to no experience here :(

If you are available over the next few days, would you be open to have a video call to further discuss this? Please let me know if you are available, and I can send you an invitation via Teams or Zoom. I'll ensure the timing aligns with my location in Leipzig, Germany and your location, to minimize any inconveniences.

Unfortunately, due to my full time job and family, I am unable to offer personal consultations. However, you may continue posting your questions here in case others may have some suggestions too.

[tomweb153]

@seanlaw

Unfortunately, due to my full time job and family, I am unable to offer personal consultations. However, you may continue posting your questions here in case others may have some suggestions too.

I fully understand that!

Thank you for the suggestion to use cloudpickle. I think storing the streaming object may be the key for solving my problem.

I'm not sure I fully understand what problem you are facing so please elaborate.

Sorry, English is not my main language, but I will try to elaborate my problem:

• I have a lot of sensors (5..7 combined are ONE system) that provide me with 2 readings/minute
• I want to check whether some of the systems are experiencing anormal states
• I am very limited on how I can implement an anomaly detection for the sensors
  • no infrastructure for continuous calculations using stumpi (so that I could access the readings as they arrive)
• I thought of using data of one week (e.g. 11. Oct to 18. Oct) for an initial calculation of all the neccessary mps
  • do anomaly detection on the mps (e.g. use the left mp as stated in your tutorial)
  • save the mps into a database table (seems I have to save the streaming object!)
  • wait 1 day, so I can access the data that was collected for the last day (e.g. 19. Oct)
  • get the streaming object out of the database
  • update my streaming object with the newly available data (streaming object with egress=True )
  • check the updated mps for anomalies
  • save updated streaming object back to database
  • -> repeat this process every day (I'd call it "pseudo-streaming-anomaly-detection")

I hope, I could clarify what I'm trying to do.

[seanlaw]

Frankly, at some point, you need to ask yourself if it takes MORE total time to:

1. pickle the stream object
2. upload/download the stream object
3. upload/download the new data
4. load the stream object into memory
5. update the stream and potentially incur a numba warm up cost

than it would take to just re-compute the matrix profile each time.

[kaganilter]

Hi Sean,

I have been using the STUMPY for a while, and your work is very impressive—thank you! Similar to a previous question, we are planning to use STUMPY libraries in our streaming device time series data application. The data stream period is one hour.

Our plan is to calculate one year of historical data matrix profiles for each device using stumpy.stumpi and create daily inferences for upcoming data streams and merge/update into old matrix profiles that can be used an anomaly detection tasks. It looks like recalculating matrix profiles using stumpy.stump is faster than incremental updates by stumpy.stumpi. Since we have more than 100K devices, this performance comparison is critical for us.

Below is a sample code snippet of our use case. We are looking forward to hearing your comments and guidance on this:

import numpy as np
import time
from stumpy import stump, stumpi

# Generate a random dataset with 10,240 points
data = np.random.rand(10240)
m = 50  # window size

# Approach 1: Using stumpi for incremental updates
# Step 1: Calculate stumpi for the first 10,000 points
initial_data = data[:10000]
stream_data = data[10000:]

start_time = time.time()
mp_stream = stumpi(initial_data, m)

# Step 2: Incrementally add the rest of the 240 points
for point in stream_data:
    mp_stream.update(point)

stumpi_time = time.time() - start_time
print(f"Time taken using stumpi: {stumpi_time:.4f} seconds")

# Approach 2: Using stumpy.stump for the entire dataset
start_time = time.time()
mp_full = stump(data, m)
stumpy_time = time.time() - start_time
print(f"Time taken using stumpy: {stumpy_time:.4f} seconds")

# Compare the times
print(f"Incremental update time: {stumpi_time:.4f} seconds")
print(f"Direct computation time: {stumpy_time:.4f} seconds")

and this is the output of the above code:

Time taken using stumpi: 7.0623 seconds
Time taken using stumpy: 0.0818 seconds
Incremental update time: 7.0623 seconds
Direct computation time: 0.0818 seconds

Looking forward to your insights.

Best regards,

Khan

[NimaSarajpoor]

Great job spotting that!

Thanks! I found it by just checking different things. I should have done profiling ! will do it.

Anyways, I forget if the order in the two for-loops matter. That is, is each the "searchsorted" and "shift-and-insert" step completely independent (i.e., has no effect) of the others? If so, then I wonder if we may be able to use prange there?

In the first for-loop, the iterations are independent from each other (each iteration corresponds to one row of 2D array). So, we can use prange. One thing that got my attention is that each row in P has length k. So, in the default setting, when k=1, each row has length 1. Therefore, np.searchsorted searches value in an array of length 1 when we have the default setting k=1.

In the second for-loop, however, we just update the last row of a 2D array in each iteration. So, I think we cannot use prange there.

Additionally, since _shift_insert_at_index is already using njit, I wonder if we're incurring ? Realistically, I wouldn't expect np.searchsorted to be that slow?? Where is the bottleneck?

Didn't come to my mind! Good point! Will do profiling and share my findings.

[NimaSarajpoor]

@seanlaw

Where is the bottleneck?

After profiling, it turns out the second for-loop is responsible for the majority of the running time.

stumpy/stumpy/stumpi.py

Lines 372 to 376 in 5892869

	for i, d in enumerate(D):
	if d < self._P[-1, -1]:
	idx = np.searchsorted(self._P[-1], d, side="right")
	core._shift_insert_at_index(self._P[-1], idx, d)
	core._shift_insert_at_index(self._I[-1], idx, i + self._n_appended)

The reason is that, in the second for-loop, the number of iterations that leads to execution of "searchsorted + two shift-insert" is much higher (~ 10x) than the number of iterations in the first-loop. In case you are interested in the "why": I think it is because in the first-loop we are just updating the NN of all previous subsequences by considering the last, new subsequence. It should be fair to assume that the chance of having many subsequences whose NN is that last, new subsequence is low (particularly when the input consists of random values). In the second for-loop, however, we are trying to find NN of the last subsequence. So, it is fair to say that this can take more effort as we need to update the best-so-far NN as we traverse the subsgeuqneces from the start till the very end.

For this input

seed = 0
np.random.seed(seed)
T = np.random.rand(10000)
m = 50 

stream = stumpi(T, m)
r = np.random.rand(1000)
# for i in range(len(r)):
#    out = stream._update_egress(r[i]) 

On average, the first for-loop in _update_egress has 1 iteration, and the second for-loop has 10 iterations (that pass the inner if condition).

[seanlaw]

In the second for-loop, however, we just update the last row of a 2D array in each iteration. So, I think we cannot use prange there.

Actually, I wonder if we can split the distance profile up into n_thread chunks and have each thread find the NN within its chunk and then perform a reduction across threads thereafter?

[NimaSarajpoor]

Actually, I wonder if we can split the distance profile up into n_thread chunks and have each thread find the NN within its chunk and then perform a reduction across threads thereafter?

I noticed that the non-parallel version is fast(er) for finding top-k in array with length upto 1M.

---

The following is the current version of finding top-k smallest values in D (to fill P[-1, :])

 for i, d in enumerate(D): 
     if d < self._P[-1, -1]: 
         idx = np.searchsorted(self._P[-1], d, side="right") 
         core._shift_insert_at_index(self._P[-1], idx, d) 
         core._shift_insert_at_index(self._I[-1], idx, i + self._n_appended) 

Let's ignore the self._n_appended for now, and wrap a njit-decorated function around it.

@njit
def func1(D, k):
    P = np.full(k, np.inf, dtype=np.float64)
    I = np.full(k, -1, dtype=np.int64)
    for i, d in enumerate(D):
        if d < P[-1]:
            idx = np.searchsorted(P, d, side="right")
            core._shift_insert_at_index(P, idx, d)
            core._shift_insert_at_index(I, idx, i)


    return P, I

In the parallel version, we find top-k in each chunk, and then merge top-k of chunks.

@njit(parallel=True)
def func2(D, k, P_threads, I_threads, tmp_P, tmp_I):
    n_threads = P_threads.shape[0]
    step = int(math.ceil(len(D) / n_threads))
    
    P_threads[:, :] = np.inf
    I_threads[:, :] = -1
    for i in prange(n_threads):
        start = i * step
        stop = min(start + step, D.shape[0])
        for j in range(start, stop):
            if D[j] < P_threads[i, -1]:
                idx = np.searchsorted(P_threads[i], D[j], side="right")
                core._shift_insert_at_index(P_threads[i], idx, D[j])
                core._shift_insert_at_index(I_threads[i], idx, j)

    P = P_threads[0]
    I = I_threads[0]
    tmp_P[:] = np.inf 
    tmp_I[:] = -1 

    # using the idea used in `core._merge_topk_PI`
    for i in range(1, n_threads):
        overlap = set(I_threads[i]).intersection(set(I))
        aj, bj = 0, 0
        idx = 0
        tmp_P[:] = np.inf
        tmp_I[:] = -1
        for _ in range(2 * k):
            if idx >= k:
                break
            if bj < k and P_threads[i, bj] < P[aj]:
                if I_threads[i, bj] not in overlap:
                    tmp_P[idx] = P_threads[i, bj]
                    tmp_I[idx] = I_threads[i, bj]
                    idx += 1
                bj += 1
            else:
                tmp_P[idx] = P[aj]
                tmp_I[idx] = I[aj]
                idx += 1
                aj += 1

        P[:] = tmp_P
        I[:] = tmp_I

    return P, I

I tested their performance.

k = 1
n_threads = numba.config.NUMBA_NUM_THREADS

P_threads = np.full((n_threads, k), np.inf, dtype=np.float64)
I_threads = np.full((n_threads, k), -1, dtype=np.int64)

tmp_P = np.full(k, np.inf, dtype=np.float64)
tmp_I = np.full(k, -1, dtype=np.int64)

for n in [1000, 10000, 100000, 1000000, 1000000, 10000000]:
    print(f'length of D --> {n}')
    
    lst_func1 = []
    lst_func2 = []
    
    for _ in range(100):
        D = np.random.rand(n)
    
        t1 = time.time()
        P_ref, I_ref = func1(D, k)
        t2 = time.time()
        P_comp, I_comp = func2(D, k, P_threads, I_threads, tmp_P, tmp_I)
        t3 = time.time()
        
        lst_func1.append(t2 - t1)
        lst_func2.append(t3 - t2)
    
        np.testing.assert_allclose(P_ref, P_comp)
        np.testing.assert_allclose(I_ref, I_comp)


    print(f"njit non-parallel version --> {np.mean(lst_func1[1:]):.6f} seconds")
    print(f"njit parallel version --> {np.mean(lst_func2[1:]):.6f} seconds")
    print('----------------------------------------------------------')

And got the following result:

# k = 1

length of D --> 1000
njit non-parallel version --> 0.000004 seconds
njit parallel version --> 0.000337 seconds
----------------------------------------------------------
length of D --> 10000
njit non-parallel version --> 0.000007 seconds
njit parallel version --> 0.000559 seconds
----------------------------------------------------------
length of D --> 100000
njit non-parallel version --> 0.000036 seconds
njit parallel version --> 0.000281 seconds
----------------------------------------------------------
length of D --> 1000000
njit non-parallel version --> 0.000335 seconds
njit parallel version --> 0.000485 seconds
----------------------------------------------------------
length of D --> 1000000
njit non-parallel version --> 0.000338 seconds
njit parallel version --> 0.000492 seconds
----------------------------------------------------------
length of D --> 10000000
njit non-parallel version --> 0.003280 seconds
njit parallel version --> 0.001895 seconds
----------------------------------------------------------

And, for k=100, got the following result:

# k = 100

length of D --> 1000
njit non-parallel version --> 0.000120 seconds
njit parallel version --> 0.000506 seconds
----------------------------------------------------------
length of D --> 10000
njit non-parallel version --> 0.000166 seconds
njit parallel version --> 0.000584 seconds
----------------------------------------------------------
length of D --> 100000
njit non-parallel version --> 0.000243 seconds
njit parallel version --> 0.000662 seconds
----------------------------------------------------------
length of D --> 1000000
njit non-parallel version --> 0.000611 seconds
njit parallel version --> 0.000912 seconds
----------------------------------------------------------
length of D --> 1000000
njit non-parallel version --> 0.000601 seconds
njit parallel version --> 0.000917 seconds
----------------------------------------------------------
length of D --> 10000000
njit non-parallel version --> 0.003660 seconds
njit parallel version --> 0.002362 seconds
----------------------------------------------------------

---

Next step can be profiling _update_egress and see what parts contribute the most to the computing time AFTER wrapping a njit-decorator function around the for-loops.

[seanlaw]

I'm getting the sense that simply njit-ing things is already "fast enough" and that parallelizing the second for-loop will add unnecessary complexity to the code. @NimaSarajpoor if you'd like to propose a PR for simply njit-ing the function (no parallelization) then I can review it