From 2c61e321e555dbd825156aefcb744f6e43991909 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Aug 2024 08:59:52 +0200 Subject: [PATCH] more --- examples/directed/directed_example.py | 157 ------------------ examples/non_abelian/closed/abelian_closed.py | 67 ++++++++ examples/non_abelian/closed/level_spacing.py | 59 +++++++ examples/non_abelian/closed/make_graph.py | 24 +++ examples/non_abelian/closed/so3_closed.py | 111 +++++++++++++ .../non_abelian/closed/so3_closed_uniform.py | 104 ++++++++++++ examples/transfer/boundary_example.py | 78 --------- examples/transfer/transfer_example.py | 119 ------------- netsalt/modes.py | 71 +++++--- netsalt/non_abelian.py | 18 +- netsalt/physics.py | 2 +- 11 files changed, 420 insertions(+), 390 deletions(-) delete mode 100644 examples/directed/directed_example.py create mode 100644 examples/non_abelian/closed/abelian_closed.py create mode 100644 examples/non_abelian/closed/level_spacing.py create mode 100644 examples/non_abelian/closed/make_graph.py create mode 100644 examples/non_abelian/closed/so3_closed.py create mode 100644 examples/non_abelian/closed/so3_closed_uniform.py delete mode 100644 examples/transfer/boundary_example.py delete mode 100644 examples/transfer/transfer_example.py diff --git a/examples/directed/directed_example.py b/examples/directed/directed_example.py deleted file mode 100644 index 83c02b0..0000000 --- a/examples/directed/directed_example.py +++ /dev/null @@ -1,157 +0,0 @@ -import numpy as np -from scipy.sparse import linalg -from tqdm import tqdm -from copy import copy -import matplotlib.pyplot as plt -from netsalt.transfer import get_edge_transfer_matrix, get_static_boundary_flow -from netsalt.quantum_graph import ( - get_total_inner_length, - create_quantum_graph, - construct_weight_matrix, - construct_incidence_matrix, - laplacian_quality, - mode_quality, - construct_laplacian, -) - -from netsalt.modes import mode_on_nodes - -from netsalt.physics import dispersion_relation_linear, set_dispersion_relation -from netsalt.physics import dispersion_relation_resistance -from netsalt.modes import scan_frequencies -from netsalt.plotting import plot_scan -import networkx as nx - - -def make_graph(): - n = 30 - graph = nx.Graph() - graph = nx.cycle_graph(n) - - graph.add_edge(0, 8) - graph.add_edge(1, 4) - graph.add_edge(2, 23) - graph.add_edge(5, 20) - graph.add_edge(8, 15) - 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) - - # graph.add_edge(0, n) - # graph.add_edge(14, n + 1) - # graph.add_edge(16, n + 2) - # graph.add_edge(3, n + 3) - # pos.append(np.array([1.4, 0])) - # pos.append(np.array([-1.4, 0.0])) - # pos.append(np.array([-1.4, -0.5])) - # pos.append(np.array([1.4, 0])) - - return graph, pos - - -if __name__ == "__main__": - np.random.seed(42) - graph, pos = make_graph() - graph_open, pos = make_graph() - graph_reversed, pos = make_graph() - params = { - "open_model": "directed", - "n_workers": 7, - "k_n": 500, - "k_min": 10.0, - "k_max": 15.0, - "alpha_n": 100, - "alpha_min": -0.8, - "alpha_max": 0.8, - } - params["c"] = 1.0 - - nx.draw(graph, pos=pos) - nx.draw_networkx_labels(graph, pos=pos) - - create_quantum_graph(graph, params=copy(params), positions=pos) - set_dispersion_relation(graph, dispersion_relation_linear) - - params["open_model"] = "open" - create_quantum_graph(graph_open, params=copy(params), positions=pos) - set_dispersion_relation(graph_open, dispersion_relation_linear) - - params["open_model"] = "directed_reversed" - create_quantum_graph(graph_reversed, params=copy(params), positions=pos) - set_dispersion_relation(graph_reversed, dispersion_relation_linear) - - qualities = scan_frequencies(graph) - plot_scan(graph, qualities) - - qualities = scan_frequencies(graph_reversed) - plot_scan(graph, qualities) - plt.show() - - k = 1.0 - ks = np.linspace(10.0, 15, 1000) - qs = [] - qs_reversed = [] - qs_open = [] - for k in tqdm(ks): - imk = 0.05 - qs.append(mode_quality([k, imk], graph)) - qs_reversed.append(mode_quality([k, imk], graph_reversed)) - qs_open.append(mode_quality([k, imk], graph_open)) - plt.figure() - plt.plot(ks, qs, label="directed") - plt.plot(ks, qs_reversed, label="directed reversed") - plt.plot(ks, qs_open, label="open") - plt.legend() - plt.yscale("log") - plt.show() - - -def lkj(): - # print(L.dot(f), f.sum()) - # print(out_flow, f) - - plt.figure() - # pp=nx.draw(graph, pos=pos, node_color=p) - nx.draw_networkx_nodes(graph, pos=pos, node_color=p, node_size=50) - f = abs(B.T.dot(p)) - pp = nx.draw_networkx_edges(graph, pos=pos, edge_color=f, width=5) - plt.colorbar(pp) - - r = get_edge_transfer_matrix(_k, graph, input_flow) - L = construct_laplacian(_k, graph) - K = np.append(graph.graph["ks"], graph.graph["ks"]) - BT, B = construct_incidence_matrix(graph) - _input_flow = BT.dot(K * input_flow) - r = np.abs(linalg.spsolve(L, _input_flow)) - - plt.figure() - nx.draw_networkx_nodes(graph, pos=pos, node_size=50, node_color=r) - # pp=nx.draw_networkx_edges(graph, pos=pos, edge_color=ff, width=5) - plt.colorbar(pp) - plt.figure() - plt.plot(r, p, "+") - plt.show() - - plt.figure() - ax = plt.gca() - f = lambda x: np.abs(x) - ax.plot(ks, f(np.array(res)[:, 3]), label="input") - ax.plot(ks, f(np.array(res)[:, 0]), label="output") - # ax.plot(ks,- np.abs(np.array(res)[:, 1]) - np.real(np.array(res)[:, 2]), label="sum") - # ax.plot( - # ks, - # np.real(np.array(res)[:, 1]) + np.real(np.array(res)[:, 2]) + np.real(np.array(res)[:, 0]), - # label="sum2", - # ) - ax.plot(ks, f(np.array(res)[:, 4]), label="input 2") - ax.plot(ks, f(np.array(res)[:, 1]), label="output 2") - ax.plot(ks, f(np.array(res)[:, 5]), label="input 3") - ax.plot(ks, f(np.array(res)[:, 2]), label="output 3") - ax.axvline(_k) - # ax2 = ax.twinx() - # qualities = scan_frequencies(graph) - # plot_scan(graph, qualities, ax=ax2) - # ax.set_ylim(-0.1, 1.1) - plt.legend(loc="upper right") - - plt.show() diff --git a/examples/non_abelian/closed/abelian_closed.py b/examples/non_abelian/closed/abelian_closed.py new file mode 100644 index 0000000..5236a45 --- /dev/null +++ b/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() diff --git a/examples/non_abelian/closed/level_spacing.py b/examples/non_abelian/closed/level_spacing.py new file mode 100644 index 0000000..0aaf77d --- /dev/null +++ b/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() diff --git a/examples/non_abelian/closed/make_graph.py b/examples/non_abelian/closed/make_graph.py new file mode 100644 index 0000000..5603c23 --- /dev/null +++ b/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 diff --git a/examples/non_abelian/closed/so3_closed.py b/examples/non_abelian/closed/so3_closed.py new file mode 100644 index 0000000..3142fe3 --- /dev/null +++ b/examples/non_abelian/closed/so3_closed.py @@ -0,0 +1,111 @@ +import numpy as np +from pathlib import Path +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 make_graph import make_graph +from netsalt.non_abelian import construct_so3_laplacian + + +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], + "laplacian_constructor": construct_so3_laplacian, + } + + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + + create_quantum_graph(graph, params=params, positions=pos) + + chi = np.array([0.0, 0.0, 1.0]) + chis = np.array(len(graph.edges) * [chi]) + for ei, (u, v) in enumerate(graph.edges): + if graph[u][v]["length"] > 0.5: + graph[u][v]["chi"] = np.array([0.0, 1.0, 0.0]) + else: + graph[u][v]["chi"] = np.array([0.0, 0.0, 1.0]) + + ks = np.linspace(5, 7, 2000) + if not Path("so3_qualities_non_uniform.csv").exists(): + res = [] + for k in tqdm(ks): + res.append(mode_quality([k, 0], graph)) + + modes = ks[peak_local_max(1 / (1e-10 + np.array(res))).flatten()] + + np.savetxt("so3_modes_non_uniform.csv", modes) + np.savetxt("so3_qualities_non_uniform.csv", res) + else: + res = np.loadtxt("so3_qualities_non_uniform.csv") + modes = np.loadtxt("so3_modes_non_uniform.csv") + 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("so3_close_scan.pdf") + + modes_df = pd.DataFrame() + modes_df.loc[:, "passive"] = modes + over_graph = oversample_graph(graph, 0.01) + + i = 2 + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_{i}.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(f"so3_close_mode_2_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_x") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_x_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_x") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_x_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_y") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_y_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_y") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_y_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_z_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_z") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_z_{i}.pdf") diff --git a/examples/non_abelian/closed/so3_closed_uniform.py b/examples/non_abelian/closed/so3_closed_uniform.py new file mode 100644 index 0000000..507d69d --- /dev/null +++ b/examples/non_abelian/closed/so3_closed_uniform.py @@ -0,0 +1,104 @@ +import numpy as np +from pathlib import Path +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 make_graph import make_graph +from netsalt.non_abelian import construct_so3_laplacian + + +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], + "laplacian_constructor": construct_so3_laplacian, + } + + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + + create_quantum_graph(graph, params=params, positions=pos) + ks = np.linspace(5, 7, 2000) + if not Path("so3_qualities.csv").exists(): + res = [] + for k in tqdm(ks): + res.append(mode_quality([k, 0], graph)) + + modes = ks[peak_local_max(1 / (1e-10 + np.array(res))).flatten()] + np.savetxt("so3_modes_uniform.csv", modes) + np.savetxt("so3_qualities.csv", res) + else: + res = np.loadtxt("so3_qualities.csv") + modes = np.loadtxt("so3_modes_uniform.csv") + + 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("so3_close_scan_uniform.pdf") + + modes_df = pd.DataFrame() + modes_df.loc[:, "passive"] = modes + i = 4 + over_graph = oversample_graph(graph, 0.01) + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_uniform_{i}.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(f"so3_close_mode_2_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_x") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_x_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_x") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_x_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_y") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_y_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_y") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_y_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") + plt.tight_layout() + plt.savefig(f"so3_close_mode_1_z_uniform_{i}.pdf") + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_z") + plt.tight_layout() + plt.savefig(f"so3_close_mode_2_z_uniform_{i}.pdf") + + #plt.show() diff --git a/examples/transfer/boundary_example.py b/examples/transfer/boundary_example.py deleted file mode 100644 index 70790d4..0000000 --- a/examples/transfer/boundary_example.py +++ /dev/null @@ -1,78 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from netsalt.modes import get_edge_transfer -from netsalt.quantum_graph import create_quantum_graph - - -from netsalt.physics import set_dispersion_relation -from netsalt.physics import dispersion_relation_resistance -import networkx as nx -from netsalt.modes import estimate_boundary_flow - - -def make_graph2(): - n = 30 - graph = nx.Graph() - graph = nx.cycle_graph(n) - - graph.add_edge(0, 8) - graph.add_edge(1, 4) - graph.add_edge(2, 23) - graph.add_edge(5, 20) - graph.add_edge(8, 15) - graph.add_edge(2, 10) - graph.add_edge(3, 11) - graph.add_edge(4, 12) - 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) - - graph.add_edge(0, n) - graph.add_edge(3, n + 1) - graph.add_edge(15, n + 2) - graph.add_edge(17, n + 3) - pos.append(np.array([1.4, 0])) - pos.append(np.array([-1.4, 0.0])) - pos.append(np.array([-1.4, -0.5])) - pos.append(np.array([1.4, 0])) - - return graph, pos - - -if __name__ == "__main__": - - params = {"open_model": "open", "c": 1.0, "R": 1000.0} - np.random.seed(42) - - graph, pos = make_graph2() - - create_quantum_graph(graph, params=params, positions=pos) - set_dispersion_relation(graph, dispersion_relation_resistance) - - n_deg = np.array([len(graph[u]) for u in graph.nodes]) - e_deg = np.array([len(graph[v]) for u, v in graph.edges]) - n_ids = list(np.argwhere(n_deg == 1).flatten()) - e_ids = list(np.argwhere(e_deg == 1).flatten()) - e_ids += list(2 * np.argwhere(e_deg == 1).flatten() + 1) - input_flow = np.zeros(2 * len(graph.edges)) - input_flow[e_ids[0]] = 1.0 - input_flow[e_ids[3]] = 1.0 - - out_flow, _k = estimate_boundary_flow(graph, input_flow) - - res_edges = [] - ks = np.logspace(-5, 2, 1000) - for k in ks: - r_edge = get_edge_transfer(k, graph, input_flow)[e_ids] - res_edges.append(r_edge) - - plt.figure() - ax = plt.gca() - f = lambda x: np.abs(x) - for i in range(len(np.array(res_edges)[0])): - ax.semilogx(ks, f(np.array(res_edges)[:, i]), label=i) - ax.axvline(_k) - ax.set_ylim(0, 2) - plt.legend(loc="upper right") - - plt.show() diff --git a/examples/transfer/transfer_example.py b/examples/transfer/transfer_example.py deleted file mode 100644 index 5b3417a..0000000 --- a/examples/transfer/transfer_example.py +++ /dev/null @@ -1,119 +0,0 @@ -import numpy as np -from scipy.sparse import linalg -from tqdm import tqdm -import matplotlib.pyplot as plt -from netsalt.modes import get_node_transfer, get_edge_transfer -from netsalt.quantum_graph import ( - get_total_inner_length, - create_quantum_graph, - construct_weight_matrix, - construct_incidence_matrix, - laplacian_quality, - mode_quality, - construct_laplacian, -) - -from netsalt.modes import mode_on_nodes - -from netsalt.physics import dispersion_relation_linear, set_dispersion_relation -from netsalt.physics import dispersion_relation_resistance -from netsalt.modes import scan_frequencies -from netsalt.plotting import plot_scan -import networkx as nx - - -def make_graph(): - n = 30 - graph = nx.Graph() - graph = nx.cycle_graph(n) - - # graph.add_edge(0, 8) - # graph.add_edge(1, 4) - # graph.add_edge(2, 23) - # graph.add_edge(5, 20) - # graph.add_edge(8, 15) - 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) - - graph.add_edge(0, n) - graph.add_edge(15, n + 1) - # graph.add_edge(16, n + 2) - # graph.add_edge(3, n + 3) - pos.append(np.array([1.4, 0])) - pos.append(np.array([-1.4, 0.0])) - # pos.append(np.array([-1.4, -0.5])) - # pos.append(np.array([1.4, 0])) - print(len(graph.edges)) - - return graph, pos - - -if __name__ == "__main__": - graph, pos = make_graph() - params = { - "open_model": "open", - "n_workers": 7, - "k_n": 10000, - "k_min": 10.0001, - "k_max": 15.0, - "alpha_n": 20, - "alpha_min": 0.00, - "alpha_max": 0.2, - } - np.random.seed(42) - a = 3 + 0.0 * np.random.uniform(0.0, 1.0, len(graph.edges)) - - e_deg = np.array([len(graph[v]) for u, v in graph.edges]) - # a[e_deg == 1] = 1.0 - params["c"] = 1.0 # np.sqrt(a) - params["R"] = 0.0 # 0 / a - - 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_resistance) - n_deg = np.array([len(graph[u]) for u in graph.nodes]) - n_ids = list(np.argwhere(n_deg == 1).flatten()) - e_ids = list(np.argwhere(e_deg == 1).flatten()) - e_ids += list(2 * np.argwhere(e_deg == 1).flatten() + 1) - input_flow = np.zeros(len(graph.nodes)) - input_flow[n_ids[0]] = 1.0 - - resonance = 2 * np.pi / get_total_inner_length(graph) - res_nodes = [] - res_edges = [] - ks = np.linspace(params["k_min"], params["k_max"], params["k_n"]) - for k in tqdm(ks): - r_node = get_node_transfer(k, graph, input_flow)[n_ids] - r_edge = get_edge_transfer(k, graph, input_flow)[e_ids] - res_nodes.append(r_node) - res_edges.append(r_edge) - - plt.figure() - ax = plt.gca() - f = lambda x: np.abs(x) - plt.axvline(resonance, c="k") - plt.axvline(resonance * 2, c="k") - plt.axvline(resonance * 3, c="k") - - ax.plot(ks, f(np.array(res_nodes)[:, 0]), label="0") - ax.plot(ks, f(np.array(res_nodes)[:, 1]), label="1") - ax.set_xlim(params["k_min"], params["k_max"]) - plt.legend(loc="upper right") - - plt.figure() - ax = plt.gca() - f = lambda x: np.abs(x) - plt.axvline(resonance, c="k") - plt.axvline(resonance * 2, c="k") - plt.axvline(resonance * 3, c="k") - - ax.plot(ks, f(np.array(res_edges)[:, 0]), label="0") - ax.plot(ks, f(np.array(res_edges)[:, 1]), label="1") - ax.plot(ks, f(np.array(res_edges)[:, 2]), label="2") - ax.plot(ks, f(np.array(res_edges)[:, 3]), label="3") - ax.set_xlim(params["k_min"], params["k_max"]) - plt.legend(loc="upper right") - - plt.show() diff --git a/netsalt/modes.py b/netsalt/modes.py index f60e67d..8a79394 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -267,9 +267,10 @@ def mode_on_nodes(mode, graph): to_complex(mode), graph )[0] min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=1, sigma=0, v0=np.ones(laplacian.shape[0]), which="LM" + laplacian, k=4, sigma=0, v0=np.ones(laplacian.shape[0]), which="LM" ) quality_thresh = graph.graph["params"].get("quality_threshold", 1e-2) + print("eigenvalues:", np.abs(min_eigenvalue)) # [np.abs(min_eigenvalue) < quality_thresh]) if abs(min_eigenvalue[0]) > quality_thresh: raise Exception( "Not a mode, as quality is too high: " @@ -285,12 +286,10 @@ def mode_on_nodes(mode, graph): def flux_on_edges(mode, graph): """Compute the flux on each edge (in both directions).""" - node_solution = mode_on_nodes(mode, graph) _, _, B, Winv = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( to_complex(mode), graph, with_k=False ) - return Winv.dot(B).dot(node_solution) @@ -318,23 +317,19 @@ def mean_on_edges(edge_flux, graph, norm_type="abs", mode=None): z[0, 0] = (np.exp(length * (k + np.conj(k))) - 1.0) / (length * (k + np.conj(k))) else: z[0, 0] = 1.0 - z[1, 1] = 1.0 z[0, 1] = (np.exp(length * k) - np.exp(length * np.conj(k))) / ( length * (k - np.conj(k)) ) z[1, 0] = z[0, 1] z[1, 1] = z[0, 0] - mean_edge_solution[ei] = np.abs( + mean_edge_solution[ei] = np.real( edge_flux[2 * ei : 2 * ei + 2].T.dot(z.dot(np.conj(edge_flux[2 * ei : 2 * ei + 2]))) ) if norm_type.startswith("real"): if mode is None: raise Exception("We need the mode for norm_type='real'") - _, _, _, Winv = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( - to_complex(mode), graph, with_k=False - ) if len(edge_flux) > 2 * len(graph.edges): DIM = 3 else: @@ -348,29 +343,55 @@ def _ext(i, shift=1): length = graph.graph["lengths"][ei] z = np.zeros([2 * DIM, 2 * DIM], dtype=np.complex128) if DIM == 3: - w_perp = Ad(2.0 * length * chi).dot(proj_perp(chi)) - w_paral = np.exp(2.0j * length * norm(chi)) * proj_paral(chi) + axis = None if norm_type.endswith("x"): - _z = np.diag([1, 0, 0]) - elif norm_type.endswith("y"): - _z = np.diag([0, 1, 0]) - elif norm_type.endswith("z"): - _z = np.diag([0, 0, 1]) - else: - _z = ( - 0 * chi.T.dot(w_perp) / norm(chi) ** 2 - + w_paral / norm(chi) - - np.eye(3) / norm(chi) - ) / (2.0 * length) + axis = 0 + if norm_type.endswith("y"): + axis = 1 + if norm_type.endswith("z"): + axis = 2 + + def sol(x): + s_paral = np.exp(1.0j * x * norm(chi)) * proj_paral(chi).dot( + edge_flux[_ext(2 * ei)] + ) + s_paral += np.exp(1.0j * (length - x) * norm(chi)) * proj_paral(chi).dot( + edge_flux[_ext(2 * ei + 1)] + ) + s_perp = Ad(x * chi).dot(proj_perp(chi)).dot(edge_flux[_ext(2 * ei)]) + s_perp += ( + Ad((length - x) * chi).dot(proj_perp(chi)).dot(edge_flux[_ext(2 * ei + 1)]) + ) + s = np.real((s_paral + s_perp) ** 2) / length + # find the equation equivalent to this integration!!! + # return np.linalg.norm(s) + if axis is None: + return np.linalg.norm(s) + return s[axis] if DIM == 1: - _z = Winv[_ext(2 * ei), _ext(2 * ei)].toarray() - z[_ext(0), _ext(0)] = z[_ext(1), _ext(1)] = _z + z[0, 0] = z[1, 1] = (np.exp(2.0j * length * chi) - 1) / (2.0j * length * chi) + z[0, 1] = z[1, 0] = np.exp(1.0j * length * chi) + + def sol(x): + return np.real( + ( + np.exp(1.0j * x * chi) * edge_flux[2 * ei] + + np.exp(1.0j * (length - x) * chi) * edge_flux[2 * ei + 1] + ) + ** 2 + / length + ) + + import scipy.integrate as integrate + + res = integrate.quad(sol, 0, length)[0] mean_edge_solution[ei] = np.real( - edge_flux[_ext(2 * ei, shift=2)].T.dot(z.dot(edge_flux[_ext(2 * ei, shift=2)])) + edge_flux[_ext(2 * ei, shift=2)].T.dot(z).dot(edge_flux[_ext(2 * ei, shift=2)]) ) - + # print("check:", mean_edge_solution[ei], res) + mean_edge_solution[ei] = res return mean_edge_solution diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 2ee812f..1675b60 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -26,21 +26,19 @@ def hat(xi): return xi_vec -def proj_perp(chi): +def proj_perp(chi_mat): """Perpendicular projection.""" - return np.eye(3) - proj_paral(chi) + return chi_mat.dot(chi_mat.T) / norm(chi_mat) ** 2 def proj_paral(chi_mat): """Paralell projection.""" - chi_vec = hat(chi_mat) - return np.outer(chi_vec, chi_vec) / np.linalg.norm(chi_vec) ** 2 + return np.eye(3) - proj_perp(chi_mat) def norm(chi_mat): """Norm of chi""" - chi_vec = hat(chi_mat) - return np.linalg.norm(chi_vec) + return np.sqrt(np.trace(0.5 * chi_mat.dot(chi_mat.T))) def Ad(chi_mat): @@ -74,8 +72,8 @@ def _ext(i): expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) expl += ( abelian_scale - * proj_paral(graph.graph["ks"][ei]) * np.exp(1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei])) + * proj_paral(graph.graph["ks"][ei]) ) B[_ext(2 * ei), _ext(u)] = -one @@ -117,9 +115,9 @@ def _ext(i): chi = graph.graph["ks"][ei] length = graph.graph["lengths"][ei] - w_perp = (Ad(2.0 * length * chi) - np.eye(3)).dot(proj_perp(chi)) - w_paral = abelian_scale * (np.exp(2.0j * length * norm(chi)) - np.eye(3)) * proj_paral(chi) - w = w_perp + w_paral + w_perp = Ad(2.0 * length * chi).dot(proj_perp(chi)) + w_paral = abelian_scale * np.exp(2.0j * length * norm(chi)) * proj_paral(chi) + w = w_perp + w_paral - np.eye(3) winv = linalg.inv(w) diff --git a/netsalt/physics.py b/netsalt/physics.py index 7878f62..6bbd4ec 100644 --- a/netsalt/physics.py +++ b/netsalt/physics.py @@ -183,7 +183,7 @@ def update_params_dielectric_constant(graph, params): graph (networkx graph): current graph params (dict): parameters, must include 'gamma_perp' and 'k_a' """ - params["dielectric_constant"] = [graph[u][v]["dielectric_constant"] for u, v in graph.edges] + params["dielectric_constant"] = [graph[u][v].get("dielectric_constant", None) for u, v in graph.edges] def q_value(mode):