polished plots

src/timing_profile.py

+# %%
 import timeit
 from my_ext import array_sum, array_sum_nocopy, array_sum_pyobj  # noqa
 import numpy as np
@@ -5,25 +6,39 @@ import matplotlib.pyplot as plt
 from seaborn import despine
 
 
-# %%
-M = 500
-N = np.logspace(1, 5, 10).astype(int)
+# %% Utils
+def bs_mean_ci(x, alpha=0.95, n_bootstrap=1000):
+    """Bootstrap mean and confidence interval."""
+    n, _ = x.shape
+    mu = np.zeros(n)
+    ci = np.zeros((n, 2))
+    for k, xi in enumerate(x):
+        bs = np.random.choice(xi, size=(n_bootstrap,), replace=True)
+        mu[k] = bs.mean()
+        ci[k] = np.percentile(bs, [(1 - alpha) * 100 / 2, 100 - (1 - alpha) * 100 / 2])
+    return mu, ci
+
+
+# %% Perform timing experiments
+R = 20  # Number of repeated timing experiments
+M = 500  # Number of iterations in each timing measurement
+N = np.logspace(1, 5, 10).astype(int)  # Sizes to test
 
 res = {"np": [], "array_sum": [], "array_sum_nocopy": [], "array_sum_pyobj": []}
 for Ni in N:
     print(f"Ni = {Ni}")
     x = np.random.normal(0, 1, int(Ni))
-    ri = timeit.repeat("x.sum()", globals=globals(), number=M, repeat=5)
-    res["np"].append(np.median(ri) / M)
+    ri = np.array(timeit.repeat("x.sum()", globals=globals(), number=M, repeat=R)) / M
+    res["np"].append(ri)
 
-    ri = timeit.repeat("array_sum(x)", globals=globals(), number=M, repeat=5)
-    res["array_sum"].append(np.median(ri) / M)
+    ri = np.array(timeit.repeat("array_sum(x)", globals=globals(), number=M, repeat=R)) / M
+    res["array_sum"].append(ri)
 
-    ri = timeit.repeat("array_sum_nocopy(x)", globals=globals(), number=M, repeat=5)
-    res["array_sum_nocopy"].append(np.median(ri) / M)
+    ri = np.array(timeit.repeat("array_sum_nocopy(x)", globals=globals(), number=M, repeat=R)) / M
+    res["array_sum_nocopy"].append(ri)
 
-    ri = timeit.repeat("array_sum_pyobj(x)", globals=globals(), number=M, repeat=5)
-    res["array_sum_pyobj"].append(np.median(ri) / M)
+    ri = np.array(timeit.repeat("array_sum_pyobj(x)", globals=globals(), number=M, repeat=R)) / M
+    res["array_sum_pyobj"].append(ri)
 
 for k in res:
     res[k] = np.array(res[k])
@@ -31,10 +46,10 @@ for k in res:
 
 # %% Plot results
 fig, ax = plt.subplots(num=10, clear=True)
-ax.loglog(N, res["np"] * 1e9, label="np")
-ax.loglog(N, res["array_sum"] * 1e9, label="array_sum")
-ax.loglog(N, res["array_sum_nocopy"] * 1e9, label="array_sum_nocopy")
-ax.loglog(N, res["array_sum_pyobj"] * 1e9, label="array_sum_pyobj")
+for k in res:
+    mu, ci = bs_mean_ci(res[k])
+    ax.loglog(N, mu, label=k)
+    ax.fill_between(N, ci[:, 0], ci[:, 1], alpha=0.2)
 ax.legend(frameon=False)
 ax.set_xlabel("N")
 ax.set_ylabel("Time (ns)")