Skip to content

Commit

Permalink
fix test
Browse files Browse the repository at this point in the history
  • Loading branch information
arnaudon committed Nov 14, 2024
1 parent c3be8d7 commit 52fa5c0
Show file tree
Hide file tree
Showing 21 changed files with 3,745 additions and 1,639 deletions.
51 changes: 13 additions & 38 deletions src/pygenstability/constructors.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,9 @@ 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):
Expand Down Expand Up @@ -71,9 +67,7 @@ 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

Expand All @@ -86,9 +80,7 @@ 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:
Expand All @@ -113,9 +105,7 @@ def __init__(
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
Expand Down Expand Up @@ -197,9 +187,7 @@ 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]
Expand Down Expand Up @@ -244,9 +232,7 @@ 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)

Expand Down Expand Up @@ -403,17 +389,11 @@ 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])

Expand Down Expand Up @@ -467,10 +447,9 @@ 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)
)

Expand All @@ -479,11 +458,7 @@ 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])

Expand Down
16 changes: 4 additions & 12 deletions src/pygenstability/optimal_scales.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,7 @@ def _pool2d_nvi(A, kernel_size, stride, padding=0):
return np.nanmean(A_w, axis=(2, 3))


def identify_optimal_scales(
results, kernel_size=3, window_size=3, max_nvi=1, basin_radius=1
):
def identify_optimal_scales(results, kernel_size=3, window_size=3, max_nvi=1, basin_radius=1):
"""Identifies optimal scales in Markov Stability [1]_.
Robust scales are found in a sequential way. We first search for large diagonal blocks
Expand Down Expand Up @@ -82,9 +80,7 @@ def identify_optimal_scales(

# smooth diagonal with moving window
block_nvi = np.roll(
np.asarray(
pd.Series(diagonal).rolling(window=window_size, win_type="triang").mean()
),
np.asarray(pd.Series(diagonal).rolling(window=window_size, win_type="triang").mean()),
-int(window_size / 2),
)
results["block_nvi"] = block_nvi
Expand All @@ -97,19 +93,15 @@ def identify_optimal_scales(

if (
np.count_nonzero(
np.around(
block_nvi[not_nan_ind[0] : not_nan_ind[0] + 2 * basin_radius + 1], 5
)
np.around(block_nvi[not_nan_ind[0] : not_nan_ind[0] + 2 * basin_radius + 1], 5)
)
== 0
):
basin_centers = np.insert(basin_centers, 0, not_nan_ind[0] + basin_radius)

if (
np.count_nonzero(
np.around(
block_nvi[not_nan_ind[-1] - 2 * basin_radius : not_nan_ind[-1] + 1], 5
)
np.around(block_nvi[not_nan_ind[-1] - 2 * basin_radius : not_nan_ind[-1] + 1], 5)
)
== 0
):
Expand Down
32 changes: 8 additions & 24 deletions src/pygenstability/plotting.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,7 @@
try:
import networkx as nx
except ImportError: # pragma: no cover
print(
'Please install networkx via pip install "pygenstability[networkx]" for full plotting.'
)
print('Please install networkx via pip install "pygenstability[networkx]" for full plotting.')

import numpy as np
from matplotlib import gridspec
Expand Down Expand Up @@ -52,9 +50,7 @@ def plot_scan(
plotly_filename (str): filename of .html figure from plotly
"""
if len(all_results["scales"]) == 1: # pragma: no cover
L.info(
"Cannot plot the results if only one scale point, we display the result instead:"
)
L.info("Cannot plot the results if only one scale point, we display the result instead:")
L.info(all_results)
return None

