Fix Precision in Unit Test

[seanlaw]

This unit test causes a failure:

import numpy as np
import numpy.testing as npt
import naive
from stumpy import stump, config
import pandas as pd

def test_stump_identical_subsequence_self_join():
    seed = 27343
    np.random.seed(seed)

    identical = np.random.rand(8)
    T_A = np.random.rand(20)
    T_A[1 : 1 + identical.shape[0]] = identical
    T_A[11 : 11 + identical.shape[0]] = identical
    m = 3
    zone = int(np.ceil(m / 4))
    ref_mp = naive.stump(T_A, m, exclusion_zone=zone, row_wise=True)
    comp_mp = stump(T_A, m, ignore_trivial=True)
    naive.replace_inf(ref_mp)
    naive.replace_inf(comp_mp)
    npt.assert_almost_equal(
        ref_mp[:, 0], comp_mp[:, 0], decimal=config.STUMPY_TEST_PRECISION
    )  # ignore indices

    comp_mp = stump(pd.Series(T_A), m, ignore_trivial=True)
    naive.replace_inf(comp_mp)
    npt.assert_almost_equal(
        ref_mp[:, 0], comp_mp[:, 0], decimal=config.STUMPY_TEST_PRECISION
    )  # ignore indices

[NimaSarajpoor]

@seanlaw
I think I found the source of error. According to my investigation, the problem is in the precision of $\rho$. The discrepancy is in index 2 and 12. In fact S1 = T[2:2+m] and S2=T[12:12+m] are exactly the same and the $\rho$ between them must be 1.

However, according to the calculation of stumpy.stump, we can see that the highest value in $\rho[:, 2, 0]$ (print it in this line) is 0.9999999999483034.

idx = 2 # the start index of subsequence of interest
arg = np.argmax(ρ[:, idx, 0]) 
>>> I[arg, 2, 0]
12 # this is NN to the subsequence at index idx=2

>>> ρ[arg, idx, 0]
0.9999999999483034

A potential solution:
We may prevent this by creating config.STUMPY_RHO_PRECISION and use it in the line of code I mentioned above before comparing the $\rho$ values, calculated in each thread, with each other (and then set it to 1 if the calculated value is too close to 1 (?!))

[NimaSarajpoor]

@seanlaw
I do not know if this discrepancy (due to the numerical error) can happen in other cases as well or not... like when the true (exact) $\rho$ is not 1 (not identical) but we have difference between stumpy.stump and naive.py/stump due to numerical error.

---

How to handle identical case without using something like config.STUMPY_RHO_PRECISION? I mean... if we try to add STUMPY_RHO_PRECISION to config.py and use it to replace close-to-1 $\rho$ values to 1, then it may wrongly impact other cases.

So, how about just handling identical cases? The following solution works for me BUT it is not efficient for massive data, say 10M data points. I think large-size data has higher chance of hash collision and thus increase the computing time (?)

# in core.py 

def preprocess_diagonal(T, m):
   ...
   # add a new part to find duplicates
   identical = {}
   for i in range(T.shape[0] - m + 1):
       identical[tuple(T[i : i + m])] = []  
   for i in range(T.shape[0] - m + 1): 
       identical[tuple(T[i : i + m])].append(i) 
    
   # O(N) time complexity
   T_subseq_IdenticalCode = np.full(T.shape[0] - m + 1, -1, np.int64)
   for i, (key, item) in enumerate(identical.items()): 
       for idx in item:
           T_subseq_IdenticalCode[idx] = i


   # return T_subseq_IdenticalCode in addition to the other outputs

And, then we check if two subsequences are identical in stumpy._compute_diagonal

# in stumpy._compute_diagonal
if T_A_subseq_IdenticalCode[i + k] == T_A_subseq_IdenticalCode[i]:
     pearson = 1.0

# or we can do it at the end to correct the values of matrix profile

---

An alternative way is to let collisions happen to make the process faster, and check it later.

# I think it O(N) time complexity

l = T.shape[0]  - m + 1
subsequence_hash = np.empty(l , dtype=np.int64)
for i in range(l ):
    subsequence_hash[i] = hash(tuple(T[i:i+m])) 

# there is a chance that we might have collision but that is okay as we will check it out later.

Then, we can correct the matrix profile (i.e. matrix profile/ matrix profile left / matrix profile right) at the end:

# I think it O(N) time complexity having assumed that the z-norm distance between two subsequences is O(1).

for row in range(l):
   idx = row
   for col in range(3): 
       idx_nn = I[row, col] # use matrix profile index to find NN of idx
       if subsequence_hash[idx] == subsequence_hash[idx_nn]:
             # two subsequences might be identical since their hash is the same
             if np.all(Tznorm[idx : idx + m] - Tznorm[idx_nn : idx_nn + m]) == 0): # check if two subsequences are actually identical
                  P[row, col] == 0

Is it worth it? we are increasing time in pre-processing and post-processing steps. (I believe both of these steps can be parallelized via prange)

[seanlaw]

