update model output return

Author: KOVALW

experiments/phase1/reports/calibration.qmd

@@ -52,12 +52,12 @@ print(
 )
 ```
 
-Histograms of posterior samples compared to the probability density of each prior. To generate the histogram, we sample $n$ particles from the posterior population, where $n$ is the effective sample size of the population, in order to reflect the weight distribution of the posterior. For each plot, the prior is plotted as a blue line and the histogram of the posterior samples is plotted in orange with a kernel density estimator overlay.
+Histograms of posterior samples compared to the probability density of each prior. For each plot, the prior is plotted as a blue line and the histogram of the posterior samples is plotted in orange with a kernel density estimator overlay. To generate the histogram, we sample $n$ particles from the posterior population, where $n$ is the effective sample size of the population. Re-sampling particles in this manner reflects the posterior weight distribution, the probability mass function over accepted particles determined by the ABC-SMC algorithm perturbation kernel and prior distributions.
 
 ```{python}
 #| echo: false
 #| eval: false
-# Preferred API ----------
+# Preferred API for development in calibrationtools----------
 for param in results.fitted_params:
     # get values and use weights to alter histogram instead of sampling with ESS
     values = results.posterior_particles.get_values_of(param)
@@ -120,7 +120,7 @@ for param in results.fitted_params:
 ```
 ```{python}
 #| echo: false
-## Obtaining the importation time series and death incidence data frmaes for a random smaple form the posterior
+## Obtaining the importation time series and death incidence data frames for a random smaple form the posterior
 # Re-generating a random sample of parameter sets from posterior
 particle_count = int(min(100, results.ess))
 particles = results.sample_posterior_particles(n=particle_count)
@@ -129,18 +129,22 @@ default_params_file = wd / 'input' / 'default_params.json'
 with open(default_params_file, "rb") as fp:
     default_params = json.load(fp)
 
-default_params['epimodel.GlobalParams']['max_time'] = 200
 ixa_overrides = {
-    "synth_population_file": "/mnt/S_CFA_Predict/team-CMEI/synthetic_populations/cbsa_all_work_school_household_2020-04-24/cbsa_all_work_school_household/IN/Bloomington IN.csv"
+    "synth_population_file": "/mnt/S_CFA_Predict/team-CMEI/synthetic_populations/cbsa_all_work_school_household_2020-04-24/cbsa_all_work_school_household/IN/Bloomington IN.csv",
+    "imported_cases_timeseries": {
+        "filename": "./experiments/phase1/projection/imported_cases_timeseries.csv"
+    },
+    "max_time": 200
 }
 default_params = apply_dict_overrides(default_params, {'epimodel.GlobalParams': ixa_overrides})
 
 mrp_defaults = {
     'ixa_inputs': default_params,
     "config_inputs": {
         "exe_file": "./target/release/ixa-epi-covid",
-        "output_dir": "./experiments/phase1/calibration/output",
+        "output_dir": "./experiments/phase1/projection/output",
         "force_overwrite": True,
+        "outputs_to_read": ['prevalence_report', 'imported_cases_timeseries']
     },
     "importation_inputs": {
         "state": "Indiana",
@@ -153,41 +157,16 @@ uniq_id = 0
 model = CovidModel()
 importation_curves = []
 prevalence_data = []
+os.makedirs(mrp_defaults["config_inputs"]["output_dir"], exist_ok=True)
+
+reader = ParticleReader(results.priors.params, mrp_defaults)
+for p in particles:
+    model_inputs = reader.read_particle(p)
+    outputs = model.simulate(model_inputs)
 
-with tempfile.TemporaryDirectory() as tmpdir:
-    reader = ParticleReader(results.priors.params, mrp_defaults)
-    for p in particles:
-        model_inputs = reader.read_particle(p)
-
-        model_inputs["config_inputs"]["output_dir"] = str(
-            Path(tmpdir, f"{uniq_id}")
-        )
-        os.makedirs(model_inputs["config_inputs"]["output_dir"], exist_ok=True)
-        importation_path = Path(
-            tmpdir, f"{uniq_id}", "importation_timeseries.csv"
-        )
-
-        model_inputs["ixa_inputs"]["epimodel.GlobalParams"][
-            "imported_cases_timeseries"
-        ]["filename"] = str(importation_path)
-        model.simulate(model_inputs)
-        prevalence_data.append(
-            pl.read_csv(
-                Path(
-                    tmpdir,
-                    f"{uniq_id}",
-                    model_inputs["ixa_inputs"]["epimodel.GlobalParams"][
-                        "prevalence_report"
-                    ]["filename"],
-                )
-            ).with_columns(pl.lit(uniq_id).alias("id"))
-        )
-        importation_curves.append(
-            pl.read_csv(importation_path).with_columns(
-                pl.lit(uniq_id).alias("id")
-            )
-        )
-        uniq_id += 1
+    importation_curves.append(outputs["imported_cases_timeseries"].with_columns(pl.lit(uniq_id).alias("id")))
+    prevalence_data.append(outputs["prevalence_report"].with_columns(pl.lit(uniq_id).alias("id")))
+    uniq_id += 1
 
 importations = pl.concat(importation_curves)
 all_prevalence_data = pl.concat(prevalence_data)

scripts/phase_1_calibration.py

@@ -26,26 +26,32 @@
 )
 ixa_overrides = {
     "synth_population_file": "/mnt/S_CFA_Predict/team-CMEI/synthetic_populations/cbsa_all_work_school_household_2020-04-24/cbsa_all_work_school_household/IN/Bloomington IN.csv",
-    "first_death_terminates_run": True
+    "first_death_terminates_run": True,
 }
 force_overwrite = False
+outputs_to_read = ["incidence_report"]
 
 # State importation model declaration parameters
 state = "Indiana"
 year = 2020
+symptomatic_reporting_prob_default = 0.5
 
 # Calibration inputs
 priors_file = Path("experiments", "phase1", "input", "priors.json")
-tolerance_values = [2.0, 0.1]  # , 2.0, 0.01]
+tolerance_values = [2.0, 0.1]
 generation_particle_count = 500
 target_data = 75
 
 
 # Output processing function for calibration
-def outputs_to_distance(model_output: pl.DataFrame, target_data: int):
-    first_death_observed = model_output.filter(
-        (pl.col("event") == "Dead") & (pl.col("count") > 0)
-    ).filter(pl.col("t_upper") == pl.min("t_upper"))
+def outputs_to_distance(
+    model_output: dict[str, pl.DataFrame], target_data: int
+):
+    first_death_observed = (
+        model_output["incidence_report"]
+        .filter((pl.col("event") == "Dead") & (pl.col("count") > 0))
+        .filter(pl.col("t_upper") == pl.min("t_upper"))
+    )
     if first_death_observed.height > 0:
         return abs(target_data - first_death_observed.item(0, "t_upper"))
     else:
@@ -74,11 +80,12 @@ def outputs_to_distance(model_output: pl.DataFrame, target_data: int):
         "exe_file": str(exe_file),
         "output_dir": str(output_dir),
         "force_overwrite": force_overwrite,
+        "outputs_to_read": outputs_to_read,
     },
     "importation_inputs": {
-        "state": "Indiana",
-        "year": 2020,
-        "symptomatic_reporting_prob": 0.5,
+        "state": state,
+        "year": year,
+        "symptomatic_reporting_prob": symptomatic_reporting_prob_default,
     },
 }
 

