#!/usr/bin/env python
# coding: utf-8

# # Model Predictive Control for Closed-Loop Insulin Delivery
# 
# **(Or: how I learned to stop worrying and replaced my daughter's pancreas with TensorFlow)**
# 
# Marius Eriksen ([@marius](http://twitter.com/marius), [marius@monkey.org](mailto:marius@monkey.org)), March 2018.
# 
# People with Type-1 Diabetes do not have a functioning pancreas,
# the organ responsible for regulating insulin levels in your body.
# Insulin is a hormone that regulates the metabolism of glucose (and also fats and proteins);
# it works by signaling, among others, fat and liver cells to absorb glucose.
# Blood glucose needs to remain in a fairly narrow "Goldilocks" zone:
# high levels of blood glucose can have severe long-term health consequences;
# low levels of blood glucose can lead to unconsciousness or even death.
# 
# Without a functioning pancreas,
# people living with Type-1 Diabetes must assume its job.
# Insulin must be provided when blood glucose is high 
# (it is usually injected subcutaneously with a syringe or an [insulin pump](https://www.medtronicdiabetes.com/treatments/insulin-pump-therapy)),
# but not too much or blood glucose levels will go too low.
# (If glucose levels do go low, fast-acting carbohydrates must be quickly consumed
# to counteract the effects of the excess insulin.)
# Typically, people with Type-1 Diabetes deliver insulin through an insulin pump or through subcutaneous injections.
# 
# What's more,
# whereas the pancreas can secrete insulin directly into the blood stream,
# so that cells are signaled with little delay,
# externally administered insulin takes much longer to reach its destination.
# [FiASP](https://www.myfiasp.com/) is the fastest acting subcutaneously delivered insulin on the market today, 
# and while it enters the blood stream quite quickly,
# its maximum effect on blood glucose is between 1 and 3 hours after injection.
# This hints at yet another complication: the pharmacokinetics of insulin analogs like used to treat Type-1 Diabetes.
# After a unit of insulin is delivered, it trickles into the blood stream over the course of many hours.
# The following shows how the concentration of FiASP 
# evolves in blood for the first two hours after an injection.
# 
# <img src="https://s3-us-west-2.amazonaws.com/tinyap/fiasp.png" alt="FiASP pharmacokinetic curve" style="width: 400px;"/>
# 
# Carbohydrates consumed in food, by the way, has its own kinetics.
# Depending on the amount of carbohydrates and the food's glycaemic index,
# a meal can take anywhere between an hour or six to fully digest:
# the carbohydrates in the meal enter the blood stream over the same course of time.
# 
# And these are far from the only factors that affect blood glucose!
# Stress levels, exercise, diurnal rhythms, and illness 
# all have a significant effect on how your body metabolizes sugar,
# or how your cells respond to insulin.
# So too do imperfections in insulin manufacturing,
# the site at which insulin is injected,
# the temperature at which insulin is kept,
# and even the mechanical design of an insulin pump.
# 
# It's easy to see how Type-1 (T1D) diabetes is a therapy intense disorder.
# People living with T1D must constantly adjust insulin delivery (or drink some juice),
# depending on any of the above factors.
# People with T1D must also anticipate meals 
# and estimate the amount of carbohydrates that they will eat.
# Ideally, they inject insulin in good time before starting a meal
# so that the insulin has time to enter the blood stream in sufficient
# quantities.
# (But don't wait too long, or your blood glucose will go too low!)
# 
# Closed-loop insulin delivery systems, 
# or "artificial pancreases"
# are designed to assume much of this burden.
# Such systems automate insulin delivery 
# in order to maintain "Goldilocks-zone" blood glucose levels.
# These systems read glucose levels from a
# [continuous glucose monitor](https://www.dexcom.com/)
# and make regular insulin deliveries through an [insulin pump](https://www.animas.com/diabetes-insulin-pump-and-bloog-glucose-meter/onetouch-ping-blood-glucose-monitor)
# based on current and predicted glucose levels.
# 
# We describe a model for predicting future blood glucose values,
# and a technique for using this model to determine how much insulin should be delivered
# to maintain a good blood glucose range.
# 
# Our initial model will be quite simple,
# but the framework we set up is flexible:
# the model can be improved or even replaced later.
# A follow-up notebook will explore some additional ideas
# to improve the model.

