Fitting ODE system to one single data

New to Julia
,
1

Hi everyone,
I have an ODE system with 5 equations and I want to fit it to a data. However, the data is related to only one of the equations.

I am trying to use the functions in Parameters Estimation of Differential equations, but all the examples it has are fitting N equations to N-column data
https://docs.sciml.ai/stable/analysis/parameter_estimation/#parameter_estimation-1

here is the system

I am not sure if the function is fitting the first eq of the system with the data (deaths) that I have, because I am getting some terrible results. (dont forget that “deaths” is a vector)

Anyone knows if this is correct? Or if it is not, how could I use those functions to estimate the parameter, making sure it is fitting the data only to the first equation of the system?

I thank you all in advance,
Thomas

code here

function f(du,u,p,t)
Θ = (ip.ρ/(1-ip.ρ))*(ip.μ+ip.γ)
    #Θ = (p[2]/(1-p[2]))*(ip.μ+ip.γ)
    σ = isolation(t)
    du[1] = Θ*u[3]     ##here is the eq. I want the results
    du[5] = ip.ν-ip.μ*u[5]-σ*p[1]*u[5]*u[3]
    du[2] = σ*p[1]*u[5]*u[3]-ip.α*u[2]-ip.μ*u[2]
    du[3] = ip.α*u[2]-(ip.μ+Θ+ip.γ)*u[3]
    du[4] = ip.γ*u[3]-ip.μ*u[4]
end
u0 = [0.0;0.0;1.0;0.0;ip.pop_est]
tspan = (0.0,150.0)
prob1 = ODEProblem(f,u0,tspan,p)
#here is the optmization function.
cost_function = build_loss_objective(prob1,Tsit5(),L2Loss(time_d,deaths),maxiters=10000,verbose=false)
result = optimize(cost_function,lower,upper,p)
p = result.minimizer
2

Use save_idxs, I’d recommend you use DiffEqFlux.jl which is a much more flexible parameter estimation system.

3

Thanks for your answer.
The selection of the equation worked. Now I am trying to use the package you suggested.
I copied the example in
https://diffeqflux.sciml.ai/dev/examples/LV-ODE/
and I am getting a error. Sorry to bother you

