unit testing failure for snippets

[NimaSarajpoor]

With python=3.10.11, numpy=2.0.2, numba=0.60.0 and llvm=0.43.0, the test function test_mpdist_snippets_s_with_isconstant fails for the following input:

seed = 1615
np.random.seed(seed)
T = np.random.uniform(-1000.0, 1000.0, [64])
m = 10
s = 3
k = 3

[NimaSarajpoor]

A similar error was observed for m=9, s=3, k=3.

OS: MacOS

python=3.12
numpy=2.1.3
numba= 0.61.0
llvm=0.44.0

And the error is:

_______________ test_mpdist_snippets_s_with_isconstant[3-3-9-T0] _______________

T = array([ 293.67128677,  308.19248823,  -90.46204887, -484.08056803,
       -815.73526712, -404.95991025, -591.06946033,...15185, -230.50234171, -360.47856789,  -41.299737  ,
       -955.80732529,   47.51877429,  372.92107172, -390.49869207])
m = 9, k = 3, s = 3

    @pytest.mark.parametrize("T", test_data)
    @pytest.mark.parametrize("m", m)
    @pytest.mark.parametrize("k", k)
    @pytest.mark.parametrize("s", s)
    def test_mpdist_snippets_s_with_isconstant(T, m, k, s):
        isconstant_custom_func = functools.partial(
            naive.isconstant_func_stddev_threshold, quantile_threshold=0.05
        )
        (
            ref_snippets,
            ref_indices,
            ref_profiles,
            ref_fractions,
            ref_areas,
            ref_regimes,
        ) = naive.mpdist_snippets(
            T, m, k, s=s, mpdist_T_subseq_isconstant=isconstant_custom_func
        )
        (
            cmp_snippets,
            cmp_indices,
            cmp_profiles,
            cmp_fractions,
            cmp_areas,
            cmp_regimes,
        ) = snippets(T, m, k, s=s, mpdist_T_subseq_isconstant=isconstant_custom_func)
    
>       npt.assert_almost_equal(
            ref_snippets, cmp_snippets, decimal=config.STUMPY_TEST_PRECISION
        )

tests/test_snippets.py:211: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/contextlib.py:81: in inner
    return func(*args, **kwds)
/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/contextlib.py:81: in inner
    return func(*args, **kwds)
