1#!/usr/bin/env julia
2# Offline benchmark of the Pac-Man problem in Julia.
3#
4# A port of the two best methods from the in-browser benchmark, the harmonic
5# basis and the graded-panel double layer, on the same problem, the same error
6# lattice, and the same timing protocol as the chapter driver. Run
7#
8# julia bench/julia/bench_pacman.jl
9#
10# from anywhere; it writes bench-pacman-julia.json into the current directory.
11# Load that file from the benchmark chapter's "load offline results" button.
12# Because it ran on the same machine as your browser, the timings share a
13# clock and can share an axis. Uses only the standard library.
14#
15# Each port is written the way its own language wants to be written, since the
16# comparison is meant to be between languages and not between the author's
17# habits in each. Here that means scalar loops with @inbounds for the matrix
18# assembly, which is both the fastest and the most readable shape available;
19# the same loop in NumPy would be hopeless, which is why the Python port
20# assembles with whole-array expressions instead. Where the work is one
21# LAPACK call, all three ports do the same thing and should agree.
23using LinearAlgebra
24using Dates
25using Printf
27const KIND = "both"
28const Z0 = 1.6 + 1.3im
30function exact(x, y)
31 th = atan(y, x)
32 th < -1e-13 && (th += 2pi)
33 r = hypot(x, y)
34 sing = r^(2 / 3) * sin(2 * th / 3)
35 smooth = 0.5 * log((x - real(Z0))^2 + (y - imag(Z0))^2)
36 KIND == "sing" && return sing
37 KIND == "smooth" && return smooth
38 return sing + smooth
39end
40exact(x::AbstractVector, y::AbstractVector) = [exact(x[k], y[k]) for k in eachindex(x)]
42"""The driver's error lattice: spacing 1/8, margins from the boundary."""
43function lattice()
44 px = Float64[]
45 py = Float64[]
46 for x in -1.0:0.125:1.0, y in -1.0:0.125:1.0
47 r2 = x^2 + y^2
48 keep = r2 <= 0.7225 + 1e-9 && r2 > 1e-12 &&
49 !(x > 1e-12 && y < -1e-12) &&
50 !(x > 1e-12 && y < 0.25 - 1e-9) &&
51 !(y < -1e-12 && x > -0.25 + 1e-9)
52 keep || continue
53 push!(px, x)
54 push!(py, y)
55 end
56 return px, py
57end
59# ------------------------------------------- the method of particular sols --
61function pac_basis(x, y, M, K)
62 np = length(x)
63 A = zeros(np, 2M + 1 + K)
64 z = complex.(x, y)
65 col = 1
66 p = ones(ComplexF64, np)
67 for mdeg in 0:M
68 A[:, col] = real.(p)
69 col += 1
70 if mdeg > 0
71 A[:, col] = imag.(p)
72 col += 1
73 end
74 p .*= z
75 end
76 th = [t < -1e-13 ? t + 2pi : t for t in atan.(y, x)]
77 r = abs.(z)
78 k, q = 0, 0
79 while q < K
80 k += 1
81 k % 3 == 0 && continue
82 q += 1
83 A[:, col] = r .^ (2k / 3) .* sin.(2k .* th ./ 3)
84 col += 1
85 end
86 return A
87end
89function pac_bpoints(ntot)
90 ne = max(6, round(Int, ntot / 6.7))
91 na = max(12, ntot - 2ne)
92 t = ((1:ne) .- 0.5) ./ ne
93 c = t .^ 2
94 s = 1.5pi .* ((1:na) .- 0.5) ./ na
95 bx = vcat(c, cos.(s), zeros(ne))
96 by = vcat(zeros(ne), sin.(s), -c)
97 return bx, by
98end
100function mps_solve(n, px, py)
101 nb = 3n + 1
102 bx, by = pac_bpoints(3nb)
103 A = pac_basis(bx, by, n, n)
104 c = A \ exact(bx, by) # pivoted QR least squares
105 return pac_basis(px, py, n, n) * c
106end
108# --------------------------------------------- boundary integral, panels ----
110"""Gauss-Legendre nodes and weights on [-1, 1], by Newton on P_p."""
111function gauss(p)
112 x = [cos(pi * (k - 0.25) / (p + 0.5)) for k in 1:p]
113 dp = ones(p)
114 for _ in 1:50
115 p0 = ones(p)
116 p1 = copy(x)
117 for k in 2:p
118 p2 = ((2k - 1) .* x .* p1 .- (k - 1) .* p0) ./ k
119 p0 = p1
120 p1 = p2
121 end
122 dp = p .* (x .* p1 .- p0) ./ (x .^ 2 .- 1)
123 dx = p1 ./ dp
124 x .-= dx
125 maximum(abs, dx) < 1e-15 && break
126 end
127 w = 2 ./ ((1 .- x .^ 2) .* dp .^ 2)
128 return x, w
129end
131function graded_breaks(T, d, base)
132 half = T / 2
133 left = vcat(0.0, half .* 2.0 .^ (-d:0))
134 bp0 = vcat(left, T .- reverse(left[1:end-1]))
135 bp = [bp0[1]]
136 for k in 2:length(bp0)
137 a, b = bp0[k-1], bp0[k]
138 ns = max(1, ceil(Int, (b - a) / base))
139 for q in 1:ns
140 push!(bp, a + (b - a) * q / ns)
141 end
142 end
143 return bp
144end
146function panel_nodes(bp, gx, gw)
147 t = Float64[]
148 w = Float64[]
149 for k in 1:length(bp)-1
150 a, b = bp[k], bp[k+1]
151 append!(t, (a + b) / 2 .+ (b - a) / 2 .* gx)
152 append!(w, (b - a) / 2 .* gw)
153 end
154 return t, w
155end
157function pac_panels(d)
158 gx, gw = gauss(12)
159 be = graded_breaks(1.0, d, 0.35)
160 ba = graded_breaks(1.5pi, d, 0.35)
161 t1, w1 = panel_nodes(be, gx, gw)
162 t2, w2 = panel_nodes(ba, gx, gw)
163 t3, w3 = panel_nodes(be, gx, gw)
164 n1, n3 = length(t1), length(t3)
165 qx = vcat(t1, cos.(t2), zeros(n3))
166 qy = vcat(zeros(n1), sin.(t2), t3 .- 1)
167 qnx = vcat(zeros(n1), cos.(t2), ones(n3))
168 qny = vcat(-ones(n1), sin.(t2), zeros(n3))
169 qw = vcat(w1, w2, w3)
170 qkap = vcat(zeros(n1), ones(length(t2)), zeros(n3))
171 return qx, qy, qnx, qny, qw, qkap
172end
174function bie_solve(d, px, py)
175 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
176 n = length(qx)
177 # column by column, so the writes run down memory, and with the bounds
178 # checks elided: every index here is a loop variable over 1:n
179 A = Matrix{Float64}(undef, n, n)
180 c = 1 / (2pi)
181 for j in 1:n
182 xj = qx[j]; yj = qy[j]
183 nxj = qnx[j]; nyj = qny[j]; wj = qw[j]
184 for i in 1:n
185 dx = qx[i] - xj
186 dy = qy[i] - yj
187 A[i, j] = c * (nxj * dx + nyj * dy) / (dx * dx + dy * dy) * wj
188 end
189 # the smooth limit of the kernel on the diagonal, and the jump relation
190 A[j, j] = -0.5 - qkap[j] * wj / (4pi)
191 end
193 sigma = A \ exact(qx, qy)
195 v = zeros(length(px))
196 for k in eachindex(px)
197 xk = px[k]; yk = py[k]
198 s = 0.0
199 for j in 1:n
200 dx = xk - qx[j]
201 dy = yk - qy[j]
202 s += (qnx[j] * dx + qny[j] * dy) / (dx * dx + dy * dy) * qw[j] * sigma[j]
203 end
204 v[k] = c * s
205 end
206 return v
207end
209# ------------------------------------------------------------- the sweep ----
211function sweep(solve, ns, px, py, ve)
212 # Warm up, then time each resolution the way the chapter driver does.
213 # Keep the fastest of several runs. Three runs at least whenever a solve
214 # is cheap enough to allow it, many more when it is very cheap, and one
215 # when it is expensive: the budget below caps the wall clock per point.
216 # The minimum over runs is what makes these numbers comparable at all,
217 # here and in the offline scripts, which apply the same rule. A single
218 # measurement of a fifty-millisecond solve reports whatever else the
219 # machine happened to be doing at the time.
220 E = Float64[]
221 T = Float64[]
222 local v
223 solve(ns[1], px[1:1], py[1:1]) # warm-up, which in Julia is compilation
224 for n in ns
225 el = Inf
226 tot = 0.0
227 rep = 0
228 while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
229 dt = v = solve(n, px, py)
230 el = min(el, dt)
231 tot += dt
232 rep += 1
233 end
234 push!(E, maximum(abs.(v .- ve)))
235 push!(T, el * 1000)
236 end
237 return E, T
238end
240function json_num(v)
241 isfinite(v) || return "null"
242 return string(v)
243end
244json_arr(v) = "[" * join(json_num.(float.(v)), ",") * "]"
246function main()
247 px, py = lattice()
248 ve = exact(px, py)
249 entries = String[]
250 for (name, slot, solve, ns) in [
251 ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
252 ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
253 ]
254 E, T = sweep(solve, ns, px, py, ve)
255 println(name)
256 for k in eachindex(ns)
257 (" n=%3d err=%.3e %8.2f ms\n", ns[k], E[k], T[k])
258 end
259 push!(entries,
260 """{"name":"$(name)","slot":"$(slot)","ns":$(json_arr(ns)),"E":$(json_arr(E)),"T":$(json_arr(T))}""")
261 end
262 label = "Julia $(VERSION.major).$(VERSION.minor), $(BLAS.get_num_threads()) threads"
263 open("bench-pacman-julia.json", "w") do f
264 write(f,
265 """{"schema":"laplace-dirichlet-bench-1","label":"$(label)","kind":"$(KIND)","date":"$(Dates.today())","methods":[$(join(entries, ","))]}""")
266 end
267 println("wrote bench-pacman-julia.json")
268end
270main()