Expand Down Expand Up @@ -310,16 +306,12 @@ def plot_communities(
plot_single_partition(
graph, all_results, scale_id, edge_color=edge_color, edge_width=edge_width
)
plt.savefig(
os.path.join(folder, "scale_" + str(scale_id) + ext), bbox_inches="tight"
)
plt.savefig(os.path.join(folder, "scale_" + str(scale_id) + ext), bbox_inches="tight")
plt.close()
matplotlib.use(mpl_backend)


def plot_communities_matrix(
graph, all_results, folder="communities_matrix", ext=".pdf"
):
def plot_communities_matrix(graph, all_results, folder="communities_matrix", ext=".pdf"):
"""Plot communities at all scales in matrix form.
Args:
Expand Down Expand Up @@ -348,9 +340,7 @@ def plot_communities_matrix(
plt.plot((lines[i + 1], lines[i + 1]), (lines[i + 1], lines[i]), c="k")
plt.plot((lines[i + 1], lines[i]), (lines[i + 1], lines[i + 1]), c="k")

plt.savefig(
os.path.join(folder, "scale_" + str(scale_id) + ext), bbox_inches="tight"
)
plt.savefig(os.path.join(folder, "scale_" + str(scale_id) + ext), bbox_inches="tight")


def _get_scales(all_results, scale_axis=True):
Expand All @@ -364,19 +354,15 @@ def _get_scales(all_results, scale_axis=True):

def _plot_number_comm(all_results, ax, scales):
"""Plot number of communities."""
ax.plot(
scales, all_results["number_of_communities"], "-", c="C3", label="size", lw=2.0
)
ax.plot(scales, all_results["number_of_communities"], "-", c="C3", label="size", lw=2.0)
ax.set_ylim(0, 1.1 * max(all_results["number_of_communities"]))
ax.set_ylabel("# clusters", color="C3")
ax.tick_params("y", colors="C3")


def _plot_ttprime(all_results, ax, scales):
"""Plot ttprime."""
contourf_ = ax.contourf(
scales, scales, all_results["ttprime"], cmap="YlOrBr_r", extend="min"
)
contourf_ = ax.contourf(scales, scales, all_results["ttprime"], cmap="YlOrBr_r", extend="min")
ax.set_ylabel(r"$log_{10}(t^\prime)$")
ax.yaxis.tick_left()
ax.yaxis.set_label_position("left")
Expand Down Expand Up @@ -447,9 +433,7 @@ def _plot_optimal_scales(all_results, ax, scales, ax1, ax2):
ax2.axvline(scale, ls="--", color="C4")


def plot_scan_plt(
all_results, figsize=(6, 5), scale_axis=True, figure_name="scan_results.svg"
):
def plot_scan_plt(all_results, figsize=(6, 5), scale_axis=True, figure_name="scan_results.svg"):
"""Plot results of pygenstability with matplotlib."""
scales = _get_scales(all_results, scale_axis=scale_axis)
plt.figure(figsize=figsize)
Expand Down
44 changes: 12 additions & 32 deletions src/pygenstability/pygenstability.py
Original file line number Diff line number Diff line change
Expand Up @@ -120,9 +120,7 @@ def _get_constructor_data(constructor, scales, pool, tqdm_disable=False):

def _check_method(method): # pragma: no cover
if _NO_LEIDEN and _NO_LOUVAIN:
raise Exception(
"Without Louvain or Leiden solver, we cannot run PyGenStability"
)
raise Exception("Without Louvain or Leiden solver, we cannot run PyGenStability")

if method == "louvain" and _NO_LOUVAIN:
print("Louvain is not available, we fallback to leiden.")
Expand Down Expand Up @@ -255,9 +253,7 @@ def run(
communities = _process_runs(t, results, all_results)

if with_NVI:
_compute_NVI(
communities, all_results, pool, n_partitions=min(n_NVI, n_tries)
)
_compute_NVI(communities, all_results, pool, n_partitions=min(n_NVI, n_tries))

if with_all_tries:
all_results["all_tries"].append(results)
Expand All @@ -266,9 +262,7 @@ def run(

if with_postprocessing:
L.info("Apply postprocessing...")
_apply_postprocessing(
all_results, pool, constructor_data, tqdm_disable, method=method
)
_apply_postprocessing(all_results, pool, constructor_data, tqdm_disable, method=method)

if with_ttprime or with_optimal_scales:
L.info("Compute ttprimes...")
Expand All @@ -282,9 +276,7 @@ def run(
"window_size": max(2, int(0.1 * n_scale)),
"basin_radius": max(1, int(0.01 * n_scale)),
}
all_results = identify_optimal_scales(
all_results, **optimal_scales_kwargs
)
all_results = identify_optimal_scales(all_results, **optimal_scales_kwargs)

save_results(all_results, filename=result_file)

Expand Down Expand Up @@ -313,9 +305,7 @@ def _compute_NVI(communities, all_results, pool, n_partitions=10):
worker = partial(evaluate_NVI, partitions=selected_partitions)
index_pairs = [[i, j] for i in range(n_partitions) for j in range(n_partitions)]
chunksize = _get_chunksize(len(index_pairs), pool)
all_results["NVI"].append(
np.mean(list(pool.imap(worker, index_pairs, chunksize=chunksize)))
)
all_results["NVI"].append(np.mean(list(pool.imap(worker, index_pairs, chunksize=chunksize))))


def evaluate_NVI(index_pair, partitions):
Expand Down Expand Up @@ -360,9 +350,7 @@ def _to_indices(matrix, directed=False):


@_timing
def _optimise(
_, quality_indices, quality_values, null_model, global_shift, method="louvain"
):
def _optimise(_, quality_indices, quality_values, null_model, global_shift, method="louvain"):
"""Worker for generalized Markov Stability optimisation runs."""
if method == "louvain":
stability, community_id = generalized_louvain.run_louvain(
Expand Down Expand Up @@ -402,9 +390,7 @@ def _optimise(
return stability + global_shift, community_id


def _evaluate_quality(
partition_id, qualities_index, null_model, global_shift, method="louvain"
):
def _evaluate_quality(partition_id, qualities_index, null_model, global_shift, method="louvain"):
"""Worker for generalized Markov Stability optimisation runs."""
if method == "louvain":
quality = generalized_louvain.evaluate_quality(
Expand Down Expand Up @@ -461,18 +447,14 @@ def _compute_ttprime(all_results, pool):
chunksize = _get_chunksize(len(index_pairs), pool)
ttprime_list = pool.map(worker, index_pairs, chunksize=chunksize)

all_results["ttprime"] = np.zeros(
[len(all_results["scales"]), len(all_results["scales"])]
)
all_results["ttprime"] = np.zeros([len(all_results["scales"]), len(all_results["scales"])])
for i, ttp in enumerate(ttprime_list):
all_results["ttprime"][index_pairs[i][0], index_pairs[i][1]] = ttp
all_results["ttprime"] += all_results["ttprime"].T


@_timing
def _apply_postprocessing(
all_results, pool, constructors, tqdm_disable=False, method="louvain"
):
def _apply_postprocessing(all_results, pool, constructors, tqdm_disable=False, method="louvain"):
"""Apply postprocessing."""
all_results_raw = all_results.copy()

Expand All @@ -494,10 +476,8 @@ def _apply_postprocessing(
)
)

all_results["community_id"][i] = all_results_raw["community_id"][
all_results["community_id"][i] = all_results_raw["community_id"][best_quality_id]
all_results["stability"][i] = all_results_raw["stability"][best_quality_id]
all_results["number_of_communities"][i] = all_results_raw["number_of_communities"][
best_quality_id
]
all_results["stability"][i] = all_results_raw["stability"][best_quality_id]
all_results["number_of_communities"][i] = all_results_raw[
"number_of_communities"
][best_quality_id]
Binary file added tests/data/partitions/scale_1.pdf
Binary file not shown.
Loading

0 comments on commit 52fa5c0

Please sign in to comment.