Merge pull request #62 from yihuicai/dev-alan-fix-report-issue

Fix Risk Index calculation and add data sanitization

simglucose/analysis/report.py

@@ -92,17 +92,22 @@ def percent_stats(BG, ax=None):
     return p_stats, fig, ax
 
 
-def risk_index_trace(df_BG, visualize=False):
-    chunk_BG = [df_BG.iloc[i:i + 60, :] for i in range(0, len(df_BG), 60)]
+def risk_index_trace(df_BG, sample_time=3, window_length=60, visualize=False):
+    step_size = int(window_length / sample_time)  # window size set to 1 hour for calculating Risk Index
+    chunk_BG = [df_BG.iloc[i:i + step_size, :] for i in range(0, len(df_BG), step_size)]
+
+    if len(chunk_BG[-1]) != step_size:  # Remove the last chunk which is not full
+        chunk_BG.pop()
 
     fBG = [
-        np.mean(1.509 * (np.log(BG[BG > 0])**1.084 - 5.381)) for BG in chunk_BG
+        1.509 * (np.log(BG[BG > 0]) ** 1.084 - 5.381) for BG in chunk_BG
     ]
 
-    fBG_df = pd.concat(fBG, axis=1).transpose()
+    rl = [(10 * (fbg * (fbg < 0)) ** 2).mean() for fbg in fBG]
+    rh = [(10 * (fbg * (fbg > 0)) ** 2).mean() for fbg in fBG]
 
-    LBGI = 10 * (fBG_df * (fBG_df < 0))**2
-    HBGI = 10 * (fBG_df * (fBG_df > 0))**2
+    LBGI = pd.concat(rl, axis=1).transpose()
+    HBGI = pd.concat(rh, axis=1).transpose()
     RI = LBGI + HBGI
 
     ri_per_hour = pd.concat(
@@ -234,7 +239,7 @@ def CVGA(BG_list, label=None):
             edgecolors='k',
             zorder=4,
             label='%s (A: %d%%, B: %d%%, C: %d%%, D: %d%%, E: %d%%)' %
-            (l, 100 * A, 100 * B, 100 * C, 100 * D, 100 * E))
+                  (l, 100 * A, 100 * B, 100 * C, 100 * D, 100 * E))
         zone_stats.append((A, B, C, D, E))
 
     zone_stats = pd.DataFrame(zone_stats, columns=['A', 'B', 'C', 'D', 'E'])
@@ -243,12 +248,15 @@ def CVGA(BG_list, label=None):
     return zone_stats, fig, ax
 
 
-def report(df, save_path=None):
+def report(df, cgm_sensor=None, save_path=None):
     BG = df.unstack(level=0).BG
 
     fig_ensemble, ax1, ax2, ax3 = ensemblePlot(df)
     pstats, fig_percent, ax4 = percent_stats(BG)
-    ri_per_hour, ri_mean, fig_ri, ax5 = risk_index_trace(BG, visualize=False)
+    if cgm_sensor is not None:
+        ri_per_hour, ri_mean, fig_ri, ax5 = risk_index_trace(BG, sample_time=cgm_sensor.sample_time, visualize=False)
+    else:
+        ri_per_hour, ri_mean, fig_ri, ax5 = risk_index_trace(BG, visualize=False)
     zone_stats, fig_cvga, ax6 = CVGA(BG, label='')
     axes = [ax1, ax2, ax3, ax4, ax5, ax6]
     figs = [fig_ensemble, fig_percent, fig_ri, fig_cvga]

simglucose/simulation/user_interface.py

@@ -361,9 +361,10 @@ def simulate(sim_time=None,
     if controller is None:
         controller = pick_controller()
 
+    cgm_sensor = CGMSensor.withName(cgm_name, seed=cgm_seed)
+
     def local_build_env(pname):
         patient = T1DPatient.withName(pname)
-        cgm_sensor = CGMSensor.withName(cgm_name, seed=cgm_seed)
         insulin_pump = InsulinPump.withName(insulin_pump_name)
         scen = copy.deepcopy(scenario)
         env = T1DSimEnv(patient, cgm_sensor, insulin_pump, scen)
@@ -380,7 +381,7 @@ def local_build_env(pname):
     results = batch_sim(sim_instances, parallel=parallel)
 
     df = pd.concat(results, keys=[s.env.patient.name for s in sim_instances])
-    results, ri_per_hour, zone_stats, figs, axes = report(df, save_path)
+    results, ri_per_hour, zone_stats, figs, axes = report(df, cgm_sensor, save_path)
 
     return results
 

tests/test_report.py

@@ -0,0 +1,41 @@
+from simglucose.analysis.report import risk_index_trace
+from simglucose.sensor.cgm import CGMSensor
+from simglucose.simulation.rendering import Viewer
+from datetime import datetime
+import pandas as pd
+import unittest
+import logging
+import os
+
+logger = logging.getLogger(__name__)
+TESTDATA_FILENAME = os.path.join(os.path.dirname(__file__), 'sim_results.csv')
+
+
+class TestReport(unittest.TestCase):
+    def setUp(self):
+        self.df = pd.concat([pd.read_csv(TESTDATA_FILENAME, index_col=0)],  keys=['test'])
+
+    def test_risk_index_trace(self):
+        BG = self.df.unstack(level=0).BG
+        sample_time = CGMSensor.withName("Dexcom").sample_time
+        ri_per_hour, ri_mean, fig, axes = risk_index_trace(BG, sample_time=sample_time)
+
+        LBGI = ri_per_hour.transpose().LBGI
+        HBGI = ri_per_hour.transpose().HBGI
+        RI = ri_per_hour.transpose()["Risk Index"]
+
+        self.assertEqual(LBGI.size, 48)
+        self.assertEqual(round(LBGI.iloc[-1].test, 3), 0.843)
+        self.assertEqual(round(LBGI.iloc[0].test, 3), 0.0)
+
+        self.assertEqual(HBGI.size, 48)
+        self.assertEqual(round(HBGI.iloc[-1].test,3), 0.0)
+        self.assertEqual(round(HBGI.iloc[0].test,3), 2.755)
+
+        self.assertEqual(RI.size, 48)
+        self.assertEqual(round(RI.iloc[-1].test,3), 0.843)
+        self.assertEqual(round(RI.iloc[0].test,3), 2.755)
+
+
+if __name__ == '__main__':
+    unittest.main()