ODE Parameter Estimation. Simple Local Optimization Example

New to Julia
,
1

Hi,

I am using Julia 1.4.0 in Atom

I am trying to reproduce the results in the Docs Simple Local Optimization

The first example optimizing the the Lotka-Volterra equation with one parameter works fine, but then I cannot reproduce the next step:

“Lastly, we can use the same tools to estimate multiple parameters simultaneously. Let’s use the Lotka-Volterra equation with all parameters free:”

I have used the following code

# Simple Local Optimization
# https://docs.sciml.ai/stable/analysis/parameter_estimation/
using DifferentialEquations, DiffEqParamEstim, Plots, Optim

# We choose to optimize the parameters on the Lotka-Volterra equation. We do so
# by defining the function as a ParameterizedFunction:
function f2(du,u,p,t)
   du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
   du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f2,u0,tspan,p)

# We create data using the numerical result with a=1.5 = p[1]:
sol = solve(prob,Tsit5())
t = collect(range(0,stop=10,length=200))
using RecursiveArrayTools # for VectorOfArray
noiseAmp = 0.1 # 0.01
randomized = VectorOfArray([(sol(t[i]) + noiseAmp*randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
# Here we used VectorOfArray from RecursiveArrayTools.jl to turn the result of
# an ODE into a matrix.
# If we plot the solution with the parameter at a=1.42, we get the following:
sol2 = solve(ODEProblem(f,u0,tspan,[1.4,1.2,3.1,1.1]),Tsit5())
dat2 = convert(Array, VectorOfArray([sol2(t[i]) for i in 1:length(t)]))
scatter(t,data[1,:])
scatter!(t,data[2,:])
plot!(t,dat2[1,:])
plot!(t,dat2[2,:])

# Notice that after one period this solution begins to drift very far off: this
# problem is sensitive to the choice of a.
#
# To build the objective function for Optim.jl, we simply call the
#   build_loss_objective function:
# We can build an objective function and solve the multiple parameter version just as before:
cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)

# Before optimizing, let's visualize our cost function by plotting it for a
# range of parameter values:
vals = 0.0:0.1:10.0
# using Plots; plotly()
plot(vals,[cost_function(i) for i in vals],yscale=:log10,
     xaxis = "Parameter", yaxis = "Cost", title = "1-Parameter Cost Function",
     lw = 3)

# Here we see that there is a very well-defined minimum in our cost function at
# the real parameter (because this is where the solution almost exactly fits the
# dataset).
#
# Now this cost function can be used with Optim.jl in order to get the
# parameters. For example, we can use Brent's algorithm to search for the best
# solution on the interval [0,10] by:
result = optimize(cost_function, 0.0, 10.0)

# This returns result.minimizer[1]==1.5 as the best parameter to match the data.
# When we plot the fitted equation on the data, we receive the following:

sol3 = solve(ODEProblem(f,u0,tspan,[result.minimizer[1]]),Tsit5())
dat3 = convert(Array, VectorOfArray([sol3(t[i]) for i in 1:length(t)]))
scatter(t,data[1,:])
scatter!(t,data[2,:])
plot!(t,dat3[1,:])
plot!(t,dat3[2,:])

# Thus we see that after fitting, the lines match up with the generated data and
# receive the right parameter value.

Everything works fine until Julia tries to plot the cost_function and gives the error:


LoadError: BoundsError
Stacktrace:
 [1] getindex at ./number.jl:78 [inlined]
 [2] f2 at /home/s/Dropbox/Code/Julia/Tutorials/DifferentialEquations/Simple Local Optimization2.jl:8 [inlined]
 [3] ODEFunction at /home/s/.julia/packages/DiffEqBase/bWm8j/src/diffeqfunction.jl:248 [inlined]
 [4] initialize!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{Float64,1},Nothing,Float64,Float64,Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Tuple{},Array{Float64,1},Tuple{}},Array{Float64,1},Float64,Nothing,OrdinaryDiffEq.DefaultInit}, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}) at /home/s/.julia/packages/OrdinaryDiffEq/RYCSN/src/perform_step/low_order_rk_perform_step.jl:623
 [5] __init(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Type{Val{true}}; saveat::Array{Float64,1}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/s/.julia/packages/OrdinaryDiffEq/RYCSN/src/solve.jl:403
 [6] #__solve#358 at /home/s/.julia/packages/OrdinaryDiffEq/RYCSN/src/solve.jl:4 [inlined]
 [7] solve_call(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Float64,ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Tsit5; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:saveat, :save_everystep, :dense, :maxiters, :verbose),Tuple{Array{Float64,1},Bool,Bool,Int64,Bool}}}) at /home/s/.julia/packages/DiffEqBase/bWm8j/src/solve.jl:60
 [8] #solve#443 at /home/s/.julia/packages/DiffEqBase/bWm8j/src/solve.jl:78 [inlined]
 [9] (::DiffEqParamEstim.var"#43#48"{Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Integer,Tuple{Symbol,Symbol},NamedTuple{(:maxiters, :verbose),Tuple{Int64,Bool}}},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},Array{Float64,2},Nothing,Nothing,Nothing},Nothing})(::Float64) at /home/s/.julia/packages/DiffEqParamEstim/TLBTz/src/build_loss_objective.jl:40
 [10] DiffEqObjective at /home/s/.julia/packages/DiffEqParamEstim/TLBTz/src/build_loss_objective.jl:24 [inlined]
 [11] iterate at ./generator.jl:47 [inlined]
 [12] collect(::Base.Generator{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},DiffEqObjective{DiffEqParamEstim.var"#43#48"{Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Integer,Tuple{Symbol,Symbol},NamedTuple{(:maxiters, :verbose),Tuple{Int64,Bool}}},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},Array{Float64,2},Nothing,Nothing,Nothing},Nothing},DiffEqParamEstim.var"#47#53"{DiffEqParamEstim.var"#43#48"{Nothing,Bool,Int64,typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR),Base.Iterators.Pairs{Symbol,Integer,Tuple{Symbol,Symbol},NamedTuple{(:maxiters, :verbose),Tuple{Int64,Bool}}},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f2),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,L2Loss{Array{Float64,1},Array{Float64,2},Nothing,Nothing,Nothing},Nothing}}}}) at ./array.jl:665
 [13] top-level scope at /home/s/Dropbox/Code/Julia/Tutorials/DifferentialEquations/Simple Local Optimization2.jl:47
