fragility.py

1import collections2import itertools3import sys4import time5from fractions import Fraction67import numpy as np89SPAN = 1310SIZE = 611ORDERS = (2, 3)12DIAL = ([0, 1, 2], [0, 1, 4, 6], [0, 1, 2, 3, 7], [0, 1, 4, 7, 9], [0, 1, 2, 3, 7, 11])13CAP = 1700014LAGS = 8_000_0001516def coarray(S):17    return {x - y for x in S for y in S}1819def weights(S):20    w = collections.Counter()21    for x in S:22        for y in S:23            w[x - y] += 124    return w2526def design(G, b, r):27    F = [0]28    for i in range(r):29        F = sorted({f + g * b ** i for f in F for g in G})30    return F3132def digits(s, b, r):33    return [(s // b ** i) % b for i in range(r)]3435def brute(S):36    D = coarray(S)37    return {s for s in S if coarray([x for x in S if x != s]) != D}3839def fast(S):40    pairs = collections.defaultdict(list)41    for x in S:42        for y in S:43            if x != y and len(pairs[x - y]) < 3:44                pairs[x - y].append((x, y))45    E = set()46    for p in pairs.values():47        if len(p) == 1:48            E.update(p[0])49        elif len(p) == 2:50            E.update(set(p[0]) & set(p[1]))51    return E5253def paired(G):54    w = weights(G)55    return {g for g in G if any(w[g - h] == 1 for h in G if h != g)}5657def holefree(G):58    a = max(G) - min(G)59    return coarray(G) == set(range(-a, a + 1))6061def generators():62    for a in range(1, SPAN + 1):63        for L in range(2, SIZE + 1):64            for mid in itertools.combinations(range(1, a), L - 2):65                G = [0, *mid, a]66                if tuple(a - x for x in reversed(G)) < tuple(G):67                    continue68                if holefree(G):69                    yield G7071def exact():72    t0 = time.time()73    gens = list(generators())74    cases = bad = tight_bad = 075    loose = []76    for G in gens:77        L, a = len(G), max(G)78        M = 2 * a + 179        U = paired(G)80        E = brute(G)81        assert U <= E82        for r in ORDERS:83            F = design(G, M, r)84            assert len(F) == L ** r and coarray(F) == set(range(-(M ** r - 1) // 2, (M ** r + 1) // 2))85            got = brute(F)86            want = {s for s in F if all(d in U for d in digits(s, M, r))}87            cases += 188            bad += got != want89            tight_bad += (len(got) == len(E) ** r) != (E == U)90            assert fast(F) == got91        if E != U:92            loose.append((G, L, len(E), len(U)))93    econ = [x for x in loose if x[2] == x[1]]94    print(f"generators {len(gens)} cases {cases} mismatches {bad} tightness-criterion failures {tight_bad}")95    print(f"loose generators (F_G > u/L) {len(loose)}, of them maximally economic (F_G = 1) {len(econ)}")96    byL = collections.Counter(x[1] for x in econ)97    print(f"  maximally economic loose by L {dict(sorted(byL.items()))}; card E(G) - u over the loose {sorted({e - u for _, _, e, u in loose})}")98    for G, L, e, u in loose:99        print(f"  {G} L {L} essential {e} paired {u}{' economic' if e == L else ''}")100    for name, G in (("ula3", [0, 1, 2]), ("mra4", [0, 1, 4, 6]), ("ula4", [0, 1, 2, 3])):101        M = 2 * max(G) + 1102        row = [(len(brute(design(G, M, r))), len(G) ** r) for r in (1, 2, 3)]103        print(f"  {name} {G} M {M} (essential, sensors) at r = 1, 2, 3: {row}")104    print(f"exact {time.time() - t0:.1f}s")105106def law(row, u):107    if all(e == u ** (i + 1) for i, e in enumerate(row) if i):108        return "u^r"109    for k in (2, 3):110        tail = row[k - 1:]111        if len(tail) < 3:112            continue113        d0, d1 = tail[1] - tail[0], tail[2] - tail[1]114        lam = Fraction(d1, d0) if d0 else Fraction(1)115        c = tail[1] - lam * tail[0]116        if lam >= 1 and all(tail[j + 1] == lam * tail[j] + c for j in range(len(tail) - 1)):117            return f"e_(r+1) = {lam} e_r + {c} from r = {k}, fitted on 2 steps, checked on {len(tail) - 3}"118    return "no affine law"119120def ess(F, w, A):121    F = np.asarray(F, dtype=np.int64)122    inside = np.zeros(A + 1, dtype=bool)123    inside[F] = True124    out = np.zeros(len(F), dtype=bool)125    ch = max(1, 4_000_000 // len(F))126    for i in range(0, len(F), ch):127        s = F[i:i + ch, None]128        t = s - F[None, :]129        W = w[t + A]130        k = s + t131        two = (W == 2) & (t != 0) & (k >= 0) & (k <= A) & inside[np.clip(k, 0, A)]132        out[i:i + ch] = ((W == 1) & (t != 0)).any(axis=1) | two.any(axis=1)133    return set(F[out].tolist())134135def dial():136    t0 = time.time()137    for G in DIAL:138        a, L, u = max(G), len(G), len(paired(G))139        wG = weights(G)140        for b in range(a + 1, 2 * a + 2):141            row, lags, w, A = [], 0, np.ones(1, dtype=np.int64), 0142            for r in range(1, 12):143                An = a * (b ** r - 1) // (b - 1)144                if L ** r > CAP or 2 * An + 1 > LAGS:145                    break146                wn = np.zeros(2 * An + 1, dtype=np.int64)147                for d, m in wG.items():148                    lo = An + d * b ** (r - 1) - A149                    wn[lo:lo + 2 * A + 1] += m * w150                w, A = wn, An151                F = design(G, b, r)152                assert len(F) == L ** r and F[-1] == A and (w > 0).all() and w.sum() == len(F) ** 2153                E = ess(F, w, A)154                assert len(F) > 250 or (E == brute(F) == fast(F) and weights(F) == {t - A: int(x) for t, x in enumerate(w) if x})155                row.append(len(E))156                lags = 2 * A + 1157            print(f"  {G} b {b} e_r {row} {law(row, u)}; {L ** len(row)} sensors, {lags} lags at the last r")158    print(f"dial {time.time() - t0:.1f}s")159160if __name__ == "__main__":161    verbs = sys.argv[1:] or ["exact", "dial"]162    for v in verbs:163        {"exact": exact, "dial": dial}[v]()