PFun CMA Model -- Modular Demos¶

In [1]:
import pfun_cma_model as pfun

Instantiate the model¶

In [2]:
import json
cma = pfun.CMASleepWakeModel()
cma.to_dict()
Out[2]:
{'_DEFAULT_PARAMS_MODEL': {'N': 1024,
  'd': 0.0,
  'taup': 1.0,
  'taug': 1.0,
  'B': 0.05,
  'Cm': 0.0,
  'toff': 0.0,
  'tM': [7.0, 11.0, 17.5],
  'seed': None,
  'eps': 1e-18},
 'bounds': {'lb': [-12.0, 0.5, 0.1, 0.0, 0.0, -3.0],
  'ub': [14.0, 3.0, 3.0, 1.0, 2.0, 3.0],
  'keep_feasible': [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]},
 'rng': None,
 '_DEFAULT_PARAMS': {'N': 1024,
  'd': 0.0,
  'taup': 1.0,
  'taug': 1.0,
  'B': 0.05,
  'Cm': 0.0,
  'toff': 0.0,
  'tM': [7.0, 11.0, 17.5],
  'seed': None,
  'eps': 1e-18},
 'param_keys': ['N',
  'd',
  'taup',
  'taug',
  'B',
  'Cm',
  'toff',
  'tM',
  'seed',
  'eps']}
In [3]:
cma.params
Out[3]:
CMAModelParams(N=1024, d=0.0, taup=1.0, taug=1.0, B=0.05, Cm=0.0, toff=0.0, tM=[7.0, 11.0, 17.5], seed=None, eps=1e-18)

Run the model (generate 24 hour-long bins)¶

In [4]:
cma.run()
Out[4]:
t c m a I_S I_E L g_0 g_1 g_2 G is_meal
0 days 00:00:00 0.000000 0.083045 0.853553 0.250723 0.152953 0.038349 1.106922e-11 0.100000 0.100000 0.100000 0.300000 False
0 days 00:01:24.457478005 0.023460 0.081543 0.855718 0.250761 0.151199 0.037915 1.354846e-11 0.100000 0.100000 0.100000 0.300000 False
0 days 00:02:48.914956010 0.046921 0.080067 0.857870 0.250800 0.149451 0.037482 1.656335e-11 0.100000 0.100000 0.100000 0.300000 False
0 days 00:04:13.372434019 0.070381 0.078618 0.860007 0.250840 0.147711 0.037052 2.022528e-11 0.100000 0.100000 0.100000 0.300000 False
0 days 00:05:37.829912024 0.093842 0.077194 0.862132 0.250881 0.145977 0.036623 2.466780e-11 0.100000 0.100000 0.100000 0.300000 False
... ... ... ... ... ... ... ... ... ... ... ... ...
0 days 23:54:22.170087975 23.906158 0.021855 0.844762 0.249999 0.175555 0.043888 2.466780e-11 0.100005 0.100032 0.101864 0.301901 False
0 days 23:55:46.627565980 23.929619 0.021434 0.846979 0.249999 0.173500 0.043375 2.022528e-11 0.100005 0.100032 0.101830 0.301867 False
0 days 23:57:11.085043989 23.953079 0.021021 0.849184 0.249999 0.171457 0.042864 1.656335e-11 0.100005 0.100031 0.101797 0.301834 False
0 days 23:58:35.542521994 23.976540 0.020616 0.851375 0.249999 0.169424 0.042356 1.354846e-11 0.100005 0.100031 0.101765 0.301801 False
1 days 00:00:00 24.000000 0.020218 0.853553 0.249999 0.167403 0.041851 1.106922e-11 0.100005 0.100031 0.101734 0.301769 False

1024 rows × 12 columns

In [5]:
cma.run().plot(x='t', y='G', kind='line')
Out[5]:
<Axes: xlabel='t'>
No description has been provided for this image

Load & Process an example dataset¶

In [6]:
from pfun_cma_model.engine.fit import fit_model
from pfun_cma_model.engine.data_utils import format_data
from pfun_cma_model.misc.pathdefs import PFunDataPaths

df_original = PFunDataPaths().read_sample_data()

df = df_original.copy()

print('(raw) sample data')
print(df.tail().to_markdown())
print()

print('formatted sample data')
df_formatted = format_data(df)
print(df_formatted.tail().to_markdown())
print()

ax = df_formatted.plot(
    x='time',
    y=['G'],
    title='Sample Glucose Dataset',
    figsize=(10, 5),
    linestyle='',
    marker='o',
)

