Assimulo and sensitivity

gistfile1.py

from assimulo.solvers import CVode
from assimulo.problem import Explicit_Problem

from scipy.integrate import odeint
from scipy.optimize import fmin_l_bfgs_b

import numpy as np
import nlopt
import matplotlib.pyplot as plt

outside_variable = 'JIMMY'

def mm_scipy(y, t, *args):
    Vmax = args[0][0]
    km = args[0][1]
    St = args[0][2]
    P = y[0]
    S = St - P
    
    dP = Vmax * (S / (S+km))
    return dP

def mm_assimulo(t, y, p):
    Vmax = p[0]
    km = p[1]
    St = p[2]
    
    P = y
    S = St - P
    
    dP = Vmax * (S / (S+km))
    dP = np.array(dP)
    return dP

# Parameters MM
Vmax = 1
km = 3
St = 10
p0 = np.array([Vmax, km, St])
n_params = len(p0)

# Initial Conditions MM
y0 = [0.1]
yS0 = np.array([[0.0], [0.0], [0.0]])


exp_mod = Explicit_Problem(mm_assimulo, y0 = y0, p0 = p0)
exp_mod.yS0 = yS0
exp_sim = CVode(exp_mod)

exp_sim.iter = 'Newton'
exp_sim.discr = 'BDF'
exp_sim.rtol = 1e-7
exp_sim.atol = 1e-6
exp_sim.pbar = [1,1,1] #pbar is used to estimate the tolerances for the parameters
exp_sim.report_continuously = True #Need to be able to store the result using the interpolate methods
exp_sim.sensmethod = 'SIMULTANEOUS' #Defines the sensitvity method used
exp_sim.suppress_sens = False            #Dont suppress the sensitivity variables in the error test.

t, assi_y = exp_sim.simulate(100) #Simulate 100 seconds



sci_y = odeint(mm_scipy, y0, t, args = (p0,))

plt.plot(t, assi_y, 'bo')
plt.plot(t, sci_y, 'r--')
# Ok - exactly same results.

n_steps = len(t)
    
# Generate random data:
noisy_y = (assi_y + np.random.randn(n_steps, 1))[::15]
noisy_t = t[::15]



# Create a function that takes as an input a vector of parameters,
# and returns RSS and gradient, as a function of that vector:

exp_mod = Explicit_Problem(mm_assimulo, y0 = y0,
                           p0 = np.array([0.0, 0, 0]))
exp_mod.yS0 = yS0
exp_sim = CVode(exp_mod)

exp_sim.iter = 'Newton'
exp_sim.discr = 'BDF'
exp_sim.rtol = 1e-7
exp_sim.atol = 1e-6
exp_sim.pbar = [1,1,1] #pbar is used to estimate the tolerances for the parameters
exp_sim.report_continuously = True #Need to be able to store the result using the interpolate methods
exp_sim.sensmethod = 'SIMULTANEOUS' #Defines the sensitvity method used
exp_sim.suppress_sens = False
exp_sim.verbosity = 50

def calc_rss_and_grad(p):
    
    exp_sim.reset()
    #exp_sim.p = p
    exp_sim.p0 = p

    t, assi_y = exp_sim.simulate(100) #Simulate 100 seconds
    n_steps = len(t)
    # Sensitivities:
    sens = np.empty((n_steps, n_params))
    for i in range(n_params):
        sens[:, i] = np.array(exp_sim.p_sol[i]).flatten()
        
    t_idx = np.searchsorted(t, noisy_t)
    t_idx[t_idx>len(t)-1] = len(t)-1
    
    res = assi_y[t_idx] - noisy_y
    
    
    # Calculate gradient at every timepoint:
    # It's (f(p)_i - y_i) * df(p)/dp_j - last term is sensitivity
    gradient = res * sens[t_idx, :]

    gradient[np.isnan(gradient)] = 0

    # Sum across timepoints
    gradient = np.sum(gradient, axis = 0)
    
    rss = np.sum(res**2)
    
    return rss, gradient
    
def nlopt_calc_rss_and_grad(p, grad):
    exp_sim.reset()
    #exp_sim.p = p
    exp_sim.p0 = p
    
    t, assi_y = exp_sim.simulate(100) #Simulate 100 seconds
    n_steps = len(t)
    # Sensitivities:
    sens = np.empty((n_steps, n_params))
    for i in range(n_params):
        sens[:, i] = np.array(exp_sim.p_sol[i]).flatten()
        
    t_idx = np.searchsorted(t, noisy_t)
    t_idx[t_idx>len(t)-1] = len(t)-1
    
    res = assi_y[t_idx] - noisy_y
    
    
    # Calculate gradient at every timepoint:
    # It's (f(p)_i - y_i) * df(p)/dp_j - last term is sensitivity
    gradient = res * sens[t_idx, :]

    gradient[np.isnan(gradient)] = 0

    # Sum across timepoints
    gradient = np.sum(gradient, axis = 0)
    
    rss = np.sum(res**2)
    
    if grad.size > 0:
        grad[:] = gradient
    
    return rss
    
x_bfgs, f, d = fmin_l_bfgs_b(calc_rss_and_grad,
                        x0 = np.array([0.0, 0, 0]), fprime = None)
                        
opt = nlopt.opt(nlopt.LD_MMA, 3)
opt.set_lower_bounds([0.0, 0.0, 0.0])
opt.set_upper_bounds([10.0, 5.0, 15.0])
opt.set_min_objective(nlopt_calc_rss_and_grad)
opt.set_maxtime(10)
x_nlopt = opt.optimize(np.array([0.0, 0, 0]))
minf = opt.last_optimum_value()