Merge pull request #40 from lcossu/CGM-abstraction

Abstraction of CGM sensor error model

Author: gcappon

docs/.vuepress/config.js

@@ -1,75 +1,79 @@
-import { defineUserConfig } from 'vuepress/cli'
-import { viteBundler } from '@vuepress/bundler-vite'
-import { hopeTheme } from "vuepress-theme-hope";
+import {defineUserConfig} from 'vuepress/cli'
+import {viteBundler} from '@vuepress/bundler-vite'
+import {hopeTheme} from "vuepress-theme-hope";
 
 export default defineUserConfig({
-  lang: 'en-US',
+    lang: 'en-US',
 
   title: 'ReplayBG',
   base: '/py_replay_bg/',
   description: 'A digital twin based framework for the development and assessment of new algorithms for type 1 ' +
     'diabetes management',
 
-  theme: hopeTheme({
-    logo: 'https://i.postimg.cc/gJn8Sy0X/replay-bg-logo.png',
-    navbar: ['/', '/documentation/get_started', '/documentation/'],
-    repo: 'https://github.com/gcappon/py_replay_bg',
-    docsDir: 'docs/',
-    docsBranch: 'master',
-    markdown: {
-      highlighter: {
-        type: "shiki",
-        langs: ['python', 'json', 'md', 'bash', 'diff'],
-        themes: {
-          dark: 'catppuccin-mocha',
-          light: 'catppuccin-latte'
-        }
-      },
-      math: {
-        type: "mathjax",
-      },
-    },
-    sidebar: [
-      {
-        text: 'Get Started',
-        link: 'documentation/get_started.md'
-      },
-      {
-        text: 'Data Requirements',
-        link: 'documentation/data_requirements.md'
-      },
-      {
-        text: 'The ReplayBG Object',
-        link: 'documentation/replaybg_object.md'
-      },
-      {
-        text: 'Choosing Blueprint',
-        link: 'documentation/choosing_blueprint.md'
-      },
-      {
-        text: 'Twinning Procedure',
-        link: 'documentation/twinning_procedure.md'
-      },
-      {
-        text: 'Replaying',
-        link: 'documentation/replaying.md'
-      },
-      {
-        text: 'The _results/_ Folder',
-        link: 'documentation/results_folder.md'
-      },
-      {
-        text: 'Visualizing Replay Results',
-        link: 'documentation/visualizing_replay_results.md'
-      },
-      {
-        text: 'Analyzing Replay Results',
-        link: 'documentation/analyzing_replay_results.md'
-      }
-    ],
-  }),
+    theme: hopeTheme({
+        logo: 'https://i.postimg.cc/gJn8Sy0X/replay-bg-logo.png',
+        navbar: ['/', '/documentation/get_started', '/documentation/'],
+        repo: 'https://github.com/gcappon/py_replay_bg',
+        docsDir: 'docs/',
+        docsBranch: 'master',
+        markdown: {
+            highlighter: {
+                type: "shiki",
+                langs: ['python', 'json', 'md', 'bash', 'diff'],
+                themes: {
+                    dark: 'catppuccin-mocha',
+                    light: 'catppuccin-latte'
+                }
+            },
+            math: {
+                type: "mathjax",
+            },
+        },
+        sidebar: [
+            {
+                text: 'Get Started',
+                link: 'documentation/get_started.md'
+            },
+            {
+                text: 'Data Requirements',
+                link: 'documentation/data_requirements.md'
+            },
+            {
+                text: 'The ReplayBG Object',
+                link: 'documentation/replaybg_object.md'
+            },
+            {
+                text: 'Choosing Blueprint',
+                link: 'documentation/choosing_blueprint.md'
+            },
+            {
+                text: 'The CGM error model',
+                link: 'documentation/cgm_model.md'
+            },
+            {
+                text: 'Twinning Procedure',
+                link: 'documentation/twinning_procedure.md'
+            },
+            {
+                text: 'Replaying',
+                link: 'documentation/replaying.md'
+            },
+            {
+                text: 'The _results/_ Folder',
+                link: 'documentation/results_folder.md'
+            },
+            {
+                text: 'Visualizing Replay Results',
+                link: 'documentation/visualizing_replay_results.md'
+            },
+            {
+                text: 'Analyzing Replay Results',
+                link: 'documentation/analyzing_replay_results.md'
+            }
+        ],
+    }),
 
-  bundler: viteBundler(),
+    bundler: viteBundler(),
 
-  plugins: [],
+    plugins: [],
 })

