Skip to content

Commit 51e18a2

Browse files
codegithubkacodegithubka
committed
decoupled plotting utilities from hpc script entirely for our sanity
1 parent 7fed52b commit 51e18a2

2 files changed

Lines changed: 787 additions & 295 deletions

File tree

‎scripts/pp_analysis.py‎

Lines changed: 1 addition & 295 deletions
Original file line numberDiff line numberDiff line change
@@ -571,7 +571,6 @@ def run_single_simulation_fss(
571571

572572
sample_counter += 1
573573

574-
# <-- FIX: Build result dict AFTER the loop, not inside it
575574
result = {
576575
"prey_birth": prey_birth,
577576
"prey_death": prey_death,
@@ -814,298 +813,6 @@ def run_debug_mode(cfg: Config, logger: logging.Logger):
814813
logger.info("Simulation complete.")
815814
input("Press Enter to exit...")
816815

817-
# =============================================================================
818-
# PLOTTING
819-
# =============================================================================
820-
821-
def generate_plots(cfg: Config, output_dir: Path, logger: logging.Logger):
822-
"""Generate all analysis plots from saved data."""
823-
import matplotlib.pyplot as plt
824-
from collections import defaultdict
825-
826-
plt.rcParams["figure.figsize"] = (14, 10)
827-
plt.rcParams["font.size"] = 11
828-
829-
prey_births = cfg.get_prey_births()
830-
prey_deaths = cfg.get_prey_deaths()
831-
n_pb, n_pd = len(prey_births), len(prey_deaths)
832-
extent = [prey_births[0], prey_births[-1], prey_deaths[0], prey_deaths[-1]]
833-
834-
# Load sweep results
835-
sweep_file = output_dir / "sweep_results.npz"
836-
if not sweep_file.exists():
837-
# Try JSON fallback
838-
sweep_file = output_dir / "sweep_results.json"
839-
if not sweep_file.exists():
840-
logger.error(f"Sweep results not found")
841-
return
842-
with open(sweep_file, "r") as f:
843-
results = json.load(f)
844-
else:
845-
results = load_sweep_binary(sweep_file)
846-
847-
logger.info(f"Loaded {len(results)} results")
848-
849-
# Initialize grids
850-
grids = {
851-
"prey_pop_no_evo": np.full((n_pd, n_pb), np.nan),
852-
"prey_pop_evo": np.full((n_pd, n_pb), np.nan),
853-
"pred_pop_no_evo": np.full((n_pd, n_pb), np.nan),
854-
"pred_pop_evo": np.full((n_pd, n_pb), np.nan),
855-
"survival_prey_no_evo": np.full((n_pd, n_pb), np.nan),
856-
"survival_prey_evo": np.full((n_pd, n_pb), np.nan),
857-
"tau_prey": np.full((n_pd, n_pb), np.nan),
858-
"evolved_prey_death": np.full((n_pd, n_pb), np.nan),
859-
"segregation_index": np.full((n_pd, n_pb), np.nan),
860-
"prey_clustering_index": np.full((n_pd, n_pb), np.nan),
861-
}
862-
863-
# Group by parameters
864-
grouped = defaultdict(list)
865-
for r in results:
866-
key = (round(r["prey_birth"], 4), round(r["prey_death"], 4), r["with_evolution"])
867-
grouped[key].append(r)
868-
869-
# Aggregate into grids
870-
for i, pd in enumerate(prey_deaths):
871-
for j, pb in enumerate(prey_births):
872-
pd_r, pb_r = round(pd, 4), round(pb, 4)
873-
874-
# No evolution
875-
no_evo = grouped.get((pb_r, pd_r, False), [])
876-
if no_evo:
877-
grids["prey_pop_no_evo"][i, j] = np.mean([r["prey_mean"] for r in no_evo])
878-
grids["pred_pop_no_evo"][i, j] = np.mean([r["pred_mean"] for r in no_evo])
879-
grids["survival_prey_no_evo"][i, j] = np.mean([r["prey_survived"] for r in no_evo]) * 100
880-
881-
taus = [r["prey_tau"] for r in no_evo if not np.isnan(r.get("prey_tau", np.nan))]
882-
if taus:
883-
grids["tau_prey"][i, j] = np.mean(taus)
884-
885-
seg = [r.get("segregation_index", np.nan) for r in no_evo]
886-
seg = [s for s in seg if not np.isnan(s)]
887-
if seg:
888-
grids["segregation_index"][i, j] = np.mean(seg)
889-
890-
clust = [r.get("prey_clustering_index", np.nan) for r in no_evo]
891-
clust = [c for c in clust if not np.isnan(c)]
892-
if clust:
893-
grids["prey_clustering_index"][i, j] = np.mean(clust)
894-
895-
# With evolution
896-
evo = grouped.get((pb_r, pd_r, True), [])
897-
if evo:
898-
grids["prey_pop_evo"][i, j] = np.mean([r["prey_mean"] for r in evo])
899-
grids["pred_pop_evo"][i, j] = np.mean([r["pred_mean"] for r in evo])
900-
grids["survival_prey_evo"][i, j] = np.mean([r["prey_survived"] for r in evo]) * 100
901-
902-
evolved = [r.get("evolved_prey_death_mean", np.nan) for r in evo]
903-
evolved = [e for e in evolved if not np.isnan(e)]
904-
if evolved:
905-
grids["evolved_prey_death"][i, j] = np.mean(evolved)
906-
907-
# Compute Hydra derivative
908-
dd = prey_deaths[1] - prey_deaths[0]
909-
dN_dd_no_evo = np.zeros_like(grids["prey_pop_no_evo"])
910-
dN_dd_evo = np.zeros_like(grids["prey_pop_evo"])
911-
912-
for j in range(n_pb):
913-
pop_smooth = gaussian_filter1d(grids["prey_pop_no_evo"][:, j], sigma=0.8)
914-
dN_dd_no_evo[:, j] = np.gradient(pop_smooth, dd)
915-
pop_smooth = gaussian_filter1d(grids["prey_pop_evo"][:, j], sigma=0.8)
916-
dN_dd_evo[:, j] = np.gradient(pop_smooth, dd)
917-
918-
# =========================================================================
919-
# PLOT 1: Phase Diagrams
920-
# =========================================================================
921-
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
922-
923-
ax = axes[0, 0]
924-
im = ax.imshow(grids["prey_pop_no_evo"], origin="lower", aspect="auto",
925-
extent=extent, cmap="YlGn")
926-
ax.contour(prey_births, prey_deaths, grids["survival_prey_no_evo"],
927-
levels=[50], colors="black", linewidths=2)
928-
plt.colorbar(im, ax=ax, label="Population")
929-
ax.set_xlabel("Prey Birth Rate")
930-
ax.set_ylabel("Prey Death Rate")
931-
ax.set_title("Prey Pop (No Evolution)")
932-
933-
ax = axes[0, 1]
934-
im = ax.imshow(grids["prey_pop_evo"], origin="lower", aspect="auto",
935-
extent=extent, cmap="YlGn")
936-
ax.contour(prey_births, prey_deaths, grids["survival_prey_evo"],
937-
levels=[50], colors="black", linewidths=2)
938-
plt.colorbar(im, ax=ax, label="Population")
939-
ax.set_xlabel("Prey Birth Rate")
940-
ax.set_ylabel("Prey Death Rate")
941-
ax.set_title("Prey Pop (With Evolution)")
942-
943-
ax = axes[0, 2]
944-
advantage = np.where(
945-
grids["prey_pop_no_evo"] > 10,
946-
(grids["prey_pop_evo"] - grids["prey_pop_no_evo"]) / grids["prey_pop_no_evo"] * 100,
947-
np.where(grids["prey_pop_evo"] > 10, 500, 0),
948-
)
949-
im = ax.imshow(np.clip(advantage, -50, 200), origin="lower", aspect="auto",
950-
extent=extent, cmap="RdYlGn", vmin=-50, vmax=200)
951-
plt.colorbar(im, ax=ax, label="Advantage (%)")
952-
ax.set_xlabel("Prey Birth Rate")
953-
ax.set_ylabel("Prey Death Rate")
954-
ax.set_title("Evolution Advantage")
955-
956-
ax = axes[1, 0]
957-
im = ax.imshow(grids["tau_prey"], origin="lower", aspect="auto",
958-
extent=extent, cmap="coolwarm", vmin=1.5, vmax=2.5)
959-
ax.contour(prey_births, prey_deaths, grids["tau_prey"],
960-
levels=[2.05], colors="green", linewidths=2)
961-
plt.colorbar(im, ax=ax, label="τ")
962-
ax.set_xlabel("Prey Birth Rate")
963-
ax.set_ylabel("Prey Death Rate")
964-
ax.set_title("Prey τ (Green: τ=2.05)")
965-
966-
ax = axes[1, 1]
967-
im = ax.imshow(grids["evolved_prey_death"], origin="lower", aspect="auto",
968-
extent=extent, cmap="viridis")
969-
plt.colorbar(im, ax=ax, label="Evolved d")
970-
ax.set_xlabel("Prey Birth Rate")
971-
ax.set_ylabel("Initial Prey Death Rate")
972-
ax.set_title("Evolved Prey Death Rate")
973-
974-
ax = axes[1, 2]
975-
im = ax.imshow(dN_dd_no_evo, origin="lower", aspect="auto",
976-
extent=extent, cmap="RdBu_r", vmin=-5000, vmax=5000)
977-
ax.contour(prey_births, prey_deaths, dN_dd_no_evo,
978-
levels=[0], colors="black", linewidths=2)
979-
plt.colorbar(im, ax=ax, label="dN/dd")
980-
ax.set_xlabel("Prey Birth Rate")
981-
ax.set_ylabel("Prey Death Rate")
982-
ax.set_title("HYDRA: dN/dd (Red: Prey ↑ with mortality)")
983-
984-
plt.tight_layout()
985-
plt.savefig(output_dir / "phase_diagrams.png", dpi=150, bbox_inches="tight")
986-
plt.close()
987-
logger.info("Saved phase_diagrams.png")
988-
989-
# =========================================================================
990-
# PLOT 2: Hydra Analysis
991-
# =========================================================================
992-
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
993-
994-
ax = axes[0]
995-
im = ax.imshow(dN_dd_no_evo, origin="lower", aspect="auto",
996-
extent=extent, cmap="RdBu_r", vmin=-5000, vmax=5000)
997-
ax.contour(prey_births, prey_deaths, dN_dd_no_evo,
998-
levels=[0], colors="black", linewidths=2)
999-
plt.colorbar(im, ax=ax, label="dN/dd")
1000-
ax.set_xlabel("Prey Birth Rate")
1001-
ax.set_ylabel("Prey Death Rate")
1002-
ax.set_title("Hydra (No Evolution)")
1003-
1004-
ax = axes[1]
1005-
im = ax.imshow(dN_dd_evo, origin="lower", aspect="auto",
1006-
extent=extent, cmap="RdBu_r", vmin=-5000, vmax=5000)
1007-
ax.contour(prey_births, prey_deaths, dN_dd_evo,
1008-
levels=[0], colors="black", linewidths=2)
1009-
plt.colorbar(im, ax=ax, label="dN/dd")
1010-
ax.set_xlabel("Prey Birth Rate")
1011-
ax.set_ylabel("Prey Death Rate")
1012-
ax.set_title("Hydra (With Evolution)")
1013-
1014-
ax = axes[2]
1015-
mid_pb_idx = n_pb // 2
1016-
target_pb = prey_births[mid_pb_idx]
1017-
no_evo_slice = grids["prey_pop_no_evo"][:, mid_pb_idx]
1018-
evo_slice = grids["prey_pop_evo"][:, mid_pb_idx]
1019-
1020-
ax.plot(prey_deaths, no_evo_slice, 'b-o',
1021-
label=f'No Evo (pb={target_pb:.2f})', markersize=4)
1022-
ax.plot(prey_deaths, evo_slice, 'g-s',
1023-
label=f'With Evo (pb={target_pb:.2f})', markersize=4)
1024-
1025-
ax.set_xlabel("Prey Death Rate")
1026-
ax.set_ylabel("Prey Population")
1027-
ax.set_title(f"Prey Pop vs Death Rate (pb={target_pb:.2f})") # Dynamic title
1028-
ax.legend()
1029-
ax.grid(True, alpha=0.3)
1030-
1031-
plt.tight_layout()
1032-
plt.savefig(output_dir / "hydra_analysis.png", dpi=150, bbox_inches="tight")
1033-
plt.close()
1034-
logger.info("Saved hydra_analysis.png")
1035-
1036-
# =========================================================================
1037-
# PLOT 3: PCF Analysis
1038-
# =========================================================================
1039-
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
1040-
1041-
ax = axes[0]
1042-
im = ax.imshow(grids["segregation_index"], origin="lower", aspect="auto",
1043-
extent=extent, cmap="RdBu", vmin=0.5, vmax=1.5)
1044-
ax.contour(prey_births, prey_deaths, grids["segregation_index"],
1045-
levels=[1.0], colors="black", linewidths=2)
1046-
plt.colorbar(im, ax=ax, label="C_cr")
1047-
ax.set_xlabel("Prey Birth Rate")
1048-
ax.set_ylabel("Prey Death Rate")
1049-
ax.set_title("Segregation Index")
1050-
1051-
ax = axes[1]
1052-
im = ax.imshow(grids["prey_clustering_index"], origin="lower", aspect="auto",
1053-
extent=extent, cmap="Greens", vmin=1.0, vmax=3.0)
1054-
plt.colorbar(im, ax=ax, label="C_rr")
1055-
ax.set_xlabel("Prey Birth Rate")
1056-
ax.set_ylabel("Prey Death Rate")
1057-
ax.set_title("Prey Clustering")
1058-
1059-
ax = axes[2]
1060-
im = ax.imshow(grids["segregation_index"], origin="lower", aspect="auto",
1061-
extent=extent, cmap="RdBu", vmin=0.5, vmax=1.5, alpha=0.8)
1062-
ax.contour(prey_births, prey_deaths, dN_dd_no_evo,
1063-
levels=[0], colors="lime", linewidths=3)
1064-
ax.contour(prey_births, prey_deaths, grids["survival_prey_no_evo"],
1065-
levels=[50], colors="black", linewidths=2, linestyles='--')
1066-
plt.colorbar(im, ax=ax, label="C_cr")
1067-
ax.set_xlabel("Prey Birth Rate")
1068-
ax.set_ylabel("Prey Death Rate")
1069-
ax.set_title("Segregation + Hydra Boundary")
1070-
1071-
plt.tight_layout()
1072-
plt.savefig(output_dir / "pcf_analysis.png", dpi=150, bbox_inches="tight")
1073-
plt.close()
1074-
logger.info("Saved pcf_analysis.png")
1075-
1076-
# =========================================================================
1077-
# Summary Statistics
1078-
# =========================================================================
1079-
summary = {
1080-
"coexistence_no_evo": int(np.sum((grids["survival_prey_no_evo"] > 80))),
1081-
"hydra_region_size": int(np.sum((dN_dd_no_evo > 0) & (grids["prey_pop_no_evo"] > 50))),
1082-
"max_hydra_strength": float(np.nanmax(dN_dd_no_evo)),
1083-
"hydra_region_size_evo": int(np.sum((dN_dd_evo > 0) & (grids["prey_pop_evo"] > 50))),
1084-
"mean_segregation_index": float(np.nanmean(grids["segregation_index"])),
1085-
"mean_prey_clustering": float(np.nanmean(grids["prey_clustering_index"])),
1086-
}
1087-
1088-
# Find critical point
1089-
dist_crit = np.abs(grids["tau_prey"] - 2.05)
1090-
if not np.all(np.isnan(dist_crit)):
1091-
min_idx = np.unravel_index(np.nanargmin(dist_crit), dist_crit.shape)
1092-
summary["critical_prey_birth"] = float(prey_births[min_idx[1]])
1093-
summary["critical_prey_death"] = float(prey_deaths[min_idx[0]])
1094-
summary["critical_tau_prey"] = float(grids["tau_prey"][min_idx])
1095-
1096-
with open(output_dir / "summary.json", "w") as f:
1097-
json.dump(summary, f, indent=2)
1098-
1099-
logger.info("=" * 60)
1100-
logger.info("ANALYSIS SUMMARY")
1101-
logger.info("=" * 60)
1102-
logger.info(f"Hydra region: {summary['hydra_region_size']} combinations")
1103-
logger.info(f"Max Hydra strength: {summary['max_hydra_strength']:.1f}")
1104-
logger.info(f"Mean segregation index: {summary['mean_segregation_index']:.3f}")
1105-
if "critical_prey_birth" in summary:
1106-
logger.info(f"Critical point: pb={summary['critical_prey_birth']:.3f}, pd={summary['critical_prey_death']:.3f}")
1107-
1108-
1109816
# =============================================================================
1110817
# MAIN
1111818
# =============================================================================
@@ -1182,8 +889,7 @@ def main():
1182889
run_fss(cfg, output_dir, logger)
1183890

1184891
if args.mode in ["full", "plot"]:
1185-
generate_plots(cfg, output_dir, logger)
1186-
892+
pass #NOTE: Decoupled plots into a separate script for clarity
1187893
elapsed = time.time() - start_time
1188894
logger.info(f"Total runtime: {elapsed/3600:.2f} hours")
1189895
logger.info("Done!")

0 commit comments

Comments
 (0)