Skip to content

Commit f2075e8

Browse files
Abou1308Abou1308
authored
Merge pull request #12 from codegithubka/Sary
Bifurcation Analysis and Warmup Study
2 parents c283ce2 + 164d66a commit f2075e8

3 files changed

Lines changed: 840 additions & 7 deletions

File tree

‎docs/sary_prompts.md‎

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
# CA Stochastic Bifurcation Diagram:
2+
Mutations (evolutions) parameter OFF
3+
Control parameter: prey death rate
4+
5+
Possible statistical observables:
6+
- Fraction of prey cells at equilibrium
7+
- Measure of entropy of the generated pattern.
8+
- Prey population count
9+
- Predator population count
10+
11+
Run simulation:
12+
- Let the system run until a steady state is observed
13+
- For each death rate value, let the CA run for a specified number of iterations after warmp up, show distribution (scatters) for each sim run at a given prey death rate, and the average line
14+
15+
16+
# Phase 1: finding the critical point
17+
- Create bifurcation diagram of mean population count, varying prey death rate
18+
- Look for critical transition
19+
- Create log-log plot of cluster size distribution, varying prey death rate
20+
- Look for power-law
21+
22+
# Experiment Phase: CA Stochastic Bifurcation Diagram:
23+
24+
1) Write a Config Object specific to that experiment
25+
2) Make sure the experiment running on the cluster is running 15 reps of each runs at all sweeped values.
26+
3) Make sure the outputs of the experiment are a 1D and 2D array (explained below)
27+
28+
# Bifurcation Diagram Prompts:
29+
1) Help me write a function for creating a stochastic bifurcation diagram, of the population count at equilibrium, varying the prey death rate (as the control variable).
30+
2) At each sweeped value of the prey death control variable, we should be measuring the population count at equilibrium for at least 15 simulation runs.
31+
3) Which means that the two inputs for my function should be a 1D Array for the sweep parameter, and a 2D array for the experiment results at each sweep for the rows, and the results for each iteration for the columns.
32+
4) When running my function, using the argparse module, my command-line arguments specifies which analysis to do, in this case the analysis is the bifurcation diagram.
33+
34+
35+
# Output:
36+
def load_bifurcation_results(results_dir: Path) -> Tuple[np.ndarray, np.ndarray]:
37+
"""
38+
Load bifurcation analysis results.
39+
40+
Returns
41+
-------
42+
sweep_params : np.ndarray
43+
1D array of control parameter values (prey death rates).
44+
results : np.ndarray
45+
2D array of shape (n_sweep, n_replicates) with population counts
46+
at equilibrium.
47+
"""
48+
npz_file = results_dir / "bifurcation_results.npz"
49+
json_file = results_dir / "bifurcation_results.json"
50+
51+
if npz_file.exists():
52+
logging.info(f"Loading bifurcation results from {npz_file}")
53+
data = np.load(npz_file)
54+
return data['sweep_params'], data['results']
55+
elif json_file.exists():
56+
logging.info(f"Loading bifurcation results from {json_file}")
57+
with open(json_file, 'r') as f:
58+
data = json.load(f)
59+
return np.array(data['sweep_params']), np.array(data['results'])
60+
else:
61+
raise FileNotFoundError(f"Bifurcation results not found in {results_dir}")
62+
63+
64+
def plot_bifurcation_diagram(sweep_params: np.ndarray, results: np.ndarray,
65+
output_dir: Path, dpi: int = 150,
66+
control_label: str = "Prey Death Rate",
67+
population_label: str = "Population at Equilibrium"):
68+
"""
69+
Generate a stochastic bifurcation diagram.
70+
71+
Shows the distribution of equilibrium population counts as a function of
72+
a control parameter (e.g., prey death rate), with scatter points for each
73+
replicate run overlaid on summary statistics.
74+
75+
Parameters
76+
----------
77+
sweep_params : np.ndarray
78+
1D array of control parameter values (e.g., prey death rates).
79+
Shape: (n_sweep,)
80+
results : np.ndarray
81+
2D array of population counts at equilibrium.
82+
Shape: (n_sweep, n_replicates) where rows correspond to sweep_params
83+
and columns are replicate simulation runs.
84+
output_dir : Path
85+
Directory to save the output figure.
86+
dpi : int
87+
Output resolution (default: 150).
88+
control_label : str
89+
Label for x-axis (control parameter).
90+
population_label : str
91+
Label for y-axis (population count).
92+
"""
93+
n_sweep, n_replicates = results.shape
94+
95+
fig, ax = plt.subplots(figsize=(12, 7))
96+
97+
# Scatter all individual replicates with transparency
98+
for i, param in enumerate(sweep_params):
99+
ax.scatter(
100+
np.full(n_replicates, param),
101+
results[i, :],
102+
alpha=0.3, s=15, c='steelblue', edgecolors='none'
103+
)
104+
105+
# Compute summary statistics
106+
means = np.mean(results, axis=1)
107+
medians = np.median(results, axis=1)
108+
q25 = np.percentile(results, 25, axis=1)
109+
q75 = np.percentile(results, 75, axis=1)
110+
111+
# Plot median line and IQR envelope
112+
ax.fill_between(sweep_params, q25, q75, alpha=0.25, color='coral',
113+
label='IQR (25th-75th percentile)')
114+
ax.plot(sweep_params, medians, 'o-', color='darkred', linewidth=2,
115+
markersize=5, label='Median')
116+
ax.plot(sweep_params, means, 's--', color='black', linewidth=1.5,
117+
markersize=4, alpha=0.7, label='Mean')
118+
119+
ax.set_xlabel(control_label)
120+
ax.set_ylabel(population_label)
121+
ax.set_title(f"Stochastic Bifurcation Diagram\n({n_replicates} replicates per parameter value)")
122+
ax.legend(loc='best')
123+
ax.grid(True, alpha=0.3)
124+
125+
# Add rug plot at bottom showing parameter sampling density
126+
ax.plot(sweep_params, np.zeros_like(sweep_params), '|', color='gray',
127+
markersize=10, alpha=0.5)
128+
129+
plt.tight_layout()
130+
output_file = output_dir / "bifurcation_diagram.png"
131+
plt.savefig(output_file, dpi=dpi)
132+
plt.close()
133+
logging.info(f"Saved {output_file}")
134+
135+
return output_file