(raw) sample data
|      |   Unnamed: 0.1 |   Unnamed: 0 |   user_id | ts_utc              | ts_local            |   is_sg |   sg |   is_fl |   meal_tag | tag_bef_meal   | tag_after_meal   | tag_after_snack   | tag_seems_high   | tag_seems_low   |
|-----:|---------------:|-------------:|----------:|:--------------------|:--------------------|--------:|-----:|--------:|-----------:|:---------------|:-----------------|:------------------|:-----------------|:----------------|
| 2874 |           2874 |         2874 |  10130489 | 2021-01-28 07:08:28 | 2021-01-27 23:08:28 |       1 |  140 |       0 |        nan | True           | False            | False             | True             | False           |
| 2875 |           2875 |         2875 |  10130489 | 2021-01-28 07:13:28 | 2021-01-27 23:13:28 |       1 |  142 |       0 |        nan | True           | False            | False             | True             | False           |
| 2876 |           2876 |         2876 |  10130489 | 2021-01-28 07:18:28 | 2021-01-27 23:18:28 |       1 |  142 |       0 |        nan | True           | False            | False             | True             | False           |
| 2877 |           2877 |         2877 |  10130489 | 2021-01-28 07:23:28 | 2021-01-27 23:23:28 |       1 |  142 |       0 |        nan | True           | False            | False             | True             | False           |
| 2878 |           2878 |         2878 |  10130489 | 2021-01-28 07:28:28 | 2021-01-27 23:28:28 |       1 |  121 |       0 |        nan | True           | False            | False             | False            | False           |

formatted sample data
|       t | time                                |   value |     tod |       t |       G |
|--------:|:------------------------------------|--------:|--------:|--------:|--------:|
| 23.8043 | 2021-01-21 23:43:09.470185456+00:00 | 126.333 | 23.8044 | 23.8043 | 1.37625 |
| 23.8066 | 2021-01-25 23:40:15.322580100+00:00 | 126     | 23.8066 | 23.8066 | 1.36254 |
| 23.8866 | 2021-01-20 23:47:22.102639092+00:00 | 115.333 | 23.8867 | 23.8866 | 1.13052 |
| 23.8891 | 2021-01-24 23:44:27.955033736+00:00 | 123.333 | 23.8892 | 23.8891 | 1.29308 |
| 23.93   | 2021-01-23 23:48:40.587487372+00:00 | 117.5   | 23.9301 | 23.93   | 1.1675  |
No description has been provided for this image

Zoom on Last day (raw data)¶

In [7]:
ax.set_xlim('2021-01-27 16:00:00', df_formatted['time'].max())
ax.get_figure()
Out[7]:
No description has been provided for this image

Fit the model to sample data¶

Raw fit result (solution)¶

In [49]:
print('model fit result')
result = fit_model(data=df_formatted, tcol='time', tm_freq='45min')
print(result.soln.shape)
print(result)
print()

model fit result
(1024, 13)
soln=                                   t         c             m         a  \
0 days 00:00:00             0.000000  0.000936  3.181676e-17  0.250121   
0 days 00:01:24.457478005   0.023460  0.000936  3.166205e-17  0.250123   
0 days 00:02:48.914956010   0.046921  0.000935  3.150475e-17  0.250126   
0 days 00:04:13.372434019   0.070381  0.000935  3.134493e-17  0.250129   
0 days 00:05:37.829912024   0.093842  0.000934  3.118266e-17  0.250132   
...                              ...       ...           ...       ...   
0 days 23:54:22.170087975  23.906158  0.000304  3.292593e-22  0.474829   
0 days 23:55:46.627565980  23.929619  0.000303  3.475400e-22  0.473340   
0 days 23:57:11.085043989  23.953079  0.000303  3.667490e-22  0.471854   
0 days 23:58:35.542521994  23.976540  0.000302  3.869292e-22  0.470372   
1 days 00:00:00            24.000000  0.000302  4.081252e-22  0.468895   

                                I_S       I_E         L       g_0       g_1  \
0 days 00:00:00            0.999785  0.250067  0.999995  0.288256  0.288256   
0 days 00:01:24.457478005  0.999785  0.250070  0.999995  0.288256  0.288256   
0 days 00:02:48.914956010  0.999785  0.250072  0.999995  0.288256  0.288256   
0 days 00:04:13.372434019  0.999785  0.250075  0.999995  0.288256  0.288256   
0 days 00:05:37.829912024  0.999785  0.250078  0.999995  0.288256  0.288256   
...                             ...       ...       ...       ...       ...   
0 days 23:54:22.170087975  0.999930  0.474796  1.000000  0.288282  0.288444   
0 days 23:55:46.627565980  0.999930  0.473307  1.000000  0.288282  0.288440   
0 days 23:57:11.085043989  0.999930  0.471821  1.000000  0.288281  0.288436   
0 days 23:58:35.542521994  0.999931  0.470340  1.000000  0.288281  0.288432   
1 days 00:00:00            0.999931  0.468862  1.000000  0.288280  0.288428   

                                g_2       g_3         G  is_meal  
