From 7fe81680f253f6ce951275da08221055f388daf4 Mon Sep 17 00:00:00 2001 From: Yihui Cai Date: Wed, 23 Aug 2023 14:55:39 +0800 Subject: [PATCH] Parameterized sample rate for risk index caculation; Added UT --- simglucose/analysis/report.py | 22 ++++++++----- simglucose/simulation/user_interface.py | 5 +-- tests/test_report.py | 42 +++++++++++++++++++++++++ 3 files changed, 59 insertions(+), 10 deletions(-) create mode 100644 tests/test_report.py diff --git a/simglucose/analysis/report.py b/simglucose/analysis/report.py index f0bec0e8..34f808ff 100644 --- a/simglucose/analysis/report.py +++ b/simglucose/analysis/report.py @@ -8,6 +8,8 @@ # from pandas.plotting import lag_plot import logging +from simglucose.sensor.cgm import CGMSensor + logger = logging.getLogger(__name__) @@ -92,18 +94,18 @@ 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 + 20, :] for i in range(0, len(df_BG), 20)] +def risk_index_trace(df_BG, sample_rate, visualize=False): + chunk_BG = [df_BG.iloc[i:i + sample_rate, :] for i in range(0, len(df_BG), sample_rate)] - if len(chunk_BG[-1]) != 20: # Remove the last chunk which is not full + if len(chunk_BG[-1]) != sample_rate: # Remove the last chunk which is not full chunk_BG.pop() fBG = [ 1.509 * (np.log(BG[BG > 0]) ** 1.084 - 5.381) for BG in chunk_BG ] - rl = [(10 * (fbg * (fbg < 0))**2).mean() for fbg in fBG] - rh = [(10 * (fbg * (fbg > 0))**2).mean() for fbg in fBG] + rl = [(10 * (fbg * (fbg < 0)) ** 2).mean() for fbg in fBG] + rh = [(10 * (fbg * (fbg > 0)) ** 2).mean() for fbg in fBG] LBGI = pd.concat(rl, axis=1).transpose() HBGI = pd.concat(rh, axis=1).transpose() @@ -238,7 +240,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']) @@ -247,12 +249,16 @@ 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 + sample_rate = 20 # Use 20 as default sample rate in one hour + if cgm_sensor is not None: + sample_rate = int(60 / cgm_sensor.sample_time) + 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) + ri_per_hour, ri_mean, fig_ri, ax5 = risk_index_trace(BG, sample_rate, 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] diff --git a/simglucose/simulation/user_interface.py b/simglucose/simulation/user_interface.py index 6826406f..a6f460fc 100644 --- a/simglucose/simulation/user_interface.py +++ b/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 diff --git a/tests/test_report.py b/tests/test_report.py new file mode 100644 index 00000000..64792516 --- /dev/null +++ b/tests/test_report.py @@ -0,0 +1,42 @@ +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 + sample_rate = int(60 / sample_time) + ri_per_hour, ri_mean, fig, axes = risk_index_trace(BG, sample_rate) + + 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(LBGI.iloc[-1].test, 0.8429957158900777) + self.assertEqual(LBGI.iloc[0].test, 0.0) + + self.assertEqual(HBGI.size, 48) + self.assertEqual(HBGI.iloc[-1].test, 0.0) + self.assertEqual(HBGI.iloc[0].test, 2.755277346918188) + + self.assertEqual(RI.size, 48) + self.assertEqual(RI.iloc[-1].test, 0.8429957158900777) + self.assertEqual(RI.iloc[0].test, 2.755277346918188) + + +if __name__ == '__main__': + unittest.main()