# # A simple model for predicting blood glucose
# 
# Our first model is going to be a linear approximation of 
# [Bergman's minimal model](https://www.ncbi.nlm.nih.gov/pubmed/443421).
# Bergman's is a compartmental model, which presents some difficulty in 
# parameter estimation ("training").
# Instead, we're going to replace Bergman's insulin and glucose compartments
# (i.e., the parts of the model that deal with the pharmacokinetics of insulin and glucose)
# with commonplace approximations 
# (which are themselves models: later revisions of our model will incorporate these directly).
# 
# We begin by assuming that the _change_ in blood glucose over a 5-minute period
# is linearly related to 
# (1) the amount of insulin absorbed into blood; 
# (2) the amount of carbohydates that are absorbed into blood;
# and (3) the amount of endogenous glucose production —
# i.e., glucose that is secreted into the blood by your liver.
# 
# We can write this as:
# 
# $$
# \delta_t = c_e + c_C*C_t - c_I*I_t
# $$
# 
# where:
# - $C_t$: the amount of carbohydates absorbed in period $t$;
# - $I_t$: the amount of insulin absorbed in period $t$; and
# - $c_e$: the amount of endogenous production (independent of $t$).
# 
# It turns out that this model is _too simple_.
# In order to compensate for this, 
# we observe that blood glucose typically moves consistently:
# that is, $\delta_t$, the change in blood glucose at time $t$ 
# is highly dependent on the change in blood glucose at time $t-1$, $\delta_{t-1}$.
# 
# So, we'll add this as another term in our model, i.e.:
# 
# $$
# \delta_t = c_e + c_C*C_t - c_I*I_t + c_{\delta}*\delta_{t-1} 
# $$
# 
# In order to make it simpler to  computing with this model,
# we sample all of our data at 5-minute intervals.
# Every 5 minutes, we have a row of data that includes:
# the amount of insulin delivered in the 5-minute period,
# the blood glucose change in the 5-minute period,
# the amount of carbohydrates delivered in the 5-minute period,
# and so on.
# 
# Thus when we talk about time,
# we are talking about a single period, or row, of data.
# 
# The coefficients ($c_e$, $c_C$, $c_I$, and $c_{\delta}$)
# are learned (as we'll see) by the model to weigh
# the relative importance of the various inputs.
# 
# This model is too simple to be useful in practice,
# but it's a good starting point that we can build on.

# ## Building the simple model

# Let's implement this model with actual code,
# and see how it does.
# 
# First, we import some Python packages that can read and clean up raw data
# that comes from our artificial pancreas system (which already samples
# data every 5 minutes).
# We also import some other utilities.

# In[1]:


from IPython.display import HTML, display
import tabulate

import numpy as np
import sys
sys.path.extend(["/Users/marius/src/tinyap.org/model"])
from tapmodel import data
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
import dateutil
tz = dateutil.tz.gettz('US/Pacific')


# We then read some recent data,
# and sample a few rows to show what's there.
# 
# The columns are:
# 
# - `deltaipid`: pump-delivered insulin (mU/5m)
# - `sgv`: sensor glucose value (mg/dL)
# - `ubi`: user bolus input (externally administered insulin, as we often do with FiASP for meals) (mU/5m)
# - `uciXS`, `uciS`, `uciM`, `uciL`: user carb input (entered carb values and their estimated "t-shirt sized" glycaemic index) (g/5m)
# 

# In[2]:


frame = data.read_raw_csv("/Users/marius/Documents/tinyAP/frame.2018-02-10.csv", tz)
frame = frame[['sgv', 'deltaipid', 'ubi', 'uciXS', 'uciS', 'uciM', 'uciL']]
frame[np.isfinite(frame['sgv'])].tail(5)


# It's useful to look at glucose values as _deltas_:
# i.e., the change in blood glucose over a 5-minute period.
# In fact, this is what our model is going to be predicting.
# We compute these and add them as new columns to the frame.
# 
# We also compute a rolling mean of glucose changes;
# this is used to smooth the effect of glucose changes.

# In[3]:


# Maxdelta is the maximum allowable glucose delta.
# We clip deltas greater than maxdelta or smaller than -maxdelta.
maxdelta = 15

frame['deltasgv'] = frame['sgv'] - frame['sgv'].shift(1)
frame['deltasgv'] = np.maximum(np.minimum(frame['deltasgv'], maxdelta), -maxdelta)

frame['meandeltasgv'] = frame['deltasgv'].rolling(window=4, min_periods=1).mean()
frame['prevmeandeltasgv'] = frame['meandeltasgv'].shift(1)
frame[np.isfinite(frame['sgv'])].tail(5)


