DDE parameter estimation

blackboxoptim.jl

using DelayDiffEq, DiffEqParamEstim, BlackBoxOptim, DataFrames, LsqFit
using Plots
gr()

include("importData.jl")
include("plot.jl")

# import data from the path, in which:
# pop: population data
# g1, g2: g1 and g2 data
# initial: initial number of cells in g1 and in g2 at time 0
pop, g2, g1, g2_0, g1_0 = get_data(joinpath("..", "data", "lap.csv"),
                                   joinpath("..", "data", "lap_pop.csv"))

# time points
times = range(0.0; stop = 95.5, length = 192)

# model
function DDEmodel(du, u, h, p, t)
    du[1] = -p[1]*(h(p, t-p[3])[1]) + 2*p[2]*(h(p, t-p[4])[2]) - p[5]*u[1]
    du[2] = p[1]*(h(p, t-p[3])[1]) - p[2]*(h(p, t-p[4])[2]) - p[6]*u[2]
end

# estimate history function
exp_model(t, p) = @. p[1] * exp(t * p[2]) # exponential model
fit1 = curve_fit(exp_model, times, g1[:, 1], [1.0, 0.5])
fit2 = curve_fit(exp_model, times, g2[:, 1], [1.0, 0.5])
h(p, t) = [exp_model(t, coef(fit1)); exp_model(t, coef(fit2))]

# initial guess
initial_guess = [0.02798, 0.025502, 21.3481, 10.2881, 0.0001, 0.0001]

# problem
prob = DDEProblem(DDEmodel, [g1_0[3], g2_0[3]], h, extrema(times), initial_guess;
                  constant_lags = [initial_guess[3], initial_guess[4]])

# DDE algorithm:
# switch automatically between non-stiff solver Tsit5 and stiff solver Rosenbrock23
alg = MethodOfSteps(AutoTsit5(Rosenbrock23()))

# plot initial guess
plot(solve(prob, alg), label = ["u1(t) (guess)", "u2(t) (guess)"], legend = :topleft,
     show = true)

