insysbio / LikelihoodProfiler.jl

LikelihoodProfiler is a Julia package for practical identifiability analysis and confidence intervals evaluation.
https://insysbio.github.io/LikelihoodProfiler.jl/latest
MIT License
16 stars 4 forks source link

Error when `scan_bouns` and `theta_bounds` are identical #10

Open TorkelE opened 5 months ago

TorkelE commented 5 months ago

I am trying to perform likelihood profiling for a problem from parameter fitting to data. Here, I know bounds for all parameters, and want to incorporate this both into theta_bounds and scan_bounds. Turns out that if scan_bounds is identical to the bounds for that parameter in theta_bounds you get an error. Wouldn't it be natural to allow this?

In practice this is solvable by making scan_bounds ever so slightly narrower on both sides (and maybe this can be done internally if thus turns out to be the case).

Just a thought, that this might be a nice improvement to make.

Minimal example:

using LikelihoodProfiler

# testing profile function
g(x) = 5.0 + (x[1]-3.0)^2 + (x[1]-x[2]-1.0)^2 + 0*x[3]^2

# Calculate parameters intervals for first parameter component, x[1]
res_1 = get_interval(
  [3., 2., 2.1], # starting point
  1,             # parameter component to analyze
  g,             # profile function
  :LIN_EXTRAPOL; # method
  loss_crit = 9., # critical level of loss function
  theta_bounds = [(0.0,10.0), (0.0,10.0), (0.0,10.0)],
  scan_bounds = (0.0,10.0)
  )
#

gives:

ERROR: ArgumentError: scan_bound are outside of the theta_bounds (0.0, 10.0)
Stacktrace:
 [1] get_endpoint(theta_init::Vector{…}, theta_num::Int64, loss_func::typeof(g), method::Symbol, direction::Symbol; loss_crit::Float64, scale::Vector{…}, theta_bounds::Vector{…}, scan_bound::Float64, scan_tol::Float64, loss_tol::Float64, local_alg::Symbol, silent::Bool, kwargs::@Kwargs{})
   @ LikelihoodProfiler ~/.julia/packages/LikelihoodProfiler/Qi97K/src/get_endpoint.jl:145
 [2] get_endpoint
   @ ~/.julia/packages/LikelihoodProfiler/Qi97K/src/get_endpoint.jl:112 [inlined]
 [3] #27
   @ ./none:0 [inlined]
 [4] iterate(g::Base.Generator, s::Vararg{Any})
   @ Base ./generator.jl:47 [inlined]
 [5] collect(itr::Base.Generator{UnitRange{…}, LikelihoodProfiler.var"#27#29"{…}})
   @ Base ./array.jl:834
 [6] get_interval(theta_init::Vector{…}, theta_num::Int64, loss_func::typeof(g), method::Symbol; loss_crit::Float64, scale::Vector{…}, theta_bounds::Vector{…}, scan_bounds::Tuple{…}, scan_tol::Float64, loss_tol::Float64, local_alg::Symbol, kwargs::@Kwargs{})
   @ LikelihoodProfiler ~/.julia/packages/LikelihoodProfiler/Qi97K/src/get_interval.jl:119
 [7] top-level scope
   @ ~/Desktop/Julia Playground/Catalyst playground/Envirionment - PEtab_profile_liklihood_test/profile_test_playground.jl:157
Some type information was truncated. Use `show(err)` to see complete types.a NaN value in the state, parameter
ivborissov commented 5 months ago

Thanks for bringing this to our attention. Yes, it is natural to allow scan_bounds to be identical with theta_bounds, however we have seen certain use-cases, where it led to unexpected results due to the differences in the way scan_bounds are introduced in LikelihoodProfiler code and the way bounds (theta_bounds) are interpreted in NLopt. This question is part of a more general one: rewriting augmented loss function and moving away from NLopt to Optimization.jl.

TorkelE commented 5 months ago

Got it.

Yes, I think that if we could simply pass SciML's OptimizationFunction (or OptimizationProblem) into get_interval that would be very useful. E.g. stuff like theta_bounds is already defined there, and SciML/Optimization algorithms would be easier to use.

One could also pass a OptimizationSolution as the starting point.

I think it could produce a really neat interface in the end, but yes, it also requires some work to carry out the changes.