moire_local_limit.py
16.0 kB · python · 392 lines
1import time2from fractions import Fraction as Q3from math import ceil, floor, gcd, log10, sqrt45import numpy as np67WINDOWS = 328CODE = 79LADDER = [(21, 23), (101, 103), (2321, 2323), (23231, 23233), (232321, 232323)]10GAPS = [(2321, 2323), (2320, 2322), (2321, 2322), (2321, 2325), (2319, 2325)]11LEVEL2 = [201, 1001, 2001, 200, 1000, 2000]12FINE = 102413FINE_PAIR = (232321, 232323)1415# DESIGNS1617def up(x, digits=4):18 scale = 10 ** (digits - 1 - floor(log10(x)))19 return ceil(x * scale) / scale2021def corners(code):22 return [[(code >> (2 * a + b)) & 1 for b in (0, 1)] for a in (0, 1)]2324def walsh(code):25 t = corners(code)26 tau = [[1 - 2 * t[a][b] for b in (0, 1)] for a in (0, 1)]27 out = {}28 for sa in (0, 1):29 for sb in (0, 1):30 out[(sa, sb)] = Q(sum(tau[a][b] * (-1) ** (sa * a + sb * b) for a in (0, 1) for b in (0, 1)), 4)31 return out3233def fill(code):34 return Q(bin(code).count("1"), 4)3536# TRIANGLE WAVE3738def dist2(x):39 r = x - 2 * (x // 2)40 return r if r <= 1 else 2 - r4142def tri(x):43 return 1 - 2 * dist2(x)4445def dist2_integral(t):46 whole = t // 247 r = t - 2 * whole48 part = r * r / 2 if r <= 1 else 2 * r - r * r / 2 - 149 return whole + part5051def tri_mean(h, a, b):52 return 1 - 2 * (dist2_integral(h * b) - dist2_integral(h * a)) / (h * (b - a))5354# LEVEL ONE5556def law1(n, m, w):57 big = n * m * w58 cuts = np.unique(np.concatenate([np.arange(0, big + 1, m * w), np.arange(0, big + 1, n * w), np.arange(0, big + 1, n * m)]))59 left, size = cuts[:-1], np.diff(cuts)60 state = (left // (m * w)) % 2 + 2 * ((left // (n * w)) % 2)61 window = left // (n * m)62 counts = np.bincount(window * 4 + state, weights=size.astype(np.float64), minlength=4 * w)63 return counts.reshape(w, 4) / float(n * m), np.rint(counts).astype(np.int64).reshape(w, 4)6465def kernel1(code):66 t = corners(code)67 k = np.zeros((4, 4))68 for sx in range(4):69 for sy in range(4):70 k[sx, sy] = float(t[sx & 1][sy & 1] != t[sx >> 1][sy >> 1])71 return k7273def lag_law(rho):74 same, diff = (1 + rho) / 4, (1 - rho) / 475 return [same, diff, diff, same]7677def limit_law1(h, w):78 rows = []79 for i in range(w):80 rows.append([float(x) for x in lag_law(tri_mean(h, Q(i, w), Q(i + 1, w)))])81 return np.array(rows)8283def walsh_limit(code, rs, rt):84 c = walsh(code)85 return (1 - c[(0, 0)] ** 2 - c[(1, 0)] ** 2 * rs - c[(0, 1)] ** 2 * rt - c[(1, 1)] ** 2 * rs * rt) / 28687def kernel_limit(code, rs, rt):88 t = corners(code)89 ps, pt = lag_law(rs), lag_law(rt)90 return sum(ps[sx] * pt[sy] for sx in range(4) for sy in range(4) if t[sx & 1][sy & 1] != t[sx >> 1][sy >> 1])9192def code7_h(u, v):93 return (1 - abs(1 - 2 * u) * abs(1 - 2 * v)) / 29495def check_walsh():96 points = [Q(k, 12) for k in range(13)]97 bad = 098 for code in range(16):99 for u in points:100 for v in points:101 for h in (1, 2, 3, 4, 6):102 if walsh_limit(code, tri(h * u), tri(h * v)) != kernel_limit(code, tri(h * u), tri(h * v)):103 bad += 1104 classes = {}105 for code in range(16):106 c = walsh(code)107 key = (c[(1, 0)] ** 2, c[(0, 1)] ** 2, c[(1, 1)] ** 2)108 classes.setdefault(key, []).append(code)109 exact7 = all(walsh_limit(CODE, tri(2 * u), tri(2 * v)) == code7_h(u, v) for u in points for v in points)110 means = all(walsh_limit(code, tri_mean(h, Q(0), Q(1)), tri_mean(h, Q(0), Q(1))) == 2 * fill(code) * (1 - fill(code)) for code in range(16) for h in range(1, 9))111 print(f"walsh form against the lag kernel: {bad} mismatches over 16 codes, gaps 1 2 3 4 6, 169 points")112 print(f"code 7 at gap 2 is (1 - |1-2u| |1-2v|)/2 on all 169 points: {exact7}")113 print(f"mean over the square is 2 p (1 - p) for all 16 codes and gaps 1..8: {means}")114 for key, codes in sorted(classes.items()):115 print(f" squared walsh (first, second, both) = {tuple(str(x) for x in key)}: codes {codes}")116117def direct2d(code, n, m, w):118 side = n * m119 i = np.arange(side)120 a, b = (i // m) % 2, (i // n) % 2121 t = np.array(corners(code))122 ink_n = t[a[:, None], a[None, :]]123 ink_m = t[b[:, None], b[None, :]]124 x = (ink_n != ink_m).astype(np.int64)125 step = side // w126 return x.reshape(w, step, w, step).sum(axis=(1, 3)), step * step127128def check_direct():129 n, m, w = 21, 23, 7130 bad = 0131 for code in range(16):132 counts, cells = direct2d(code, n, m, w)133 _, exact = law1(n, m, w)134 k = kernel1(code).astype(np.int64)135 fact = exact @ k @ exact.T136 bad += int(np.sum(fact != counts * (n * m) ** 2 // cells))137 print(f"2d raster of side 483 against the factorised box mean, 16 codes, 7 by 7 windows: {bad} mismatches")138139def ladder():140 print(f"code 7, gap 2, {WINDOWS} by {WINDOWS} windows of side 1/{WINDOWS}, error = max |box mean - H(centre)|")141 k = kernel1(CODE)142 us = (np.arange(WINDOWS) + 0.5) / WINDOWS143 hc = (1 - np.abs(1 - 2 * us)[:, None] * np.abs(1 - 2 * us)[None, :]) / 2144 lim = limit_law1(2, WINDOWS)145 exact_centre = np.max(np.abs(lim @ k @ lim.T - hc))146 print(f" box average of H equals H at the box centre: {exact_centre:.1e}")147 for n, m in LADDER:148 t0 = time.time()149 law, _ = law1(n, m, WINDOWS)150 means = law @ k @ law.T151 err = np.max(np.abs(means - hc))152 bound = (20 * (m - n) + 36 * WINDOWS) / n153 print(f" {n}/{m}: mean {means.mean():.6f}, error {err:.3e}, error N l {err * n / WINDOWS:.4f}, proved bound {up(bound)}, {time.time() - t0:.2f}s")154155def gaps():156 print("other gaps and shared factors, code 7, error against the gap-h limit, global covariance of the two bits")157 k = kernel1(CODE)158 for n, m in GAPS:159 h = m - n160 law, exact = law1(n, m, WINDOWS)161 lim = limit_law1(h, WINDOWS)162 err = np.max(np.abs(law @ k @ law.T - lim @ k @ lim.T))163 tot = exact.sum(axis=0)164 whole = n * m * WINDOWS165 p11, pa, pb = Q(int(tot[3]), whole), Q(int(tot[1] + tot[3]), whole), Q(int(tot[2] + tot[3]), whole)166 cov = p11 - pa * pb167 print(f" {n}/{m}: gap {h}, gcd {gcd(n, m)}, error {err:.3e}, global covariance {cov}")168169def obstruction():170 print("local correlation against global covariance, code 7, 2321/2323")171 n, m = 2321, 2323172 law, exact = law1(n, m, WINDOWS)173 a = law[:, 1] + law[:, 3]174 b = law[:, 2] + law[:, 3]175 pearson = (law[:, 3] - a * b) / np.sqrt(a * (1 - a) * b * (1 - b))176 lim = np.array([float(tri_mean(2, Q(i, WINDOWS), Q(i + 1, WINDOWS))) for i in range(WINDOWS)])177 print(f" local correlation runs {pearson.min():.4f} to {pearson.max():.4f}, max distance to the window mean of the triangle wave {np.max(np.abs(pearson - lim)):.2e}")178 tot = exact.sum(axis=0)179 whole = n * m * WINDOWS180 p11, pa, pb = Q(int(tot[3]), whole), Q(int(tot[1] + tot[3]), whole), Q(int(tot[2] + tot[3]), whole)181 qn, qm = Q(n - 1, 2 * n), Q(m - 1, 2 * m)182 print(f" global covariance {p11 - pa * pb}, marginals {pa == qn} {pb == qm}")183 fill_n, fill_m = 1 - qn * qn, 1 - qm * qm184 glob = fill_n + fill_m - 2 * fill_n * fill_m185 k = kernel1(CODE)186 joint = [Q(int(x), whole) for x in tot]187 direct = sum(joint[sx] * joint[sy] for sx in range(4) for sy in range(4) if k[sx, sy])188 print(f" global overlay mean from the window integrals {direct} = {float(direct):.9f}, independent value {glob}, equal {direct == glob}, limit 3/8")189190# OCTAGON191192def octagon():193 kappa = sqrt(2) - 1194 c = (1 - kappa) / 2195 sag = (1 - sqrt(kappa)) ** 2 / sqrt(2)196 chord = sqrt(2) * (1 - kappa)197 print(f"regular chord octagon at kappa = sqrt 2 - 1, level c = 1 - 1/sqrt 2 = {c:.6f}")198 print(f" chord {chord:.6f}, sag {sag:.6f}, sag over chord {sag / chord:.6f} in centred units")199 n, m = FINE_PAIR200 law, _ = law1(n, m, FINE)201 k = kernel1(CODE)202 chord_u = (1 - (1 + kappa) / 2) / 2203 arc_u = (1 - sqrt(kappa)) / 2204 out = []205 for u in (chord_u, arc_u):206 i = int(u * FINE)207 val = law[i] @ k @ law[i]208 lim = lag_law(tri_mean(2, Q(i, FINE), Q(i + 1, FINE)))209 lv = sum(float(lim[sx]) * float(lim[sy]) * k[sx, sy] for sx in range(4) for sy in range(4))210 out.append((u, val, lv))211 print(f" {n}/{m}, window side 1/{FINE} on the diagonal:")212 print(f" chord midpoint u = {out[0][0]:.6f}: box mean {out[0][1]:.6f}, box average of H {out[0][2]:.6f}, H there 0.25")213 print(f" arc point u = {out[1][0]:.6f}: box mean {out[1][1]:.6f}, box average of H {out[1][2]:.6f}, level {c:.6f}")214215# LEVEL TWO216217def law2(n, m, w):218 big = n * n * m * m * w219 cuts = np.unique(np.concatenate([np.arange(0, big + 1, m * m * w), np.arange(0, big + 1, n * n * w), np.arange(0, big + 1, n * n * m * m)]))220 left, size = cuts[:-1], np.diff(cuts)221 fn, fm = left // (m * m * w), left // (n * n * w)222 a, c = (fn // n) % 2, (fn % n) % 2223 b, d = (fm // m) % 2, (fm % m) % 2224 state = a + 2 * b + 4 * c + 8 * d225 window = left // (n * n * m * m)226 counts = np.bincount(window * 16 + state, weights=size.astype(np.float64), minlength=16 * w)227 return counts.reshape(w, 16) / float(n * n * m * m)228229def kernel2(code):230 t = corners(code)231 k = np.zeros((16, 16))232 for sx in range(16):233 for sy in range(16):234 ax, bx, cx, dx = sx & 1, (sx >> 1) & 1, (sx >> 2) & 1, (sx >> 3) & 1235 ay, by, cy, dy = sy & 1, (sy >> 1) & 1, (sy >> 2) & 1, (sy >> 3) & 1236 ink_n = t[ax][ay] & t[cx][cy]237 ink_m = t[bx][by] & t[dx][dy]238 k[sx, sy] = float(ink_n != ink_m)239 return k240241def par(x):242 return int((x // 1) % 2)243244def limit_law2(u, odd):245 s = 2 * u246 cuts = {Q(0), Q(1), Q(2)}247 for k in range(-4, 5):248 cuts.add(k - s)249 for floor_phi in (0, 1):250 for j in (0, 1, 2):251 for k in range(-12, 24):252 cuts.add(floor_phi + Q(k, 4) - s / 2 - Q(odd * j, 4))253 pts = sorted(x for x in cuts if 0 <= x <= 2)254 law = [Q(0)] * 16255 for lo, hi in zip(pts, pts[1:]):256 if hi == lo:257 continue258 phi = (lo + hi) / 2259 f = phi - (phi // 1)260 j = (f + s) // 1261 lam = 4 * f + 2 * s + odd * j262 rho = tri(lam)263 a, b = par(phi), par(phi + s)264 wgt = (hi - lo) / 2265 for c in (0, 1):266 for d in (0, 1):267 law[a + 2 * b + 4 * c + 8 * d] += wgt * (1 + (1 if c == d else -1) * rho) / 4268 return law269270def simpson_law2(lo, hi, odd):271 pa, pm, pb = limit_law2(lo, odd), limit_law2((lo + hi) / 2, odd), limit_law2(hi, odd)272 return [(x + 4 * y + z) / 6 for x, y, z in zip(pa, pm, pb)]273274def quadratic_check(odd):275 bad = 0276 for q in range(4):277 lo = Q(q, 4)278 xs = [lo + Q(k, 20) for k in range(6)]279 vals = [limit_law2(x, odd) for x in xs]280 for s in range(16):281 y = [v[s] for v in vals]282 d2 = [y[i] - 2 * y[i + 1] + y[i + 2] for i in range(4)]283 if len(set(d2)) != 1:284 bad += 1285 return bad286287def h2(u, v, odd, k):288 pu, pv = limit_law2(u, odd), limit_law2(v, odd)289 return sum(pu[sx] * pv[sy] for sx in range(16) for sy in range(16) if k[sx, sy])290291def level2():292 k = kernel2(CODE)293 print("level 2, code 7, gap 2")294 for odd in (1, 0):295 print(f" limit law piecewise quadratic on each quarter, {'odd' if odd else 'even'} sides: {quadratic_check(odd)} breaches over 16 states")296 lims = {}297 for odd in (1, 0):298 rows = [simpson_law2(Q(i, WINDOWS), Q(i + 1, WINDOWS), odd) for i in range(WINDOWS)]299 lims[odd] = np.array([[float(x) for x in r] for r in rows])300 for n in LEVEL2:301 t0 = time.time()302 odd = n % 2303 law = law2(n, n + 2, WINDOWS)304 means = law @ k @ law.T305 lim = lims[odd] @ k @ lims[odd].T306 other = lims[1 - odd] @ k @ lims[1 - odd].T307 us = (np.arange(WINDOWS) + 0.5) / WINDOWS308 g = np.abs(1 - 2 * us)309 gg = g[:, None] * g[None, :]310 na = 2 * 9 / 16 - 2 * 9 / 16 * (0.5 + gg / 4)311 nb = 2 * 9 / 16 - 2 * (0.5 + gg / 4) ** 2312 print(f" {n}/{n + 2}: error {np.max(np.abs(means - lim)):.3e}, against the other parity {np.max(np.abs(means - other)):.3e}, naive {np.max(np.abs(means - na)):.3e} and {np.max(np.abs(means - nb)):.3e}, {time.time() - t0:.2f}s")313 for odd in (1, 0):314 whole = [sum(x) for x in zip(*[simpson_law2(Q(q, 4), Q(q + 1, 4), odd) for q in range(4)])]315 whole = [x / 4 for x in whole]316 top = all(sum(whole[s] for s in range(16) if (s & 3) == ab) == Q(1, 4) for ab in range(4))317 indep = all(whole[s] == Q(1, 16) for s in range(16))318 mean = sum(whole[sx] * whole[sy] for sx in range(16) for sy in range(16) if k[sx, sy])319 print(f" {'odd' if odd else 'even'} sides: global law uniform on 16 states {indep}, top pair uniform {top}, mean of the limit {mean} = {float(mean):.6f}, 2 p (1 - p) at p = 9/16 is {2 * Q(9, 16) * Q(7, 16)}")320 closed_form(k)321322def tee(y):323 r = y - (y // 1)324 return (1 if (y // 1) % 2 == 0 else -1) * r * (1 - r)325326def weight(u):327 return -tee(4 * min(u, 1 - u)) / 2328329def alpha(u):330 return (1 + tri(2 * u)) / 4331332def gamma(u, odd):333 return (1 + 2 * weight(u)) / 4 if odd else Q(1, 4)334335def beta(u):336 return alpha(u) / 4 + weight(u) / 8337338def h2_closed(u, v, odd):339 return Q(5, 8) - alpha(u) * alpha(v) - gamma(u, odd) * gamma(v, odd) - 2 * beta(u) * beta(v)340341def h2_naive(u, v):342 return Q(9, 16) - Q(9, 8) * alpha(u) * alpha(v)343344def closed_form(k):345 pts = [Q(i, 96) for i in range(97)]346 for odd in (1, 0):347 bad = 0348 for u in pts:349 law = limit_law2(u, odd)350 a = sum(law[s] for s in range(16) if s & 3 == 3)351 c = sum(law[s] for s in range(16) if s >> 2 == 3)352 b = law[15]353 bad += int(a != alpha(u)) + int(c != gamma(u, odd)) + int(b != beta(u))354 grid = [Q(i, 12) for i in range(13)]355 bad2 = sum(int(h2(u, v, odd, k) != h2_closed(u, v, odd)) for u in grid for v in grid)356 print(f" {'odd' if odd else 'even'} sides: alpha, beta, gamma against the lag law at 97 points {bad} mismatches, H2 closed form against the kernel at 169 points {bad2} mismatches")357 pts = [Q(0), Q(1, 8), Q(1, 4)]358 for odd in (1, 0):359 mat = [[h2_closed(u, v, odd) for v in pts] for u in pts]360 det = (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) + mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]))361 p1, p2 = h2_closed(Q(0), Q(3, 8), odd), h2_closed(Q(1, 4), Q(1, 4), odd)362 print(f" {'odd' if odd else 'even'}: H2 corner {h2_closed(Q(0), Q(0), odd)}, centre {h2_closed(Q(1, 2), Q(1, 2), odd)}, at (3/8, 3/8) {h2_closed(Q(3, 8), Q(3, 8), odd)}")363 print(f" 3 by 3 minor at u, v in 0, 1/8, 1/4: {det}")364 print(f" g(u) g(v) = 1/4 at both (0, 3/8) and (1/4, 1/4): H2 reads {p1} and {p2}")365 fine = [Q(i, 256) for i in range(257)]366 parts = {odd: [(alpha(u), beta(u), gamma(u, odd)) for u in fine] for odd in (1, 0)}367 def at(p, q):368 return Q(5, 8) - p[0] * q[0] - p[2] * q[2] - 2 * p[1] * q[1]369 gap = max(abs(at(p, q) - at(r, t)) for p, r in zip(parts[1], parts[0]) for q, t in zip(parts[1], parts[0]))370 na = {odd: max(abs(at(p, q) - Q(9, 16) + Q(9, 8) * p[0] * q[0]) for p in parts[odd] for q in parts[odd]) for odd in (1, 0)}371 print(f" max |H2 odd - H2 even| on the 1/256 grid {gap} = {float(gap):.6f}, attained at (3/8, 3/8): {abs(h2_closed(Q(3, 8), Q(3, 8), 1) - h2_closed(Q(3, 8), Q(3, 8), 0)) == gap}")372 print(f" max |H2 - naive| on the 1/256 grid, a lower bound on the sup, truncated: odd {floor(float(na[1]) * 1e6) / 1e6:.6f}, even {floor(float(na[0]) * 1e6) / 1e6:.6f}")373 grid = [Q(i, 64) + Q(1, 128) for i in range(64)]374 for odd in (1, 0):375 mat = np.array([[float(h2_closed(u, v, odd)) for v in grid] for u in grid])376 sv = np.linalg.svd(mat, compute_uv=False)377 rank = int(np.sum(sv > 1e-12 * sv[0]))378 print(f" kernel rank on a 64 by 64 grid, {'odd' if odd else 'even'}: {rank}, singular values {', '.join(f'{x:.3e}' for x in sv[:rank + 1])}")379380def main():381 t0 = time.time()382 check_walsh()383 check_direct()384 ladder()385 gaps()386 obstruction()387 octagon()388 level2()389 print(f"total {time.time() - t0:.1f}s")390391if __name__ == "__main__":392 main()