# plot data
data = vcat(g1[:, 3]', g2[:, 3]')
scatter!(times, data', show = true)

# define problem generator (optimization in log space)
function prob_generator(prob, p)
    exp_p = exp.(p)
    remake(prob; p = exp_p, constant_lags = [exp_p[3], exp_p[4]])
end

# define objective function (L2 loss)
obj = build_loss_objective(prob, alg, L2Loss(times, data);
                           prob_generator = prob_generator)

# perform global optimization
results_dde = bboptimize(obj;
                         SearchRange = [(-6.0, 6.0), (-6.0, 6.0), (0.0, 6.0),
                                        (0.0, 6.0), (-10.0, 0.0), (-10.0, 0.0)],
                         NumDimensions = 6)

# plot solution with best candidate
min_p = exp.(best_candidate(results_dde))
plot!(solve(remake(prob; p = min_p, constant_lags = [min_p[3], min_p[4]]), alg), show = true)
println("minimum = " , best_fitness(results_dde), " with argmin = ", min_p)

nlopt.jl

using DelayDiffEq, DiffEqParamEstim, NLopt, DataFrames, LsqFit
using Plots
gr()

include("importData.jl")
include("plot.jl")

# import data from the path, in which:
# pop: population data
# g1, g2: g1 and g2 data
# initial: initial number of cells in g1 and in g2 at time 0
pop, g2, g1, g2_0, g1_0 = get_data(joinpath("..", "data", "lap.csv"),
                                   joinpath("..", "data", "lap_pop.csv"))

# time points
times = range(0.0; stop = 95.5, length = 192)

# model
function DDEmodel(du, u, h, p, t)
    du[1] = -p[1]*(h(p, t-p[3])[1]) + 2*p[2]*(h(p, t-p[4])[2]) - p[5]*u[1]
    du[2] = p[1]*(h(p, t-p[3])[1]) - p[2]*(h(p, t-p[4])[2]) - p[6]*u[2]
end

# estimate history function
exp_model(t, p) = @. p[1] * exp(t * p[2]) # exponential model
fit1 = curve_fit(exp_model, times, g1[:, 1], [1.0, 0.5])
fit2 = curve_fit(exp_model, times, g2[:, 1], [1.0, 0.5])
h(p, t) = [exp_model(t, coef(fit1)); exp_model(t, coef(fit2))]

# initial guess
initial_guess = [0.02798, 0.025502, 21.3481, 10.2881, 0.0001, 0.0001]

# problem
prob = DDEProblem(DDEmodel, [g1_0[3], g2_0[3]], h, extrema(times), initial_guess;
                  constant_lags = [initial_guess[3], initial_guess[4]])

# DDE algorithm:
# switch automatically between non-stiff solver Tsit5 and stiff solver Rosenbrock23
alg = MethodOfSteps(AutoTsit5(Rosenbrock23()))

# plot initial guess
plot(solve(prob, alg), label = ["u1(t) (guess)", "u2(t) (guess)"], legend = :topleft,
     show = true)

# plot data
data = vcat(g1[:, 3]', g2[:, 3]')
scatter!(times, data', show = true)

# define problem generator (optimization in log space)
function prob_generator(prob, p)
    exp_p = exp.(p)
    remake(prob; p = exp_p, constant_lags = [exp_p[3], exp_p[4]])
end
# function prob_generator(prob, p)
#     exp_p = exp.(p)
#     T = eltype(exp_p)
#     remake(prob; p = exp_p, u0 = T.(prob.u0), tspan = T.(prob.tspan),
#            constant_lags = [exp_p[3], exp_p[4]])
# end

# define objective function (L2 loss)
obj = build_loss_objective(prob, alg, L2Loss(times, data);
                           prob_generator = prob_generator,
                           verbose_opt = true)
# obj = build_loss_objective(prob, alg, L2Loss(times, data);
#                            prob_generator = prob_generator,
#                            verbose_opt = true, mpg_autodiff = true)

# perform global optimization
# opt = Opt(:GD_STOGO, 6)
opt = Opt(:GN_ESCH, 6)
min_objective!(opt, obj)
lower_bounds!(opt, [-6.0, -6.0, 0.0, 0.0, -10.0, -10.0])
upper_bounds!(opt, [6.0, 6.0, 6.0, 6.0, 0.0, 0.0])
maxeval!(opt, 10_000)
minf, minx, ret = NLopt.optimize(opt, log.(initial_guess))

# plot solution with best candidate
min_p = exp.(minx)
plot!(solve(remake(prob; p = min_p, constant_lags = [min_p[3], min_p[4]]), alg), show = true)
println("minimum = " , minf, " with argmin = ", min_p)

optim.jl

using DelayDiffEq, DiffEqParamEstim, Optim, DataFrames, LsqFit
using Plots
gr()

include("importData.jl")
include("plot.jl")

# import data from the path, in which:
# pop: population data
# g1, g2: g1 and g2 data
# initial: initial number of cells in g1 and in g2 at time 0
pop, g2, g1, g2_0, g1_0 = get_data(joinpath("..", "data", "lap.csv"),
                                   joinpath("..", "data", "lap_pop.csv"))

# time points
times = range(0.0; stop = 95.5, length = 192)

# model
function DDEmodel(du, u, h, p, t)
    du[1] = -p[1]*(h(p, t-p[3])[1]) + 2*p[2]*(h(p, t-p[4])[2]) - p[5]*u[1]
    du[2] = p[1]*(h(p, t-p[3])[1]) - p[2]*(h(p, t-p[4])[2]) - p[6]*u[2]
end

# estimate history function
exp_model(t, p) = @. p[1] * exp(t * p[2]) # exponential model
fit1 = curve_fit(exp_model, times, g1[:, 1], [1.0, 0.5])
fit2 = curve_fit(exp_model, times, g2[:, 1], [1.0, 0.5])
h(p, t) = [exp_model(t, coef(fit1)); exp_model(t, coef(fit2))]

# initial guess
initial_guess = [0.02798, 0.025502, 21.3481, 10.2881, 0.0001, 0.0001]

# problem
prob = DDEProblem(DDEmodel, [g1_0[3], g2_0[3]], h, extrema(times), initial_guess;
                  constant_lags = [initial_guess[3], initial_guess[4]])

# DDE algorithm:
# switch automatically between non-stiff solver Tsit5 and stiff solver Rosenbrock23
alg = MethodOfSteps(AutoTsit5(Rosenbrock23()))

# plot initial guess
plot(solve(prob, alg), label = ["u1(t) (guess)", "u2(t) (guess)"], legend = :topleft,
     show = true)

# plot data
data = vcat(g1[:, 3]', g2[:, 3]')
scatter!(times, data', show = true)

# define problem generator (optimization in log space)
function prob_generator(prob, p)
    exp_p = exp.(p)
    remake(prob; p = exp_p, constant_lags = [exp_p[3], exp_p[4]])
end
# function prob_generator(prob, p)
#     exp_p = exp.(p)
#     T = eltype(exp_p)
#     remake(prob; p = exp_p, u0 = T.(prob.u0), tspan = T.(prob.tspan),
#            constant_lags = [exp_p[3], exp_p[4]])
# end

# define objective function (L2 loss)
obj = build_loss_objective(prob, alg, L2Loss(times, data);
                           prob_generator = prob_generator,
                           verbose_opt = true)
# obj = build_loss_objective(prob, alg, L2Loss(times, data);
#                            prob_generator = prob_generator,
#                            verbose_opt = true, mpg_autodiff = true)

# perform global optimization
results_dde = optimize(obj, [-6.0, -6.0, 0.0, 0.0, -10.0, -10.0],
                       [6.0, 6.0, 6.0, 6.0, 0.0, 0.0],
                       log.(initial_guess),
                       Fminbox(Optim.LBFGS()))
# results_dde = optimize(obj, obj, [-6.0, -6.0, 0.0, 0.0, -10.0, -10.0],
#                        [6.0, 6.0, 6.0, 6.0, 0.0, 0.0],
#                        log.(initial_guess),
#                        Fminbox(Optim.LBFGS()))

# plot solution with best candidate
min_p = exp.(Optim.minimizer(results_dde))
plot!(solve(remake(prob; p = min_p, constant_lags = [min_p[3], min_p[4]]), alg), show = true)
println("minimum = " , Optim.minimum(results_dde), " with argmin = ", min_p)