# The model we proposed above is parameterized in terms of _active_ insulin and glucose, 
# so before we proceed we have to convert glucose and insulin deliveries 
# (e.g., at time $t$) into _active_ glucose and insulin values
# (at times greater than $t$).
# 
# This requires us to model the pharmacodynamics of insulin.
# For the time being we'll use the model worked by [Dragan Maksimovic for Loop](https://github.com/ps2/LoopIOB/blob/master/ScalableExp.ipynb)
# with some standard parameters;
# we'll see later how we can learn its parameters from our data.
# 
# Since our data contains two types of insulin deliveries
# (pump deliveries are Humalog, externally administerd insulin is FiASP),
# we'll use two different sets of parameters.
# 
# Once we're done we have a new column, `ia`, in our frame
# that computes insulin activity (insulin absorption per 5 minutes)
# per our chosen pharmacodynamic model.

# In[4]:


def expia(t, tp, td):
    """Exponential insulin curve, parameterized by peak and duration, 
    due to Dragan Maksimovic (@dm61).
    Worked by Pete (@ps2) in the notebook
    	https://github.com/ps2/LoopIOB/blob/master/ScalableExp.ipynb
    Reworked here to be vectorized by numpy.
    """
    tau = tp*(1-tp/td)/(1-2*tp/td)
    a = 2*(tau/td)
    S = 1/(1-a+(1+a)*np.exp(-td/tau))
    return np.maximum(0., (S/np.power(tau,2))*t*(1-t/td)*np.exp(-t/tau))

# W is the window of history that we consider. 
# Our periods are 5 minutes, so 12*6 = 6 hours.
W = 12*6
Wtime = np.linspace(0., 1.*W, W, endpoint=False, dtype='float32')
humalog_coeffs = expia(Wtime, 13., 48.)  # 1:05h peak, 4h duration
fiasp_coeffs = expia(Wtime, 7., 45.)  # 35m peak, 3:45h duration

# We flip them because we're going to be applying these to "trailing" data.
humalog_coeffs = np.flip(humalog_coeffs, 0)
fiasp_coeffs = np.flip(fiasp_coeffs, 0)

# Ia indicates insulin activity. This is computed by adding up the 
# contributions of each delivery over a rolling window. Conveniently,
# this is equivalent to taking the dot product of the deliveries in the
# window with the coefficients computed above.
frame['ia'] = frame['deltaipid'].rolling(window=W).apply(lambda pids: np.dot(pids, humalog_coeffs)) + \
    frame['ubi'].rolling(window=W).apply(lambda pids: np.dot(pids, fiasp_coeffs))
frame[np.isfinite(frame['sgv'])].tail(5)


# Now we'll use the same technique to also compute carbohydrate absorption.
# For this, we'll use the classic [Walsh curves](https://www.amazon.com/Pumping-Insulin-Everything-Need-Success/dp/1884804845).
# We use a set of parameters for each glycaemic category: 
# these are added together to compute the aggregate carbohydrate absorption, `ca` (g/5m).

# In[5]:


def walshca(t, tdel, tdur):
    """Walsh carb absorption curves with provided delays and duration."""
    return ((t >= tdel) & (t <= tdel+tdur/2)) * (4*(t-tdel)/np.square(tdur)) + \
        ((t > tdel+tdur/2) & (t <= tdel+tdur)) * (4/tdur * (1-(t-tdel)/tdur))

# Carbohydrate parameters (delay, duration): XS, S, M, L.
carb_params = {
    'uciXS': (2., 6.),
    'uciS': (3., 12.), 
    'uciM': (3., 24.),
    'uciL': (3., 36.),
}
carb_coeffs = {which: walshca(Wtime, tdel, tdur) for (which, (tdel, tdur)) in carb_params.items()}
carb_coeffs = {which: np.flip(coeffs, 0) for (which, coeffs) in carb_coeffs.items()}

frame['ca'] = np.zeros(len(frame))
# Add up the contributions from each carbohydrate type.
for column, coeffs in carb_coeffs.items():
    frame['ca'] += frame[column].rolling(window=W).apply(lambda ucis: np.dot(ucis, coeffs))

frame[np.isfinite(frame['sgv'])].tail(5)


# We can now use [Tensorflow](https://www.tensorflow.org)
# to define and fit a model to these data.
# While Tensorflow is overkill here,
# we'll use some of its more advanced features soon.

# In[6]:


