STIMP result mechanism

[aycabayraktar]

Hi community,
I have a question about the class STIMP which is implementation of SKIMP algorithm. For better understanding, I used the dataset in the official SKIMP Dataset page: Termitte DNA https://sites.google.com/view/pan-matrix-profile/datasets
The dataset length is 16326

1. I used stimp and get pmp results.
2. I checked the length of outputs respectively. len(pmp.M_) len(pmp.PAN_). They have same length that is what I expected.
3. I want to control the M_values. pmp.M_[:100]. The largest value is 13060. Why this value is greater than half of the dataset, it is closed to the length of the dataset.
   
4. My other question is, let's say I picked the first element of M_ array (M_[0]. As I understood, this is subsequence length. Then I investigate it in PAN_. I would expect to see the length of 1 values in the related PAN_ array as M_[0]*2.
   

[seanlaw]

Hi @aycabayraktar!

Firstly, if you haven't had a chance, I would recommend reading the original SKIMP paper as it contains a lot of details that have been captured in stumpy.stimp.

I want to control the M_values. pmp.M_[:100]. The largest value is 13060. Why this value is greater than half of the dataset, it is closed to the length of the dataset.

There is a lot to unpack here so I'll try to summarize the important points!

For a time series, T with a length n = 16326, the stimp function will attempt to compute the matrix profiles for 3 <= m <= 13060. You may ask, "where does the upper bound of 13060come from?!? It's actually computed from thecore.get_max_window_size` function:

https://github.com/TDAmeritrade/stumpy/blob/3cb836ba88d942314035ac74410e6128b7dd33ab/stumpy/core.py#L528C1-L553C17

So, the one thing to keep in mind is that the maximum (theoretical) allowable window size actually can't be the full length of the time series due to the fact that we need to always account for an (left and right) exclusion zone. To demonstrate this, we can naively set m = 13060 and directly compute the matrix profile:

import numpy as np
import stumpy

T = np.random.rand(16326)
m = 13060
print(stumpy.stump(T, m))

mparray([[161.58358094149827, 3266, -1, 3266],
         [inf, -1, -1, -1],
         [inf, -1, -1, -1],
         ...,
         [inf, -1, -1, -1],
         [inf, -1, -1, -1],
         [161.58358094149827, 0, 0, -1]], dtype=object)

Notice that only the FIRST and LAST subsequences has "enough space" to have a nearest neighbor. If we simply increase our window size by one to m = 13061, we'll get:

mparray([[inf, -1, -1, -1],
         [inf, -1, -1, -1],
         [inf, -1, -1, -1],
         ...,
         [inf, -1, -1, -1],
         [inf, -1, -1, -1],
         [inf, -1, -1, -1]], dtype=object)

And no subsequence has any nearest neighbors because they have all been excluded since all nearest neighbors (subsequences) would lie within their respective exclusion zones!

My other question is, let's say I picked the first element of M_ array (M_[0]. As I understood, this is subsequence length. Then I investigate it in PAN_. I would expect to see the length of 1 values in the related PAN_ array as M_[0]*2.

Perhaps, there is a bit of misunderstanding here. I am wondering if this tutorial (still a work-in-progress) might help?

[aycabayraktar]

Hi @seanlaw ,
Based on paper, I am trying to verify PMP return. In section 3, the author describes both calculation and heat map creation (with DNA example). Also in the next section, for example, brute force approach addressing to create PMP with looping over a range and calculating MP's and storing them in the rows of PMP. My point is, how we achieve to the heat map with matrix profiles described in here.
In SCIMP pmp.PAN_ is a 0-1 matrix so we can plot based on it, but how we achieved pmp.PAN_? I hope, I could explain.

III. COMPUTING THE PAN MATRIX PROFILE

[seanlaw]

I think all you need to do is:

import numpy as np
import stumpy

T = np.random.rand(1000)
pmp = stumpy.stimp(T, percentage=1.0)

for _ in range(len(pmp.M_)):
    pmp.update()

Then, you will find all of your UNALTERED matrix profiles in pmp._PAN (and all of the transformed matrix profiles in pmp.PAN_). You can do whatever you want with all of the raw matrix profiles in pmp._PAN (i.e., you can plot them).

[aycabayraktar]

According to stimp documentation:

1. This parts generates the final shape inf full _PAN matrix
   self._PAN = np.full( (self._M.shape[0], self._T.shape[0]), fill_value=np.inf, dtype=np.float64 )

2. Then in the update section, _PAN is updated based on the parameters. For %100 percentage, the else part, uses the stump. For example, in the first step of the loop: self._PAN[ self._bfs_indices[0], : out[:, 0].shape[0]] = out[:, 0]
   This is what I expect. But my point is, why I can't get the same results for the selection index by using stump itself.
   ` See Table 2

        if self._n_processed < self._M.shape[0]:
            m = self._M[self._n_processed]
            if self._percentage < 1.0:
                approx = scrump(
                    self._T,
                    m,
                    ignore_trivial=True,
                    percentage=self._percentage,
                    pre_scrump=self._pre_scrump,
                    k=1,
                    T_A_subseq_isconstant=self._T_subseq_isconstant_func,
                )
                approx.update()
                self._PAN[
                    self._bfs_indices[self._n_processed], : approx.P_.shape[0]
                ] = approx.P_
            else:
                out = self._partial_mp_func(
                    self._T,
                    m,
                    ignore_trivial=True,
                    T_A_subseq_isconstant=self._T_subseq_isconstant_func,
                )
                self._PAN[
                    self._bfs_indices[self._n_processed], : out[:, 0].shape[0]
                ] = out[:, 0]
            self._n_processed += 1

1. I mean based on the 1000 length T that you have provided, why pmp._PAN[0] differs from m = 401 stump = stumpy.stump(T, m) stump[:,0]

[seanlaw]

I mean based on the 1000 length T that you have provided, why pmp._PAN[0] differs from m = 401 stump = stumpy.stump(T, m) stump[:,0]

I think I see the problem. pmp._PAN stores the matrix profiles ordered from smallest m to largest m. However, the matrix profiles are NOT computed in this order and, instead, are computed in "breadth-first-search" (BFS) order. Thus, in order to retrieve the correct BFS order, we must follow the procedure in the pan transformation function (shown in the permalink above) but this example should work:

import numpy as np
import stumpy

T = np.random.rand(1000)
pmp = stumpy.stimp(T, percentage=1.0)
pmp.update()

# retrieve the BFS indices for the matrix profiles that have already been processed
idx = pmp._bfs_indices[: pmp._n_processed]

# retrieve the specific raw matrix profiles from the pan matrix profile and make a copy
# Note that `_PAN` is a private variable and should NEVER be accessed directly as there is no guarantee that it will not change!!
processed_mps = pmp._PAN[idx].copy()  # Don't do this!!

# Compute the first matrix profile by hand
mp = stumpy.stump(T, pmp.M_[0])

# Compare `mp[:, 0]` and `processed_mps[0, :len(T)-pmp.M_[0]+1]`

I think what you've helped me identify is that it is somewhat difficult to retrieve the raw matrix profiles from the pmp object. I will take a look at whether there is a more convenient way to do this. Perhaps, it would make sense to add another class attribute for something like pmp.P_ that represents the raw matrix profiles but ordered by BFS?

[seanlaw]

I've created a new issue to address this. Thank you for your patience @aycabayraktar!

[aycabayraktar]

You're welcome, I really thank you for your effort @seanlaw