0 days 00:00:00            0.288256  0.288256  1.153024    False  
0 days 00:01:24.457478005  0.288256  0.288256  1.153024    False  
0 days 00:02:48.914956010  0.288256  0.288256  1.153024    False  
0 days 00:04:13.372434019  0.288256  0.288256  1.153024    False  
0 days 00:05:37.829912024  0.288256  0.288256  1.153024    False  
...                             ...       ...       ...      ...  
0 days 23:54:22.170087975  0.290648  0.295124  1.162498    False  
0 days 23:55:46.627565980  0.290578  0.294893  1.162192    False  
0 days 23:57:11.085043989  0.290509  0.294672  1.161898    False  
0 days 23:58:35.542521994  0.290444  0.294458  1.161614    False  
1 days 00:00:00            0.290380  0.294253  1.161341    False  

[1024 rows x 13 columns] formatted_data=                                         time       value        tod  \
t                                                                      
0.021667  2021-01-21 00:01:18.484847970+00:00  113.666667   0.021667   
0.043056  2021-01-24 00:02:36.969695940+00:00  123.333333   0.043333   
0.065278  2021-01-27 00:03:55.454543910+00:00  132.000000   0.065278   
0.091667  2021-01-20 00:05:31.117301709+00:00  120.333333   0.091944   
0.113333  2021-01-23 00:06:49.602149679+00:00  161.333333   0.113611   
...                                       ...         ...        ...   
23.881389 2021-01-22 23:52:53.219940492+00:00  182.333333  23.881389   
23.903056 2021-01-25 23:54:11.704788462+00:00  128.000000  23.903056   
23.929444 2021-01-18 23:55:47.367546261+00:00  112.333333  23.929722   
23.951389 2021-01-21 23:57:05.852394231+00:00  118.333333  23.951389   
23.973056 2021-01-24 23:58:24.337242201+00:00  126.000000  23.973333   

                   t         G  
t                               
0.021667    0.021667  1.105198  
0.043056    0.043056  1.292240  
0.065278    0.065278  1.540036  
0.091667    0.091667  1.222375  
0.113333    0.113333  1.994650  
...              ...       ...  
23.881389  23.881389  1.999995  
23.903056  23.903056  1.419478  
23.929444  23.929444  1.087873  
23.951389  23.951389  1.181691  
23.973056  23.973056  1.362543  

[1024 rows x 5 columns] cma=<pfun_cma_model.engine.cma.CMASleepWakeModel object at 0x7f1b0036f7d0> popt=array([ 5.25189771e+00,  3.25261494e+03,  8.05064622e-01,  1.44128012e-01,
        5.36199106e-06, -2.21182508e-02]) pcov=array([[ 3.20850760e+02,  7.28704681e+04,  4.55994585e-01,
        -1.01618765e-01, -1.26469592e+01, -4.93003338e-01],
       [ 7.28704681e+04,  3.70697577e+10, -4.27225411e+05,
        -1.16332174e+04, -2.24840479e+05, -2.53833986e+05],
       [ 4.55994585e-01, -4.27225411e+05,  6.06027153e+00,
         1.27139962e-01,  2.88119606e+00,  2.92547162e+00],
       [-1.01618765e-01, -1.16332174e+04,  1.27139962e-01,
         6.78302164e-03,  7.95648285e-02,  7.96633438e-02],
       [-1.26469592e+01, -2.24840479e+05,  2.88119606e+00,
         7.95648285e-02,  1.95481373e+00,  1.53938468e+00],
       [-4.93003338e-01, -2.53833986e+05,  2.92547162e+00,
         7.96633438e-02,  1.53938468e+00,  2.73812026e+00]]) infodict={'message': 'Optimization terminated successfully.', 'error': None, 'ier': 0, 'result':   message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH
  success: True
   status: 0
      fun: 71.67492014774733
        x: [ 5.252e+00  3.253e+03  8.051e-01  1.441e-01  5.362e-06
            -2.212e-02]
      nit: 70
      jac: [ 1.450e-04  0.000e+00  1.262e-02  8.066e-02  3.604e-03
             0.000e+00]
     nfev: 686
     njev: 98
 hess_inv: <6x6 LbfgsInvHessProduct with dtype=float64>} mesg='Optimization terminated successfully.' ier=0 popt_named={'d': np.float64(5.251897714211664), 'taup': np.float64(3252.6149427248142), 'taug': np.float64(0.8050646215473204), 'B': np.float64(0.14412801161625127), 'Cm': np.float64(5.361991056700414e-06), 'toff': np.float64(-0.02211825075866171)} cond=np.float64(81671288781510.69) diag=array([3.20850760e+02, 3.70697577e+10, 6.06027153e+00, 6.78302164e-03,
       1.95481373e+00, 2.73812026e+00])

