# Eigenvalue art: the eigenvalues of a complex matrix that depends on two angles (alpha, beta),
# drawn for many random angle pairs. Two families: "rand" (random sparse matrices) and "ssh" (a ring with gain and loss).
# Needs: pip install numpy matplotlib
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import hsv_to_rgb

FAMILY = "rand"      # "rand" or "ssh"
N = 100_000          # number of angle pairs; try 10_000 first
SIZE = 1000
PAPER = (0.957, 0.941, 0.902)

# Colour scheme: "paper" (translucent pastel ink on a light background), "script" (hue = angle of the first
# parameter, saturation 0.8, brightness 0.2 ... 1 from the second angle, on black) or "gradient" (your own colours below)
SCHEME = "paper"
GRADIENT = [(0.0, "#05030f"), (0.25, "#3a0f6b"), (0.5, "#c2185b"), (0.75, "#ff8a3d"), (1.0, "#fff3c4")]   # (position 0..1, colour)
GRADIENT_FOR = "density"   # "density" (dark = few roots) or "source" (the gradient is looped over the colour source; make the last colour equal the first)

rng = np.random.default_rng(7)
al, be = rng.random(N) * 2 * np.pi, rng.random(N) * 2 * np.pi

if FAMILY == "rand":
    # M(alpha, beta) = A0 + e^{i alpha} A1 + e^{i beta} A2 with three fixed random sparse 6 x 6 matrices
    n = 6
    A = [(rng.standard_normal((n, n)) + 1j * rng.standard_normal((n, n))) * (rng.random((n, n)) < 0.6) for _ in range(3)]
    M = A[0] + np.exp(1j * al)[:, None, None] * A[1] + np.exp(1j * be)[:, None, None] * A[2]
    HALF = 3.4
else:
    # non-Hermitian SSH ring: hopping v inside a cell, w e^{i beta} between cells, gain/loss +-i*gamma, twist alpha
    cells, v, w, gam = 12, 1.0, 1.6, 0.8
    n = 2 * cells
    M = np.zeros((N, n, n), dtype=complex)
    for j in range(n - 1):
        h = v if j % 2 == 0 else w * np.exp(1j * be)
        M[:, j, j + 1] = h
        M[:, j + 1, j] = h
    M[:, 0, n - 1] += w * np.exp(1j * (be - al))
    M[:, n - 1, 0] += w * np.exp(1j * (be + al))
    M[:, np.arange(n), np.arange(n)] += np.where(np.arange(n) % 2 == 0, 1j * gam, -1j * gam)
    HALF = v + w + gam + 0.2

ev = np.linalg.eigvals(M)                        # shape (N, n)
hue = np.repeat(al / (2 * np.pi), n)             # hue = the angle alpha
x, y = ev.real.ravel(), ev.imag.ravel()

kw = dict(bins=SIZE, range=[[-HALF, HALF], [-HALF, HALF]])
cnt, _, _ = np.histogram2d(y, x, **kw)
cs, _, _ = np.histogram2d(y, x, weights=np.cos(2 * np.pi * hue), **kw)
sn, _, _ = np.histogram2d(y, x, weights=np.sin(2 * np.pi * hue), **kw)
un, _, _ = np.histogram2d(y, x, weights=np.repeat(be / (2 * np.pi), n), **kw)   # mean of the second angle, for the script colours

dens = np.log1p(cnt / max(np.percentile(cnt[cnt > 0], 6), 1e-9))
v_ = np.clip(dens / np.percentile(dens[cnt > 0], 99.7), 0, 1)
alpha = np.clip(0.9 * v_**0.55, 0, 1)[..., None]
h = (np.arctan2(sn, cs) / (2 * np.pi)) % 1
conc = np.clip(np.hypot(cs, sn) / np.maximum(cnt, 1) * 1.5 + 0.3, 0, 1)
ink = hsv_to_rgb(np.dstack([h, 0.78 * conc, np.full_like(h, 0.93)]))
if SCHEME == "paper":
    img = np.array(PAPER) * (1 - alpha) + ink * alpha
elif SCHEME == "script":
    m = np.maximum(cnt, 1)
    bright = 0.2 + 0.8 * np.clip(un / m, 0, 1)
    col = hsv_to_rgb(np.dstack([h, np.full_like(h, 0.8), bright]))
    a2_ = np.clip(v_**0.6, 0, 1)[..., None]
    img = col * a2_                                  # on black
else:
    from matplotlib.colors import to_rgb
    pos = np.array([p for p, _ in GRADIENT]); cols = np.array([to_rgb(c) for _, c in GRADIENT])
    lut = lambda t: np.stack([np.interp(t, pos, cols[:, k]) for k in range(3)], axis=-1)
    if GRADIENT_FOR == "density":
        img = lut(v_)
    else:
        img = lut(h) * alpha + np.array(PAPER) * (1 - alpha)

plt.figure(figsize=(8, 8))
plt.imshow(img[::-1], extent=[-HALF, HALF, -HALF, HALF])
plt.axis("off")
plt.savefig("eigen_art.png", dpi=150, bbox_inches="tight", pad_inches=0)
print("saved eigen_art.png")