docs/documentation/README.md

@@ -16,6 +16,8 @@ description: Docs index
 
 ### [Choosing Blueprint](./choosing_blueprint.md)
 
+### [The CGM error model](./cgm_model.md)
+
 ### [Twinning Procedure](./twinning_procedure.md)
 
 ### [Replaying](./replaying.md)

docs/documentation/cgm_model.md

@@ -0,0 +1,82 @@
+---
+sidebar: heading
+---
+
+# The CGM error model
+
+The CGM error model is a key component of the ReplayBG framework, and it is designed to simulate the errors associated with
+these devices.
+
+## Model abstract class
+
+The CGM error model is implemented as an abstract class `CGM`, which defines the interface for
+all CGM error models. This class can be extended by specific CGM error model implementations.
+
+The constructor of the abstract class defines the following attributes:
+
+- `ts`: The sample time of the cgm sensor in minutes (default: 5).
+- `max_lifetime`: The maximum lifetime of the CGM sensor in minutes (default: 1440).
+
+The abstract class defines the following methods:
+
+- `connect_new_cgm(connected_at: int=0)`: Connects a new CGM device, setting the time at which it is connected. Used
+  internally to reset state across multiple replays.
+- `add_offset(to_add: float)`: Adjusts the model's internal offset. Rarely needs overriding.
+- `measure(ig: float, past_ig: list[float], t: float)`:  Function that provides a CGM measure from the past interstitial
+  glucose values. It takes as input the interstitial glucose values `ig`, all the past interstitial glucose values
+  `past_ig`, and the relative time `t` in days since the startup of the sensor.
+
+## Default CGM error model
+
+The default CGM error model available in ReplayBG is the Vettoretti19CGM model as published in Vettoretti et al., Sensors,
+2019, which is a factory-calibrated device with a lifespan of 10 days and a sample time of 5 minutes.
+
+Model equations are:
+
+$
+\begin{cases}
+IG_S(t) = (a_0 + a_1 \cdot t + a_2 \cdot t^2) \cdot IG(t) + b_0 \\
+CGM(t) = IG_S(t) + v(t)
+\end{cases}
+$
+
+with $a_0$, $a_1$, $a_2$, and $b_0$ (min$^{-1}$) model coefficients; and $v(t) \sim N(0, \epsilon_v)$ random white
+noise. Particularly, in order to mimic real CGM systems, $CGM(t)$ is saturated between 40 and 400 mg/dL.
+
+## How to choose the CGM error model to use during replays
+
+The replay method of the `ReplayBG` class accepts a `sensor_cgm` parameter, which can be set to the class of the CGM to
+be used.
+```python
+import MyCustomCGM
+results = rbg.replay(..., sensor_cgm=MyCustomCGM)
+```
+If no `sensor_cgm` is provided, the default CGM error model (Vettoretti19CGM) is used.
+
+### Example of use
+
+This example shows how to use a custom CGM error model in a ReplayBG simulation. The example uses a fake CGM model
+that doubles the interstitial glucose values for demonstration purposes. The sensor has a maximum lifetime of 10
+minutes.
+
+The full example is available in the `example/code/` folder in the `fake_CGM.py` and `replay_map_fakecgm.py` files.
+
+```
+class FakeCGM(CGM):
+    def __init__(self):
+        super().__init__()
+        self.max_lifetime = 10
+
+    def measure(self, ig: float, past_ig: list[float], t):
+        return ig * 2
+
+    def connect_new_cgm(self, connected_at=0):
+        super().connect_new_cgm()
+        print(f"New CGM sensor connected at {connected_at}")
+        
+def main():
+  ...
+  # Replay the twin with the same input data used for twinning
+  replay_results = rbg.replay(..., sensor_cgm=FakeCGM)
+  ...
+```

