diff --git a/.pylintrc b/.pylintrc index 67c4326..cbdd470 100644 --- a/.pylintrc +++ b/.pylintrc @@ -29,7 +29,7 @@ max-locals=20 # Maximum number of return / yield for function / method body max-returns=6 # Maximum number of branch for function / method body -max-branches=12 +max-branches=20 # Maximum number of statements in function / method body max-statements=50 # Maximum number of parents for a class (see R0901). 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_random.py b/examples/non_abelian/closed/so3_closed_random.py new file mode 100644 index 0000000..fd7d7bc --- /dev/null +++ b/examples/non_abelian/closed/so3_closed_random.py @@ -0,0 +1,114 @@ +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]) + random = [np.random.normal(0, 1, 3) for _ in range(len(graph.edges))] + chis = np.array(len(graph.edges) * [chi]) + for ei, (u, v) in enumerate(graph.edges): + graph[u][v]["chi"] = random[ei] + graph[u][v]["chi"] /= np.linalg.norm(graph[u][v]["chi"]) + + # 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) + + plt.figure(figsize=(4, 3)) + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig("so3_close_mode_1.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("so3_close_mode_1_x.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("so3_close_mode_1_y.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("so3_close_mode_1_z.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("so3_close_mode_2.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("so3_close_mode_2_x.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("so3_close_mode_2_y.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("so3_close_mode_2_z.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..4c66f2f --- /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 + 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("so3_close_mode_1_uniform.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("so3_close_mode_1_x_uniform.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("so3_close_mode_1_y_uniform.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("so3_close_mode_1_z_uniform.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("so3_close_mode_2_uniform.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("so3_close_mode_2_x_uniform.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("so3_close_mode_2_y_uniform.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("so3_close_mode_2_z_uniform.pdf") diff --git a/examples/non_abelian/open/abelian_open.py b/examples/non_abelian/open/abelian_open.py new file mode 100644 index 0000000..cc02607 --- /dev/null +++ b/examples/non_abelian/open/abelian_open.py @@ -0,0 +1,68 @@ +import numpy as np +import matplotlib.pyplot as plt +import networkx as nx +from pathlib import Path + +from netsalt.quantum_graph import oversample_graph, create_quantum_graph +from netsalt.plotting import plot_single_mode +from netsalt.physics import dispersion_relation_linear, set_dispersion_relation +from netsalt.modes import scan_frequencies, find_modes +from netsalt.plotting import plot_scan +from netsalt.io import save_modes, load_modes, save_qualities, load_qualities + +from make_graph import make_graph + +if __name__ == "__main__": + graph, pos = make_graph(with_leads=True) + params = { + "open_model": "open", + "n_workers": 7, + "k_n": 200, + "k_min": 5.0, + "k_max": 7.0, + "alpha_n": 100, + "alpha_min": -0.01, + "alpha_max": 0.3, + "c": len(graph.edges) * [1.0], + } + np.random.seed(42) + + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + plt.savefig("graph.pdf") + + create_quantum_graph(graph, params=params, positions=pos) + set_dispersion_relation(graph, dispersion_relation_linear) + + if not Path("modes_open.h5").exists(): + qualities = scan_frequencies(graph) + save_qualities(qualities, "qualities_open.h5") + modes_df = find_modes(graph, qualities) + save_modes(modes_df, filename="modes_open.h5") + else: + qualities = load_qualities("qualities_open.h5") + modes_df = load_modes("modes_open.h5") + print(modes_df) + + plt.figure(figsize=(5, 3)) + ax = plt.gca() + plot_scan(graph, qualities, modes_df=modes_df, ax=ax) + plt.savefig("open_scan.pdf") + + over_graph = oversample_graph(graph, 0.01) + over_graph.graph["params"]["c"] = len(over_graph.edges) * [1.0] + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig("open_mode_1.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig("open_mode_2.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 5, ax=plt.gca(), norm_type="real") + plt.tight_layout() + plt.savefig("open_mode_3.pdf") diff --git a/examples/non_abelian/open/make_graph.py b/examples/non_abelian/open/make_graph.py new file mode 100644 index 0000000..5603c23 --- /dev/null +++ b/examples/non_abelian/open/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/open/so3_open_random.py b/examples/non_abelian/open/so3_open_random.py new file mode 100644 index 0000000..ec45da0 --- /dev/null +++ b/examples/non_abelian/open/so3_open_random.py @@ -0,0 +1,130 @@ +import numpy as np +from pathlib import Path +import matplotlib.pyplot as plt +import networkx as nx + +from netsalt.quantum_graph import create_quantum_graph, oversample_graph +from netsalt.modes import find_modes, from_complex +from netsalt.modes import scan_frequencies +from netsalt.plotting import plot_scan, plot_single_mode +from netsalt.io import save_modes, load_modes, save_qualities, load_qualities + +from netsalt.non_abelian import construct_so3_laplacian + +from make_graph import make_graph + +if __name__ == "__main__": + np.random.seed(42) + graph, pos = make_graph(with_leads=True) + params = { + "open_model": "open", + "n_workers": 7, + "quality_threshold": 1e-7, + "max_steps": 200, + "k_n": 200, + "k_min": 5.0, + "k_max": 7.0, + "alpha_n": 100, + "alpha_min": -0.01, + "alpha_max": 0.3, + "c": len(graph.edges) * [1.0], + "laplacian_constructor": construct_so3_laplacian, + } + np.random.seed(42) + + 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): + graph[u][v]["chi"] = np.random.normal(0, 1, 3) + graph[u][v]["chi"] /= np.linalg.norm(graph[u][v]["chi"]) + pos = [graph.nodes[u]["position"] for u in graph] + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + + if not Path("qualities_non_uniform.h5").exists(): + qualities = scan_frequencies(graph) + save_qualities(qualities, "qualities_non_uniform.h5") + else: + qualities = load_qualities("qualities_non_uniform.h5") + + if not Path("modes_non_uniform.h5").exists(): + modes_df = find_modes(graph, qualities) + save_modes(modes_df, filename="modes_non_uniform.h5") + else: + modes_df = load_modes("modes_non_uniform.h5") + print(modes_df) + + plt.figure(figsize=(5, 3)) + ax = plt.gca() + plot_scan(graph, qualities, modes_df=modes_df, ax=ax) + plt.savefig("so3_open_scan_non_uniform.pdf") + + over_graph = oversample_graph(graph, 0.01) + # modes_df = find_modes(over_graph, qualities) + print(modes_df) + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_x_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_y_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_z_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 1, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 1, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_x_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 1, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_y_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 1, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_z_non.pdf") + + print("lkjlkj") + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 10, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 10, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_x_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 10, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_y_non.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 10, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_z_non.pdf") diff --git a/examples/non_abelian/open/so3_open_uniform.py b/examples/non_abelian/open/so3_open_uniform.py new file mode 100644 index 0000000..72e694e --- /dev/null +++ b/examples/non_abelian/open/so3_open_uniform.py @@ -0,0 +1,122 @@ +import numpy as np +from pathlib import Path +import matplotlib.pyplot as plt +import networkx as nx + +from netsalt.quantum_graph import create_quantum_graph, oversample_graph +from netsalt.modes import find_modes, from_complex +from netsalt.modes import scan_frequencies +from netsalt.plotting import plot_scan, plot_single_mode +from netsalt.io import save_modes, load_modes, save_qualities, load_qualities + +from netsalt.non_abelian import construct_so3_laplacian + +from make_graph import make_graph + +if __name__ == "__main__": + graph, pos = make_graph(with_leads=True) + params = { + "open_model": "open", + "n_workers": 7, + "quality_threshold": 1e-7, + "max_steps": 200, + "k_n": 200, + "k_min": 5.0, + "k_max": 7.0, + "alpha_n": 100, + "alpha_min": -0.01, + "alpha_max": 0.3, + "c": len(graph.edges) * [1.0], + "laplacian_constructor": construct_so3_laplacian, + } + np.random.seed(42) + + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + + create_quantum_graph(graph, params=params, positions=pos) + + pos = [graph.nodes[u]["position"] for u in graph] + nx.draw(graph, pos=pos) + nx.draw_networkx_labels(graph, pos=pos) + + if not Path("qualities_uniform.h5").exists(): + qualities = scan_frequencies(graph) + save_qualities(qualities, "qualities_uniform.h5") + else: + qualities = load_qualities("qualities_uniform.h5") + + if not Path("modes_uniform.h5").exists(): + modes_df = find_modes(graph, qualities) + save_modes(modes_df, filename="modes_uniform.h5") + else: + modes_df = load_modes("modes_uniform.h5") + # modes_df = load_modes("modes_open.h5") + print(modes_df) + + plt.figure(figsize=(5, 3)) + ax = plt.gca() + plot_scan(graph, qualities, modes_df=modes_df, ax=ax) + plt.savefig("so3_open_scan_uniform.pdf") + + over_graph = oversample_graph(graph, 0.01) + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_x.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_y.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 0, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_1_z.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 2, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 2, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_x.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 2, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_y.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 2, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_2_z.pdf") + print('lkjlkj') + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 5, norm_type="real", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 5, norm_type="real_x", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_x.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 5, norm_type="real_y", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_y.pdf") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 5, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_z.pdf") 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 bdcca35..7eb8bcf 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -9,20 +9,20 @@ import scipy as sc from tqdm import tqdm -from .algorithm import ( +from netsalt.algorithm import ( clean_duplicate_modes, find_rough_modes_from_scan, refine_mode_brownian_ratchet, ) -from .physics import gamma, q_value -from .quantum_graph import ( +from netsalt.physics import gamma, q_value +from netsalt.quantum_graph import ( construct_incidence_matrix, construct_laplacian, construct_weight_matrix, - set_wavenumber, mode_quality, ) -from .utils import from_complex, get_scan_grid, to_complex +from netsalt.utils import from_complex, get_scan_grid, to_complex +from netsalt.non_abelian import norm, Ad, proj_paral, proj_perp warnings.filterwarnings("ignore") warnings.filterwarnings("error", category=np.ComplexWarning) @@ -187,8 +187,7 @@ def _graph_norm(BT, Bout, Winv, z_matrix, node_solution, mask): """Compute the norm of the node solution on the graph.""" weight_matrix = Winv.dot(z_matrix).dot(Winv) inner_matrix = BT.dot(weight_matrix).dot(mask).dot(Bout) - norm = node_solution.T.dot(inner_matrix.dot(node_solution)) - return norm + return node_solution.T.dot(inner_matrix.dot(node_solution)) def compute_z_matrix(graph): @@ -264,11 +263,13 @@ def pump_linear(mode_0, graph, D0_0, D0_1): def mode_on_nodes(mode, graph): """Compute the mode solution on the nodes of the graph.""" - laplacian = construct_laplacian(to_complex(mode), graph) + laplacian = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + to_complex(mode), graph + )[0] min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=1, sigma=0, v0=np.ones(len(graph)), which="LM" + laplacian, k=5, sigma=0, v0=np.ones(laplacian.shape[0]), which="LM" ) - quality_thresh = graph.graph["params"].get("quality_threshold", 1e-4) + quality_thresh = graph.graph["params"].get("quality_threshold", 1e-2) if abs(min_eigenvalue[0]) > quality_thresh: raise Exception( "Not a mode, as quality is too high: " @@ -278,44 +279,119 @@ def mode_on_nodes(mode, graph): + ", mode: " + str(mode) ) + # this 1.5 is fairly arbitrary, it is so get similar igenvalues, close to 0, regardless of the + # choice of quality_threshold, which may pick up other small ones, not 0 if set to large values + mask = np.logical_or( + np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue)), + np.abs(min_eigenvalue) < 1e-4, + ) - return node_solution[:, 0] + n_eig = len([1 for a in mask if a]) + print(min_eigenvalue, np.abs(min_eigenvalue)) + if n_eig > 1: + L.info(f"We found {n_eig} vanishing eigenvalues, we will use the sum of their eigenvectors") + print(f"We found {n_eig} vanishing eigenvalues, we will use the sum of their eigenvectors") + return node_solution[:, mask].sum(axis=1) def flux_on_edges(mode, graph): """Compute the flux on each edge (in both directions).""" - node_solution = mode_on_nodes(mode, graph) - - _, B = construct_incidence_matrix(graph) - Winv = construct_weight_matrix(graph, with_k=False) - + _, _, B, Winv = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + to_complex(mode), graph, with_k=False + ) return Winv.dot(B).dot(node_solution) -def mean_mode_on_edges(mode, graph): - r"""Compute the average :math:`Real(E^2)` on each edge.""" +def mean_mode_on_edges(mode, graph, norm_type="abs"): + r"""Compute the average edge value.""" edge_flux = flux_on_edges(mode, graph) + return mean_on_edges(edge_flux, graph, norm_type=norm_type, mode=mode) + +def mean_on_edges(edge_flux, graph, norm_type="abs", mode=None): + r"""Mean edge values from any edge vector. + + With norm_type='abs', we use \int |u|^2 dx (does not work for non-abelian) + With norm_type='real', we use Real \int u^2 dx, needs mode value + """ mean_edge_solution = np.zeros(len(graph.edges)) - for ei in range(len(graph.edges)): - k = 1.0j * graph.graph["ks"][ei] - length = graph.graph["lengths"][ei] - z = np.zeros([2, 2], dtype=np.complex128) - if abs(np.real(k)) > 0: # in case we deal with closed graph, we have 0 / 0 - 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] + if norm_type == "abs": + for ei in range(len(graph.edges)): + k = 1.0j * graph.graph["ks"][ei] + length = graph.graph["lengths"][ei] + z = np.zeros([2, 2], dtype=np.complex128) - mean_edge_solution[ei] = np.abs( - edge_flux[2 * ei : 2 * ei + 2].T.dot(z.dot(np.conj(edge_flux[2 * ei : 2 * ei + 2]))) - ) + if abs(np.real(k)) > 0: # in case we deal with closed graph, we have 0 / 0 + z[0, 0] = (np.exp(length * (k + np.conj(k))) - 1.0) / (length * (k + np.conj(k))) + else: + z[0, 0] = 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.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'") + if len(edge_flux) > 2 * len(graph.edges): + DIM = 3 + else: + DIM = 1 + + def _ext(i, shift=1): + return slice(DIM * i, DIM * (i + shift)) + + for ei in range(len(graph.edges)): + chi = graph.graph["ks"][ei] + length = graph.graph["lengths"][ei] + z = np.zeros([2 * DIM, 2 * DIM], dtype=np.complex128) + if DIM == 3: + axis = None + if norm_type.endswith("x"): + 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] + + import scipy.integrate as integrate + + res = integrate.quad(sol, 0, length)[0] + mean_edge_solution[ei] = res + + if DIM == 1: + + 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) + + mean_edge_solution[ei] = np.real( + edge_flux[_ext(2 * ei, shift=2)].T.dot(z).dot(edge_flux[_ext(2 * ei, shift=2)]) + ) return mean_edge_solution @@ -898,19 +974,48 @@ def lasing_threshold_linear(mode, graph, D0): def get_node_transfer(k, graph, input_flow): """Compute node transfer from a given input flow.""" - return sc.sparse.linalg.spsolve(construct_laplacian(k, graph), graph.graph["ks"] * input_flow) + laplacian, BT, _, _ = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + k, graph + ) + bt = np.clip(np.ceil(np.real(BT.toarray())), -1, 0) + ks = np.empty(2 * len(graph.graph["ks"])) + ks[::2] = graph.graph["ks"] + ks[1::2] = graph.graph["ks"] + K = bt.dot(ks) + return sc.sparse.linalg.spsolve(laplacian, K * input_flow) def get_edge_transfer(k, graph, input_flow): """Compute edge transfer from a given input flow.""" - set_wavenumber(graph, k) - BT, B = construct_incidence_matrix(graph) - _r = get_node_transfer(k, graph, BT.dot(input_flow)) - Winv = construct_weight_matrix(graph, with_k=False) + laplacian, BT, B, _ = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + k, graph, with_k=True + ) + + _, _, _, Winv = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + k, graph, with_k=False + ) + + shape_k = list(np.shape(graph.graph["ks"])) + + if len(shape_k) > 1: + ks = np.zeros((6 * shape_k[0], 6 * shape_k[0]), dtype=np.complex) + for ei in range(shape_k[0]): + k = graph.graph["ks"][ei] + ks[2 * 3 * ei : 3 * (2 * ei + 1), 2 * 3 * ei : 3 * (2 * ei + 1)] = k + ks[3 * (2 * ei + 1) : 3 * (2 * ei + 2), 3 * (2 * ei + 1) : 3 * (2 * ei + 2)] = k + else: + shape_k[0] *= 2 + ks = np.empty(shape_k, dtype=np.complex) + ks[::2] = graph.graph["ks"] + ks[1::2] = graph.graph["ks"] + ks = np.diag(ks) + + _in = BT.dot(ks.dot(input_flow)) + _r = sc.sparse.linalg.spsolve(laplacian, _in) return Winv.dot(B).dot(_r) -def estimate_boundary_flow(graph, input_flow, k_frac=1e-2): +def estimate_boundary_flow(graph, input_flow, k_frac=1e-3): """Estimate boundary flow for static simulations. Arguments: @@ -926,7 +1031,7 @@ def estimate_boundary_flow(graph, input_flow, k_frac=1e-2): ) e_deg = np.array([len(graph[v]) for u, v in graph.edges]) - output_ids = list(np.argwhere(e_deg == 1).flatten()) + output_ids = list(2 * np.argwhere(e_deg == 1).flatten()) output_ids += list(2 * np.argwhere(e_deg == 1).flatten() + 1) # get the flows on all nodes diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py new file mode 100644 index 0000000..b6745bd --- /dev/null +++ b/netsalt/non_abelian.py @@ -0,0 +1,146 @@ +"""Module for non-abelian quantum graphs.""" +import numpy as np +from scipy import sparse, linalg + +DIM = 3 + + +def hat_inv(xi_vec): + """Convert vector Lie algebra to matrix Lie algebra element.""" + xi = np.zeros((3, 3), dtype=np.complex128) + xi[1, 2] = -xi_vec[0] + xi[2, 1] = xi_vec[0] + xi[0, 2] = xi_vec[1] + xi[2, 0] = -xi_vec[1] + xi[0, 1] = -xi_vec[2] + xi[1, 0] = xi_vec[2] + return xi + + +def hat(xi): + """Convert matrix Lie algebra to vector Lie algebra element.""" + xi_vec = np.zeros(3, dtype=np.complex128) + xi_vec[0] = xi[2, 1] + xi_vec[1] = xi[0, 2] + xi_vec[2] = xi[1, 0] + return xi_vec + + +def proj_perp(chi_mat): + """Perpendicular projection.""" + return chi_mat.dot(chi_mat.T) / norm(chi_mat) ** 2 + + +def proj_paral(chi_mat): + """Paralell projection.""" + return np.eye(3) - proj_perp(chi_mat) + + +def norm(chi_mat): + """Norm of chi""" + return np.sqrt(np.trace(0.5 * chi_mat.dot(chi_mat.T))) + + +def Ad(chi_mat): + """Adjoint action.""" + return linalg.expm(chi_mat) + + +def set_so3_wavenumber(graph, wavenumber): + """Set so3 matrix wavenumber.""" + chis = np.array([graph[u][v].get("chi", None) for u, v in graph.edges]) + if chis[0] is None: + chi = hat_inv(np.array([0.0, 0.0, 1.0])) + chis = np.array(len(graph.edges) * [chi]) + else: + if len(np.shape(chis[0])) == 1: + chis = np.array([hat_inv(chi) for chi in chis]) + chi_loss = np.array([proj_perp(chi) for chi in chis]) + graph.graph["ks"] = chis * np.real(wavenumber) - chi_loss * np.imag(wavenumber) + graph.graph["chis"] = chis + graph.graph["wavenumber"] = wavenumber + + +def construct_so3_incidence_matrix(graph, abelian_scale=1.0): + """Construct SO3 incidence matrix.""" + + def _ext(i): + return slice(DIM * i, DIM * (i + 1)) + + B = sparse.lil_matrix((len(graph.edges) * 2 * DIM, len(graph) * DIM), dtype=np.complex128) + BT = sparse.lil_matrix((len(graph) * DIM, len(graph.edges) * 2 * DIM), dtype=np.complex128) + for ei, (u, v) in enumerate(graph.edges): + + k = graph.graph["ks"][ei] + chi = graph.graph["chis"][ei] + length = graph.graph["lengths"][ei] + + one = np.eye(DIM) + expl = Ad(length * k) + expl = np.array(expl.dot(proj_perp(chi)), dtype=np.complex128) + expl += abelian_scale * np.exp(1.0j * length * graph.graph["wavenumber"]) * proj_paral(chi) + + B[_ext(2 * ei), _ext(u)] = -one + B[_ext(2 * ei), _ext(v)] = expl + B[_ext(2 * ei + 1), _ext(u)] = expl + B[_ext(2 * ei + 1), _ext(v)] = -one + + BT[_ext(u), _ext(2 * ei)] = -one + BT[_ext(v), _ext(2 * ei)] = expl + BT[_ext(u), _ext(2 * ei + 1)] = expl + BT[_ext(v), _ext(2 * ei + 1)] = -one + + if graph.graph["params"]["open_model"] == "open": + if len(graph[u]) == 1 or len(graph[v]) == 1: + BT[_ext(v), _ext(2 * ei)] = 0 + BT[_ext(u), _ext(2 * ei + 1)] = 0 + + if graph.graph["params"]["open_model"] == "directed": + BT[_ext(u), _ext(2 * ei + 1)] = 0 + BT[_ext(v), _ext(2 * ei + 1)] = 0 + + if graph.graph["params"]["open_model"] == "directed_reversed": + B[_ext(2 * ei + 1), _ext(u)] = 0 + B[_ext(2 * ei + 1), _ext(v)] = 0 + + return BT, B + + +def construct_so3_weight_matrix(graph, with_k=True, abelian_scale=1.0): + """Construct SO3 weight matrix.""" + + def _ext(i): + return slice(DIM * i, DIM * (i + 1)) + + Winv = sparse.lil_matrix( + (len(graph.edges) * 2 * DIM, len(graph.edges) * 2 * DIM), dtype=np.complex128 + ) + for ei, _ in enumerate(graph.edges): + k = graph.graph["ks"][ei] + chi = graph.graph["chis"][ei] + length = graph.graph["lengths"][ei] + + w_perp = Ad(2.0 * length * k).dot(proj_perp(chi)) + w_paral = ( + abelian_scale * np.exp(2.0j * length * graph.graph["wavenumber"]) * proj_paral(chi) + ) + w = w_perp + w_paral - np.eye(3) + + winv = linalg.inv(w) + + if with_k: + winv = (k.dot(proj_perp(chi)) + 1.0j * graph.graph["wavenumber"] * proj_paral(chi)).dot( + winv + ) + + Winv[_ext(2 * ei), _ext(2 * ei)] = winv + Winv[_ext(2 * ei + 1), _ext(2 * ei + 1)] = winv + return Winv + + +def construct_so3_laplacian(wavenumber, graph, abelian_scale=1.0, with_k=True): + """Construct quantum laplacian from a graph.""" + set_so3_wavenumber(graph, wavenumber) + BT, B = construct_so3_incidence_matrix(graph, abelian_scale=abelian_scale) + Winv = construct_so3_weight_matrix(graph, abelian_scale=abelian_scale, with_k=with_k) + return BT.dot(Winv).dot(B), BT, B, Winv diff --git a/netsalt/physics.py b/netsalt/physics.py index 47e536a..6bbd4ec 100644 --- a/netsalt/physics.py +++ b/netsalt/physics.py @@ -62,7 +62,7 @@ def dispersion_relation_resistance(freq, params=None): .. math:: - k(\omega) = \sqrt{\frac{\omega^2}{c^2} - i R C \omega} + k(\omega) = \sqrt{\frac{\omega^2}{c^2} + i R C \omega} Args: freq (float): frequency @@ -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): diff --git a/netsalt/plotting.py b/netsalt/plotting.py index 2e5cda0..722eafb 100644 --- a/netsalt/plotting.py +++ b/netsalt/plotting.py @@ -422,13 +422,21 @@ def plot_single_mode( edge_vmin=None, edge_vmax=None, cmap="coolwarm", + norm_type="abs", ): """Plot single mode on the graph.""" mode = modes_df[df_entry][index] if df_entry == "threshold_lasing_modes": graph.graph["params"]["D0"] = modes_df["lasing_thresholds"][index] ax = _plot_single_mode( - graph, mode, ax=ax, colorbar=colorbar, edge_vmin=edge_vmin, edge_vmax=edge_vmax, cmap=cmap + graph, + mode, + ax=ax, + colorbar=colorbar, + edge_vmin=edge_vmin, + edge_vmax=edge_vmax, + cmap=cmap, + norm_type=norm_type, ) ax.set_title( "mode " @@ -439,26 +447,60 @@ def plot_single_mode( def _plot_single_mode( - graph, mode, ax=None, colorbar=True, edge_vmin=None, edge_vmax=None, cmap="coolwarm" + graph, + mode, + ax=None, + colorbar=True, + edge_vmin=None, + edge_vmax=None, + cmap="coolwarm", + norm_type="abs", ): - positions = [graph.nodes[u]["position"] for u in graph] - edge_solution = mean_mode_on_edges(mode, graph) + edge_solution = mean_mode_on_edges(mode, graph, norm_type=norm_type) + return plot_on_graph( + graph, + edge_solution, + edge_vmin=edge_vmin, + edge_vmax=edge_vmax, + ax=ax, + cmap=cmap, + colorbar=colorbar, + norm_type=norm_type, + ) + + +def plot_on_graph( + graph, + edge_data, + edge_vmin=None, + edge_vmax=None, + ax=None, + cmap="coolwarm", + colorbar=True, + norm_type="abs", +): + """Plot edge data on graph.""" if ax is None: plt.figure(figsize=(5, 4)) # 14,3 ax = plt.gca() + positions = [graph.nodes[u]["position"] for u in graph] nx.draw(graph, pos=positions, node_size=0, width=0, ax=ax) cmap = plt.get_cmap(cmap) if edge_vmax is None: - edge_vmax = max(abs(edge_solution)) + edge_vmax = max(abs(edge_data)) if edge_vmin is None: - edge_vmin = -max(abs(edge_solution)) + if norm_type == "abs": + edge_vmin = 0 + if norm_type.startswith("real"): + edge_vmin = -max(abs(edge_data)) + nx.draw_networkx_edges( graph, pos=positions, - edge_color=edge_solution, + edge_color=edge_data, width=2, edge_cmap=cmap, edge_vmin=edge_vmin, @@ -468,7 +510,10 @@ def _plot_single_mode( if colorbar: sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=edge_vmin, vmax=edge_vmax)) sm.set_array([]) - plt.colorbar(sm, ax=ax, label=r"$|E|^2$ (a.u)", shrink=0.5) + if norm_type == "abs": + plt.colorbar(sm, ax=ax, label=r"$|E|^2$ (a.u)", shrink=0.5) + if norm_type.startswith("real"): + plt.colorbar(sm, ax=ax, label=r"Real$(E^2)$ (a.u)", shrink=0.5) return ax diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index f59aca4..a7b14a7 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -211,6 +211,9 @@ def oversample_graph(graph, edge_size): # pylint: disable=too-many-locals n_nodes = int(graph[u][v]["length"] / edge_size) if n_nodes > 1: dielectric_constant = graph[u][v].get("dielectric_constant", None) + chi = None + if "chi" in graph[u][v]: + chi = graph[u][v]["chi"] pump = graph[u][v]["pump"] oversampled_graph.remove_edge(u, v) @@ -233,6 +236,7 @@ def oversample_graph(graph, edge_size): # pylint: disable=too-many-locals first, last, dielectric_constant=dielectric_constant, + chi=chi, pump=pump, edgelabel=ei, ) @@ -241,6 +245,7 @@ def oversample_graph(graph, edge_size): # pylint: disable=too-many-locals last_node + node_index, v, dielectric_constant=dielectric_constant, + chi=chi, pump=pump, edgelabel=ei, ) @@ -255,7 +260,7 @@ def oversample_graph(graph, edge_size): # pylint: disable=too-many-locals return oversampled_graph -def construct_laplacian(wavenumber, graph): +def construct_laplacian(wavenumber, graph, with_k=True): """Construct quantum laplacian from a graph. The quantum laplacian is L(k) = B^T(k) W^{-1}(k) B(k), with quantum incidence and weight matrix. @@ -266,7 +271,7 @@ def construct_laplacian(wavenumber, graph): """ set_wavenumber(graph, wavenumber) BT, B = construct_incidence_matrix(graph) - Winv = construct_weight_matrix(graph) + Winv = construct_weight_matrix(graph, with_k=with_k) laplacian = BT.dot(Winv).dot(B) node_loss = graph.graph["params"].get("node_loss", 0) @@ -275,7 +280,7 @@ def construct_laplacian(wavenumber, graph): [graph[u].get("node_loss", node_loss) for u in graph.nodes()] ) - return laplacian + return BT.dot(Winv).dot(B), BT, B, Winv def set_wavenumber(graph, wavenumber): @@ -326,7 +331,7 @@ def construct_incidence_matrix(graph): def construct_weight_matrix(graph, with_k=True): """Construct the quantum matrix W^{-1}(k). - The with_k argument is needed for the graph laplcian, not for computing the edge amplitudes. + The with_k argument is needed for the graph laplacian, not for computing the edge amplitudes. Args: graph (graph): quantum graph @@ -455,5 +460,7 @@ def mode_quality(mode, graph, quality_method="eigenvalue"): graph (graph): quantum graph quality_method (str): method for quality evaluation (eig, singular value or det) """ - laplacian = construct_laplacian(to_complex(mode), graph) + laplacian = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + to_complex(mode), graph + )[0] return laplacian_quality(laplacian, method=quality_method)