6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 1#!/usr/bin/env python3
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 2"""Offline benchmark of the Pac-Man problem in Python (NumPy and SciPy).
4A port of the two best methods from the in-browser benchmark, the harmonic
5basis and the graded-panel double layer, on the same problem, the same error
6lattice, and the same timing protocol as the chapter driver. Run
8 python3 bench/python/bench_pacman.py
10from anywhere; it writes bench-pacman-python.json into the current directory.
11Load that file from the benchmark chapter's "load offline results" button.
12Because it ran on the same machine as your browser, the timings share a clock
13and can share an axis.
15The point of the comparison is the language and its libraries, not the
16author's fluency in each, so each port is written the way its own language
17wants to be written. For NumPy that means whole-array expressions, with two
18details that matter more than they look:
20 * the kernel matrix is assembled in place. Written as the obvious sequence
21 of broadcast expressions it allocates five n-by-n temporaries and spends
22 more time on memory traffic than the dense solve itself spends on
23 arithmetic; reusing two buffers cuts that stage by nearly half. Writing
24 it as a scalar double loop, which is the fastest shape in MATLAB and
25 Julia, would be catastrophic here, and that asymmetry is the honest
26 reason the ports do not look alike.
28 * the least-squares fit goes through a QR driver. numpy.linalg.lstsq is an
29 SVD (LAPACK gelsd), while MATLAB's backslash and Julia's backslash are
30 both pivoted QR, so using lstsq would have Python solving a harder
31 problem than the others and calling it the same benchmark.
34import json
35import platform
36import sys
37import time
38from datetime import date
40import numpy as np
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 41from scipy.linalg import lstsq as qr_lstsq
43KIND = "both"
44Z0 = 1.6 + 1.3j
47def exact(x, y, kind=KIND):
48 """lp_lexact: the exact solution, singular and/or smooth part."""
49 th = np.arctan2(y, x)
50 th = np.where(th < -1e-13, th + 2 * np.pi, th)
51 r = np.hypot(x, y)
52 sing = r ** (2.0 / 3.0) * np.sin(2.0 * th / 3.0)
53 smooth = 0.5 * np.log((x - Z0.real) ** 2 + (y - Z0.imag) ** 2)
54 if kind == "sing":
55 return sing
56 if kind == "smooth":
57 return smooth
58 return sing + smooth
61def lattice():
62 """The driver's error lattice: spacing 1/8, margins from the boundary."""
63 s = np.arange(-8, 9) / 8.0
64 X, Y = np.meshgrid(s, s)
65 R2 = X**2 + Y**2
66 keep = (
67 (R2 <= 0.7225 + 1e-9)
68 & (R2 > 1e-12)
69 & ~((X > 1e-12) & (Y < -1e-12))
70 & ~((X > 1e-12) & (Y < 0.25 - 1e-9))
71 & ~((Y < -1e-12) & (X > -0.25 + 1e-9))
72 )
73 return X[keep], Y[keep]
76# ------------------------------------------- the method of particular sols --
79def pac_basis(x, y, M, K):
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 80 """Re z^m and Im z^m up to degree M, then K corner functions.
82 Written into a preallocated Fortran-ordered array: LAPACK wants columns
83 contiguous, and building a list of columns for column_stack would copy the
84 whole matrix once more than it has to.
85 """
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 86 z = x + 1j * y
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 87 A = np.empty((x.size, 2 * M + 1 + K), order="F")
88 A[:, 0] = 1.0
89 col = 1
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 90 p = z.copy()
91 for _ in range(M):
93 A[:, col + 1] = p.imag
94 col += 2
95 p *= z
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 96 th = np.arctan2(y, x)
97 th = np.where(th < -1e-13, th + 2 * np.pi, th)
98 r = np.abs(z)
99 k, q = 0, 0
100 while q < K:
101 k += 1
102 if k % 3 != 0:
103 q += 1
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 104 A[:, col] = r ** (2.0 * k / 3.0) * np.sin(2.0 * k * th / 3.0)
105 col += 1
106 return A
109def pac_bpoints(ntot):
110 ne = max(6, round(ntot / 6.7))
111 na = max(12, ntot - 2 * ne)
112 t = (np.arange(1, ne + 1) - 0.5) / ne
113 c = t**2
114 s = 1.5 * np.pi * (np.arange(1, na + 1) - 0.5) / na
115 bx = np.concatenate([c, np.cos(s), np.zeros(ne)])
116 by = np.concatenate([np.zeros(ne), np.sin(s), -c])
117 return bx, by
120def mps_solve(n, px, py):
121 nb = 3 * n + 1
122 bx, by = pac_bpoints(3 * nb)
123 A = pac_basis(bx, by, n, n)
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 124 # gelsy is pivoted QR: what MATLAB's and Julia's backslash do on a
125 # rectangular system. The default lstsq driver is an SVD and does more
126 # work than either.
127 c = qr_lstsq(A, exact(bx, by), lapack_driver="gelsy")[0]
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 128 return pac_basis(px, py, n, n) @ c
131# --------------------------------------------- boundary integral, panels ----
134def graded_breaks(T, d, base):
135 half = T / 2.0
136 left = np.concatenate([[0.0], half * 2.0 ** np.arange(-d, 1)])
137 bp0 = np.concatenate([left, T - left[-2::-1]])
138 bp = [bp0[0]]
139 for a, b in zip(bp0[:-1], bp0[1:]):
140 ns = max(1, int(np.ceil((b - a) / base)))
141 for q in range(1, ns + 1):
142 bp.append(a + (b - a) * q / ns)
143 return np.array(bp)
146def panel_nodes(bp, gx, gw):
147 a = bp[:-1]
148 b = bp[1:]
149 t = ((a + b) / 2)[:, None] + ((b - a) / 2)[:, None] * gx[None, :]
150 w = ((b - a) / 2)[:, None] * gw[None, :]
151 return t.ravel(), w.ravel()
154def pac_panels(d):
155 gx, gw = np.polynomial.legendre.leggauss(12)
156 be = graded_breaks(1.0, d, 0.35)
157 ba = graded_breaks(1.5 * np.pi, d, 0.35)
158 t1, w1 = panel_nodes(be, gx, gw)
159 t2, w2 = panel_nodes(ba, gx, gw)
160 t3, w3 = panel_nodes(be, gx, gw)
161 z1 = np.zeros_like(t1)
162 z3 = np.zeros_like(t3)
163 qx = np.concatenate([t1, np.cos(t2), z3])
164 qy = np.concatenate([z1, np.sin(t2), t3 - 1.0])
165 qnx = np.concatenate([z1, np.cos(t2), np.ones_like(t3)])
166 qny = np.concatenate([-np.ones_like(t1), np.sin(t2), z3])
167 qw = np.concatenate([w1, w2, w3])
168 qkap = np.concatenate([z1, np.ones_like(t2), z3])
169 return qx, qy, qnx, qny, qw, qkap
172def bie_solve(d, px, py):
173 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
174 n = qx.size
176 # A is built in two buffers and finished in place: every operation below
177 # writes into memory that already exists. The obvious version, one
178 # broadcast expression per line, needs five n-by-n arrays and is about
179 # 1.6 times slower for that reason alone.
180 A = qx[:, None] - qx[None, :] # dx
181 B = qy[:, None] - qy[None, :] # dy
182 R2 = A * A
183 R2 += B * B
184 np.fill_diagonal(R2, 1.0) # keep 0/0 off the diagonal
185 A *= qnx[None, :]
186 B *= qny[None, :]
187 A += B
188 A /= R2
189 A *= qw[None, :] / (2 * np.pi)
190 # the diagonal is the smooth limit of the kernel, plus the jump relation
191 diag = np.arange(n)
192 A[diag, diag] = -0.5 - qkap * qw / (4 * np.pi)
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 194 sigma = np.linalg.solve(A, exact(qx, qy))
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 196 EX = px[:, None] - qx[None, :]
197 EY = py[:, None] - qy[None, :]
198 KE = (qnx[None, :] * EX + qny[None, :] * EY) / (2 * np.pi * (EX**2 + EY**2))
199 return (KE * qw[None, :]) @ sigma
202# ------------------------------------------------------------- the sweep ----
205def sweep(solve, ns, px, py, ve):
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 206 """Warm up, then time each resolution the way the chapter driver does.
208 Keep the fastest of several runs. Three runs at least whenever a solve
209 is cheap enough to allow it, many more when it is very cheap, and one
210 when it is expensive: the budget below caps the wall clock per point.
211 The minimum over runs is what makes these numbers comparable at all,
212 here and in the offline scripts, which apply the same rule. A single
213 measurement of a fifty-millisecond solve reports whatever else the
214 machine happened to be doing at the time.
215 """
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 216 E, T = [], []
217 solve(ns[0], px[:1], py[:1]) # warm-up
218 for n in ns:
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 219 el, tot, rep = float("inf"), 0.0, 0
220 while rep == 0 or (rep < 50 and tot < 0.6 and (rep < 3 or tot < 0.1)):
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 221 t0 = time.perf_counter()
222 v = solve(n, px, py)
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 223 dt = time.perf_counter() - t0
224 el = min(el, dt)
225 tot += dt
226 rep += 1
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 227 E.append(float(np.max(np.abs(v - ve))))
228 T.append(el * 1000.0)
229 return E, T
233 """How many threads the BLAS behind NumPy will use, if we can find out."""
234 try:
235 from threadpoolctl import threadpool_info
237 n = max((p["num_threads"] for p in threadpool_info()), default=None)
238 if n:
239 return n
240 except ImportError:
241 pass
242 return None
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 245def main():
246 px, py = lattice()
247 ve = exact(px, py)
248 methods = []
249 for name, slot, solve, ns in [
250 ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
251 ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
252 ]:
253 E, T = sweep(solve, ns, px, py, ve)
254 methods.append({"name": name, "slot": slot, "ns": ns, "E": E, "T": T})
255 print(name)
256 for n, e, t in zip(ns, E, T):
257 print(f" n={n:3d} err={e:.3e} {t:8.2f} ms")
260 label = f"Python {sys.version_info.major}.{sys.version_info.minor} + NumPy {np.__version__}"
261 if nt:
262 label += f", {nt} threads"
264 "schema": "laplace-dirichlet-bench-1",
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 266 "kind": KIND,
267 "date": date.today().isoformat(),
268 "machine": platform.machine(),
269 "methods": methods,
270 }
271 with open("bench-pacman-python.json", "w") as f:
272 json.dump(out, f)
273 print("wrote bench-pacman-python.json")
276if __name__ == "__main__":
277 main()