From c9332b8868b984713a2e080d7255a52676a9f41a Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 11 Aug 2023 13:54:42 +0200 Subject: [PATCH] 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 d74b150..f6cd252 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 ff4991e..3572e8d 100644 --- a/netsalt/quantum_graph.py +++ b/netsalt/quantum_graph.py @@ -267,7 +267,7 @@ def construct_laplacian(wavenumber, graph): set_wavenumber(graph, wavenumber) BT, B = construct_incidence_matrix(graph) Winv = construct_weight_matrix(graph) - return BT.dot(Winv).dot(B) + return BT.dot(Winv).dot(B), BT, B, Winv def set_wavenumber(graph, wavenumber): @@ -447,7 +447,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)