import tensorflow as tf
import os

# Quelch annoying Tensorflow warnings.
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"

class model:
    # The following are the columns we're using in our model.
    # The use of 'None' here indicates that the value of the dimension
    # is unknown until runtime. In our case, it's the number of samples
    # that we run through each iteration.
    car = tf.placeholder(tf.float32, [None, 1], name='car')
    ia = tf.placeholder(tf.float32, [None, 1], name='ia')
    prevdelta = tf.placeholder(tf.float32, [None, 1], name='prevdelta')

    # And the following are the (learned) parameters:
    egp_effect = tf.Variable(0., name='egp_effect')
    carb_coeff = tf.Variable(0., name='carb_coeff')
    insulin_coeff = tf.Variable(0., name='insulin_coeff')
    prevdelta_coeff = tf.Variable(0., name='delta_coeff')

    # Our model, as defined above:
    predicted_delta = egp_effect + carb_coeff*car - insulin_coeff*ia + prevdelta_coeff*prevdelta


# Now we're ready to train: 
# we need to fit the coefficients (variables)
# to values that best fit our observed data.
# 
# We'll do this with one of the optimizers that are bundled with Tensorflow. 
# These perform [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) on a _loss function_
# in order to find the values of the variables that minimize the loss.
# We'll begin by using the standard [mean squared error](https://en.wikipedia.org/wiki/Mean_squared_error) loss function.

# In[7]:


actual_delta = tf.placeholder(tf.float32, [None, 1], name='delta')
# m is the batch size. We have to determine this dynamically
# since our input sizes may vary.
m = tf.cast(tf.shape(actual_delta)[0], tf.float32)
loss = tf.reduce_sum(tf.square(actual_delta - model.predicted_delta)) / (2*m)


# Then we train the model using a gradient descent optimizer.
# First we have to create a subset of our dataset that contain only defined values.
# (We'll later apply other filters to the data as well,
# to clean it up in various ways.)

# In[8]:


# Select all the rows that are defined for the columns we use in our model.
# Note that this can make our data frame noncontiguous, but that's ok since
# we've already computed the values (ca and ia) that depend on contiguous data.
rows = np.isfinite(frame['meandeltasgv']) \
    & np.isfinite(frame['prevmeandeltasgv']) \
    & np.isfinite(frame['ca']) \
    & np.isfinite(frame['ia'])

filtered = frame[rows]
filtered.tail(5)


# We now shuffle the data and split it into a training and test set.
# The training set will be used to perform model training (parameter estimation);
# the test set is used to measure the performance of the model.
# We reserve 5% of the data (randomly) sampled for testing.

# In[9]:


shuffled = filtered.sample(frac=1)
n = int(len(shuffled)*0.95)
train, test = shuffled[:n], shuffled[n:]


# In[10]:


def train_feed_dict(batch):
    """Construct a data dictionary from a batch (a pandas frame)."""
    return {
        model.car: batch['ca'].values.reshape([-1, 1]),
        model.ia: batch['ia'].values.reshape([-1, 1]),
        model.prevdelta: batch['prevmeandeltasgv'].values.reshape([-1, 1]),
        actual_delta: batch['meandeltasgv'].values.reshape([-1, 1]),
    }

# The number of epochs to run through. Each epoch sees the entirety
# of the (training) data set.
nepoch = 1000
learning_step = tf.train.AdagradOptimizer(1.0).minimize(loss)

# batch_size is the number of samples we run through each mini-batch
# when training.
batch_size = 5000

nbatch = int(len(train)/batch_size)
if nbatch < 1:
    raise ValueError('not enough data')

feed_dicts = []
for i in range(nbatch):
    batch = train[i*batch_size:(i+1)*batch_size]
    feed_dicts.append(train_feed_dict(batch))

test_feed_dict = train_feed_dict(test)

# Tensorflow needs a session within which to run our 
# computation graphs.
sess = tf.InteractiveSession()

sess.run(tf.global_variables_initializer())
for e in range(nepoch):
    for i in range(nbatch):
        sess.run(learning_step, feed_dict=feed_dicts[i])
    # Print the loss periodically so we can monitor progress.
    if e%500 == 0:
        est_loss = sess.run(loss, test_feed_dict)
        print('epoch {} loss {}'.format(e, est_loss))

est_loss = sess.run(loss, test_feed_dict)
print('final loss: {}'.format(est_loss))
est_egp_effect, est_carb_coeff, est_insulin_coeff, est_prevdelta_coeff = \
    sess.run([model.egp_effect, model.carb_coeff, model.insulin_coeff, model.prevdelta_coeff])
