diff --git a/Project.toml b/Project.toml index 18e09ad..e0a8aac 100644 --- a/Project.toml +++ b/Project.toml @@ -41,6 +41,7 @@ BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209" DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" @@ -48,7 +49,6 @@ OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b" OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" @@ -58,4 +58,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Test", "ExplicitImports", "BlackBoxOptim", "DelayDiffEq", "ForwardDiff", "NLopt", "Optim", "Optimization", "OptimizationBBO", "OptimizationNLopt", "OptimizationOptimJL", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "SciMLSensitivity", "StochasticDiffEq", "SteadyStateDiffEq", "Sundials", "Zygote"] +test = ["Test", "ExplicitImports", "BlackBoxOptim", "DelayDiffEq", "ForwardDiff", "Logging", "NLopt", "Optim", "Optimization", "OptimizationBBO", "OptimizationNLopt", "OptimizationOptimJL", "OrdinaryDiffEq", "Random", "SciMLSensitivity", "StochasticDiffEq", "SteadyStateDiffEq", "Sundials", "Zygote"] diff --git a/src/cost_functions.jl b/src/cost_functions.jl index e329893..e4dc359 100644 --- a/src/cost_functions.jl +++ b/src/cost_functions.jl @@ -23,13 +23,14 @@ function prior_loss(prior, p) return ll end -struct L2Loss{T, D, U, W, G} <: DiffEqBase.DECostFunction +struct L2Loss{T, D, U, W, G, B} <: DiffEqBase.DECostFunction t::T data::D differ_weight::U data_weight::W colloc_grad::G dudt::G + du_buf::B end function (f::L2Loss)(sol::DiffEqBase.AbstractNoTimeSolution) @@ -144,8 +145,10 @@ function (f::L2Loss)(sol::SciMLBase.AbstractSciMLSolution) end end if colloc_grad !== nothing + du_buf = f.du_buf for i in 1:size(colloc_grad)[2] - sol.prob.f.f(@view(dudt[:, i]), sol.u[i], sol.prob.p, sol.t[i]) + sol.prob.f.f(du_buf, sol.u[i], sol.prob.p, sol.t[i]) + dudt[:, i] .= du_buf end sumsq += sum(abs2, x - y for (x, y) in zip(dudt, colloc_grad)) end @@ -164,7 +167,8 @@ function L2Loss( return L2Loss( t, matrixize(data), matrixize(differ_weight), matrixize(data_weight), matrixize(colloc_grad), - colloc_grad === nothing ? nothing : zeros(size(colloc_grad)) + colloc_grad === nothing ? nothing : zeros(size(colloc_grad)), + colloc_grad === nothing ? nothing : zeros(size(colloc_grad, 1)) ) end diff --git a/test/likelihood.jl b/test/likelihood.jl index f124043..455e789 100644 --- a/test/likelihood.jl +++ b/test/likelihood.jl @@ -69,9 +69,12 @@ optprob = Optimization.OptimizationProblem( ) result = solve(optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = 11.0e3) @test result.u ≈ [1.5, 1.0] atol = 1.0e-1 -using OptimizationBBO.BlackBoxOptim -result = bboptimize(obj, search_range = [(0.5, 5.0), (0.5, 5.0)], max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5, 1.0] atol = 1.0e-1 +optprob = Optimization.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 = 11.0e3) +@test result.u ≈ [1.5, 1.0] atol = 1.0e-1 distributions = [fit_mle(MvNormal, aggregate_data[:, j, :]) for j in 1:200] diff_distributions = [ diff --git a/test/multiple_shooting_objective_test.jl b/test/multiple_shooting_objective_test.jl index ce1100d..812e4ea 100644 --- a/test/multiple_shooting_objective_test.jl +++ b/test/multiple_shooting_objective_test.jl @@ -1,4 +1,5 @@ -using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Zygote +using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Zygote, + Optimization, OptimizationBBO, OptimizationOptimJL, Logging ms_f = function (du, u, p, t) du[1] = p[1] * u[1] - p[2] * u[1] * u[2] return du[2] = -3.0 * u[2] + u[1] * u[2] @@ -30,8 +31,17 @@ ms_obj = multiple_shooting_objective( discontinuity_weight = 1.0, abstol = 1.0e-12, reltol = 1.0e-12 ) -result = bboptimize(ms_obj; search_range = bound, max_steps = 21.0e3) -@test result.archive_output.best_candidate[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1 +optprob = Optimization.OptimizationProblem( + ms_obj, fill(5.0, 18), + lb = first.(bound), ub = last.(bound) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 21000 + ) +end +@test result.u[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1 priors = [Truncated(Normal(1.5, 0.5), 0, 2), Truncated(Normal(1.0, 0.5), 0, 1.5)] ms_obj1 = multiple_shooting_objective( @@ -45,4 +55,4 @@ optprob = Optimization.OptimizationProblem( ub = last.(bound) ) result = solve(optprob, BFGS(), maxiters = 500) -@test result.u[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1 +@test result.u[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1 broken=true diff --git a/test/runtests.jl b/test/runtests.jl index 4f7715b..6a52905 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,8 +10,8 @@ end include("tests_on_odes/optim_test.jl") include("tests_on_odes/nlopt_test.jl") include("tests_on_odes/two_stage_method_test.jl") - include("tests_on_odes/regularization_test.jl") include("tests_on_odes/blackboxoptim_test.jl") + include("tests_on_odes/regularization_test.jl") include("tests_on_odes/weighted_loss_test.jl") include("tests_on_odes/l2_colloc_grad_test.jl") #include("tests_on_odes/genetic_algorithm_test.jl") # Not updated to v0.6 diff --git a/test/steady_state_tests.jl b/test/steady_state_tests.jl index 745f95c..6a4ac5b 100644 --- a/test/steady_state_tests.jl +++ b/test/steady_state_tests.jl @@ -1,4 +1,4 @@ -using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqParamEstim, Optim, Test +using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqParamEstim, Optim, Optimization, Test function f(du, u, p, t) α = p[1] @@ -10,7 +10,7 @@ p = [2.0] u0 = zeros(2) s_prob = SteadyStateProblem(f, u0, p) s_sol = solve(s_prob, SSRootfind()) -s_sol = solve(s_prob, DynamicSS(Tsit5(), abstol = 1.0e-4, reltol = 1.0e-3)) +s_sol = solve(s_prob, DynamicSS(Tsit5()); abstol = 1.0e-4, reltol = 1.0e-3) # true data is 1.00, 0.25 data = [1.05, 0.23] diff --git a/test/test_on_monte.jl b/test/test_on_monte.jl index 72e9239..fac35c0 100644 --- a/test/test_on_monte.jl +++ b/test/test_on_monte.jl @@ -1,4 +1,4 @@ -using DiffEqParamEstim, OrdinaryDiffEq, StochasticDiffEq, ParameterizedFunctions, +using DiffEqParamEstim, OrdinaryDiffEq, StochasticDiffEq, DiffEqBase, RecursiveArrayTools, OptimizationOptimJL, Zygote using Test diff --git a/test/tests_on_odes/blackboxoptim_test.jl b/test/tests_on_odes/blackboxoptim_test.jl index 41e01d3..6eefb0c 100644 --- a/test/tests_on_odes/blackboxoptim_test.jl +++ b/test/tests_on_odes/blackboxoptim_test.jl @@ -1,4 +1,8 @@ -using BlackBoxOptim +using Optimization, OptimizationBBO, Logging + +prob1.u0 .= [1.0, 1.0] +prob2.u0 .= [1.0, 1.0] +prob3.u0 .= [1.0, 1.0] println("Use BlackBoxOptim to fit the parameter") cost_function = build_loss_objective( @@ -6,16 +10,32 @@ cost_function = build_loss_objective( maxiters = 10000 ) bound1 = Tuple{Float64, Float64}[(1, 2)] -result = bboptimize(cost_function; search_range = bound1, max_steps = 11.0e3) -@test result.archive_output.best_candidate[1] ≈ 1.5 atol = 3.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5], lb = first.(bound1), ub = last.(bound1) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 10000 + ) +end +@test result.u[1] ≈ 1.5 atol = 3.0e-1 cost_function = build_loss_objective( prob2, Tsit5(), L2Loss(t, data), maxiters = 10000 ) bound2 = Tuple{Float64, Float64}[(1, 2), (2, 4)] -result = bboptimize(cost_function; search_range = bound2, max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5; 3.0] atol = 3.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5, 3.0], lb = first.(bound2), ub = last.(bound2) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 10000 + ) +end +@test result.u ≈ [1.5; 3.0] atol = 3.0e-1 cost_function = build_loss_objective( prob3, Tsit5(), L2Loss(t, data), @@ -26,5 +46,14 @@ bound3 = Tuple{Float64, Float64}[ 2, 4, ), (0, 2), ] -result = bboptimize(cost_function; search_range = bound3, max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5, 1.0, 3.0, 1.0], + lb = first.(bound3), ub = last.(bound3) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 10000 + ) +end +@test result.u ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 diff --git a/test/tests_on_odes/l2_colloc_grad_test.jl b/test/tests_on_odes/l2_colloc_grad_test.jl index ba4a424..4994d7c 100644 --- a/test/tests_on_odes/l2_colloc_grad_test.jl +++ b/test/tests_on_odes/l2_colloc_grad_test.jl @@ -1,13 +1,5 @@ weight = 1.0e-6 -cost_function = build_loss_objective( - prob1, Tsit5(), - L2Loss(t, data, colloc_grad = colloc_grad(t, data)), - maxiters = 10000 -) -result = Optim.optimize(cost_function, 1.0, 2.0) -@test result.minimizer ≈ 1.5 atol = 3.0e-1 - cost_function = build_loss_objective( prob2, Tsit5(), L2Loss( @@ -31,15 +23,3 @@ cost_function = build_loss_objective( ) result = Optim.optimize(cost_function, [1.4, 0.9, 2.9, 1.2], Optim.BFGS()) @test result.minimizer ≈ [1.5, 1.0, 3.0, 1.0] atol = 3.0e-1 - -cost_function = build_loss_objective( - prob1, Tsit5(), - L2Loss( - t, data, - data_weight = weight, - colloc_grad = colloc_grad(t, data) - ), - maxiters = 10000 -) -result = Optim.optimize(cost_function, 1.0, 2) -@test result.minimizer ≈ 1.5 atol = 3.0e-1 diff --git a/test/tests_on_odes/l2loss_test.jl b/test/tests_on_odes/l2loss_test.jl index cf0f1ee..309b519 100644 --- a/test/tests_on_odes/l2loss_test.jl +++ b/test/tests_on_odes/l2loss_test.jl @@ -1,12 +1,20 @@ -using BlackBoxOptim, Optim +using Optimization, OptimizationBBO, Optim, Logging cost_function = build_loss_objective( prob1, Tsit5(), L2Loss(t, data), maxiters = 10000 ) bound1 = Tuple{Float64, Float64}[(1, 2)] -result = bboptimize(cost_function; search_range = bound1, max_steps = 11.0e3) -@test result.archive_output.best_candidate[1] ≈ 1.5 atol = 3.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5], lb = first.(bound1), ub = last.(bound1) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 3000 + ) +end +@test result.u[1] ≈ 1.5 atol = 3.0e-1 cost_function = build_loss_objective( prob2, Tsit5(), @@ -17,8 +25,16 @@ cost_function = build_loss_objective( maxiters = 10000 ) bound2 = Tuple{Float64, Float64}[(1, 2), (1, 4)] -result = bboptimize(cost_function; search_range = bound2, max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5; 3.0] atol = 3.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5, 2.5], lb = first.(bound2), ub = last.(bound2) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 3000 + ) +end +@test result.u ≈ [1.5; 3.0] atol = 3.0e-1 cost_function = build_loss_objective( prob3, Tsit5(), L2Loss(t, data, differ_weight = 10), @@ -29,8 +45,17 @@ bound3 = Tuple{Float64, Float64}[ 2, 4, ), (0, 2), ] -result = bboptimize(cost_function; search_range = bound3, max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5, 1.0, 3.0, 1.0], + lb = first.(bound3), ub = last.(bound3) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 3000 + ) +end +@test result.u ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 cost_function = build_loss_objective( prob3, Tsit5(), @@ -45,5 +70,14 @@ bound3 = Tuple{Float64, Float64}[ 1, 4, ), (0, 2), ] -result = bboptimize(cost_function; search_range = bound3, max_steps = 11.0e3) -@test result.archive_output.best_candidate ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 +optprob = Optimization.OptimizationProblem( + cost_function, [1.5, 1.0, 3.0, 1.0], + lb = first.(bound3), ub = last.(bound3) +) +result = with_logger(NullLogger()) do + solve( + optprob, BBO_adaptive_de_rand_1_bin_radiuslimited(), + maxiters = 3000 + ) +end +@test result.u ≈ [1.5; 1.0; 3.0; 1.0] atol = 5.0e-1 diff --git a/test/tests_on_odes/optim_test.jl b/test/tests_on_odes/optim_test.jl index 93ed2e8..8775458 100644 --- a/test/tests_on_odes/optim_test.jl +++ b/test/tests_on_odes/optim_test.jl @@ -6,10 +6,6 @@ obj = build_loss_objective( ### Optim Method -println("Use Optim Brent to fit the parameter") -result = Optim.optimize(obj, 1.0, 10.0) -@test_broken result.minimizer[1] ≈ 1.5 atol = 3.0e-1 - println("Use Optim BFGS to fit the parameter") result = Optim.optimize(obj, [1.0], Optim.BFGS()) @test result.minimizer[1] ≈ 1.5 atol = 3.0e-1 diff --git a/test/tests_on_odes/test_problems.jl b/test/tests_on_odes/test_problems.jl index 530b5fa..589fc1d 100644 --- a/test/tests_on_odes/test_problems.jl +++ b/test/tests_on_odes/test_problems.jl @@ -1,27 +1,27 @@ -using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools +using OrdinaryDiffEq, RecursiveArrayTools # Here are the problems to solve -f1 = @ode_def begin - dx = a * x - x * y - dy = -3y + x * y -end a +f1 = function (du, u, p, t) + du[1] = p[1] * u[1] - u[1] * u[2] + du[2] = -3 * u[2] + u[1] * u[2] +end u0 = [1.0; 1.0] tspan = (0.0, 10.0) p = [1.5] prob1 = ODEProblem(f1, u0, tspan, [1.5]) -f2 = @ode_def begin - dx = a * x - x * y - dy = -c * y + x * y -end a c +f2 = function (du, u, p, t) + du[1] = p[1] * u[1] - u[1] * u[2] + du[2] = -p[2] * u[2] + u[1] * u[2] +end p = [1.5, 3.0] prob2 = ODEProblem(f2, u0, tspan, p) -f3 = @ode_def begin - dx = a * x - b * x * y - dy = -c * y + d * x * y -end a b c d +f3 = function (du, u, p, t) + du[1] = p[1] * u[1] - p[2] * u[1] * u[2] + du[2] = -p[3] * u[2] + p[4] * u[1] * u[2] +end p = [1.5, 1.0, 3.0, 1.0] prob3 = ODEProblem(f3, u0, tspan, p)