geodesics.py
11.8 kB · python · 329 lines
1import heapq2import math3import sys4import time56import numpy as np78R2 = math.sqrt(2.0)9R5 = math.sqrt(5.0)10LEVEL1 = 2 * R5 - 3 * R21112# DESIGN1314def code(n):15 return sum(1 << (n * i + j) for i in range(n) for j in range(n) if not (i % 2 and j % 2))1617def voids(n, level):18 out = []19 def walk(x0, y0, size, depth):20 if depth == level:21 return22 c = size // n23 for i in range(n):24 for j in range(n):25 if i % 2 and j % 2:26 out.append((x0 + i * c, y0 + j * c, c))27 else:28 walk(x0 + i * c, y0 + j * c, c, depth + 1)29 walk(0, 0, n ** level, 0)30 return np.array(out, dtype=float).reshape(-1, 3)3132# VISIBILITY3334def blocked(p, Q, H):35 if len(H) == 0:36 return np.zeros(len(Q), bool)37 out = np.zeros(len(Q), bool)38 dx = Q[:, 0] - p[0]39 dy = Q[:, 1] - p[1]40 bx0 = np.minimum(p[0], Q[:, 0])41 bx1 = np.maximum(p[0], Q[:, 0])42 by0 = np.minimum(p[1], Q[:, 1])43 by1 = np.maximum(p[1], Q[:, 1])44 step = max(1, 4_000_000 // max(1, len(Q)))45 for s in range(0, len(H), step):46 a0 = H[s:s + step, 0][None, :]47 b0 = H[s:s + step, 1][None, :]48 a1 = a0 + H[s:s + step, 2][None, :]49 b1 = b0 + H[s:s + step, 2][None, :]50 near = (a1 > bx0[:, None]) & (a0 < bx1[:, None]) & (b1 > by0[:, None]) & (b0 < by1[:, None])51 if not near.any():52 continue53 with np.errstate(divide="ignore", invalid="ignore"):54 t1 = (a0 - p[0]) / dx[:, None]55 t2 = (a1 - p[0]) / dx[:, None]56 u1 = (b0 - p[1]) / dy[:, None]57 u2 = (b1 - p[1]) / dy[:, None]58 zx = dx[:, None] == 059 zy = dy[:, None] == 060 inx = (a0 < p[0]) & (p[0] < a1)61 iny = (b0 < p[1]) & (p[1] < b1)62 lox = np.where(zx, np.where(inx, -np.inf, np.inf), np.minimum(t1, t2))63 hix = np.where(zx, np.where(inx, np.inf, -np.inf), np.maximum(t1, t2))64 loy = np.where(zy, np.where(iny, -np.inf, np.inf), np.minimum(u1, u2))65 hiy = np.where(zy, np.where(iny, np.inf, -np.inf), np.maximum(u1, u2))66 lo = np.maximum(np.maximum(lox, loy), 0.0)67 hi = np.minimum(np.minimum(hix, hiy), 1.0)68 out |= (near & (hi - lo > 1e-12)).any(1)69 return out7071def shortest(H, A, B, bound):72 A = np.array(A, float)73 B = np.array(B, float)74 c = H[:, :2] + H[:, 2:3] / 275 r = H[:, 2] * R2 / 276 keep = np.hypot(*(c - A).T) + np.hypot(*(c - B).T) - 2 * r <= bound77 H = H[keep]78 C = np.concatenate([H[:, :2], H[:, :2] + H[:, 2:3] * [1, 0], H[:, :2] + H[:, 2:3] * [0, 1], H[:, :2] + H[:, 2:3]])79 C = np.unique(C, axis=0)80 C = C[np.hypot(*(C - A).T) + np.hypot(*(C - B).T) <= bound + 1e-9]81 P = np.concatenate([[A, B], C])82 n = len(P)83 tail = np.hypot(*(P - B).T)84 dist = np.full(n, np.inf)85 dist[0] = 0.086 done = np.zeros(n, bool)87 heap = [(tail[0], 0)]88 while heap:89 f, u = heapq.heappop(heap)90 if done[u]:91 continue92 done[u] = True93 if u == 1:94 return dist[1], len(H), n95 cand = np.where(~done)[0]96 step = np.hypot(*(P[cand] - P[u]).T)97 cand = cand[dist[u] + step + tail[cand] <= bound + 1e-9]98 if len(cand) == 0:99 continue100 ok = ~blocked(P[u], P[cand], H)101 v = cand[ok]102 nd = dist[u] + np.hypot(*(P[v] - P[u]).T)103 for w, d in zip(v[nd < dist[v] - 1e-13], nd[nd < dist[v] - 1e-13]):104 dist[w] = d105 heapq.heappush(heap, (d + tail[w], w))106 return np.inf, len(H), n107108def corner(n, level):109 S = n ** level110 H = voids(n, level)111 slack = 0.05112 while True:113 bound = (R2 + slack) * S114 d, nh, nv = shortest(H, (0, 0), (S, S), bound)115 if d < np.inf:116 return d / S, len(H), nv117 slack *= 2118119def envelope(n, level):120 return 2 - (2 - R2) * (1 - 1 / n) ** level121122# VERBS123124def verb_code():125 print("code of the side-N tile, bit N*i + j, void iff both digits odd")126 for n in (3, 5, 7, 9, 11):127 print(f" N = {n:<2} code {code(n)} fill {n * n - ((n - 1) // 2) ** 2}")128 assert code(3) == 495129130def verb_corner():131 print("corner distance D(N, L) = d((0,0), (1,1)), exact visibility graph")132 cases = [(3, 1), (3, 2), (3, 3), (3, 4), (5, 1), (5, 2), (7, 1), (7, 2), (9, 2)]133 if len(sys.argv) > 2:134 cases = [tuple(map(int, a.split(","))) for a in sys.argv[2:]]135 last = {}136 for n, level in cases:137 t = time.time()138 d, nh, nv = corner(n, level)139 env = envelope(n, level)140 assert R2 <= d <= env + 1e-12141 assert d >= last.get(n, 0.0) - 1e-12142 last[n] = d143 if n == 3:144 assert abs(d - 2 * R5 / 3) < 1e-12145 print(f" N = {n:<2} L = {level} D = {d:.12f} N(D - sqrt2) = {n * (d - R2):.6f} 2 - D = {2 - d:.6f} envelope {env:.6f} voids {nh} vertices {nv} {time.time() - t:.1f} s")146147def verb_level1():148 print("level 1: D(N, 1) against sqrt2 + (2 sqrt5 - 3 sqrt2)/N")149 worst = 0.0150 for n in range(3, 42, 2):151 d, _, _ = corner(n, 1)152 f = R2 + LEVEL1 / n153 worst = max(worst, abs(d - f))154 print(f" N = {n:<2} D = {d:.12f} formula {f:.12f} diff {d - f:+.1e}")155 print(f" largest difference {worst:.1e}; 2 sqrt5 - 3 sqrt2 = {LEVEL1:.9f}")156157# MAP158159def periodic_edges(K=3):160 corners = np.array([(1.0, 1.0), (2.0, 1.0), (1.0, 2.0), (2.0, 2.0)])161 H = np.array([(1.0 + 2 * i, 1.0 + 2 * j, 1.0) for i in range(-K - 2, K + 3) for j in range(-K - 2, K + 3)])162 ea, eb, ev = [], [], []163 for a, pa in enumerate(corners):164 Q, idx = [], []165 for b, pb in enumerate(corners):166 for i in range(-K, K + 1):167 for j in range(-K, K + 1):168 q = pb + 2 * np.array([i, j])169 if np.any(q != pa):170 Q.append(q)171 idx.append(b)172 Q = np.array(Q)173 ok = ~blocked(pa, Q, H)174 for q, b in zip(Q[ok], np.array(idx)[ok]):175 ea.append(a)176 eb.append(b)177 ev.append(q - pa)178 return np.array(ea), np.array(eb), np.array(ev)179180def gauge(h, P, V):181 return np.max((V @ P.T) / h[None, :], axis=1)182183def homogenise(h, P, E):184 ea, eb, ev = E185 w = gauge(h, P, ev)186 g0 = P @ ev.T187 lo = np.zeros(len(P))188 hi = np.full(len(P), 2.0)189 flat = ea * 4 + eb190 for _ in range(48):191 lam = (lo + hi) / 2192 g = g0 - lam[:, None] * w[None, :]193 G = np.full((len(P), 16), -np.inf)194 for k in range(16):195 sel = flat == k196 if sel.any():197 G[:, k] = g[:, sel].max(1)198 G = G.reshape(len(P), 4, 4)199 for k in range(4):200 G = np.maximum(G, G[:, :, k:k + 1] + G[:, k:k + 1, :])201 pos = (np.diagonal(G, axis1=1, axis2=2) > 1e-12).any(1)202 lo = np.where(pos, lam, lo)203 hi = np.where(pos, hi, lam)204 return (lo + hi) / 2205206def octagon(V):207 a = np.abs(V)208 big = np.maximum(a[:, 0], a[:, 1])209 small = np.minimum(a[:, 0], a[:, 1])210 return big - small + R2 * small211212def node_cycles():213 from itertools import combinations, permutations214 out = [[(a, a)] for a in range(4)]215 for k in (2, 3, 4):216 for sub in combinations(range(4), k):217 for rest in permutations(sub[1:]):218 seq = (sub[0],) + rest219 if k == 2 and rest[0] < sub[0]:220 continue221 out.append([(seq[i], seq[(i + 1) % k]) for i in range(k)])222 return out223224def hull_gauge(W, V):225 from scipy.spatial import ConvexHull226 eq = ConvexHull(W).equations227 return np.max((V @ eq[:, :2].T) / (-eq[:, 2])[None, :], axis=1)228229def cycle_points(cost, P, E, cycles):230 ea, eb, ev = E231 g0 = P @ ev.T232 pair = ea * 4 + eb233 groups = [np.where(pair == k)[0] for k in range(16)]234 lo = np.zeros(len(P))235 hi = np.full(len(P), 2.0)236 def best(lam):237 g = g0 - lam[:, None] * cost[None, :]238 G = np.full((len(P), 16), -np.inf)239 A = np.zeros((len(P), 16), int)240 for k, idx in enumerate(groups):241 if len(idx):242 j = g[:, idx].argmax(1)243 G[:, k] = g[np.arange(len(P)), idx[j]]244 A[:, k] = idx[j]245 S = np.stack([sum(G[:, a * 4 + b] for a, b in c) for c in cycles], 1)246 return S, A247 for _ in range(48):248 lam = (lo + hi) / 2249 S, _ = best(lam)250 pos = S.max(1) > 1e-12251 lo = np.where(pos, lam, lo)252 hi = np.where(pos, hi, lam)253 S, A = best(lo)254 pick = S.argmax(1)255 W = np.zeros((len(P), 2))256 for t in range(len(P)):257 idx = [A[t, a * 4 + b] for a, b in cycles[pick[t]]]258 W[t] = ev[idx].sum(0) / cost[idx].sum()259 return W260261def verb_map():262 ndir = 720263 t = time.time()264 th = np.linspace(0, 2 * np.pi, ndir, endpoint=False)265 P = np.stack([np.cos(th), np.sin(th)], 1)266 U = np.stack([np.cos(th + np.pi / ndir), np.sin(th + np.pi / ndir)], 1)267 probe = np.array([[1, 0], [math.cos(math.pi / 8), math.sin(math.pi / 8)], [2 / R5, 1 / R5], [1 / R2, 1 / R2]])268 names = ["0", "22.5", "atan(1/2)", "45"]269 cycles = node_cycles()270 print(f"homogenisation map: upper bounds on nu_L from cycle points, {ndir} support directions, {len(cycles)} node cycles")271 for K in (3, 6):272 E = periodic_edges(K)273 cost = np.hypot(*E[2].T)274 print(f" window {K} periods, {len(E[0])} edges")275 print(" L " + " ".join(f"{s:>10}" for s in names) + " gap to octagon at least")276 for level in range(1, 21):277 W = cycle_points(cost, P, E, cycles)278 cost = hull_gauge(W, E[2])279 nu = hull_gauge(W, probe)280 gap = np.max(1 - hull_gauge(W, np.concatenate([P, U])) / octagon(np.concatenate([P, U])))281 if level <= 4 or level % 5 == 0:282 print(f" {level:<4} " + " ".join(f"{x:10.6f}" for x in nu) + f" {gap:.2e}")283 print(" oct " + " ".join(f"{x:10.6f}" for x in octagon(probe)))284 two = np.array([[4.0, 2.0]])285 print(f" explicit path (0,0)-(2,1)-(3,2)-(4,2): {R5 + R2 + 1:.6f} against octagon {octagon(two)[0]:.6f} for displacement (4,2)")286 print(f" {time.time() - t:.1f} s")287288def verb_bridge():289 E = periodic_edges()290 ndir = 1440291 th = np.linspace(0, 2 * np.pi, ndir, endpoint=False)292 P = np.stack([np.cos(th), np.sin(th)], 1)293 h = homogenise(np.ones(ndir), P, E)294 print("bridge at level 1: d_N((0,0), (1, (N-1)/(2N))) against the level-1 norm of the same vector")295 for n in (11, 21, 31, 41):296 H = voids(n, 1)297 B = (n, (n - 1) // 2)298 nu = gauge(h, P, np.array([[1.0, (n - 1) / (2 * n)]]))[0]299 d, _, _ = shortest(H, (0, 0), B, 1.2 * n)300 print(f" N = {n:<2} d = {d / n:.6f} level-1 norm {nu:.6f} euclid {math.hypot(1, (n - 1) / (2 * n)):.6f} N(d - norm) = {n * (d / n - nu):.4f}")301302def verb_hull():303 from itertools import combinations, product304 from scipy.spatial import ConvexHull305 V = []306 for i in range(3):307 for s in (1, -1):308 e = np.zeros(3)309 e[i] = s310 V.append(e)311 for i, j in combinations(range(3), 2):312 for s, t in product((1, -1), repeat=2):313 e = np.zeros(3)314 e[i], e[j] = s, t315 V.append(e / R2)316 hull = ConvexHull(np.array(V))317 faces = len(np.unique(np.round(hull.equations, 9), axis=0))318 edges = len({tuple(sorted(p)) for s in hull.simplices for p in combinations(s, 2)})319 print(f"dim 3 ball: hull of the 18 free unit vectors has {len(hull.vertices)} vertices, {edges} edges, {faces} faces, all triangles: {faces == len(hull.simplices)}")320 assert (len(hull.vertices), edges, faces) == (18, 48, 32)321322VERBS = {"code": verb_code, "corner": verb_corner, "level1": verb_level1, "map": verb_map, "bridge": verb_bridge, "hull": verb_hull}323324if __name__ == "__main__":325 names = sys.argv[1:2] or list(VERBS)326 for name in names:327 t = time.time()328 VERBS[name]()329 print(f"[{name}] {time.time() - t:.1f} s", flush=True)