diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..c6763e4 --- /dev/null +++ b/python/README.md @@ -0,0 +1,91 @@ +## Numba Implementation + +A high-performance Numba JIT-compiled implementation is available in `medcouple_numba.py`. + +### Installation +```bash +pip install numba +``` +Example usage: +```python +from medcouple_numba import medcouple_matrix_full +result = medcouple_matrix_full(data) +``` +### Performance: + +The Numba implementation with parallel execution achieves sensitive speedup compared to the pure NumPy implementation. + +I have also added more descriptive tests, `test_df.py` and `test_numba.py`, that require pandas and tqdm. + +If you want to print the table in markdown format, you need to install `tabulate` and `print(summary.to_markdown())`. + +#### Numba implementation: + +**Summary of the results ordered by mean timings:** + +| Test case | mean_results | mean_timings | min_timings | max_timings | all_passed | +|:-----------------------------------------|---------------:|---------------:|--------------:|--------------:|:-------------| +| naive_loop_kernel_same_partial_False_1 | 0.246938 | 0.00855717 | 0.000555452 | 0.034805 | True | +| naive_loop_kernel_same_partial_True_0 | 0.246938 | 0.00961762 | 0.000663259 | 0.0375037 | True | +| naive_loop_kernel_same_partial_False_0 | 0.246938 | 0.0102791 | 0.000552418 | 0.0488931 | True | +| naive_loop_kernel_same_partial_True_1 | 0.246938 | 0.0105872 | 0.000536229 | 0.045166 | True | +| naive_loop_kernel_same_monotonic_False_0 | 0.246938 | 0.0164787 | 0.000534889 | 0.0921133 | True | +| naive_loop_kernel_same_full_True_1 | 0.246938 | 0.0167137 | 0.000599015 | 0.0928224 | True | +| naive_loop_kernel_same_full_False_0 | 0.246938 | 0.0169012 | 0.000556725 | 0.0990537 | True | +| naive_loop_kernel_same_full_True_0 | 0.246938 | 0.017097 | 0.000521972 | 0.0945251 | True | +| naive_loop_kernel_same_monotonic_True_0 | 0.246938 | 0.0174738 | 0.00050252 | 0.0996008 | True | +| naive_loop_kernel_same_monotonic_False_1 | 0.246938 | 0.0178723 | 0.000503608 | 0.0993816 | True | +| naive_loop_kernel_same_monotonic_True_1 | 0.246938 | 0.0181392 | 0.000615125 | 0.10721 | True | +| naive_loop_kernel_same_full_False_1 | 0.246938 | 0.0183009 | 0.000682849 | 0.10718 | True | +| naive_matrix_gen_matrix_same_partial | 0.246938 | 0.0484474 | 0.000206347 | 0.256039 | True | +| statsmodels_medcouple | 0.246938 | 0.0578549 | 0.0058599 | 0.27967 | True | +| naive_matrix_gen_matrix_same_full | 0.246938 | 0.0579054 | 0.000199287 | 0.300115 | True | +| naive_matrix_gen_matrix_same_monotonic | 0.246938 | 0.0595919 | 0.000248522 | 0.304371 | True | + +Group averages: + + Statsmodels: 0.057855s\ + naive_loop average: 0.014835s\ + naive_matrix average: 0.055315s\ + Matrix is 3.73x slower than loop\ + Loop is 3.90x faster than statsmodels\ + Matrix is 1.05x faster than statsmodels + +#### Numpy implementation: + +**Summary of the results ordered by mean timings:** + +| Test case | mean_results | mean_timings | min_timings | max_timings | all_passed | +|:-----------------------------------------|---------------:|---------------:|--------------:|--------------:|:-------------| +| statsmodels_medcouple | 0.246938 | 0.0661967 | 0.00758035 | 0.332696 | True | +| naive_matrix_gen_matrix_same_partial | 0.246938 | 0.0706036 | 0.00720444 | 0.345422 | True | +| naive_matrix_gen_matrix_same_full | 0.246938 | 0.0778316 | 0.00530281 | 0.426505 | True | +| naive_matrix_gen_matrix_same_monotonic | 0.246938 | 0.0892744 | 0.00724118 | 0.507823 | True | +| naive_loop_kernel_same_full_True_1 | 0.246938 | 0.714831 | 0.00451493 | 2.38954 | True | +| naive_loop_kernel_same_partial_False_0 | 0.246938 | 0.734988 | 0.0053756 | 2.49929 | True | +| naive_loop_kernel_same_monotonic_False_1 | 0.246938 | 0.746236 | 0.00525514 | 2.64932 | True | +| naive_loop_kernel_same_partial_True_1 | 0.246938 | 0.747932 | 0.00583358 | 2.52636 | True | +| naive_loop_kernel_same_full_False_0 | 0.246938 | 0.750571 | 0.00578721 | 2.48884 | True | +| naive_loop_kernel_same_monotonic_True_1 | 0.246938 | 0.755513 | 0.0054866 | 2.6988 | True | +| naive_loop_kernel_same_partial_True_0 | 0.246938 | 0.759996 | 0.00592637 | 2.65849 | True | +| naive_loop_kernel_same_full_False_1 | 0.246938 | 0.772812 | 0.00591261 | 2.73194 | True | +| naive_loop_kernel_same_monotonic_False_0 | 0.246938 | 0.77304 | 0.00468273 | 2.93801 | True | +| naive_loop_kernel_same_partial_False_1 | 0.246938 | 0.778947 | 0.00507722 | 2.84921 | True | +| naive_loop_kernel_same_monotonic_True_0 | 0.246938 | 0.782877 | 0.0050175 | 2.94109 | True | +| naive_loop_kernel_same_full_True_0 | 0.246938 | 0.792352 | 0.00480242 | 2.97309 | True | + +Group averages: + + Statsmodels: 0.066197s\ + naive_loop average: 0.759175s\ + naive_matrix average: 0.079237s\ + Matrix is 9.58x faster than loop\ + Loop is 11.47x slower than statsmodels\ + Matrix is 1.20x slower than statsmodels + +#### Numpy vs Numba + +Naive Loop is **51.17x** faster in Numba.\ +Matrix is 1.43x faster in Numba. + +**Loops are extremely fast in Numba!** \ No newline at end of file diff --git a/python/medcouple_numba.py b/python/medcouple_numba.py new file mode 100644 index 0000000..8cc6ff4 --- /dev/null +++ b/python/medcouple_numba.py @@ -0,0 +1,421 @@ +from typing import Callable +import numpy as np +from numba import njit, prange + + +@njit(cache=True, fastmath=False, nogil=True) +def kernel_same_partial(i, j, min_k, max_k, start_index): + # partial matrix, shortcut for loop implementation + # ensures we only add 0 for k times and ignore the other pairs as per paper + if i == j: # this only happens k times, as the antidiagonal of the square of ties + return 0.0 + else: # will be ignored + return np.nan + + +@njit(cache=True, fastmath=False, nogil=True) +def kernel_same_full(i, j, min_k, max_k, start_index): + # full matrix, necessary for matrix implementation + # simpler than paper version as it does not bother with values of k, mi, mj (troublesome as it's using a different i/j compared to the rest) + # however this is theoretically equivalent since the result is just the median as we also end up with 0 for k times, and equal numbers of 1 and -1 to offset each other, just in different positions + # also empirically equivalent based on tests + return np.sign(j - i) + + +@njit(cache=True, fastmath=False, nogil=True) +def kernel_same_monotonic(i, j, min_k, max_k, start_index): + # full matrix preserving monotonicity, necessary for fast matrix implementation + """ + # uses values of k, mi, mj according to the paper formula just to make the matrix monotonic + k = max_k - min_k + start_index + mi = i - min_k + start_index + mj = j - min_k + start_index + res = mi + mj - start_index - k + """ + # simplify above + res = i + j - min_k - max_k + return np.sign(res) + + +# Not optimised, just testing my understanding first +# May be logically easier to understand than optimised versions, especially if you just assume either value for each argument as it's all the same +@njit(cache=True, fastmath=False, nogil=True, parallel=True) +def medcouple_naive_loop(xs, kernel_same, med_centred=False, start_index=0): + # here the formulas are the same regardless of whether indexing starts at 0 or 1 due to our kernel_same implementations + xs = np.asarray(xs) # ensure xs is a numpy array + xs_sorted = np.sort(xs) + n = len(xs_sorted) + if n % 2 == 0: + med = np.float64((xs_sorted[n // 2 - 1] + xs_sorted[n // 2]) / 2.0) + else: + med = np.float64(xs_sorted[(n - 1) // 2]) + + # Build xis list + xis_list = [] + for i in range(len(xs_sorted)): + x = xs_sorted[i] + if x <= med: + if med_centred: + shifted_x = x - med + else: + shifted_x = x + xis_list.append((shifted_x, i + start_index)) + + # Build xjs list + xjs_list = [] + for j in range(len(xs_sorted)): + x = xs_sorted[j] + if x >= med: + if med_centred: + shifted_x = x - med + else: + shifted_x = x + xjs_list.append((shifted_x, j + start_index)) + + # Build ks list + ks_list = [] + for k in range(len(xs_sorted)): + x = xs_sorted[k] + if x == med: + ks_list.append(k + start_index) + + if len(ks_list) > 0: + min_k = ks_list[0] + max_k = ks_list[0] + for k in ks_list: + if k < min_k: + min_k = k + if k > max_k: + max_k = k + else: # these values will not be used + min_k = np.nan + max_k = np.nan + + # Pre-allocate buffer array for parallel computation + n_xjs = len(xjs_list) + n_xis = len(xis_list) + buffer = np.empty(n_xjs * n_xis, dtype=np.float64) + + # Parallel iteration through xjs and xis + for idx_j in prange(n_xjs): + # Access in reverse order + xj, j = xjs_list[n_xjs - 1 - idx_j] + for idx_i in range(n_xis): + # Access in reverse order + xi, i = xis_list[n_xis - 1 - idx_i] + + # Calculate linear index for buffer + buffer_idx = idx_j * n_xis + idx_i + + # if xi == med and xj == med: + if xi == xj: # simplified as this is only possible when they equal the median + item = kernel_same(i, j, min_k, max_k, start_index) + else: + # Calculate kernel_diff inline + if med_centred: + num = xj + xi + else: + num = (xj - med) - (med - xi) + den = xj - xi + item = num / den + buffer[buffer_idx] = item + + return np.nanmedian(buffer) + + +@njit(cache=True, fastmath=False, nogil=True) +def gen_matrix_same_partial(i_idx, j_idx, z_index, zs): + same_ref = j_idx == i_idx + return np.where( + same_ref, + np.zeros_like(same_ref, dtype=np.float64), + np.full_like(same_ref, np.nan, dtype=np.float64), + ) + + +@njit(cache=True, fastmath=False, nogil=True) +def gen_matrix_same_full(i_idx, j_idx, z_index, zs): + same_ref = j_idx - i_idx + return np.where( + same_ref == 0, + np.zeros_like(same_ref, dtype=np.float64), + np.where(same_ref > 0, np.ones_like(same_ref, dtype=np.float64), np.full_like(same_ref, -1.0, dtype=np.float64)), + ) + + +@njit(cache=True, fastmath=False, nogil=True) +def gen_matrix_same_monotonic(i_idx, j_idx, z_index, zs): + k_idx = z_index[zs == 0] + max_k = np.max(k_idx) + min_k = np.min(k_idx) + same_ref = i_idx + j_idx - min_k - max_k + return np.where( + same_ref == 0, + np.zeros_like(same_ref, dtype=np.float64), + np.where(same_ref > 0, np.ones_like(same_ref, dtype=np.float64), np.full_like(same_ref, -1.0, dtype=np.float64)), + ) + + +# slightly slower than statsmodels but hopefully more intuitive as logic is close to the loop version +@njit(cache=True, fastmath=False, nogil=True) +def medcouple_naive_matrix( + xs: np.ndarray, + gen_matrix_same: Callable[[np.ndarray, np.ndarray, np.ndarray], np.ndarray], +): # kernel_same=gen_matrix_same, med_centred=True, start_index=0 + xs = np.asarray(xs) # ensure xs is a numpy array + xs = np.sort(xs) # sorts and turns it into numpy array + n = len(xs) + if n % 2 == 0: + med = np.float64((xs[n // 2 - 1] + xs[n // 2]) / 2.0) + else: + med = np.float64(xs[(n - 1) // 2]) + zs = np.flip(xs - med) + zis = zs[zs <= 0] + zjs = zs[zs >= 0] + zjs = zjs[:, None] # transpose to vertical + num = zjs + zis + den = zjs - zis + if len(zs[zs == 0]): # any ties with median + z_index = np.flip(np.indices(np.shape(zs))[0]) + zis_idx = z_index[zs <= 0] + zjs_idx = z_index[zs >= 0] + zjs_idx = zjs_idx[:, None] # transpose to vertical + i_idx = np.ones_like(zjs_idx) * zis_idx + j_idx = np.ones_like(zis_idx) * zjs_idx + same = gen_matrix_same(i_idx, j_idx, z_index, zs) + else: + same = np.full_like(den, np.nan, dtype=np.float64) + # with np.errstate(invalid='ignore'): # since we are handling the divide by zero + res = np.where(den == 0, same, num / den) + # print(res) + return np.nanmedian(res) + +# Specialized versions with inlined gen_matrix_same logic (can be cached properly) +@njit(cache=True, fastmath=False, nogil=True) +def medcouple_matrix_partial(xs): + xs = np.asarray(xs) + xs = np.sort(xs) + n = len(xs) + if n % 2 == 0: + med = np.float64((xs[n // 2 - 1] + xs[n // 2]) / 2.0) + else: + med = np.float64(xs[(n - 1) // 2]) + zs = np.flip(xs - med) + zis = zs[zs <= 0] + zjs = zs[zs >= 0] + zjs = zjs[:, None] + num = zjs + zis + den = zjs - zis + if len(zs[zs == 0]): + z_index = np.flip(np.indices(np.shape(zs))[0]) + zis_idx = z_index[zs <= 0] + zjs_idx = z_index[zs >= 0] + zjs_idx = zjs_idx[:, None] + i_idx = np.ones_like(zjs_idx) * zis_idx + j_idx = np.ones_like(zis_idx) * zjs_idx + # Inline gen_matrix_same_partial + same_ref = j_idx == i_idx + same = np.where(same_ref, np.zeros_like(same_ref, dtype=np.float64), np.full_like(same_ref, np.nan, dtype=np.float64)) + else: + same = np.full_like(den, np.nan, dtype=np.float64) + res = np.where(den == 0, same, num / den) + return np.nanmedian(res) + +@njit(cache=True, fastmath=False, nogil=True) +def medcouple_matrix_full(xs): + xs = np.asarray(xs) + xs = np.sort(xs) + n = len(xs) + if n % 2 == 0: + med = np.float64((xs[n // 2 - 1] + xs[n // 2]) / 2.0) + else: + med = np.float64(xs[(n - 1) // 2]) + zs = np.flip(xs - med) + zis = zs[zs <= 0] + zjs = zs[zs >= 0] + zjs = zjs[:, None] + num = zjs + zis + den = zjs - zis + if len(zs[zs == 0]): + z_index = np.flip(np.indices(np.shape(zs))[0]) + zis_idx = z_index[zs <= 0] + zjs_idx = z_index[zs >= 0] + zjs_idx = zjs_idx[:, None] + i_idx = np.ones_like(zjs_idx) * zis_idx + j_idx = np.ones_like(zis_idx) * zjs_idx + # Inline gen_matrix_same_full + same_ref = j_idx - i_idx + same = np.where(same_ref == 0, np.zeros_like(same_ref, dtype=np.float64), + np.where(same_ref > 0, np.ones_like(same_ref, dtype=np.float64), + np.full_like(same_ref, -1.0, dtype=np.float64))) + else: + same = np.full_like(den, np.nan, dtype=np.float64) + res = np.where(den == 0, same, num / den) + return np.nanmedian(res) + +@njit(cache=True, fastmath=False, nogil=True) +def medcouple_matrix_monotonic(xs): + xs = np.asarray(xs) + xs = np.sort(xs) + n = len(xs) + if n % 2 == 0: + med = np.float64((xs[n // 2 - 1] + xs[n // 2]) / 2.0) + else: + med = np.float64(xs[(n - 1) // 2]) + zs = np.flip(xs - med) + zis = zs[zs <= 0] + zjs = zs[zs >= 0] + zjs = zjs[:, None] + num = zjs + zis + den = zjs - zis + if len(zs[zs == 0]): + z_index = np.flip(np.indices(np.shape(zs))[0]) + zis_idx = z_index[zs <= 0] + zjs_idx = z_index[zs >= 0] + zjs_idx = zjs_idx[:, None] + i_idx = np.ones_like(zjs_idx) * zis_idx + j_idx = np.ones_like(zis_idx) * zjs_idx + # Inline gen_matrix_same_monotonic + k_idx = z_index[zs == 0] + max_k = np.max(k_idx) + min_k = np.min(k_idx) + same_ref = i_idx + j_idx - min_k - max_k + same = np.where(same_ref == 0, np.zeros_like(same_ref, dtype=np.float64), + np.where(same_ref > 0, np.ones_like(same_ref, dtype=np.float64), + np.full_like(same_ref, -1.0, dtype=np.float64))) + else: + same = np.full_like(den, np.nan, dtype=np.float64) + res = np.where(den == 0, same, num / den) + return np.nanmedian(res) + +# Cacheable Numba-compiled loop implementations (with literal kernel types) +@njit(cache=True, fastmath=False, nogil=True, parallel=True) +def _medcouple_loop_impl(xs, kernel_type_code, med_centred, start_index): + """ + Internal cacheable implementation with kernel_type as integer code: + 0 = partial, 1 = full, 2 = monotonic + """ + xs = np.asarray(xs) + xs_sorted = np.sort(xs) + n = len(xs_sorted) + if n % 2 == 0: + med = np.float64((xs_sorted[n // 2 - 1] + xs_sorted[n // 2]) / 2.0) + else: + med = np.float64(xs_sorted[(n - 1) // 2]) + + n_total = len(xs_sorted) + xis_vals = np.empty(n_total, dtype=np.float64) + xis_idx = np.empty(n_total, dtype=np.int64) + n_xis = 0 + for i in range(n_total): + x = xs_sorted[i] + if x <= med: + if med_centred: + xis_vals[n_xis] = x - med + else: + xis_vals[n_xis] = x + xis_idx[n_xis] = i + start_index + n_xis += 1 + xis_vals = xis_vals[:n_xis] + xis_idx = xis_idx[:n_xis] + + xjs_vals = np.empty(n_total, dtype=np.float64) + xjs_idx = np.empty(n_total, dtype=np.int64) + n_xjs = 0 + for j in range(n_total): + x = xs_sorted[j] + if x >= med: + if med_centred: + xjs_vals[n_xjs] = x - med + else: + xjs_vals[n_xjs] = x + xjs_idx[n_xjs] = j + start_index + n_xjs += 1 + xjs_vals = xjs_vals[:n_xjs] + xjs_idx = xjs_idx[:n_xjs] + + min_k = np.nan + max_k = np.nan + for k in range(n_total): + x = xs_sorted[k] + if x == med: + if np.isnan(min_k): + min_k = k + start_index + max_k = k + start_index + else: + max_k = k + start_index + + buffer = np.empty(n_xjs * n_xis, dtype=np.float64) + + for idx_j in prange(n_xjs): + j_pos = n_xjs - 1 - idx_j + xj = xjs_vals[j_pos] + j = xjs_idx[j_pos] + + for idx_i in range(n_xis): + i_pos = n_xis - 1 - idx_i + xi = xis_vals[i_pos] + i = xis_idx[i_pos] + + buffer_idx = idx_j * n_xis + idx_i + + if xi == xj: + # Inline kernel_same logic based on kernel_type_code + if kernel_type_code == 0: # partial + if i == j: + item = 0.0 + else: + item = np.nan + elif kernel_type_code == 1: # full + item = np.sign(j - i) + else: # 2 = monotonic + res = i + j - min_k - max_k + item = np.sign(res) + else: + # Calculate kernel_diff inline + if med_centred: + num = xj + xi + else: + num = (xj - med) - (med - xi) + den = xj - xi + item = num / den + buffer[buffer_idx] = item + + return np.nanmedian(buffer) + +# Python wrapper functions that dispatch to the cacheable Numba function +def medcouple_loop_partial_False_0(xs): + return _medcouple_loop_impl(xs, 0, False, 0) + +def medcouple_loop_partial_False_1(xs): + return _medcouple_loop_impl(xs, 0, False, 1) + +def medcouple_loop_partial_True_0(xs): + return _medcouple_loop_impl(xs, 0, True, 0) + +def medcouple_loop_partial_True_1(xs): + return _medcouple_loop_impl(xs, 0, True, 1) + +def medcouple_loop_full_False_0(xs): + return _medcouple_loop_impl(xs, 1, False, 0) + +def medcouple_loop_full_False_1(xs): + return _medcouple_loop_impl(xs, 1, False, 1) + +def medcouple_loop_full_True_0(xs): + return _medcouple_loop_impl(xs, 1, True, 0) + +def medcouple_loop_full_True_1(xs): + return _medcouple_loop_impl(xs, 1, True, 1) + +def medcouple_loop_monotonic_False_0(xs): + return _medcouple_loop_impl(xs, 2, False, 0) + +def medcouple_loop_monotonic_False_1(xs): + return _medcouple_loop_impl(xs, 2, False, 1) + +def medcouple_loop_monotonic_True_0(xs): + return _medcouple_loop_impl(xs, 2, True, 0) + +def medcouple_loop_monotonic_True_1(xs): + return _medcouple_loop_impl(xs, 2, True, 1) \ No newline at end of file diff --git a/python/test_df.py b/python/test_df.py new file mode 100644 index 0000000..2a62787 --- /dev/null +++ b/python/test_df.py @@ -0,0 +1,205 @@ +from timeit import timeit +from functools import partial +import pandas as pd +from tqdm import tqdm + +import numpy as np +from statsmodels.stats.stattools import medcouple as stms_medcouple + +from medcouple import ( + medcouple_naive_loop, + kernel_same_monotonic, + kernel_same_full, + kernel_same_partial, + medcouple_naive_matrix, + gen_matrix_same_monotonic, + gen_matrix_same_full, + gen_matrix_same_partial, +) + + +def build_test_dict() -> dict: + """ + Builds a dictionary of functions to test for calculating the median couple. + + The dictionary will contain the following functions: + - statsmodels_medcouple: the medcouple function from statsmodels + - naive_matrix_{gen_matrix_same.__name__}: the naive matrix based medcouple function using the given gen_matrix_same + - naive_loop_{kernel_same.__name__}_{med_centred}_{start_index}: the naive loop based medcouple function using + the given kernel_same, med_centred, and start_index + + Returns + ------- + dict + A dictionary of functions to test. + """ + algos = {} + algos["statsmodels_medcouple"] = lambda x: stms_medcouple(x).flat[0] + for gen_matrix_same in [ + gen_matrix_same_monotonic, + gen_matrix_same_full, + gen_matrix_same_partial, + ]: + algos[f"naive_matrix_{gen_matrix_same.__name__}"] = ( + lambda x: medcouple_naive_matrix(x, gen_matrix_same) + ) + + for kernel_same in [kernel_same_monotonic, kernel_same_full, kernel_same_partial]: + for med_centred in [False, True]: + for start_index in [0, 1]: + algos[ + f"naive_loop_{kernel_same.__name__}_{med_centred}_{start_index}" + ] = lambda x: medcouple_naive_loop( + x, kernel_same, med_centred, start_index + ) + + return algos + + +def run_tests( + algos: dict, arr: list, test_name: str, num_iterations: int = 10 +) -> pd.DataFrame: + """ + Runs a test of the medcouple implementations on the given array. + + The test will calculate the median of each array using each of the implementations + and then calculate the mean of the results and mean of the timings of each implementation. + The results will be returned in a DataFrame sorted by the mean timings. + + Parameters + ---------- + algos : dict + A dictionary of functions to test. + arr : list + The array to test on. + test_name : str + The name of the test. + num_iterations : int, default=10 + The number of iterations to run each test. + + Returns + ------- + pd.DataFrame + A DataFrame containing the results and timings of each test. + """ + results = {} + timings = {} + for name, f in tqdm(algos.items()): + results[name] = f(arr) + timings[name] = timeit(partial(f, arr), number=num_iterations) + df = pd.concat([pd.Series(results), pd.Series(timings)], axis=1) + df.columns = [f"results_{test_name}", f"timings_{test_name}"] + base_result = df.loc["statsmodels_medcouple", f"results_{test_name}"] + df[f"Passed_{test_name}"] = df[f"results_{test_name}"] == base_result + return df + + +def create_summary(algos_df: pd.DataFrame) -> pd.DataFrame: + """ + Creates a summary DataFrame from the results of the tests. + + The DataFrame will contain three columns: + - mean_results: the mean of the results of each test + - mean_timings: the mean of the timings of each test + - all_passed: a boolean indicating whether all tests passed for each implementation + + The DataFrame will be sorted by the mean timings in ascending order. + + Parameters + ---------- + algos_df : pd.DataFrame + A DataFrame containing the results and timings of each test. + + Returns + ------- + pd.DataFrame + A DataFrame containing the summary of the tests. + """ + mean_results = algos_df.filter(regex="results").mean(axis=1).rename("mean_results") + mean_timings = algos_df.filter(regex="timings").mean(axis=1).rename("mean_timings") + min_timings = algos_df.filter(regex="timings").min(axis=1).rename("min_timings") + max_timings = algos_df.filter(regex="timings").max(axis=1).rename("max_timings") + all_passed = algos_df.filter(regex="Passed").all(axis=1).rename("all_passed") + df = pd.concat( + [mean_results, mean_timings, min_timings, max_timings, all_passed], axis=1 + ) + df = df.sort_values("mean_timings", ascending=True) + return df + + +def test(test_arrs: list) -> None: + """ + Runs a test of the medcouple implementations on the given list of arrays. + + The test will calculate the median of each array using each of the implementations + and then calculate the mean of the results and mean, min, and max of the timings of each implementation. + The results will be printed in a DataFrame sorted by the mean timings. + + Parameters + ---------- + test_arrs : list + A list of arrays to test on. + + Returns + ------- + None + """ + algos = build_test_dict() + algos_df = pd.Series(algos).to_frame(name="algos") + print(f"Starting test on {len(test_arrs)} arrays") + for i, arr in enumerate(test_arrs, start=1): + name = f"arr_{i}" + print(f"Test array {i} of length {len(arr)}") + res_df = run_tests(algos, arr, name) + algos_df = pd.concat([algos_df, res_df], axis=1) + + algos_df.index.name = "Test case" + summary = create_summary(algos_df) + print("Summary of the results ordered by mean timings:", end="\n\n") + print(summary) + # Comparison of the two algorithms + print("\nGroup averages:") + loop_avg = summary[summary.index.str.startswith("naive_loop_")][ + "mean_timings" + ].mean() + matrix_avg = summary[summary.index.str.startswith("naive_matrix_")][ + "mean_timings" + ].mean() + statsmodels_avg = summary.loc["statsmodels_medcouple", "mean_timings"] + print(f" Statsmodels: {statsmodels_avg:.6f}s") + print(f" naive_loop average: {loop_avg:.6f}s") + print(f" naive_matrix average: {matrix_avg:.6f}s") + if matrix_avg > loop_avg: + print(f" Matrix is {matrix_avg / loop_avg:.2f}x slower than loop") + else: + print(f" Matrix is {loop_avg / matrix_avg:.2f}x faster than loop") + if loop_avg > statsmodels_avg: + print(f" Loop is {loop_avg / statsmodels_avg:.2f}x slower than statsmodels") + else: + print(f" Loop is {statsmodels_avg / loop_avg:.2f}x faster than statsmodels") + if matrix_avg > statsmodels_avg: + print(f" Matrix is {matrix_avg / statsmodels_avg:.2f}x slower than statsmodels") + else: + print(f" Matrix is {statsmodels_avg / matrix_avg:.2f}x faster than statsmodels") + + return None + + +def main(seed:int=1968) -> None: + np.random.seed(seed) + test_arrs = [ + [0, 1, 2, 2, 3], + [1, 2, 2, 2, 3, 4], + [1, 2, 2, 2, 2, 3, 4], + [0.2, 0.17, 0.08, 0.16, 0.88, 0.86, 0.09, 0.54, 0.27, 0.14], + [1] * 10 + [5] * 3, + [10.0] * 480 + + [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], + np.random.poisson(1, size=500), + np.random.rand(500), + ] + test(test_arrs) + + +if __name__ == "__main__": + main() diff --git a/python/test_numba.py b/python/test_numba.py new file mode 100644 index 0000000..d032133 --- /dev/null +++ b/python/test_numba.py @@ -0,0 +1,226 @@ +from timeit import timeit +from functools import partial +import pandas as pd +import statsmodels +from tqdm import tqdm + +import numpy as np +from statsmodels.stats.stattools import medcouple as stms_medcouple + +from medcouple_numba import ( + medcouple_matrix_partial, + medcouple_matrix_full, + medcouple_matrix_monotonic, + medcouple_loop_partial_False_0, + medcouple_loop_partial_False_1, + medcouple_loop_partial_True_0, + medcouple_loop_partial_True_1, + medcouple_loop_full_False_0, + medcouple_loop_full_False_1, + medcouple_loop_full_True_0, + medcouple_loop_full_True_1, + medcouple_loop_monotonic_False_0, + medcouple_loop_monotonic_False_1, + medcouple_loop_monotonic_True_0, + medcouple_loop_monotonic_True_1, +) + + +def build_test_dict() -> dict: + """ + Builds a dictionary of functions to test for calculating the median couple. + + The dictionary will contain the following functions: + - statsmodels_medcouple: the medcouple function from statsmodels + - naive_matrix_{gen_matrix_same.__name__}: the naive matrix based medcouple function using the given gen_matrix_same + - naive_loop_{kernel_same.__name__}_{med_centred}_{start_index}: the naive loop based medcouple function using + the given kernel_same, med_centred, and start_index + + Returns + ------- + dict + A dictionary of functions to test. + """ + algos = {} + algos["statsmodels_medcouple"] = lambda x: stms_medcouple(x).flat[0] + + # Matrix-based implementations + algos["naive_matrix_gen_matrix_same_monotonic"] = medcouple_matrix_monotonic + algos["naive_matrix_gen_matrix_same_full"] = medcouple_matrix_full + algos["naive_matrix_gen_matrix_same_partial"] = medcouple_matrix_partial + + # Loop-based implementations + algos["naive_loop_kernel_same_monotonic_False_0"] = medcouple_loop_monotonic_False_0 + algos["naive_loop_kernel_same_monotonic_False_1"] = medcouple_loop_monotonic_False_1 + algos["naive_loop_kernel_same_monotonic_True_0"] = medcouple_loop_monotonic_True_0 + algos["naive_loop_kernel_same_monotonic_True_1"] = medcouple_loop_monotonic_True_1 + algos["naive_loop_kernel_same_full_False_0"] = medcouple_loop_full_False_0 + algos["naive_loop_kernel_same_full_False_1"] = medcouple_loop_full_False_1 + algos["naive_loop_kernel_same_full_True_0"] = medcouple_loop_full_True_0 + algos["naive_loop_kernel_same_full_True_1"] = medcouple_loop_full_True_1 + algos["naive_loop_kernel_same_partial_False_0"] = medcouple_loop_partial_False_0 + algos["naive_loop_kernel_same_partial_False_1"] = medcouple_loop_partial_False_1 + algos["naive_loop_kernel_same_partial_True_0"] = medcouple_loop_partial_True_0 + algos["naive_loop_kernel_same_partial_True_1"] = medcouple_loop_partial_True_1 + + return algos + + +def run_tests( + algos: dict, arr: list, test_name: str, num_iterations: int = 10 +) -> pd.DataFrame: + """ + Runs a test of the medcouple implementations on the given array. + + The test will calculate the median of each array using each of the implementations + and then calculate the mean of the results and mean of the timings of each implementation. + The results will be returned in a DataFrame sorted by the mean timings. + + Parameters + ---------- + algos : dict + A dictionary of functions to test. + arr : list + The array to test on. + test_name : str + The name of the test. + num_iterations : int, default=10 + The number of iterations to run each test. + + Returns + ------- + pd.DataFrame + A DataFrame containing the results and timings of each test. + """ + results = {} + timings = {} + # Ensure array is numpy array with consistent dtype (float64) + if not isinstance(arr, np.ndarray): + work_arr = np.array(arr, dtype=np.float64) + else: + work_arr = arr.astype(np.float64) + + for name, f in tqdm(algos.items()): + results[name] = f(work_arr) + timings[name] = timeit(partial(f, work_arr), number=num_iterations) + df = pd.concat([pd.Series(results), pd.Series(timings)], axis=1) + df.columns = [f"results_{test_name}", f"timings_{test_name}"] + base_result = df.loc["statsmodels_medcouple", f"results_{test_name}"] + # Use np.isclose for floating-point comparison with tolerance + df[f"Passed_{test_name}"] = df[f"results_{test_name}"].apply( + lambda x: np.isclose(x, base_result, rtol=1e-9, atol=1e-12) + ) + return df + + +def create_summary(algos_df: pd.DataFrame) -> pd.DataFrame: + """ + Creates a summary DataFrame from the results of the tests. + + The DataFrame will contain five columns: + - mean_results: the mean of the results of each test + - mean_timings: the mean of the timings of each test + - min_timings: the min of the timings of each test + - max_timings: the max of the timings of each test + - all_passed: a boolean indicating whether all tests passed for each implementation + + The DataFrame will be sorted by the mean timings in ascending order. + + Parameters + ---------- + algos_df : pd.DataFrame + A DataFrame containing the results and timings of each test. + + Returns + ------- + pd.DataFrame + A DataFrame containing the summary of the tests. + """ + mean_results = algos_df.filter(regex="results").mean(axis=1).rename("mean_results") + mean_timings = algos_df.filter(regex="timings").mean(axis=1).rename("mean_timings") + min_timings = algos_df.filter(regex="timings").min(axis=1).rename("min_timings") + max_timings = algos_df.filter(regex="timings").max(axis=1).rename("max_timings") + all_passed = algos_df.filter(regex="Passed").all(axis=1).rename("all_passed") + df = pd.concat( + [mean_results, mean_timings, min_timings, max_timings, all_passed], axis=1 + ) + df = df.sort_values("mean_timings", ascending=True) + return df + + +def test(test_arrs: list) -> None: + """ + Runs a test of the medcouple implementations on the given list of arrays. + + The test will calculate the median of each array using each of the implementations + and then calculate the mean of the results and mean, min, and max of the timings of each implementation. + The results will be printed in a DataFrame sorted by the mean timings. + + Parameters + ---------- + test_arrs : list + A list of arrays to test on. + + Returns + ------- + None + """ + algos = build_test_dict() + algos_df = pd.Series(algos).to_frame(name="algos") + print(f"Starting test on {len(test_arrs)} arrays") + for i, arr in enumerate(test_arrs, start=1): + name = f"arr_{i}" + print(f"Test array {i} of length {len(arr)}") + res_df = run_tests(algos, arr, name) + algos_df = pd.concat([algos_df, res_df], axis=1) + + algos_df.index.name = "Test case" + summary = create_summary(algos_df) + print("Summary of the results ordered by mean timings:", end="\n\n") + print(summary) + # Comparison of the two algorithms + print("\nGroup averages:") + loop_avg = summary[summary.index.str.startswith("naive_loop_")][ + "mean_timings" + ].mean() + matrix_avg = summary[summary.index.str.startswith("naive_matrix_")][ + "mean_timings" + ].mean() + statsmodels_avg = summary.loc["statsmodels_medcouple", "mean_timings"] + print(f" Statsmodels: {statsmodels_avg:.6f}s") + print(f" naive_loop average: {loop_avg:.6f}s") + print(f" naive_matrix average: {matrix_avg:.6f}s") + if matrix_avg > loop_avg: + print(f" Matrix is {matrix_avg / loop_avg:.2f}x slower than loop") + else: + print(f" Matrix is {loop_avg / matrix_avg:.2f}x faster than loop") + if loop_avg > statsmodels_avg: + print(f" Loop is {loop_avg / statsmodels_avg:.2f}x slower than statsmodels") + else: + print(f" Loop is {statsmodels_avg / loop_avg:.2f}x faster than statsmodels") + if matrix_avg > statsmodels_avg: + print(f" Matrix is {matrix_avg / statsmodels_avg:.2f}x slower than statsmodels") + else: + print(f" Matrix is {statsmodels_avg / matrix_avg:.2f}x faster than statsmodels") + + return None + + +def main(seed:int=1968) -> None: + np.random.seed(seed) + test_arrs = [ + [0, 1, 2, 2, 3], + [1, 2, 2, 2, 3, 4], + [1, 2, 2, 2, 2, 3, 4], + [0.2, 0.17, 0.08, 0.16, 0.88, 0.86, 0.09, 0.54, 0.27, 0.14], + [1] * 10 + [5] * 3, + [10.0] * 480 + + [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], + np.random.poisson(1, size=500), + np.random.rand(500), + ] + test(test_arrs) + + +if __name__ == "__main__": + main()