result_ode = DiffEqFlux.sciml_train(loss_adjoint, p, BFGS(initial_stepnorm = 0.0001),cb = callback)
ERROR: ArgumentError: tuple must be non-empty
Stacktrace:
 [1] first(::Tuple{}) at ./tuple.jl:77
 [2] _unapply(::Nothing, ::Tuple{}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:151
 [3] _unapply(::Tuple{Nothing}, ::Tuple{}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:155 (repeats 2 times)
 [4] _unapply(::Tuple{NTuple{6,Nothing},Tuple{Nothing}}, ::Tuple{Nothing,Nothing,Nothing,Nothing,Array{Float64,1},Array{Float64,1}}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:156
 [5] unapply(::Tuple{NTuple{6,Nothing},Tuple{Nothing}}, ::Tuple{Nothing,Nothing,Nothing,Nothing,Array{Float64,1},Array{Float64,1}}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:165
 [6] (::Zygote.var"#172#173"{DiffEqBase.var"#512#back#459"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#93"{Tsit5,InterpolatingAdjoint{0,true,Val{:central},Bool},Array{Float64,1},Array{Float64,1},Tuple{}}},Tuple{NTuple{6,Nothing},Tuple{Nothing}}})(::Array{Float64,2}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:172
 [7] (::Zygote.var"#337#back#174"{Zygote.var"#172#173"{DiffEqBase.var"#512#back#459"{DiffEqSensitivity.var"#adjoint_sensitivity_backpass#93"{Tsit5,InterpolatingAdjoint{0,true,Val{:central},Bool},Array{Float64,1},Array{Float64,1},Tuple{}}},Tuple{NTuple{6,Nothing},Tuple{Nothing}}}})(::Array{Float64,2}) at /home/thomas/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [8] #solve#445 at /home/thomas/.julia/packages/DiffEqBase/KnYSY/src/solve.jl:69 [inlined]
 [9] (::typeof(∂(#solve#445)))(::Array{Float64,2}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [10] (::Zygote.var"#172#173"{typeof(∂(#solve#445)),Tuple{NTuple{6,Nothing},Tuple{Nothing}}})(::Array{Float64,2}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:171
 [11] (::Zygote.var"#337#back#174"{Zygote.var"#172#173"{typeof(∂(#solve#445)),Tuple{NTuple{6,Nothing},Tuple{Nothing}}}})(::Array{Float64,2}) at /home/thomas/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [12] #solve at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0 [inlined]
 [13] (::typeof(∂(#solve)))(::Array{Float64,2}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [14] predict_adjoint at ./REPL[24]:2 [inlined]
 [15] (::typeof(∂(predict_adjoint)))(::Array{Float64,2}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [16] loss_adjoint at ./REPL[26]:2 [inlined]
 [17] (::typeof(∂(loss_adjoint)))(::Tuple{Int64,Nothing}) at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [18] #172 at /home/thomas/.julia/packages/Zygote/YeCEW/src/lib/lib.jl:171 [inlined]
 [19] (::Zygote.var"#337#back#174"{Zygote.var"#172#173"{typeof(∂(loss_adjoint)),Tuple{Tuple{Nothing},Int64}}})(::Tuple{Int64,Nothing}) at /home/thomas/.julia/packages/ZygoteRules/6nssF/src/adjoint.jl:49
 [20] #34 at /home/thomas/.julia/packages/DiffEqFlux/6oM6z/src/train.jl:176 [inlined]
 [21] (::typeof(∂(λ)))(::Int64) at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface2.jl:0
 [22] #38 at /home/thomas/.julia/packages/Zygote/YeCEW/src/compiler/interface.jl:46 [inlined]
 [23] (::DiffEqFlux.var"#37#50"{DiffEqFlux.var"#34#47"{typeof(loss_adjoint)}})(::Array{Float64,1}, ::Array{Float64,1}) at /home/thomas/.julia/packages/DiffEqFlux/6oM6z/src/train.jl:199
 [24] value_gradient!!(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}) at /home/thomas/.julia/packages/NLSolversBase/mGaJg/src/interface.jl:82
 [25] initial_state(::BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Float64,Flat}, ::Optim.Options{Float64,DiffEqFlux.var"#_cb#46"{var"#7#8",Base.Iterators.Cycle{Tuple{DiffEqFlux.NullData}}}}, ::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}) at /home/thomas/.julia/packages/Optim/L5T76/src/multivariate/solvers/first_order/bfgs.jl:66
 [26] optimize(::TwiceDifferentiable{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Float64,Flat}, ::Optim.Options{Float64,DiffEqFlux.var"#_cb#46"{var"#7#8",Base.Iterators.Cycle{Tuple{DiffEqFlux.NullData}}}}) at /home/thomas/.julia/packages/Optim/L5T76/src/multivariate/optimize/optimize.jl:33
 [27] #sciml_train#33(::Function, ::Int64, ::DiffEqFlux.ZygoteDiffMode, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqFlux.sciml_train), ::Function, ::Array{Float64,1}, ::BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Float64,Flat}, ::Base.Iterators.Cycle{Tuple{DiffEqFlux.NullData}}) at /home/thomas/.julia/packages/DiffEqFlux/6oM6z/src/train.jl:269
 [28] #sciml_train at ./none:0 [inlined] (repeats 2 times)
 [29] top-level scope at REPL[32]:1
4

Are you on the latest versions? Just to double check, I ran it on the latest releases and it was fine.

5

Worked!!!

yes, the problem was the julia version.
Thank you!
I will try to apply in my problem now

1 Like
6

Hi, which version of Julia are you using? I am having the same error with Julia v1.4 and DiffEqFlux v1.12.0.

1 Like
7

What version of Zygote and DiffEqSensitivity?

8

Zygote v0.4.22 and DiffEqSensitivity v6.14.1

9

Could you update DiffEqSensitivity? You’re before save_idxs was fixed (which was quite recent).

10

It works now after changing DiffEqSensitivity from 6.14.1 to 6.22.0. Thank you, @ChrisRackauckas