Raw fitted parameters¶

Displays { param_name: type(param_value) } for all fitted parameters in the model result.

In [50]:
import pandas as pd
raw_params = result.cma.params.model_dump()

print('types of raw parameters')
params_types = {k: ( str(type(v)), repr(v)[:100] + '...' ) for k, v in raw_params.items()}

df_params_types = pd.DataFrame.from_dict(
    params_types, orient='index', columns=['type', 'repr_head'])
df_params_types.sort_values(by='type', inplace=False)

types of raw parameters
Out[50]:
type repr_head
seed <class 'NoneType'> None...
d <class 'float'> 5.251897714211664...
taup <class 'float'> 3252.6149427248142...
B <class 'float'> 0.14412801161625127...
Cm <class 'float'> 5.361991056700414e-06...
toff <class 'float'> -0.02211825075866171...
eps <class 'float'> 1e-18...
N <class 'int'> 1024...
tM <class 'list'> [15.721666666666668, 17.971666666666664, 20.22...
taug <class 'numpy.float64'> np.float64(0.8050646215473204)...

CMAModelParams: Pretty Markdown Table¶

Includes qualitative descriptions (e.g., 'Normal', 'Low', 'High')

In [51]:
from pfun_cma_model import CMAModelParams

params = result.cma.params
print(params.generate_markdown_table())

Parameter    Type             Value    Default    Lower Bound    Upper Bound  Description
-----------  ------  --------------  ---------  -------------  -------------  -------------------------------------------------------
d            float      5.2519            0             -12               14  Time zone offset (hours) (High)
taup         float   3252.61              1               0.5              3  Photoperiod length (hours) (Very High)
taug         float      0.805065          1               0.1              3  Glucose response time constant (Normal)
B            float      0.144128          0.05            0                1  Glucose Bias constant (baseline glucose level) (Normal)
Cm           float      5.36199e-06       0               0                2  Cortisol temporal sensitivity coefficient (Normal)
toff         float     -0.0221183         0              -3                3  Solar noon offset (effects of latitude) (Normal)

Fitted Solution Plotted¶

In [52]:
ax = result.soln.plot(
    x='t',
    y=['G'],
    title='Model Fit Result',
    figsize=(10, 5),
    linestyle='--',
    marker='x',
    logy=True
)
result.formatted_data.plot(
    x='t',
    y=['G'],
    ax=ax,
    linestyle='',
    marker='x',
    color='red',
    title='Data',
)
Out[52]:
<Axes: title={'center': 'Data'}, xlabel='t'>
No description has been provided for this image

Zoomed on dinner estimate (fitted)¶

In [56]:
ax = result.soln.iloc[-200:, :].plot(x='t', y=['G'], figsize=(10, 5), linestyle='-', marker='o')
result.formatted_data.plot(
    x='t',
    y=['G'],
    ax=ax,
    linestyle='',
    marker='x',
    color='red',
    title='Model Fit (Dinner Time Detail)',
)
ax.set_xlim(18, 24)
Out[56]:
(18.0, 24.0)
No description has been provided for this image
In [54]:
from pfun_cma_model.engine.cma_plot import CMAPlotConfig

cma_plot = CMAPlotConfig()

cma_plot.plot_model_results(
    result.formatted_data,
    result.soln,
    as_blob=False,
    separate2subplots=False,
)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[54], line 5
      1 from pfun_cma_model.engine.cma_plot import CMAPlotConfig
      3 cma_plot = CMAPlotConfig()
----> 5 cma_plot.plot_model_results(
      6     result.formatted_data,
      7     result.soln,
      8     as_blob=False,
      9     separate2subplots=False,
     10 )

AttributeError: 'CMAPlotConfig' object has no attribute 'plot_model_results'
No description has been provided for this image
In [ ]:
from pfun_cma_model.engine.cma_plot import CMAPlotConfig

cma_plot = CMAPlotConfig()

cma_plot.plot_model_results(
    result.formatted_data,
    result.soln,
    as_blob=False,
    separate2subplots=True,
)

Experiments with N=...¶

In [ ]:
for n in [12, 24, 288, 1024]:
    print(f'Running model with {n} time points')
    soln = pfun.CMASleepWakeModel(N=n).run()
    print(soln.shape)
    print()

CMA Model Parameters (Curve fitting)¶

Experiments¶

./cma-model-tuning.ipynb

Export notebook to HTML¶

./cma-model-modular-demos.html

In [ ]:
%%sh

uv run jupyter nbconvert --to html cma-model-modular-demos.ipynb --output cma-model-modular-demos.html
cp cma-model-modular-demos.html ../../pfun_cma_model/static/notebooks/cma-model-modular-demos.html
In [ ]: