Source
scripts/numpy_compare.py
1
#!/usr/bin/env python32
"""HONEST cheatah `linalg` vs NumPy comparison.4
NumPy's array ops dispatch to BLAS/LAPACK (blocked, vectorized, often multi-threaded5
Fortran), so this is the hard comparison — and the honest one. We do NOT tilt it in6
cheatah's favor: the SAME matrix (a fixed-seed, well-conditioned one) is fed to both7
libraries, each runs the SAME operation N times with the result consumed, and we8
print both answers so you can see they agree. Expect:10
* small matrices → cheatah often WINS (no Python call + NumPy dispatch overhead);11
* large matrices → NumPy WINS (BLAS is a hand-tuned, threaded kernel — we don't12
pretend otherwise);13
* the crossover is the interesting part.15
The cheatah side compiles a real .purr that loops the op (linalg routines live in a16
separate .a, so the optimizer can't hoist or delete the opaque call — verified by an17
empty-loop guard); the NumPy side loops the op in CPython. Build `release` first.19
python3 scripts/numpy_compare.py20
"""21
import atexit22
import glob23
import os24
import platform25
import re26
import shutil27
import subprocess28
import sys29
import tempfile30
import time32
import numpy as np34
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))35
PURRC = os.path.join(ROOT, "build", "release", "bin", "purrc")36
CHEATAH = os.path.join(ROOT, "build", "release", "bin", "cheatah")37
# 7 striated rounds replacing 5 best-of-N trials per side. See scripts/perf_suite.py for38
# the reasoning: a minimum has no dispersion, and taking each side's minimum independently39
# pairs one side's luckiest run against the other's.40
ROUNDS = 741
rng = np.random.default_rng(0)44
def lit(M):45
"""A cheatah list literal of M's row-major elements (full float precision)."""46
return "[" + ", ".join(repr(float(v)) for v in np.asarray(M).flatten()) + "]"49
def spd(n):50
"""A symmetric positive-definite, well-conditioned n×n matrix (good for solve/inv/51
det/eig). Same matrix used on both sides."""52
A = rng.standard_normal((n, n))53
return A @ A.T + n * np.eye(n)56
def gen(n):57
"""A general (non-symmetric) well-conditioned n×n matrix — for the general58
eigenvalues and matrix_power (which must not route through the symmetric path)."""59
return rng.standard_normal((n, n)) * 0.3 + np.eye(n)62
def run(argv):63
return subprocess.run(argv, capture_output=True, text=True)66
def parse(stdout):67
lines = [l for l in stdout.strip().splitlines() if l.strip()]68
if len(lines) < 2:69
return None, None70
try:71
return float(lines[-1]), lines[-2]72
except ValueError:73
return None, None76
_TEMPDIRS = []77
atexit.register(lambda: [shutil.rmtree(d, ignore_errors=True) for d in _TEMPDIRS])80
def build_cheatah(setup, expr, iters, extract="[0]"):81
"""setup builds the operands; `expr` is the op; we accumulate a scalar from it82
(the bare scalar when extract=='', else element extract like [0] or [0, 0]) so the83
opaque library call is genuinely consumed and can't be optimized away."""84
consume = expr if extract == "" else f"ndarray.get({expr}, {extract})"85
body = f"acc = acc + {consume}"86
src = (87
"import io\nimport time\nimport ndarray\nimport linalg\n" + setup + "\n"88
"let acc = 0.0\n"89
"let t0 = time.monotonic()\n"90
f"for i in range(0, {iters}) {{\n {body}\n}}\n"91
"let t1 = time.monotonic()\n"92
"io.print(acc)\nio.print(t1 - t0)\n"93
)94
# mkdtemp, not TemporaryDirectory: the .so must outlive this call so the rounds below can95
# re-run it without recompiling. Registered for cleanup at exit rather than leaked.96
d = tempfile.mkdtemp(prefix="numpy_compare.")97
_TEMPDIRS.append(d)98
purr, so = os.path.join(d, "b.purr"), os.path.join(d, "b.so")99
open(purr, "w").write(src)100
if run([PURRC, purr, "-o", so]).returncode != 0:101
return None102
return so105
def _median(xs):106
ys = sorted(xs)107
n = len(ys)108
return ys[n // 2] if n % 2 else 0.5 * (ys[n // 2 - 1] + ys[n // 2])111
def np_once(operands, op, iters):112
"""One timed NumPy pass. op(operands, i) -> scalar contribution; `i` lets a case vary113
its input per iteration (to match a cheatah loop that does the same)."""114
import time as _t115
a = 0.0116
t0 = _t.monotonic()117
for i in range(iters):118
a += op(operands, i)119
return _t.monotonic() - t0, repr(a)122
ROWS = [] # (label, n, cheatah_us, numpy_us, ratio, lo, hi) for the generated table125
def bench(label, n, iters, ch_setup, ch_expr, np_build, np_op, extract="[0]"):126
# Compile ONCE, outside the rounds: compilation is not the thing being measured, and127
# re-running purrc between rounds would put seconds between the two sides of a pair.128
so = build_cheatah(ch_setup, ch_expr, iters, extract)129
if so is None:130
print(f"{label:<22}{n:>5} (cheatah compile/run failed)")131
return132
operands = np_build()134
# ROUNDS striated rounds: cheatah then NumPy, adjacently, before either repeats. The135
# headline is the median of the per-round PAIRED ratios — not the ratio of two medians,136
# which the two only agree on when the machine holds perfectly still.137
ch_ts, np_ts, ratios = [], [], []138
c_acc = n_acc = None139
for _ in range(ROUNDS):140
ct, ca = parse(run([CHEATAH, so]).stdout)141
if ct is None:142
print(f"{label:<22}{n:>5} (cheatah run failed)")143
return144
nt, na = np_once(operands, np_op, iters)145
ch_ts.append(ct)146
np_ts.append(nt)147
ratios.append(nt / ct)148
c_acc, n_acc = ca, na150
ct, nt = _median(ch_ts), _median(np_ts)151
cu, nu = ct / iters * 1e6, nt / iters * 1e6 # µs per op152
r = _median(ratios) # >1 = cheatah faster153
winner = "cheatah" if r > 1.0 else "numpy"154
ratio = r if r > 1.0 else 1.0 / r155
lo, hi = min(ratios), max(ratios)156
agree = ""157
try:158
if abs(float(c_acc) - float(n_acc)) > 1e-3 * max(1.0, abs(float(n_acc))):159
agree = f" ⚠ disagree (cheatah {c_acc} vs numpy {n_acc})"160
except ValueError:161
pass162
# The raw band is numpy/cheatah in both directions, so a reader can see the swing without163
# it being folded through the winner flip.164
band = f" [{lo:.2f}–{hi:.2f} raw]"165
ROWS.append((label, n, cu, nu, r, lo, hi))166
print(f"{label:<22}{n:>5}{cu:>11.2f}{nu:>11.2f} {winner:>7} {ratio:>5.1f}×{band}{agree}")170
def _capture(cmd):171
try:172
return subprocess.run(cmd, shell=True, capture_output=True, text=True).stdout.strip()173
except Exception:174
return ""177
def _blas():178
"""WHICH BLAS, and how many threads. Where the crossovers land depends far more on the179
linked BLAS (reference vs OpenBLAS vs MKL) and its thread count than on the NumPy version,180
so a stamp naming only the version cannot be reproduced from.182
np.show_config() is consulted first but on a system BLAS it honestly answers183
"blas / unknown", which identifies nothing. In that case resolve the shared object NumPy184
actually links, which does."""185
name = ""186
try:187
blas = np.show_config(mode="dicts").get("Build Dependencies", {}).get("blas", {})188
n, v = blas.get("name", ""), blas.get("version", "")189
if n and n != "blas" and v and v != "unknown":190
name = f"{n} {v}"191
except Exception:192
pass193
if not name:194
# ldd the extension module that carries the BLAS dependency. The subdirectory moved195
# between NumPy versions (core/ -> _core/), so glob rather than hardcode it.196
cands = glob.glob(os.path.join(os.path.dirname(np.__file__), "**",197
"_multiarray_umath*.so"), recursive=True)198
if cands:199
libs = _capture(f"ldd {cands[0]} 2>/dev/null | grep -iE 'blas|mkl' | head -2")200
resolved = []201
for line in libs.splitlines():202
soname = line.split("=>")[0].strip()203
path = line.split("=>")[1].split("(")[0].strip() if "=>" in line else ""204
# On Debian libblas.so.3 is an alternatives symlink; the target is the thing205
# that actually determines the numbers, so report what it points at.206
real = _capture(f"readlink -f {path}") if path else ""207
resolved.append(f"{soname} -> {os.path.basename(real)}" if real else soname)208
if resolved:209
name = "; ".join(resolved)210
if not name:211
name = "unidentified system BLAS"212
threads = (os.environ.get("OPENBLAS_NUM_THREADS") or os.environ.get("OMP_NUM_THREADS")213
or "unset (BLAS default)")214
return f"{name}, threads={threads}"217
def write_md(path, suite_name):218
"""Emit the generated region body scripts/bench_table.purr expects.220
ONE HARNESS, ONE TABLE. The linalg README used to place a cheatah column measured here221
beside an Eigen column measured by Google Benchmark, with a prose warning not to read222
across. A warning is a worse fix than a structure: this table carries only the columns223
this harness measured, and the Eigen comparison is its own generated table from its own224
harness (docs/bench/linalg-vs-eigen.md)."""225
watch = ("stdlib/linalg/, stdlib/ndarray/, scripts/numpy_compare.py"226
if suite_name.startswith("linalg")227
else "stdlib/ndarray/, scripts/numpy_compare.py")228
commit = _capture("git rev-parse --short HEAD") or "unknown"229
# Judged against THIS suite's watched sources only, with Markdown excluded — see the long230
# note in scripts/bench_run.sh. A whole-tree check means any edit anywhere invalidates231
# every later measurement, which turns one prose fix into a full re-measure.232
subprocess.run("git update-index --refresh", shell=True, capture_output=True)233
spec = " ".join(f"'{p.strip()}'" for p in watch.split(",") if p.strip())234
if spec and subprocess.run(f"git diff --quiet -- {spec} ':!*.md' ':!docs/bench'",235
shell=True).returncode != 0:236
commit += " (dirty)"237
host = _capture("awk -F': ' '/^model name/{print $2; exit}' /proc/cpuinfo") or platform.machine()238
gov = _capture("cat /sys/devices/system/cpu/cpu0/cpufreq/scaling_governor")239
if gov:240
host += f" (governor={gov})"241
lines = [242
"<!-- cheatah-bench-stamp v1",243
f" suite: {suite_name}",244
f" generated: {time.strftime('%Y-%m-%d')}",245
f" commit: {commit}",246
f" host: {host}, {os.cpu_count()} CPUs",247
" cpu-scaling: enabled",248
" build: purrc -> -O3 -march=native",249
f" competitors: NumPy {np.__version__} on {_blas()}, CPython "250
f"{sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}",251
f" harness: rounds={ROUNDS}, striated (cheatah and NumPy adjacent in each round)",252
" statistic: median of per-round PAIRED ratios; [lo-hi] is the range of those",253
f" watch: {watch}",254
" publishable: true",255
"",256
" PRODUCED BY:",257
f" python3 scripts/numpy_compare.py --suite "258
f"{suite_name.replace('-vs-numpy', '')} --md {path}",259
"-->",260
"",261
"| op | operand dimensions | cheatah | NumPy | winner | band |",262
"|----|--------------------|--------:|------:|--------|------|",263
]264
for label, n, cu, nu, r, lo, hi in ROWS:265
win = f"**cheatah {r:.1f}x**" if r > 1.0 else f"NumPy {1 / r:.1f}x"266
lines.append(f"| `{label}` | {n} | {cu:.2f} | {nu:.2f} | {win} | {lo:.2f}-{hi:.2f} |")267
with open(path, "w") as f:268
f.write("\n".join(lines) + "\n")269
print(f"\nwrote {path}")272
def main():273
if not (os.path.exists(PURRC) and os.path.exists(CHEATAH)):274
sys.exit("numpy_compare: build the `release` preset first (need purrc + cheatah).")275
suite = "linalg"276
if "--suite" in sys.argv:277
i = sys.argv.index("--suite")278
if i + 1 >= len(sys.argv):279
sys.exit("numpy_compare: --suite needs a name (linalg|ndarray)")280
suite = sys.argv[i + 1]281
if suite not in ("linalg", "ndarray"):282
sys.exit(f"numpy_compare: unknown suite {suite!r} — expected linalg or ndarray")283
md_out = None284
if "--md" in sys.argv:285
i = sys.argv.index("--md")286
if i + 1 >= len(sys.argv):287
sys.exit("numpy_compare: --md needs an output path")288
md_out = sys.argv[i + 1]289
print(f"# cheatah linalg vs NumPy {np.__version__} (BLAS/LAPACK) — µs per op, "290
f"same matrix, result consumed\n")291
print(f"{'operation':<22}{'n':>5}{'cheatah':>11}{'numpy':>11} {'winner':>13}")292
print("-" * 70)294
if suite == "ndarray":295
run_ndarray()296
else:297
run_linalg()299
print("\nNote: NumPy calls BLAS/LAPACK (often threaded). cheatah's kernels are "300
"single-threaded\nauto-vectorized C++. Small n favors cheatah (no Python/"301
"dispatch overhead); large n favors BLAS / vectorized ufuncs.")303
if md_out is not None:304
write_md(md_out, f"{suite}-vs-numpy")307
def run_linalg():308
# ---- matmul A·B (both n×n) ----309
for n, iters in [(4, 200000), (16, 50000), (32, 20000), (64, 5000), (96, 2000)]:310
A, B = spd(n), spd(n)311
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])\n" \312
f"let B = ndarray.reshape(ndarray.array({lit(B)}), [{n}, {n}])"313
bench("matmul", n, iters, setup, "linalg.matmul(A, B)",314
lambda A=A, B=B: (A, B), lambda o, i: float((o[0] @ o[1])[0, 0]),315
extract="[0, 0]")317
# ---- solve A·x = b ----318
for n, iters in [(4, 200000), (16, 50000), (32, 20000), (64, 4000)]:319
A, b = spd(n), rng.standard_normal(n)320
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])\n" \321
f"let b = ndarray.array({lit(b)})"322
bench("solve", n, iters, setup, "linalg.solve(A, b)",323
lambda A=A, b=b: (A, b), lambda o, i: float(np.linalg.solve(o[0], o[1])[0]))325
# ---- det(A) ----326
for n, iters in [(4, 200000), (16, 50000), (32, 20000), (64, 5000)]:327
A = spd(n)328
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"329
bench("det", n, iters, setup, "linalg.det(A)",330
lambda A=A: (A,), lambda o, i: float(np.linalg.det(o[0])), extract="")332
# ---- inv(A) ----333
for n, iters in [(4, 100000), (16, 30000), (32, 10000), (64, 3000)]:334
A = spd(n)335
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"336
bench("inv", n, iters, setup, "linalg.inv(A)",337
lambda A=A: (A,), lambda o, i: float(np.linalg.inv(o[0])[0, 0]),338
extract="[0, 0]")340
# ---- eigvalsh(A) (symmetric eigenvalues) ----341
# Small n is the physicist's common case (few-level systems, spin Hamiltonians,342
# parameter sweeps that solve the same-shape problem millions of times) — measure343
# finely there to find where cheatah's no-overhead Jacobi beats LAPACK dispatch.344
for n, iters in [(2, 200000), (3, 200000), (4, 100000), (6, 80000), (8, 60000),345
(16, 20000), (32, 5000), (64, 1000)]:346
A = spd(n)347
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"348
# cheatah returns eigenvalues DESCENDING ([0] = largest); NumPy returns them349
# ASCENDING ([-1] = largest) — read the same end so the cross-check agrees.350
bench("eigvalsh", n, iters, setup, "linalg.eigvalsh(A)",351
lambda A=A: (A,), lambda o, i: float(np.linalg.eigvalsh(o[0])[-1]))353
# ---- dot (vectors) ----354
for n, iters in [(64, 200000), (1024, 50000), (16384, 5000)]:355
u, v = rng.standard_normal(n), rng.standard_normal(n)356
setup = f"let u = ndarray.array({lit(u)})\nlet v = ndarray.array({lit(v)})"357
bench("dot", n, iters, setup, "linalg.dot(u, v)",358
lambda u=u, v=v: (u, v), lambda o, i: float(o[0] @ o[1]), extract="")360
# ---- element-wise math ufuncs (cheatah ndarray.sqrt etc. vs numpy np.sqrt) ----361
# cheatah's ndarray ufuncs are header templates (inlined), so a loop-invariant362
# input would be hoisted out; we add a tiny i-dependent scalar each iteration to363
# force the work to run. NumPy does the same add+ufunc (Python never hoists), so364
# the comparison stays fair — both compute "add a scalar, then the ufunc, over n".365
for fn in ["sqrt", "exp", "sin"]:366
npfn = getattr(np, fn)367
for n, iters in [(64, 200000), (1024, 50000), (16384, 4000)]:368
v = rng.random(n) + 0.1 # positive (valid for sqrt)369
setup = f"let v = ndarray.array({lit(v)})"370
expr = f"ndarray.{fn}(ndarray.add(v, ndarray.scalar(0.000001 * i)))"371
bench(f"ndarray.{fn}", n, iters, setup, expr,372
lambda v=v: (v,),373
lambda o, i, f=npfn: float(f(o[0] + 0.000001 * i)[0]), extract="[0]")375
# ---- element-wise add of a broadcast scalar — a memory-bandwidth-bound op (one read,376
# one write). This is where allocating the result buffer UNINITIALIZED instead of377
# zero-filling it before the overwrite (see ndarray::buffer_t) matters most. ----378
for n, iters in [(64, 200000), (16384, 4000)]:379
v = rng.random(n) + 0.1380
setup = f"let v = ndarray.array({lit(v)})"381
expr = "ndarray.add(v, ndarray.scalar(0.000001 * i))"382
bench("ndarray.add", n, iters, setup, expr,383
lambda v=v: (v,),384
lambda o, i: float((o[0] + 0.000001 * i)[0]), extract="[0]")386
# ---- Cholesky factor (SPD) ----387
for n, iters in [(8, 50000), (32, 8000), (64, 2000)]:388
A = spd(n)389
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"390
bench("cholesky", n, iters, setup, "linalg.cholesky(A)",391
lambda A=A: (A,), lambda o, i: float(np.linalg.cholesky(o[0])[0, 0]),392
extract="[0, 0]")394
# ---- QR (R[0,0]) ----395
for n, iters in [(8, 40000), (32, 6000), (64, 1500)]:396
A = spd(n)397
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"398
bench("qr", n, iters, setup, "linalg.qr(A).r",399
lambda A=A: (A,), lambda o, i: float(np.linalg.qr(o[0])[1][0, 0]),400
extract="[0, 0]")402
# ---- SVD — compared FAIRLY on both sides ----403
# values-only: cheatah `svdvals` vs numpy `svd(compute_uv=False)`404
for n, iters in [(8, 20000), (32, 3000), (64, 800)]:405
A = spd(n)406
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"407
bench("svdvals", n, iters, setup, "linalg.svdvals(A)",408
lambda A=A: (A,), lambda o, i: float(np.linalg.svd(o[0], compute_uv=False)[0]),409
extract="[0]")410
# full decomposition: cheatah `svd` vs numpy `svd` (both compute U and Vᵀ)411
for n, iters in [(8, 20000), (32, 2000), (64, 500)]:412
A = spd(n)413
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"414
bench("svd (full)", n, iters, setup, "linalg.svd(A).s",415
lambda A=A: (A,), lambda o, i: float(np.linalg.svd(o[0])[1][0]), extract="[0]")417
# ---- pseudo-inverse (SVD-based) ----418
for n, iters in [(8, 20000), (32, 3000), (64, 800)]:419
A = spd(n)420
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"421
bench("pinv", n, iters, setup, "linalg.pinv(A)",422
lambda A=A: (A,), lambda o, i: float(np.linalg.pinv(o[0])[0, 0]),423
extract="[0, 0]")425
# ---- condition number / matrix_rank (SVD-based) ----426
for n, iters in [(8, 20000), (32, 3000), (64, 800)]:427
A = spd(n)428
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"429
bench("cond", n, iters, setup, "linalg.cond(A)",430
lambda A=A: (A,), lambda o, i: float(np.linalg.cond(o[0])), extract="")431
for n, iters in [(8, 20000), (32, 3000), (64, 800)]:432
A = spd(n)433
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"434
bench("matrix_rank", n, iters, setup, "linalg.matrix_rank(A)",435
lambda A=A: (A,), lambda o, i: float(np.linalg.matrix_rank(o[0])), extract="")437
# ---- slogdet (log|det|, LU-based) ----438
for n, iters in [(8, 80000), (32, 15000), (64, 4000)]:439
A = spd(n)440
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"441
bench("slogdet", n, iters, setup, "linalg.slogdet(A).logabsdet",442
lambda A=A: (A,), lambda o, i: float(np.linalg.slogdet(o[0])[1]), extract="")444
# ---- eigh (symmetric, eigenvalues AND eigenvectors) ----445
for n, iters in [(8, 40000), (32, 4000), (64, 1000)]:446
A = spd(n)447
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"448
bench("eigh", n, iters, setup, "linalg.eigh(A).values",449
lambda A=A: (A,), lambda o, i: float(np.linalg.eigh(o[0])[0][-1]), extract="[0]")451
# ---- eigvals (GENERAL, non-symmetric -> Hessenberg + shifted QR) ----452
for n, iters in [(8, 8000), (16, 2000), (32, 400)]:453
G = gen(n)454
setup = f"let A = ndarray.reshape(ndarray.array({lit(G)}), [{n}, {n}])"455
bench("eigvals", n, iters, setup, "ndarray.real(linalg.eigvals(A))",456
lambda G=G: (G,), lambda o, i: float(np.sort(np.linalg.eigvals(o[0]).real)[-1]),457
extract="[0]")459
# ---- matrix_power A³ (general) ----460
for n, iters in [(8, 40000), (32, 5000), (64, 1500)]:461
G = gen(n)462
setup = f"let A = ndarray.reshape(ndarray.array({lit(G)}), [{n}, {n}])"463
bench("matrix_power", n, iters, setup, "linalg.matrix_power(A, 3)",464
lambda G=G: (G,), lambda o, i: float(np.linalg.matrix_power(o[0], 3)[0, 0]),465
extract="[0, 0]")467
# ---- trace / Frobenius norm (cheap O(n²) reductions) ----468
for n, iters in [(32, 100000), (256, 20000)]:469
A = spd(n)470
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"471
bench("trace", n, iters, setup, "linalg.trace(A)",472
lambda A=A: (A,), lambda o, i: float(np.trace(o[0])), extract="")473
for n, iters in [(32, 100000), (256, 10000)]:474
A = spd(n)475
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])"476
bench("norm(matrix)", n, iters, setup, "linalg.norm(A)",477
lambda A=A: (A,), lambda o, i: float(np.linalg.norm(o[0])), extract="")479
# ---- outer product / Kronecker product ----480
for n, iters in [(64, 50000), (256, 4000)]:481
u, v = rng.standard_normal(n), rng.standard_normal(n)482
setup = f"let u = ndarray.array({lit(u)})\nlet v = ndarray.array({lit(v)})"483
bench("outer", n, iters, setup, "linalg.outer(u, v)",484
lambda u=u, v=v: (u, v), lambda o, i: float(np.outer(o[0], o[1])[0, 0]),485
extract="[0, 0]")486
for n, iters in [(8, 40000), (16, 8000), (32, 1500)]:487
A, B = spd(n), spd(n)488
setup = f"let A = ndarray.reshape(ndarray.array({lit(A)}), [{n}, {n}])\n" \489
f"let B = ndarray.reshape(ndarray.array({lit(B)}), [{n}, {n}])"490
bench("kron", n, iters, setup, "linalg.kron(A, B)",491
lambda A=A, B=B: (A, B), lambda o, i: float(np.kron(o[0], o[1])[0, 0]),492
extract="[0, 0]")496
def run_ndarray():497
"""The elementwise ufuncs the ndarray README publishes. Kept separate from the linalg498
cases because they are a different claim on different operands — one table holding both499
would sit a 64-element sqrt beside a 96x96 matmul as though they were comparable."""500
for n, iters, fn in [(64, 200000, "sqrt"), (16384, 3000, "sqrt"),501
(16384, 3000, "exp"), (16384, 3000, "sin")]:502
x = np.abs(rng.standard_normal(n)) + 0.5503
setup = f"let X = ndarray.array({lit(x)})"504
bench(f"ndarray.{fn}", n, iters, setup, f"ndarray.{fn}(X)",505
lambda x=x: (x,), lambda o, i, fn=fn: float(getattr(np, fn)(o[0])[0]))507
# Array + scalar is the OPERATOR, not add(): ndarray.add takes two arrays (ndarray.hpp:993),508
# while the scalar broadcast is operator+ (ndarray.hpp:1365). The ndarray README billed this509
# row as "16384-element array + scalar" against `ndarray.add`, which is not a form that510
# exists — measure what the code actually provides and name it accordingly.511
x = rng.standard_normal(16384)512
setup = f"let X = ndarray.array({lit(x)})"513
bench("X + scalar", 16384, 3000, setup, "X + 1.5",514
lambda x=x: (x,), lambda o, i: float((o[0] + 1.5)[0]))516
# And the two-array form, which is what ndarray.add really is.517
y = rng.standard_normal(16384)518
setup2 = f"let X = ndarray.array({lit(x)})\nlet Y = ndarray.array({lit(y)})"519
bench("ndarray.add", 16384, 3000, setup2, "ndarray.add(X, Y)",520
lambda x=x, y=y: (x, y), lambda o, i: float((o[0] + o[1])[0]))523
if __name__ == "__main__":524
main()