From fe20adcdb3bc1a9bcc735a5a3072888547af92c5 Mon Sep 17 00:00:00 2001 From: d-schindler <60650591+d-schindler@users.noreply.github.com> Date: Fri, 8 Nov 2024 09:46:26 +0000 Subject: [PATCH] use linalg.eigh instead of linalg.eig --- src/pygenstability/constructors.py | 59 ++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 16 deletions(-) diff --git a/src/pygenstability/constructors.py b/src/pygenstability/constructors.py index 743d2ca..7b7852b 100644 --- a/src/pygenstability/constructors.py +++ b/src/pygenstability/constructors.py @@ -31,9 +31,13 @@ def load_constructor(constructor, graph, **kwargs): """Load a constructor from its name, or as a custom Constructor class.""" if isinstance(constructor, str): if graph is None: - raise Exception(f"No graph was provided with a generic constructor {constructor}") + raise Exception( + f"No graph was provided with a generic constructor {constructor}" + ) try: - return getattr(sys.modules[__name__], f"constructor_{constructor}")(graph, **kwargs) + return getattr(sys.modules[__name__], f"constructor_{constructor}")( + graph, **kwargs + ) except AttributeError as exc: raise Exception(f"Could not load constructor {constructor}") from exc if not isinstance(constructor, Constructor): @@ -53,8 +57,8 @@ def limit(*args, **kwargs): def _compute_spectral_decomp(matrix): - """Solve eigenalue problem.""" - lambdas, v = la.eig(matrix.toarray()) + """Solve eigenalue problem for symmetric matrix.""" + lambdas, v = la.eigh(matrix.toarray()) vinv = la.inv(v.real) return lambdas.real, v.real, vinv @@ -67,7 +71,9 @@ def _check_total_degree(degrees): def _get_spectral_gap(laplacian): """Compute spectral gap.""" - spectral_gap = np.round(max(np.real(sp.linalg.eigs(laplacian, which="SM", k=2)[0])), 8) + spectral_gap = np.round( + max(np.real(sp.linalg.eigs(laplacian, which="SM", k=2)[0])), 8 + ) L.info("Spectral gap = 10^%s", np.around(np.log10(spectral_gap), 2)) return spectral_gap @@ -80,7 +86,9 @@ class Constructor: to return quality matrix, null model, and possible global shift. """ - def __init__(self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **kwargs): + def __init__( + self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **kwargs + ): """The constructor calls the prepare method upon initialisation. Args: @@ -105,7 +113,9 @@ def __init__(self, graph, with_spectral_gap=False, exp_comp_mode="spectral", **k def _get_exp(self, scale): """Compute matrix exponential at a given scale.""" if self.exp_comp_mode == "expm": - exp = sp.linalg.expm(-scale * self.partial_quality_matrix.toarray().astype(_DTYPE)) + exp = sp.linalg.expm( + -scale * self.partial_quality_matrix.toarray().astype(_DTYPE) + ) if self.exp_comp_mode == "spectral": lambdas, v, vinv = self.spectral_decomp exp = v @ np.diag(np.exp(-scale * lambdas)) @ vinv @@ -181,7 +191,9 @@ class constructor_continuous_combinatorial(Constructor): @_limit_numpy def prepare(self, **kwargs): """Prepare the constructor with non-scale dependent computations.""" - laplacian, degrees = sp.csgraph.laplacian(self.graph, return_diag=True, normed=False) + laplacian, degrees = sp.csgraph.laplacian( + self.graph, return_diag=True, normed=False + ) _check_total_degree(degrees) laplacian /= degrees.mean() pi = np.ones(self.graph.shape[0]) / self.graph.shape[0] @@ -221,7 +233,9 @@ class constructor_continuous_normalized(Constructor): @_limit_numpy def prepare(self, **kwargs): """Prepare the constructor with non-scale dependent computations.""" - laplacian, degrees = sp.csgraph.laplacian(self.graph, return_diag=True, normed=False) + laplacian, degrees = sp.csgraph.laplacian( + self.graph, return_diag=True, normed=False + ) _check_total_degree(degrees) normed_laplacian = sp.diags(1.0 / degrees).dot(laplacian) @@ -354,7 +368,9 @@ class constructor_directed(Constructor): @_limit_numpy def prepare(self, **kwargs): """Prepare the constructor with non-scale dependent computations.""" - assert self.exp_comp_mode == "expm" + assert ( + self.exp_comp_mode == "expm" + ), 'exp_comp_mode="expm" is required for "constructor_directed"' alpha = kwargs.get("alpha", 0.8) n_nodes = self.graph.shape[0] @@ -366,11 +382,17 @@ def prepare(self, **kwargs): self.partial_quality_matrix = sp.csr_matrix( alpha * np.diag(dinv).dot(self.graph.toarray()) - + ((1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0))).dot(ones) + + ( + (1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0)) + ).dot(ones) - np.eye(n_nodes) ) - pi = abs(sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][:, 0]) + pi = abs( + sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][ + :, 0 + ] + ) pi /= pi.sum() self.partial_null_model = np.array([pi, pi]) @@ -424,9 +446,10 @@ def prepare(self, **kwargs): self.partial_quality_matrix = sp.csr_matrix( alpha * np.diag(dinv).dot(self.graph.toarray()) - + ((1 - alpha) * np.diag(np.ones(n_nodes)) + np.diag(alpha * (dinv == 0.0))).dot( - ones - ) + + ( + (1 - alpha) * np.diag(np.ones(n_nodes)) + + np.diag(alpha * (dinv == 0.0)) + ).dot(ones) - np.eye(n_nodes) ) @@ -435,7 +458,11 @@ def prepare(self, **kwargs): sp.diags(dinv).dot(self.graph) - sp.diags(np.ones(n_nodes)) ) - pi = abs(sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][:, 0]) + pi = abs( + sp.linalg.eigs(self.partial_quality_matrix.transpose(), which="SM", k=1)[1][ + :, 0 + ] + ) pi /= pi.sum() self.partial_null_model = np.array([pi, pi])