Merge pull request #73 from hetalang/opt_prob

adds OptimizationProblem interface

Author: ivborissov

.github/workflows/ci.yml

@@ -16,9 +16,6 @@ jobs:
       fail-fast: false
       matrix:
         include:
-          - {version: '1.9',  os: windows-latest, arch: x64}
-          - {version: '1.9',  os: ubuntu-latest,  arch: x64}
-          - {version: '1.9',  os: macos-latest,   arch: aarch64}
           - {version: 'lts', os: windows-latest, arch: x64}
           - {version: 'lts', os: ubuntu-latest,  arch: x64}
           - {version: 'lts', os: macos-latest,   arch: aarch64}

Project.toml

@@ -11,11 +11,11 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
-NLopt_jll = "079eb43e-fd8e-5478-9966-2cf3e3edb778"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@@ -31,18 +31,18 @@ CSV = "0.10"
 DataFrames = "1"
 DataStructures = "0.18"
 Distributions = "0.25"
+ForwardDiff = "0.10.38, 1.0.1"
 LabelledArrays = "1.6.1"
-NLopt = "1.0.3"
-NLopt_jll = "=2.8.0"
 NaNMath = "1"
+OptimizationNLopt = "0.2.2, 0.3.2"
 OrdinaryDiffEq = "6"
 ProgressMeter = "1"
 RecipesBase = "1"
 Reexport = "1"
 SciMLBase = "2.15"
 Sundials = "4.24"
 XLSX = "0.7, 0.8"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/basics/overview.md

@@ -564,8 +564,8 @@ fit_res = fit(p, to_fit)
 ```
 
 ```julia
-FitResult with status :XTOL_REACHED
-   Status: XTOL_REACHED
+FitResult with status :Success
+   Status: Success
    Optimal values: [:kabs => 18.868605026704916, :Q => 4.043662480774219, :kel => 0.17104243648378176, :sigma1 => 0.020347955494158528, :sigma2 => 0.31561050699802246, :sigma3 => 0.5716026958426483]
    OF value: 140.96503722972034
    OF count: 8612

docs/src/tutorial/fit-files/fit-run.jl

@@ -67,19 +67,17 @@ res_optim = fit(p, params_df)
 res_optim = fit(
     p, 
     params_df, 
-    fit_alg = :LN_SBPLX, 
+    fit_alg = NLopt.LN_SBPLX(), 
     ftol_abs = 1e-5, 
     ftol_rel = 0, 
-    maxeval = 10^6
+    maxiters = 10^6
 )
 optim(res_optim)
 
 #= see API docs
 fit_alg : fitting algorithm. Default is :LN_NELDERMEAD
-ftol_abs : absolute tolerance on function value. Default is 0.0
-ftol_rel : relative tolerance on function value. Default is 1e-4
-xtol_rel : relative tolerance on optimization parameters. Default is 0.0
-xtol_rel : absolute tolerance on optimization parameters. Default is 0.0
-maxeval : maximum number of function evaluations. Default is 1e4
+ftol_abs : absolute tolerance on objective value. Default is 0.0
+ftol_rel : relative tolerance on objective value. Default is 1e-4
+maxiters : maximum number of objective evaluations. Default is 1e4
 maxtime : maximum optimization time (in seconds). Default is 0
 =#

docs/src/tutorial/fit.md

@@ -172,8 +172,8 @@ res_optim = fit(p, to_fit) # default fitting
 ```
 ┌ Info: Scenario ":multiple_15" has no measurements. It will be excluded from fitting.
 └ @ HetaSimulator y:\HetaSimulator.jl\src\fit.jl:77
-FitResult with status :XTOL_REACHED
-   Status: XTOL_REACHED
+FitResult with status :Success
+   Status: Success
    Optimal values: [:kabs => 18.868605026704916, :Q => 4.043662480774219, :kel => 0.17104243648378176, :sigma1 => 0.020347955494158528, :sigma2 => 0.31561050699802246, :sigma3 => 0.5716026958426483]
    OF value: 140.96503722971997
    OF count: 8612
@@ -246,38 +246,36 @@ res_optim = fit(p, params_df)
 ```
 ┌ Info: Scenario ":multiple_15" has no measurements. It will be excluded from fitting.
 └ @ HetaSimulator 
