diff --git a/examples/bioma/fig4.py b/examples/bioma/fig4.py index 70587f8..2d96d7d 100644 --- a/examples/bioma/fig4.py +++ b/examples/bioma/fig4.py @@ -1,7 +1,7 @@ from concurrent.futures import ProcessPoolExecutor from pathlib import Path -from plotly.subplots import make_subplots +import matplotlib.pyplot as plt from problems import ackley4, griewank, schwefel2_26 from utils import IMAGES_DIR, FunctionConfig, build_cmlon @@ -31,56 +31,94 @@ N_VAR = 5 +def render_3d_cmlons(func_names, cmlons): + standard_camera = dict( + up=dict(x=0, y=0, z=1), center=dict(x=0, y=0, z=0), eye=dict(x=1.55, y=1.55, z=0.4) + ) + + axis_config = dict( + visible=True, + showgrid=True, + gridcolor="lightgray", + showline=True, + linecolor="black", + showbackground=True, + backgroundcolor="rgb(250, 250, 250)", + zeroline=True, + zerolinecolor="gray", + showticklabels=True, + ) + + plot_paths = [] + + for func_name in func_names: + viz = LONVisualizer() + fig = viz.plot_3d(cmlons[func_name]) + + fig.update_layout( + scene=dict( + xaxis=dict(**axis_config, title="X"), + yaxis=dict(**axis_config, title="Y"), + zaxis=dict(**axis_config, title="Fitness"), + camera=dict(**standard_camera), + aspectmode="cube", + ), + title=dict( + text=f"{func_name}", + x=0.5, + y=0.95, + xanchor="center", + font=dict(size=20), + ), + showlegend=False, + width=600, + height=600, + margin=dict(l=20, r=20, t=60, b=20), + paper_bgcolor="rgb(255, 255, 255)", + plot_bgcolor="rgb(255, 255, 255)", + ) + + path = Path(IMAGES_DIR) / f"fig4_{func_name.replace(' ', '_')}.png" + fig.write_image(path, scale=2) + plot_paths.append((func_name, path)) + print(f"Saved 3D panel to {path}") + + return plot_paths + + +def merge_3d_plots(plot_paths): + merged_fig, axes = plt.subplots(1, len(plot_paths), figsize=(6 * len(plot_paths), 6)) + if len(plot_paths) == 1: + axes = [axes] + + for ax, (_, image_path) in zip(axes, plot_paths): + img = plt.imread(image_path) + ax.imshow(img) + ax.axis("off") + + merged_fig.tight_layout() + + final_path = Path(IMAGES_DIR) / "fig4.png" + merged_fig.savefig(final_path, dpi=200, bbox_inches="tight") + plt.close(merged_fig) + print(f"Successfully saved merged figure to {final_path}") + + def main() -> None: Path(IMAGES_DIR).mkdir(parents=True, exist_ok=True) - - viz = LONVisualizer() func_names = list(FUNCTIONS.keys()) - # Build CMLONs in parallel, then create 3D figures with ProcessPoolExecutor() as executor: futures = { func_name: executor.submit(build_cmlon, FUNCTIONS[func_name], N_VAR) for func_name in func_names } - cmlons = {name: fut.result() for name, fut in futures.items()} - figures = [] - for func_name in func_names: - fig = viz.plot_3d(cmlons[func_name]) - figures.append(fig) - - # Combine into a single 1x3 subplot figure - combined = make_subplots( - rows=1, - cols=3, - specs=[[{"type": "scene"}, {"type": "scene"}, {"type": "scene"}]], - subplot_titles=[ - f"(a) {func_names[0]}", - f"(b) {func_names[1]}", - f"(c) {func_names[2]}", - ], - horizontal_spacing=0.02, - ) + cmlons = {name: fut.result() for name, fut in futures.items()} - for idx, fig in enumerate(figures, start=1): - scene_name = f"scene{idx}" if idx > 1 else "scene" - for trace in fig.data: - trace.scene = scene_name - combined.add_trace(trace, row=1, col=idx) - # Copy camera / axis settings from original figure - if "scene" in fig.layout: - combined.layout[scene_name].update(fig.layout.scene) - - combined.update_layout( - showlegend=False, - width=1800, - height=600, - margin=dict(l=0, r=0, t=60, b=0), - ) + plot_paths = render_3d_cmlons(func_names, cmlons) - combined.write_image(f"{IMAGES_DIR}/fig4.png", scale=2) - combined.write_html(f"{IMAGES_DIR}/fig4.html") + merge_3d_plots(plot_paths) if __name__ == "__main__": diff --git a/examples/bioma/utils.py b/examples/bioma/utils.py index d367008..a6a44c3 100644 --- a/examples/bioma/utils.py +++ b/examples/bioma/utils.py @@ -230,7 +230,7 @@ def save_metrics_figure( fontsize=10, ) - fig.tight_layout(rect=[0, 0.06, 1, 1]) + fig.tight_layout(rect=(0, 0.06, 1, 1)) fig.savefig(str(output_path), dpi=150, bbox_inches="tight", facecolor="white") plt.close(fig) print(f"Saved {output_path}") diff --git a/examples/number_partitioning/npp_paths.py b/examples/number_partitioning/npp_paths.py new file mode 100644 index 0000000..349d1e2 --- /dev/null +++ b/examples/number_partitioning/npp_paths.py @@ -0,0 +1 @@ +IMAGES_DIR = "images/number_partitioning" diff --git a/examples/number_partitioning/plot_lons.py b/examples/number_partitioning/plot_lons.py new file mode 100644 index 0000000..4463d46 --- /dev/null +++ b/examples/number_partitioning/plot_lons.py @@ -0,0 +1,116 @@ +from pathlib import Path + +import matplotlib.pyplot as plt +from npp_paths import IMAGES_DIR + +from lonkit import ( + CMLON, + ILSSampler, + ILSSamplerConfig, + LONConfig, + LONVisualizer, + NumberPartitioning, +) + +N = 20 +INSTANCE_SEED = 1 +N_RUNS = 100 +N_ITER = 500 +RANDOM_SEED = 42 + +K_VALUES = [0.3, 0.7, 0.95] + + +def render_3d_lons(cmlon_by_k: dict[float, CMLON], output_dir: Path = Path(IMAGES_DIR)) -> None: + output_dir.mkdir(parents=True, exist_ok=True) + + standard_camera = dict( + up=dict(x=0, y=0, z=1), + center=dict(x=0, y=0, z=0), + eye=dict(x=1.55, y=1.55, z=0.4), + ) + + axis_config = dict( + visible=True, + showgrid=True, + gridcolor="lightgray", + showline=True, + linecolor="black", + showbackground=True, + backgroundcolor="rgb(250, 250, 250)", + zeroline=True, + zerolinecolor="gray", + showticklabels=True, + ) + + for k, cmlon in cmlon_by_k.items(): + vis = LONVisualizer(min_edge_width=0.5, max_edge_width=1, min_node_size=2.5, arrow_size=0.1) + fig = vis.plot_3d(cmlon) + fig.update_layout( + scene=dict( + xaxis=dict(**axis_config, title="X"), + yaxis=dict(**axis_config, title="Y"), + zaxis=dict(**axis_config, title="Fitness"), + camera=dict(**standard_camera), + aspectmode="cube", + ), + showlegend=False, + width=900, + height=700, + margin=dict(l=20, r=20, t=40, b=20), + ) + + fig.write_image(output_dir / f"NPP_{k}_3d.png", scale=2) + + +def render_merged_lon_grid(k_values: list[float], output_dir: Path = Path(IMAGES_DIR)) -> None: + fig, axes = plt.subplots(2, len(k_values), figsize=(5 * len(k_values), 10)) + if len(k_values) == 1: + axes = [[axes[0]], [axes[1]]] + + for idx, k in enumerate(k_values): + img_2d = plt.imread(output_dir / f"NPP_{k}_2d.png") + img_3d = plt.imread(output_dir / f"NPP_{k}_3d.png") + + ax_top = axes[0][idx] + ax_top.imshow(img_2d) + ax_top.set_title(f"2D CMLON (k={k})", fontsize=12) + ax_top.axis("off") + + ax_bottom = axes[1][idx] + ax_bottom.imshow(img_3d) + ax_bottom.set_title(f"3D CMLON (k={k})", fontsize=12) + ax_bottom.axis("off") + + fig.suptitle("Number Partitioning CMLON Views", fontsize=16, y=0.98) + fig.tight_layout(rect=(0, 0, 1, 0.96)) + fig.savefig(output_dir / "NPP_merged_cmlon_views.png", dpi=200) + plt.close(fig) + + +def main(): + Path(IMAGES_DIR).mkdir(parents=True, exist_ok=True) + + sampler_config = ILSSamplerConfig(n_runs=N_RUNS, n_iter_no_change=N_ITER, seed=RANDOM_SEED) + + lon_config = LONConfig(eq_atol=1e-8) + cmlon_by_k = {} + + for k in K_VALUES: + problem = NumberPartitioning(n=N, k=k, instance_seed=INSTANCE_SEED) + sampler = ILSSampler(sampler_config) + result = sampler.sample(problem) + + lon = sampler.sample_to_lon(result, lon_config) + cmlon = lon.to_cmlon() + cmlon_by_k[k] = cmlon + + vis = LONVisualizer(0.5, 1, arrow_size=0.1) + vis.plot_2d(cmlon, f"{IMAGES_DIR}/NPP_{k}_2d.png") + + render_3d_lons(cmlon_by_k) + render_merged_lon_grid(K_VALUES) + + +if __name__ == "__main__": + main() diff --git a/examples/number_partitioning/plot_metrics.py b/examples/number_partitioning/plot_metrics.py new file mode 100644 index 0000000..f094630 --- /dev/null +++ b/examples/number_partitioning/plot_metrics.py @@ -0,0 +1,285 @@ +""" +Number Partitioning Problem: Phase Transition Sweep +=================================================== +The example script sweeps the phase-transition parameter k across its full +range [0, 1] and plots LON / CMLON network metrics, ILS success, and the +transition-rate derivative. + +The transition centre k* is located at the steepest descent of the ILS success +rate - i.e. the k value where the smoothed curve of -d(success)/dk is maximal. + +References +---------- +Mertens, S. (1998). Phase transition in the number partitioning problem. + Physical Review Letters, 81(20), 4281-4284. +Borgs, C., Chayes, J., Mertens, S., & Nair, C. (2001). Phase transition and + finite-size scaling for the integer partitioning problem. + Random Structures & Algorithms, 19(3-4), 261-294. +Gent, I. P., & Walsh, T. (1998). Analysis of heuristics for number partitioning. + Computational Intelligence, 14(2), 430-451. + +Outputs one figure: ``NPP_phase_transition_sweep.png`` +""" + +from pathlib import Path + +import matplotlib.patches as mpatches +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.lines import Line2D +from npp_paths import IMAGES_DIR +from scipy.ndimage import uniform_filter1d + +from lonkit import ILSSampler, ILSSamplerConfig, LONConfig, NumberPartitioning + +SMOOTH_WIDTH = 3 # uniform moving-average window (points); odd integer + + +def detect_via_derivative(ks, success, smooth_width=SMOOTH_WIDTH): + """ + Locate k* from the ILS success curve: k at the maximum of smoothed -d(success)/dk + (positive-clipped). Returns the derivative series used for the bottom panel. + """ + smoothed = uniform_filter1d(success.astype(float), size=smooth_width, mode="nearest") + deriv = -np.gradient(smoothed, ks) + deriv = np.clip(deriv, 0, None) + + peak_idx = int(np.argmax(deriv)) + k_centre = float(ks[peak_idx]) + + return k_centre, deriv + + +def plot_npp_metrics( + ks, + n_optima, + n_funnels, + global_strength, + global_funnel_prop, + ils_success, + k_centre, + deriv_signal, +): + C_LINE = "#1a1a2e" + C_GLOBAL = "#2563eb" + C_SINK = "#dc2626" + C_SUCCESS = "#16a34a" + C_DERIV = "#7c3aed" + C_BAND = "#fef08a" + C_REF = "#6b7280" + + # Mark the "hard" region in the figure + + _X_K_RIGHT = float(ks.max()) + 0.02 + + _YLABEL_KW = { + "rotation": 90, + "ha": "center", + "va": "center", + "labelpad": 20, + "fontsize": 8, + } + + SUCCESS_YLABEL = "ILS success rate\nvs. best fitness in\nthe sampled network" + + PANELS = [ + ("Global to local\n funnel proportion", global_funnel_prop, C_GLOBAL, "^"), + ("Number of CMLON\nlocal optima", n_optima, C_LINE, "o"), + ("Number of CMLON\nfunnels", n_funnels, C_SINK, "s"), + (SUCCESS_YLABEL, ils_success, C_SUCCESS, "o"), + ("Global CMLON\nstrength", global_strength, C_GLOBAL, "D"), + ] + + def _draw_background(ax): + ax.axvspan(k_centre, _X_K_RIGHT, color=C_BAND, alpha=0.55, zorder=0) + ax.axvline(k_centre, color=C_REF, linewidth=0.9, linestyle="--", zorder=1) + + def _plot_metric_ax(ax, ylabel, series, colour, marker): + _draw_background(ax) + ax.plot( + ks, + series, + color=colour, + marker=marker, + markersize=4, + linewidth=1.4, + markeredgewidth=0.5, + markeredgecolor="white", + zorder=3, + ) + ax.set_ylabel(ylabel, **_YLABEL_KW) + ax.tick_params(axis="both", labelsize=9) + ax.spines[["top", "right"]].set_visible(False) + ax.grid(axis="y", linewidth=0.4, alpha=0.4) + if series.min() >= 0: + rng = series.max() - series.min() + ax.set_ylim(bottom=-0.05 * rng if rng > 0 else -0.1) + + def _plot_derivative_ax(ax): + _draw_background(ax) + ax.plot( + ks, + deriv_signal, + color=C_DERIV, + linewidth=1.5, + zorder=3, + label=r"$-\,d(\mathrm{success})/dk$ (smoothed)", + ) + ax.set_ylabel("Transition\nrate", **_YLABEL_KW) + ax.set_xlabel("Phase-transition parameter k", fontsize=10.5, labelpad=8) + ax.tick_params(axis="both", labelsize=9) + ax.spines[["top", "right"]].set_visible(False) + ax.legend( + fontsize=7.5, + loc="upper center", + bbox_to_anchor=(0.5, -0.38), + ncol=1, + framealpha=0.85, + edgecolor="#d1d5db", + ) + ax.set_xlim(ks.min() - 0.02, ks.max() + 0.02) + + # Legend + + legend_handles = [ + mpatches.Patch( + color=C_BAND, + alpha=0.55, + label=(f"hard region k \u2208 [{k_centre:.2f},\u202f{ks.max():.2f}] "), + ), + Line2D( + [0], + [0], + color=C_REF, + linestyle="--", + linewidth=0.9, + label=(f"transition centre k*\u202f=\u202f{k_centre:.2f} "), + ), + ] + + # Grid arrengment + + output_path_grid = Path(IMAGES_DIR) / "NPP_phase_transition_sweep.png" + fig_g = plt.figure(figsize=(14, 8)) + gs = fig_g.add_gridspec( + 3, + 3, + hspace=0.38, + wspace=0.36, + height_ratios=[1, 1, 0.36], + ) + + ax_g00 = fig_g.add_subplot(gs[0, 0]) + _plot_metric_ax(ax_g00, *PANELS[0]) + ax_g01 = fig_g.add_subplot(gs[0, 1], sharex=ax_g00) + _plot_metric_ax(ax_g01, *PANELS[1]) + ax_g02 = fig_g.add_subplot(gs[0, 2], sharex=ax_g00) + _plot_metric_ax(ax_g02, *PANELS[2]) + + ax_g10 = fig_g.add_subplot(gs[1, 0], sharex=ax_g00) + _plot_metric_ax(ax_g10, *PANELS[3]) + ax_g11 = fig_g.add_subplot(gs[1, 1], sharex=ax_g00) + _plot_metric_ax(ax_g11, *PANELS[4]) + ax_g12 = fig_g.add_subplot(gs[1, 2], sharex=ax_g00) + _plot_derivative_ax(ax_g12) + + ax_g_leg = fig_g.add_subplot(gs[2, :]) + ax_g_leg.set_axis_off() + ax_g_leg.legend( + handles=legend_handles, + fontsize=10, + loc="center", + ncol=1, + framealpha=0.95, + edgecolor="#d1d5db", + borderpad=0.45, + labelspacing=0.65, + ) + + for ax in (ax_g00, ax_g01, ax_g02, ax_g10, ax_g11): + ax.tick_params(axis="x", labelbottom=False) + + # Titles and saving + + fig_g.suptitle( + f"NPP phase transition - LON metric analysis (N={N}, {N_RUNS} ILS runs per k)", + fontsize=11, + fontweight="bold", + y=0.98, + ) + fig_g.subplots_adjust(left=0.07, right=0.98, top=0.90, bottom=0.11) + fig_g.savefig(output_path_grid, dpi=300, bbox_inches="tight", pad_inches=0.12) + plt.close(fig_g) + print(f"Saved \u2192 {output_path_grid}") + + +N = 20 +INSTANCE_SEED = 1 +N_RUNS = 100 +N_ITER = 500 +RANDOM_SEED = 42 +EQ_ATOL = 1e-8 + +K_VALUES = np.linspace(0.1, 1.0, 50) + + +def main(): + Path(IMAGES_DIR).mkdir(parents=True, exist_ok=True) + + print(f"Sweeping k across {len(K_VALUES)} values (n={N}, n_runs={N_RUNS})\n") + + records = [] + sampler_config = ILSSamplerConfig(n_runs=N_RUNS, n_iter_no_change=N_ITER, seed=RANDOM_SEED) + lon_config = LONConfig(eq_atol=EQ_ATOL) + + for k in K_VALUES: + problem = NumberPartitioning(n=N, k=k, instance_seed=INSTANCE_SEED) + sampler = ILSSampler(sampler_config) + result = sampler.sample(problem) + + lon = sampler.sample_to_lon(result, lon_config) + cmlon = lon.to_cmlon() + + m = cmlon.compute_metrics() + + records.append( + { + "k": k, + "n_optima": m["n_optima"], + "n_funnels": m["n_funnels"], + "global_strength": m["global_strength"], + "global_funnel_prop": m["global_funnel_proportion"], + "ils_success": m["success"], + } + ) + print( + f" k={k:.3f} optima={m['n_optima']:>3} funnels={m['n_funnels']:>2} " + f"success={m['success']:.0%}" + ) + + ks = np.array([r["k"] for r in records]) + n_optima = np.array([r["n_optima"] for r in records]) + n_funnels = np.array([r["n_funnels"] for r in records]) + global_strength = np.array([r["global_strength"] for r in records]) + global_funnel_prop = np.array([r["global_funnel_prop"] for r in records]) + ils_success = np.array([r["ils_success"] for r in records]) + + k_centre, deriv_signal = detect_via_derivative(ks, ils_success) + + print("\nHard-region detection results") + print(f" Transition centre k*: {k_centre:.3f}") + + plot_npp_metrics( + ks, + n_optima, + n_funnels, + global_strength, + global_funnel_prop, + ils_success, + k_centre, + deriv_signal, + ) + + +if __name__ == "__main__": + main() diff --git a/images/fig3.png b/images/fig3.png index d339323..433752f 100644 Binary files a/images/fig3.png and b/images/fig3.png differ diff --git a/images/fig4.png b/images/fig4.png index e009fe0..d6d12ef 100644 Binary files a/images/fig4.png and b/images/fig4.png differ diff --git a/images/fig4_Ackley_4.png b/images/fig4_Ackley_4.png new file mode 100644 index 0000000..cdb1ccc Binary files /dev/null and b/images/fig4_Ackley_4.png differ diff --git a/images/fig4_Griewank.png b/images/fig4_Griewank.png new file mode 100644 index 0000000..5d1915e Binary files /dev/null and b/images/fig4_Griewank.png differ diff --git a/images/fig4_Schwefel_2.26.png b/images/fig4_Schwefel_2.26.png new file mode 100644 index 0000000..0b90843 Binary files /dev/null and b/images/fig4_Schwefel_2.26.png differ diff --git a/images/number_partitioning/NPP_0.3_2d.png b/images/number_partitioning/NPP_0.3_2d.png new file mode 100644 index 0000000..4da200a Binary files /dev/null and b/images/number_partitioning/NPP_0.3_2d.png differ diff --git a/images/number_partitioning/NPP_0.3_3d.png b/images/number_partitioning/NPP_0.3_3d.png new file mode 100644 index 0000000..9bd5d6d Binary files /dev/null and b/images/number_partitioning/NPP_0.3_3d.png differ diff --git a/images/number_partitioning/NPP_0.7_2d.png b/images/number_partitioning/NPP_0.7_2d.png new file mode 100644 index 0000000..a86bc6d Binary files /dev/null and b/images/number_partitioning/NPP_0.7_2d.png differ diff --git a/images/number_partitioning/NPP_0.7_3d.png b/images/number_partitioning/NPP_0.7_3d.png new file mode 100644 index 0000000..55ad1bd Binary files /dev/null and b/images/number_partitioning/NPP_0.7_3d.png differ diff --git a/images/number_partitioning/NPP_0.95_2d.png b/images/number_partitioning/NPP_0.95_2d.png new file mode 100644 index 0000000..5cb6f15 Binary files /dev/null and b/images/number_partitioning/NPP_0.95_2d.png differ diff --git a/images/number_partitioning/NPP_0.95_3d.png b/images/number_partitioning/NPP_0.95_3d.png new file mode 100644 index 0000000..55c22e1 Binary files /dev/null and b/images/number_partitioning/NPP_0.95_3d.png differ diff --git a/images/number_partitioning/NPP_merged_cmlon_views.png b/images/number_partitioning/NPP_merged_cmlon_views.png new file mode 100644 index 0000000..6055cf7 Binary files /dev/null and b/images/number_partitioning/NPP_merged_cmlon_views.png differ diff --git a/images/number_partitioning/NPP_phase_transition_sweep.png b/images/number_partitioning/NPP_phase_transition_sweep.png new file mode 100644 index 0000000..a06c49b Binary files /dev/null and b/images/number_partitioning/NPP_phase_transition_sweep.png differ