/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/numpy/_utils/__init__.py:85: in wrapper
    return fun(*args, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

args = (<function assert_array_almost_equal.<locals>.compare at 0x10ce03a60>, array([[-336.77194671, -854.63885204, -784.0031....1080978 ,   92.0933612 ,
        -213.1082644 ,  684.15146167,  639.03341797,  221.1559157 ,
          70.8424083 ]]))
kwds = {'err_msg': '', 'header': 'Arrays are not almost equal to 5 decimals', 'precision': 5, 'verbose': True}

    @wraps(func)
    def inner(*args, **kwds):
        with self._recreate_cm():
>           return func(*args, **kwds)
E           AssertionError: 
E           Arrays are not almost equal to 5 decimals
E           
E           Mismatched elements: 9 / 27 (33.3%)
E           Max absolute difference among violations: 1230.1028783
E           Max relative difference among violations: 7.87042798
E            ACTUAL: array([[-336.77195, -854.63885, -784.00315, -230.50234, -360.47857,
E                    -41.29974, -955.80733,   47.51877,  372.92107],
E                  [ 552.12966,  652.35062,  -73.57151,  601.59769, -180.18051,...
E            DESIRED: array([[-336.77195, -854.63885, -784.00315, -230.50234, -360.47857,
E                    -41.29974, -955.80733,   47.51877,  372.92107],
E                  [ 552.12966,  652.35062,  -73.57151,  601.59769, -180.18051,...

/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/contextlib.py:81: AssertionError
======================== 1 failed, 146 passed in 35.82s ========================
Error: Test execution encountered exit code 1

---

[Update]
I was able to get a similar error with the following inputs:

seed=2636
np.random.seed(seed)
T = np.random.uniform(-1000, 1000, [64]).astype(np.float64)

m=9
k=3
s=3

[seanlaw]

Do we know why this is happening?

[NimaSarajpoor]

Do we know why this is happening?

According to my investigation, I think the issue is coming from naive. In the naive snippet, a loss of precision occurs when we do np.sum in this line:

https://github.com/TDAmeritrade/stumpy/blob/423c679fdfcdd547744d2b9a2aed7565124fd6e1/tests/naive.py#L1617

To be more precise, there are cases where the non-zero values of the 1D array np.minimum(D[i], Q) and np.minimum(D[j], Q) are the same BUT their sum is different! Although the error is negligible, it can result in a completely different outcome as the algorithm looks for the index with minimum value. The former points to index i and the latter points to index j. Hence, their corresponding snippets will be different.

---

For instance, for this input, one can follow the naive snippet to see that loss of precision. For the convenience, I am going to share those 1D arrays here to demonstrate the loss of precision in np.sum:

import numpy as np

non_zero_elems = """
0.024861743494183196  0.024861743494183196  0.024861743494183196
 0.021792470064670937  0.021792470064670937  0.018949939845733486
 0.018949939845733486  0.0015757541842979532 0.0015757541842979532
 0.0015757541842979532 0.0015757541842979532 0.0015757541842979532
 0.0015757541842979532 0.0015757541842979532 0.0015757541842979532
 0.04033997886552635   0.04033997886552635   0.04033997886552635
 0.04033997886552635  
 """
non_zero_elems = np.array(non_zero_elems.split()).astype(np.float64)
IDX1 = np.array([28, 29, 30, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
       53, 54])

IDX2 = np.array([ 8,  9, 10, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
       53, 54])

# `x` and `y` are two arrays with the same non-zero elements but at different positions
x = np.zeros(55, dtype=np.float64)
y = np.zeros(55, dtype=np.float64)

x[IDX1] = non_zero_elems.copy()
y[IDX2] = non_zero_elems.copy()

print(np.sum(x) == np.sum(y))  # False!! 
print(np.sum(x[x!=0]) == np.sum(y[y!=0]))  # True

---

Regarding the other case, I was able to produce a similar error with this input:

seed=2636
np.random.seed(seed)
T = np.random.uniform(-1000, 1000, [64]).astype(np.float64)

m=9
k=3
s=3

And those two 1D arrays with different sum are shared below:

a = np.array([
     0.005604770709214046, 0.005604770709214046, 0.005604770709214046,
     0.0352694999331615,   0.0352694999331615,   0.06102325948161525,
     0.06102325948161525,  0.04229368621667206,  0.004366737610721767,
     0.004366737610721767, 0.004366737610721767, 0.004366737610721767,
     0.004366737610721767, 0.004366737610721767, 0.004366737610721767,
     0.0658140591851339,   0.0658140591851339,  0.05934110672062023,
]).astype(np.float64)

IDX1 = np.array([16, 17, 18, 19, 20, 34, 35, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
       55])

IDX2 = np.array([ 7,  8,  9, 10, 11, 34, 35, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
       55])


# `x` and `y` are two arrays with the same non-zero elements but at different positions
x = np.zeros(56, dtype=np.float64)
x[IDX1] = a

y = np.zeros(56, dtype=np.float64)
y[IDX2] = a

x_sum = np.sum(x)
y_sum = np.sum(y)
assert x_sum == y_sum    # THIS FAILS!

[seanlaw]

Here's another one https://github.com/TDAmeritrade/stumpy/actions/runs/14322499583/job/40141974811

[NimaSarajpoor]

I updated my previous comment and expanded the response. In short, the issue is a loss of precision that occurs in np.sum in this line:

https://github.com/TDAmeritrade/stumpy/blob/423c679fdfcdd547744d2b9a2aed7565124fd6e1/tests/naive.py#L1617

A small loss of precision can result in drastic error in the snippet index.

A few notes:
(1) Solution? We can use math.fsum in the naive version which should give more accurate floating points. I tried it locally and the tests with those seeds are passing. My concern is that this creates an inconsistency in naive implementations unless we use math.fsum instead of np.sum in other naive implementations.

(2) Although some of the non-zero elements of the 1D array in the performant version can be different than their corresponding values in the naive version (due to loss of precision), that np.sum issue was not observed in the performant version for those two cases shared in my previous comment.

---

[Update]
Reported this in an issue on Numpy's GitHub repo

[seanlaw]

We can use math.fsum in the naive version which should give more accurate floating points.

I think using math.fsum in the naive version is reasonable (for all naive functions)