# 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 *

init = State(y = 2)
system = System(init=init, t_0=1, t_end=3)

def slope_func(state, t, system):
    [y] = state
    dydt = y + t
    return [dydt]

results, details = run_euler(system, slope_func)

get_last_value(results.y)

"""
Explanation: Modeling and Simulation in Python
Experiments with different ODE solvers
Copyright 2019 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
"""


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

"""
Explanation: Glucose minimal model
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

%time results = run_simulation(system, update_func);

results

"""
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_euler...
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
"""


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=1)

%time results2, details = run_euler(system, slope_func)

"""
Explanation: Here's how we run the ODE solver.
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_euler, we can see that they are not very different.
End of explanation
"""


diff = results.G - results2.G
percent_diff = diff / results2.G * 100

max(abs(percent_diff.dropna()))

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


m = UNITS.meter
s = UNITS.second

"""
Explanation: Dropping pennies
I'll start by getting the units we need from Pint.
End of explanation
"""


init = State(y=381 * m, 
             v=0 * m/s)

"""
Explanation: And defining the initial state.
End of explanation
"""


g = 9.8 * m/s**2

"""
Explanation: Acceleration due to gravity is about 9.8 m / s$^2$.
End of explanation
"""


t_end = 10 * s

"""
Explanation: When we call odeint, we need an array of timestamps where we want to compute the solution.
I'll start with a duration of 10 seconds.
End of explanation
"""


system = System(init=init, g=g, t_end=t_end)

"""
Explanation: Now we make a System object.
End of explanation
"""


def slope_func(state, t, system):
    """Compute derivatives of the state.
    
    state: position, velocity
    t: time
    system: System object containing `g`
    
    returns: derivatives of y and v
    """
    y, v = state
    g = system.g    

    dydt = v
    dvdt = -g
    
    return dydt, dvdt

"""
Explanation: And define the slope function.
End of explanation
"""


dydt, dvdt = slope_func(system.init, 0, system)
print(dydt)
print(dvdt)

"""
Explanation: It's always a good idea to test the slope function with the initial conditions.
End of explanation
"""


system.set(dt=0.1*s)
results, details = run_euler(system, slope_func, max_step=0.5)
details.message

results

def crossings(series, value):
    """Find the labels where the series passes through value.

    The labels in series must be increasing numerical values.

    series: Series
    value: number

    returns: sequence of labels
    """
    units = get_units(series.values[0])
    values = magnitudes(series - value)
    interp = InterpolatedUnivariateSpline(series.index, values)
    return interp.roots()

t_crossings = crossings(results.y, 0)

system.set(dt=0.1*s)
results, details = run_ralston(system, slope_func, max_step=0.5)
details.message

t_crossings = crossings(results.y, 0)

"""
Explanation: Now we're ready to call run_euler
End of explanation
"""


results

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


def plot_position(results):
    plot(results.y, label='y')
    decorate(xlabel='Time (s)',
             ylabel='Position (m)')

plot_position(results)
savefig('figs/chap09-fig01.pdf')

"""
Explanation: And here's position as a function of time:
End of explanation
"""


def crossings(series, value):
    """Find the labels where the series passes through value.

    The labels in series must be increasing numerical values.

    series: Series
    value: number

    returns: sequence of labels
    """
    units = get_units(series.values[0])
    values = magnitudes(series - value)
    interp = InterpolatedUnivariateSpline(series.index, values)
    return interp.roots()

t_crossings = crossings(results.y, 0)

"""
Explanation: Onto the sidewalk
To figure out when the penny hit the sidewalk, we can use crossings, which finds the times where a Series passes through a given value.
End of explanation
"""


t_sidewalk = t_crossings[0] * s

"""
Explanation: For this example there should be just one crossing, the time when the penny hits the sidewalk.
End of explanation
"""


sqrt(2 * init.y / g)

"""
Explanation: We can compare that to the exact result.  Without air resistance, we have
$v = -g t$
and
$y = 381 - g t^2 / 2$
Setting $y=0$ and solving for $t$ yields
$t = \sqrt{\frac{2 y_{init}}{g}}$
End of explanation
"""


def event_func(state, t, system):
    """Return the height of the penny above the sidewalk.
    """
    y, v = state
    return y

"""
Explanation: The estimate is accurate to about 10 decimal places.
Events
Instead of running the simulation until the penny goes through the sidewalk, it would be better to detect the point where the penny hits the sidewalk and stop.  run_ralston provides exactly the tool we need, event functions.
Here's an event function that returns the height of the penny above the sidewalk:
End of explanation
"""


results, details = run_ralston(system, slope_func, events=event_func)
details

"""
Explanation: And here's how we pass it to run_ralston.  The solver should run until the event function returns 0, and then terminate.
End of explanation
"""


results

"""
Explanation: The message from the solver indicates the solver stopped because the event we wanted to detect happened.
Here are the results:
End of explanation
"""


t_sidewalk = get_last_label(results) * s

"""
Explanation: With the events option, the solver returns the actual time steps it computed, which are not necessarily equally spaced. 
The last time step is when the event occurred:
End of explanation
"""


v_sidewalk = get_last_value(results.v)

"""
Explanation: The result is accurate to about 15 decimal places.
We can also check the velocity of the penny when it hits the sidewalk:
End of explanation
"""


km = UNITS.kilometer
h = UNITS.hour
v_sidewalk.to(km / h)

"""
Explanation: And convert to kilometers per hour.
End of explanation
"""


