diff --git a/.gitignore b/.gitignore index 14e1f5b..3c358cc 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,6 @@ out/sim-*/plots out/sim-*/dataset-*.hdf5 lib venv -*~ \ No newline at end of file +*~ +out/*/.DS_Store +.DS_Store \ No newline at end of file diff --git a/config/simulation-charmeleon-alt.json b/config/simulation-charmeleon-alt.json index 40ba792..a422ddd 100644 --- a/config/simulation-charmeleon-alt.json +++ b/config/simulation-charmeleon-alt.json @@ -4,12 +4,12 @@ "seed": 42, "remaster_xml": "src/remaster-template.xml", "num_simulations": 2000, - "num_workers": 16, + "num_workers": 5, "simulation_hyperparameters": { "duration_range": { "dist": "uniform_int", - "lower_bound": 30, - "upper_bound": 50 + "lower_bound": 10, + "upper_bound": 30 }, "num_changes": { "dist": "uniform_int", diff --git a/config/simulation-wartortle.json b/config/simulation-wartortle.json new file mode 100644 index 0000000..ffeba01 --- /dev/null +++ b/config/simulation-wartortle.json @@ -0,0 +1,37 @@ +{ + "simulation_name": "sim-wartortle", + "output_hdf5": "dataset-wartortle.hdf5", + "seed": 42, + "remaster_xml": "src/remaster-template-contemporaneous-lenient.xml", + "num_simulations": 2000, + "num_workers": 5, + "simulation_hyperparameters": { + "duration_range": { + "dist": "uniform_int", + "lower_bound": 10, + "upper_bound": 30 + }, + "num_changes": { + "dist": "uniform_int", + "lower_bound": 1, + "upper_bound": 2 + }, + "r0": { + "dist": "lognormal", + "LN_mean": 1.17, + "LN_sigma": 0.38 + }, + "net_removal_rate": { + "dist": "lognormal", + "LN_mean": -1.81, + "LN_sigma": 0.2 + }, + "sampling_prop": { + "dist": "beta", + "alpha": 2, + "beta": 10 + }, + "num_temp_measurements": 10, + "contemporaneous_sample": true + } +} diff --git a/main.py b/main.py index 2281341..b32d9c9 100644 --- a/main.py +++ b/main.py @@ -643,6 +643,11 @@ def create_database(pickle_files): "present_cumulative", data=sim["simulation_results"]["present_cumulative"], ) + + #put rho in database if doing contemporaneous sampling + if "rho" in sim["parameters"].keys(): + params_grp.create_dataset("rho", data=sim["parameters"]["rho"]["values"]) + db_conn.attrs["num_simulations"] = num_sims db_conn.attrs["creation_date"] = datetime.datetime.now().isoformat() db_conn.close() diff --git a/src/get_output_distributions.py b/src/get_output_distributions.py new file mode 100644 index 0000000..85446d9 --- /dev/null +++ b/src/get_output_distributions.py @@ -0,0 +1,125 @@ +import h5py +import numpy as np + +#look at an hdf5 file, pull out present prevalence and cumulative incidence + + + +DB_FNAME = "out/sim-wartortle/dataset-wartortle.hdf5" +#DB_FNAME = "out/sim-charmeleon-alt/dataset-charmeleon-alt.hdf5" + +print(f"\nReading from {DB_FNAME}\n") + + +#read in +db_conn = h5py.File(DB_FNAME, "r") + +epidemic_durations = np.array([db_conn[f"{key}/output/parameters/epidemic_duration"][()] for key in db_conn.keys()]) + +r0_values = [db_conn[f"{key}/output/parameters/r0/values"][()] for key in db_conn.keys()] +r0_ct = [db_conn[f"{key}/output/parameters/r0/change_times"][()] for key in db_conn.keys()] + +net_removal_rates = [db_conn[f"{key}/output/parameters/net_removal_rate/values"][()] for key in db_conn.keys()] +net_removal_ct = [db_conn[f"{key}/output/parameters/net_removal_rate/change_times"][()] for key in db_conn.keys()] + +sampling_prop_values = [db_conn[f"{key}/output/parameters/sampling_prop/values"][()] for key in db_conn.keys()] +sampling_prop_ct = [db_conn[f"{key}/output/parameters/sampling_prop/change_times"][()] for key in db_conn.keys()] + +if "rho" in db_conn[f"record_000001/output/parameters/"].keys(): + rho_values = np.array([db_conn[f"{key}/output/parameters/rho/values"][()] for key in db_conn.keys()]) + +num_param_changes = np.array([len(r0_ct[ix]) for ix in range(len(r0_ct))]) +all_change_times = np.concatenate(r0_ct) + +present_prevalence = np.array([db_conn[f"{key}/output/present_prevalence"][()] for key in db_conn.keys()]) +present_cumulative = np.array([db_conn[f"{key}/output/present_cumulative"][()] for key in db_conn.keys()]) + + + + + +#process parameter samples +all_r0_samples = np.concatenate(r0_values) +r0_weights = np.concatenate([ + np.diff(np.concatenate((np.array([0.0]), change_times, np.array([epidemic_durations[ix]])))) for (ix, change_times) in enumerate(r0_ct) +]) + +all_net_removal_samples = np.concatenate(net_removal_rates) +net_removal_weights = np.concatenate([ + np.diff(np.concatenate((np.array([0.0]), change_times, np.array([epidemic_durations[ix]])))) for (ix, change_times) in enumerate(net_removal_ct) +]) + +all_sampling_prop_samples = np.concatenate(sampling_prop_values) +sampling_prop_weights = np.concatenate([ + np.diff(np.concatenate((np.array([0.0]), change_times, np.array([epidemic_durations[ix]])))) for (ix, change_times) in enumerate(sampling_prop_ct) +]) + + + + + + + + +# MODEL PARAMS +print("\nModel Parameter Distributions:") + +#r0 +r0_median = np.median(all_r0_samples) +r0_interval = np.percentile(all_r0_samples, [2.5, 97.5], weights=r0_weights, method='inverted_cdf') +print(f"R0: Median={r0_median:.4f}, 95% interval=({r0_interval[0]:.4f}, {r0_interval[1]:.4f})") + +#net removal rate +net_removal_median = np.median(all_net_removal_samples) +net_removal_interval = np.percentile(all_net_removal_samples, [2.5, 97.5], weights=net_removal_weights, method='inverted_cdf') +print(f"Net Removal Rate: Median={net_removal_median:.4f}, 95% interval=({net_removal_interval[0]:.4f}, {net_removal_interval[1]:.4f})") + +#sampling proportion +sampling_prop_median = np.median(all_sampling_prop_samples) +sampling_prop_interval = np.percentile(all_sampling_prop_samples, [2.5, 97.5], weights=sampling_prop_weights, method='inverted_cdf') +print(f"Sampling Proportion: Median={sampling_prop_median:.4f}, 95% interval=({sampling_prop_interval[0]:.4f}, {sampling_prop_interval[1]:.4f})") + +#rho if present +if "rho" in db_conn[f"record_000001/output/parameters/"].keys(): + rho_median = np.median(rho_values) + rho_interval = np.percentile(rho_values, [2.5, 97.5]) + print(f"Rho (Contemporaneous Sampling Proportion): Median={rho_median:.4f}, 95% interval=({rho_interval[0]:.4f}, {rho_interval[1]:.4f})") + + + + +#HYPERPARAMS +print("\nHyperparameter Distributions:") + +#epidemic duration +epidemic_duration_median = np.median(epidemic_durations) +epidemic_duration_interval = np.percentile(epidemic_durations, [2.5, 97.5]) +print(f"Epidemic Duration: Median={epidemic_duration_median:.4f}, 95% interval=({epidemic_duration_interval[0]:.4f}, {epidemic_duration_interval[1]:.4f})") + +#number of parameter changes +num_param_changes_median = np.median(num_param_changes) +num_param_changes_interval = np.percentile(num_param_changes, [2.5, 97.5]) +print(f"Number of Parameter Changes: Median={num_param_changes_median:.4f}, 95% interval=({num_param_changes_interval[0]:.4f}, {num_param_changes_interval[1]:.4f})") + +#change times +all_change_times_median = np.median(all_change_times) +all_change_times_interval = np.percentile(all_change_times, [2.5, 97.5]) +print(f"Change Times: Median={all_change_times_median:.4f}, 95% interval=({all_change_times_interval[0]:.4f}, {all_change_times_interval[1]:.4f})") + + + + +#OUTPUT +print("\nOutput Distributions:") + +#prevalence +prevalence_median = np.median(present_prevalence) +prevalence_interval = np.percentile(present_prevalence, [2.5, 97.5]) +print(f"Present Prevalence: Median={prevalence_median:.4f}, 95% interval=({prevalence_interval[0]:.4f}, {prevalence_interval[1]:.4f})") + + +#cumulative incidence +cumulative_median = np.median(present_cumulative) +cumulative_interval = np.percentile(present_cumulative, [2.5, 97.5]) +print(f"Present Cumulative Incidence: Median={cumulative_median:.4f}, 95% interval=({cumulative_interval[0]:.4f}, {cumulative_interval[1]:.4f})") +print("\n") \ No newline at end of file diff --git a/src/remaster-template-contemporaneous-lenient.xml b/src/remaster-template-contemporaneous-lenient.xml new file mode 100644 index 0000000..adcfe97 --- /dev/null +++ b/src/remaster-template-contemporaneous-lenient.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + X -> 2X + X -> Psi + X -> Mu + + + + + + + + + + + + + diff --git a/src/remaster-template-contemporaneous.xml b/src/remaster-template-contemporaneous.xml index baf16a2..5ed353e 100644 --- a/src/remaster-template-contemporaneous.xml +++ b/src/remaster-template-contemporaneous.xml @@ -15,7 +15,7 @@ X -> 2X - X -> Psi + X -> Psi X -> Mu