Add adtype arg, rename to two_stage_objective, updates tests

Project.toml

@@ -6,16 +6,13 @@ version = "1.28.0"
 [deps]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PenaltyFunctions = "06bb1623-fdd5-5ca2-a01c-88eae3ea319e"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 
 [compat]
 Calculus = "0.5"
@@ -43,6 +40,16 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
 test = ["Test", "BlackBoxOptim", "DelayDiffEq", "LeastSquaresOptim", "NLopt", "Optim", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "SteadyStateDiffEq"]

src/build_loss_objective.jl

@@ -1,12 +1,12 @@
 export build_loss_objective
 
 function build_loss_objective(prob::SciMLBase.AbstractSciMLProblem, alg, loss,
+ adtype = SciMLBase.NoAD(),
  regularization = nothing, args...;
  priors = nothing,
  prob_generator = STANDARD_PROB_GENERATOR,
  kwargs...)
-
- cost_function = function (p, nothing)
+ cost_function = function (p, _p = nothing)
  tmp_prob = prob_generator(prob, p)
  if typeof(loss) <: Union{L2Loss, LogLikeLoss}
  sol = solve(tmp_prob, alg, args...; saveat = loss.t, save_everystep = false,
@@ -23,9 +23,8 @@ function build_loss_objective(prob::SciMLBase.AbstractSciMLProblem, alg, loss,
  if regularization !== nothing
  loss_val += regularization(p)
  end
-
  loss_val
  end
 
- OptimizationFunction(cost_function)
+ return OptimizationFunction(cost_function, adtype)
 end

src/multiple_shooting_objective.jl

@@ -22,12 +22,12 @@ struct Merged_Solution{T1, T2, T3}
 end
 
 function multiple_shooting_objective(prob::DiffEqBase.DEProblem, alg, loss,
+ adtype = SciMLBase.NoAD(),
  regularization = nothing; priors = nothing,
  discontinuity_weight = 1.0,
  prob_generator = STANDARD_MS_PROB_GENERATOR,
  kwargs...)
-
- cost_function = function (p, nothing)
+ cost_function = function (p, _p = nothing)
  t0, tf = prob.tspan
  P, N = length(prob.p), length(prob.u0)
  K = Int((length(p) - P) / N)
@@ -59,9 +59,10 @@ function multiple_shooting_objective(prob::DiffEqBase.DEProblem, alg, loss,
  if priors !== nothing
  loss_val += prior_loss(priors, p[(end - length(priors)):end])
  end
- if regularization !== nothing
+ if !isnothing(regularization)
  loss_val += regularization(p)
  end
+
  for k in 2:K
  if typeof(discontinuity_weight) <: Real
  loss_val += discontinuity_weight *
@@ -74,5 +75,5 @@ function multiple_shooting_objective(prob::DiffEqBase.DEProblem, alg, loss,
  loss_val
  end
 
- return OptimizationFunction(cost_function)
+ return OptimizationFunction(cost_function, adtype)
 end

src/two_stage_method.jl

@@ -1,12 +1,12 @@
-export TwoStageCost, two_stage_method
+export TwoStageCost, two_stage_objective
 
 struct TwoStageCost{F, D} <: Function
  cost_function::F
  estimated_solution::D
  estimated_derivative::D
 end
 
-(f::TwoStageCost)(p, _ = nothing) = f.cost_function(p, _)
+(f::TwoStageCost)(p, _p = nothing) = f.cost_function(p, _p)
 
 decide_kernel(kernel::CollocationKernel) = kernel
 function decide_kernel(kernel::Symbol)
@@ -70,8 +70,36 @@ function construct_estimated_solution_and_derivative!(data, kernel, tpoints)
  estimated_derivative, estimated_solution
 end
 
-function two_stage_objective(prob::DiffEqBase.DEProblem, tpoints, data;
- kernel = EpanechnikovKernel())
+function construct_iip_cost_function(f, du, preview_est_sol, preview_est_deriv, tpoints)
+ function (p, nothing)
+ _du = PreallocationTools.get_tmp(du, p)
+ vecdu = vec(_du)
+ cost = zero(first(p))
+ for i in 1:length(preview_est_sol)
+ est_sol = preview_est_sol[i]
+ f(_du, est_sol, p, tpoints[i])
+ vecdu .= vec(preview_est_deriv[i]) .- vec(_du)
+ cost += sum(abs2, vecdu)
+ end
+ cost
+ end
+end
+
+function construct_oop_cost_function(f, du, preview_est_sol, preview_est_deriv, tpoints)
+ function (p, nothing)
+ cost = zero(first(p))
+ for i in 1:length(preview_est_sol)
+ est_sol = preview_est_sol[i]
+ _du = f(est_sol, p, tpoints[i])
+ cost += sum(abs2, vec(preview_est_deriv[i]) .- vec(_du))
+ end
+ cost
+ end
+end
+
+function two_stage_objective(prob::DiffEqBase.DEProblem, tpoints, data,
+ adtype = SciMLBase.NoAD();
+ kernel = EpanechnikovKernel())
  f = prob.f
  kernel_function = decide_kernel(kernel)
  estimated_derivative, estimated_solution = construct_estimated_solution_and_derivative!(data,
@@ -84,16 +112,12 @@ function two_stage_objective(prob::DiffEqBase.DEProblem, tpoints, data;
  preview_est_deriv = [@view estimated_derivative[:, i]
  for i in 1:size(estimated_solution, 2)]
 
- function cost_function(p, nothing)
- cost = zero(first(p))
- for i in 1:length(preview_est_sol)
- est_sol = preview_est_sol[i]
- _du = f(est_sol, p, tpoints[i])
- cost += sum(abs2, vec(preview_est_deriv[i]) .- vec(_du))
- end
- cost
+ cost_function = if isinplace(prob)
+ construct_oop_cost_function(f, prob.u0, preview_est_sol, preview_est_deriv, tpoints)
+ else
+ construct_iip_cost_function(f, prob.u0, preview_est_sol, preview_est_deriv, tpoints)
  end
 
  return OptimizationFunction(TwoStageCost(cost_function, estimated_solution,
- estimated_derivative))
+  estimated_derivative), adtype)
 end

test/dae_tests.jl

@@ -1,31 +1,36 @@
 using DelayDiffEq, OrdinaryDiffEq, RecursiveArrayTools, Test,
  Sundials
 
-function resrob(r, yp, y, p, tres)
- r[1] = -p[1] * y[1] + 1.0e4 * y[2] * y[3]
- r[2] = -r[1] - 3.0e7 * y[2] * y[2] - yp[2]
- r[1] -= yp[1]
- r[3] = y[1] + y[2] + y[3] - 1.0
+function f(out, du, u, p, t)
+ out[1] = -p[1] * u[1] + 1e4 * u[2] * u[3] - du[1]
+ out[2] = +p[1] * u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
+ out[3] = u[1] + u[2] + u[3] - 1.0
 end
-u0 = [1.0, 0, 0]
-du0 = [-0.04, 0.04, 0.0]
-p = [0.04]
-prob = DAEProblem(pf, u0, du0, (0.0, 100000.0), p, differential_vars = [3])
-sol = solve(prob, IDA())
+u₀ = [1.0, 0, 0]
+du₀ = [-0.04, 0.04, 0.0]
+tspan = (0.0, 100000.0)
+differential_vars = [true, true, false]
+prob = DAEProblem(f, du₀, u₀, tspan, [0.04], differential_vars = differential_vars)
+sol = solve(prob, DFBDF())
 
 t = collect(range(0, stop = 10, length = 30))
 randomized = VectorOfArray([(sol(t[i]) + 0.003randn(3)) for i in 1:length(t)])
 data = convert(Array, randomized)
 
-using DiffEqParamEstim, NLopt
-cost_function = build_loss_objective(prob, IDA(), L2Loss(t, data), maxiter = 10000,
- abstol = 1e-8, reltol = 1e-8, verbose = false)
+using DiffEqParamEstim, OptimizationNLopt, OptimizationOptimJL, ForwardDiff, FiniteDiff,
+ Zygote, Optimization, SciMLSensitivity
+cost_function = build_loss_objective(prob, DFBDF(), L2Loss(t, data),
+ Optimization.AutoZygote(), abstol = 1e-8,
+ reltol = 1e-8, verbose = false)
+optprob = OptimizationProblem(cost_function, [0.01]; lb = [0.0], ub = [1.0])
+res = solve(optprob, OptimizationOptimJL.BFGS())
 
+cost_function = build_loss_objective(prob, DFBDF(), L2Loss(t, data),
+ Optimization.AutoForwardDiff(), abstol = 1e-8,
+ reltol = 1e-8, verbose = false)
+optprob = OptimizationProblem(cost_function, [0.01]; lb = [0.0], ub = [1.0])
+res = solve(optprob, OptimizationOptimJL.BFGS())
 opt = Opt(:GN_ESCH, 1)
-min_objective!(opt, cost_function.cost_function2)
-lower_bounds!(opt, [0.0])
-upper_bounds!(opt, [1.0])
-maxeval!(opt, 10000)
 
-(minf, minx, ret) = NLopt.optimize(opt, [0.2])
-@test minx[1]≈0.04 atol=5e-3
+res = solve(optprob, opt)
+# @test minx[1]≈0.04 atol=5e-3

test/dde_tests.jl

@@ -17,7 +17,7 @@ t = collect(range(0, stop = 10, length = 30))
 randomized = VectorOfArray([(sol(t[i]) + 0.03randn(2)) for i in 1:length(t)])
 data = convert(Array, randomized)
 
-using DiffEqParamEstim, NLopt
+using DiffEqParamEstim, OptimizationNLopt
 
 function f_lotka2(du, u, h, p, t)
  du[1] = 0.5 * u[1] - p[1] * u[1] * u[2]
@@ -29,14 +29,12 @@ p = [0.5]
 
 prob_opt = DDEProblem(f_lotka2, u0, h, tspan, p, constant_lags = [0.5])
 cost_function = build_loss_objective(prob_opt, MethodOfSteps(Tsit5()),
- L2Loss(t, data), maxiter = 10000, abstol = 1e-8,
- reltol = 1e-8, verbose = false)
+ L2Loss(t, data), Optimization.AutoZygote(),
+ abstol = 1e-8,
+ reltol = 1e-8)
 
+optprob = OptimizationProblem(cost_function, [1.0], lb = [0.0], ub = [1.0])
 opt = Opt(:GN_ESCH, 1)
-min_objective!(opt, cost_function.cost_function2)
-lower_bounds!(opt, [0.0])
-upper_bounds!(opt, [1.0])
-maxeval!(opt, 10000)
+res = solve(optprob, BFGS())
 
-(minf, minx, ret) = NLopt.optimize(opt, [0.2])
 @test minx[1]≈0.5 atol=5e-3

test/likelihood.jl

@@ -1,4 +1,5 @@
-using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Test, RecursiveArrayTools
+using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Test, RecursiveArrayTools,
+ OptimizationBBO
 # ,BlackBoxOptim,
 
 pf_func = function (du, u, p, t)
@@ -18,32 +19,37 @@ end
 aggregate_data = convert(Array, VectorOfArray([generate_data(sol, t) for i in 1:100]))
 
 distributions = [fit_mle(Normal, aggregate_data[i, j, :]) for i in 1:2, j in 1:200]
-obj = build_loss_objective(prob1, Tsit5(), LogLikeLoss(t, distributions),
- maxiters = 10000, verbose = false)
-bound1 = Tuple{Float64, Float64}[(0.5, 5), (0.5, 5)]
-result = bboptimize(obj; SearchRange = bound1, MaxSteps = 11e3)
-@test result.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1
+obj = build_loss_objective(prob1, Tsit5(), LogLikeLoss(t, distributions), maxiters = 10000,
+ verbose = false)
+
+optprob = OptimizationProblem(obj, [2.0, 2.0], lb = [0.5, 0.5], ub = [5.0, 5.0])
+result = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 11e3)
+@test result.original.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1
 
 data_distributions = [fit_mle(Normal, aggregate_data[i, j, :]) for i in 1:2, j in 1:200]
 diff_distributions = [fit_mle(Normal,
  aggregate_data[i, j, :] - aggregate_data[i, j - 1, :])
  for j in 2:200, i in 1:2]
 obj = build_loss_objective(prob1, Tsit5(),
  LogLikeLoss(t, data_distributions, diff_distributions),
- maxiters = 10000, verbose = false)
-bound1 = Tuple{Float64, Float64}[(0.5, 5), (0.5, 5)]
-result = bboptimize(obj; SearchRange = bound1, MaxSteps = 11e3)
-@test result.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1
+ Optimization.AutoForwardDiff(), maxiters = 10000,
+ verbose = false)
+optprob = OptimizationProblem(obj, [2.0, 2.0], lb = [0.5, 0.5], ub = [5.0, 5.0])
+result = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 11e3)
+@test result.original.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1
 
 data_distributions = [fit_mle(Normal, aggregate_data[i, j, :]) for i in 1:2, j in 1:200]
 diff_distributions = [fit_mle(Normal,
  aggregate_data[i, j, :] - aggregate_data[i, j - 1, :])
  for j in 2:200, i in 1:2]
 obj = build_loss_objective(prob1, Tsit5(),
  LogLikeLoss(t, data_distributions, diff_distributions, 0.3),
- maxiters = 10000, verbose = false)
-bound1 = Tuple{Float64, Float64}[(0.5, 5), (0.5, 5)]
-result = bboptimize(obj; SearchRange = bound1, MaxSteps = 11e3)
+ Optimization.AutoForwardDiff(), maxiters = 10000,
+ verbose = false)
+optprob = OptimizationProblem(obj, [2.0, 2.0], lb = [0.5, 0.5], ub = [5.0, 5.0])
+result = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 11e3)
+using OptimizationBBO.BlackBoxOptim
+bboptimize((x) -> obj(x, nothing), SearchRange = bound1)
 @test result.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1
 
 distributions = [fit_mle(MvNormal, aggregate_data[:, j, :]) for j in 1:200]
@@ -53,7 +59,8 @@ diff_distributions = [fit_mle(MvNormal,
 priors = [Truncated(Normal(1.5, 0.1), 0, 2), Truncated(Normal(1.0, 0.1), 0, 1.5)]
 obj = build_loss_objective(prob1, Tsit5(),
  LogLikeLoss(t, distributions, diff_distributions),
- maxiters = 10000, verbose = false, priors = priors)
-bound1 = Tuple{Float64, Float64}[(0.5, 5), (0.5, 5)]
-result = bboptimize(obj; SearchRange = bound1, MaxSteps = 11e3)
+ Optimization.AutoForwardDiff(), maxiters = 10000,
+ verbose = false, priors = priors)
+optprob = OptimizationProblem(obj, [2.0, 2.0], lb = [0.5, 0.5], ub = [5.0, 5.0])
+result = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 11e3)
 @test result.archive_output.best_candidate≈[1.5, 1.0] atol=1e-1