-FitResult with status :FTOL_REACHED
-   Status: FTOL_REACHED
+FitResult with status :Success
+   Status: Success
    Optimal values: [:kabs => 8.669590504032879, :kel => 0.2299120380231296, :Q => 3.386457652767808, :sigma1 => 0.010105725225267037, :sigma2 => 0.09951673713071268, :sigma3 => 0.6024808584834973]
    OF value: -101.7645013649068
    OF count: 417
 ```
 
 ## Additional optimization-specific options
 
-Internally `HetaSimulator` uses `NLopt` library. One can choose the optimization algorithm as well as additional options.
+Internally, `HetaSimulator` uses the `Optimization.jl` interface, which supports various optimization packages. By default, `HetaSimulator` includes the NLopt library, allowing you to choose an NLopt-based optimization algorithm without installing additional packages. To access more optimizers, you can add any optimization package supported by `Optimization.jl` and use the preferred optimizer algorithm in `fit()`.
 
 Read more about NLopt algorithms choice: <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/>
 
 ```julia
 res_optim = fit(
     p, 
     params_df, 
-    fit_alg = :LN_SBPLX, 
+    fit_alg = NLopt.LN_SBPLX(), 
     ftol_abs = 1e-5, 
     ftol_rel = 0, 
-    maxeval = 10^6
+    maxiters = 10^6
 )
 optim(res_optim)
 ```
 
 There are several optimization related arguments, which are available for a user. 
 To learn more read about [`fit`](@ref) method in API documentation.
 
-- `fit_alg` : fitting algorithm. Default is `:LN_NELDERMEAD`
-- `ftol_abs` : absolute tolerance on function value. Default is `0.0`
-- `ftol_rel` : relative tolerance on function value. Default is `1e-4`
-- `xtol_rel` : relative tolerance on optimization parameters. Default is `0.0`
-- `xtol_rel` : absolute tolerance on optimization parameters. Default is `0.0`
-- `maxeval` : maximum number of function evaluations. Default is `1e4`
+- `fit_alg` : fitting algorithm. Default is `NLopt.LN_NELDERMEAD()`
+- `ftol_abs` : absolute tolerance on objective value. Default is `0.0`
+- `ftol_rel` : relative tolerance on objective value. Default is `1e-4`
+- `maxiters` : maximum number of objective evaluations. Default is `1e4`
 - `maxtime` : maximum optimization time (in seconds). Default is `0`

src/HetaSimulator.jl

@@ -28,9 +28,9 @@ module HetaSimulator
   @reexport using SciMLBase.EnsembleAnalysis
   @reexport using OrdinaryDiffEq
   using Sundials
-
+  using ForwardDiff
   # fitting
-  using NLopt
+  @reexport using OptimizationNLopt
 
   # utils
   using LabelledArrays
@@ -66,6 +66,7 @@ module HetaSimulator
   include("solution_interface.jl")
   include("plots.jl")
   include("loss.jl")
+  include("optprob.jl")
   include("fit.jl")
   include("estimator.jl")
   include("monte_carlo.jl")
@@ -86,7 +87,7 @@ module HetaSimulator
   export CVODE_BDF, CVODE_Adams
   export optim, obj
   export sim, mc, mc!
-  export fit, loss, estimator
+  export fit, loss, estimator, generate_optimization_problem
   export HetaSimulatorDir
   export update
   export times, vals, status, status_summary

src/estimator.jl

@@ -3,7 +3,7 @@
 # aux function to allow user write both `Pair{Symbol,Float64}` and `Pair{Symbol,Real}`
 # convert to regular `Vector{Pair{Symbol,Float64}}` for further processing
 normalize_params(v::AbstractVector{<:Pair{Symbol,<:Real}}) =
-    [k => Float64(val) for (k,val) in v]
+    [k => float(val) for (k,val) in v]
 
 """
     function estimator(
@@ -35,7 +35,7 @@ normalize_params(v::AbstractVector{<:Pair{Symbol,<:Real}}) =
 
   Returns:
 
-    function(x:Vector{Float64}=last.(parameters_fitted))
+    function(x=last.(parameters_fitted))
   
   The method returns anonimous function which depends on parameters vector in the same order as in `parameters_fitted`.
   This function is ready to be passed to optimizer routine or identifiability analysis.
@@ -90,13 +90,17 @@ function estimator(
     function(prob,i,repeat) # internal_prob_func
       #update_init_values(selected_prob[i], last(selected_scenario_pairs[i]).init_func, x)
       scn = last(selected_scenario_pairs[i])
-      params_total = merge_strict(scn.parameters, x)
+      T = eltype(x)
+      params_total = (T <: ForwardDiff.Dual) ? merge_strict(NamedTuple{keys(scn.parameters)}(T.(collect(scn.parameters))), x) : merge_strict(scn.parameters, x)
       remake_prob(selected_prob[i], scn.init_func, params_total; safetycopy=true) #?safetycopy=false
     end
   end
 
   function _output(sol, i)
-    sol.retcode != :Success && error("Simulated scenario $i returned $(sol.retcode) status")
+    if !SciMLBase.successful_retcode(sol.retcode)
+      @warn "Simulated scenario $i returned $(sol.retcode) status"
+      return (Inf, false)
+    end
     sv = sol.prob.kwargs[:callback].discrete_callbacks[1].affect!.saved_values
     simulation = Simulation(sv, NamedTuple(parameters_fitted_norm), sol.retcode)
     result = SimResult(simulation, last(selected_scenario_pairs[i]))
@@ -118,7 +122,7 @@ function estimator(
 
   ### function ready for fitting
 
-  function out(x::Vector{Float64}=last.(parameters_fitted_norm))
+  function out(x=last.(parameters_fitted_norm))
     x_nt = NamedTuple{Tuple(parameters_fitted_names)}(x)
     solution = solve(
       prob(x_nt),
@@ -234,3 +238,4 @@ function estimator(
 
   return estimator(scenario_pairs, parameters_fitted; kwargs...)
 end
+