in expression starting at /home/s/Dropbox/Code/Julia/Tutorials/DifferentialEquations/Simple Local Optimization2.jl:47

julia>

I am probably doing a very basic mistake in the code but I cannot find where it is.

Thanks.

2

It expects all parameter values as a vector but you are giving a scalar to cost_function because prob uses f2 opposed to linked document it uses f with one parameter.

3

Thanks very much tomakluftu for your prompt response. But I still do not know how to fix the problem.

in the code the parameters are given by,

p = [1.5,1.0,3.0,1.0]

and this is used by prob (I think).

Also prob is using f2 (and no f) because I just copied the code in the documentation.
However you have to scroll down a little bit from the linked document to get to the relevant part that starts at the quoted text:

" Lastly, we can use the same tools to estimate multiple parameters simultaneously. Let’s use the Lotka-Volterra equation with all parameters free:"

Just below that line in the document is where f2 is defined.
(Sorry, I do not know how to link directly to the exact position in the document).

function f2(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f2,u0,tspan,p)

And just after the definition of “prob” says
“We can build an objective function and solve the multiple parameter version just as before:”
And there is the code defining cost_function that I just also copied and paste

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                      maxiters=10000,verbose=false)

How should I change the call to prob or to cost_function?

Thanks again

4

Here, the plot cannot be produced because there would be 4 parameters plus cost value(5D plot?) and there is no such plot functionality available. You can just skip this part / comment out.

5

Thanks Lutfullah,

I got it now. The key was, as you said at the beginning, that I was passing a scalar to the cost function, and it was expecting a vector.