tab = [
    ["EGP",       est_egp_effect],
    ["Carb",      est_carb_coeff],
    ["Insulin",   est_insulin_coeff],             
    ["Prevdelta", est_prevdelta_coeff],
]
display(HTML(tabulate.tabulate(tab, tablefmt='html')))


# We can interpret our estimated coeffients as follows: 
# 
# 1. The EGP and insulin coefficients can give a _basal rate_.
# 2. The carb and insulin coefficients can give a _carb ratio_.
# 3. Prevdelta as the amount which blood glucose tends to move in the direction it's already moving.
#   Since this is generally true, you can also think of this as glucose changes that probably should
#   be explained by something that isn't modeled explicitly in the model.
# 
# The estimated loss of $1.7$ is the mean squared error we expect from our model.
# In our case, this means that the mean square of the error in our predictions is $1.7$ every five minutes.

# ## Predicting future blood glucose values

# Now we can use our model to visualize how it fares in a real-world scenario
# by feeding the model past data (insulin and carb history for the last 6 hours),
# and then asking for its prediction.
# 
# It's not terribly meaningful to infer _too_ far into the future,
# since this prediction doesn't account for _future_ carb or insulin deliveries
# (including what the AP controller will deliver; more on this soon).

# In[11]:


def plot_glucose(time, actual, predicted):
    """Plot actual vs. predicted glucose."""
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.set_ylabel('blood glucose (mg/dL)')
    ax.set_xlabel('time')
    ax.plot(time, actual, label='actual')
    ax.plot(time, predicted, label='predicted')
    ax.xaxis.set_major_locator(plt.MaxNLocator(10))
    fmt = DateFormatter('%H:%M', tz=tz)
    ax.xaxis.set_major_formatter(fmt)
    plt.legend()
    plt.show()


# ## Predicting a single time step
# 
# The following illustrates how our model performs 
# when predicting a single time step: the next glucose delta.
# 
# As expected, it's generally pretty good at this.
# We can see that it follows the trends quite well,
# though it may miss some of the "spikyness" of the actual data. 
# (This isn't a big problem in practice.)
# 
# The first plot shows the _deltas_ predicted directly,
# the second shows these deltas accumulated over a 3 hour period
# (and thus errors are also accumulated).
# Again, this isn't necessarily a super meaningful way to visualize
# how the model is doing, 
# but it gives you a good idea of cumulative error.
# 
# You'll also note that the timeframe chosen was during the night time,
# which is easier for any model:
# it is a more controlled environment as
# there is (usually) little to no carbohydrate activity,
# and also little physical activity.

# In[12]:


testframe = frame[(frame.index >= '2018-02-06 00:00 PST') & (frame.index <= '2018-02-06 03:00 PST')]
predicted_deltas = sess.run(model.predicted_delta, feed_dict=train_feed_dict(testframe))
plot_glucose(testframe.index, testframe.deltasgv, predicted_deltas)


# In[13]:


predicted_sgvs = testframe.iloc[0].sgv + np.cumsum(predicted_deltas)
plot_glucose(testframe.index, testframe.sgv, predicted_sgvs)


# We can also see how the model performs if we ask if to infer 3 hours in the future.
# This is not a meaningful thing to do in this context,
# as it does not account for insulin (or carb) deliveries
# during this time period,
# but it's nevertheless a fun visualization.

# In[14]:


def feed_dict_row_delta(row, prevdelta):
    """Construct a data dictionary for inference 
    from a single row and a prevdelta."""
    return {
        model.car: np.array([row['ca']]).reshape([-1, 1]),
        model.ia: np.array([row['ia']]).reshape([-1, 1]),
        model.prevdelta: np.array([prevdelta]).reshape([-1, 1]),
    }

predicted_deltas = np.zeros(len(testframe))
for i in range(len(testframe)):
    row = testframe.iloc[i]
    if i == 0:
        # Start off with the initial mean delta.
        prevdelta = row['prevmeandeltasgv']
    else:
        prevdelta = predicted_deltas[i-1]
    predicted_deltas[i] = \
        sess.run(model.predicted_delta, feed_dict=feed_dict_row_delta(row, prevdelta))
predicted_sgvs = testframe.iloc[0].sgv + np.cumsum(predicted_deltas)
plot_glucose(testframe.index, testframe.sgv, predicted_sgvs)


