diff --git a/README.md b/README.md index 000ed9d..5b1e77e 100644 --- a/README.md +++ b/README.md @@ -208,26 +208,26 @@ All implementations share the calling syntax and produce the same results up to ## Performance -Best of five timed runs, in milliseconds, for 500,000 values in 1,000 groups (lower is better): +Best of three passes of five timed runs each, in milliseconds, for 500,000 values in 1,000 groups (lower is better): -| function | numpy | numba | numbagg | pandas | -| -------- | -----: | -----: | ------: | -----: | -| `sum` | 0.92 | 0.43 | 0.49 | 15.89 | -| `mean` | 1.56 | 0.51 | 0.56 | 15.99 | -| `max` | 1.71 | 0.64 | 0.88 | 15.15 | -| `std` | 3.09 | 0.67 | 0.68 | 16.82 | -| `nansum` | 4.17 | 0.44 | 1.76 | 22.26 | -| `median` | 48.41 | 11.49 | --- | 27.70 | -| `sort` | 135.15 | 25.99 | --- | --- | -| `cumsum` | 44.84 | 0.58 | --- | 7.32 | -| custom callable | 144.10 | 42.97 | --- | 108.59 | +| function | numpy | numba | numbagg | pandas | +| --------------- | -----: | -----: | ------: | -----: | +| `sum` | 0.95 | 0.44 | 0.47 | 15.96 | +| `mean` | 1.52 | 0.51 | 0.55 | 16.98 | +| `max` | 1.68 | 0.63 | 0.85 | 16.26 | +| `std` | 3.09 | 0.66 | 0.69 | 17.02 | +| `nansum` | 4.24 | 0.45 | 1.78 | 21.71 | +| `median` | 48.77 | 11.52 | --- | 28.57 | +| `sort` | 136.16 | 25.96 | --- | --- | +| `cumsum` | 45.88 | 0.58 | --- | 7.57 | +| custom callable | 142.54 | 43.16 | --- | 109.89 | -In short: the NumPy implementation is about 5–24× faster than pandas on common reductions, and Numba adds another 2–10× on top — and about 77× for `cumsum`. numbagg is par with us on the light scans; we lead on the nan-heavy accumulations (up to 4×), and numbagg has no `median`, cumulatives, `trapezoid` or `sort`. `median`, `sort` and custom callables are the expensive cases everywhere. Absolute numbers depend on your machine, so treat the ratios as the takeaway. +In short: the NumPy implementation is about 5–18× faster than pandas on common reductions, and Numba adds another 2–10× on top — and about 79× for `cumsum`. numbagg is par with us on the light scans; we lead on the nan-heavy accumulations (up to 4×), and numbagg has no `median`, cumulatives, `trapezoid` or `sort`. `median`, `sort` and custom callables are the expensive cases everywhere. Absolute numbers depend on your machine, so treat the ratios as the takeaway. -The complete table and the benchmark setup are in [docs/benchmarks.md](https://github.com/ml31415/numpy-groupies/blob/master/docs/benchmarks.md). To run it yourself, from the repository root (add `--pandas --numbagg` for those columns, `--purepy` for the pure-Python one): +The complete table and the benchmark setup are in [docs/benchmarks.md](https://github.com/ml31415/numpy-groupies/blob/master/docs/benchmarks.md). To run it yourself, from the repository root — `--help` explains which columns to ask for: ```sh -python -m numpy_groupies.benchmarks.generic --pandas --numbagg +python -m numpy_groupies.benchmarks --numbagg --pandas ``` The `numpy` and `numba` columns above are not the whole story: NumPy 1.25 brought major [speed improvements to ufuncs](https://numpy.org/doc/stable/release/1.25.0-notes.html), which narrowed the gap between the NumPy and Numba implementations considerably. The authors hope that `ufunc.at` or an equivalent in NumPy or SciPy will eventually become fast enough to make this package redundant. diff --git a/docs/benchmarks.md b/docs/benchmarks.md index 18db7b4..15b65ad 100644 --- a/docs/benchmarks.md +++ b/docs/benchmarks.md @@ -7,70 +7,68 @@ Complete benchmark results for `aggregate`. A short summary is in the [README](. - 500,000 indices, uniformly drawn from `[0, 1000)`. - The values of `a` are uniform on `[0, 1)`, and everything above `0.8` is then set to `0`, so that there are falsy values for the boolean functions. - For the `nan…` functions, of the remaining values those below `0.2` are set to `nan` (about 20%), leaving the rest on the interval `[0.2, 0.8]`. -- Times in **milliseconds**, taking the minimum over 5 timed runs after one untimed warm-up call (`repeat=5` in `numpy_groupies/benchmarks/generic.py`). Lower is better. +- Times in **milliseconds**, the fastest of five timed executions after one untimed warm-up call, over three complete passes of the table (`--repeat=5 --runs=3` of `numpy_groupies.benchmarks`). Lower is better. - Software: Python 3.14.7, NumPy 2.5.4, Numba 0.68.0, pandas 3.0.6, numbagg 0.9.6 (Linux, x86_64). -Absolute numbers depend heavily on the machine, Python and library versions, so rely on the ratios between columns rather than on the values. On shared machines single runs wobble - `--runs=3` takes the per-cell minimum of three complete runs. Re-run on your own hardware, from the repository root: +Absolute numbers depend heavily on the machine, Python and library versions, so rely on the ratios between columns rather than on the values. On shared machines single runs wobble, so `--runs=3` measures everything three times and reports the fastest measurement of each cell. Re-run on your own hardware, from the repository root: ```sh -python -m numpy_groupies.benchmarks.generic --pandas --numbagg +python -m numpy_groupies.benchmarks --ufunc --numbagg --pandas --runs=3 ``` -Without `--pandas` the pandas column is omitted, and the pure-Python implementation only appears with `--purepy`. The optional numbagg column (from the `numbagg` package, included in the dev extras) appears with `--numbagg`. numbagg ships nan-guarded grouped reductions only, so its plain-func numbers run through the nan twins - on the benchmark's nan-free plain data the results coincide. - -Columns: `ufunc` is the `ufunc.at`-based implementation (incomplete, for benchmarking only), `numpy` the default NumPy implementation, `numba` the Numba implementation, `numbagg` the numbagg grouped-reduction wrapper, and `pandas` the pandas `groupby` wrapper. `----` means the function is not available in that implementation. The row `arbitrary` is a custom Python callable passed as `func`. +Which columns are measured by default, what the optional ones are and how the output arrives is what `--help` explains. `----` in the table below means the function is not available in that implementation, the row `arbitrary` is a custom Python callable passed as `func`. ## Results -| function | ufunc | numpy | numba | numbagg | pandas | -| ------------ | -----: | ------: | -----: | ------: | ------: | -| sum | 0.866 | 0.924 | 0.434 | 0.485 | 15.891 | -| prod | 0.920 | 0.928 | 0.468 | 0.444 | 15.191 | -| min | 1.730 | 1.687 | 0.611 | 0.815 | 15.018 | -| max | 1.712 | 1.713 | 0.638 | 0.876 | 15.154 | -| len | 0.942 | 0.830 | 0.373 | 0.471 | 14.100 | -| all | 32.574 | 1.805 | 0.823 | 1.567 | 15.725 | -| any | 30.603 | 2.506 | 0.644 | 1.568 | 15.675 | -| anynan | 5.221 | 0.890 | 0.623 | 2.674 | 15.258 | -| allnan | 7.607 | 2.122 | 0.826 | 3.255 | 15.315 | -| mean | ---- | 1.559 | 0.511 | 0.556 | 15.985 | -| median | ---- | 48.409 | 11.493 | ---- | 27.700 | -| trapezoid | ---- | 3.895 | 0.667 | ---- | ---- | -| std | ---- | 3.094 | 0.673 | 0.677 | 16.823 | -| var | ---- | 3.187 | 0.675 | 0.675 | 16.696 | -| first | ---- | 1.765 | 0.362 | 0.313 | 15.390 | -| last | ---- | 0.952 | 0.373 | 0.431 | 15.573 | -| argmax | ---- | 2.671 | 0.606 | 0.888 | 15.260 | -| argmin | ---- | 4.992 | 0.586 | 0.838 | 14.827 | -| nansum | ---- | 4.169 | 0.443 | 1.755 | 22.259 | -| nanprod | ---- | 4.115 | 0.479 | 1.770 | 20.854 | -| nanmin | ---- | 4.875 | 0.977 | 2.011 | 21.029 | -| nanmax | ---- | 4.945 | 1.010 | 2.085 | 21.752 | -| nanlen | ---- | 2.477 | 0.465 | 1.712 | 20.736 | -| nanall | ---- | 4.985 | 1.058 | 3.249 | 21.618 | -| nanany | ---- | 5.664 | 0.634 | 2.659 | 21.525 | -| nanmean | ---- | 4.678 | 1.634 | 1.758 | 22.228 | -| nanmedian | ---- | 42.432 | 10.681 | ---- | 30.505 | -| nantrapezoid | ---- | 6.558 | 1.772 | ---- | ---- | -| nanvar | ---- | 6.255 | 1.860 | 1.921 | 22.432 | -| nanstd | ---- | 6.216 | 1.893 | 1.929 | 22.450 | -| nanfirst | ---- | 5.077 | 0.393 | 0.368 | 22.225 | -| nanlast | ---- | 4.370 | 0.434 | 1.627 | 21.577 | -| nanargmin | ---- | 6.385 | 0.774 | 1.958 | 16.259 | -| nanargmax | ---- | 3.998 | 0.807 | 2.028 | 16.571 | -| cumsum | ---- | 44.838 | 0.583 | ---- | 7.322 | -| cumprod | ---- | ---- | 0.627 | ---- | 13.258 | -| cummax | ---- | ---- | 0.758 | ---- | 13.085 | -| cummin | ---- | ---- | 0.786 | ---- | 12.785 | -| arbitrary | ---- | 144.095 | 42.973 | ---- | 108.594 | -| sort | ---- | 135.146 | 25.992 | ---- | ---- | +| function | numpy | numba | ufunc | numbagg | pandas | +| ------------ | -----------: | -----------: | -----------: | -----------: | -----------: | +| sum | 0.953 | 0.444 | 0.903 | 0.472 | 15.961 | +| prod | 0.886 | 0.456 | 0.909 | 0.443 | 15.274 | +| min | 1.675 | 0.618 | 1.687 | 0.837 | 16.073 | +| max | 1.684 | 0.632 | 1.767 | 0.849 | 16.259 | +| len | 0.845 | 0.382 | 0.929 | 0.488 | 15.132 | +| all | 1.756 | 0.819 | 32.878 | 1.575 | 16.028 | +| any | 2.465 | 0.650 | 30.555 | 1.562 | 16.047 | +| anynan | 0.920 | 0.614 | 5.350 | 2.658 | 15.452 | +| allnan | 2.109 | 0.831 | 7.601 | 3.289 | 15.752 | +| mean | 1.524 | 0.511 | ---- | 0.545 | 16.978 | +| median | 48.770 | 11.517 | ---- | ---- | 28.566 | +| trapezoid | 3.814 | 0.653 | ---- | ---- | ---- | +| std | 3.093 | 0.656 | ---- | 0.693 | 17.020 | +| var | 3.133 | 0.667 | ---- | 0.685 | 16.854 | +| first | 1.790 | 0.293 | ---- | 0.308 | 15.966 | +| last | 0.938 | 0.366 | ---- | 0.425 | 15.951 | +| argmax | 2.625 | 0.624 | ---- | 0.908 | 15.760 | +| argmin | 5.136 | 0.583 | ---- | 0.847 | 15.562 | +| nansum | 4.237 | 0.446 | ---- | 1.778 | 21.705 | +| nanprod | 4.159 | 0.449 | ---- | 1.741 | 21.129 | +| nanmin | 4.844 | 0.936 | ---- | 1.996 | 21.081 | +| nanmax | 4.923 | 1.005 | ---- | 2.083 | 21.219 | +| nanlen | 2.495 | 0.448 | ---- | 1.713 | 20.246 | +| nanall | 4.929 | 1.057 | ---- | 3.270 | 21.970 | +| nanany | 5.698 | 0.642 | ---- | 2.656 | 21.667 | +| nanmean | 4.705 | 1.650 | ---- | 1.797 | 21.271 | +| nanmedian | 42.596 | 10.618 | ---- | ---- | 30.263 | +| nantrapezoid | 6.610 | 1.771 | ---- | ---- | ---- | +| nanvar | 6.111 | 1.860 | ---- | 1.922 | 21.491 | +| nanstd | 6.135 | 1.863 | ---- | 1.934 | 21.645 | +| nanfirst | 4.872 | 0.348 | ---- | 0.317 | 21.347 | +| nanlast | 4.224 | 0.449 | ---- | 1.610 | 21.495 | +| nanargmin | 6.409 | 0.758 | ---- | 1.939 | 16.622 | +| nanargmax | 3.947 | 0.794 | ---- | 2.019 | 16.535 | +| cumsum | 45.878 | 0.583 | ---- | ---- | 7.567 | +| cumprod | ---- | 0.563 | ---- | ---- | 13.840 | +| cummax | ---- | 0.750 | ---- | ---- | 14.388 | +| cummin | ---- | 0.761 | ---- | ---- | 13.083 | +| arbitrary | 142.543 | 43.163 | ---- | ---- | 109.890 | +| sort | 136.157 | 25.961 | ---- | ---- | ---- | ## Reading the table -- **numba vs numpy:** Numba is faster for every function both provide - by 1.4–13× across the reductions (`sum` 2.1×, `argmin` 8.5×, `nanfirst` 13×), by 5.2× on `sort` (26 ms against 135 ms), and by far the most on `cumsum` (about 77×). -- **numba vs numbagg:** par on the light scans (`sum` 0.43 vs 0.49 ms, `nanfirst` 0.39 vs 0.37 ms, and `first` even a sixth behind); numba leads wherever the nan handling branches per value (`anynan` 4.3×, `nansum` 4.0×, `nanlast` 3.8×, `nanargmax` 2.5×), and is the only one of the two with `median`, the cumulatives, `trapezoid` and `sort`. -- **numpy vs pandas:** the NumPy implementation is about 5× (`var`, `std`) to 17× (`sum`, `anynan`, `len`; `prod` 16×) faster than pandas on common reductions. The exceptions are `median`, where pandas (28 ms) beats NumPy (48 ms) but not Numba (11 ms), and the cumulatives, where pandas is six times faster than NumPy (`cumsum` 7.3 against 44.8 ms) and still twelve times slower than Numba (0.58 ms). -- **ufunc:** `ufunc.at` was historically slow; NumPy 1.25 narrowed the gap considerably, but this implementation remains incomplete and `all`/`any` are still very slow. +- **numba vs numpy:** Numba is faster for every function both provide - by 1.5-6.1× across the plain reductions (`sum` 2.1×, `first` 6.1×), by 2.9-14× on the nan ones (`nanmean` 2.9×, `nanfirst` 14×), by 4.2-8.8× on the argreductions (`argmin` 8.8×), by 5.2× on `sort` (26 ms against 136 ms), and by far the most on `cumsum` (about 79×). +- **numba vs numbagg:** par on the light scans (`sum` 0.44 against 0.47 ms, `first` 0.29 against 0.31 ms, while `prod` and `nanfirst` go to numbagg by a few per cent); numba leads wherever the nan handling branches per value (`anynan` 4.3×, `nansum` 4.0×, `nanlen` 3.8×, `nanlast` 3.6×, `nanargmax` 2.5×), and is the only one of the two with `median`, the cumulatives, `trapezoid` and `sort`. +- **numpy vs pandas:** the NumPy implementation is 5.4× (`var`) to 17.9× (`len`) faster on the plain reductions, up to 8.1× on the nan ones (`nanlen`) and 2.6-6.0× on the argreductions. The exceptions are `median` and `nanmedian`, where pandas (28.6 ms and 30.3 ms) beats NumPy (48.8 ms and 42.6 ms) but not Numba (11.5 ms and 10.6 ms), and the cumulatives, where pandas is six times faster than NumPy (`cumsum` 7.6 against 45.9 ms) and still thirteen times slower than Numba (0.58 ms). +- **ufunc:** `ufunc.at` was historically slow; NumPy 1.25 brought [major speed improvements](https://numpy.org/doc/stable/release/1.25.0-notes.html) and on the arithmetic scans the column now stands level with NumPy (`sum` 0.90 against 0.95 ms, `prod`, `min`, `max` and `len` within a few per cent). It is six times behind on `anynan` and twelve to nineteen times on `all` and `any`, and the implementation remains incomplete. - **Expensive everywhere:** `median`, `sort` and custom callables. ## Sorting across group counts @@ -78,14 +76,14 @@ Columns: `ufunc` is the `ufunc.at`-based implementation (incomplete, for benchma `sort` emits one value per input item and orders every group on its own, so its cost turns on the group *size* - a few large blocks to order against many small ones. `--sort` measures that function alone over a hundredth, one, ten and a hundred times the group count of the table above, on the same 500,000 values: ```sh -python -m numpy_groupies.benchmarks.generic --sort --runs=3 +python -m numpy_groupies.benchmarks --sort --runs=3 ``` | groups | values per group | numpy | numba | | -----: | ---------------: | ------: | -----: | -| 100 | 5000 | 124.895 | 34.421 | -| 1000 | 500 | 141.278 | 26.088 | -| 10000 | 50 | 159.103 | 21.794 | -| 100000 | 5 | 167.581 | 11.088 | +| 100 | 5000 | 126.200 | 34.335 | +| 1000 | 500 | 143.293 | 26.167 | +| 10000 | 50 | 157.599 | 21.341 | +| 100000 | 5 | 167.970 | 12.095 | -NumPy slows down monotonically as the groups multiply. The Numba column dips in the middle: with few large groups the comparison sort of a whole block is what costs, while many tiny groups are ordered by insertion sort in the flat array, without the per-group view that a library sort has to build - which is three times the speed it had before at five values per group. The other implementations have no column here because none of them provides `sort`. +NumPy slows down monotonically as the groups multiply, from 126 to 168 ms. The Numba column falls the whole way down, to 12 ms: with few large groups the comparison sort of a whole block is what costs, while many tiny groups are ordered by insertion sort in the flat array, without the per-group view that a library sort has to build - which is almost three times the speed of the five-thousand-value groups. The other implementations have no column here because none of them provides `sort`. diff --git a/numpy_groupies/benchmarks/__init__.py b/numpy_groupies/benchmarks/__init__.py index e69de29..1a0895d 100644 --- a/numpy_groupies/benchmarks/__init__.py +++ b/numpy_groupies/benchmarks/__init__.py @@ -0,0 +1,420 @@ +"""The benchmark suite of the aggregation implementations. + +Run it with ``python -m numpy_groupies.benchmarks`` - the module measures every +function of ``func_list`` once per implementation and prints the milliseconds +per call as a Markdown table. NumPy and Numba are always measured, ``--help`` +explains the columns to ask for, ``--all`` takes them all. A single run, the +default, prints its rows as they are measured; a repeated one reports a counter +and then the fastest measurement of all runs. +""" + +import argparse +import platform +import sys +import timeit + +import numpy as np + +from numpy_groupies import aggregate_numpy, aggregate_purepy +from numpy_groupies.benchmarks import _grouploop, _ufunc + +try: + from numpy_groupies import aggregate_numba +except ImportError: + aggregate_numba = None +try: + from numpy_groupies.benchmarks import _numbagg +except ImportError: + _numbagg = None +try: + from numpy_groupies.benchmarks import _pandas +except ImportError: + _pandas = None + + +def arbitrary(iterator): + tmp = 0 + for i, x in enumerate(iterator, 1): + tmp += x**i + return tmp + + +func_list = ( + np.sum, + np.prod, + np.min, + np.max, + len, + np.all, + np.any, + "anynan", + "allnan", + np.mean, + np.median, + np.trapezoid, + np.std, + np.var, + "first", + "last", + "argmax", + "argmin", + np.nansum, + np.nanprod, + np.nanmin, + np.nanmax, + "nanlen", + "nanall", + "nanany", + np.nanmean, + np.nanmedian, + "nantrapezoid", + np.nanvar, + np.nanstd, + "nanfirst", + "nanlast", + "nanargmin", + "nanargmax", + "cumsum", + "cumprod", + "cummax", + "cummin", + arbitrary, + "sort", +) + + +def benchmark_data(size=5e5, seed=100, groups=1000): + rnd = np.random.RandomState(seed=seed) + group_idx = rnd.randint(0, int(groups), int(size)) + a = rnd.random_sample(group_idx.size) + a[a > 0.8] = 0 + nana = a.copy() + nana[(nana < 0.2) & (nana != 0)] = np.nan + nan_share = np.mean(np.isnan(nana)) + assert 0.15 < nan_share < 0.25, f"{nan_share * 100:3f}% nans" + return a, nana, group_idx + + +def single_run(implementations, group_idx, a, nana, *, repeat, raise_errors, funcs=func_list, on_row=None): + """One complete measurement pass - returns the table as + ``{func_name: [cell, ...]}``, where a cell is a float in + milliseconds, None for unavailable combinations or a marker string + ("ERROR"/"FAIL"). + + ``on_row(label, cells)`` is called after every function, so the caller can + report the row as it is measured instead of waiting for the pass. ``repeat`` + is the number of timed executions of a cell, capped per column by + ``MAX_REPEAT``.""" + table = {} + for func in funcs: + func_name = getattr(func, "__name__", func) + used_a = nana if "nan" in func_name else a + cells = [] + results = [] + for name, aggregate_fn in implementations: + try: + res = aggregate_fn(group_idx, used_a, func=func) + except NotImplementedError: + cells.append(None) + continue + except Exception: + if raise_errors: + raise + cells.append("ERROR") + else: + results.append(res) + try: + np.testing.assert_array_almost_equal(res, results[0]) + except AssertionError: + cells.append("FAIL") + else: + timer = timeit.Timer(lambda: aggregate_fn(group_idx, used_a, func=func)) + t0 = min(timer.repeat(repeat=min(repeat, MAX_REPEAT.get(name, repeat)), number=1)) + cells.append(t0 * 1000) + table[func_name] = cells + if on_row is not None: + on_row(func_name, cells) + return table + + +def _merge_runs(tables): + """Per-cell minimum over the complete runs - filters the system noise + (cpu frequency, background load) that makes single runs wobble.""" + if len(tables) == 1: + return tables[0] + table = {} + for func_name in tables[0]: + merged = [] + for col in zip(*(t[func_name] for t in tables)): + numbers = [c for c in col if isinstance(c, float)] + if numbers: + merged.append(min(numbers)) + else: + merged.append(next((c for c in col if c is not None), None)) + table[func_name] = merged + return table + + +# Width of every measured column: a Markdown table reads just as well in a +# terminal when the cells are padded, and padding does not hurt the rendering. +CELL_WIDTH = 12 +NO_IMPL = "----" +# The reference paths are timed once whatever --repeat asks for: their number +# is a ratio rather than a result this package competes on, and every further +# execution only makes the whole measurement wait for them. +MAX_REPEAT = {"purepy": 1, "grouploop": 1} + + +def _cell_text(cell): + """The cell as text: the milliseconds, ``NO_IMPL`` for a combination the + implementation does not offer, or the ERROR/FAIL marker.""" + if cell is None: + return NO_IMPL + if isinstance(cell, float): + return f"{cell:.3f}" + return str(cell) + + +def _table_row(label, cells, label_width): + return "| " + " | ".join([label.ljust(label_width)] + [_cell_text(cell).rjust(CELL_WIDTH) for cell in cells]) + " |" + + +def _table_rule(header, implementations, label_width): + """The two Markdown lines that introduce a table - the numeric columns + carry a trailing colon, which is what right-aligns them.""" + names = [name for name, _ in implementations] + return [ + "| " + " | ".join([header.ljust(label_width)] + [name.rjust(CELL_WIDTH) for name in names]) + " |", + "| " + " | ".join(["-" * label_width] + ["-" * (CELL_WIDTH - 1) + ":" for _ in names]) + " |", + ] + + +def _table_lines(header, rows, implementations, label_width): + """``rows`` are (row label, cells) pairs - the complete Markdown table.""" + return [ + *_table_rule(header, implementations, label_width), + *(_table_row(label, cells, label_width) for label, cells in rows), + ] + + +def _open_report(header, implementations, labels, runs): + """How the measurement reports itself while it works. + + A single run is its own report: the rows stream out as Markdown, header + first, and nothing is repeated at the end. Every further run would print + the same table again, so a repeated measurement announces itself with a + counter alone and the merged minimum is the only output.""" + label_width = max(len(header), *(len(label) for label in labels)) + if runs > 1: + return None, label_width + print(*_table_rule(header, implementations, label_width), sep="\n", flush=True) + return sys.stdout, label_width + + +def _versions_line(implementations): + """The setup of this run as an italic Markdown paragraph.""" + implementation_names = [name for name, _ in implementations] + postfix = "" + if "numba" in implementation_names: + import numba + + postfix += f", Numba {numba.__version__}" + if "numbagg" in implementation_names: + import numbagg + + postfix += f", numbagg {numbagg.__version__}" + if "pandas" in implementation_names: + import pandas + + postfix += f", Pandas {pandas.__version__}" + return ( + f"*{platform.system()}({platform.machine()}), Python {sys.version.split()[0]}, " + f"Numpy {np.version.version}{postfix}*" + ) + + +def benchmark(implementations, repeat=5, size=5e5, seed=100, raise_errors=False, runs=1): + a, nana, group_idx = benchmark_data(size=size, seed=seed) + labels = [getattr(func, "__name__", func) for func in func_list] + out, label_width = _open_report("function", implementations, labels, runs) + + def report(label, cells): + if out is not None: + print(_table_row(label, cells, label_width), file=out, flush=True) + + tables = [] + for run in range(runs): + if runs > 1: + print(f"run {run + 1}/{runs}", file=sys.stderr, flush=True) + tables.append( + single_run( + implementations, + group_idx, + a, + nana, + repeat=repeat, + raise_errors=raise_errors, + on_row=report, + ) + ) + + if runs > 1: + table = _merge_runs(tables) + print( + *_table_lines("function", [(label, table[label]) for label in labels], implementations, label_width), + _versions_line(implementations), + sep="\n", + ) + else: + print(_versions_line(implementations)) + + +# The group counts of the ``--sort`` series - the default of the main table and +# a hundredth and a hundred times it. ``sort`` orders every group on its own, +# so its cost turns on the group SIZE (a few large blocks against many small +# ones) and the main table would only ever show one point of that curve. +SORT_GROUPS = (100, 1_000, 10_000, 100_000) + + +def sort_benchmark(implementations, *, repeat=5, size=5e5, seed=100, raise_errors=False, runs=1): + """Measure ``func='sort'`` alone, over the group counts of ``SORT_GROUPS``.""" + labels = [str(groups) for groups in SORT_GROUPS] + out, label_width = _open_report("groups", implementations, labels, runs) + + def report(label, cells): + if out is not None: + print(_table_row(label, cells, label_width), file=out, flush=True) + + rows = [] + for groups, label in zip(SORT_GROUPS, labels): + a, nana, group_idx = benchmark_data(size=size, seed=seed, groups=min(groups, int(size))) + tables = [] + for run in range(runs): + if runs > 1: + print(f"groups {label}: run {run + 1}/{runs}", file=sys.stderr, flush=True) + tables.append( + single_run( + implementations, + group_idx, + a, + nana, + repeat=repeat, + raise_errors=raise_errors, + funcs=("sort",), + on_row=lambda _func_name, cells: report(label, cells), + ) + ) + rows.append((label, _merge_runs(tables)["sort"])) + + if runs > 1: + print(*_table_lines("groups", rows, implementations, label_width), _versions_line(implementations), sep="\n") + else: + print(_versions_line(implementations)) + + +# The columns in the order they appear in the table. +COLUMN_ORDER = ("numpy", "numba", "purepy", "grouploop", "ufunc", "numbagg", "pandas") +# of those, the ones that have to be asked for +OPTIONAL_COLUMNS = COLUMN_ORDER[2:] + + +def _build_implementations(args): + """the participants as (column name, aggregate callable) pairs, in + ``COLUMN_ORDER``, the optional ones only when they were asked for""" + if aggregate_numba is None: + raise SystemExit("this benchmark needs the numba package installed") + aggregates = { + "numpy": aggregate_numpy.aggregate, + "numba": aggregate_numba.aggregate, + "purepy": aggregate_purepy.aggregate, + "grouploop": _grouploop.aggregate, + "ufunc": _ufunc.aggregate, + # the optional dependencies are None when they are not installed + "numbagg": _numbagg.aggregate if _numbagg is not None else None, + "pandas": _pandas.aggregate if _pandas is not None else None, + } + impls = [] + for name in COLUMN_ORDER: + if name in OPTIONAL_COLUMNS and not getattr(args, name): + continue + aggregate = aggregates[name] + if aggregate is None: + raise SystemExit(f"--{name} needs the optional {name} package installed") + impls.append((name, aggregate)) + return impls + + +def main(argv=None): + parser = argparse.ArgumentParser( + prog="python -m numpy_groupies.benchmarks", + description="Benchmark the aggregate implementations against each other.\n" + "The numpy and numba columns are always measured, the rest is asked for.", + epilog="""The table is Markdown: milliseconds per call, '----' where an implementation +has no such function, and a line naming the versions it was measured with. A +single run prints its rows as they are measured; a repeated run replaces them +with a counter and prints the fastest measurement of all runs once at the end. + +numbagg ships nan-guarded reductions only, so its numbers for the plain +functions are those of the nan twins - on the nan-free data of this benchmark +the two coincide.""", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + parser.add_argument("--all", action="store_true", help="include every optional column listed below") + parser.add_argument( + "--purepy", + action="store_true", + help="the pure-Python implementation, the reference the others are checked against", + ) + parser.add_argument( + "--grouploop", + action="store_true", + help="numpy with its fast path disabled, what a func it does not recognise costs", + ) + parser.add_argument( + "--ufunc", action="store_true", help="the ufunc.at-based implementation, incomplete and for comparison only" + ) + parser.add_argument( + "--numbagg", + action="store_true", + help="the numbagg grouped reductions (optional dependency, part of the dev extras)", + ) + parser.add_argument("--pandas", action="store_true", help="the pandas groupby wrapper (optional dependency)") + parser.add_argument( + "--sort", + action="store_true", + help=f"measure 'sort' only, over the group counts {'-'.join(str(g) for g in SORT_GROUPS)} instead of the full table", + ) + parser.add_argument( + "--runs", type=int, default=1, help="complete measurement passes, the per-cell minimum is reported (default: 1)" + ) + parser.add_argument( + "--repeat", + type=int, + default=5, + help=( + "timed executions per cell, the fastest counts, except for the pure-Python and " + "the group-loop column (default: 5)" + ), + ) + parser.add_argument("--size", type=float, default=5e5, help="number of values in the benchmark data (default: 5e5)") + parser.add_argument( + "--raise-errors", + action="store_true", + help="abort on the first failing implementation instead of printing ERROR", + ) + args = parser.parse_args(argv) + if args.all: + for column in OPTIONAL_COLUMNS: + setattr(args, column, True) + measure = sort_benchmark if args.sort else benchmark + measure( + _build_implementations(args), + repeat=args.repeat, + size=args.size, + raise_errors=args.raise_errors, + runs=max(1, args.runs), + ) + + +if __name__ == "__main__": + main() diff --git a/numpy_groupies/benchmarks/__main__.py b/numpy_groupies/benchmarks/__main__.py new file mode 100644 index 0000000..2a9216a --- /dev/null +++ b/numpy_groupies/benchmarks/__main__.py @@ -0,0 +1,15 @@ +"""Entry point of ``python -m numpy_groupies.benchmarks``.""" + +import signal + +from . import main + +if hasattr(signal, "SIGPIPE"): + # a reader that went away (``| head``) ends the run quietly, as it does for + # any other command line filter; left alone, python prints a traceback + signal.signal(signal.SIGPIPE, signal.SIG_DFL) + +try: + main() +except BrokenPipeError: + pass diff --git a/numpy_groupies/benchmarks/_grouploop.py b/numpy_groupies/benchmarks/_grouploop.py new file mode 100644 index 0000000..a82f641 --- /dev/null +++ b/numpy_groupies/benchmarks/_grouploop.py @@ -0,0 +1,33 @@ +"""The numpy implementation with its fast path disabled, imported with +``--grouploop`` in :mod:`numpy_groupies.benchmarks`. + +Wrapping the function in a lambda keeps ``aggregate_numpy`` from recognising it +as one of its own, so it groups the values and loops over the groups instead - +the column therefore measures what a ``func`` that the implementation does not +know costs, which is the price of the recognition itself.""" + +from operator import itemgetter + +from .. import aggregate_numpy +from ..utils import allnan, anynan, nanfirst, nanlast + +# the string names of func_list that a per-group callable can stand in for - +# everything else (the cumulatives, the argreductions, sort) needs the +# implementation's own handling and has no group-loop equivalent +extrafuncs = { + "allnan": allnan, + "anynan": anynan, + "first": itemgetter(0), + "last": itemgetter(-1), + "nanfirst": nanfirst, + "nanlast": nanlast, +} + + +def aggregate(*args, **kwargs): + """``aggregate_numpy.aggregate`` on a function it cannot recognise.""" + func = kwargs.pop("func") + func = extrafuncs.get(func, func) + if isinstance(func, str): + raise NotImplementedError("Grouploop needs to be called with a function") + return aggregate_numpy.aggregate(*args, func=lambda x: func(x), **kwargs) diff --git a/numpy_groupies/benchmarks/generic.py b/numpy_groupies/benchmarks/generic.py deleted file mode 100755 index d523a40..0000000 --- a/numpy_groupies/benchmarks/generic.py +++ /dev/null @@ -1,331 +0,0 @@ -#!/usr/bin/python -B - -import argparse -import platform -import sys -import timeit -from operator import itemgetter - -import numpy as np - -from numpy_groupies import aggregate_numpy, aggregate_purepy -from numpy_groupies.benchmarks import _ufunc -from numpy_groupies.utils import allnan, anynan, nanfirst, nanlast - -try: - from numpy_groupies import aggregate_numba -except ImportError: - aggregate_numba = None -try: - from numpy_groupies.benchmarks import _numbagg -except ImportError: - _numbagg = None -try: - from numpy_groupies.benchmarks import _pandas -except ImportError: - _pandas = None - - -def aggregate_grouploop(*args, **kwargs): - """wraps func in lambda which prevents aggregate_numpy from - recognising and optimising it. Instead it groups and loops.""" - extrafuncs = { - "allnan": allnan, - "anynan": anynan, - "first": itemgetter(0), - "last": itemgetter(-1), - "nanfirst": nanfirst, - "nanlast": nanlast, - } - func = kwargs.pop("func") - func = extrafuncs.get(func, func) - if isinstance(func, str): - raise NotImplementedError("Grouploop needs to be called with a function") - return aggregate_numpy.aggregate(*args, func=lambda x: func(x), **kwargs) - - -def arbitrary(iterator): - tmp = 0 - for i, x in enumerate(iterator, 1): - tmp += x**i - return tmp - - -func_list = ( - np.sum, - np.prod, - np.min, - np.max, - len, - np.all, - np.any, - "anynan", - "allnan", - np.mean, - np.median, - np.trapezoid, - np.std, - np.var, - "first", - "last", - "argmax", - "argmin", - np.nansum, - np.nanprod, - np.nanmin, - np.nanmax, - "nanlen", - "nanall", - "nanany", - np.nanmean, - np.nanmedian, - "nantrapezoid", - np.nanvar, - np.nanstd, - "nanfirst", - "nanlast", - "nanargmin", - "nanargmax", - "cumsum", - "cumprod", - "cummax", - "cummin", - arbitrary, - "sort", -) - - -def benchmark_data(size=5e5, seed=100, groups=1000): - rnd = np.random.RandomState(seed=seed) - group_idx = rnd.randint(0, int(groups), int(size)) - a = rnd.random_sample(group_idx.size) - a[a > 0.8] = 0 - nana = a.copy() - nana[(nana < 0.2) & (nana != 0)] = np.nan - nan_share = np.mean(np.isnan(nana)) - assert 0.15 < nan_share < 0.25, f"{nan_share * 100:3f}% nans" - return a, nana, group_idx - - -def single_run(implementations, group_idx, a, nana, *, repeat, raise_errors, verbose, funcs=func_list): - """One complete measurement pass - returns the table as - ``{func_name: [cell, ...]}``, where a cell is a float in - milliseconds, None for unavailable combinations or a marker string - ("ERROR"/"FAIL").""" - table = {} - for func in funcs: - func_name = getattr(func, "__name__", func) - used_a = nana if "nan" in func_name else a - cells = [] - results = [] - if verbose: - print(func_name.ljust(13), end="", file=sys.stderr) - for name, aggregate_fn in implementations: - try: - res = aggregate_fn(group_idx, used_a, func=func) - except NotImplementedError: - cells.append(None) - if verbose: - print("----".rjust(14), end="", file=sys.stderr) - continue - except Exception: - if raise_errors: - raise - cells.append("ERROR") - if verbose: - print("ERROR".rjust(14), end="", file=sys.stderr) - else: - results.append(res) - try: - np.testing.assert_array_almost_equal(res, results[0]) - except AssertionError: - cells.append("FAIL") - if verbose: - print("FAIL".rjust(14), end="", file=sys.stderr) - else: - t0 = min( - timeit.Timer(lambda: aggregate_fn(group_idx, used_a, func=func)).repeat(repeat=repeat, number=1) - ) - cells.append(t0 * 1000) - if verbose: - print(f"{t0 * 1000:.3f}".rjust(14), end="", file=sys.stderr) - if verbose: - sys.stderr.flush() - if verbose: - print(file=sys.stderr) - table[func_name] = cells - return table - - -def _merge_runs(tables): - """Per-cell minimum over the complete runs - filters the system noise - (cpu frequency, background load) that makes single runs wobble.""" - if len(tables) == 1: - return tables[0] - table = {} - for func_name in tables[0]: - merged = [] - for col in zip(*(t[func_name] for t in tables)): - numbers = [c for c in col if isinstance(c, float)] - if numbers: - merged.append(min(numbers)) - else: - merged.append(next((c for c in col if c is not None), None)) - table[func_name] = merged - return table - - -def _print_table(header, rows, implementations): - """``rows`` are (row label, cells) pairs as produced by ``single_run``.""" - print(f"{header:13}" + "".join(name.rjust(14) for name, _ in implementations)) - print("-" * (13 + 14 * len(implementations))) - for label, cells in rows: - print(label.ljust(13), end="") - for cell in cells: - text = "----" if cell is None else (f"{cell:.3f}" if isinstance(cell, float) else cell) - print(text.rjust(14), end="") - print() - - -def _print_versions(implementations): - implementation_names = [name for name, _ in implementations] - postfix = "" - if "numba" in implementation_names: - import numba - - postfix += f", Numba {numba.__version__}" - if "numbagg" in implementation_names: - import numbagg - - postfix += f", numbagg {numbagg.__version__}" - if "pandas" in implementation_names: - import pandas - - postfix += f", Pandas {pandas.__version__}" - print( - f"{platform.system()}({platform.machine()}), Python {sys.version.split()[0]}, Numpy {np.version.version}" - f"{postfix}" - ) - - -def benchmark(implementations, repeat=5, size=5e5, seed=100, raise_errors=False, runs=1): - a, nana, group_idx = benchmark_data(size=size, seed=seed) - - tables = [] - for r in range(runs): - if runs > 1: - print(f"run {r + 1}/{runs}", file=sys.stderr) - tables.append( - single_run(implementations, group_idx, a, nana, repeat=repeat, raise_errors=raise_errors, verbose=runs == 1) - ) - - table = _merge_runs(tables) - _print_table( - "function", - [(name, table[name]) for name in (getattr(func, "__name__", func) for func in func_list)], - implementations, - ) - _print_versions(implementations) - - -# The group counts of the ``--sort`` series - the default of the main table and -# a hundredth and a hundred times it. ``sort`` orders every group on its own, -# so its cost turns on the group SIZE (a few large blocks against many small -# ones) and the main table would only ever show one point of that curve. -SORT_GROUPS = (100, 1_000, 10_000, 100_000) - - -def sort_benchmark(implementations, *, repeat=5, size=5e5, seed=100, raise_errors=False, runs=1): - """Measure ``func='sort'`` alone, over the group counts of ``SORT_GROUPS``.""" - rows = [] - for groups in SORT_GROUPS: - a, nana, group_idx = benchmark_data(size=size, seed=seed, groups=min(groups, int(size))) - tables = [] - for r in range(runs): - if runs > 1: - print(f"groups {groups}: run {r + 1}/{runs}", file=sys.stderr) - tables.append( - single_run( - implementations, - group_idx, - a, - nana, - repeat=repeat, - raise_errors=raise_errors, - verbose=False, - funcs=("sort",), - ) - ) - rows.append((str(groups), _merge_runs(tables)["sort"])) - - _print_table("groups", rows, implementations) - _print_versions(implementations) - - -def _build_implementations(args): - """the participants as (column name, aggregate callable) pairs, in - column order - purepy first, numbagg before the right-most pandas""" - if aggregate_numba is None: - raise SystemExit("this benchmark needs the numba package installed") - impls = [] - if args.purepy: - impls.append(("purepy", aggregate_purepy.aggregate)) - impls.extend( - ( - ("ufunc", _ufunc.aggregate), - ("numpy", aggregate_numpy.aggregate), - ("numba", aggregate_numba.aggregate), - ) - ) - if args.numbagg: - if _numbagg is None: - raise SystemExit("--numbagg needs the optional numbagg package installed") - impls.append(("numbagg", _numbagg.aggregate)) - if args.pandas: - if _pandas is None: - raise SystemExit("--pandas needs the optional pandas package installed") - impls.append(("pandas", _pandas.aggregate)) - return impls - - -def main(argv=None): - parser = argparse.ArgumentParser( - prog="python -m numpy_groupies.benchmarks.generic", - description="Benchmark the aggregate implementations against each other " - "(default columns: ufunc, numpy, numba).", - ) - parser.add_argument("--purepy", action="store_true", help="include the pure-Python implementation") - parser.add_argument("--pandas", action="store_true", help="include the pandas implementation") - parser.add_argument( - "--numbagg", action="store_true", help="include the numbagg implementation (optional dependency, dev extras)" - ) - parser.add_argument( - "--sort", - action="store_true", - help=f"measure 'sort' only, over the group counts {'-'.join(str(g) for g in SORT_GROUPS)} instead of the full table", - ) - parser.add_argument( - "--runs", type=int, default=1, help="complete measurement passes, the per-cell minimum is reported (default: 1)" - ) - parser.add_argument( - "--repeat", type=int, default=5, help="timed executions per cell, the fastest counts (default: 5)" - ) - parser.add_argument("--size", type=float, default=5e5, help="number of values in the benchmark data (default: 5e5)") - parser.add_argument( - "--raise-errors", - action="store_true", - help="abort on the first failing implementation instead of printing ERROR", - ) - args = parser.parse_args(argv) - measure = sort_benchmark if args.sort else benchmark - measure( - _build_implementations(args), - repeat=args.repeat, - size=args.size, - raise_errors=args.raise_errors, - runs=max(1, args.runs), - ) - - -if __name__ == "__main__": - main() diff --git a/numpy_groupies/benchmarks/simple.py b/numpy_groupies/benchmarks/simple.py deleted file mode 100755 index de5d5ac..0000000 --- a/numpy_groupies/benchmarks/simple.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/python -B - -import timeit - -import numpy as np - -from numpy_groupies import aggregate_np, aggregate_py -from numpy_groupies.benchmarks._pandas import aggregate as aggregate_pd -from numpy_groupies.benchmarks._ufunc import aggregate as aggregate_ufunc -from numpy_groupies.utils import aliasing - - -def aggregate_group_loop(*args, **kwargs): - """wraps func in lambda which prevents aggregate_numpy from - recognising and optimising it. Instead it groups and loops.""" - func = kwargs["func"] - del kwargs["func"] - return aggregate_np(*args, func=lambda x: func(x), **kwargs) - - -print("-----simple examples----------") -test_a = np.array([12.0, 3.2, -15, 88, 12.9]) -test_group_idx = np.array([1, 0, 1, 4, 1]) -print("test_a: ", test_a) -print("test_group_idx: ", test_group_idx) -print("aggregate(test_group_idx, test_a):") -print(aggregate_np(test_group_idx, test_a)) # group vals by idx and sum -# array([3.2, 9.9, 0., 0., 88.]) -print("aggregate(test_group_idx, test_a, sz=8, func='min', fill_value=np.nan):") -print(aggregate_np(test_group_idx, test_a, size=8, func="min", fill_value=np.nan)) -# array([3.2, -15., nan, 88., nan, nan, nan, nan]) -print("aggregate_py(test_group_idx, test_a, sz=5, func=lambda x: ' + '.join(str(xx) for xx in x),fill_value='')") -print( - aggregate_py( - test_group_idx, - test_a, - size=5, - func=lambda x: " + ".join(str(xx) for xx in x), - fill_value="", - ) -) - - -print() -print("---------testing--------------") -print("compare against group-and-loop with numpy") -testable_funcs = {aliasing[f]: f for f in (np.sum, np.prod, np.any, np.all, np.min, np.max, np.std, np.var, np.mean)} -test_group_idx = np.random.randint(0, int(1e3), int(1e5)) -test_a = np.random.rand(int(1e5)) * 100 - 50 -test_a[test_a > 25] = 0 # for use with bool functions -for name, f in testable_funcs.items(): - numpy_loop_group = aggregate_group_loop(test_group_idx, test_a, func=f) - - for acc_func, acc_name in [ - (aggregate_np, "np-optimised"), - (aggregate_ufunc, "np-ufunc-at"), - (aggregate_py, "purepy"), - (aggregate_pd, "pandas"), - ]: - try: - test_out = acc_func(test_group_idx, test_a, func=name) - test_out = np.asarray(test_out) - if not np.allclose(test_out, numpy_loop_group.astype(test_out.dtype)): - print( - name, - acc_name, - "FAILED test, output: [" + acc_name + "; correct]...", - ) - print(np.vstack((test_out, numpy_loop_group))) - else: - print(name, acc_name, "PASSED test") - except NotImplementedError: - print(name, acc_name, "NOT IMPLEMENTED") - -print() -print("----------benchmarking-------------") -print("Note that the actual observed speedup depends on a variety of properties of the input.") -print("Here we are using 100,000 indices uniformly picked from [0, 1000).") -print("Specifically, about 25% of the values are 0 (for use with bool operations),") -print("the remainder are uniformly distributed on [-50,25).") -print("Times are scaled to 10 repetitions (actual number of reps used may not be 10).") - -print( - "".join( - [ - "function".rjust(8), - "pure-py".rjust(14), - "np-grouploop".rjust(14), - "np-ufuncat".rjust(14), - "np-optimised".rjust(14), - "pandas".rjust(14), - "ratio".rjust(15), - ] - ) -) - -for name, f in testable_funcs.items(): - print(name.rjust(8), end="") - times: list[float | None] = [None] * 5 - for ii, acc_func in enumerate( - [ - aggregate_py, - aggregate_group_loop, - aggregate_ufunc, - aggregate_np, - aggregate_pd, - ] - ): - try: - func = f if acc_func is aggregate_group_loop else name - reps = 3 if acc_func is aggregate_py else 20 - times[ii] = ( - timeit.Timer(lambda: acc_func(test_group_idx, test_a, func=func)).timeit(number=reps) / reps * 10 - ) - print(f"{times[ii] * 1000:.1f}ms".rjust(13), end="") # type: ignore[operator] # always set unless the raise above - except NotImplementedError: - print("no-impl".rjust(13), end="") - - denom = min(t for t in times if t is not None) - ratios = [("-".center(4) if t is None else str(round(t / denom, 1))).center(5) for t in times] - print(" ", (":".join(ratios)))