src/ixa_epi_covid/covid_model.py

@@ -85,7 +85,7 @@ def simulate(model_inputs: dict[str, Any]) -> pl.DataFrame:
             str(input_file_path),
             "--output",
             config_inputs["output_dir"],
-            "-f",
+            "--force-overwrite",
             "--no-stats",
         ]
 
@@ -99,11 +99,21 @@ def simulate(model_inputs: dict[str, Any]) -> pl.DataFrame:
             raise e
 
         # Read the model incidence report from the specified location and return as a DataFrame
-        incidence_report_filename = ixa_inputs["epimodel.GlobalParams"][
-            "incidence_report"
-        ]["filename"]
-        incidence_report_path = Path(
-            config_inputs["output_dir"], incidence_report_filename
-        )
-
-        return pl.read_csv(incidence_report_path)
+        outputs = {}
+        for output in config_inputs["outputs_to_read"]:
+            fp = ixa_inputs["epimodel.GlobalParams"][output]["filename"]
+            if Path(config_inputs["output_dir"], fp).exists():
+                outputs.update(
+                    {
+                        output: pl.read_csv(
+                            Path(config_inputs["output_dir"], fp)
+                        )
+                    }
+                )
+            elif Path(fp).exists():
+                outputs.update({output: pl.read_csv(Path(fp))})
+            else:
+                raise FileNotFoundError(
+                    f"Expected output file {fp} not found. Looked in {config_inputs['output_dir']}"
+                )
+        return outputs

src/model.rs

@@ -1,8 +1,8 @@
 use ixa::{ExecutionPhase, prelude::*};
 
 use crate::{
-    infection_importation, infection_propagation_loop, population_loader, reports, settings,
-    symptom_status_manager, abort_run
+    abort_run, infection_importation, infection_propagation_loop, population_loader, reports,
+    settings, symptom_status_manager,
 };
 
 pub fn initialize_model(context: &mut Context, seed: u64, max_time: f64) -> Result<(), IxaError> {