From 68db0134f33d7a614f277525b7604f4d8e7f58b7 Mon Sep 17 00:00:00 2001 From: adalseno Date: Wed, 5 Nov 2025 00:23:17 +0100 Subject: [PATCH 1/2] Add Numba-accelerated implementation with parallel execution - Add medcouple_numba.py: Numba JIT-compiled version with parallel support - Add test_numba.py: Comprehensive test suite comparing implementations - Add test_df.py: DataFrame-based testing utilities - Achieves ~4x speedup over matrix implementations - All implementations produce identical results (within floating-point tolerance) --- python/README.md | 48 +++++ python/medcouple_numba.py | 421 ++++++++++++++++++++++++++++++++++++++ python/test_df.py | 205 +++++++++++++++++++ python/test_numba.py | 226 ++++++++++++++++++++ 4 files changed, 900 insertions(+) create mode 100644 python/README.md create mode 100644 python/medcouple_numba.py create mode 100644 python/test_df.py create mode 100644 python/test_numba.py diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..d4b08fd --- /dev/null +++ b/python/README.md @@ -0,0 +1,48 @@ +### General + +Calculates medcouple for Python 3 and Scala (native or Spark) using the naive algorithm (does not scale well to large n). + +Implemented based on +- https://wis.kuleuven.be/stat/robust/papers/2004/medcouple.pdf +- https://en.wikipedia.org/wiki/Medcouple +- https://github.com/statsmodels/statsmodels/issues/5395 + +### Python + +Recommend using [statsmodels](https://github.com/statsmodels/statsmodels) instead for small n, but you can use this code for reference to understand the logic of the medcouple algorithm. See the bottom of `python/test.py` for an example of how it's used. Tests show that the results match up but statsmodels is slightly faster. + +As for scalable libraries using the fast algorithm, haven't found any Python libraries that are 100% accurate yet, existing ones have a problem with handling ties to the median like in the example in the GitHub issue above, and the results do not match those of the naive algorithms and R's fast algorithm. + +### Scala + +Can use as a library project/jar, or just copy the relevant files over to your project. + +To run without Spark: +``` +cd scala +# ensure `scalaVersion` in `build.sbt` matches your version or you might get a runtime error like: Exception in thread "main" java.lang.NoSuchMethodError: scala.Predef$.doubleArrayOps([D)[D +sbt assembly +scala -classpath target/medcouple.jar stats.medcouple.TestAlgos +``` +To run the Spark code with UDF, do the above but instead of the last line, spark-submit that jar with class `stats.medcouple.TestSpark`. + +#### Advanced, safe to ignore + +If you see these lines during runtime, it means that the Medcouple implementations using Breeze matrices will be unoptimised and slower than it should be: +``` +WARN BLAS: Failed to load implementation from: com.github.fommil.netlib.NativeSystemBLAS +WARN BLAS: Failed to load implementation from: com.github.fommil.netlib.NativeRefBLAS +``` +To fix that, ensure these libraries are included in `build.sbt`: +``` + "org.scalanlp" %% "breeze-natives" % "1.0", + "com.github.fommil.netlib" % "all" % "1.1.2" pomOnly() +``` +Shade so they don't conflict with spark's version of breeze: +``` +assemblyShadeRules in assembly := Seq( + ShadeRule.rename("breeze.**" -> "shaded.breeze.@1").inAll +) +``` +And get libraries like BLAS natively installed on your system. +If you manage to do this, you can check if `Medcouple.calcMatrix` is faster. Otherwise, just stick with the default `Medcouple.calcLoop`. 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() From 057b25673fecbc91217fab5cb632dce9b4cfce06 Mon Sep 17 00:00:00 2001 From: adalseno Date: Thu, 6 Nov 2025 01:34:50 +0100 Subject: [PATCH 2/2] fIXed README.md in python folder --- python/README.md | 119 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 81 insertions(+), 38 deletions(-) diff --git a/python/README.md b/python/README.md index d4b08fd..c6763e4 100644 --- a/python/README.md +++ b/python/README.md @@ -1,48 +1,91 @@ -### General +## Numba Implementation -Calculates medcouple for Python 3 and Scala (native or Spark) using the naive algorithm (does not scale well to large n). +A high-performance Numba JIT-compiled implementation is available in `medcouple_numba.py`. -Implemented based on -- https://wis.kuleuven.be/stat/robust/papers/2004/medcouple.pdf -- https://en.wikipedia.org/wiki/Medcouple -- https://github.com/statsmodels/statsmodels/issues/5395 +### Installation +```bash +pip install numba +``` +Example usage: +```python +from medcouple_numba import medcouple_matrix_full +result = medcouple_matrix_full(data) +``` +### Performance: -### Python +The Numba implementation with parallel execution achieves sensitive speedup compared to the pure NumPy implementation. -Recommend using [statsmodels](https://github.com/statsmodels/statsmodels) instead for small n, but you can use this code for reference to understand the logic of the medcouple algorithm. See the bottom of `python/test.py` for an example of how it's used. Tests show that the results match up but statsmodels is slightly faster. +I have also added more descriptive tests, `test_df.py` and `test_numba.py`, that require pandas and tqdm. -As for scalable libraries using the fast algorithm, haven't found any Python libraries that are 100% accurate yet, existing ones have a problem with handling ties to the median like in the example in the GitHub issue above, and the results do not match those of the naive algorithms and R's fast algorithm. +If you want to print the table in markdown format, you need to install `tabulate` and `print(summary.to_markdown())`. -### Scala +#### Numba implementation: -Can use as a library project/jar, or just copy the relevant files over to your project. +**Summary of the results ordered by mean timings:** -To run without Spark: -``` -cd scala -# ensure `scalaVersion` in `build.sbt` matches your version or you might get a runtime error like: Exception in thread "main" java.lang.NoSuchMethodError: scala.Predef$.doubleArrayOps([D)[D -sbt assembly -scala -classpath target/medcouple.jar stats.medcouple.TestAlgos -``` -To run the Spark code with UDF, do the above but instead of the last line, spark-submit that jar with class `stats.medcouple.TestSpark`. +| 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 | -#### Advanced, safe to ignore +Group averages: -If you see these lines during runtime, it means that the Medcouple implementations using Breeze matrices will be unoptimised and slower than it should be: -``` -WARN BLAS: Failed to load implementation from: com.github.fommil.netlib.NativeSystemBLAS -WARN BLAS: Failed to load implementation from: com.github.fommil.netlib.NativeRefBLAS -``` -To fix that, ensure these libraries are included in `build.sbt`: -``` - "org.scalanlp" %% "breeze-natives" % "1.0", - "com.github.fommil.netlib" % "all" % "1.1.2" pomOnly() -``` -Shade so they don't conflict with spark's version of breeze: -``` -assemblyShadeRules in assembly := Seq( - ShadeRule.rename("breeze.**" -> "shaded.breeze.@1").inAll -) -``` -And get libraries like BLAS natively installed on your system. -If you manage to do this, you can check if `Medcouple.calcMatrix` is faster. Otherwise, just stick with the default `Medcouple.calcLoop`. + 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