# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

"""
Explanation: Modeling and Simulation in Python
Chapter 18
Copyright 2017 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
"""


data = pd.read_csv('data/glucose_insulin.csv', index_col='time');

"""
Explanation: Code from the previous chapter
Read the data.
End of explanation
"""


I = interpolate(data.insulin)

"""
Explanation: Interpolate the insulin data.
End of explanation
"""


G0 = 290
k1 = 0.03
k2 = 0.02
k3 = 1e-05

"""
Explanation: Initialize the parameters
End of explanation
"""


Gb = data.glucose[0]
Ib = data.insulin[0]

"""
Explanation: To estimate basal levels, we'll use the concentrations at t=0.
End of explanation
"""


init = State(G=G0, X=0)

"""
Explanation: Create the initial condtions.
End of explanation
"""


t_0 = get_first_label(data)
t_end = get_last_label(data)

system = System(G0=G0, k1=k1, k2=k2, k3=k3,
                init=init, Gb=Gb, Ib=Ib, I=I,
                t_0=t_0, t_end=t_end, dt=2)

def update_func(state, t, system):
    """Updates the glucose minimal model.
    
    state: State object
    t: time in min
    system: System object
    
    returns: State object
    """
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    dt = system.dt
    
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    G += dGdt * dt
    X += dXdt * dt

    return State(G=G, X=X)

def run_simulation(system, update_func):
    """Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    """
    t_0, t_end, dt = system.t_0, system.t_end, system.dt
    
    frame = TimeFrame(columns=init.index)
    frame.row[t_0] = init
    ts = linrange(t_0, t_end, dt)
    
    for t in ts:
        frame.row[t+dt] = update_func(frame.row[t], t, system)
    
    return frame

results = run_simulation(system, update_func);

"""
Explanation: Make the System object.
End of explanation
"""


def slope_func(state, t, system):
    """Computes derivatives of the glucose minimal model.
    
    state: State object
    t: time in min
    system: System object
    
    returns: derivatives of G and X
    """
    G, X = state
    k1, k2, k3 = system.k1, system.k2, system.k3 
    I, Ib, Gb = system.I, system.Ib, system.Gb
    
    dGdt = -k1 * (G - Gb) - X*G
    dXdt = k3 * (I(t) - Ib) - k2 * X
    
    return dGdt, dXdt

"""
Explanation: Numerical solution
In the previous chapter, we approximated the differential equations with difference equations, and solved them using run_simulation.
In this chapter, we solve the differential equation numerically using run_ode_solver, which is a wrapper for the SciPy ODE solver.
Instead of an update function, we provide a slope function that evaluates the right-hand side of the differential equations.  We don't have to do the update part; the solver does it for us.
End of explanation
"""


slope_func(init, 0, system)

"""
Explanation: We can test the slope function with the initial conditions.
End of explanation
"""


results2, details = run_ode_solver(system, slope_func, t_eval=data.index);

"""
Explanation: Here's how we run the ODE solver.
End of explanation
"""


details

"""
Explanation: details is a ModSimSeries object with information about how the solver worked.
End of explanation
"""


results2

"""
Explanation: results is a TimeFrame with one row for each time step and one column for each state variable:
End of explanation
"""


plot(results.G, '-')
plot(results2.G, '-')
plot(data.glucose, 'bo')

"""
Explanation: Plotting the results from run_simulation and run_ode_solver, we can see that they are not very different.
End of explanation
"""


diff = results.G - results2.G
percent_diff = diff / results2.G * 100
percent_diff.dropna()

"""
Explanation: The differences in G are less than 1%.
End of explanation
"""


params = Params(G0 = 290,
                k1 = 0.03,
                k2 = 0.02,
                k3 = 1e-05)

"""
Explanation: Optimization
Now let's find the parameters that yield the best fit for the data.  
We'll use these values as an initial estimate and iteratively improve them.
End of explanation
"""


def make_system(params, data):
    """Makes a System object with the given parameters.
    
    params: sequence of G0, k1, k2, k3
    data: DataFrame with `glucose` and `insulin`
    
    returns: System object
    """
    # params might be a Params object or an array,
    # so we have to unpack it like this
    G0, k1, k2, k3 = params
    
    Gb = data.glucose[0]
    Ib = data.insulin[0]
    I = interpolate(data.insulin)
    
    t_0 = get_first_label(data)
    t_end = get_last_label(data)

    init = State(G=G0, X=0)
    
    return System(G0=G0, k1=k1, k2=k2, k3=k3,
                  init=init, Gb=Gb, Ib=Ib, I=I,
                  t_0=t_0, t_end=t_end, dt=2)

system = make_system(params, data)

"""
Explanation: make_system takes the parameters and actual data and returns a System object.
End of explanation
"""


system = make_system(params, data)
results, details = run_ode_solver(system, slope_func, t_eval=data.index)
details

def error_func(params, data):
    """Computes an array of errors to be minimized.
    
    params: sequence of parameters
    data: DataFrame of values to be matched
    
    returns: array of errors
    """
    print(params)
    
    # make a System with the given parameters
    system = make_system(params, data)
    
    # solve the ODE
    results, details = run_ode_solver(system, slope_func, t_eval=data.index)
    
    # compute the difference between the model
    # results and actual data
    errors = results.G - data.glucose
    return errors

"""
Explanation: error_func takes the parameters and actual data, makes a System object, and runs odeint, then compares the results to the data.  It returns an array of errors.
End of explanation
"""


error_func(params, data)

"""
Explanation: When we call error_func, we provide a sequence of parameters as a single object.
Here's how that works:
End of explanation
"""


best_params, fit_details = leastsq(error_func, params, data)

"""
Explanation: leastsq is a wrapper for scipy.optimize.leastsq
Here's how we call it.
End of explanation
"""


best_params

"""
Explanation: The first return value is a Params object with the best parameters:
End of explanation
"""


fit_details

"""
Explanation: The second return value is a ModSimSeries object with information about the results.
End of explanation
"""


system = make_system(best_params, data)
results, details = run_ode_solver(system, slope_func, t_eval=data.index)
details.message

"""
Explanation: Now that we have best_params, we can use it to make a System object and run it.
End of explanation
"""


plot(results.G, label='simulation')
plot(data.glucose, 'bo', label='glucose data')

decorate(xlabel='Time (min)',
         ylabel='Concentration (mg/dL)')

savefig('figs/chap08-fig04.pdf')

"""
Explanation: Here are the results, along with the data.  The first few points of the model don't fit the data, but we don't expect them to.
End of explanation
"""


def indices(params):
    """Compute glucose effectiveness and insulin sensitivity.
    
    params: sequence of G0, k1, k2, k3
    data: DataFrame with `glucose` and `insulin`
    
    returns: State object containing S_G and S_I
    """
    G0, k1, k2, k3 = params
    return State(S_G=k1, S_I=k3/k2)

"""
Explanation: Interpreting parameters
Based on the parameters of the model, we can estimate glucose effectiveness and insulin sensitivity.
End of explanation
"""


indices(best_params)

"""
Explanation: Here are the results.
End of explanation
"""


source_code(run_ode_solver)

source_code(leastsq)

"""
Explanation: Under the hood
Here's the source code for run_ode_solver and leastsq, if you'd like to know how they work.
End of explanation
"""


