detector.py
13.0 kB · python · 384 lines
1import math2import sys3import time45NGRID = 60006WINDOW_FLOOR = 10007OMEGA_LO, OMEGA_HI, OMEGA_N = 0.5, 14.0, 13518TOL = 0.019RATIO = 10.010COMB = 811OM3 = 2 * math.pi / math.log(3)12OM5 = 2 * math.pi / math.log(5)1314INDEPENDENT = (3, (0, 1), 5, (0, 1, 2))15DEPENDENT = (3, (0, 1), 9, (0, 1, 3))1617# ENUMERATION1819def walk(L, p, pd, q, qd, sink):20 pw_p = [p ** i for i in range(L + 1)]21 pw_q = [1]22 while pw_q[-1] <= pw_p[L]:23 pw_q.append(pw_q[-1] * q)24 allowed = frozenset(qd)25 stack = [(0, 0, len(pw_q) - 1)]26 nodes = 027 while stack:28 j, a, s = stack.pop()29 nodes += 130 hi = a + pw_p[L - j] - 131 t = s32 while t > 0 and a // pw_q[t - 1] == hi // pw_q[t - 1]:33 t -= 134 ok = True35 for u in range(t, s):36 if (a // pw_q[u]) % q not in allowed:37 ok = False38 break39 if not ok:40 continue41 if j == L:42 sink(a)43 continue44 step = pw_p[L - j - 1]45 for d in reversed(pd):46 stack.append((j + 1, a + d * step, t))47 return nodes484950def in_base(n, b, ds):51 if n == 0:52 return 0 in ds53 while n:54 if n % b not in ds:55 return False56 n //= b57 return True585960def brute(L, p, pd, q, qd):61 return [n for n in range(p ** L) if in_base(n, p, pd) and in_base(n, q, qd)]626364def digit_count(N, q, ds):65 if N <= 0:66 return 067 k = len(ds)68 dig = []69 m = N70 while m:71 dig.append(m % q)72 m //= q73 total = 074 for i in range(len(dig) - 1, -1, -1):75 total += sum(1 for d in ds if d < dig[i]) * k ** i76 if dig[i] not in ds:77 break78 return total - (1 if 0 in ds else 0)798081def histogram(L, p, pd, q, qd):82 h = L * math.log(p) / (NGRID - 1)83 cnt = [0] * (NGRID + 2)84 log = math.log8586 def sink(n, cnt=cnt, h=h, log=log):87 if n:88 i = int(log(n) / h) + 189 if i <= NGRID:90 cnt[i] += 19192 nodes = walk(L, p, pd, q, qd, sink)93 run = 094 out = []95 for i in range(NGRID):96 run += cnt[i]97 out.append(run)98 return h, out, nodes99100101def one_base_ladder(L, p, b, ds):102 h = L * math.log(p) / (NGRID - 1)103 return h, [digit_count(int(math.exp(i * h)), b, ds) for i in range(NGRID)]104105106def crop(h, counts):107 start = 0108 while start < len(counts) and counts[start] < WINDOW_FLOOR:109 start += 1110 return [i * h for i in range(start, len(counts))], list(counts[start:])111112# DETECTOR113114def blackman(m):115 return [0.42 - 0.5 * math.cos(2 * math.pi * j / (m - 1)) + 0.08 * math.cos(4 * math.pi * j / (m - 1)) for j in range(m)]116117118def power_at(y, h, w):119 z = complex(math.cos(w * h), -math.sin(w * h))120 acc = 0j121 for v in reversed(y):122 acc = acc * z + v123 return abs(acc) ** 2124125126def golden(f, lo, hi):127 phi = (math.sqrt(5) - 1) / 2128 a, b = lo, hi129 c, d = b - phi * (b - a), a + phi * (b - a)130 fc, fd = f(c), f(d)131 for _ in range(60):132 if fc < fd:133 b, d, fd = d, c, fc134 c = b - phi * (b - a)135 fc = f(c)136 else:137 a, c, fc = c, d, fd138 d = a + phi * (b - a)139 fd = f(d)140 return (a + b) / 2141142143def detrend(us, cs, deg):144 ys = [math.log(c) for c in cs]145 mid = (us[0] + us[-1]) / 2146 half = (us[-1] - us[0]) / 2147 ts = [(x - mid) / half for x in us]148 basis = []149 for k in range(deg + 1):150 v = [t ** k for t in ts]151 for b in basis:152 d = sum(a * c for a, c in zip(v, b))153 v = [a - d * c for a, c in zip(v, b)]154 nrm = math.sqrt(sum(a * a for a in v))155 basis.append([a / nrm for a in v])156 g = list(ys)157 for b in basis:158 d = sum(a * c for a, c in zip(g, b))159 g = [a - d * c for a, c in zip(g, b)]160 n = len(ts)161 mt = sum(ts) / n162 my = sum(ys) / n163 slope = sum((t - mt) * (y - my) for t, y in zip(ts, ys)) / sum((t - mt) ** 2 for t in ts) / half164 return slope, g165166167def spectrum(g, h):168 w = blackman(len(g))169 y = [a * b for a, b in zip(g, w)]170 m = sum(y) / len(y)171 y = [v - m for v in y]172 step = (OMEGA_HI - OMEGA_LO) / (OMEGA_N - 1)173 ws = [OMEGA_LO + i * step for i in range(OMEGA_N)]174 return y, ws, [power_at(y, h, x) for x in ws]175176177def maxima(y, h, ws, ps):178 out = []179 for i in range(1, len(ps) - 1):180 if ps[i] > ps[i - 1] and ps[i] >= ps[i + 1]:181 x = golden(lambda w: -power_at(y, h, w), ws[i - 1], ws[i + 1])182 out.append((x, power_at(y, h, x)))183 return out184185186def median(xs):187 s = sorted(xs)188 n = len(s)189 return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2190191192CANDIDATES = [193 ("2pi/ln3", 2 * math.pi / math.log(3)),194 ("2pi/ln5", 2 * math.pi / math.log(5)),195 ("2pi/ln9", 2 * math.pi / math.log(9)),196 ("2pi/ln15", 2 * math.pi / math.log(15)),197 ("beat", 2 * math.pi / math.log(3) - 2 * math.pi / math.log(5)),198 ("sum", 2 * math.pi / math.log(3) + 2 * math.pi / math.log(5)),199]200201# READING202203def nearest_comb(w):204 best = None205 for m in range(-COMB, COMB + 1):206 for n in range(-COMB, COMB + 1):207 if m == 0 and n == 0:208 continue209 v = m * OM3 + n * OM5210 if v <= 0:211 continue212 e = abs(v - w) / v * 100213 if best is None or e < best[0]:214 best = (e, m, n, v)215 return best216217218def analyse(us, cs, asked, deg):219 h = us[1] - us[0]220 slope, g = detrend(us, cs, deg)221 y, ws, ps = spectrum(g, h)222 peaks = maxima(y, h, ws, ps)223 med = median(ps)224 loud = [t for t in peaks if t[1] >= RATIO * med]225 cover = sum(2 * TOL * w for w, _ in loud) / (OMEGA_HI - OMEGA_LO)226 top = sorted(peaks, key=lambda t: -t[1])[:4]227 verdict = {}228 for label, target in CANDIDATES:229 if label not in asked:230 continue231 best = min(peaks, key=lambda t: abs(t[0] - target))232 err = abs(best[0] - target) / target * 100233 ratio = best[1] / med234 verdict[label] = (err < TOL * 100 and ratio >= RATIO, err, ratio, best[0], target)235 return slope, max(g) - min(g), len(loud), cover, [(w, p / med) for w, p in top], verdict236237238def reading(us, cs, asked, deg):239 slope, swing, nloud, cover, top, verdict = analyse(us, cs, asked, deg)240 print(" detrend degree %d fitted exponent %.6f swing %.6f maxima above %.0fx median %d chance of a hit by position alone %.3f" % (deg, slope, swing, RATIO, nloud, cover))241 print(" strongest maxima " + " ".join("%.4f (%.3g x)" % t for t in top))242 e, m, n, v = nearest_comb(top[0][0])243 print(" strongest maximum %.6f: nearest m 2pi/ln3 + n 2pi/ln5 with |m|, |n| <= %d is (%d, %d) = %.6f, error %.3f%%" % (top[0][0], COMB, m, n, v, e))244 for label, _ in CANDIDATES:245 if label not in verdict:246 continue247 ok, err, ratio, got, target = verdict[label]248 print(" %-8s predicted %.6f nearest max %.6f error %.3f%% power %.3g x median present %s" % (label, target, got, err, ratio, ok))249 return slope, verdict250251252def report(name, us, cs, asked):253 print(" %s" % name)254 print(" window u in [%.4f, %.4f], span %.4f, %d points, count %s to %s, resolution 2pi/span %.4f" % (us[0], us[-1], us[-1] - us[0], len(us), fmt(cs[0]), fmt(cs[-1]), 2 * math.pi / (us[-1] - us[0])))255 reading(us, cs, asked, 1)256 reading(us, cs, asked, 3)257258259def fmt(c):260 return "%d" % c if float(c).is_integer() else "%.6g" % c261262# VERBS263264def verb_control(M):265 p, pd, q, qd = INDEPENDENT266 print("CONTROL: the pruned walk against direct digit filtering, base %d digits %s and base %d digits %s" % (p, list(pd), q, list(qd)))267 print("L height direct walk nodes agree")268 for L in range(M + 1):269 ref = brute(L, p, pd, q, qd)270 got = []271 nodes = walk(L, p, pd, q, qd, got.append)272 print("%-2d %-13d %-7d %-5d %-6d %s" % (L, p ** L, len(ref), len(got), nodes, ref == got))273 sys.stdout.flush()274 terms = []275 walk(10, p, pd, q, qd, terms.append)276 print("first members of the cell below %d^10: %s" % (p, ", ".join(str(n) for n in terms[1:13])))277278279def pair_counts(L, cell):280 p, pd, q, qd = cell281 t0 = time.time()282 h, counts, nodes = histogram(L, p, pd, q, qd)283 budget = math.log(len(pd), p) + math.log(len(qd), q) - 1284 print(" base %d digits %s against base %d digits %s: budget %.6f rounded up, %d nodes, %d hits, %.1f s" % (p, list(pd), q, list(qd), math.ceil(budget * 1e6) / 1e6, nodes, counts[-1], time.time() - t0))285 sys.stdout.flush()286 return crop(h, counts)287288289def verb_cell(L):290 p, pd, q, qd = INDEPENDENT291 print("CELL: height %d^%d = %d" % (p, L, p ** L))292 print(" criterion: a local maximum of the Blackman periodogram of ln C(e^u) detrended in u, lying within %.0f%% of the prediction and carrying at least %.0fx the median power of the band [%.1f, %.1f]" % (TOL * 100, RATIO, OMEGA_LO, OMEGA_HI))293 print("INDEPENDENT PAIR")294 us, cs = pair_counts(L, INDEPENDENT)295 report("C(N) of the intersection", us, cs, ("2pi/ln3", "2pi/ln5", "2pi/ln15", "beat", "sum"))296 half = len(us) // 2297 report("C(N) on the upper half of the window", us[half:], cs[half:], ("2pi/ln3", "2pi/ln5"))298 print("ONE-BASE CONTROLS")299 hp, cp = one_base_ladder(L, p, p, pd)300 up, vp = crop(hp, cp)301 report("base %d digits %s alone" % (p, list(pd)), up, vp, ("2pi/ln3", "2pi/ln5"))302 hq, cq = one_base_ladder(L, p, q, qd)303 uq, vq = crop(hq, cq)304 report("base %d digits %s alone" % (q, list(qd)), uq, vq, ("2pi/ln3", "2pi/ln5"))305306307def verb_ladder(A, B):308 p, pd, q, qd = INDEPENDENT309 print("LADDER: the criterion at every height from %d^%d to %d^%d, full window, both detrends" % (p, A, p, B))310 print(" a height whose new hits are 0 is a decade of %d containing no member of the cell, on which C is exactly constant" % p)311 print("L hits new span pts loud cover deg 2pi/ln3 err power present 2pi/ln5 err power present")312 prev = None313 for L in range(A, B + 1):314 h, counts, nodes = histogram(L, p, pd, q, qd)315 us, cs = crop(h, counts)316 new = "" if prev is None else "%d" % (counts[-1] - prev)317 prev = counts[-1]318 for deg in (1, 3):319 _, _, nloud, cover, _, v = analyse(us, cs, ("2pi/ln3", "2pi/ln5"), deg)320 a, b = v["2pi/ln3"], v["2pi/ln5"]321 print("%-2d %-9d %-9s %6.2f %5d %4d %6.3f %3d %8.3f%% %7.3g %-8s %8.3f%% %7.3g %-8s" % (L, counts[-1], new, us[-1] - us[0], len(us), nloud, cover, deg, a[1], a[2], a[0], b[1], b[2], b[0]))322 new = ""323 sys.stdout.flush()324325326def top_band(b, k, ds):327 return b ** k, max(ds) * (b ** (k + 1) - 1) // (b - 1)328329330def verb_blocks(L):331 p, pd, q, qd = INDEPENDENT332 print("BLOCKS: the applicability hypothesis of the detector, height %d^%d" % (p, L))333 print(" a member of the base %d design with k+1 digits lies in [%d^k, (%d^(k+1)-1)/2], so the design occupies the fraction ln(%d/2)/ln %d = %.6f of every decade of %d and the rest of the decade is empty" % (p, p, p, p, p, math.log(p / 2) / math.log(p), p))334 print(" the base %d design occupies ln(%d/2)/ln %d = %.6f of every decade of %d" % (q, q, q, math.log(q / 2) / math.log(q), q))335 empty = []336 for j in range(L):337 lo3, hi3 = top_band(p, j, pd)338 hit = False339 i = 0340 while q ** i <= hi3:341 lo5, hi5 = top_band(q, i, qd)342 if max(lo3, lo5) <= min(hi3, hi5):343 hit = True344 break345 i += 1346 if not hit:347 empty.append(j)348 print(" decades [%d^j, %d^(j+1)) whose two design bands do not meet, so the cell has no member there and C is constant across the whole decade, j = %s" % (p, p, ", ".join(str(j) for j in empty)))349 print(" the fraction of the exponents j < %d carrying no member is %d / %d = %.4f" % (L, len(empty), L, len(empty) / L))350 print("BLOCK MODEL: C3(N) C5(N) / N, the two band structures multiplied with no joint arithmetic")351 h = L * math.log(p) / (NGRID - 1)352 model = []353 for i in range(NGRID):354 N = int(math.exp(i * h))355 model.append(digit_count(N, p, pd) * digit_count(N, q, qd) / N if N > 0 else 0)356 um, cm = crop(h, model)357 report("C3(N) C5(N) / N", um, cm, ("2pi/ln3", "2pi/ln5", "2pi/ln15", "beat", "sum"))358359360def verb_collapse(L):361 p, pd, q, qd = DEPENDENT362 print("COLLAPSE: height %d^%d = %d" % (p, L, p ** L))363 print(" the dependent pair collapses to the base 9 design {0,1,3} of exponent log_9 3 = 0.5, one lattice of frequency 2pi/ln 9 whose first harmonic is 2pi/ln 3")364 ud, cd = pair_counts(L, DEPENDENT)365 report("C(N) of the dependent intersection", ud, cd, ("2pi/ln3", "2pi/ln5", "2pi/ln9"))366367368def main():369 args = sys.argv[1:]370 if len(args) == 3 and args[0] == "ladder":371 verb_ladder(int(args[1]), int(args[2]))372 elif len(args) == 2 and args[0] == "blocks":373 verb_blocks(int(args[1]))374 elif len(args) == 2 and args[0] == "collapse":375 verb_collapse(int(args[1]))376 elif len(args) == 2 and args[0] == "control":377 verb_control(int(args[1]))378 elif len(args) == 2 and args[0] == "cell":379 verb_cell(int(args[1]))380 else:381 print("verbs: control M | cell L | collapse M | ladder A B | blocks L")382383384main()