Precision Unit Tests Failing

[seanlaw]

@NimaSarajpoor Is there any chance that you can take a look at these failing precision-related unit tests:

https://github.com/stumpy-dev/automate/actions/runs/12051839571/job/33655731907

The numba release candidate (rc) has been updated/changed to Numba 0.61.0rc1 and llvmlite 0.44.0rc1.

Note that the last numba release candidate that had passing tests were for:

numba-0.60.0rc1 was llvmlite-0.43.0rc1 (see here)

[seanlaw]

Tagging @esc and @kc611 for awareness. This is the first time an RC has caused a failure so maybe you can suggest a path forward for assessing this?

[esc]

@seanlaw thank you for flagging this! I've added a link to the thread tracking issues with the RC on Numba's discourse:

https://numba.discourse.group/t/ann-numba-0-61-0rc1-and-llvmlite-0-44-0rc1-are-out/2816/3

which is a good start and makes sure that the Numba dev community is aware of the issue.

Ideally, the next step would be to convert the failing unit-test into a minimal reproducer and then open an issue on the Numba issue tracker with the reproducer. I would suggest to try to narrow things down in this thread here until a reproducer can be established and to then port that over to the Numba issue tracker.

Also pinging @stuartarchibald and @sklam for visibility.

[esc]

The failure is this:

============================= test session starts ==============================
platform linux -- Python 3.11.10, pytest-8.3.3, pluggy-1.5.0
rootdir: /home/runner/work/automate/automate
configfile: pytest.ini
collected 5 items

tests/test_precision.py ...Fs                                            [100%]

=================================== FAILURES ===================================
________________________________ test_snippets _________________________________

    def test_snippets():
        # This test function raises an error if there is a considerable loss of precision
        # that violates the symmetry property of a distance measure.
        m = 10
        k = 3
        s = 3
        seed = [332](https://github.com/stumpy-dev/automate/actions/runs/12051839571/job/33655731907#step:10:333)
        np.random.seed(seed)
        T = np.random.uniform(-1000.0, 1000.0, [64])
    
        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,
        ) = stumpy.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_precision.py:149: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/contextlib.py:81: in inner
    return func(*args, **kwds)
/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/contextlib.py:81: in inner
    return func(*args, **kwds)
