degeneracy.py

1import numpy as np2from scipy.sparse import coo_matrix3from scipy.sparse.csgraph import connected_components45CODE = 76BASE = 27TOL = 1e-98LEVELS = list(range(1, 9))9SWEEP = [1e-12, 1e-11, 1e-10, 1e-09, 1e-08, 1e-07, 1e-06, 1e-05]10SECOND = 1.0 - np.sqrt(30.0) / 6.011THIRD = 1.0 - 0.98833242156612NEAR = 1e-6131415def tile(code, q):16    cells = np.array([(code >> i) & 1 for i in range(q * q)], dtype=np.uint8)17    return cells.reshape(q, q)181920def fractal(code, q, level):21    unit = tile(code, q)22    out = unit23    for _ in range(1, level):24        out = np.kron(out, unit)25    return out262728def and_form(grid):29    side = grid.shape[0]30    rows = np.arange(side).reshape(-1, 1)31    cols = np.arange(side).reshape(1, -1)32    return np.array_equal(grid != 0, (rows & cols) == 0)333435def graph(grid):36    flat = grid.reshape(-1)37    index = np.full(flat.size, -1, dtype=np.int64)38    filled = np.flatnonzero(flat)39    index[filled] = np.arange(filled.size)40    index = index.reshape(grid.shape)41    pairs = []42    for axis in range(grid.ndim):43        low = np.take(index, np.arange(index.shape[axis] - 1), axis=axis)44        high = np.take(index, np.arange(1, index.shape[axis]), axis=axis)45        keep = (low >= 0) & (high >= 0)46        pairs.append(np.stack([low[keep], high[keep]], axis=1))47    return filled.size, np.concatenate(pairs)484950def one_component(nodes, edges):51    ones = np.ones(edges.shape[0])52    adj = coo_matrix((ones, (edges[:, 0], edges[:, 1])), shape=(nodes, nodes))53    count, _ = connected_components(adj, directed=False)54    return count == 1555657def laplacian(nodes, edges):58    degree = np.zeros(nodes)59    np.add.at(degree, edges[:, 0], 1.0)60    np.add.at(degree, edges[:, 1], 1.0)61    root = 1.0 / np.sqrt(degree)62    out = np.zeros((nodes, nodes))63    np.fill_diagonal(out, 1.0)64    weight = root[edges[:, 0]] * root[edges[:, 1]]65    out[edges[:, 0], edges[:, 1]] = -weight66    out[edges[:, 1], edges[:, 0]] = -weight67    return out686970def spectrum(level):71    grid = fractal(CODE, BASE, level)72    if not and_form(grid):73        raise SystemExit("kronecker fractal is not the AND set at level %d" % level)74    nodes, edges = graph(grid)75    if not one_component(nodes, edges):76        raise SystemExit("graph is disconnected at level %d" % level)77    return np.sort(np.linalg.eigvalsh(laplacian(nodes, edges)))787980def classes(mu, tol):81    starts = np.concatenate(([0], np.flatnonzero(np.diff(mu) > tol) + 1))82    ends = np.concatenate((starts[1:], [mu.size]))83    counts = ends - starts84    values = np.add.reduceat(mu, starts) / counts85    spreads = mu[ends - 1] - mu[starts]86    return values, counts, spreads878889def locate(values, target):90    at = int(np.argmin(np.abs(values - target)))91    if abs(values[at] - target) < NEAR:92        return at93    return -1949596def multiplicity(values, counts, target):97    at = locate(values, target)98    return int(counts[at]) if at >= 0 else 099100101def fib(n):102    a, b = 0, 1103    for _ in range(n):104        a, b = b, a + b105    return a106107108def main():109    print("object: design code %d, base %d, dimension 2, normalised Laplacian" % (CODE, BASE))110    print("clustering tolerance: %g" % TOL)111    print()112    print("level nodes distinct degenerate repeated mult_1 mult_second max_spread")113    table = {}114    top = None115    for level in LEVELS:116        mu = spectrum(level)117        values, counts, spreads = classes(mu, TOL)118        nodes = mu.size119        big = counts > 1120        distinct = int(values.size)121        degenerate = int(big.sum())122        repeated = float(counts[big].sum()) / nodes123        unit = multiplicity(values, counts, 1.0)124        second = multiplicity(values, counts, SECOND)125        spread = float(spreads[big].max()) if degenerate else 0.0126        print(127            "%d %d %d %d %.4f %d %d %.2e"128            % (level, nodes, distinct, degenerate, repeated, unit, second, spread)129        )130        third_at = locate(values, THIRD)131        table[level] = {132            "nodes": nodes,133            "distinct": distinct,134            "degenerate": degenerate,135            "unit": unit,136            "second": second,137            "second_value": values[locate(values, SECOND)] if second else None,138            "third": int(counts[third_at]) if third_at >= 0 else 0,139            "third_value": float(values[third_at]) if third_at >= 0 else None,140            "spread": spread,141        }142        if level == LEVELS[-1]:143            top = (mu, values, counts)144    print()145146    print("multiplicity of eigenvalue 1, levels 1 to 8")147    print("  measured: %s" % ", ".join(str(table[l]["unit"]) for l in LEVELS))148    print("  3^(L-1):  %s" % ", ".join(str(3 ** (l - 1)) for l in LEVELS))149    print()150151    print("multiplicity of the 1 -/+ sqrt(30)/6 family, levels 3 to 8")152    print("  measured:   %s" % ", ".join(str(table[l]["second"]) for l in LEVELS[2:]))153    print("  3^(L-3)+1:  %s" % ", ".join(str(3 ** (l - 3) + 1) for l in LEVELS[2:]))154    worst = 0.0155    for level in LEVELS[1:]:156        worst = max(worst, abs(table[level]["second_value"] - SECOND))157    print("  1 - sqrt(30)/6 = %.12f" % SECOND)158    print("  1 + sqrt(30)/6 = %.12f" % (2.0 - SECOND))159    print("  worst deviation over levels 2 to 8: %.2e" % worst)160    print()161162    print("third family, levels 4 to 8")163    print("  eigenvalue sought: 1 - %.12f" % (1.0 - THIRD))164    print("  measured at level 8: 1 - %.12f" % (1.0 - table[8]["third_value"]))165    print("  measured:   %s" % ", ".join(str(table[l]["third"]) for l in LEVELS[3:]))166    print("  3^(L-4)+1:  %s" % ", ".join(str(3 ** (l - 4) + 1) for l in LEVELS[3:]))167    _, values8, counts8 = top168    print("  classes at level 8 with multiplicity 82: %d" % int((counts8 == 82).sum()))169    print()170171    print("counting fits")172    print("  distinct:        %s" % ", ".join(str(table[l]["distinct"]) for l in LEVELS))173    print("  2*Fib(2L)+1:     %s" % ", ".join(str(2 * fib(2 * l) + 1) for l in LEVELS))174    print("  degenerate:      %s" % ", ".join(str(table[l]["degenerate"]) for l in LEVELS[1:]))175    print("  2*Fib(2L-3)-1:   %s" % ", ".join(str(2 * fib(2 * l - 3) - 1) for l in LEVELS[1:]))176    print()177178    mu8, values8, counts8 = top179    print("reading the integers off a floating-point spectrum, level 8")180    print("  widest degenerate class spans %.2e" % table[8]["spread"])181    at = locate(values8, 1.0)182    print("  nearest distinct eigenvalue below 1: %.10f" % values8[at - 1])183    print("  nearest distinct eigenvalue above 1: %.10f" % values8[at + 1])184    print("  isolation of the class at 1: %.3f below, %.3f above"185          % (1.0 - values8[at - 1], values8[at + 1] - 1.0))186    print()187188    stop = 3 ** 9189    print("why level 9 stops under this method")190    print("  nodes at level 9: %d" % stop)191    print("  dense matrix before eigensolver workspace: %.2f GB" % (stop * stop * 8 / 1e9))192    print()193194    print("tolerance sweep at level 8")195    print("  tol distinct mult_1 mult_second")196    for tol in SWEEP:197        values, counts, _ = classes(mu8, tol)198        print(199            "  %.0e %d %d %d"200            % (tol, values.size, multiplicity(values, counts, 1.0),201               multiplicity(values, counts, SECOND))202        )203204205main()