search

tpapp / MultistartOptimization.jl

Multistart optimization methods in Julia.
Other
64 stars 10 forks source link

TaskFailedException when solving a problem depending on an ODE #30

Closed DanielVandH closed 2 years ago

DanielVandH commented 2 years ago

Whenever I try to apply these methods to a problem which depends on an ODE, I get a TaskFailedException. Could this be related to #5? An example is shown below, where I also commented out some code which removes a for loop to see if it relates to #5 directly, but it shows no difference. Maybe there's a weird interaction between the integrator interface and the interface used here.

using MultistartOptimization
using NLopt
using DifferentialEquations
using Random
using Distributions

Random.seed!(2992999)
λ = -0.5
y₀ = 15.0
σ = 0.1
T = 5.0
n = 200
Δt = T / n
t = [j * Δt for j in 0:n]
y = y₀ * exp.(λ * t)
yᵒ = y .+ [0.0, rand(Normal(0, σ), n)...]
function ode_fnc(u, p, t)
    λ = p
    du = λ * u
    return du
end
function loglik(θ, data, integrator)
    yᵒ, n = data
    λ, σ, u0 = θ
    integrator.p = λ
    reinit!(integrator, u0)
    solve!(integrator)
    ℓℓ = -0.5n * log(2π * σ^2)
    for i in eachindex(yᵒ)
        ℓℓ = ℓℓ - 0.5 / σ^2 * (yᵒ[i] - integrator.sol.u[i])^2
    end
    #ℓℓ = -0.5n * log(2π * σ^2) - 0.5 / σ^2 * sum((yᵒ - integrator.sol.u).^2)
end

θ₀ = [-0.5, 0.1, 15.0]
data = (yᵒ, n)
lb = [1e-12, -10.0, -10.0, -10.0, -10.0]
ub = [10.0, 10.0, 10.0, 10.0, 10.0]
ode_problem = ODEProblem(ode_fnc, y₀, (0.0, T), -0.5)
integrator = init(ode_problem, Tsit5(), saveat = t)

P = MinimizationProblem(θ -> loglik(θ, data, integrator), lb, ub)
local_method = NLoptLocalMethod(NLopt.LN_NELDERMEAD)
multistart_method = TikTak(100)
p = multistart_minimization(multistart_method, local_method, P)
ERROR: TaskFailedException
Stacktrace:
  [1] wait
    @ .\task.jl:345 [inlined]
  [2] fetch
    @ .\task.jl:360 [inlined]
  [3] iterate
    @ .\generator.jl:47 [inlined]
  [4] _collect(c::Vector{Task}, itr::Base.Generator{Vector{Task}, typeof(fetch)}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:807
  [5] collect_similar(cont::Vector{Task}, itr::Base.Generator{Vector{Task}, typeof(fetch)})
    @ Base .\array.jl:716
  [6] map(f::Function, A::Vector{Task})
    @ Base .\abstractarray.jl:2931
  [7] sobol_starting_points(minimization_problem::MinimizationProblem{var"#7#8", Vector{Float64}}, N::Int64, use_threads::Bool)
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:62
  [8] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#7#8", Vector{Float64}}; use_threads::Bool)   
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:88
  [9] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#7#8", Vector{Float64}})
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:85
 [10] top-level scope
    @ c:\Users\licer\.julia\dev\ProfileLikelihood\test\refinement.jl:82

    nested task error: BoundsError: attempt to access 201-element Vector{Float64} at index [138]
    Stacktrace:
     [1] getindex
       @ .\array.jl:924 [inlined]
     [2] solution_endpoint_match_cur_integrator!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})   
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:139
     [3] _postamble!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, 
OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, 
Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:123
     [4] postamble!
       @ C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:120 [inlined]
     [5] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, 
Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\solve.jl:485
     [6] loglik(θ::Vector{Float64}, data::Tuple{Vector{Float64}, Int64}, integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ Main c:\Users\licer\.julia\dev\ProfileLikelihood\test\refinement.jl:64
     [7] (::var"#7#8")(θ::Vector{Float64})
       @ Main c:\Users\licer\.julia\dev\ProfileLikelihood\test\refinement.jl:79
     [8] _initial
       @ C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:60 [inlined]
     [9] (::MultistartOptimization.var"#3#6"{Vector{Float64}, MultistartOptimization.var"#_initial#4"{var"#7#8"}})()
       @ MultistartOptimization .\threadingconstructs.jl:258
tpapp commented 2 years ago

To me it looks like an error within loglik. I would debug this by simply calling it at the relevant θ, eg

loglik(θ, data, integrator)

You can print the θ first thing in the function to see where it fails.

DanielVandH commented 2 years ago

That's what I had initially assumed, but I wasn't able to identify anything to do with loglik. Here is a modified form which doesn't actually use the integrator result (also fixed the bounds, they had extra entries by accident).

using MultistartOptimization
using NLopt
using DifferentialEquations
using Random
using Distributions

Random.seed!(2992999)
λ = -0.5
y₀ = 15.0
σ = 0.1
T = 5.0
n = 200
Δt = T / n
t = [j * Δt for j in 0:n]
y = y₀ * exp.(λ * t)
yᵒ = y .+ [0.0, rand(Normal(0, σ), n)...]
function ode_fnc(u, p, t)
    λ = p
    du = λ * u
    return du
end
function loglik(θ, data, integrator)
    solve!(integrator)
    2.0
end

θ₀ = [-0.5, 0.1, 15.0]
data = (yᵒ, n)
lb = [1e-12, 1e-12, 1e-12]
ub = [10.0, 10.0, 10.0]
ode_problem = ODEProblem(ode_fnc, y₀, (0.0, T), -0.5)
integrator = init(ode_problem, Tsit5(), saveat = t)

P = MinimizationProblem(θ -> loglik(θ, data, integrator), lb, ub)
local_method = NLoptLocalMethod(NLopt.LN_NELDERMEAD)
multistart_method = TikTak(100)
p = multistart_minimization(multistart_method, local_method, P)

This leads to the same error:

ERROR: TaskFailedException
Stacktrace:
  [1] wait
    @ .\task.jl:345 [inlined]
  [2] fetch
    @ .\task.jl:360 [inlined]
  [3] iterate
    @ .\generator.jl:47 [inlined]
  [4] _collect(c::Vector{Task}, itr::Base.Generator{Vector{Task}, typeof(fetch)}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:807
  [5] collect_similar(cont::Vector{Task}, itr::Base.Generator{Vector{Task}, typeof(fetch)})
    @ Base .\array.jl:716
  [6] map(f::Function, A::Vector{Task})
    @ Base .\abstractarray.jl:2931
  [7] sobol_starting_points(minimization_problem::MinimizationProblem{var"#123#124", Vector{Float64}}, N::Int64, use_threads::Bool)
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:62
  [8] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#123#124", Vector{Float64}}; use_threads::Bool)
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:88
  [9] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#123#124", Vector{Float64}})
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:85
 [10] top-level scope
    @ c:\Users\licer\.julia\dev\ProfileLikelihood\dev\rvm.jl:39

    nested task error: BoundsError: attempt to access 200-element Vector{Float64} at index [202]
    Stacktrace:
     [1] getindex
       @ .\array.jl:924 [inlined]
     [2] solution_endpoint_match_cur_integrator!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:139
     [3] _postamble!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, 
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:123
     [4] postamble!
       @ C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\integrators\integrator_utils.jl:120 [inlined]
     [5] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\UG9Mz\src\solve.jl:485
     [6] loglik(θ::Vector{Float64}, data::Tuple{Vector{Float64}, Int64}, integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, false, Float64, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Float64}, ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, Float64, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 
Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, ODEFunction{false, typeof(ode_fnc), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, 
DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, Float64, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
       @ Main c:\Users\licer\.julia\dev\ProfileLikelihood\dev\rvm.jl:25
     [7] (::var"#123#124")(θ::Vector{Float64})
       @ Main c:\Users\licer\.julia\dev\ProfileLikelihood\dev\rvm.jl:36
     [8] _initial
       @ C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:60 [inlined]
     [9] (::MultistartOptimization.var"#3#6"{Vector{Float64}, MultistartOptimization.var"#_initial#4"{var"#123#124"}})()
       @ MultistartOptimization .\threadingconstructs.jl:258

That error message for the nested task error, BoundsError: attempt to access 200-element Vector{Float64} at index [202], seems to change the number of elements pretty often. Interestingly, if I add @show θ to the loglik function, it works:

function loglik(θ, data, integrator)
    @show θ
    solve!(integrator)
    2.0
end
...
(location = [0.07812500000099219, 6.640625000000336, 5.546875000000445], value = 2.0)

It also doesn't break if I simply used solve(ode_problem) rather than work with the integrator interface:

function ode_fnc(u, p, t)
    λ = p
    du = λ * u
    return du
end
function loglik(θ, data, integrator)
    @show θ
    solve(integrator)
    2.0
end

θ₀ = [-0.5, 0.1, 15.0]
data = (yᵒ, n)
lb = [1e-12, 1e-12, 1e-12]
ub = [10.0, 10.0, 10.0]
ode_problem = ODEProblem(ode_fnc, y₀, (0.0, T), -0.5)
integrator = init(ode_problem, Tsit5(), saveat = t)

P = MinimizationProblem(θ -> loglik(θ, data, ode_problem), lb, ub)
local_method = NLoptLocalMethod(NLopt.LN_NELDERMEAD)
multistart_method = TikTak(100)
p = multistart_minimization(multistart_method, local_method, P)

Finally, returnig to the original loglik, adding @show θ leads to a MethodError with isless:

function loglik(θ, data, integrator)
    @show θ
    yᵒ, n = data
    λ, σ, u0 = θ
    integrator.p = λ
    reinit!(integrator, u0)
    solve!(integrator)
    ℓℓ = -0.5n * log(2π * σ^2)
    for i in eachindex(yᵒ)
        ℓℓ = ℓℓ - 0.5 / σ^2 * (yᵒ[i] - integrator.sol.u[i])^2
    end
    #ℓℓ = -0.5n * log(2π * σ^2) - 0.5 / σ^2 * sum((yᵒ - integrator.sol.u).^2)
end

θ₀ = [-0.5, 0.1, 15.0]
data = (yᵒ, n)
lb = [1e-12, 1e-12, 1e-12]
ub = [10.0, 10.0, 10.0]
ode_problem = ODEProblem(ode_fnc, y₀, (0.0, T), -0.5)
integrator = init(ode_problem, Tsit5(), saveat = t)

P = MinimizationProblem(θ -> loglik(θ, data, integrator), lb, ub)
local_method = NLoptLocalMethod(NLopt.LN_NELDERMEAD)
multistart_method = TikTak(100)
p = multistart_minimization(multistart_method, local_method, P)
ERROR: MethodError: no method matching isless(::Nothing, ::Nothing)
Closest candidates are:
  isless(::Any, ::Missing) at missing.jl:88
  isless(::Missing, ::Any) at missing.jl:87
Stacktrace:
  [1] lt(o::Base.Order.ForwardOrdering, a::Nothing, b::Nothing)
    @ Base.Order .\ordering.jl:117
  [2] lt(o::Base.Order.By{MultistartOptimization.var"#7#8", Base.Order.ForwardOrdering}, a::NamedTuple{(:location, :value), Tuple{Vector{Float64}, Nothing}}, b::NamedTuple{(:location, :value), Tuple{Vector{Float64}, Nothing}})
    @ Base.Order .\ordering.jl:119
  [3] selectpivot!
    @ .\sort.jl:533 [inlined]
  [4] partition!(v::Vector{NamedTuple{(:location, :value), Tuple{Vector{Float64}, Nothing}}}, lo::Int64, hi::Int64, o::Base.Order.By{MultistartOptimization.var"#7#8", Base.Order.ForwardOrdering})
    @ Base.Sort .\sort.jl:555
  [5] sort!
    @ .\sort.jl:633 [inlined]
  [6] partialsort!
    @ .\sort.jl:98 [inlined]
  [7] #partialsort!#2
    @ .\sort.jl:157 [inlined]
  [8] partialsort(v::Vector{NamedTuple{(:location, :value), Tuple{Vector{Float64}, Nothing}}}, k::UnitRange{Int64}; kws::Base.Pairs{Symbol, MultistartOptimization.var"#7#8", Tuple{Symbol}, NamedTuple{(:by,), Tuple{MultistartOptimization.var"#7#8"}}})
    @ Base.Sort .\sort.jl:166
  [9] _keep_lowest(xs::Vector{NamedTuple{(:location, :value), Tuple{Vector{Float64}, Nothing}}}, N::Int64)
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:75
 [10] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#143#144", Vector{Float64}}; use_threads::Bool)
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:90
 [11] multistart_minimization(multistart_method::TikTak, local_method::NLoptLocalMethod{Set{Symbol}}, minimization_problem::MinimizationProblem{var"#143#144", Vector{Float64}})
    @ MultistartOptimization C:\Users\licer\.julia\packages\MultistartOptimization\vzoMo\src\tiktak.jl:85
 [12] top-level scope
    @ c:\Users\licer\.julia\dev\ProfileLikelihood\dev\rvm.jl:46

Maybe this last one is a separate issue? If it is (and if this isn't just me making a mistake) I could open an issue for this error.

This seems really weird. I'm not sure why adding @show θ would have anything to do with making it work? Any ideas?

DanielVandH commented 2 years ago

Something I didn't realise I could do is to set use_threads = false. If I do this, the TaskFailedException goes away, although that isless error returns. So I guess there's a thread-safe issue just because I'm using some mutables in the objective function. Probably up to the user (me) to be responsible for locking that / using use_threads = false.

Adding @show removes the isless error though. I'm clueless here.

(For the TaskFailedException, based on the error message I imagine the issue is that the integrator is being accessed by separate threads at the same time, and the adaptive timestepping in solving the ODE makes the vectors different sizes, hence why it presents as a BoundsError.)

tpapp commented 2 years ago

First, I would try to get code running with use_threads = false.

I think that you get nothings in the initial points. But I don't understand how the code posted above can do that. Can you do a 

@show xs

in _keep_lowest and post the output here?

DanielVandH commented 2 years ago

Ahh... my loglik,

function loglik(θ, data, integrator)
    yᵒ, n = data
    λ, σ, u0 = θ
    integrator.p = λ
    reinit!(integrator, u0)
    solve!(integrator)
    ℓℓ = -0.5n * log(2π * σ^2)
    for i in eachindex(yᵒ)
        ℓℓ = ℓℓ - 0.5 / σ^2 * (yᵒ[i] - integrator.sol.u[i])^2
    end
    #ℓℓ = -0.5n * log(2π * σ^2) - 0.5 / σ^2 * sum((yᵒ - integrator.sol.u).^2)
end

wasn't actually returning anything. All the initial objective values were nothing. No wonder it worked when I used @show - it made it be the returned value. Very sorry for that dumb mistake..

So, no problems with the isless error - that's on me. The TaskFailedException still persists as given in the original question (with use_threads = true), but as mentioned maybe it's on me to make sure the objective function is thread-safe (maybe that's worth documenting, at least?). This can probably be closed.

tpapp commented 2 years ago

I think the error message could be improved --- thanks for reporting and tracking down the issue. I am closing it for now and made a note about it in #31. If you still experience problems, please feel free to comment here and I will reopen.