1"""Figures from results/*.npz (written by experiments.py). Run: python plots.py"""
3from __future__ import annotations
5import os
7import matplotlib
9matplotlib.use("Agg")
10import matplotlib.pyplot as plt # noqa: E402
11import numpy as np # noqa: E402
13from model import spectral_entropy # noqa: E402
15HERE = os.path.dirname(os.path.abspath(__file__))
16RES = os.path.join(HERE, "results")
19def scan():
20 d = np.load(os.path.join(RES, "scan.npz"))
21 energies = d["energies"]
22 fig, ax = plt.subplots(1, 3, figsize=(17, 4.6))
23 rows = []
24 for mu, col in ((0.0, "C0"), (1.0, "C3")):
25 lam = np.array([np.median(d[f"mu{mu:g}_E{i}_lyap"]) for i in range(len(energies))])
26 nfin = np.array([np.median(d[f"mu{mu:g}_E{i}_neff"][-1]) for i in range(len(energies))])
27 ax[0].loglog(energies, lam, "o-", color=col, label=f"mu = {mu:g}")
28 ax[1].semilogx(energies, nfin, "o-", color=col, label=f"mu = {mu:g}")
29 for i in range(0, len(energies), 3):
30 ne = np.median(d[f"mu{mu:g}_E{i}_neff"], axis=1)
31 t = 20.0 * (np.arange(len(ne)) + 1) # records every 20 time units
32 ax[2].semilogx(t, ne, color=col, alpha=0.3 + 0.7 * i / len(energies), lw=1)
33 rows.append((mu, lam, nfin))
34 ax[0].set(xlabel="strike energy E (N = 32, beta = 1)", ylabel="Lyapunov exponent (4000 time units)",
35 title="Rate of separation")
36 ax[1].set(xlabel="strike energy E", ylabel="effective fraction of modes at t = 4000",
37 title="Equipartition (about 0.7 for random equipartition)")
38 ax[2].set(xlabel="time", ylabel="effective fraction of modes", title="Spreading over time (darker: more energy)")
39 for a in ax:
40 a.grid(True, which="both", alpha=0.3)
41 ax[0].legend()
42 ax[1].legend()
43 fig.tight_layout()
44 fig.savefig(os.path.join(RES, "scan.png"), dpi=100)
45 return energies, rows
48def cascade():
49 d = np.load(os.path.join(RES, "cascade.npz"))
50 Es = (0.5, 2.0, 10.0)
51 fig, ax = plt.subplots(2, 3, figsize=(17, 7.5), sharex=True)
52 for i, mu in enumerate((0.0, 1.0)):
53 for j, E in enumerate(Es):
54 k = f"mu{mu:g}_E{E:g}"
55 t, Ek, om = d[k + "_t"], d[k + "_Ek"], d[k + "_omega"]
56 frac = Ek / Ek.sum(1, keepdims=True)
57 im = ax[i, j].pcolormesh(t, np.arange(1, Ek.shape[1] + 1), np.log10(frac.T + 1e-12),
58 vmin=-6, vmax=0, cmap="magma", shading="auto", rasterized=True)
59 ax[i, j].set(title=f"mu = {mu:g}, E = {E:g}", ylabel="mode k" if j == 0 else None)
60 if i == 1:
61 ax[i, j].set_xlabel("time")
62 fig.colorbar(im, ax=ax, label="log10 share of energy in mode k", shrink=0.8)
63 fig.savefig(os.path.join(RES, "cascade.png"), dpi=100, bbox_inches="tight")
66def damped(mu):
67 d = np.load(os.path.join(RES, "damped.npz"))
68 E0s = (0.5, 5.0, 30.0)
69 fig, ax = plt.subplots(4, 3, figsize=(17, 14))
70 summary = []
71 for j, E0 in enumerate(E0s):
72 k = f"mu{mu:g}_E{E0:g}"
73 t, Ek, E, logsum = d[k + "_t"], d[k + "_Ek"], d[k + "_E"], d[k + "_logsum"]
74 neff = spectral_entropy(Ek) # (T, M)
75 ax[0, j].semilogy(t, E[:, 0], "k")
76 ax[0, j].set(title=f"mu = {mu:g}, strike energy E0 = {E0:g}", ylabel="energy")
77 ax[1, j].plot(t, neff, lw=0.6)
78 ax[1, j].set(ylabel="effective fraction of modes\n(16 strikes differing by 1e-10)", ylim=(0, 1))
79 frac = Ek[-1] / Ek[-1].sum(1, keepdims=True) # (M, N)
80 ks = np.arange(1, frac.shape[1] + 1)
81 for m in range(frac.shape[0]):
82 ax[2, j].semilogy(ks, frac[m], lw=0.7)
83 ax[2, j].set(xlabel="mode k", ylabel="final share of energy", ylim=(1e-8, 1))
84 w = max(1, int(round(100 / (t[1] - t[0]))))
85 loc = (logsum[w:, 0] - logsum[:-w, 0]) / (t[w:] - t[:-w])
86 ax[3, j].plot(t[w:], loc, "k", lw=0.8)
87 ax[3, j].axhline(0, color="0.6", lw=0.5)
88 ax[3, j].set(xlabel="time", ylabel="local Lyapunov exponent (100-unit windows)")
89 top = frac.mean(0) > 1e-3
90 spread = np.log10(frac[:, top]).std(0)
91 summary.append((E0, float(np.median(spread)), float(spread.max()), int(top.sum())))
92 fig.tight_layout()
93 fig.savefig(os.path.join(RES, f"damped_mu{mu:g}.png"), dpi=95)
94 return summary
97if __name__ == "__main__":
98 if os.path.exists(os.path.join(RES, "scan.npz")):
99 energies, rows = scan()
100 for mu, lam, nfin in rows:
101 print(f"mu {mu:g}: " + " ".join(f"E {E:.2g}: lam {l:.4f} n {n:.2f}" for E, l, n in zip(energies, lam, nfin)))
102 if os.path.exists(os.path.join(RES, "cascade.npz")):
103 cascade()
104 if os.path.exists(os.path.join(RES, "damped.npz")):
105 for mu in (0.0, 1.0):
106 for E0, med, mx, n in damped(mu):
107 print(f"damped mu {mu:g}, E0 {E0:g}: spread of final log10 mode shares across strikes over the {n} modes "
108 f"holding > 0.1%: median {med:.3f}, max {mx:.3f} (decades)")
109 print("done")