implement plotting functionality

Author: staskh

README.md

@@ -101,6 +101,7 @@ IGLU_PYTHON extends beyond the capabilities of the original IGLU-R package by of
 | load_dexcom()     | Load Timeseries from Dexcom device file (CGM reading converted into mg/dL)
 | **PLOT/VISUALISE CGM **
 | plot_daily()      | Plot daily Glucose values for each day |
+| plot_statistics() | Plot median + quantile daily statistics |
 
 # Installation
 

iglu_python/__init__.py

@@ -10,7 +10,7 @@
 from .ea1c import ea1c
 from .episode_calculation import episode_calculation
 from .extension.load_data import load_dexcom, load_libre
-from .extension.plots import plot_daily
+from .extension.plots import plot_daily, plot_statistics
 from .gmi import gmi
 from .grade import grade
 from .grade_eugly import grade_eugly
@@ -87,6 +87,7 @@
     "modd",
     "pgs",
     "plot_daily",
+    "plot_statistics",
     "process_data",
     "quantile_glu",
     "range_glu",

iglu_python/extension/plots.py

@@ -3,6 +3,7 @@
 """
 
 import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
 
 
@@ -38,7 +39,7 @@ def plot_daily(cgm_timeseries: pd.Series, lower: int = 70, upper: int = 140) ->
         axes[i].plot(cgm_one_day.index, cgm_one_day.values)
         axes[i].set_title(f"Day: {day.strftime('%Y-%m-%d')}")
         axes[i].set_ylabel("Glucose (mg/dL)")
-        axes[i].set_ylim(0, max(max(cgm_one_day.values), 300))
+        axes[i].set_ylim(0, max(np.nanmax(cgm_one_day.values), 300))
 
         # Fill area above upper limit and plot it in orange
         upper_array = [upper] * len(cgm_one_day.values)
@@ -66,3 +67,82 @@ def plot_daily(cgm_timeseries: pd.Series, lower: int = 70, upper: int = 140) ->
 
     fig.tight_layout()
     return fig
+
+def plot_statistics(cgm_timeseries: pd.Series, lower: int = 70, upper: int = 140) -> plt.Figure:
+    """
+    Plot statistical representation of daily trends
+    in the single 24h timeline, this will plot mean sample trends, 10%, +25% and 75% and 90% quantiles
+    """
+    # check if cgm_timeseries is a pandas series
+    if not isinstance(cgm_timeseries, pd.Series):
+        raise AttributeError("cgm_timeseries must be a pandas series")
+
+    # check if cgm_timeseries is not a datetime index
+    if not isinstance(cgm_timeseries.index, pd.DatetimeIndex):
+        raise AttributeError("cgm_timeseries must have a datetime index")
+
+    # check if cgm_timeseries is not empty
+    if len(cgm_timeseries) < 16:
+        raise ValueError("cgm_timeseries is too short to plot statistics")
+
+    # get sampling frequency
+    time_diffs = cgm_timeseries.index.to_series().diff()
+    dt0 = int(time_diffs.mode().iloc[0].total_seconds() / 60)
+
+
+    # Create time grid
+    start_time = cgm_timeseries.index.min().floor("D")
+    end_time = cgm_timeseries.index.max().ceil("D")
+    time_grid = pd.date_range(start=start_time, end=end_time, freq=f"{dt0}min")
+    # remove the last time point
+    time_grid = time_grid[:-1]
+
+    # interpolate
+    cgm_timeseries_interpolated = np.interp(
+        (time_grid - start_time).total_seconds() / 60,
+        (cgm_timeseries.index - start_time).total_seconds() / 60,
+        cgm_timeseries.values,
+        left=np.nan,
+        right=np.nan,
+    )
+
+    # reorganise as 2d array with rows as timepoints and columns as days
+    # Reshape to days
+    n_days = (end_time - start_time).days
+    n_points_per_day = 24 * 60 // dt0
+    cgm_timeseries_2d = cgm_timeseries_interpolated.reshape(n_days, n_points_per_day)
+
+    # one day time grid
+    time_grid_one_day = time_grid[0:n_points_per_day]
+    # get mean sample trends
+    mean_sample_trends = np.nanmean(cgm_timeseries_2d, axis=0)
+
+    # get 10%, +25% and 75% and 90% quantiles
+    quantiles = np.nanpercentile(cgm_timeseries_2d, [10, 25, 75, 90], axis=0)
+
+    # create figure and axes
+    fig, ax = plt.subplots(figsize=(12, 6))
+
+    # plot mean sample trends
+    ax.plot(time_grid_one_day,mean_sample_trends, color="orange", alpha=1, linewidth=3, label='Mean sample trends')
+
+    # plot quantiles
+    ax.fill_between(time_grid_one_day, quantiles[0], quantiles[1], alpha=0.25, color="blue",label='10% quantile')
+    ax.fill_between(time_grid_one_day, quantiles[1], mean_sample_trends, alpha=0.50, color="blue",label='25% quantile')
+    ax.fill_between(time_grid_one_day, mean_sample_trends, quantiles[2], alpha=0.50, color="blue",label='75% quantile')
+    ax.fill_between(time_grid_one_day, quantiles[2], quantiles[3], alpha=0.25, color="blue",label='90% quantile')
+
+    ax.axhline(y=upper, color="orange", linestyle="--", alpha=0.7, label=f"Hyper threshold ({upper} mg/dL)")
+    ax.axhline(y=lower, color="green", linestyle="--", alpha=0.7, label=f"Hypo threshold ({lower} mg/dL)")
+
+    ax.set_ylim(min(30, np.nanmin(cgm_timeseries.values)), max(np.nanmax(cgm_timeseries.values), 300))
+    ax.set_xlabel("Time (hours)")
+    time_grid_one_day = pd.date_range(start=start_time, periods=24, freq="1h")
+    ax.set_xticks(time_grid_one_day)  # Show every hour from 0 to 24
+    ax.set_xticklabels([f"{h.hour}" for h in time_grid_one_day])  # Format as HH:00
+    ax.grid(True, alpha=0.3, linestyle="--")
+    ax.legend()
+    fig.tight_layout()
+
+    # plot the results
+    return fig