hospital usage

src/covid19sim/log/track.py

@@ -286,6 +286,7 @@ def __init__(self, env, city, conf, logfile):
 
         # testing & hospitalization
         self.hospitalization_per_day = [0]
+        self.hospital_usage_per_day = [0]
         self.critical_per_day = [0]
         self.deaths_per_day = [0]
         self.test_results_per_day = defaultdict(lambda :{'positive':0, 'negative':0})
@@ -654,6 +655,7 @@ def increment_day(self):
         self.r_per_day.append(sum(h.is_removed for h in self.city.humans))
         self.ei_per_day.append(self.e_per_day[-1] + self.i_per_day[-1])
 
+        num_humans_in_hospital = 0
         for human in self.city.humans:
             if human.is_susceptible:
                 state = 'S'
@@ -665,7 +667,8 @@ def increment_day(self):
                 state = 'R'
             else:
                 raise ValueError(f"{human} is not in any of SEIR states")
-
+            if human.mobility_planner.hospitalization_timestamp is not None:
+                num_humans_in_hospital += 1
             self.humans_quarantined_state[human.name].append(human.intervened_behavior.is_under_quarantine)
             self.humans_state[human.name].append(state)
             self.humans_rec_level[human.name].append(human.rec_level)
@@ -674,6 +677,7 @@ def increment_day(self):
         # test_per_day
         self.tested_per_day.append(0)
         self.hospitalization_per_day.append(0)
+        self.hospital_usage_per_day.append(num_humans_in_hospital)
         self.critical_per_day.append(0)
         self.deaths_per_day.append(0)
 
@@ -1834,7 +1838,7 @@ def write_for_training(self, humans, outfile, conf):
         """ Writes some data out for the ML predictor """
         data = dict()
         data['hospitalization_per_day'] = self.hospitalization_per_day # how many new people get admitted to the hospital
-        import pdb; pdb.set_trace()
+        data['hospital_usage_per_day'] = self.hospital_usage_per_day # how many people are in the hospital for covid
         # TODO: write total number of _currently hospitalized people_, not how many people _go into the hospital_.
 
         # parse test results

src/covid19sim/plotting/compare_real_and_sim.py

@@ -1,4 +1,4 @@
-import os
+import os, sys
 import pickle
 import numpy as np
 import pandas as pd
@@ -10,6 +10,8 @@
 quebec_population = 8485000
 csv_path = "path/to/csv"
 sims_dir_path = "path/to/simulations/"
+if len(sys.argv) > 1:
+    sims_dir_path = sys.argv[1]
 
 # Load data
 qc_data = pd.read_csv(csv_path)
@@ -62,28 +64,28 @@ def parse_tracker(sim_tracker_data):
     sim_prior_data = pickle.load(open(sim_priors_path, "rb"))
     # Parse data
     sim_dates, sim_deaths, sim_tests, sim_cases = parse_tracker(sim_tracker_data)
-    sim_hospitalizations = [float(x)*100/sim_tracker_data['n_humans'] for x in sim_prior_data['hospitalization_per_day']]
+    sim_hospitalizations = [float(x)*100/sim_tracker_data['n_humans'] for x in sim_prior_data['hospital_usage_per_day']]
     # change key above in sim_prior_data['hospital_usage_per_day']
 
     all_sim_cases.append(sim_cases)
     all_sim_hospitalizations.append(sim_hospitalizations)
     all_sim_deaths.append(sim_deaths)
 
-avg_sim_cases = [sum(elem)/len(elem) for elem in zip(*all_sim_cases)]
-avg_sim_hospitalizations = [sum(elem)/len(elem) for elem in zip(*all_sim_hospitalizations)]
-avg_sim_deaths = [sum(elem)/len(elem) for elem in zip(*all_sim_deaths)]
+avg_sim_cases = np.array([sum(elem)/len(elem) for elem in zip(*all_sim_cases)])
+avg_sim_hospitalizations = np.array([sum(elem)/len(elem) for elem in zip(*all_sim_hospitalizations)])
+avg_sim_deaths = np.array([sum(elem)/len(elem) for elem in zip(*all_sim_deaths)])
 
 # Plot Quebec Data
 last_index = 63 # index of last day of Quebec data, use 63 for 30 days
 fig, ax = plt.subplots(figsize=(7.5, 7.5))
 real_dates = qc_data.loc[34:last_index, 'dates'].to_numpy()
 real_cases = [100 * float(x if str(x) != "nan" else 0) / quebec_population for x in qc_data.loc[34:last_index, 'change_cases']]
 # change key to 'total_hospitalizations' if needed
-real_hospitalizations = [100 * float(x if str(x) != "nan" else 0) / quebec_population for x in qc_data.loc[34:last_index, 'change_hospitalizations']]
+real_hospitalizations = [100 * float(x if str(x) != "nan" else 0) / quebec_population for x in qc_data.loc[34:last_index, 'total_hospitalizations']]
 real_deaths = [100 * float(x if str(x) != "nan" else 0) / quebec_population for x in qc_data.loc[34:last_index,'change_fatalities']]
 
 ax.plot(real_dates, real_cases, label="Diagnoses per day")
-ax.plot(real_dates, real_hospitalizations, label="Hospital admissions")
+ax.plot(real_dates, real_hospitalizations, label="Hospital utilization")
 ax.plot(real_dates, real_deaths, label="Mortalities per day")
 
 ax.legend()
@@ -94,20 +96,18 @@ def parse_tracker(sim_tracker_data):
 plt.title("Real Quebec COVID Statistics")
 plt.savefig("qc_stats.png")
 
-# Plot average across simulations data
+# Plot average across simulations
 fig, ax = plt.subplots(figsize=(7.5, 7.5))
 # ax.set_ylim([0.,0.03])
 avg_sim_cases = np.array(sim_cases[1:])*.1 # adjust for unknown cases
 
-'''
-ax.errorbar(np.array(sim_dates), avg_sim_cases.mean(axis=0), yerr=avg_sim_cases.std(axis=0), label="Diagnoses per day")
+ax.errorbar(sim_dates, avg_sim_cases, yerr=avg_sim_cases.std(axis=0), label="Diagnoses per day")
 # change caption below if using 'total_hospitalizations'
-ax.errorbar(sim_dates, avg_sim_hospitalizations.mean(axis=0)[1:], yerr=avg_sim_hospitalizations.std(axis=0)[1:], label="Hospital admissions")
-ax.errorbar(sim_dates, avg_sim_deaths.mean(axis=0), yerr=avg_sim_deaths.std(axis=0), label="Mortalities per day")
-'''
+ax.errorbar(sim_dates, avg_sim_hospitalizations[1:], yerr=avg_sim_hospitalizations.std(axis=0), label="Hospital admissions")
+ax.errorbar(sim_dates, avg_sim_deaths, yerr=avg_sim_deaths.std(axis=0), label="Mortalities per day")
 
 ax.plot(sim_dates, avg_sim_cases, label="Diagnoses per day")
-ax.plot(sim_dates, avg_sim_hospitalizations[1:], label="Hospital admissions")
+ax.plot(sim_dates, avg_sim_hospitalizations[1:], label="Hospital utilization")
 ax.plot(sim_dates, avg_sim_deaths, label="Mortalities per day")
 
 ax.legend()