Skip to content

Commit

Permalink
renamed block_detection_curve to block_nvi
Browse files Browse the repository at this point in the history
  • Loading branch information
d-schindler committed Nov 8, 2024
1 parent fe20adc commit d3b728c
Show file tree
Hide file tree
Showing 5 changed files with 57 additions and 25 deletions.
29 changes: 19 additions & 10 deletions src/pygenstability/optimal_scales.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,14 +33,21 @@ def _pool2d_nvi(A, kernel_size, stride, padding=0):
)
shape_w = (output_shape[0], output_shape[1], kernel_size, kernel_size)
# pylint: disable=unsubscriptable-object
strides_w = (stride * A.strides[0], stride * A.strides[1], A.strides[0], A.strides[1])
strides_w = (
stride * A.strides[0],
stride * A.strides[1],
A.strides[0],
A.strides[1],
)
A_w = as_strided(A, shape_w, strides_w)

# Return the result of pooling
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 All @@ -57,7 +64,7 @@ def identify_optimal_scales(results, kernel_size=3, window_size=3, max_nvi=1, ba
basin_radius (int): radius of basin around local minima of the pooled diagonal
Returns:
result dictionary with two new keys: 'selected_partitions' and 'block_detection_curve'
result dictionary with two new keys: 'selected_partitions' and 'block_nvi'
References:
.. [1] D. J. Schindler, J. Clarke, and M. Barahona, 'Multiscale Mobility Patterns and
Expand All @@ -74,22 +81,24 @@ def identify_optimal_scales(results, kernel_size=3, window_size=3, max_nvi=1, ba
diagonal = np.diag(nvi_tt_pooled)[: len(nvi_t)]

# smooth diagonal with moving window
block_detection_curve = np.roll(
np.asarray(pd.Series(diagonal).rolling(window=window_size, win_type="triang").mean()),
block_nvi = np.roll(
np.asarray(
pd.Series(diagonal).rolling(window=window_size, win_type="triang").mean()
),
-int(window_size / 2),
)
results["block_detection_curve"] = block_detection_curve
results["block_nvi"] = block_nvi

# find local minima on diagonal of pooled NVI(s,s')
basin_centers, _ = find_peaks(-block_detection_curve, height=-max_nvi)
basin_centers, _ = find_peaks(-block_nvi, height=-max_nvi)

# add robust scales located in large 0 margins
not_nan_ind = np.argwhere(~np.isnan(block_detection_curve)).flatten()
not_nan_ind = np.argwhere(~np.isnan(block_nvi)).flatten()

if (
np.count_nonzero(
np.around(
block_detection_curve[not_nan_ind[0] : not_nan_ind[0] + 2 * basin_radius + 1], 5
block_nvi[not_nan_ind[0] : not_nan_ind[0] + 2 * basin_radius + 1], 5
)
)
== 0
Expand All @@ -99,7 +108,7 @@ def identify_optimal_scales(results, kernel_size=3, window_size=3, max_nvi=1, ba
if (
np.count_nonzero(
np.around(
block_detection_curve[not_nan_ind[-1] - 2 * basin_radius : not_nan_ind[-1] + 1], 5
block_nvi[not_nan_ind[-1] - 2 * basin_radius : not_nan_ind[-1] + 1], 5
)
)
== 0
Expand Down
47 changes: 35 additions & 12 deletions src/pygenstability/plotting.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,9 @@
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 @@ -48,7 +50,9 @@ 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 @@ -263,15 +267,24 @@ def plot_optimal_partitions(

for optimal_scale_id in selected_scales:
plot_single_partition(
graph, all_results, optimal_scale_id, edge_color=edge_color, edge_width=edge_width
graph,
all_results,
optimal_scale_id,
edge_color=edge_color,
edge_width=edge_width,
)
plt.savefig(f"{folder}/scale_{optimal_scale_id}{ext}", bbox_inches="tight")
if show: # pragma: no cover
plt.show()


def plot_communities(
graph, all_results, folder="communities", edge_color="0.5", edge_width=0.5, ext=".pdf"
graph,
all_results,
folder="communities",
edge_color="0.5",
edge_width=0.5,
ext=".pdf",
):
"""Plot the community structures at each scale in a folder.
Expand All @@ -293,12 +306,16 @@ 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 @@ -327,7 +344,9 @@ def plot_communities_matrix(graph, all_results, folder="communities_matrix", ext
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 @@ -341,15 +360,19 @@ 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 @@ -394,15 +417,15 @@ def _plot_optimal_scales(all_results, ax, scales, ax1, ax2):
"""Plot stability."""
ax.plot(
scales,
all_results["block_detection_curve"],
all_results["block_nvi"],
"-",
lw=2.0,
c="C4",
label="Block NVI",
)
ax.plot(
scales[all_results["selected_partitions"]],
all_results["block_detection_curve"][all_results["selected_partitions"]],
all_results["block_nvi"][all_results["selected_partitions"]],
"o",
lw=2.0,
c="C4",
Expand Down Expand Up @@ -456,7 +479,7 @@ def plot_scan_plt(all_results, scale_axis=True, figure_name="scan_results.svg"):
_plot_number_comm(all_results, ax=ax3, scales=scales)
axes.append(ax3)

if "block_detection_curve" in all_results:
if "block_nvi" in all_results:
ax4 = plt.subplot(gs[2, 0])
_plot_optimal_scales(all_results, ax=ax4, scales=scales, ax1=ax1, ax2=ax2)
axes.append(ax4)
Expand Down
2 changes: 1 addition & 1 deletion src/pygenstability/pygenstability.py
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ def run(
optimisation evaluations at each scale (included if with_all_tries==True)
- 'NVI': NVI(t) at each scale
- 'ttprime': NVI(t,tprime) matrix
- 'block_detection_curve': block detection curve (included if with_optimal_scales==True)
- 'block_nvi': block NVI curve (included if with_optimal_scales==True)
- 'selected_partitions': selected partitions (included if with_optimal_scales==True)
"""
Expand Down
2 changes: 1 addition & 1 deletion tests/data/test_dataclustering_default.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -2400,7 +2400,7 @@ labels:
- 31.0
- 31.0
results:
block_detection_curve:
block_nvi:
- .nan
- 0.0680975542371894
- 0.05125455197832432
Expand Down
2 changes: 1 addition & 1 deletion tests/data/test_dataclustering_knn.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -1800,7 +1800,7 @@ labels:
- 25.0
- 25.0
results:
block_detection_curve:
block_nvi:
- .nan
- 0.10668153950024195
- 0.08542249390900095
Expand Down

0 comments on commit d3b728c

Please sign in to comment.