@NimaSarajpoor I don't think hash is a good idea as it is slow and pushing subsequences into a Dict may cause a lot of memory to be used as the length of the time series increases. Additionally, we'd need to account for z-normalization in these cases. Just because two subsequences aren't identical doesn't mean that their z-normalized subsequences aren't identical.

I don't know if it's worth addressing or what the best option is. Setting a config.STUMPY_RHO_PRECISION (or some other name like config.STUMPY_PERFECT_CORRELATION = 0.9999999) might not be the worst idea. At least, people can set it it to 1.0 if they dislike the precision that we've chosen and then we'd only need to compare if $\rho$ > config.STUMPY_PERFECT_CORRELATION: and then shift it to 1.0` or something.

Alternatively, I wonder if we can check $\rho$ and if it is above the config.STUMPY_PERFECT_CORRELATION then we should compute it directly using the dot product instead of using the distance from the previous diagonal.

[NimaSarajpoor]

I was wondering if what you have in mind is to check each single $\rho$ that is being calculated in stumpy.stump._compute_diagonal. If yes, then in that case, the alternative approach

Alternatively, I wonder if we can check and if it is above the config.STUMPY_PERFECT_CORRELATION then we should compute it directly using the dot product instead of using the distance from the previous diagonal.

might (hugely ?) increase the computing time if the data mostly show periodical behavior. (because, we will have many identical subsequences)

---

On the other hand, I think we also want to make sure the values in final matrix profile are precise. I mean...if we say idx_nn = I[idx], then dist(T[idx:idx+m], T[idx_nn, idx_nn+m]) should be exact.

So, I think we can just re-calculate distance at the end after checking $\rho > config.STUMPY\_PERFECT\_CORRELATION$.
(I know we should better re-calculate it throughout the stumpy.stump._compute_diagonal but I think this is too much for a rare case like this and, as I mentioned earlier, we may increase computing time when there are many identical subsequences. So, just correcting matrix profile at the end should be okay)

[seanlaw]

I think the only thing we'd need to change is add a line after L188 in _compute_diagonal:

https://github.com/TDAmeritrade/stumpy/blob/576de379994896063be71d264e48a08ae17476cf/stumpy/stump.py#L183-L188

if T_B_subseq_isfinite[i + k] and T_A_subseq_isfinite[i]: 
    # Neither subsequence contains NaNs 
    if T_B_subseq_isconstant[i + k] or T_A_subseq_isconstant[i]: 
        pearson = 0.5 
    else: 
        pearson = cov * Σ_T_inverse[i + k] * σ_Q_inverse[i] 
        if pearson > config.STUMPY_PERFECT_CORRELATION:
            pearson = 1.0

I think that this is the simplest and should be a negligible cost as there is an if comparison (which should be cheap??) and then a value is set (not computed).

This should be simple, easy to maintain, easy to interpret, and if users don't like it then config.STUMPY_PERFECT_CORRELATION is easy to understand

I think we also want to make sure the values in final matrix profile are precise

I don't know. Unfortunately, I don't think we'll get much more "precise". It's a trade off between speed and a slight loss in precision (but not much, actually). In the worst case, if you have a long sine wave (time series) OR when each subsequence has a perfect match, you'll end up computing all of the pairwise distances twice (once during _stump and then once during post-processing if you do it your way! I believe that in my example above, it should be reasonably efficient but one should verify.

[seanlaw]

As I think about it more, if we were to do this as a post-processing step, we would write a separate function to stumpy.stump. I would not add the post-processing step to the existing function. I think my proposal above is the way to go

[NimaSarajpoor]

I think that this is the simplest and should be a negligible cost as there is an if comparison (which should be cheap??) and then a value is set (not computed).

Correct. As you said, I think is better and enough for such rare case (you do not have to recalculate distances and you modify the values throughout the process).

[seanlaw]

As a cross reference to #655 this related test is also failing:

def test_stumpi_identical_subsequence_self_join_egress():
    m = 3

    seed = np.random.randint(100000)
    seed = 84451
    np.random.seed(seed)

[NimaSarajpoor]

@seanlaw
I think you may want to see this! Apparently, it is not just for the identical subsequences. So, I created a test where I randomly scale up (* 1000) or scale down (0.001) chunks of time series.

Here is the test function:

def test_stump_imprecision():
    seed = np.random.randint(100000)
    np.random.seed(seed)

    T = np.random.rand(64)
    
    # randomly choose some chunks of time series and scale them up or down 
    T_mask = np.full(len(T), 0, dtype=bool)
    k = np.random.randint(len(T) // 4, len(T))
    IDX = np.random.choice(np.arange(len(T)), k, replace=False)
    T_mask[IDX] = True
    slices = core._get_mask_slices(T_mask)
    scale = np.random.choice(np.array([0.001, 1000]), len(slices), replace=True)
    for i, (start, stop) in enumerate(slices):
        T[start:stop] = scale[i] * T[start:stop]

    m = 3
    zone = int(np.ceil(m / 4))
    ref_mp = naive.stump(T, m, exclusion_zone=zone)
    comp_mp = stump(T, m, ignore_trivial=True)
    naive.replace_inf(ref_mp)
    naive.replace_inf(comp_mp)

    ref_P = ref_mp[:, 0].astype(np.float64)
    comp_P = comp_mp[:, 0].astype(np.float64)

    npt.assert_almost_equal(ref_P, comp_P)

And, here is the error:

 Mismatched elements: 37 / 62 (59.7%)
E       Max absolute difference: 0.22458425
E       Max relative difference: 0.78388078
E        x: array([2.6611068e-02, 7.4557520e-05, 5.2291006e-02, 3.3167208e-02,
E              1.0766382e-04, 4.5937140e-02, 4.4588591e-02, 3.0992438e-04,
E              1.8665326e-01, 4.5937140e-02, 7.9446847e-02, 1.1810141e-01,...
E        y: array([2.6611068e-02, 8.7222889e-05, 7.4846320e-02, 3.3167208e-02,
E              1.3051016e-04, 4.5937140e-02, 7.4861289e-02, 3.0992438e-04,
E              1.8665326e-01, 4.5937140e-02, 7.9446847e-02, 1.1810142e-01,...

[seanlaw]

Hmm, so what do you think is the issue? I'm assuming that naive.stump is "correct"?

[NimaSarajpoor]

I'm assuming that naive.stump is "correct"?

Yeah...I think so. I took a look at the naive.stump and it just do z-norm on each single subsequence and then calculate the pair-wise distances. So, I assume naive.stump should be correct.

---

According to my investigation, I think one problem is the imprecision in the calculation of sliding standard deviation Σ_T, particularly when we have large numbers. (Just for the record: Even if we take std of each subseq (instead of rolling std), we can reduce the error significantly, but we still have error!)

---

One way to make it a little bit better is standardization on the whole time series T (so, we scale std of each subseq by std(T)), and then rescale back.

Another way that I think should work (I tested it to some extent) is that we refine the welford nanvar/nanstd. We can do it by calculating the rolling nanvar/nanstd for the subsequences each scaled by their minmax, and then scale back at the end.

# in stumpy/core.py
@njit
def rolling_minmax_scale(a, m):
    """
    Docstring
    """
    l = len(a) - m + 1
    scale = np.empty(l, dtype=np.float64)
    for i in range(l):
        scale[i] = np.max(a[i : i + m]) - np.min(a[i : i + m])
       if scale[i] == 0:
          scale[i] = 1
     
    return scale

Then, we use this to modify the welford as follows:

def _welford_nanvar(a, w, a_subseq_isfinite):
    """
    Compute the rolling variance 
    """
    all_variances = np.empty(a.shape[0] - w + 1, dtype=np.float64)
    prev_mean = 0.0
    prev_var = 0.0

    scale = rolling_minmax_scale(a, w) # NEW

    for start_idx in range(a.shape[0] - w + 1):
        prev_start_idx = start_idx - 1
        stop_idx = start_idx + w  # Exclusive index value
        last_idx = start_idx + w - 1  # Last inclusive index value

        if (
            start_idx == 0
            or not a_subseq_isfinite[prev_start_idx]
            or not a_subseq_isfinite[start_idx]
        ):
            curr_mean = np.nanmean(a[start_idx:stop_idx])
            curr_var = np.nanvar(a[start_idx:stop_idx] / scale[start_idx]) # MODIFIED: minmax-scaled
        else:
            curr_mean = prev_mean + (a[last_idx] - a[prev_start_idx]) / w

            curr_var = (
                prev_var * np.square(scale[prev_start_idx] / scale[start_idx])
                + ((a[last_idx] - a[prev_start_idx]) / scale[start_idx])
                * (
                    (a[last_idx] - curr_mean + a[prev_start_idx] - prev_mean)
                    / scale[start_idx]
                )
                / w
            ) # MODIFIED: reversing scale and scale again! # So, our curr_var is `var` of minmax-scaled subseq

        all_variances[start_idx] = curr_var

        prev_mean = curr_mean
        prev_var = curr_var

    return all_variances # this is now var of subsequences, where each subseq is scaled by minmax

And finally, reverse everything back to normal at the end:

def welford_nanstd(a, w=None):
    """
    This a convenience wrapper around `welford_nanvar`.
     """
    if w is None:
        w = a.shape[0]

    a_subseq_isfinite = rolling_isfinite(a, w)
    scale = rolling_minmax_scale(a, w) # NEW

    out = np.sqrt(np.clip(welford_nanvar(a, w), a_min=0, a_max=None))
    out[a_subseq_isfinite] = scale[a_subseq_isfinite] * out[a_subseq_isfinite] # scale back to normal

    return out

---

For instance, I am comparing the proposed one with the existing version:

• imprecision test mentioned above: test_stump_imprecision (seed = 0) --> existing version: max abs diff: 0.36033988, proposed welford approach: 0.09287705
• for the test function for identical case (but not scaling random chunks by 1000 or 0.001): --> existing version: resolved by new config variable, proposed welford approach: no need to add config. It is all good!

[NimaSarajpoor]

@seanlaw
In another test that I designed, the proposed welford performed worse. So, you may ignore the idea (I just did not delete it for our future ref)