/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/numpy/_utils/__init__.py:85: in wrapper
    return fun(*args, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

args = (<function assert_array_almost_equal.<locals>.compare at 0x7f86b45da020>, array([[  61.66844566,  600.71676057,   93.4....299013  ,
         703.22299054, -807.06416843,  962.8838992 ,   38.81042428,
         -83.86261799,   18.55095966]]))
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: 10 / 30 (33.3%)
E           Max absolute difference among violations: 1546.55831117
E           Max relative difference among violations: 20.34471955
E            ACTUAL: array([[  61.66845,  600.71676,   93.43324,  380.43982, -[341](https://github.com/stumpy-dev/automate/actions/runs/12051839571/job/33655731907#step:10:342).5808 ,
E                    455.06518,  -92.4037 , -321.6817 ,  731.41598, -694.49969],
E                  [ 743.18048, -887.71494, -286.9[345](https://github.com/stumpy-dev/automate/actions/runs/12051839571/job/33655731907#step:10:346)7, -562.81888,  306.91077,...
E            DESIRED: array([[  61.66845,  600.71676,   93.43324,  380.43982, -341.5808 ,
E                    455.06518,  -92.4037 , -321.6817 ,  731.41598, -694.49969],
E                  [ 743.18048, -887.71494, -286.93457, -562.81888,  306.91077,...

/opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/contextlib.py:81: AssertionError
=========================== short test summary info ============================
SKIPPED [1] tests/test_precision.py:169: Skipping Tests No GPUs Available
==================== 1 failed, 3 passed, 1 skipped in 7.01s ====================
Error: Test execution encountered exit code 1

[NimaSarajpoor]

@seanlaw @esc
I will work on it to provide a minimal reproducer.

[seanlaw]

Thank you, @NimaSarajpoor! I really appreciate it. 🙏 Please let me know if there's anything I can do to help.

In case it matters, I actually re-ran that failed test 2 times (i.e., 3 times in total) and they all failed at the same spot and with the same array mismatch error:

Attempt #1
Attempt #2
Attempt #3

[kc611]

The failing tests are probably a symptom of changes in fastmath across LLVM versions. This can be demonstrated by turning off all fastmath operations within the library.

The exact combination / location of flags causing the issue is still unknown, however the fact still remains that the precision tests are highly sensitive to fastmath changes.

[NimaSarajpoor]

@seanlaw
I was able to produce the failure in Colab using numba==0.61.0rc1. After some investigation, I was able to narrow it down to the precision issue that comes from core._mass. IMO, this is not MWE, and I need to dig deeper.

inputs

s = 3  

seed = 332
np.random.seed(seed)
T = np.random.uniform(-1000.0, 1000.0, [64])
T = np.pad(T, pad_width=(0, 6), mode="constant", constant_values=np.nan)

M_T, Σ_T = stumpy.core.compute_mean_std(T, s)

T_subseq_isconstant = np.zeros(len(T) -s + 1, dtype=bool)
IDX = np.array([ 8, 17, 26, 51])
T_subseq_isconstant[IDX] = True

outputs

QT_T2 = core._sliding_dot_product(T[2 : 2 + s], T)
dist_prof_T2 = core._mass(
    T[2 : 2 + s],
    T,
    QT_T2,
    M_T[2],
    Σ_T[2],
    M_T,
    Σ_T,
    T_subseq_isconstant[2],
    T_subseq_isconstant,
)


QT_T10 = core._sliding_dot_product(T[10 : 10 + s], T)
dist_prof_T10 = core._mass(
    T[10 : 10 + s],
    T,
    QT_T10,
    M_T[10],
    Σ_T[10],
    M_T,
    Σ_T,
    T_subseq_isconstant[10],
    T_subseq_isconstant,
)

we expect to see dist_prof_T2[10] == dist_prof_T10[2] since they are just the distance between the subsequence T[2:2+s] and T[10:10+s]. However, they are not equal.

---

[Update 1] Going deeper
So, I checked the core._mass and went deeper to see what's happening down the road...

s = 3
QT = 1453258.774722079
μ_Q = -55.35781464461123
σ_Q = 650.912209452633
μ_T = -264.9540642227162
σ_T = 722.0717285148526
isconstant_Q = False
isconstant_T = False


ref = core._calculate_squared_distance(
    s, 
    QT, 
    μ_Q, 
    σ_Q, 
    μ_T, 
    σ_T, 
    isconstant_Q, 
    isconstant_T,
)


comp = core._calculate_squared_distance(
    s, 
    QT,  
    μ_T, 
    σ_T,
    μ_Q, 
    σ_Q, 
    isconstant_T,
    isconstant_Q
)

If we look at the function _calculate_squared_distance , we can see that ref and comp are expected to be the same value. However, we can see ref != comp

---

[Update 2] Going deeper again!

So, I tried to check the function.

from numba import njit
@njit(
    # "f8(i8, f8, f8, f8, f8, f8)",
    fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}
)
def func(
    m, QT, μ_Q, σ_Q, M_T, Σ_T, Q_subseq_isconstant, T_subseq_isconstant
):
    if np.isinf(M_T) or np.isinf(μ_Q):
        D_squared = np.inf
    elif Q_subseq_isconstant and T_subseq_isconstant:
        D_squared = 0
    elif Q_subseq_isconstant or T_subseq_isconstant:
        D_squared = m
    else:
        denom = (σ_Q * Σ_T) * m
        denom = max(denom, 1e-14)

        ρ = (QT - (μ_Q * M_T) * m) / denom
        ρ = min(ρ, 1.0)

        D_squared = np.abs(2 * m * (1.0 - ρ))

    return D_squared

For the inputs provided in my previous update above, only else part is executed. Now, if I revise the function to just keep the else part, then:

from numba import njit
@njit(
    fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}
)
def func(
    m, QT, μ_Q, σ_Q, M_T, Σ_T, Q_subseq_isconstant, T_subseq_isconstant
):   
    denom = (σ_Q * Σ_T) * m
    denom = max(denom, 1e-14)

    ρ = (QT - (μ_Q * M_T) * m) / denom
    ρ = min(ρ, 1.0)

    D_squared = np.abs(2 * m * (1.0 - ρ))

    return D_squared


s = 3
QT = 1453258.774722079
μ_Q = -55.35781464461123
σ_Q = 650.912209452633
μ_T = -264.9540642227162
σ_T = 722.0717285148526
isconstant_Q = False
isconstant_T = False


ref = func(
    s, 
    QT, 
    μ_Q, 
    σ_Q, 
    μ_T, 
    σ_T, 
    isconstant_Q, 
    isconstant_T,
)


comp = func(
    s, 
    QT,  
    μ_T, 
    σ_T,
    μ_Q, 
    σ_Q, 
    isconstant_T,
    isconstant_Q
)

assert ref == comp

and it is all good now!

The is strange! So, the function with those if/elif parts results in assertion failure BUT the second function in which I just kept the else part, does not result in assertion failure!!

[NimaSarajpoor]

Minimal Working Example (MWE)

# numba version: 0.61.0rc1

from numba import njit
@njit(
    fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}
)
def func(
    a, b
):  
    m = 3
    out = (a * b) * m
        
    return out

##################
a = 650.912209452633
b = 722.0717285148526

ref = func(a, b)


comp = func(b, a)

assert ref == comp

[Potential Solution]
Removing reassoc seems to fix the issue here.

[Note]
I do not know why previously, as mentioned in my previous comment, the issue was resolved by just removing those if / elif parts from the function core._calculate_squared_distance and keeping the else part only.

[seanlaw]

@NimaSarajpoor and I had a discussion over here but repeating here for transparency:

@NimaSarajpoor As I think about it more, I am starting to question whether turning off reassoc is a good idea because we're getting divergance in how we use fastmath. Sometimes we have:

fastmath=True
fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}

and now we are adding a third variant fastmath={"nsz", "arcp", "contract", "afn"} (i.e., removing "reassoc:). This is bothersome to maintain and smells fishy IMHO. I feel like our North Star should probably be to match our fast implementation with our naive implementations (i.e., what would a user reasonably expect when they try to calculate things themselves using numpy or base Python)? What do you think?

@kc611 Would you happen to have any guidance as to what to do under these circumstances? Consistency is key and we want to avoid/minimize having one-off exceptions like this as it will make it hard to maintain long term.

[seanlaw]

@kc611 and @esc I just wanted to provide a quick comment that @NimaSarajpoor is actively crafting up a solution for this issue in TDAmeritrade/stumpy#1048

[esc]

@seanlaw any updates here? We are about to tag the finals today.

[seanlaw]

@seanlaw any updates here? We are about to tag the finals today.

@esc Yes! @NimaSarajpoor is wrapping up a PR that will work around this issue.

[seanlaw]

Closed by TDAmeritrade/stumpy#1048

@NimaSarajpoor Thanks again for all of your persistence and hard work!