<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Modeling and Simulation in Python
Chapter 18
Copyright 2017 Allen Downey
License
Step1:  Code from the previous chapter
Read the data.
Step2:  Interpolate the insulin data.
Step3:  The glucose minimal model
I'll cheat by starting with parameters that fit the data roughly; then we'll see how to improve them.
Step5:  Here's a version of make_system that takes the parameters and data
Step7:  And here's the update function.
Step8:  Before running the simulation, it is always a good idea to test the update function using the initial conditions.  In this case we can veryify that the results are at least qualitatively correct.
Step10:  Now run_simulation is pretty much the same as it always is.
Step11:  And here's how we run it.
Step12:  The results are in a TimeFrame object with one column per state variable.
Step13:  The following plot shows the results of the simulation along with the actual glucose data.
Step15:  Numerical solution
Now let's solve the differential equation numerically using run_ode_solver, which is an implementation of Ralston's method.
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.
Step16:  We can test the slope function with the initial conditions.
Step17:  Here's how we run the ODE solver.
Step18:  details is a ModSimSeries object with information about how the solver worked.
Step19:  results is a TimeFrame with one row for each time step and one column for each state variable
Step20:  Plotting the results from run_simulation and run_ode_solver, we can see that they are not very different.
Step21:  The differences in G are less than 2%.
Step22:  Exercises
Exercise
Step23:  Under the hood
Here's the source code for run_ode_solver if you'd like to know how it works.
Notice that run_ode_solver is another name for run_ralston, which implements Ralston's method.
<ASSISTANT_TASK:>
Python Code:
# 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
params = Params(G0 = 290,
                k1 = 0.03,
                k2 = 0.02,
                k3 = 1e-05)
Explanation: The glucose minimal model
I'll cheat by starting with parameters that fit the data roughly; then we'll see how to 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
    
    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(params,
                  init=init, Gb=Gb, Ib=Ib, I=I,
                  t_0=t_0, t_end=t_end, dt=2)
system = make_system(params, data)
Explanation: Here's a version of make_system that takes the parameters and data:
End of explanation
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)
Explanation: And here's the update function.
End of explanation
update_func(system.init, system.t_0, system)
Explanation: Before running the simulation, it is always a good idea to test the update function using the initial conditions.  In this case we can veryify that the results are at least qualitatively correct.
End of explanation
def run_simulation(system, update_func):
    Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    
    init = system.init
    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
Explanation: Now run_simulation is pretty much the same as it always is.
End of explanation
results = run_simulation(system, update_func);
Explanation: And here's how we run it.
End of explanation
results
Explanation: The results are in a TimeFrame object with one column per state variable.
End of explanation
subplot(2, 1, 1)
plot(results.G, 'b-', label='simulation')
plot(data.glucose, 'bo', label='glucose data')
decorate(ylabel='Concentration (mg/dL)')
subplot(2, 1, 2)
plot(results.X, 'C1', label='remote insulin')
decorate(xlabel='Time (min)', 
         ylabel='Concentration (arbitrary units)')
savefig('figs/chap18-fig01.pdf')
Explanation: The following plot shows the results of the simulation along with the actual glucose data.
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
Now let's solve the differential equation numerically using run_ode_solver, which is an implementation of Ralston's method.
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(system.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)
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, 'C0', label='run_simulation')
plot(results2.G, 'C2--', label='run_ode_solver')
decorate(xlabel='Time (min)', ylabel='Concentration (mg/dL)')
savefig('figs/chap18-fig02.pdf')
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
max(abs(percent_diff))
Explanation: The differences in G are less than 2%.
End of explanation
# Solution
system3 = System(system, dt=1)
results3, details = run_ode_solver(system3, slope_func)
details
# Solution
plot(results2.G, 'C2--', label='run_ode_solver (dt=2)')
plot(results3.G, 'C3:', label='run_ode_solver (dt=1)')
decorate(xlabel='Time (m)', ylabel='mg/dL')
# Solution
diff = (results2.G - results3.G).dropna()
percent_diff = diff / results2.G * 100
# Solution
max(abs(percent_diff))
Explanation: Exercises
Exercise:  Our solution to the differential equations is only approximate because we used a finite step size, dt=2 minutes.
If we make the step size smaller, we expect the solution to be more accurate.  Run the simulation with dt=1 and compare the results.  What is the largest relative error between the two solutions?
End of explanation
source_code(run_ode_solver)
Explanation: Under the hood
Here's the source code for run_ode_solver if you'd like to know how it works.
Notice that run_ode_solver is another name for run_ralston, which implements Ralston's method.
End of explanation
<END_TASK>