# ## Visualizing the model by synthesizing scenarios

# Another useful "sanity check" on our model is 
# to see how it responds to certain common scenarios.
# We have certain ideas for how your body should respond
# to, say, eating carbohydrates without insulin,
# or giving a correction bolus.
# 
# The following sets up code to synthesize (future) events,
# and then visualize the model's predictions for
# how the body responds to these events.

# In[15]:


from numpy.lib.stride_tricks import as_strided

def rolldot(events, coeffs, w=W):
    """Helper function to compute a rolling dot-product,
    used to compute current insulin and carbohydrate activity
    according to our curves."""
    events = np.append(np.zeros(w), events)
    stride, = events.strides
    strides = as_strided(events, shape=[len(events) - w + 1, w], strides=[stride, stride])
    strides = strides[1:]
    return np.dot(strides, coeffs)

def feed_dict_synthetic(carbs, boluses, prevdelta):
    """Construct a data dictionary with the 
    provided carb and bolus events. The returned
    feed dict contains only the last entry."""
    return {
        model.car: rolldot(carbs, carb_coeffs['uciM'])[-1].reshape([1, 1]),
        model.ia: rolldot(boluses, humalog_coeffs)[-1].reshape([-1, 1]),
        model.prevdelta: np.array([prevdelta]).reshape([-1, 1]),
    }

def predict_glucose(carbs, boluses, npad=12*3):
    """Infers future glucose values based on the provided
    carbohydrate and bolus (insulin) events."""
    pad = np.zeros(npad)
    carbs = np.append(carbs, pad)
    boluses = np.append(boluses, pad)
    prevdelta = 0
    deltas = np.zeros(len(carbs))
    for i in range(len(carbs)):
        prevdelta = sess.run(
            model.predicted_delta,
            feed_dict=feed_dict_synthetic(carbs[:i+1], boluses[:i+1], prevdelta),
        )
        deltas[i] = prevdelta
    return np.cumsum(deltas)

def plot_predict_glucose(title, events):
    """Plot the model's glucose predictions given 
    a set of events."""
    maxtime = 0
    for (time, _, _) in events:
        if time > maxtime:
            maxtime = time
    carbs, boluses = np.zeros(maxtime+1), np.zeros(maxtime+1)
    for (time, event, v) in events:
        if event == 'C':
            carbs[time] = v
        elif event == 'B':
            boluses[time] = v
        else:
            raise Exception("unknown event type")
    glucose = predict_glucose(carbs, boluses)
    n = len(glucose)
    time = np.linspace(0., 5.*n, n, endpoint=False)
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.set_ylabel('change in blood glucose (mg/dL)')
    ax.set_xlabel('time (m)')
    ax.set_title(title)
    ax.plot(time, glucose)
    plt.show()

plot_predict_glucose('prebolus for food', [(0, 'B', 100), (4, 'B', 2000), (6, 'C', 20)])
plot_predict_glucose('do nothing', [])
plot_predict_glucose('give big bolus', [(0, 'B', 2000)])
plot_predict_glucose('eat carbs then bolus later', [(0, 'C', 10), (4, 'B', 1000)])


# As we can see,
# our simple model follows intuition:
# 
# - when we pre-bolus for food, glucose dips a little before rising
#   and settling at a slightly higher level;
# - if we do nothing (i.e., no insulin delivery, no carb delivery),
#   endogenous glucose production ensures that glucose levels steadily
#   rise;
# - if we give a big bolus, glucose levels drop significantly;
# - if we eat carbs before bolusing, glucose levels rise and then level off.

# # Creating an insulin delivery schedule
# 
# We are now ready to use our blood glucose model to create an insulin delivery controller.
# The controller should use the predictions given by the glucose model 
# to learn a _schedule_ of insulin deliveries that 
# maximizes the amount of time where glucose is in range (near the setpoint).
# 
# Here's how it's going to work:
# given the last 6 hours of insulin and carb deliveries, 
# we use our model to predict the next 6 hours.
# While the previous 6 hours of insulin delivery are fixed immutable,
# we can _train_ the schedule for the next 6 hours.
# We use a loss function that minimizes the (squared) distance from our desired setpoint.

# In[16]:


def roll_indices(n):
    m = np.zeros([n, n], dtype=int)
    for i in range(0, n):
        m[i, :] = np.arange(i, i+n, dtype=int)
    return m

