Investigate Failed Test in test_motifs.py

[seanlaw]

STUMPY_SEED=462698499 NUMBA_DISABLE_JIT=1 NUMBA_ENABLE_CUDASIM=1 pixi exec --spec black=26.5.1 --spec coverage=7.15.2 --spec dask=2026.7.1 --spec distributed=2026.7.1 --spec flake8=7.3.0 --spec isort=8.0.1 --spec numba=0.66.0 --spec numpy=2.4.6 --spec pandas=3.0.5 --spec polars=1.43.2 --spec pytest=9.1.1 --spec scipy=1.18.0 ./test.sh custom 1 tests/test_motifs.py

This test is failing here. It would be good to understand why this is happening and how best to fix it without compromising the test

=================================== FAILURES ===================================
_________________________________ test_motifs __________________________________

    def test_motifs():
        T = rng.RNG.rand(64)
        m = 3
    
        max_motifs = 3
        max_matches = 4
        max_distance = np.inf
        cutoff = np.inf
    
        # naive
        # `max_distance` and `cutoff` are hard-coded, and set to np.inf.
        ref_distances, ref_indices = naive_motifs(T, m, max_motifs, max_matches)
    
        # performant
        mp = naive.stump(T, m, row_wise=True)
        comp_distance, comp_indices = motifs(
            T,
            mp[:, 0].astype(np.float64),
            min_neighbors=1,
            max_distance=max_distance,
            cutoff=cutoff,
            max_matches=max_matches,
            max_motifs=max_motifs,
        )
    
        npt.assert_almost_equal(ref_indices, comp_indices)
>       npt.assert_almost_equal(ref_distances, comp_distance)
E       AssertionError: 
E       Arrays are not almost equal to 7 decimals
E       
E       Mismatched elements: 1 / 12 (8.33%)
E       Mismatch at index:
E        [0, 1]: 9.35810366834369e-07 (ACTUAL), 1.3102010687964719e-06 (DESIRED)
E       Max absolute difference among violations: 3.74390702e-07
E       Max relative difference among violations: 0.28575057
E        ACTUAL: array([[0.0000000e+00, 9.3581037e-07, 2.3419014e-01, 5.2195509e-01],
E              [0.0000000e+00, 6.9339867e-04, 6.7135066e-03, 2.2566633e-02],
E              [0.0000000e+00, 8.3248216e-04, 9.2865913e-02, 1.7967654e-01]])
E        DESIRED: array([[0.0000000e+00, 1.3102011e-06, 2.3419014e-01, 5.2195509e-01],
E              [0.0000000e+00, 6.9339868e-04, 6.7135066e-03, 2.2566633e-02],
E              [0.0000000e+00, 8.3248216e-04, 9.2865913e-02, 1.7967654e-01]])

tests/test_motifs.py:619: AssertionError
========================= 1 failed, 21 passed in 0.23s =========================
Error: Test execution encountered exit code 1
Error: Process completed with exit code 1.

[jackwalkerlabs]

I reproduced the failure with the reported seed and dependency versions, both with JIT disabled and enabled. I have not changed any code yet.

What I found:

• The motif selection is correct: all 12 computed indices match the naive reference.
• The mismatch is limited to the distance between the length-3 subsequences at indices 32 and 43. They are almost perfectly correlated.
• A 100-digit calculation from the underlying float64 inputs gives a distance of 9.358103667676396e-7.
• The naive calculation returns 9.35810366834369e-7 (about 6.7e-17 absolute error).
• MASS returns 1.3102010687964719e-6 with JIT disabled and 1.2299552403865104e-6 with JIT enabled.
• Here, 1 - rho is only about 1.46e-13. The raw sliding-dot-product/correlation formula therefore amplifies normal float64 rounding through cancellation. This is not specific to Numba fastmath.

My proposed minimal approach is:

1. Keep the randomized input and the existing exact motif-index comparison.
2. Change only the distance assertion to use decimal=config.STUMPY_TEST_PRECISION, as other naive-vs-performant comparisons do throughout the suite.
3. Re-run this exact seed with JIT disabled and enabled, then run the complete formatting, lint, and test checks.

This preserves the test's behavioral signal: incorrect motif selection still fails exactly, and material distance errors still fail at the project's configured five-decimal precision. It avoids fixing the seed or adding a one-off tolerance.

Would you be comfortable with this approach?

[seanlaw]

Here, 1 - rho is only about 1.46e-13. The raw sliding-dot-product/correlation formula therefore amplifies normal float64 rounding through cancellation.

Is there a better mathematical formulation that can reduce/minimize this amplification, improve the precision of the computation, without significantly increasing the overall computational time? I think that is the solution that we're looking for

[jackwalkerlabs]

Thanks — I investigated a numerically stable hybrid formulation locally.

The approach is:

1. Keep the existing FFT/correlation calculation for ordinary candidates.
2. When the computed squared distance is positive and below 2 * m * sqrt(float64 epsilon) — equivalently, when 1 - rho < sqrt(epsilon) — recompute only that candidate as the sum of squared pointwise differences between the two z-normalized subsequences.
3. Fold the threshold check into the existing distance-profile loop, so there is no second full scan. The cost remains O(n) after the FFT plus O(c * m), where c is the number of extremely close candidates.

This avoids subtracting a correlation from one in precisely the range where cancellation dominates. For the failing subsequences:

• 100-digit reference: 9.358103667676396e-7
• current MASS, JIT disabled: 1.3102010687964719e-6
• stable calculation, JIT enabled: 9.35810366765695e-7 (about 1.9e-18 absolute error)

I added a deterministic regression using those two subsequences. The regression and the original STUMPY_SEED=462698499 motif test pass with JIT disabled and enabled. The focused core, precision, and motif suites pass in both modes (271 passed, 2 skipped each), and Black, isort, Flake8, docstring, and fastmath-policy checks pass.

A compiled-kernel microbenchmark over 99,901 candidates (m=100, one threshold candidate) measured about 114 microseconds for the existing profile calculation and 156 microseconds for the hybrid loop; end-to-end mass timing was within run-to-run noise because the FFT and preprocessing dominate.

Would you be comfortable with this formulation and threshold? I will wait for confirmation before finalizing the branch and opening a pull request.

[jackwalkerlabs]

One correction to the condition above: the threshold should include computed zero distances, not only positive ones. A stress test found a non-identical near-affine pair (Q = [1, 2, 4], T = [-4, -0.99999999, 5]) whose correlation formula rounds the distance to exactly zero, while the direct z-normalized distance is 2.1428570903602775e-9. Recomputing all values below the threshold, including zero, recovers that distance to float64 precision. I added this as a second local regression. The asymptotic cost is unchanged; it also recomputes the usual self-match once in O(m).

[seanlaw]

Would you be comfortable with this formulation and threshold? I will wait for confirmation before finalizing the branch and opening a pull request.

Frankly, I'm finding it a bit difficult to follow your logic. At this point, it's probably better to see some code so please go ahead and open up a PR.