more

examples/non_abelian/closed/abelian_closed.py

@@ -0,0 +1,67 @@
+import numpy as np
+import pandas as pd
+from tqdm import tqdm
+import matplotlib.pyplot as plt
+from skimage.feature import peak_local_max
+import networkx as nx
+
+from netsalt.quantum_graph import oversample_graph, create_quantum_graph, mode_quality
+from netsalt.plotting import plot_single_mode
+from netsalt.physics import dispersion_relation_linear, set_dispersion_relation
+
+from make_graph import make_graph
+
+
+if __name__ == "__main__":
+    graph, pos = make_graph()
+    params = {
+        "open_model": "open",
+        "n_workers": 7,
+        "k_n": 2000,
+        "k_min": 0.00001,
+        "k_max": 5.2,
+        "alpha_n": 20,
+        "alpha_min": 0.00,
+        "alpha_max": 0.2,
+        "quality_threshold": 1e-3,
+        "c": len(graph.edges) * [1.0],
+    }
+
+    nx.draw(graph, pos=pos)
+    nx.draw_networkx_labels(graph, pos=pos)
+
+    create_quantum_graph(graph, params=params, positions=pos)
+    set_dispersion_relation(graph, dispersion_relation_linear)
+
+    res = []
+    ks = np.linspace(5, 7, 2000)
+    for k in tqdm(ks):
+        res.append(mode_quality([k, 0], graph))
+
+    modes = ks[peak_local_max(1 / (1e-10 + np.array(res))).flatten()]
+    print(modes)
+    plt.figure(figsize=(4, 2))
+    for mode in modes:
+        plt.axvline(mode, c="k")
+    plt.semilogy(ks, res, "-")
+    plt.axis([ks[0], ks[-1], 1e-3, 1])
+    plt.xlabel("wavenumber")
+    plt.ylabel("mode quality")
+    plt.tight_layout()
+    plt.savefig("close_scan_abelian.pdf")
+
+    modes_df = pd.DataFrame()
+    modes_df.loc[:, "passive"] = modes
+    over_graph = oversample_graph(graph, 0.01)
+    over_graph.graph["params"]["c"] = len(over_graph.edges) * [1.0]
+    set_dispersion_relation(over_graph, dispersion_relation_linear)
+    plt.figure(figsize=(4, 3))
+    plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real")
+    plt.tight_layout()
+    plt.savefig("close_mode_1_abelian.pdf")
+
+    plt.figure(figsize=(4, 3))
+    plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real")
+    plt.tight_layout()
+    plt.savefig("close_mode_2_abelian.pdf")
+    plt.show()

examples/non_abelian/closed/level_spacing.py

@@ -0,0 +1,59 @@
+import numpy as np
+from functools import partial
+from tqdm import tqdm
+from multiprocessing import Pool
+import matplotlib.pyplot as plt
+from skimage.feature import peak_local_max
+
+from netsalt.quantum_graph import create_quantum_graph, mode_quality
+from netsalt.physics import dispersion_relation_linear, set_dispersion_relation
+
+from make_graph import make_graph
+
+
+def _qual(k, graph=None):
+    return mode_quality([k, 0], graph)
+
+
+if __name__ == "__main__":
+    graph, pos = make_graph()
+    k_max = 200
+    k_res = 500
+    params = {"open_model": "open", "quality_threshold": 1e-3, "c": len(graph.edges) * [1.0]}
+
+    create_quantum_graph(graph, params=params, positions=pos)
+    set_dispersion_relation(graph, dispersion_relation_linear)
+
+    res = []
+    ks = np.linspace(10, k_max, k_res * k_max)
+
+    with Pool() as pool:
+        res = list(tqdm(pool.imap(partial(_qual, graph=graph), ks), total=len(ks)))
+
+    modes = ks[sorted(peak_local_max(1 / (1e-10 + np.array(res))).flatten())]
+    print(len(modes))
+
+    plt.figure(figsize=(20, 2))
+    for mode in modes:
+        plt.axvline(mode, c="k")
+
+    plt.semilogy(ks, res, "-")
+    plt.axis([ks[0], ks[-1], 1e-3, 1])
+    plt.xlabel("wavenumber")
+    plt.ylabel("mode quality")
+    plt.tight_layout()
+
+    modes_inter = np.diff(modes)
+    mean_modes_inter = np.mean(modes_inter)
+
+    plt.figure(figsize=(5, 3))
+    plt.hist(modes_inter / mean_modes_inter, bins=40, histtype="step", density=True, label="data")
+    s = np.linspace(0, 4, 100)
+    plt.plot(s, np.pi * s / 2 * np.exp(-np.pi / 4 * s**2), label="GOE")
+    plt.plot(s, np.exp(-s), label="Poisson")
+    plt.xlabel("s")
+    plt.ylabel("P(s)")
+    plt.legend(loc="best")
+    plt.tight_layout()
+    plt.savefig("closed_level_spacing_abelian.pdf")
+    plt.show()

examples/non_abelian/closed/make_graph.py

@@ -0,0 +1,24 @@
+import networkx as nx
+import numpy as np
+
+
+def make_graph(with_leads=False):
+    n = 30
+    graph = nx.Graph()
+    graph = nx.cycle_graph(n)
+
+    graph.add_edge(2, 8)
+    graph.add_edge(27, 16)
+    graph.add_edge(16, 10)
+    x = np.linspace(0, 2 * np.pi * (1 - 1.0 / (len(graph) - 1)), len(graph))
+    pos = np.array([np.cos(x), np.sin(x)]).T
+    pos = list(pos)
+    if with_leads:
+        graph.add_edge(0, n)
+        graph.add_edge(14, n + 1)
+        graph.add_edge(16, n + 2)
+        pos.append(np.array([1.4, 0]))
+        pos.append(np.array([-1.4, 0.3]))
+        pos.append(np.array([-1.4, -0.3]))
+
+    return graph, pos