6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 1#!/usr/bin/env python3
2"""Offline benchmark of the Pac-Man problem in Python (NumPy).
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.
14"""
16import json
17import platform
18import sys
19import time
20from datetime import date
22import numpy as np
24KIND = "both"
25Z0 = 1.6 + 1.3j
28def exact(x, y, kind=KIND):
29 """lp_lexact: the exact solution, singular and/or smooth part."""
30 th = np.arctan2(y, x)
31 th = np.where(th < -1e-13, th + 2 * np.pi, th)
32 r = np.hypot(x, y)
33 sing = r ** (2.0 / 3.0) * np.sin(2.0 * th / 3.0)
34 smooth = 0.5 * np.log((x - Z0.real) ** 2 + (y - Z0.imag) ** 2)
35 if kind == "sing":
36 return sing
37 if kind == "smooth":
38 return smooth
39 return sing + smooth
42def lattice():
43 """The driver's error lattice: spacing 1/8, margins from the boundary."""
44 s = np.arange(-8, 9) / 8.0
45 X, Y = np.meshgrid(s, s)
46 R2 = X**2 + Y**2
47 keep = (
48 (R2 <= 0.7225 + 1e-9)
49 & (R2 > 1e-12)
50 & ~((X > 1e-12) & (Y < -1e-12))
51 & ~((X > 1e-12) & (Y < 0.25 - 1e-9))
52 & ~((Y < -1e-12) & (X > -0.25 + 1e-9))
53 )
54 return X[keep], Y[keep]
57# ------------------------------------------- the method of particular sols --
60def pac_basis(x, y, M, K):
61 z = x + 1j * y
62 cols = [np.ones_like(x)]
63 p = z.copy()
64 for _ in range(M):
65 cols.append(p.real)
66 cols.append(p.imag)
67 p = p * z
68 th = np.arctan2(y, x)
69 th = np.where(th < -1e-13, th + 2 * np.pi, th)
70 r = np.abs(z)
71 k, q = 0, 0
72 while q < K:
73 k += 1
74 if k % 3 != 0:
75 q += 1
76 cols.append(r ** (2.0 * k / 3.0) * np.sin(2.0 * k * th / 3.0))
77 return np.column_stack(cols)
80def pac_bpoints(ntot):
81 ne = max(6, round(ntot / 6.7))
82 na = max(12, ntot - 2 * ne)
83 t = (np.arange(1, ne + 1) - 0.5) / ne
84 c = t**2
85 s = 1.5 * np.pi * (np.arange(1, na + 1) - 0.5) / na
86 bx = np.concatenate([c, np.cos(s), np.zeros(ne)])
87 by = np.concatenate([np.zeros(ne), np.sin(s), -c])
88 return bx, by
91def mps_solve(n, px, py):
92 nb = 3 * n + 1
93 bx, by = pac_bpoints(3 * nb)
94 A = pac_basis(bx, by, n, n)
95 # QR-based least squares, the closest match to MATLAB's backslash
96 c, *_ = np.linalg.lstsq(A, exact(bx, by), rcond=None)
97 return pac_basis(px, py, n, n) @ c
100# --------------------------------------------- boundary integral, panels ----
103def graded_breaks(T, d, base):
104 half = T / 2.0
105 left = np.concatenate([[0.0], half * 2.0 ** np.arange(-d, 1)])
106 bp0 = np.concatenate([left, T - left[-2::-1]])
107 bp = [bp0[0]]
108 for a, b in zip(bp0[:-1], bp0[1:]):
109 ns = max(1, int(np.ceil((b - a) / base)))
110 for q in range(1, ns + 1):
111 bp.append(a + (b - a) * q / ns)
112 return np.array(bp)
115def panel_nodes(bp, gx, gw):
116 a = bp[:-1]
117 b = bp[1:]
118 t = ((a + b) / 2)[:, None] + ((b - a) / 2)[:, None] * gx[None, :]
119 w = ((b - a) / 2)[:, None] * gw[None, :]
120 return t.ravel(), w.ravel()
123def pac_panels(d):
124 gx, gw = np.polynomial.legendre.leggauss(12)
125 be = graded_breaks(1.0, d, 0.35)
126 ba = graded_breaks(1.5 * np.pi, d, 0.35)
127 t1, w1 = panel_nodes(be, gx, gw)
128 t2, w2 = panel_nodes(ba, gx, gw)
129 t3, w3 = panel_nodes(be, gx, gw)
130 z1 = np.zeros_like(t1)
131 z3 = np.zeros_like(t3)
132 qx = np.concatenate([t1, np.cos(t2), z3])
133 qy = np.concatenate([z1, np.sin(t2), t3 - 1.0])
134 qnx = np.concatenate([z1, np.cos(t2), np.ones_like(t3)])
135 qny = np.concatenate([-np.ones_like(t1), np.sin(t2), z3])
136 qw = np.concatenate([w1, w2, w3])
137 qkap = np.concatenate([z1, np.ones_like(t2), z3])
138 return qx, qy, qnx, qny, qw, qkap
141def bie_solve(d, px, py):
142 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
143 n = qx.size
144 DX = qx[:, None] - qx[None, :]
145 DY = qy[:, None] - qy[None, :]
146 R2 = DX**2 + DY**2
147 np.fill_diagonal(R2, 1.0)
148 K = (qnx[None, :] * DX + qny[None, :] * DY) / (2 * np.pi * R2)
149 np.fill_diagonal(K, -qkap / (4 * np.pi))
150 A = -0.5 * np.eye(n) + K * qw[None, :]
151 sigma = np.linalg.solve(A, exact(qx, qy))
152 EX = px[:, None] - qx[None, :]
153 EY = py[:, None] - qy[None, :]
154 KE = (qnx[None, :] * EX + qny[None, :] * EY) / (2 * np.pi * (EX**2 + EY**2))
155 return (KE * qw[None, :]) @ sigma
158# ------------------------------------------------------------- the sweep ----
161def sweep(solve, ns, px, py, ve):
162 E, T = [], []
163 solve(ns[0], px[:1], py[:1]) # warm-up
164 for n in ns:
165 t0 = time.perf_counter()
166 v = solve(n, px, py)
167 el = time.perf_counter() - t0
168 rep = 0
169 while el < 0.005 and rep < 4:
170 rep += 1
171 t0 = time.perf_counter()
172 v = solve(n, px, py)
173 el = min(el, time.perf_counter() - t0)
174 E.append(float(np.max(np.abs(v - ve))))
175 T.append(el * 1000.0)
176 return E, T
179def main():
180 px, py = lattice()
181 ve = exact(px, py)
182 methods = []
183 for name, slot, solve, ns in [
184 ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
185 ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
186 ]:
187 E, T = sweep(solve, ns, px, py, ve)
188 methods.append({"name": name, "slot": slot, "ns": ns, "E": E, "T": T})
189 print(name)
190 for n, e, t in zip(ns, E, T):
191 print(f" n={n:3d} err={e:.3e} {t:8.2f} ms")
193 out = {
194 "schema": "laplace-dirichlet-bench-1",
195 "label": f"Python {sys.version_info.major}.{sys.version_info.minor} + NumPy {np.__version__}",
196 "kind": KIND,
197 "date": date.today().isoformat(),
198 "machine": platform.machine(),
199 "methods": methods,
200 }
201 with open("bench-pacman-python.json", "w") as f:
202 json.dump(out, f)
203 print("wrote bench-pacman-python.json")
206if __name__ == "__main__":
207 main()