BUG: Fix np.inf + np.nan sum issue on groupby mean

[parthi-siva]

• closes BUG: Groupby mean differs depending on where np.inf value was introduced #50367
• Tests added and passed if fixing a bug or adding a new feature
• All code checks passed.

[parthi-siva]

@jbrockmendel @rhshadrach please review. Please let me know if it 's correct. Then I'll write test cases.

[rhshadrach]

This does not look correct to me; in particular, np.inf + np.nan should be np.nan and if I'm understanding this code it would come out to np.inf. In addition, this doesn't look like it handles -np.inf.

Instead, I would guess that when val is +/- np.inf, then compensation should be set to 0.

This additional check would make the current implementation more complex in the inner for loop, and I'm wondering if the additional numerical accuracy is worth the perf hit. But maybe branch prediction means the perf hit would be minimal in what I would expect to be the typical case (no inf).

[parthi-siva]

Thanks @rhshadrach for the review.

We have to make compensation 0 not only when val is +/- np.inf also when sumx[lab, j] is +/- np.inf. Then only the results are consistent

something like

if val == np.inf or sumx[lab, j] == np.inf:
    compensation[lab, j] = 0.0
else:
    compensation[lab, j] = t - sumx[lab, j] - y

And yeah we are introducing additional check which will impact the performance. Also you are right about fact that sacrificing performance for numerical accuracy is worth or not.

What should we do? Is it okay that we live with this bug ? or fix it ?

[rhshadrach]

Can you benchmark the proposed fix?

[parthi-siva]

sure @rhshadrach.. Can you help me with that..? direct me to any resource or docs on standardized way to do benchmark ?

%timeit will be sufficient right ?

[rhshadrach]

Yes - %timeit is a good way to approach. Generate some test data (e.g. np.random.randint and np.random.rand) of sufficient size so that the groupby operation takes at least a decent amount of time - maybe 80ms or more. Then post the timings when you run it on main vs this PR.

[parthi-siva]

instead of checking val for inf or -inf we can check for compensation is nan or not:

 1078  if not isna_entry:
                nobs[lab, j] += 1
                y = val - compensation[lab, j]
                t = sumx[lab, j] + y
                compensation[lab, j] = t - sumx[lab, j] - y
                if compensation[lab, j] != compensation[lab, j] :
                    compensation[lab, j] = 0.
                sumx[lab, j] = t

Benchmark against Main Branch

In [5]: n = 5000000
   ...: ...: df = pd.DataFrame(
   ...: ...:     {
   ...: ...:         "z": [rand.randint(0,1) for i in range(n)],
   ...: ...:         "x1": [rand.randint(0,n) for i in range(n)],
   ...: ...:         "x2": [rand.randint(0,n) for i in range(n)],
   ...: ...:      }
   ...: ...: )

In [6]: %timeit df.groupby("z").mean()
79.8 ms ± 306 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [7]: %timeit df.groupby("z").mean()
79.9 ms ± 789 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [8]: %timeit df.groupby("z").mean()
79.8 ms ± 614 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [9]: %timeit df.groupby("z").mean()
80.5 ms ± 254 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [10]: %timeit df.groupby("z").mean()
80.1 ms ± 650 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [11]: %timeit df.groupby("z").mean()
79.2 ms ± 686 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [12]: %timeit df.groupby("z").mean()
79.4 ms ± 674 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [13]: %timeit df.groupby("z").mean()
79.4 ms ± 743 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [14]: %timeit df.groupby("z").mean()
83.2 ms ± 3.34 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Benchmark against branch which contains proposed fix

In [9]: n = 5000000
   ...: ...: df = pd.DataFrame(
   ...: ...:     {
   ...: ...:         "z": [rand.randint(0,1) for i in range(n)],
   ...: ...:         "x1": [rand.randint(0,n) for i in range(n)],
   ...: ...:         "x2": [rand.randint(0,n) for i in range(n)],
   ...: ...:      }
   ...: ...: )

In [10]: %timeit df.groupby("z").mean()
81.5 ms ± 378 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [11]: %timeit df.groupby("z").mean()
81.8 ms ± 891 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [12]: %timeit df.groupby("z").mean()
81.5 ms ± 538 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [13]: %timeit df.groupby("z").mean()
82 ms ± 754 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [14]: %timeit df.groupby("z").mean()
83.8 ms ± 1.49 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [15]: %timeit df.groupby("z").mean()
81.3 ms ± 544 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [16]: %timeit df.groupby("z").mean()
81.7 ms ± 768 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [17]: %timeit df.groupby("z").mean()
81.6 ms ± 404 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [18]: %timeit df.groupby("z").mean()
81.7 ms ± 500 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

The overhead is less than 2 ms @rhshadrach

[rhshadrach]

Looks good, cc @mroeschke for a 2nd eye

[parthi-siva]

Done

[jbrockmendel]

this is to check for nan? can you use an explicit check? If not, please add a comment about what you're doing and why

[parthi-siva]

Hi @jbrockmendel, I tried do it explicitly but it is not working. I tried to use utils.is_nan since it is no gil it is not compiling. So I just used this way to determine Nan (That's what we are doing in utils.is_nan as well.)
I will add appropriate comment in the code

[jbrockmendel]

nitpick: please define the expected on a previous line

[parthi-siva]

Done

[rhshadrach]

@parthi-siva - looks like the 32-bit is failing here.

https://github.com/pandas-dev/pandas/actions/runs/4911671351/jobs/8769921340?pr=52964#step:4:1956

[parthi-siva]

@rhshadrach Linux 32-Bit test failing issue fixed.

[rhshadrach]

lgtm, cc @jbrockmendel - good here?

[mroeschke]

Thanks @parthi-siva