# rolled is a W-by-W matrix of indices used to construct
# lookback windows for carb and insulin curves.
# These cannot be computed outside of Tensorflow anymore because they 
# are part of the computation that is involved in optimization.
rolled = roll_indices(W)

class schedule:
    # future_dpids is the next 6 hours of insulin deliveries, as learned
    # by the schedule model.
    future_dpids = tf.Variable(tf.zeros([W, 1]), name='future_dpids')
    
    # These are placeholders used to feed the glucose model.
    past_dpids = tf.placeholder(tf.float32, [W, 1], name='past_dpids')
    past_ucis = tf.placeholder(tf.float32, [W, 1], name='past_ucis')
    
    # Last_delta is used to initialize the pdgsv factor from the glucose model.
    last_delta = tf.placeholder(tf.float32, [1], name='last_delta')
    # Deviation is the deviation from the desired setpoint: i.e., the amount 
    # (up or down) glucose should move and settle to, relative to the starting
    # point.
    deviation = tf.placeholder(tf.float32, [], name='deviation')
    
    dpids = tf.gather(
        tf.reshape(tf.concat([past_dpids, future_dpids], 0), [-1]), 
        rolled,
    )
    ucis = tf.gather(
        tf.reshape(tf.concat([past_ucis, tf.zeros([W, 1])], 0), [-1]),
        rolled,
    )
    ia = tf.matmul(dpids, humalog_coeffs.reshape([-1, 1]))
    car = tf.matmul(ucis, carb_coeffs['uciM'].reshape([-1, 1]))
    
    egp_effect = model.egp_effect
    carb_effect = car*model.carb_coeff
    insulin_effect = ia*model.insulin_coeff
    
    # Delta contains the predicted deltas, sans the prevdelta term.
    # Since prevdelta uses the previous SGV delta, we must use the 
    # predictions yielded by the model itself to compute this term
    # for future data points. We do this by scanning over the tensor
    # of deltas, accumulating contributions.
    delta = egp_effect + carb_effect - insulin_effect
    delta = tf.scan(lambda a, x: a*model.prevdelta_coeff+x, delta, last_delta)
    cumdelta = tf.scan(lambda a, x: a + x, delta)
    
    # Define a loss function. This is the function we will optimize.
    # Negative insulin deliveries are not permitted (they are not physically
    # possible—it would be nice if they were!), and so we must penalize
    # these heavily so that the model does not find schedules that would
    # subtract insulin. (Interestingly, if it could, it could almost always
    # correct positive deviations immediately, since it could give a large 
    # amount of insulin and then cut the tail.)
    scaleneg = 100000.
    loss = tf.reduce_sum(tf.pow(deviation+cumdelta, 2))/(2*W) \
        - scaleneg*tf.minimum(tf.reduce_min(future_dpids), 0)
        
    optimizer = tf.train.GradientDescentOptimizer(5.0)
    grads, v = zip(*optimizer.compute_gradients(loss, [future_dpids]))
    grads, _ = tf.clip_by_global_norm(grads, 1.)
    #assert len(grads) == 1
    #grads = [tf.clip_by_value(grads[0], -1., 1.)]
    learning_step = optimizer.apply_gradients(zip(grads, v))
      
    clip = future_dpids.assign(tf.maximum(future_dpids, 0.))
    
    sess.run(future_dpids.initializer)


# With a scheduling model in hand,
# we can go through the same exercise as above:
# synthesize a scenario, and see what insulin delivery schedule is learned.
# The big difference here is that the glucose values reflect how the model
# believes glucose is going to respond not just to recent history but to the
# insulin schedule that's been learned by the model itself.
# 
# We'll plot several scenarios: 
# each includes both the predicted glucose schedule 
# as well as the learned insulin delivery schedule.

# In[17]:


def schedule_insulin(ucis, dpids, deviation, niter=500, display_loss=False):
    """Create a 6-hour insulin schedule based on the past carb and insulin inputs,
    correcting for the provided deviation."""
    # Reset the future DPIDS back to their initial value.
    sess.run(schedule.future_dpids.assign(np.zeros([W, 1])))
    # Extend history to go W in the past.
    feed_dict = {
        schedule.past_ucis: ucis.reshape([-1, 1]),
        schedule.past_dpids: dpids.reshape([-1, 1]),
        schedule.last_delta: np.array([0.]).reshape([1]),
        schedule.deviation: deviation,
    }
    # Learn the insulin delivery schedule.
    for i in range(niter):
        sess.run(schedule.learning_step, feed_dict=feed_dict)
        if display_loss and i%100 == 0:
            loss = sess.run(schedule.loss, feed_dict=feed_dict)
            print("iter {} loss {}".format(i, loss))
    dpids, cumdelta = sess.run([schedule.future_dpids, schedule.cumdelta], feed_dict=feed_dict)
    return dpids, cumdelta

