From 1c0a139c301784d33926b9ccd09db3e5e0fedd50 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 21 Nov 2022 16:38:28 +0100 Subject: [PATCH 01/22] tmp test --- examples/non_abelian/test.py | 238 +++++++++++++++++++++++++++++++++++ 1 file changed, 238 insertions(+) create mode 100644 examples/non_abelian/test.py diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py new file mode 100644 index 0000000..197b972 --- /dev/null +++ b/examples/non_abelian/test.py @@ -0,0 +1,238 @@ +import numpy as np +from tqdm import tqdm +import matplotlib.pyplot as plt +import scipy as sc +from netsalt.quantum_graph import create_quantum_graph, laplacian_quality, mode_quality +from netsalt.modes import mode_on_nodes + +from netsalt.physics import dispersion_relation_linear, set_dispersion_relation +import networkx as nx +from scipy import sparse, linalg + + +def hat_inv(xi_vec): + """Convert vector Lie algebra to matrix Lie algebra element.""" + xi = np.zeros((3, 3)) + 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) + 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): + return np.eye(3) - proj_paral(chi) + + +def proj_paral(chi_mat, center_scale=1.0): + chi_vec = hat(chi_mat) + return center_scale * np.outer(chi_vec, chi_vec) / np.linalg.norm(chi_vec) ** 2 + + +def norm(chi_mat): + chi_vec = hat(chi_mat) + return np.linalg.norm(chi_vec) + + +def Ad(chi_mat): + return linalg.expm(chi_mat) + + +def set_so3_wavenumber(graph, wavenumber): + x = np.array([0.2, 0.2, 1.0]) + chi_vec = wavenumber * x / np.linalg.norm(x) + chi_mat = hat_inv(chi_vec) + graph.graph["ks"] = len(graph.edges) * [chi_mat] + + x2 = np.array([1.0, 0.2, 0.2]) + chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) + chi_mat2 = hat_inv(chi_vec2) + graph.graph["ks"][:10] = 10 * [chi_mat2] + graph.graph["ks"][20:30] = 10 * [chi_mat2] + + +def construct_so3_incidence_matrix(graph): + dim = 3 + + def _ext(i): + return slice(dim * i, dim * (i + 1)) + + Bout = 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): + + one = np.eye(dim) + expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) + expl = 0.0j + expl.dot(proj_perp(graph.graph["ks"][ei])) + expl += proj_paral(graph.graph["ks"][ei]) * np.exp( + 1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei]) + ) + out = True if (len(graph[u]) == 1 or len(graph[v]) == 1) else False + + Bout[_ext(2 * ei), _ext(u)] = -one + Bout[_ext(2 * ei), _ext(v)] = 0 if out else expl + Bout[_ext(2 * ei + 1), _ext(u)] = 0 if out else expl + Bout[_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 + + return BT, Bout + + +def construct_so3_weight_matrix(graph, with_k=True): + dim = 3 + + 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, (u, v) in enumerate(graph.edges): + chi = graph.graph["ks"][ei] + length = graph.graph["lengths"][ei] + + w_perp = Ad(2.0 * length * chi).dot(proj_perp(chi)) + w_paral = np.exp(2.0j * length * norm(chi)) * proj_paral(chi) + + w = w_perp + w_paral - np.eye(3) + + winv = linalg.inv(w) + + if with_k: + winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * 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): + """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. + + Args: + wavenumber (complex): wavenumber + graph (graph): quantum graph + """ + set_so3_wavenumber(graph, wavenumber) + BT, B = construct_so3_incidence_matrix(graph) + Winv = construct_so3_weight_matrix(graph) + return BT.dot(Winv).dot(B) + + +def so3_mode_on_nodes(laplacian): + v0 = np.random.random(laplacian.shape[0]) + min_eigenvalue, node_solution = sc.sparse.linalg.eigs( + laplacian, k=1, sigma=0, v0=v0, which="LM" + ) + quality_thresh = graph.graph["params"].get("quality_threshold", 1e2) + if abs(min_eigenvalue[0]) > quality_thresh: + raise Exception( + "Not a mode, as quality is too high: " + + str(abs(min_eigenvalue[0])) + + " > " + + str(quality_thresh) + ) + + return node_solution[:, 0] + + +if __name__ == "__main__": + + params = {"open_model": "open", "c": 1.0} + n = 70 + + graph = nx.cycle_graph(n) + # graph.add_edge(0, n) + # graph.add_edge(10, n + 1) + # graph.add_edge(20, n + 2) + + graph.add_edge(0, 2) + graph.add_edge(0, 5) + graph.add_edge(0, 7) + graph.add_edge(0, 20) + graph.add_edge(0, 50) + graph.add_edge(10, 12) + graph.add_edge(20, 12) + + graph_u1 = nx.cycle_graph(n) + # graph_u1.add_edge(0, n) + # graph_u1.add_edge(10, n + 1) + # graph_u1.add_edge(20, n + 2) + + graph_u1.add_edge(0, 2) + graph_u1.add_edge(0, 5) + graph_u1.add_edge(0, 7) + graph_u1.add_edge(0, 20) + graph_u1.add_edge(0, 50) + graph_u1.add_edge(10, 12) + graph_u1.add_edge(20, 12) + + x = np.linspace(0, 2 * np.pi / (len(graph) - 1), len(graph)) + pos = 100.1 * np.array([np.cos(x), np.sin(x)]).T + pos = list(pos) + pos.append([0, 0]) + create_quantum_graph(graph, params=params, positions=pos) + create_quantum_graph(graph_u1, params=params, positions=pos) + set_dispersion_relation(graph_u1, dispersion_relation_linear) + + ks = np.linspace(20.0, 25.0, 500) + qs = [] + qs_u1 = [] + for k in tqdm(ks): + L = construct_so3_laplacian(k, graph) + qs.append(laplacian_quality(L)) + qs_u1.append(mode_quality([k, 0], graph_u1)) + + plt.figure() + plt.plot(ks, qs_u1, "+-r") + plt.plot(ks, qs, "-") + plt.yscale("log") + plt.show() + + k = ks[np.argmin(qs)] + L = construct_so3_laplacian(k, graph) + mode = so3_mode_on_nodes(L) + + k_u1 = k # ks[np.argmin(qs_u1)] + mode_u1 = mode_on_nodes([k_u1, 0], graph_u1) + + plt.figure() + x = np.abs(mode[::3]) + y = np.abs(mode[1::3]) + z = np.abs(mode[2::3]) + plt.plot(x, label="x") + plt.plot(y, label="y") + plt.plot(z, label="z") + n = np.sqrt(x**2 + y**2 + z**2) + + plt.plot(n, "+-", label="norm") + plt.plot(np.abs(mode_u1), label="u1") + plt.legend() + + plt.figure() + plt.plot(ks, qs_u1, "+-r") + plt.plot(ks, qs, "-") + plt.axvline(k) + plt.axvline(k_u1, c="r") + plt.yscale("log") + + plt.figure() + plt.plot(np.sqrt(x**2 + y**2 + z**2) - np.abs(mode_u1)) + plt.show() From 8b2a8e03b05bbffbffc15d670848d2f00a885b92 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 22 Nov 2022 10:07:08 +0100 Subject: [PATCH 02/22] more --- examples/non_abelian/test.py | 12 ++++++------ netsalt/modes.py | 2 +- netsalt/quantum_graph.py | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py index 197b972..c169242 100644 --- a/examples/non_abelian/test.py +++ b/examples/non_abelian/test.py @@ -159,9 +159,9 @@ def so3_mode_on_nodes(laplacian): n = 70 graph = nx.cycle_graph(n) - # graph.add_edge(0, n) - # graph.add_edge(10, n + 1) - # graph.add_edge(20, n + 2) + graph.add_edge(0, n) + graph.add_edge(10, n + 1) + graph.add_edge(20, n + 2) graph.add_edge(0, 2) graph.add_edge(0, 5) @@ -172,9 +172,9 @@ def so3_mode_on_nodes(laplacian): graph.add_edge(20, 12) graph_u1 = nx.cycle_graph(n) - # graph_u1.add_edge(0, n) - # graph_u1.add_edge(10, n + 1) - # graph_u1.add_edge(20, n + 2) + graph_u1.add_edge(0, n) + graph_u1.add_edge(10, n + 1) + graph_u1.add_edge(20, n + 2) graph_u1.add_edge(0, 2) graph_u1.add_edge(0, 5) diff --git a/netsalt/modes.py b/netsalt/modes.py index bdcca35..968fff9 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -268,7 +268,7 @@ def mode_on_nodes(mode, graph): min_eigenvalue, node_solution = sc.sparse.linalg.eigs( laplacian, k=1, sigma=0, v0=np.ones(len(graph)), 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: " diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index f59aca4..f1366c3 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -326,7 +326,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 From ace90cc915c7ab37129f9222509a3b94640022fa Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 24 Nov 2022 14:59:58 +0100 Subject: [PATCH 03/22] more --- examples/non_abelian/test.py | 191 ++++------------------------------- netsalt/non_abelian.py | 145 ++++++++++++++++++++++++++ 2 files changed, 162 insertions(+), 174 deletions(-) create mode 100644 netsalt/non_abelian.py diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py index c169242..361839c 100644 --- a/examples/non_abelian/test.py +++ b/examples/non_abelian/test.py @@ -1,202 +1,45 @@ import numpy as np from tqdm import tqdm import matplotlib.pyplot as plt -import scipy as sc from netsalt.quantum_graph import create_quantum_graph, laplacian_quality, mode_quality from netsalt.modes import mode_on_nodes from netsalt.physics import dispersion_relation_linear, set_dispersion_relation import networkx as nx -from scipy import sparse, linalg -def hat_inv(xi_vec): - """Convert vector Lie algebra to matrix Lie algebra element.""" - xi = np.zeros((3, 3)) - 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) - 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): - return np.eye(3) - proj_paral(chi) - - -def proj_paral(chi_mat, center_scale=1.0): - chi_vec = hat(chi_mat) - return center_scale * np.outer(chi_vec, chi_vec) / np.linalg.norm(chi_vec) ** 2 - - -def norm(chi_mat): - chi_vec = hat(chi_mat) - return np.linalg.norm(chi_vec) - - -def Ad(chi_mat): - return linalg.expm(chi_mat) - - -def set_so3_wavenumber(graph, wavenumber): - x = np.array([0.2, 0.2, 1.0]) - chi_vec = wavenumber * x / np.linalg.norm(x) - chi_mat = hat_inv(chi_vec) - graph.graph["ks"] = len(graph.edges) * [chi_mat] - - x2 = np.array([1.0, 0.2, 0.2]) - chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) - chi_mat2 = hat_inv(chi_vec2) - graph.graph["ks"][:10] = 10 * [chi_mat2] - graph.graph["ks"][20:30] = 10 * [chi_mat2] - - -def construct_so3_incidence_matrix(graph): - dim = 3 - - def _ext(i): - return slice(dim * i, dim * (i + 1)) - - Bout = 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): - - one = np.eye(dim) - expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) - expl = 0.0j + expl.dot(proj_perp(graph.graph["ks"][ei])) - expl += proj_paral(graph.graph["ks"][ei]) * np.exp( - 1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei]) - ) - out = True if (len(graph[u]) == 1 or len(graph[v]) == 1) else False - - Bout[_ext(2 * ei), _ext(u)] = -one - Bout[_ext(2 * ei), _ext(v)] = 0 if out else expl - Bout[_ext(2 * ei + 1), _ext(u)] = 0 if out else expl - Bout[_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 - - return BT, Bout - - -def construct_so3_weight_matrix(graph, with_k=True): - dim = 3 - - 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, (u, v) in enumerate(graph.edges): - chi = graph.graph["ks"][ei] - length = graph.graph["lengths"][ei] - - w_perp = Ad(2.0 * length * chi).dot(proj_perp(chi)) - w_paral = np.exp(2.0j * length * norm(chi)) * proj_paral(chi) - - w = w_perp + w_paral - np.eye(3) - - winv = linalg.inv(w) - - if with_k: - winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * 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): - """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. - - Args: - wavenumber (complex): wavenumber - graph (graph): quantum graph - """ - set_so3_wavenumber(graph, wavenumber) - BT, B = construct_so3_incidence_matrix(graph) - Winv = construct_so3_weight_matrix(graph) - return BT.dot(Winv).dot(B) - - -def so3_mode_on_nodes(laplacian): - v0 = np.random.random(laplacian.shape[0]) - min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=1, sigma=0, v0=v0, which="LM" - ) - quality_thresh = graph.graph["params"].get("quality_threshold", 1e2) - if abs(min_eigenvalue[0]) > quality_thresh: - raise Exception( - "Not a mode, as quality is too high: " - + str(abs(min_eigenvalue[0])) - + " > " - + str(quality_thresh) - ) +from netsalt.non_abelian import construct_so3_laplacian, so3_mode_on_nodes +def make_graph(): + graph = nx.cycle_graph(n) + graph.add_edge(0, 8) + graph.add_edge(0, 20) + graph.add_edge(10, 15) + x = np.linspace(0, 2*np.pi*(1-1. / (len(graph) - 1)), len(graph)) + pos = np.array([np.cos(x), np.sin(x)]).T + pos = list(pos) - return node_solution[:, 0] + return graph, pos if __name__ == "__main__": params = {"open_model": "open", "c": 1.0} - n = 70 + n = 30 + graph, pos = make_graph() + graph_u1, pos = make_graph() - graph = nx.cycle_graph(n) - graph.add_edge(0, n) - graph.add_edge(10, n + 1) - graph.add_edge(20, n + 2) + nx.draw(graph, pos=pos) - graph.add_edge(0, 2) - graph.add_edge(0, 5) - graph.add_edge(0, 7) - graph.add_edge(0, 20) - graph.add_edge(0, 50) - graph.add_edge(10, 12) - graph.add_edge(20, 12) - - graph_u1 = nx.cycle_graph(n) - graph_u1.add_edge(0, n) - graph_u1.add_edge(10, n + 1) - graph_u1.add_edge(20, n + 2) - - graph_u1.add_edge(0, 2) - graph_u1.add_edge(0, 5) - graph_u1.add_edge(0, 7) - graph_u1.add_edge(0, 20) - graph_u1.add_edge(0, 50) - graph_u1.add_edge(10, 12) - graph_u1.add_edge(20, 12) - - x = np.linspace(0, 2 * np.pi / (len(graph) - 1), len(graph)) - pos = 100.1 * np.array([np.cos(x), np.sin(x)]).T - pos = list(pos) - pos.append([0, 0]) create_quantum_graph(graph, params=params, positions=pos) create_quantum_graph(graph_u1, params=params, positions=pos) set_dispersion_relation(graph_u1, dispersion_relation_linear) - ks = np.linspace(20.0, 25.0, 500) + ks = np.linspace(1.53, 1.58, 2000) qs = [] qs_u1 = [] for k in tqdm(ks): - L = construct_so3_laplacian(k, graph) + k += 0.0966j + L = construct_so3_laplacian(k, graph, abelian_scale=20.0) qs.append(laplacian_quality(L)) qs_u1.append(mode_quality([k, 0], graph_u1)) diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py new file mode 100644 index 0000000..6830d5d --- /dev/null +++ b/netsalt/non_abelian.py @@ -0,0 +1,145 @@ +"""Module for non-abelian quantum graphs.""" +import numpy as np +import scipy as sc + +from scipy import sparse, linalg + + +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): + return np.eye(3) - proj_paral(chi) + + +def proj_paral(chi_mat): + chi_vec = hat(chi_mat) + return np.outer(chi_vec, chi_vec) / np.linalg.norm(chi_vec) ** 2 + + +def norm(chi_mat): + chi_vec = hat(chi_mat) + return np.linalg.norm(chi_vec) + + +def Ad(chi_mat): + return linalg.expm(chi_mat) + + +def set_so3_wavenumber(graph, wavenumber, chis=None): + if chis is None: + x = np.array([0.2, 0.2, 1.0]) + chi_vec = wavenumber * x / np.linalg.norm(x) + chi_mat = hat_inv(chi_vec) + graph.graph["ks"] = len(graph.edges) * [chi_mat] + + x2 = np.array([1.0, 0.2, 0.2]) + chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) + chi_mat2 = hat_inv(chi_vec2) + graph.graph["ks"][:10] = 10 * [chi_mat2] + graph.graph["ks"][20:30] = 10 * [chi_mat2] + else: + graph.graph["ks"] = chis + + +def construct_so3_incidence_matrix(graph, abelian_scale=1.0): + dim = 3 + + def _ext(i): + return slice(dim * i, dim * (i + 1)) + + Bout = 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): + + one = np.eye(dim) + expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) + 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])) + ) + out = True if (len(graph[u]) == 1 or len(graph[v]) == 1) else False + + Bout[_ext(2 * ei), _ext(u)] = -one + Bout[_ext(2 * ei), _ext(v)] = 0 if out else expl + Bout[_ext(2 * ei + 1), _ext(u)] = 0 if out else expl + Bout[_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 + + return BT, Bout + + +def construct_so3_weight_matrix(graph, with_k=True, abelian_scale=1.0): + dim = 3 + + 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, (u, v) in enumerate(graph.edges): + chi = graph.graph["ks"][ei] + length = graph.graph["lengths"][ei] + + 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) + + if with_k: + winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * 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, chis=None): + """Construct quantum laplacian from a graph.""" + set_so3_wavenumber(graph, wavenumber, chis=chis) + BT, B = construct_so3_incidence_matrix(graph, abelian_scale=abelian_scale) + Winv = construct_so3_weight_matrix(graph, abelian_scale=abelian_scale) + return BT.dot(Winv).dot(B) + + +def so3_mode_on_nodes(laplacian, quality_thresh=1e2): + v0 = np.random.random(laplacian.shape[0]) + min_eigenvalue, node_solution = sc.sparse.linalg.eigs( + laplacian, k=1, sigma=0, v0=v0, which="LM" + ) + if abs(min_eigenvalue[0]) > quality_thresh: + raise Exception( + "Not a mode, as quality is too high: " + + str(abs(min_eigenvalue[0])) + + " > " + + str(quality_thresh) + ) + + return node_solution[:, 0] From d2f49dc9f04b1d5e72c9c53c7e52b817d1af2ad2 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Sat, 15 Jul 2023 11:03:47 +0200 Subject: [PATCH 04/22] more --- examples/non_abelian/test.py | 21 +++++++++++++-------- netsalt/non_abelian.py | 10 +++++----- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py index 361839c..544d63f 100644 --- a/examples/non_abelian/test.py +++ b/examples/non_abelian/test.py @@ -9,18 +9,20 @@ from netsalt.non_abelian import construct_so3_laplacian, so3_mode_on_nodes + + def make_graph(): graph = nx.cycle_graph(n) graph.add_edge(0, 8) graph.add_edge(0, 20) graph.add_edge(10, 15) - x = np.linspace(0, 2*np.pi*(1-1. / (len(graph) - 1)), len(graph)) + 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) - return graph, pos + if __name__ == "__main__": params = {"open_model": "open", "c": 1.0} @@ -32,23 +34,26 @@ def make_graph(): create_quantum_graph(graph, params=params, positions=pos) create_quantum_graph(graph_u1, params=params, positions=pos) + set_dispersion_relation(graph_u1, dispersion_relation_linear) - ks = np.linspace(1.53, 1.58, 2000) + ks = np.linspace(1.0, 3, 200) qs = [] qs_u1 = [] for k in tqdm(ks): - k += 0.0966j - L = construct_so3_laplacian(k, graph, abelian_scale=20.0) + kim = 0 # 0.0966j + L = construct_so3_laplacian(k - kim, graph, abelian_scale=1.0) qs.append(laplacian_quality(L)) - qs_u1.append(mode_quality([k, 0], graph_u1)) + qs_u1.append(mode_quality([k, kim], graph_u1)) plt.figure() - plt.plot(ks, qs_u1, "+-r") - plt.plot(ks, qs, "-") + plt.plot(ks, qs_u1, "+-r", label="u1") + plt.plot(ks, qs, "-", label="so3") + plt.legend(loc="best") plt.yscale("log") plt.show() + k = ks[np.argmin(qs)] L = construct_so3_laplacian(k, graph) mode = so3_mode_on_nodes(L) diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 6830d5d..db4ead1 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -46,12 +46,12 @@ def Ad(chi_mat): def set_so3_wavenumber(graph, wavenumber, chis=None): if chis is None: - x = np.array([0.2, 0.2, 1.0]) + x = np.array([0.1, 0.5, 1.0]) chi_vec = wavenumber * x / np.linalg.norm(x) chi_mat = hat_inv(chi_vec) graph.graph["ks"] = len(graph.edges) * [chi_mat] - x2 = np.array([1.0, 0.2, 0.2]) + x2 = np.array([1.0, 0.5, 0.0]) chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) chi_mat2 = hat_inv(chi_vec2) graph.graph["ks"][:10] = 10 * [chi_mat2] @@ -72,7 +72,7 @@ def _ext(i): one = np.eye(dim) expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) - expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) + #expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) expl += ( abelian_scale * proj_paral(graph.graph["ks"][ei]) @@ -108,13 +108,13 @@ def _ext(i): 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) if with_k: - winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) + # winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) + winv = (chi + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) Winv[_ext(2 * ei), _ext(2 * ei)] = winv Winv[_ext(2 * ei + 1), _ext(2 * ei + 1)] = winv From a345a8d82c2031e7281d6850d139bfe579d90f09 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Sun, 16 Jul 2023 13:26:29 +0200 Subject: [PATCH 05/22] more --- examples/non_abelian/test.py | 44 ++++++++++++++++++++++++------ netsalt/non_abelian.py | 53 ++++++++++++++++++++++++++++++++---- 2 files changed, 82 insertions(+), 15 deletions(-) diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py index 544d63f..5dcb54e 100644 --- a/examples/non_abelian/test.py +++ b/examples/non_abelian/test.py @@ -4,14 +4,16 @@ from netsalt.quantum_graph import create_quantum_graph, laplacian_quality, mode_quality from netsalt.modes import mode_on_nodes +from netsalt.modes import scan_frequencies +from netsalt.plotting import plot_scan from netsalt.physics import dispersion_relation_linear, set_dispersion_relation import networkx as nx -from netsalt.non_abelian import construct_so3_laplacian, so3_mode_on_nodes +from netsalt.non_abelian import construct_so3_laplacian, so3_mode_on_nodes, scan_frequencies_so3 -def make_graph(): +def make_graph(n): graph = nx.cycle_graph(n) graph.add_edge(0, 8) graph.add_edge(0, 20) @@ -19,16 +21,30 @@ def make_graph(): 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(1, n) + graph.add_edge(15, n + 1) + pos += [[1.2, 0]] + pos += [[-1.2, 0]] return graph, pos if __name__ == "__main__": - params = {"open_model": "open", "c": 1.0} + params = { + "open_model": "open", + "c": 1.0, + "k_min": 10.0, + "k_max": 12.0, + "k_n": 200, + "alpha_min": 0.0, + "alpha_max": 0.4, + "alpha_n": 50, + "n_workers": 7, + } n = 30 - graph, pos = make_graph() - graph_u1, pos = make_graph() + graph, pos = make_graph(n) + graph_u1, pos = make_graph(n) nx.draw(graph, pos=pos) @@ -37,12 +53,23 @@ def make_graph(): set_dispersion_relation(graph_u1, dispersion_relation_linear) - ks = np.linspace(1.0, 3, 200) + qualities_u1 = scan_frequencies(graph_u1) + plot_scan(graph_u1, qualities_u1) + plt.suptitle('u1') + + qualities = scan_frequencies_so3(graph) + plot_scan(graph, qualities) + plt.suptitle("so3") + plt.show() + + +def lkj(): + ks = np.linspace(10.0, 15, 500) qs = [] qs_u1 = [] for k in tqdm(ks): - kim = 0 # 0.0966j - L = construct_so3_laplacian(k - kim, graph, abelian_scale=1.0) + kim = 0.05 + L = construct_so3_laplacian(k + 1j * kim, graph) qs.append(laplacian_quality(L)) qs_u1.append(mode_quality([k, kim], graph_u1)) @@ -53,7 +80,6 @@ def make_graph(): plt.yscale("log") plt.show() - k = ks[np.argmin(qs)] L = construct_so3_laplacian(k, graph) mode = so3_mode_on_nodes(L) diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index db4ead1..946ea6a 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -1,8 +1,12 @@ """Module for non-abelian quantum graphs.""" import numpy as np import scipy as sc +from tqdm import tqdm +import multiprocessing from scipy import sparse, linalg +from netsalt.quantum_graph import laplacian_quality +from netsalt.utils import get_scan_grid, to_complex def hat_inv(xi_vec): @@ -46,16 +50,15 @@ def Ad(chi_mat): def set_so3_wavenumber(graph, wavenumber, chis=None): if chis is None: - x = np.array([0.1, 0.5, 1.0]) + x = np.array([0.0, 0.0, 1.0]) chi_vec = wavenumber * x / np.linalg.norm(x) chi_mat = hat_inv(chi_vec) graph.graph["ks"] = len(graph.edges) * [chi_mat] - x2 = np.array([1.0, 0.5, 0.0]) + x2 = np.array([1.0, 0.0, 0.0]) chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) chi_mat2 = hat_inv(chi_vec2) graph.graph["ks"][:10] = 10 * [chi_mat2] - graph.graph["ks"][20:30] = 10 * [chi_mat2] else: graph.graph["ks"] = chis @@ -72,12 +75,13 @@ def _ext(i): one = np.eye(dim) expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) - #expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) + 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])) ) + out = True if (len(graph[u]) == 1 or len(graph[v]) == 1) else False Bout[_ext(2 * ei), _ext(u)] = -one @@ -113,8 +117,7 @@ def _ext(i): winv = linalg.inv(w) if with_k: - # winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) - winv = (chi + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) + winv = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) Winv[_ext(2 * ei), _ext(2 * ei)] = winv Winv[_ext(2 * ei + 1), _ext(2 * ei + 1)] = winv @@ -143,3 +146,41 @@ def so3_mode_on_nodes(laplacian, quality_thresh=1e2): ) return node_solution[:, 0] + + +class WorkerScan: + """Worker to scan complex frequency.""" + + def __init__(self, graph, quality_method="eigenvalue"): + self.graph = graph + self.quality_method = quality_method + np.random.seed(42) + + def __call__(self, freq): + L = construct_so3_laplacian(to_complex(freq), self.graph, abelian_scale=1.0) + return laplacian_quality(L) + + +def scan_frequencies_so3(graph, quality_method="eigenvalue"): + """Scan a range of complex frequencies and return mode qualities.""" + ks, alphas = get_scan_grid(graph) + freqs = [[k, a] for k in ks for a in alphas] + + worker_scan = WorkerScan(graph, quality_method=quality_method) + chunksize = max(1, int(0.1 * len(freqs) / graph.graph["params"]["n_workers"])) + with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool: + qualities_list = list( + tqdm( + pool.imap(worker_scan, freqs, chunksize=chunksize), + total=len(freqs), + ) + ) + + id_k = [k_i for k_i in range(len(ks)) for a_i in range(len(alphas))] + id_a = [a_i for k_i in range(len(ks)) for a_i in range(len(alphas))] + qualities = sc.sparse.coo_matrix( + (qualities_list, (id_k, id_a)), + shape=(graph.graph["params"]["k_n"], graph.graph["params"]["alpha_n"]), + ).toarray() + + return qualities From 99082a2faa6153e7118faae8a8efdafa11faac9c Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 9 Aug 2023 13:39:15 +0200 Subject: [PATCH 06/22] more --- examples/transfer/boundary_example.py | 2 +- netsalt/non_abelian.py | 110 ++++++-------------------- netsalt/physics.py | 2 +- netsalt/quantum_graph.py | 4 +- 4 files changed, 31 insertions(+), 87 deletions(-) diff --git a/examples/transfer/boundary_example.py b/examples/transfer/boundary_example.py index 70790d4..0c62107 100644 --- a/examples/transfer/boundary_example.py +++ b/examples/transfer/boundary_example.py @@ -41,7 +41,7 @@ def make_graph2(): if __name__ == "__main__": - params = {"open_model": "open", "c": 1.0, "R": 1000.0} + params = {"open_model": "open", "c": 1.0, "R": 0.0} np.random.seed(42) graph, pos = make_graph2() diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 946ea6a..6866554 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -1,12 +1,6 @@ """Module for non-abelian quantum graphs.""" import numpy as np -import scipy as sc -from tqdm import tqdm - -import multiprocessing from scipy import sparse, linalg -from netsalt.quantum_graph import laplacian_quality -from netsalt.utils import get_scan_grid, to_complex def hat_inv(xi_vec): @@ -31,49 +25,50 @@ def hat(xi): def proj_perp(chi): + """Perpendicular projection.""" return np.eye(3) - proj_paral(chi) 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 def norm(chi_mat): + """Norm of chi""" chi_vec = hat(chi_mat) return np.linalg.norm(chi_vec) def Ad(chi_mat): + """Adjoint action.""" return linalg.expm(chi_mat) -def set_so3_wavenumber(graph, wavenumber, chis=None): +def set_so3_wavenumber(graph, wavenumber): + """Set so3 matrix wavenumber.""" + chis = graph.graph["params"].get("chis", None) if chis is None: - x = np.array([0.0, 0.0, 1.0]) - chi_vec = wavenumber * x / np.linalg.norm(x) - chi_mat = hat_inv(chi_vec) - graph.graph["ks"] = len(graph.edges) * [chi_mat] - - x2 = np.array([1.0, 0.0, 0.0]) - chi_vec2 = wavenumber * x2 / np.linalg.norm(x2) - chi_mat2 = hat_inv(chi_vec2) - graph.graph["ks"][:10] = 10 * [chi_mat2] + chi = hat_inv(np.array([0.0, 0.0, 1.0])) + chis = np.array(len(graph.edges) * [chi]) else: - graph.graph["ks"] = chis + if len(np.shape(chis[0])) == 1: + chis = np.array([hat_inv(chi) for chi in chis]) + graph.graph["ks"] = chis * wavenumber def construct_so3_incidence_matrix(graph, abelian_scale=1.0): - dim = 3 + """Construct SO3 incidence matrix.""" + DIM = 3 def _ext(i): - return slice(dim * i, dim * (i + 1)) + return slice(DIM * i, DIM * (i + 1)) - Bout = 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) + Bout = 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): - - one = np.eye(dim) + one = np.eye(DIM) expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) expl += ( @@ -82,7 +77,7 @@ def _ext(i): * np.exp(1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei])) ) - out = True if (len(graph[u]) == 1 or len(graph[v]) == 1) else False + out = len(graph[u]) == 1 or len(graph[v]) == 1 Bout[_ext(2 * ei), _ext(u)] = -one Bout[_ext(2 * ei), _ext(v)] = 0 if out else expl @@ -98,15 +93,16 @@ def _ext(i): def construct_so3_weight_matrix(graph, with_k=True, abelian_scale=1.0): - dim = 3 + """Construct SO3 weight matrix.""" + DIM = 3 def _ext(i): - return slice(dim * i, dim * (i + 1)) + return slice(DIM * i, DIM * (i + 1)) Winv = sparse.lil_matrix( - (len(graph.edges) * 2 * dim, len(graph.edges) * 2 * dim), dtype=np.complex128 + (len(graph.edges) * 2 * DIM, len(graph.edges) * 2 * DIM), dtype=np.complex128 ) - for ei, (u, v) in enumerate(graph.edges): + for ei, _ in enumerate(graph.edges): chi = graph.graph["ks"][ei] length = graph.graph["lengths"][ei] @@ -124,63 +120,9 @@ def _ext(i): return Winv -def construct_so3_laplacian(wavenumber, graph, abelian_scale=1.0, chis=None): +def construct_so3_laplacian(wavenumber, graph, abelian_scale=1.0): """Construct quantum laplacian from a graph.""" - set_so3_wavenumber(graph, wavenumber, chis=chis) + 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) return BT.dot(Winv).dot(B) - - -def so3_mode_on_nodes(laplacian, quality_thresh=1e2): - v0 = np.random.random(laplacian.shape[0]) - min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=1, sigma=0, v0=v0, which="LM" - ) - if abs(min_eigenvalue[0]) > quality_thresh: - raise Exception( - "Not a mode, as quality is too high: " - + str(abs(min_eigenvalue[0])) - + " > " - + str(quality_thresh) - ) - - return node_solution[:, 0] - - -class WorkerScan: - """Worker to scan complex frequency.""" - - def __init__(self, graph, quality_method="eigenvalue"): - self.graph = graph - self.quality_method = quality_method - np.random.seed(42) - - def __call__(self, freq): - L = construct_so3_laplacian(to_complex(freq), self.graph, abelian_scale=1.0) - return laplacian_quality(L) - - -def scan_frequencies_so3(graph, quality_method="eigenvalue"): - """Scan a range of complex frequencies and return mode qualities.""" - ks, alphas = get_scan_grid(graph) - freqs = [[k, a] for k in ks for a in alphas] - - worker_scan = WorkerScan(graph, quality_method=quality_method) - chunksize = max(1, int(0.1 * len(freqs) / graph.graph["params"]["n_workers"])) - with multiprocessing.Pool(graph.graph["params"]["n_workers"]) as pool: - qualities_list = list( - tqdm( - pool.imap(worker_scan, freqs, chunksize=chunksize), - total=len(freqs), - ) - ) - - id_k = [k_i for k_i in range(len(ks)) for a_i in range(len(alphas))] - id_a = [a_i for k_i in range(len(ks)) for a_i in range(len(alphas))] - qualities = sc.sparse.coo_matrix( - (qualities_list, (id_k, id_a)), - shape=(graph.graph["params"]["k_n"], graph.graph["params"]["alpha_n"]), - ).toarray() - - return qualities diff --git a/netsalt/physics.py b/netsalt/physics.py index 47e536a..7878f62 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 diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index f1366c3..802616f 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -455,5 +455,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 + ) return laplacian_quality(laplacian, method=quality_method) From 4482aa8a0673a1c55f4c8b5e63340201dbf6896c Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 9 Aug 2023 13:48:03 +0200 Subject: [PATCH 07/22] clean --- examples/non_abelian/test.py | 112 ------------------------ examples/transfer/boundary_example.py | 78 ----------------- examples/transfer/transfer_example.py | 119 -------------------------- 3 files changed, 309 deletions(-) delete mode 100644 examples/non_abelian/test.py delete mode 100644 examples/transfer/boundary_example.py delete mode 100644 examples/transfer/transfer_example.py diff --git a/examples/non_abelian/test.py b/examples/non_abelian/test.py deleted file mode 100644 index 5dcb54e..0000000 --- a/examples/non_abelian/test.py +++ /dev/null @@ -1,112 +0,0 @@ -import numpy as np -from tqdm import tqdm -import matplotlib.pyplot as plt -from netsalt.quantum_graph import create_quantum_graph, laplacian_quality, mode_quality -from netsalt.modes import mode_on_nodes - -from netsalt.modes import scan_frequencies -from netsalt.plotting import plot_scan -from netsalt.physics import dispersion_relation_linear, set_dispersion_relation -import networkx as nx - - -from netsalt.non_abelian import construct_so3_laplacian, so3_mode_on_nodes, scan_frequencies_so3 - - -def make_graph(n): - graph = nx.cycle_graph(n) - graph.add_edge(0, 8) - graph.add_edge(0, 20) - graph.add_edge(10, 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(1, n) - graph.add_edge(15, n + 1) - pos += [[1.2, 0]] - pos += [[-1.2, 0]] - - return graph, pos - - -if __name__ == "__main__": - - params = { - "open_model": "open", - "c": 1.0, - "k_min": 10.0, - "k_max": 12.0, - "k_n": 200, - "alpha_min": 0.0, - "alpha_max": 0.4, - "alpha_n": 50, - "n_workers": 7, - } - n = 30 - graph, pos = make_graph(n) - graph_u1, pos = make_graph(n) - - nx.draw(graph, pos=pos) - - create_quantum_graph(graph, params=params, positions=pos) - create_quantum_graph(graph_u1, params=params, positions=pos) - - set_dispersion_relation(graph_u1, dispersion_relation_linear) - - qualities_u1 = scan_frequencies(graph_u1) - plot_scan(graph_u1, qualities_u1) - plt.suptitle('u1') - - qualities = scan_frequencies_so3(graph) - plot_scan(graph, qualities) - plt.suptitle("so3") - plt.show() - - -def lkj(): - ks = np.linspace(10.0, 15, 500) - qs = [] - qs_u1 = [] - for k in tqdm(ks): - kim = 0.05 - L = construct_so3_laplacian(k + 1j * kim, graph) - qs.append(laplacian_quality(L)) - qs_u1.append(mode_quality([k, kim], graph_u1)) - - plt.figure() - plt.plot(ks, qs_u1, "+-r", label="u1") - plt.plot(ks, qs, "-", label="so3") - plt.legend(loc="best") - plt.yscale("log") - plt.show() - - k = ks[np.argmin(qs)] - L = construct_so3_laplacian(k, graph) - mode = so3_mode_on_nodes(L) - - k_u1 = k # ks[np.argmin(qs_u1)] - mode_u1 = mode_on_nodes([k_u1, 0], graph_u1) - - plt.figure() - x = np.abs(mode[::3]) - y = np.abs(mode[1::3]) - z = np.abs(mode[2::3]) - plt.plot(x, label="x") - plt.plot(y, label="y") - plt.plot(z, label="z") - n = np.sqrt(x**2 + y**2 + z**2) - - plt.plot(n, "+-", label="norm") - plt.plot(np.abs(mode_u1), label="u1") - plt.legend() - - plt.figure() - plt.plot(ks, qs_u1, "+-r") - plt.plot(ks, qs, "-") - plt.axvline(k) - plt.axvline(k_u1, c="r") - plt.yscale("log") - - plt.figure() - plt.plot(np.sqrt(x**2 + y**2 + z**2) - np.abs(mode_u1)) - plt.show() diff --git a/examples/transfer/boundary_example.py b/examples/transfer/boundary_example.py deleted file mode 100644 index 0c62107..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": 0.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() From e5d8cc2cb0f849a18a86ff40cf7a3986156bb139 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 10 Aug 2023 10:53:19 +0200 Subject: [PATCH 08/22] add back tansfer --- examples/transfer/boundary_example.py | 78 +++++++++++++++++ examples/transfer/transfer_example.py | 119 ++++++++++++++++++++++++++ 2 files changed, 197 insertions(+) create mode 100644 examples/transfer/boundary_example.py create mode 100644 examples/transfer/transfer_example.py diff --git a/examples/transfer/boundary_example.py b/examples/transfer/boundary_example.py new file mode 100644 index 0000000..70790d4 --- /dev/null +++ b/examples/transfer/boundary_example.py @@ -0,0 +1,78 @@ +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 new file mode 100644 index 0000000..5b3417a --- /dev/null +++ b/examples/transfer/transfer_example.py @@ -0,0 +1,119 @@ +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() From 6d875d5b16e2bb20abd00618e56c62818d148e86 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 10 Aug 2023 15:25:41 +0200 Subject: [PATCH 09/22] fix transfer --- netsalt/modes.py | 14 +++++++++++--- netsalt/plotting.py | 24 ++++++++++++++++++++---- 2 files changed, 31 insertions(+), 7 deletions(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 968fff9..67a5fa8 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -296,7 +296,10 @@ def flux_on_edges(mode, graph): def mean_mode_on_edges(mode, graph): r"""Compute the average :math:`Real(E^2)` on each edge.""" edge_flux = flux_on_edges(mode, graph) + return mean_on_edges(edge_flux, graph) + +def mean_on_edges(edge_flux, graph): mean_edge_solution = np.zeros(len(graph.edges)) for ei in range(len(graph.edges)): k = 1.0j * graph.graph["ks"][ei] @@ -898,14 +901,19 @@ 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) + BT, _ = construct_incidence_matrix(graph) + bt = np.clip(np.ceil(np.real(BT.toarray())), -1, 0) + K = bt.dot(np.repeat(graph.graph["ks"], 2)) + return sc.sparse.linalg.spsolve(construct_laplacian(k, graph), 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)) + _r = sc.sparse.linalg.spsolve( + construct_laplacian(k, graph), BT.dot(np.repeat(graph.graph["ks"], 2) * input_flow) + ) Winv = construct_weight_matrix(graph, with_k=False) return Winv.dot(B).dot(_r) @@ -926,7 +934,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/plotting.py b/netsalt/plotting.py index 2e5cda0..38cce10 100644 --- a/netsalt/plotting.py +++ b/netsalt/plotting.py @@ -441,24 +441,40 @@ def plot_single_mode( def _plot_single_mode( graph, mode, ax=None, colorbar=True, edge_vmin=None, edge_vmax=None, cmap="coolwarm" ): - positions = [graph.nodes[u]["position"] for u in graph] edge_solution = mean_mode_on_edges(mode, graph) + return plot_on_graph( + graph, + edge_solution, + edge_vmin=edge_vmin, + edge_vmax=edge_vmax, + ax=ax, + cmap=cmap, + colorbar=colorbar, + ) + + +def plot_on_graph( + graph, edge_data, edge_vmin=None, edge_vmax=None, ax=None, cmap="coolwarm", colorbar=True +): + """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)) + edge_vmin = -max(abs(edge_data)) + print(len(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, From fef93f829026bff3fc305d1962fef049120e73ac Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 11 Aug 2023 13:54:42 +0200 Subject: [PATCH 10/22] wip non-abelian transfer --- netsalt/modes.py | 36 +++++++++++++++++++++++++++--------- netsalt/non_abelian.py | 2 +- netsalt/quantum_graph.py | 4 ++-- 3 files changed, 30 insertions(+), 12 deletions(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 67a5fa8..14ebaca 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -264,7 +264,9 @@ 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" ) @@ -901,20 +903,36 @@ def lasing_threshold_linear(mode, graph, D0): def get_node_transfer(k, graph, input_flow): """Compute node transfer from a given input flow.""" - BT, _ = construct_incidence_matrix(graph) + laplacian, BT, _, _ = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + k, graph + ) bt = np.clip(np.ceil(np.real(BT.toarray())), -1, 0) - K = bt.dot(np.repeat(graph.graph["ks"], 2)) - return sc.sparse.linalg.spsolve(construct_laplacian(k, graph), K * input_flow) + 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 = sc.sparse.linalg.spsolve( - construct_laplacian(k, graph), BT.dot(np.repeat(graph.graph["ks"], 2) * input_flow) - ) - Winv = construct_weight_matrix(graph, with_k=False) + laplacian, BT, B, Winv = graph.graph["params"].get( + "laplacian_constructor", construct_laplacian + )(k, graph) + s = list(np.shape(graph.graph["ks"])) + s[0] *= 2 + ks = np.empty(s, dtype=np.complex) + ks[::2] = graph.graph["ks"] + ks[1::2] = graph.graph["ks"] + if len(s) > 1: + _in = np.einsum("ikl,ij", ks, np.diag(input_flow)) + _in = BT.dot(sc.linalg.block_diag(*_in)) + # WIP HERE + else: + _in = BT.dot(ks * input_flow) + + _r = sc.sparse.linalg.spsolve(laplacian, _in) return Winv.dot(B).dot(_r) diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 6866554..3db02a8 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -125,4 +125,4 @@ def construct_so3_laplacian(wavenumber, graph, abelian_scale=1.0): 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) - return BT.dot(Winv).dot(B) + return BT.dot(Winv).dot(B), BT, B, Winv diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index 802616f..fa296b8 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -275,7 +275,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): @@ -455,7 +455,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 = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( + laplacian = graph.graph["params"].get("laplacian_constructor", construct_laplacian)[0]( to_complex(mode), graph ) return laplacian_quality(laplacian, method=quality_method) From 100ddfee93342882b057a9ccaa3efafb6ef49ab4 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 18 Aug 2023 14:23:13 +0200 Subject: [PATCH 11/22] fix --- netsalt/modes.py | 1 + netsalt/quantum_graph.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 14ebaca..2bde73f 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -302,6 +302,7 @@ def mean_mode_on_edges(mode, graph): def mean_on_edges(edge_flux, graph): + """Mean edge values from any edge vector.""" mean_edge_solution = np.zeros(len(graph.edges)) for ei in range(len(graph.edges)): k = 1.0j * graph.graph["ks"][ei] diff --git a/netsalt/quantum_graph.py b/netsalt/quantum_graph.py index fa296b8..512a63c 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -455,7 +455,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 = graph.graph["params"].get("laplacian_constructor", construct_laplacian)[0]( + laplacian = graph.graph["params"].get("laplacian_constructor", construct_laplacian)( to_complex(mode), graph - ) + )[0] return laplacian_quality(laplacian, method=quality_method) From 7b0b7c4ff50554d35515302585fdebac1de5d760 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 29 Aug 2023 16:43:08 +0200 Subject: [PATCH 12/22] more --- netsalt/modes.py | 99 ++++++++++++++++++++++++++++++---------- netsalt/non_abelian.py | 28 ++++++------ netsalt/plotting.py | 43 ++++++++++++++--- netsalt/quantum_graph.py | 9 +++- 4 files changed, 131 insertions(+), 48 deletions(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 2bde73f..8f3c322 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -23,6 +23,7 @@ mode_quality, ) from .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 +188,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): @@ -268,7 +268,7 @@ 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(len(graph)), which="LM" + laplacian, k=1, sigma=0, v0=np.ones(laplacian.shape[0]), which="LM" ) quality_thresh = graph.graph["params"].get("quality_threshold", 1e-2) if abs(min_eigenvalue[0]) > quality_thresh: @@ -288,39 +288,88 @@ 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) + 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. -def mean_on_edges(edge_flux, graph): - """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) + + 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] - 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]))) + 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 norm_type.startswith("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: + 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: + + w_perp = Ad(2.0 * length * chi).dot(proj_perp(chi)) + w_paral = np.exp(2.0j * length * norm(chi)) * proj_paral(chi) + 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) + + if DIM == 1: + _z = Winv[_ext(2 * ei), _ext(2 * ei)].toarray() + + z[_ext(0), _ext(0)] = z[_ext(1), _ext(1)] = _z + 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 diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 3db02a8..3ad25ec 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -2,6 +2,8 @@ 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.""" @@ -48,8 +50,8 @@ def Ad(chi_mat): def set_so3_wavenumber(graph, wavenumber): """Set so3 matrix wavenumber.""" - chis = graph.graph["params"].get("chis", None) - if chis is None: + chis = [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: @@ -60,12 +62,11 @@ def set_so3_wavenumber(graph, wavenumber): def construct_so3_incidence_matrix(graph, abelian_scale=1.0): """Construct SO3 incidence matrix.""" - DIM = 3 def _ext(i): return slice(DIM * i, DIM * (i + 1)) - Bout = sparse.lil_matrix((len(graph.edges) * 2 * DIM, len(graph) * DIM), dtype=np.complex128) + 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): one = np.eye(DIM) @@ -79,22 +80,21 @@ def _ext(i): out = len(graph[u]) == 1 or len(graph[v]) == 1 - Bout[_ext(2 * ei), _ext(u)] = -one - Bout[_ext(2 * ei), _ext(v)] = 0 if out else expl - Bout[_ext(2 * ei + 1), _ext(u)] = 0 if out else expl - Bout[_ext(2 * ei + 1), _ext(v)] = -one + 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)] = 0 if out else expl + BT[_ext(u), _ext(2 * ei + 1)] = 0 if out else expl BT[_ext(v), _ext(2 * ei + 1)] = -one - return BT, Bout + return BT, B def construct_so3_weight_matrix(graph, with_k=True, abelian_scale=1.0): """Construct SO3 weight matrix.""" - DIM = 3 def _ext(i): return slice(DIM * i, DIM * (i + 1)) @@ -120,9 +120,9 @@ def _ext(i): return Winv -def construct_so3_laplacian(wavenumber, graph, abelian_scale=1.0): +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) + 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/plotting.py b/netsalt/plotting.py index 38cce10..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,9 +447,16 @@ 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", ): - 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, @@ -451,11 +466,19 @@ def _plot_single_mode( 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 + 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: @@ -469,8 +492,11 @@ def plot_on_graph( if edge_vmax is None: edge_vmax = max(abs(edge_data)) if edge_vmin is None: - edge_vmin = -max(abs(edge_data)) - print(len(edge_data)) + 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, @@ -484,7 +510,10 @@ def plot_on_graph( 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 512a63c..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) From 8ae70464d18cfe9d2761334638532710136ce7e1 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 29 Aug 2023 17:07:17 +0200 Subject: [PATCH 13/22] lint --- netsalt/modes.py | 1 - 1 file changed, 1 deletion(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 8f3c322..a88b32b 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -347,7 +347,6 @@ 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) if norm_type.endswith("x"): From 01fc1e4c8aa9a6e5eb1233c3efdbd5c5d43e8a0b Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 5 Sep 2023 11:25:40 +0200 Subject: [PATCH 14/22] more --- .pylintrc | 2 +- netsalt/modes.py | 28 +++++++++++++++------------- 2 files changed, 16 insertions(+), 14 deletions(-) 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/netsalt/modes.py b/netsalt/modes.py index a88b32b..96f1b49 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -9,20 +9,19 @@ 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") @@ -331,6 +330,8 @@ def mean_on_edges(edge_flux, graph, norm_type="abs", mode=None): ) 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 ) @@ -965,24 +966,25 @@ def get_node_transfer(k, graph, input_flow): def get_edge_transfer(k, graph, input_flow): """Compute edge transfer from a given input flow.""" - set_wavenumber(graph, k) laplacian, BT, B, Winv = graph.graph["params"].get( "laplacian_constructor", construct_laplacian - )(k, graph) - s = list(np.shape(graph.graph["ks"])) - s[0] *= 2 - ks = np.empty(s, dtype=np.complex) + )(k, graph, with_k=True) + + shape_k = list(np.shape(graph.graph["ks"])) + shape_k[0] *= 2 + ks = np.empty(shape_k, dtype=np.complex) ks[::2] = graph.graph["ks"] ks[1::2] = graph.graph["ks"] - if len(s) > 1: + + if len(shape_k) > 1: + # non-abelian part, WIP _in = np.einsum("ikl,ij", ks, np.diag(input_flow)) _in = BT.dot(sc.linalg.block_diag(*_in)) - # WIP HERE else: _in = BT.dot(ks * input_flow) _r = sc.sparse.linalg.spsolve(laplacian, _in) - return Winv.dot(B).dot(_r) + return Winv.dot(B).dot(_r) / ks # we divide by ks as Winv has it (with_k=True) def estimate_boundary_flow(graph, input_flow, k_frac=1e-2): From 43a4cf6882089dd015c456f6913c20c340c484e1 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 11 Sep 2023 14:41:23 +0200 Subject: [PATCH 15/22] more --- netsalt/modes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index 96f1b49..3a806e1 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -987,7 +987,7 @@ def get_edge_transfer(k, graph, input_flow): return Winv.dot(B).dot(_r) / ks # we divide by ks as Winv has it (with_k=True) -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: From d850d63ba3bd4ba0a6116e680dc03308378bfcca Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 14 Sep 2023 11:29:35 +0200 Subject: [PATCH 16/22] more --- netsalt/algorithm.py | 1 + netsalt/modes.py | 31 +++++++++++++++++++------------ 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/netsalt/algorithm.py b/netsalt/algorithm.py index 92db27e..887eca8 100644 --- a/netsalt/algorithm.py +++ b/netsalt/algorithm.py @@ -70,6 +70,7 @@ def refine_mode_brownian_ratchet( ) new_quality = mode_quality(new_mode, graph, quality_method=quality_method) + print(current_mode, current_quality, new_mode, new_quality) if disp: L.debug( diff --git a/netsalt/modes.py b/netsalt/modes.py index 3a806e1..f60e67d 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -966,25 +966,32 @@ def get_node_transfer(k, graph, input_flow): def get_edge_transfer(k, graph, input_flow): """Compute edge transfer from a given input flow.""" - laplacian, BT, B, Winv = graph.graph["params"].get( - "laplacian_constructor", construct_laplacian - )(k, graph, with_k=True) + 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"])) - shape_k[0] *= 2 - ks = np.empty(shape_k, dtype=np.complex) - ks[::2] = graph.graph["ks"] - ks[1::2] = graph.graph["ks"] if len(shape_k) > 1: - # non-abelian part, WIP - _in = np.einsum("ikl,ij", ks, np.diag(input_flow)) - _in = BT.dot(sc.linalg.block_diag(*_in)) + 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: - _in = BT.dot(ks * input_flow) + 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) / ks # we divide by ks as Winv has it (with_k=True) + return Winv.dot(B).dot(_r) def estimate_boundary_flow(graph, input_flow, k_frac=1e-3): From 0c45720cc6d011aefd8a1e6ff9d526662f00e970 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 15 Sep 2023 14:02:19 +0200 Subject: [PATCH 17/22] more --- netsalt/algorithm.py | 1 - netsalt/non_abelian.py | 19 +++++++++++++++---- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/netsalt/algorithm.py b/netsalt/algorithm.py index 887eca8..92db27e 100644 --- a/netsalt/algorithm.py +++ b/netsalt/algorithm.py @@ -70,7 +70,6 @@ def refine_mode_brownian_ratchet( ) new_quality = mode_quality(new_mode, graph, quality_method=quality_method) - print(current_mode, current_quality, new_mode, new_quality) if disp: L.debug( diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 3ad25ec..fc4e371 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -78,18 +78,29 @@ def _ext(i): * np.exp(1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei])) ) - out = len(graph[u]) == 1 or len(graph[v]) == 1 - 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)] = 0 if out else expl - BT[_ext(u), _ext(2 * ei + 1)] = 0 if out else expl + 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 From 944a9d2aa24d223e302773ec60f3db44db3ab468 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 16 May 2024 12:51:57 +0200 Subject: [PATCH 18/22] tmp --- netsalt/non_abelian.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index fc4e371..2ee812f 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -117,9 +117,9 @@ def _ext(i): chi = graph.graph["ks"][ei] length = graph.graph["lengths"][ei] - 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) + 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 winv = linalg.inv(w) From 2c61e321e555dbd825156aefcb744f6e43991909 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Aug 2024 08:59:52 +0200 Subject: [PATCH 19/22] 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): From ce939a03068b82dad65db8487fd16cc0dd65c0a8 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Aug 2024 15:03:11 +0200 Subject: [PATCH 20/22] open example --- examples/non_abelian/open/abelian_open.py | 68 ++++++++++ examples/non_abelian/open/make_graph.py | 24 ++++ examples/non_abelian/open/so3_open.py | 109 ++++++++++++++++ examples/non_abelian/open/so3_open_uniform.py | 122 ++++++++++++++++++ netsalt/modes.py | 38 +++--- 5 files changed, 339 insertions(+), 22 deletions(-) create mode 100644 examples/non_abelian/open/abelian_open.py create mode 100644 examples/non_abelian/open/make_graph.py create mode 100644 examples/non_abelian/open/so3_open.py create mode 100644 examples/non_abelian/open/so3_open_uniform.py 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.py b/examples/non_abelian/open/so3_open.py new file mode 100644 index 0000000..cf07fc0 --- /dev/null +++ b/examples/non_abelian/open/so3_open.py @@ -0,0 +1,109 @@ +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.3, + "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") 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..6b09e61 --- /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.3, + "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") + + plt.figure(figsize=(5, 3)) + plot_single_mode(over_graph, modes_df, 7, 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, 7, 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, 7, 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, 7, norm_type="real_z", ax=plt.gca()) + plt.tight_layout() + plt.savefig("so3_open_mode_3_z.pdf") diff --git a/netsalt/modes.py b/netsalt/modes.py index 8a79394..cb8220b 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -267,10 +267,9 @@ def mode_on_nodes(mode, graph): to_complex(mode), graph )[0] min_eigenvalue, node_solution = sc.sparse.linalg.eigs( - laplacian, k=4, sigma=0, v0=np.ones(laplacian.shape[0]), which="LM" + laplacian, k=5, 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: " @@ -280,8 +279,13 @@ def mode_on_nodes(mode, graph): + ", mode: " + str(mode) ) - - return node_solution[:, 0] + # 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 + n_eig = len([1 for a in np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue)) if a]) + 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[:, np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue))].sum(axis=1) def flux_on_edges(mode, graph): @@ -369,29 +373,19 @@ def sol(x): 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) - 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)]) - ) - # print("check:", mean_edge_solution[ei], res) - mean_edge_solution[ei] = res + 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 From 36c2b288142aab7a98882199d5adcfaac7fba183 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 7 Aug 2024 10:18:06 +0200 Subject: [PATCH 21/22] more --- .../{so3_closed.py => so3_closed_random.py} | 41 ++++++++++--------- .../non_abelian/closed/so3_closed_uniform.py | 34 +++++++-------- .../open/{so3_open.py => so3_open_random.py} | 23 ++++++++++- examples/non_abelian/open/so3_open_uniform.py | 12 +++--- netsalt/non_abelian.py | 23 +++++++---- 5 files changed, 82 insertions(+), 51 deletions(-) rename examples/non_abelian/closed/{so3_closed.py => so3_closed_random.py} (81%) rename examples/non_abelian/open/{so3_open.py => so3_open_random.py} (82%) diff --git a/examples/non_abelian/closed/so3_closed.py b/examples/non_abelian/closed/so3_closed_random.py similarity index 81% rename from examples/non_abelian/closed/so3_closed.py rename to examples/non_abelian/closed/so3_closed_random.py index 3142fe3..fd7d7bc 100644 --- a/examples/non_abelian/closed/so3_closed.py +++ b/examples/non_abelian/closed/so3_closed_random.py @@ -35,12 +35,16 @@ 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): - 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]) + 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(): @@ -69,43 +73,42 @@ 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.savefig("so3_close_mode_1.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real") + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_x") plt.tight_layout() - plt.savefig(f"so3_close_mode_2_{i}.pdf") + 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_x") + 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_x_{i}.pdf") + plt.savefig("so3_close_mode_1_y.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_x") + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") plt.tight_layout() - plt.savefig(f"so3_close_mode_2_x_{i}.pdf") + plt.savefig("so3_close_mode_1_z.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_y") + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real") plt.tight_layout() - plt.savefig(f"so3_close_mode_1_y_{i}.pdf") + 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_y") + 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_y_{i}.pdf") + plt.savefig("so3_close_mode_2_x.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_y") plt.tight_layout() - plt.savefig(f"so3_close_mode_1_z_{i}.pdf") + 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(f"so3_close_mode_2_z_{i}.pdf") + 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 index 507d69d..4c66f2f 100644 --- a/examples/non_abelian/closed/so3_closed_uniform.py +++ b/examples/non_abelian/closed/so3_closed_uniform.py @@ -58,47 +58,47 @@ 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.savefig("so3_close_mode_1_uniform.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real") + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_x") plt.tight_layout() - plt.savefig(f"so3_close_mode_2_uniform_{i}.pdf") + 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_x") + 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_x_uniform_{i}.pdf") + plt.savefig("so3_close_mode_1_y_uniform.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_x") + plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") plt.tight_layout() - plt.savefig(f"so3_close_mode_2_x_uniform_{i}.pdf") + plt.savefig("so3_close_mode_1_z_uniform.pdf") plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_y") + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real") plt.tight_layout() - plt.savefig(f"so3_close_mode_1_y_uniform_{i}.pdf") + 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_y") + 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_y_uniform_{i}.pdf") + plt.savefig("so3_close_mode_2_x_uniform.pdf") + plt.figure(figsize=(4, 3)) - plot_single_mode(over_graph, modes_df, 1, ax=plt.gca(), norm_type="real_z") + plot_single_mode(over_graph, modes_df, 2, ax=plt.gca(), norm_type="real_y") plt.tight_layout() - plt.savefig(f"so3_close_mode_1_z_uniform_{i}.pdf") + 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(f"so3_close_mode_2_z_uniform_{i}.pdf") - - #plt.show() + plt.savefig("so3_close_mode_2_z_uniform.pdf") diff --git a/examples/non_abelian/open/so3_open.py b/examples/non_abelian/open/so3_open_random.py similarity index 82% rename from examples/non_abelian/open/so3_open.py rename to examples/non_abelian/open/so3_open_random.py index cf07fc0..ec45da0 100644 --- a/examples/non_abelian/open/so3_open.py +++ b/examples/non_abelian/open/so3_open_random.py @@ -25,7 +25,7 @@ "k_min": 5.0, "k_max": 7.0, "alpha_n": 100, - "alpha_min": -0.3, + "alpha_min": -0.01, "alpha_max": 0.3, "c": len(graph.edges) * [1.0], "laplacian_constructor": construct_so3_laplacian, @@ -107,3 +107,24 @@ 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 index 6b09e61..72e694e 100644 --- a/examples/non_abelian/open/so3_open_uniform.py +++ b/examples/non_abelian/open/so3_open_uniform.py @@ -24,7 +24,7 @@ "k_min": 5.0, "k_max": 7.0, "alpha_n": 100, - "alpha_min": -0.3, + "alpha_min": -0.01, "alpha_max": 0.3, "c": len(graph.edges) * [1.0], "laplacian_constructor": construct_so3_laplacian, @@ -100,23 +100,23 @@ 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, 7, norm_type="real", ax=plt.gca()) + 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, 7, norm_type="real_x", ax=plt.gca()) + 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, 7, norm_type="real_y", ax=plt.gca()) + 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, 7, norm_type="real_z", ax=plt.gca()) + 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/netsalt/non_abelian.py b/netsalt/non_abelian.py index 1675b60..00d303f 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -55,7 +55,9 @@ def set_so3_wavenumber(graph, wavenumber): else: if len(np.shape(chis[0])) == 1: chis = np.array([hat_inv(chi) for chi in chis]) - graph.graph["ks"] = chis * wavenumber + graph.graph["ks"] = chis * np.real(wavenumber) - np.eye(3) * np.imag(wavenumber) + graph.graph["chis"] = chis + graph.graph["wavenumber"] = wavenumber def construct_so3_incidence_matrix(graph, abelian_scale=1.0): @@ -69,11 +71,11 @@ def _ext(i): for ei, (u, v) in enumerate(graph.edges): one = np.eye(DIM) expl = Ad(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) - expl = np.array(expl.dot(proj_perp(graph.graph["ks"][ei])), dtype=np.complex128) + expl = np.array(expl.dot(proj_perp(graph.graph["chis"][ei])), dtype=np.complex128) expl += ( abelian_scale - * np.exp(1.0j * graph.graph["lengths"][ei] * norm(graph.graph["ks"][ei])) - * proj_paral(graph.graph["ks"][ei]) + * np.exp(1.0j * graph.graph["lengths"][ei] * graph.graph["wavenumber"]) + * proj_paral(graph.graph["chis"][ei]) ) B[_ext(2 * ei), _ext(u)] = -one @@ -112,17 +114,22 @@ def _ext(i): (len(graph.edges) * 2 * DIM, len(graph.edges) * 2 * DIM), dtype=np.complex128 ) for ei, _ in enumerate(graph.edges): - chi = graph.graph["ks"][ei] + k = graph.graph["ks"][ei] + chi = graph.graph["chis"][ei] length = graph.graph["lengths"][ei] - 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_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 = (chi.dot(proj_perp(chi)) + 1.0j * norm(chi) * proj_paral(chi)).dot(winv) + 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 From 29895b6a63b525147b902b18607bca168c7c1e5e Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 30 Oct 2024 09:35:13 +0100 Subject: [PATCH 22/22] more --- netsalt/modes.py | 10 ++++++++-- netsalt/non_abelian.py | 20 +++++++++++--------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/netsalt/modes.py b/netsalt/modes.py index cb8220b..7eb8bcf 100644 --- a/netsalt/modes.py +++ b/netsalt/modes.py @@ -281,11 +281,17 @@ def mode_on_nodes(mode, graph): ) # 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 - n_eig = len([1 for a in np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue)) if a]) + mask = np.logical_or( + np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue)), + np.abs(min_eigenvalue) < 1e-4, + ) + + 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[:, np.abs(min_eigenvalue) < 1.5 * min(np.abs(min_eigenvalue))].sum(axis=1) + return node_solution[:, mask].sum(axis=1) def flux_on_edges(mode, graph): diff --git a/netsalt/non_abelian.py b/netsalt/non_abelian.py index 00d303f..b6745bd 100644 --- a/netsalt/non_abelian.py +++ b/netsalt/non_abelian.py @@ -48,14 +48,15 @@ def Ad(chi_mat): def set_so3_wavenumber(graph, wavenumber): """Set so3 matrix wavenumber.""" - chis = [graph[u][v].get("chi", None) for u, v in graph.edges] + 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]) - graph.graph["ks"] = chis * np.real(wavenumber) - np.eye(3) * np.imag(wavenumber) + 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 @@ -69,14 +70,15 @@ def _ext(i): 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(graph.graph["lengths"][ei] * graph.graph["ks"][ei]) - expl = np.array(expl.dot(proj_perp(graph.graph["chis"][ei])), dtype=np.complex128) - expl += ( - abelian_scale - * np.exp(1.0j * graph.graph["lengths"][ei] * graph.graph["wavenumber"]) - * proj_paral(graph.graph["chis"][ei]) - ) + 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