nodal.py

1import sys2import time3import numpy as np4from scipy.sparse import coo_matrix5from scipy.sparse.csgraph import connected_components67ZERO = 1e-98MULT = 1e-89SHOW = 1210BINS = np.array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])1112def carpet(level):13    tile = np.ones((3, 3), dtype=np.uint8)14    tile[1, 1] = 015    grid = np.ones((1, 1), dtype=np.uint8)16    for _ in range(level):17        grid = np.kron(grid, tile)18    assert int(grid.sum()) == 8 ** level19    return grid2021def edges(grid):22    index = -np.ones(grid.shape, dtype=np.int64)23    filled = np.argwhere(grid == 1)24    index[filled[:, 0], filled[:, 1]] = np.arange(len(filled))25    right = (grid[:, :-1] == 1) & (grid[:, 1:] == 1)26    down = (grid[:-1, :] == 1) & (grid[1:, :] == 1)27    a = np.concatenate([index[:, :-1][right], index[:-1, :][down]])28    b = np.concatenate([index[:, 1:][right], index[1:, :][down]])29    return len(filled), a, b3031def laplacian(n, a, b):32    lap = np.zeros((n, n))33    lap[a, b] = -1.034    lap[b, a] = -1.035    lap[np.arange(n), np.arange(n)] = -lap.sum(axis=1)36    return lap3738def nodal(n, a, b, v):39    sign = np.where(v > ZERO, 1, np.where(v < -ZERO, -1, 0))40    keep = (sign[a] == sign[b]) & (sign[a] != 0)41    graph = coo_matrix((np.ones(int(keep.sum())), (a[keep], b[keep])), shape=(n, n))42    parts, _ = connected_components(graph, directed=False)43    zeros = int((sign == 0).sum())44    return parts - zeros, zeros4546def clusters(values):47    n = len(values)48    last = np.zeros(n, dtype=np.int64)49    mult = np.zeros(n, dtype=np.int64)50    start = 051    for i in range(1, n + 1):52        if i == n or values[i] - values[i - 1] > MULT:53            last[start:i] = i54            mult[start:i] = i - start55            start = i56    return last, mult5758def report(name, n, a, b, values, vectors):59    k = np.arange(1, n + 1)60    nu = np.zeros(n, dtype=np.int64)61    zeros = np.zeros(n, dtype=np.int64)62    for i in range(n):63        nu[i], zeros[i] = nodal(n, a, b, vectors[:, i])64    last, mult = clusters(values)65    ratio = nu / k66    hist = np.histogram(np.minimum(ratio, 1.0), bins=BINS)[0]67    violations = nu > k68    degenerate_k = int((mult > 1).sum())69    degenerate_classes = len(set(last[mult > 1].tolist()))70    classes = len(set(last.tolist()))71    tail = ratio[1:]72    print(f"{name}: {n} nodes, {len(a)} edges, lambda_2 = {values[1]:.10f}, lambda_max = {values[-1]:.10f}")73    print(f"  nu_k, k = 1..{SHOW}: {nu[:SHOW].tolist()}")74    print(f"  r_k, k = 1..{SHOW}: {mult[:SHOW].tolist()}")75    print(f"  zero cells, k = 1..{SHOW}: {zeros[:SHOW].tolist()}")76    print(f"  max nu_k = {nu.max()} at k = {int(k[nu.argmax()])}; over k >= 2 max nu_k/k = {tail.max():.6f} at k = {int(tail.argmax()) + 2}, mean nu_k/k = {tail.mean():.6f}")77    print(f"  equality nu_k = k at k = {k[nu == k][:SHOW].tolist()}, {int((nu == k).sum())} times")78    print(f"  Courant nu_k <= k: {int((~violations).sum())} of {n}, fraction {(~violations).mean():.6f}; violations {int(violations.sum())}, of which at degenerate k {int((mult[violations] > 1).sum())}")79    print(f"  DGLS nu_k <= k + r_k - 1: {int((nu <= last).sum())} of {n}, fraction {(nu <= last).mean():.6f}")80    print(f"  eigenvalue classes {classes}, degenerate classes {degenerate_classes}, degenerate indices {degenerate_k} of {n}, fraction {degenerate_k / n:.6f}, max multiplicity {int(mult.max())}")81    print(f"  eigenvectors with a zero cell: {int((zeros > 0).sum())}")82    top = vectors[:, -1]83    sign = np.sign(top) * (np.abs(top) > ZERO)84    print(f"  top eigenvector: zero cells {int(zeros[-1])}, min |v| = {np.abs(top).min():.3e}, edges joining two nonzero cells of one sign: {int(((sign[a] == sign[b]) & (sign[a] != 0)).sum())}")85    print(f"  nu_k/k in [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1]: {hist.tolist()}, above 1: {int((ratio > 1).sum())}")86    return nu, mult, last8788def mixing(n, a, b, vectors, mult):89    found = []90    for i in np.flatnonzero(mult == 2)[::2]:91        u, v = vectors[:, i], vectors[:, i + 1]92        counts = [nodal(n, a, b, w)[0] for w in (u, v, (u + v) / np.sqrt(2), (u - v) / np.sqrt(2))]93        if len(set(counts)) > 1:94            found.append((i + 1, counts))95    print(f"  double eigenvalues where the two returned vectors, their sum and their difference disagree in nu: {len(found)} of {int((mult == 2).sum()) // 2}; first at k = {found[0][0] if found else None}, counts {found[0][1] if found else None}")9697def separable(side):98    i = np.arange(side)99    lam = 2.0 - 2.0 * np.cos(np.pi * np.arange(side) / side)100    cos = np.cos(np.pi * np.outer(np.arange(side), i + 0.5) / side)101    items = sorted(((lam[p] + lam[q], p, q) for p in range(side) for q in range(side)))102    vectors = np.empty((side * side, side * side))103    expect = np.empty(side * side, dtype=np.int64)104    for col, (_, p, q) in enumerate(items):105        vectors[:, col] = np.outer(cos[p], cos[q]).ravel()106        expect[col] = (p + 1) * (q + 1)107    return np.array([t[0] for t in items]), vectors, expect108109def run_carpet(level):110    t = time.time()111    n, a, b = edges(carpet(level))112    values, vectors = np.linalg.eigh(laplacian(n, a, b))113    assert connected_components(coo_matrix((np.ones(len(a)), (a, b)), shape=(n, n)), directed=False)[0] == 1114    nu, mult, last = report(f"carpet L = {level}", n, a, b, values, vectors)115    mixing(n, a, b, vectors, mult)116    print(f"  {time.time() - t:.1f} s")117    print()118119def run_grid(side):120    t = time.time()121    n, a, b = edges(np.ones((side, side), dtype=np.uint8))122    values, vectors = np.linalg.eigh(laplacian(n, a, b))123    nu, mult, last = report(f"grid {side} x {side}, numpy basis", n, a, b, values, vectors)124    mixing(n, a, b, vectors, mult)125    svalues, svectors, expect = separable(side)126    assert np.abs(svalues - values).max() < 1e-9127    snu, smult, slast = report(f"grid {side} x {side}, separable basis", n, a, b, svalues, svectors)128    print(f"  separable nu_k equals (p + 1)(q + 1) for all k: {bool((snu == expect).all())}")129    print(f"  {time.time() - t:.1f} s")130    print()131132def main():133    t = time.time()134    print(f"domain: carpet L = 3, 4 on 4-neighbour adjacency, combinatorial Laplacian D - A, numpy.linalg.eigh; zero tolerance {ZERO:g}, multiplicity tolerance {MULT:g}; controls grid 22 x 22 and 64 x 64")135    print("nu_k counts strong nodal domains of the k-th returned eigenvector; on a degenerate eigenvalue the count depends on the basis, so every nu_k below is a fact about numpy's returned basis only")136    print()137    for level in (3, 4):138        run_carpet(level)139    for side in (22, 64):140        run_grid(side)141    print(f"total {time.time() - t:.1f} s")142143if __name__ == "__main__":144    main()