def plot_schedule_insulin(title, deviation, events):
    """Plot the learned insulin schedule given the provided deviation and
    past events."""
    ucis, dpids = np.zeros(W), np.zeros(W)
    for (time, event, v) in events:
        if time >= 0:
            raise Exception("bad time")
        if event == 'C':
            ucis[time] = v
        elif event == 'B':
            dpids[time] = v
        else:
            raise Exception("unknown event type")

    dpids, cumdelta = schedule_insulin(ucis, dpids, deviation)
    n = len(dpids)
    time = np.linspace(0., 5.*n, n, endpoint=False)
    fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()
    ax1.plot(time, deviation+cumdelta)
    ax1.set_title(title)
    ax2.bar(time, dpids, width=5)
    ax1.set_xlabel("minutes")
    ax1.set_ylabel("deviation (mg/dL) ")
    ax1.set_ylim([-40, 40])
    ax2.set_ylabel("delivery (mU/5m)")
    ax2.set_ylim([0, 500])
    ax1.legend()
    plt.show()

plot_schedule_insulin('hold steady', 0, [])
plot_schedule_insulin('correct deviation of 30', 30, [])
plot_schedule_insulin('correct deviation of -30', -30, [])
plot_schedule_insulin('hold steady after bolus', 0, [(-1, 'B', 400)])


# We can see here that the scheduling model is able to 
# learn a schedule that minimizes deviation from our desired setpoint (0 on the y axis).
# Since the model minimizes the mean-squared-error,
# this means that it also takes into account cumulative effects.
# For example, in the second scenario, the learned insulin schedule
# frontloads most of the insulin (and, presumably, withholds delivery later)
# in order to quickly bring the blood glucose down to the setpoint.
# Thus, the traditional distinction between "basal" and "bolus" insulin deliveries is gone:
# the model simply delivers insulin that maximizes the amount 
# of time spent near the desired setpoint.
# 
# A closed-loop insulin delivery system computes a schedule at regular intervals
# in order to account for new inputs (insulin delivered, carbs consumed) and for 
# unexplained deviations in the model 
# (i.e., where the actual deviations differ from the expected ones).
# The system then delivers insulin according to the schedule.

# # Epilogue
# 
# In this notebook,
# we have derived a (very) simple model for blood glucose,
# and shown how to make use of this model in an insulin controller,
# such as may be used by a closed-loop insulin delivery system.
# 
# At a very high level,
# this is the strategy that is used in [tinyAP](https://github.com/tinyap),
# but reality demands a much more complicated model,
# to account for time-related effects, glucose-related effects,
# and other nonlinearities.
# 
# The framework set up here is a flexible one:
# any modelled effect is automatically accounted for when computing insulin deliveries.
# For example, if we make EGP a function of time (as is done in the model used by tinyAP),
# time is automatically accounted for:
# if, in an hour, the estimated EGP increases or decreases, 
# the learned schedule will anticipate this in good time and deliver (or withhold) insulin 
# so that the desired setpoint is maintained.
# This, of course, is similar to a traditional "basal schedule",
# but with a key difference:
# Basal schedules are _delivery schedules_, that is, 
# they specify the rate of insulin to be delivered at specific times,
# in anticipation of increased (or decreased) insulin needs in the future.
# With model predictive control, we instead model EGP and other effects directly;
# the insulin scheduler learns a _delivery schedule_ based instead on the directly modelled effects.
# 
# Another advantage of this approach is its composability. 
# Because we have a general framework for relating inputs
# (insulin, carbohydrates, time, glucose levels, etc.) 
# to predicted blood glucose,
# improving the model (in the sense that it accounts for more of the observed glucose level variance)
# automatically improves the AP's performance.
# As an example, the model used by tinyAP includes parameters for both
# time (to account for different insulin sensitivies and endogenous glucose production rates throughout the day) 
# and glucose levels (to account for observed effects of increased insulin resistance at higher glucose levels, presumably due to hepatic activity).
# In order to account for these observed effects, 
# only the model itself was modified;
# the rest of the system remained unchanged.
# 
# In the next installment,
# we'll build on the simple model here to account for some of the aforementioned effects.