Merge pull request #39 from insysbio/single_par

fix single parameter issue

Author: ivborissov

src/odeprob_utils.jl

@@ -25,7 +25,7 @@ function solver_init(sciml_prob::SciMLBase.AbstractODEProblem,
     function affect!(integrator)
       set_x_fixed!(solver_state_opt.reinit_cache.p, integrator.u[idx])
       solver_state_opt.reinit_cache.u0 = integrator.u[1:end-1][1:end .!= idx]
-      sol = solve!(solver_state_opt)
+      sol = solve_optcache(solver_state_opt)
       for i in 1:length(integrator.u)-1
         i == idx && continue
         integrator.u[i] = sol[i - (i>idx)]

src/optprob_utils.jl

@@ -59,7 +59,7 @@ function build_optprob_reduced(optprob::OptimizationProblem, optpars)
   else
     ub_reduced = zeros(eltype(ub_full), lenreduced)
   end
-
+  
   u0_full = optprob.u0
   u0_reduced = zeros(eltype(u0_full), lenreduced)
   remake(optprob, f=optf_reduced, u0=u0_reduced, p=FixedParamCache(optprob.p, 1, u0_full[1], 1.0), lb=lb_reduced, ub=ub_reduced)

src/problem_interface.jl

@@ -111,7 +111,6 @@ function validate_profile_bound(bound)
 end
 
 function validate_optpars(npars::Int, u)
-  npars < 2 && throw(DimensionMismatch("`optpars` length must be ≥ 2."))
   !(u isa AbstractVector && npars == length(u)) && 
     throw(ArgumentError("OptimizationProblem initial values must be of the same type and size as `optpars`."))
   return nothing

src/profiler_state.jl

@@ -155,9 +155,14 @@ end
 
 function profiler_step!(profiler::ProfilerState, method::OptimizationProfiler)
   stepper = get_stepper(method)
+  optcache = get_solver_state(profiler)
+  idx = get_idx(profiler)
 
   pars_guess = compute_next_pars!(profiler, stepper)
-  sol = solve_optcache(profiler)
+  fill_x_reduced!(optcache.reinit_cache.u0, pars_guess, idx)
+  set_x_fixed!(optcache.reinit_cache.p, pars_guess[idx])
+  
+  sol = solve_optcache(optcache)
 
   if SciMLBase.successful_retcode(sol.retcode)
     idx = get_idx(profiler)
@@ -173,15 +178,23 @@ function profiler_step!(profiler::ProfilerState, method::OptimizationProfiler)
 
 end
 
-function solve_optcache(profiler::ProfilerState)
-  optcache = get_solver_state(profiler)
-  pars_guess = get_parscache(profiler)
-  idx = get_idx(profiler)
+function solve_optcache(optcache::OptimizationCache)  
+  if isempty(optcache.reinit_cache.u0)
+    return solve_empty_optcache(optcache)
+  else
+    return solve!(optcache)
+  end
+end
 
-  fill_x_reduced!(optcache.reinit_cache.u0, pars_guess, idx)
-  set_x_fixed!(optcache.reinit_cache.p, pars_guess[idx])
+function solve_empty_optcache(optcache::OptimizationCache)
+  u = optcache.reinit_cache.u0
+  p = optcache.reinit_cache.p
+  
+  t = @elapsed obj = optcache.f(u, p)
+  stats = SciMLBase.OptimizationStats(; iterations = 1, time = t, fevals = 1)
 
-  return solve!(optcache)
+  return SciMLBase.build_solution(optcache, optcache.opt, u, obj; 
+    retcode = ReturnCode.Success, stats)
 end
 
 
@@ -200,7 +213,7 @@ function profiler_step!(profiler::ProfilerState, method::IntegrationProfiler)
     profiler.numiter += 1
   else
     @warn "Solver returned $(integrator.sol.retcode) retcode at profile point x = $(get_curx(profiler)). Profiling is interrupted."
-    profiler.solver_retcode = sol.retcode
+    profiler.solver_retcode = integrator.sol.retcode
   end
 
 end

test/problems/analytic_funcs.jl

@@ -1,7 +1,7 @@
 
 funcs_dict = Dict(
 
-#=
+
   :f_1p => Dict(
     :func => (x,p) -> 5.0 + (x[1]-3.0)^2,
     :optim => [3.],
@@ -10,7 +10,7 @@ funcs_dict = Dict(
     :profile_range => [(-20.,20.)],
     :retcode => [(:Identifiable,:Identifiable)]
   ),
-=#
+
   :f_1p_ex => Dict(
     :func => (x,p) -> 5.0 + (x[1]-1e-8)^2 + 0.0*x[2], 
     :optim => [1e-8, 2.],

test/test_prob_interface.jl

@@ -3,7 +3,6 @@ using LikelihoodProfiler, Test
 # ParameterProfile tests
 optprob1 = OptimizationProblem((x,p)->x[1]^2+x[2]^2, [1., 2.0])
 
-@test_throws DimensionMismatch PLProblem(optprob1, [0.])
 @test_throws ArgumentError PLProblem(optprob1, [0.,0]) 
 @test_throws ArgumentError PLProblem(optprob1, [0.,0, 0.]) 
 @test_throws ArgumentError PLProblem(optprob1, [0.,0], (-Inf,4)) 

test/test_profile.jl

@@ -1,15 +1,16 @@
-using LikelihoodProfiler, Test, OptimizationNLopt, Plots, OrdinaryDiffEq, ForwardDiff
+using LikelihoodProfiler, Test
+using Optimization, OptimizationNLopt, Plots, OrdinaryDiffEq, ForwardDiff
 
 ######################################### PLProblem ##########################################
 
-optf = OptimizationFunction((x,p) -> 5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2, Optimization.AutoForwardDiff())
-optprob = OptimizationProblem(optf, [0.,0.])
+optf = OptimizationFunction((x,p) -> 5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2, AutoForwardDiff())
+optprob = OptimizationProblem(optf, [0.,0.]; lb=[-3,-1], ub=[10,9])
 
 plprob = PLProblem(optprob, [3.,4.], [(-5,20), (-2,15)]; threshold=4.0)
 
 #################################### OptimizationProfiler ####################################
 
-method1 = OptimizationProfiler(optimizer=NLopt.LN_NELDERMEAD(), step=FixedStep(; initial_step=0.1))
+method1 = OptimizationProfiler(optimizer=NLopt.LN_NELDERMEAD(), stepper = FixedStep(; initial_step=0.1))
 sol1 = profile(plprob, method1; verbose=true)
 
 
@@ -22,3 +23,26 @@ sol2 = profile(plprob, method2; verbose=true)
 
 method3 = CICOProfiler(optimizer=:LN_NELDERMEAD, scan_tol=1e-3)
 sol3 = profile(plprob, method3; verbose=true)
+
+
+######################################### PLProblem w 1 parameter ##########################################
+
+optf = OptimizationFunction((x,p) -> 5.0 + (x[1]-3.0)^2, AutoForwardDiff())
+optprob = OptimizationProblem(optf, [0.]; lb=[-10.], ub=[50.])
+
+plprob = PLProblem(optprob, [3.], (-5,20); threshold=4.0)
+
+#################################### OptimizationProfiler ####################################
+
+method1 = OptimizationProfiler(optimizer=Optimization.LBFGS(), stepper = FixedStep(; initial_step=0.1))
+sol1 = profile(plprob, method1; verbose=true)
+
+#################################### IntegrationProfiler ####################################
+
+method2 = IntegrationProfiler(integrator=FBDF(autodiff = AutoFiniteDiff()), integrator_opts=(dtmax=0.1,), matrix_type=:hessian)
+sol2 = profile(plprob, method2; verbose=true)
+
+######################################## CICOProfiler #######################################
+
+method3 = CICOProfiler(optimizer=:LN_NELDERMEAD, scan_tol=1e-3)
+sol3 = profile(plprob, method3; verbose=true)