docs/documentation/choosing_blueprint.md

@@ -139,19 +139,9 @@ $
 where $G_{th}$ is the hypoglycemic threshold (set to 60 mg/dl) and $r_1$ (dimensionless) and $r_2$ (dimensionless) 
 are two model parameters, without direct physiological interpretation, set to 1.44 and 0.81, respectively.
 
-#### CGM sensor error
-The subsystem describes the CGM sensor error model of a factory-calibrated CGM sensor with 10-day lifetime.
-Model equations are:
-
-$
-\begin{cases}
-    IG_S(t) = (a_0 + a_1 \cdot t + a_2 \cdot t^2) \cdot IG(t) + b_0 \\
-    CGM(t) = IG_S(t) + v(t) 
-\end{cases}
-$
-
-with $a_0$, $a_1$, $a_2$, and $b_0$ (min$^{-1}$) model coefficients; and $v(t) \sim N(0, \epsilon_v)$ random white 
-noise. Particularly, in order to mimic real CGM systems, $CGM(t)$ is saturated between 40 and 400 mg/dL.
+#### CGM sensor error model
+The subsystem describes the CGM sensor error model. For more details on the CGM error model, please refer to the
+[CGM Error Model](./cgm_model.md) page.
 
 #### Overall model
 
@@ -306,18 +296,8 @@ where $G_{th}$ is the hypoglycemic threshold (set to 60 mg/dl) and $r_1$ (dimens
 are two model parameters, without direct physiological interpretation, set to 1.44 and 0.81, respectively.
 
 #### CGM sensor error
-The subsystem describes the CGM sensor error model of a factory-calibrated CGM sensor with 10-day lifetime.
-Model equations are:
-
-$
-\begin{cases}
-    IG_S(t) = (a_0 + a_1 \cdot t + a_2 \cdot t^2) \cdot IG(t) + b_0 \\
-    CGM(t) = IG_S(t) + v(t) 
-\end{cases}
-$
-
-with $a_0$, $a_1$, $a_2$, and $b_0$ (min$^{-1}$) model coefficients; and $v(t) \sim N(0, \epsilon_v)$ random white 
-noise. Particularly, in order to mimic real CGM systems, $CGM(t)$ is saturated between 40 and 400 mg/dL.
+The subsystem describes the CGM sensor error model. For more details on the CGM error model, please refer to the
+[CGM Error Model](./cgm_model.md) page.
 
 #### Overall model
 

docs/documentation/replaying.md

@@ -48,7 +48,8 @@ rbg.replay(data: pd.DataFrame,
    save_suffix: str = '',
    save_workspace: bool = False,
    n_replay: int = 1000,
-   sensors: list | None = None
+   sensors: list | None = None,
+   sensor_cgm: CGM = Vettoretti19CGM,
 ) -> Dict:
 ```
 ### Input parameters
@@ -117,6 +118,7 @@ in the `results/workspaces` folder or not.
 replay simulations. Ignored if twinning_method is 'map'.
 - `sensors`: , optional, default: `None`: A `list[Sensors]` to be used in each of the replay simulations. Its length 
 must coincide with the selected `n_replay`. Used when working with intervals. If `None` new sensors will be used.
+- `sensor_cgm`, optional, default: `Vettoretti19CGM`: The class of the CGM error model to be used during the replay simulation.
 
 ::: tip REMEMBER
 The total length of the simulation, `simulation_length`, is defined in minutes and determined by ReplayBG automatically 

py_replay_bg.egg-info/SOURCES.txt

@@ -9,7 +9,6 @@ py_replay_bg.egg-info/SOURCES.txt
 py_replay_bg.egg-info/dependency_links.txt
 py_replay_bg.egg-info/requires.txt
 py_replay_bg.egg-info/top_level.txt
-py_replay_bg/analyzer/__init__.py
 py_replay_bg/data/__init__.py
 py_replay_bg/dss/__init__.py
 py_replay_bg/dss/default_dss_handlers.py

py_replay_bg/example/code/fake_CGM.py

@@ -0,0 +1,20 @@
+import numpy as np
+from py_replay_bg.sensors import CGM
+
+
+class FakeCGM(CGM):
+    def __init__(self):
+        super().__init__()
+        self.max_lifetime = 10  # This sensor only lasts 10 minutes
+
+    def measure(self, ig: float, past_ig: list[float], t):
+        # Fake sensor that just doubles the ig value
+        # return ig * 2
+
+        # Fake sensor that averages the last 10 ig values (10 minutes) plus some noise
+        return np.nanmean(past_ig[-10:]) + np.random.uniform(-10,
+                                                             10)
+
+    def connect_new_cgm(self, connected_at=0):
+        super().connect_new_cgm()
+        print(f"New CGM sensor connected at {connected_at}")  # Just to show that this method is called

py_replay_bg/example/code/replay_map.py

@@ -1,6 +1,7 @@
 import os
 import numpy as np
 
+from py_replay_bg.sensors.Vettoretti19CGM import Vettoretti19CGM
 from utils import load_test_data, load_patient_info
 
 from py_replay_bg.py_replay_bg import ReplayBG
@@ -13,7 +14,7 @@
 
 # Set other parameters for twinning
 blueprint = 'multi-meal'
-save_folder = os.path.join(os.path.abspath(''),'..','..','..')
+save_folder = os.path.join(os.path.abspath(''), '..', '..', '..')
 
 # load patient_info
 patient_info = load_patient_info()

py_replay_bg/example/code/replay_map_fakecgm.py

@@ -0,0 +1,49 @@
+import os
+import numpy as np
+
+from py_replay_bg.example.code.fake_CGM import FakeCGM
+from py_replay_bg.sensors.Vettoretti19CGM import Vettoretti19CGM
+from utils import load_test_data, load_patient_info
+
+from py_replay_bg.py_replay_bg import ReplayBG
+from py_replay_bg.visualizer import Visualizer
+from py_replay_bg.analyzer import Analyzer
+
+# Set verbosity
+verbose = True
+plot_mode = False
+
+# Set other parameters for twinning
+blueprint = 'multi-meal'
+save_folder = os.path.join(os.path.abspath(''), '..', '..', '..')
+
+# load patient_info
+patient_info = load_patient_info()
+p = np.where(patient_info['patient'] == 1)[0][0]
+# Set bw and u2ss
+bw = float(patient_info.bw.values[p])
+
+# Instantiate ReplayBG
+rbg = ReplayBG(blueprint=blueprint, save_folder=save_folder,
+               yts=5, exercise=False,
+               seed=1,
+               verbose=verbose, plot_mode=plot_mode)
+
+# Load data and set save_name
+data = load_test_data(day=1)
+save_name = 'data_day_' + str(1)
+
+print("Replaying " + save_name)
+
+# Replay the twin with the same input data used for twinning
+replay_results = rbg.replay(data=data, bw=bw, save_name=save_name,
+                            twinning_method='map',
+                            save_workspace=True,
+                            save_suffix='_replay_map', sensor_cgm=FakeCGM)
+
+# Visualize and analyze results
+Visualizer.plot_replay_results(replay_results, data=data)
+analysis = Analyzer.analyze_replay_results(replay_results, data=data)
+print('Mean glucose: %.2f mg/dl' % analysis['median']['glucose']['variability']['mean_glucose'])
+print('TIR: %.2f %%' % analysis['median']['glucose']['time_in_ranges']['time_in_target'])
+print('N Days: %.2f days' % analysis['median']['glucose']['data_quality']['number_days_of_observation'])