TST: stats.cosine: modify test to silence failure - #18382

@ilayn That failure in test_cosine_logpdf_endpoints is because

x = [-np.pi, np.pi]
c = np.cos(x)
np.testing.assert_equal(c, -1)

fails. I think that sort of change would cause problems in other libraries, too. Should we report to NumPy?

I think I'm looking at the wrong place maybe but I am reading

E    x: array([-37.88153, -37.88153])
E    y: array([-inf, -inf])

as the mismatch and not the element. Is that correct?

In this PR, the failure is in the revised line 8764:

_________________________ test_cosine_logpdf_endpoints _________________________
[gw1] linux -- Python 3.11.3 /opt/hostedtoolcache/Python/3.11.3/x64/bin/python
scipy/stats/tests/test_distributions.py:8764: in test_cosine_logpdf_endpoints
    assert_equal(c, -1)
        c          = array([-1., -1.])
        x          = [-3.141592653589793, 3.141592653589793]
/opt/hostedtoolcache/Python/3.11.3/x64/lib/python3.11/contextlib.py:81: in inner
    return func(*args, **kwds)
E   AssertionError: 
E   Arrays are not equal
E   
E   Mismatched elements: 2 / 2 (100%)
E   Max absolute difference: 2.22044605e-16
E   Max relative difference: 2.22044605e-16
E    x: array([-1., -1.])
E    y: array(-1)
        args       = (<built-in function eq>, array([-1., -1.]), -1)
        func       = <function assert_array_compare at 0x7fd0d2944220>
        kwds       = {'err_msg': '', 'header': 'Arrays are not equal', 'strict': False, 'verbose': True}
        self       = <contextlib._GeneratorContextManager object at 0x7fd0d2943190>

revealing that the failure you mentioned is actually due to a change in the value returned by np.cos(-np.pi)/np.cos(np.pi).

Ah OK. this revealed a different issue. I think

E    x: array([-1., -1.])
E    y: array(-1)

indeed should invalidate the assertion. We should probably test the size too anyways.

I'm confused. It sounds like you're saying that the assertion fails because they are of different shapes, but assert_equal allows comparison of arrays to scalars.

It says that it's failing because the results of np.cos(np.pi) and np.cos(-np.pi), despite being displayed as -1, differ by 2.22044605e-16 from the integer -1. I modified the existing test to reveal that this is the underlying cause of the failure you mentioned.

My question was whether you agree this is something to report to NumPy. My question was whether you agree this is something to report to NumPy. It looks like it could be a bug in NumPy main, so if that's the case and they fix it, our failure would be resolved.

Oh maybe linked to numpy/numpy#23399 ??? I think Astropy could be seeing the same thing astropy/astropy#14705

Ah then I got what you mean now, sorry for my confusion. That value is just floating number jitter. If you look at the eps value,

np.spacing(1.)
Out[3]: 2.220446049250313e-16

which is the smallest number you can distinguish a float from another differing that amount. That's why typically you don't do equal with floats. Integer vs. float comparisons are always prone to that. If you want to tighten the equality assert_allclose with atol very small and rtol = 0. can be used. Or any other meaningful atol, rtol pair.

Why suddenly we started seeing is probably due to what @tupui mentioned, that they fiddled with some internals of trig funcs.

Re: astropy -- it might have broken some stuff in astropy.units too (for the coordinates failures) but I would need @mhvk to confirm. On the modeling side, does look like floating point differences but enough to cause some bounding box comparisons to fail which needs @WilliamJamieson to investigate.

That value is just floating number jitter

I know. However, with something as fundamental as np.cos(np.pi), it is pretty disruptive to change it.

Well I would argue that the equal shoudn't have been used for any float comparison anyways. That would fix this issue at hand but losing this convenience is a different discussion.

Following that, assert_equal should not work if the dtype is float no?

For certain float values assert_equal works OK. If you are working with decimal float values etc. So I don't think the conclusion follows from that but if your values are coming out of a float operation from a black box then you shouldn't rely on the architecture for generating pure integers.

What I mean is this type of stuff should work

x = 5.
assert_equal(x + 1, 6)
assert_equal(x+0.3, 5.3)
assert_equal(0.1+0.1+0.1, 0.3)  # this famous FP example expectedly breaks 

Think of equal as a bit-level comparison. If ULP guarantees become loose, then you can't pull the comparison off with ==. So here the our comparison quantifier is not suitable to start with, but this does not change the fact that the NumPy PR did a substantial regression for potentially high-precision computations.

@mdhaber What do you think? Should we move forward with the assert_allclose conversion?

You mean np.isclose(c, -1, atol=...) within the code itself?

That would stop the test from complaining, but it's unclear whether that is the best fix. There is no need for the result of cosine._logpdf to jump to negative infinity unless c is exactly -1 because the calculation can still be performed in the primary branch of _lazywhere, without using the fillvalue. In that case, the test_cosine_logpdf_endpoints could be modified to accept any value less than or equal to some threshold.

Then again, I suppose I would prefer for the cosine distribution to have exactly zero probability density when evaluated at np.pi , so yeah.

With cf373ac in #18402 I was able to get down to one failure:

=========================== short test summary info ============================
FAILED scipy/stats/tests/test_distributions.py::test_cosine_logpdf_endpoints - AssertionError: 
Arrays are not equal

x and y -inf location mismatch:
 x: array([-37.88153, -37.88153])
 y: array([-inf, -inf])
= 1 failed, 39630 passed, 2194 skipped, 140 xfailed, 10 xpassed in 425.35s (0:07:05) =

AIUI, this is similarly due to numpy/numpy#23399. Seems there's not much but to patch our tests to have less sharp expectations.

As the original author of that overly optimistic test, I'm happy to fix the issue. But the fix that I have in mind is more radical than just tweaking the test.

The true 64 bit floating point value for cosine.logpdf(np.pi) is actually not -inf, but -75.80844227238724, because np.pi is not equal to $\pi$. I have the cosine distribution implemented in mpsci:

In [58]: import numpy as np

In [59]: from mpmath import mp

In [60]: from mpsci.distributions import cosine

In [61]: mp.dps = 200

In [62]: float(cosine.logpdf(np.pi))
Out[62]: -75.80844227238724

The fix that I have in mind will update both the pdf and logpdf methods to provide much more accurate values in the neighborhoods of $x = \pi$ and $x = -\pi$. Note that because of the catastrophic loss of precision in the expression $1 + \cos(x)$ for $x$ near $\pm \pi$, the pdf method of stats.cosine is not very good near the endpoints:

In [77]: from scipy import stats

In [78]: x = 3.141592

In [79]: stats.cosine.pdf(x)
Out[79]: 3.4014265343429277e-14

In [80]: float(cosine.pdf(x))  # Reference value
Out[80]: 3.399387386564633e-14

What I have in mind is to implement the pdf and/or logpdf methods as a ufunc, similar to what I did for the cdf method in #10170.

I think we need something to get CI green in the meantime. (Also, personally, I would not be interested in improving the accuracy of the cosine distribution in this range of inputs at this time.)

@h-vetinari you suggested giving the test "less sharp expectations". This does that by ensuring that the resulting logpdf is below a certain threshold, a few ticks above what we might expect.

I think we need something to get CI green in the meantime

Agreed. This is a one-line fix, and a follow-up issue can be opened for a more desirable change in the longer term. Let's merge this, and follow up on it at a more leisurely pace later?

In it goes. I'm looking at the related failures in signal.windows in gh-18406.