Wow even with float64 there are many errors ...
Small mistake in the check.

Sorry for opening a PR and running away :)

LGTM

Indeed that could be problematic. LGTM

LGTM

Any ideas of pathological cases to test with quickly?

No idea, except comparing roc_auc_score(y_true, y_pred, sample_weight=sample_weight_32) and  roc_auc_score(y_true, y_pred, sample_weight=sample_weight_64), yet it does not test the bug we want to fix, but the correction we propose, so I am not fond of this solution.

I only think it's an effective test for very large vectors, hence
impractical.

On 6 September 2016 at 20:05, Tom Dupré la Tour notifications@github.com
wrote:

No idea, except comparing roc_auc_score(y_true, y_pred,
sample_weight=sample_weight_32) and roc_auc_score(y_true, y_pred,
sample_weight=sample_weight_64), yet it does not test the bug we want to
fix, but the correction we propose, so I am not fond of this solution.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#7331 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAEz6zdyrhBjU7xPTazZ20BVeBDqeRrgks5qnTrggaJpZM4JzY5A
.

how long does a test with a large-enough vector run?

lgtm

import sys
import time
import numpy as np
import pandas as pd
n_trials = 50
all_results = []
for dtype in [np.float32, np.float64]:
    for i in range(3, 8):
        n = 1 * 10 ** i
        absdiff, reldiff = [], []
        s = time.time()
        for j in range(n_trials):
            x = np.random.rand(n).astype(dtype)
            a = np.cumsum(x)[-1]
            b = np.sum(x)
            absdiff.append(np.abs(a - b))
            reldiff.append(absdiff[-1] / b)
        all_results.append((dtype.__name__, n, (time.time() - s) / len(results), np.log10(np.mean(absdiff)), np.log10(np.mean(reldiff))))

pd.DataFrame(all_results, columns=['dtype', 'n', 'time', 'log abs diff', 'log rel diff'])

np.allclose has default rtol=1e-5, atol=1e-8

Two questions:

• can we build a non-regression case where the old code would have been broken? possibly: all the absolute diffs above for float32 exceed the atol; for 1e6 samples we get relative diff just greater than allclose's rtol in 0.02s, but this would be slower embedded in the rest of roc_auc_score.
• can we get test coverage for the error message here? probably not, without having rtol as a parameter to stable_cumsum. This might be worth doing.

So I'll try build a non-regression test with 1e6 samples.

A regression test for roc_auc_score turns out to be too slow.

I've added a direct test of stable_cumsum using configurable atol, rtol

I suppose the remaining question is: are there other places in the codebase where cumsum is happening and might be over float32 data or very very large arrays?

One of the Travis error will go away if you rebase on master but there is one on Python 2.7 that seems genuine:

======================================================================
FAIL: sklearn.utils.tests.test_extmath.test_stable_cumsum
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/miniconda/envs/testenv/lib/python2.6/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/travis/sklearn_build_oldest/scikit-learn/sklearn/utils/tests/test_extmath.py", line 656, in test_stable_cumsum
    stable_cumsum, r, rtol=0, atol=0)
  File "/home/travis/sklearn_build_oldest/scikit-learn/sklearn/utils/testing.py", line 438, in assert_raise_message
    (names, function.__name__))
AssertionError: RuntimeError not raised by stable_cumsum

One of the Travis error will go away if you rebase on master but there is one on Python 2.7 that seems genuine:

Yeah, I was wondering if there'd be some platform that didn't fail my test....

I assume that means numpy used to have an unstable implementation of sum, rather than it used to have a stable implementation of cumsum

I wonder: Should I only run the test on recent numpy, or just remove it?

Hm... testing the error message seems slightly odd since it relies on numpy being broken. I guess having a test with a known correct result that wasn't correct before might be nicer, but not offer the error coverage?

Like

x = np.empty(1e7)
x.fill(1e-7)
np.isclose(stable_cumsum(x)[-1], 1, rtol=0, atol=1e-12)

It might just be my bedtime, but I'm not sure I get what test you're suggesting.

np.isclose(cumsum(x)[-1], 1, rtol=0, atol=1e-12)
fails for np.cumsum currently, but should work for stable_cumsum

No, I don't think it will work for stable_cumsum in old numpy. I might be wrong though.

ah... huh...

So just drop the test, I suppose.

And merge the PR?

fine with me

I believe the improved stability of np.sum was done in numpy 1.9: http://docs.scipy.org/doc/numpy/release.html#better-numerical-stability-for-sum-in-some-cases and I also quickly checked that in 1.8 cumsum and sum were giving the same wrong result in one of your snippet but 1.9 was fine. Maybe skip the test for numpy < 1.9?

could please run a check for np1.9, @lesteve?

could please run a check for np1.9, @lesteve?

Not sure what you exactly mean, but here is what I tried:

import numpy as np

np.random.seed(42)
n_samples = 4 * 10 ** 7
y = np.random.randint(2, size=n_samples)
prediction = np.random.normal(size=n_samples) + y * 0.01
trivial_weight = np.ones(n_samples)

print(np.cumsum(trivial_weight.astype('float32'))[-1])
print(np.sum(trivial_weight.astype('float32')))

Output for numpy 1.8:

1.67772e+07
1.67772e+07

Output for numpy 1.9:

1.67772e+07
4e+07

I also made sure that the test was run for numpy 1.9. LGTM will wait for AppVeyor and then merge.

Merging, thanks!