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.
15using LinearAlgebra
16using Dates
17using Printf
19const KIND = "both"
20const Z0 = 1.6 + 1.3im
22function exact(x, y)
23 th = atan(y, x)
24 th < -1e-13 && (th += 2pi)
25 r = hypot(x, y)
26 sing = r^(2 / 3) * sin(2 * th / 3)
27 smooth = 0.5 * log((x - real(Z0))^2 + (y - imag(Z0))^2)
28 KIND == "sing" && return sing
29 KIND == "smooth" && return smooth
30 return sing + smooth
31end
32exact(x::AbstractVector, y::AbstractVector) = [exact(x[k], y[k]) for k in eachindex(x)]
34"""The driver's error lattice: spacing 1/8, margins from the boundary."""
35function lattice()
36 px = Float64[]
37 py = Float64[]
38 for x in -1.0:0.125:1.0, y in -1.0:0.125:1.0
39 r2 = x^2 + y^2
40 keep = r2 <= 0.7225 + 1e-9 && r2 > 1e-12 &&
41 !(x > 1e-12 && y < -1e-12) &&
42 !(x > 1e-12 && y < 0.25 - 1e-9) &&
43 !(y < -1e-12 && x > -0.25 + 1e-9)
44 keep || continue
45 push!(px, x)
46 push!(py, y)
47 end
48 return px, py
49end
51# ------------------------------------------- the method of particular sols --
53function pac_basis(x, y, M, K)
54 np = length(x)
55 A = zeros(np, 2M + 1 + K)
56 z = complex.(x, y)
57 col = 1
58 p = ones(ComplexF64, np)
59 for mdeg in 0:M
60 A[:, col] = real.(p)
61 col += 1
62 if mdeg > 0
63 A[:, col] = imag.(p)
64 col += 1
65 end
66 p .*= z
67 end
68 th = [t < -1e-13 ? t + 2pi : t for t in atan.(y, x)]
69 r = abs.(z)
70 k, q = 0, 0
71 while q < K
72 k += 1
73 k % 3 == 0 && continue
74 q += 1
75 A[:, col] = r .^ (2k / 3) .* sin.(2k .* th ./ 3)
76 col += 1
77 end
78 return A
79end
81function pac_bpoints(ntot)
82 ne = max(6, round(Int, ntot / 6.7))
83 na = max(12, ntot - 2ne)
84 t = ((1:ne) .- 0.5) ./ ne
85 c = t .^ 2
86 s = 1.5pi .* ((1:na) .- 0.5) ./ na
87 bx = vcat(c, cos.(s), zeros(ne))
88 by = vcat(zeros(ne), sin.(s), -c)
89 return bx, by
90end
92function mps_solve(n, px, py)
93 nb = 3n + 1
94 bx, by = pac_bpoints(3nb)
95 A = pac_basis(bx, by, n, n)
96 c = A \ exact(bx, by) # pivoted QR least squares
97 return pac_basis(px, py, n, n) * c
98end
100# --------------------------------------------- boundary integral, panels ----
102"""Gauss-Legendre nodes and weights on [-1, 1], by Newton on P_p."""
103function gauss(p)
104 x = [cos(pi * (k - 0.25) / (p + 0.5)) for k in 1:p]
105 dp = ones(p)
106 for _ in 1:50
107 p0 = ones(p)
108 p1 = copy(x)
109 for k in 2:p
110 p2 = ((2k - 1) .* x .* p1 .- (k - 1) .* p0) ./ k
111 p0 = p1
112 p1 = p2
113 end
114 dp = p .* (x .* p1 .- p0) ./ (x .^ 2 .- 1)
115 dx = p1 ./ dp
116 x .-= dx
117 maximum(abs, dx) < 1e-15 && break
118 end
119 w = 2 ./ ((1 .- x .^ 2) .* dp .^ 2)
120 return x, w
121end
123function graded_breaks(T, d, base)
124 half = T / 2
125 left = vcat(0.0, half .* 2.0 .^ (-d:0))
126 bp0 = vcat(left, T .- reverse(left[1:end-1]))
127 bp = [bp0[1]]
128 for k in 2:length(bp0)
129 a, b = bp0[k-1], bp0[k]
130 ns = max(1, ceil(Int, (b - a) / base))
131 for q in 1:ns
132 push!(bp, a + (b - a) * q / ns)
133 end
134 end
135 return bp
136end
138function panel_nodes(bp, gx, gw)
139 t = Float64[]
140 w = Float64[]
141 for k in 1:length(bp)-1
142 a, b = bp[k], bp[k+1]
143 append!(t, (a + b) / 2 .+ (b - a) / 2 .* gx)
144 append!(w, (b - a) / 2 .* gw)
145 end
146 return t, w
147end
149function pac_panels(d)
150 gx, gw = gauss(12)
151 be = graded_breaks(1.0, d, 0.35)
152 ba = graded_breaks(1.5pi, d, 0.35)
153 t1, w1 = panel_nodes(be, gx, gw)
154 t2, w2 = panel_nodes(ba, gx, gw)
155 t3, w3 = panel_nodes(be, gx, gw)
156 n1, n3 = length(t1), length(t3)
157 qx = vcat(t1, cos.(t2), zeros(n3))
158 qy = vcat(zeros(n1), sin.(t2), t3 .- 1)
159 qnx = vcat(zeros(n1), cos.(t2), ones(n3))
160 qny = vcat(-ones(n1), sin.(t2), zeros(n3))
161 qw = vcat(w1, w2, w3)
162 qkap = vcat(zeros(n1), ones(length(t2)), zeros(n3))
163 return qx, qy, qnx, qny, qw, qkap
164end
166function bie_solve(d, px, py)
167 qx, qy, qnx, qny, qw, qkap = pac_panels(d)
168 n = length(qx)
169 A = zeros(n, n)
170 for j in 1:n, i in 1:n
171 if i == j
172 A[i, j] = -0.5 - qkap[j] / (4pi) * qw[j]
173 else
174 dx = qx[i] - qx[j]
175 dy = qy[i] - qy[j]
176 A[i, j] = (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j]
177 end
178 end
179 sigma = A \ exact(qx, qy)
180 v = zeros(length(px))
181 for k in eachindex(px)
182 s = 0.0
183 for j in 1:n
184 dx = px[k] - qx[j]
185 dy = py[k] - qy[j]
186 s += (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j] * sigma[j]
187 end
188 v[k] = s
189 end
190 return v
191end
193# ------------------------------------------------------------- the sweep ----
195function sweep(solve, ns, px, py, ve)
196 E = Float64[]
197 T = Float64[]
198 solve(ns[1], px[1:1], py[1:1]) # warm-up, which in Julia is compilation
199 for n in ns
200 el = v = solve(n, px, py)
201 # a second timed run always: stray compilation and GC pauses
202 # otherwise land in single measurements
203 el = min(el, v = solve(n, px, py))
204 rep = 0
205 while el < 0.005 && rep < 4
206 rep += 1
207 el = min(el, v = solve(n, px, py))
208 end
209 push!(E, maximum(abs.(v .- ve)))
210 push!(T, el * 1000)
211 end
212 return E, T
213end
215function json_num(v)
216 isfinite(v) || return "null"
217 return string(v)
218end
219json_arr(v) = "[" * join(json_num.(float.(v)), ",") * "]"
221function main()
222 px, py = lattice()
223 ve = exact(px, py)
224 entries = String[]
225 for (name, slot, solve, ns) in [
226 ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
227 ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
228 ]
229 E, T = sweep(solve, ns, px, py, ve)
230 println(name)
231 for k in eachindex(ns)
232 (" n=%3d err=%.3e %8.2f ms\n", ns[k], E[k], T[k])
233 end
234 push!(entries,
235 """{"name":"$(name)","slot":"$(slot)","ns":$(json_arr(ns)),"E":$(json_arr(E)),"T":$(json_arr(T))}""")
236 end
237 label = "Julia $(VERSION.major).$(VERSION.minor)"
238 open("bench-pacman-julia.json", "w") do f
239 write(f,
240 """{"schema":"laplace-dirichlet-bench-1","label":"$(label)","kind":"$(KIND)","date":"$(Dates.today())","methods":[$(join(entries, ","))]}""")
241 end
242 println("wrote bench-pacman-julia.json")
243end
245main()