‎scripts/analysis.py‎

Lines changed: 122 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -6,12 +6,13 @@
66
Designed to run locally (not on HPC) for fast iteration.
77
88
Usage:
9-
python plot_pp_results.py results/ # All plots
10-
python plot_pp_results.py results/ --phase-only # Just phase diagrams
11-
python plot_pp_results.py results/ --hydra-only # Just Hydra analysis
12-
python plot_pp_results.py results/ --pcf-only # Just PCF analysis
13-
python plot_pp_results.py results/ --fss-only # Just FSS plots
14-
python plot_pp_results.py results/ --dpi 300 # High-res for publication
9+
python plot_pp_results.py results/ # All plots
10+
python plot_pp_results.py results/ --phase-only # Just phase diagrams
11+
python plot_pp_results.py results/ --hydra-only # Just Hydra analysis
12+
python plot_pp_results.py results/ --pcf-only # Just PCF analysis
13+
python plot_pp_results.py results/ --fss-only # Just FSS plots
14+
python plot_pp_results.py results/ --bifurcation-only # Just bifurcation diagram
15+
python plot_pp_results.py results/ --dpi 300 # High-res for publication
1516
"""
1617

1718
import argparse
@@ -124,6 +125,34 @@ def load_sensitivity_results(results_dir: Path) -> List[Dict]:
124125
return json.load(f)
125126

126127

128+
def load_bifurcation_results(results_dir: Path) -> Tuple[np.ndarray, np.ndarray]:
129+
"""
130+
Load bifurcation analysis results.
131+
132+
Returns
133+
-------
134+
sweep_params : np.ndarray
135+
1D array of control parameter values (prey death rates).
136+
results : np.ndarray
137+
2D array of shape (n_sweep, n_replicates) with population counts
138+
at equilibrium.
139+
"""
140+
npz_file = results_dir / "bifurcation_results.npz"
141+
json_file = results_dir / "bifurcation_results.json"
142+
143+
if npz_file.exists():
144+
logging.info(f"Loading bifurcation results from {npz_file}")
145+
data = np.load(npz_file)
146+
return data['sweep_params'], data['results']
147+
elif json_file.exists():
148+
logging.info(f"Loading bifurcation results from {json_file}")
149+
with open(json_file, 'r') as f:
150+
data = json.load(f)
151+
return np.array(data['sweep_params']), np.array(data['results'])
152+
else:
153+
raise FileNotFoundError(f"Bifurcation results not found in {results_dir}")
154+
155+
127156
# =============================================================================
128157
# DATA PROCESSING
129158
# =============================================================================
@@ -546,6 +575,80 @@ def plot_fss_analysis(fss_results: List[Dict], output_dir: Path, dpi: int = 150)
546575
logging.info(f"Saved {output_file}")
547576

548577

578+
def plot_bifurcation_diagram(sweep_params: np.ndarray, results: np.ndarray,
579+
output_dir: Path, dpi: int = 150,
580+
control_label: str = "Prey Death Rate",
581+
population_label: str = "Population at Equilibrium"):
582+
"""
583+
Generate a stochastic bifurcation diagram.
584+
585+
Shows the distribution of equilibrium population counts as a function of
586+
a control parameter (e.g., prey death rate), with scatter points for each
587+
replicate run overlaid on summary statistics.
588+
589+
Parameters
590+
----------
591+
sweep_params : np.ndarray
592+
1D array of control parameter values (e.g., prey death rates).
593+
Shape: (n_sweep,)
594+
results : np.ndarray
595+
2D array of population counts at equilibrium.
596+
Shape: (n_sweep, n_replicates) where rows correspond to sweep_params
597+
and columns are replicate simulation runs.
598+
output_dir : Path
599+
Directory to save the output figure.
600+
dpi : int
601+
Output resolution (default: 150).
602+
control_label : str
603+
Label for x-axis (control parameter).
604+
population_label : str
605+
Label for y-axis (population count).
606+
"""
607+
n_sweep, n_replicates = results.shape
608+
609+
fig, ax = plt.subplots(figsize=(12, 7))
610+
611+
# Scatter all individual replicates with transparency
612+
for i, param in enumerate(sweep_params):
613+
ax.scatter(
614+
np.full(n_replicates, param),
615+
results[i, :],
616+
alpha=0.3, s=15, c='steelblue', edgecolors='none'
617+
)
618+
619+
# Compute summary statistics
620+
means = np.mean(results, axis=1)
621+
medians = np.median(results, axis=1)
622+
q25 = np.percentile(results, 25, axis=1)
623+
q75 = np.percentile(results, 75, axis=1)
624+
625+
# Plot median line and IQR envelope
626+
ax.fill_between(sweep_params, q25, q75, alpha=0.25, color='coral',
627+
label='IQR (25th-75th percentile)')
628+
ax.plot(sweep_params, medians, 'o-', color='darkred', linewidth=2,
629+
markersize=5, label='Median')
630+
ax.plot(sweep_params, means, 's--', color='black', linewidth=1.5,
631+
markersize=4, alpha=0.7, label='Mean')
632+
633+
ax.set_xlabel(control_label)
634+
ax.set_ylabel(population_label)
635+
ax.set_title(f"Stochastic Bifurcation Diagram\n({n_replicates} replicates per parameter value)")
636+
ax.legend(loc='best')
637+
ax.grid(True, alpha=0.3)
638+
639+
# Add rug plot at bottom showing parameter sampling density
640+
ax.plot(sweep_params, np.zeros_like(sweep_params), '|', color='gray',
641+
markersize=10, alpha=0.5)
642+
643+
plt.tight_layout()
644+
output_file = output_dir / "bifurcation_diagram.png"
645+
plt.savefig(output_file, dpi=dpi)
646+
plt.close()
647+
logging.info(f"Saved {output_file}")
648+
649+
return output_file
650+
651+
549652
def plot_sensitivity_analysis(sens_results: List[Dict], output_dir: Path, dpi: int = 150):
550653
"""Generate evolution sensitivity analysis plots."""
551654
# Group by evolve_sd
@@ -694,6 +797,8 @@ def main():
694797
help='Generate only FSS plots')
695798
parser.add_argument('--sensitivity-only', action='store_true',
696799
help='Generate only sensitivity analysis plots')
800+
parser.add_argument('--bifurcation-only', action='store_true',
801+
help='Generate only bifurcation diagram')
697802
parser.add_argument('--dpi', type=int, default=150,
698803
help='Output resolution (default: 150)')
699804
parser.add_argument('--output', type=Path, default=None,
@@ -720,7 +825,7 @@ def main():
720825

721826
# Determine what to plot
722827
plot_all = not any([args.phase_only, args.hydra_only, args.pcf_only,
723-
args.fss_only, args.sensitivity_only])
828+
args.fss_only, args.sensitivity_only, args.bifurcation_only])
724829

725830
# Main sweep plots
726831
if plot_all or args.phase_only or args.hydra_only or args.pcf_only:
@@ -779,6 +884,16 @@ def main():
779884
except FileNotFoundError as e:
780885
logging.warning(f"Sensitivity results not found: {e}")
781886

887+
# Bifurcation diagram
888+
if plot_all or args.bifurcation_only:
889+
try:
890+
sweep_params, bifurc_results = load_bifurcation_results(results_dir)
891+
logging.info(f"Loaded bifurcation results: {len(sweep_params)} sweep values, "
892+
f"{bifurc_results.shape[1]} replicates each")
893+
plot_bifurcation_diagram(sweep_params, bifurc_results, output_dir, args.dpi)
894+
except FileNotFoundError as e:
895+
logging.warning(f"Bifurcation results not found: {e}")
896+
782897
logging.info("Done!")
783898

784899

0 commit comments

Comments
 (0)