A component of the DiffEq ecosystem for enabling sensitivity analysis for scientific machine learning (SciML). Optimize-then-discretize, discretize-then-optimize, adjoint methods, and more for ODEs, SDEs, DDEs, DAEs, etc.
Home Page: https://docs.sciml.ai/SciMLSensitivity/stable/
License: Other
SciMLSensitivity.jl is a component package in the SciML Scientific Machine Learning ecosystem. It holds the sensitivity analysis utilities. Users interested in using this functionality should check out DifferentialEquations.jl.
For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.
![dependabot[bot] avatar](https://avatars.githubusercontent.com/in/29110?v=4 "dependabot[bot]")
![femtocleaner[bot] avatar](https://avatars.githubusercontent.com/in/4123?v=4 "femtocleaner[bot]")
![github-actions[bot] avatar](https://avatars.githubusercontent.com/in/15368?v=4 "github-actions[bot]")
@JuliaRegistrator register()
Another approach to address the ill-conditioning of the medium-sized parameter sensitivity matrix is the likelihood profiler method proposed here: insysbio/LikelihoodProfiler.jl#3 (comment)
I am not sure if I missed this, but if so I couldn't find anything in the docs: Is it possible to compute derivatives of the solution with respect to initial values? This is crucial for shooting methods.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3730677/pdf/rsif20121018.pdf
I think the nicest thing to do would be to return a GLM regression model along with the PRCC coefficients.
Implement a sensitivity(::ODEProblem,tspan) which applies Morris global sensitivity analysis on a ParameterizedFunction.
http://www.abe.ufl.edu/Faculty/jjones/ABE_5646/2010/Morris.1991%20SA%20paper.pdf
http://www.math.ualberta.ca/~hwang/sensitivity
http://discovery.ucl.ac.uk/19896/1/19896.pdf
It seems to me that the means and variances are calculated over the whole simulation time and not over the elementary effects for each time point.
If I run the example code from the documentation (http://docs.juliadiffeq.org/latest/analysis/global_sensitivity.html):
function f(du,u,p,t)
du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
du[2] = -3*u[2] + u[1]*u[2]
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0]
prob = ODEProblem(f,u0,tspan,p)
t = collect(range(0, stop=10, length=200))
m = DiffEqSensitivity.morris_sensitivity(prob,Tsit5(),t,[[1,5],[0.5,5]],[10,10],len_trajectory=1500,total_num_trajectory=1000,num_trajectory=150)
I would expect an output array of length 200 for each parameter and output combination (as written in the documentation). But on my installation I will get very lengthy arrays instead.
Elementary effects are correct:
size.(m.elementary_effects[1])
gives
107759-element Array{Tuple{Int64,Int64},1}:
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
⋮
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
(2, 200)
but
size.(m.means)
gives
2-element Array{Tuple{Int64,Int64,Int64},1}:
(2, 1, 107759)
(2, 1, 117091)
107759+ 117091 = 224850 is exactly the number of samples (len_trajectory*num_trajectory-150 = 1500×150−150 = 224850)!
If I calculate the means and variances by hand (just changing dims from 2 to 3) I will get following:
means = []
variances = []
for k in m.elementary_effects
push!(means,mean(convert(Array, VectorOfArray(k)),dims=3))
push!(variances,var(convert(Array, VectorOfArray(k)),dims=3))
end
means[1]
which gives
2×200×1 Array{Float64,3}:
[:, :, 1] =
0.0 0.0513992 0.106989 0.169601 … 7.24055 7.67254 7.95 7.98303
0.0 0.0011582 0.00429882 0.00909004 3.84639 3.7428 3.6888 3.68657
which seems to be correct.
Please, could someone check this! Maybe this could also be a problem on my local installation.
Thanks a lot!
It would be nice to have a backsolve function which is like solve but has the adjoints already setup with Zygote.jl, where it runs a choice gradient calculation when Zygote is used. We can set it up so that way it uses forward-mode automatically when it makes sense, switches to adjoint sensitivity analysis in the cases that are applicable (i.e. currently ODEs), and then falls back to pure Zygote reverse-mode in the case where length(p) >> length(u0). Of course, then there can just be a mode overload or something like that so you can pass down all of the requested SensitivityAlg args to override the choice, and then we should get this setup in the DiffEqParamEstim/DiffEqBayes to make all of the gradients efficient.
@YingboMa had an idea to use save points from the forward solution to reset the backsolve. We should make the current BacksolveAdjoint always checkpoint like this.
It seems to me that the effects array should be preallocated by the number of parameters e.g. effects = [[] for i=1:length(p_steps)] and not effects = [[]]. If not, elem_effect with change_index != length(effects)-1 will be assigned to the wrong parameter array.
In this case
if length(effects) >= change_index && change_index > 0 push!(effects[change_index],elem_effect) elseif change_index > 0 push!(effects,[elem_effect]) end
can also simplified to
if change_index > 0 push!(effects[change_index],elem_effect) end
These are all methods necessary for getting gradients to fit chaotic systems.
https://www.sciencedirect.com/science/article/pii/S0021999117304783
https://arxiv.org/abs/1807.05568
https://arxiv.org/abs/1902.11112
http://discovery.ucl.ac.uk/19896/1/19896.pdf
I was trying to go through the documentation for adjoint sensitivity analysis.
In addition to having to fix some function signatures that were incorrect in the documentation, the resulting sensitivities seem to be wrong when compared to autodifferentiation and finite differences.
I am following the documentation with this code:
`using DifferentialEquations
using DiffEqSensitivity
function f(du,u,p,t)
du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end
p = [1.5,1.0,3.0]
prob = ODEProblem(f,[1.0; 1.0],(0.0, 10.0),p)
sol = solve(prob,Vern9(),abstol=1e-10,reltol=1e-10);
dg(out,u,p,t,i) = (out .= 1.0 .- u) # this function signature is wrong/incomplete in the documentation
t = [0.5, 1.0] # this definition is missing from the documentation
res = adjoint_sensitivities(sol,Vern9(),dg,t,abstol=1e-14,
reltol=1e-14,iabstol=1e-14,ireltol=1e-12)`
which results in
1×3 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}: 0.616953 -0.357961 0.100422
Again, copying the code from the documentation to check against other methods:
`
using ForwardDiff, Calculus
function G(p)
tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
sol = solve(tmp_prob,Vern9(), abstol=1e-14, reltol=1e-14, saveat=t)
A = convert(Array,sol)
sum(((1 .- A).^2)./2)
end
@show G([1.5,1.0,3.0])
@show res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
@show res3 = Calculus.gradient(G,[1.5,1.0,3.0]);
`
which results in the entirely different (but consistent) gradients
G([1.5, 1.0, 3.0]) = 2.161316471053881 res2 = ForwardDiff.gradient(G, [1.5, 1.0, 3.0]) = [4.698722374012723, -2.227674237734399, 1.1820220043862988] res3 = Calculus.gradient(G, [1.5, 1.0, 3.0]) = [4.698722373289014, -2.227674241910004, 1.1820220061302733]
Posting here after asking on Gitter.
I was trying to make a SecondOrderODEProblem work in a neural ODE layer. This is what I have so far:
using DifferentialEquations
using Flux
using DiffEqFlux
model = Chain(Dense(2, 50, tanh), Dense(50, 2))
u0 = Float32[0.; 2.]
du0 = Float32[0.; 0.]
tspan = (0.0f0, 1.0f0)
t = range(tspan[1], tspan[2], length=20)
p1 = param(Flux.data(DiffEqFlux.destructure(model)))
ps = Flux.params(p1)
function f(du, u, p, t)
println("u is ", u, " du is ", du)
return DiffEqFlux.restructure(model, p)(u)
end
prob = SecondOrderODEProblem{false}(f, du0, u0, tspan, p1)
function predict_adjoint()
diffeq_adjoint(p1, prob, Tsit5(), u0=(du0, u0), saveat=t, abstol=1e-8, reltol=1e-6)
end
correct_pos = Float32.(transpose(hcat(collect(0:0.05:1)[2:end], collect(2:-0.05:1)[2:end])))
loss_n_ode() = sum(abs2, correct_pos .- predict_adjoint()[1:2, :])
data = Iterators.repeated((), 1000)
opt = ADAM(0.1)
cb = function ()
println(loss_n_ode())
loss_n_ode() < 0.01 && Flux.stop()
end
Flux.train!(loss_n_ode, ps, data, opt, cb=cb)The ODE seems to run okay, as in running diffeq_adjoint works. But then when I try to run the Flux.train! line I get errors about zero and similar not being defined for the relevant argument types, and if I try to fix these I get ERROR: type DynamicalODEFunction has no field jac.
It seems that the neural ODE layer is not implemented for this ODE type. It has been pointed out that diffeq_rd should work in this case, so I'll give it a go.
Hello, when I type
using DiffEqSensitivityreturns
[ Info: Precompiling DiffEqSensitivity [41bf760c-e81c-5289-8e54-58b1f1f8abe2]
ERROR: LoadError: LoadError: syntax: invalid function name "Morris <: GSAMethod"
Stacktrace:
[1] include at ./boot.jl:317 [inlined]
[2] include_relative(::Module, ::String) at ./loading.jl:1044
[3] include at ./sysimg.jl:29 [inlined]
[4] include(::String) at /home/jarr/.julia/packages/DiffEqSensitivity/xlGSs/src/DiffEqSensitivity.jl:3
[5] top-level scope at none:0
[6] include at ./boot.jl:317 [inlined]
[7] include_relative(::Module, ::String) at ./loading.jl:1044
[8] include(::Module, ::String) at ./sysimg.jl:29
[9] top-level scope at none:2
[10] eval at ./boot.jl:319 [inlined]
[11] eval(::Expr) at ./client.jl:393
[12] top-level scope at ./none:3
in expression starting at /home/jarr/.julia/packages/DiffEqSensitivity/xlGSs/src/morris_sensitivity.jl:1
in expression starting at /home/jarr/.julia/packages/DiffEqSensitivity/xlGSs/src/DiffEqSensitivity.jl:17
ERROR: Failed to precompile DiffEqSensitivity [41bf760c-e81c-5289-8e54-58b1f1f8abe2] to /home/jarr/.julia/compiled/v1.0/DiffEqSensitivity/02xYn.ji.
Stacktrace:
[1] error(::String) at ./error.jl:33
[2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1203
[3] _require(::Base.PkgId) at ./loading.jl:960
[4] require(::Base.PkgId) at ./loading.jl:858
[5] require(::Module, ::Symbol) at ./loading.jl:853
I think it is better to pass in f_p, g_p and g_u explicitly for adjoint sensitivity analysis instead of relying on the parameter p in f(du,u,p,t) and dg(out,u,p,t). Because p can be used for other purposes like caching and storing configurations which are not for differentiation.
It needs to be a little different since it won't have a forward solution first, so we will want to drop down to a helper function and handle it there.
It's only used for linear regression, e.g. https://github.com/JuliaDiffEq/DiffEqSensitivity.jl/blob/a565c73581d88b29b84eee6b5ee0feb0eb1faa04/src/regression_sensitivity.jl#L113 and, in this application, it would be simpler just to use \. A nice side effect would be that you avoid Meta.parse and eval.
There's still a lot to do in terms of simplifying the implementations. Essentially things like Jv and v'J should just call a helper function in SparseDiffTools. auto_jacvec and num_jacvec already exist. We should add some for v'J as well, maybe with Tracker, ReverseDiff, and Zygote for benchmarking. Then we just call those helper functions here, with dropdowns to generating Jacobians if necessary.
@JuliaRegistrator register()
When the running-cost function, g, depends explicitly on the system parameters, p. I get a mismatch between the gradients reported by adjoint_sensitivities and Calculus.gradient.
I suspect this discrepancy is due to the fact that I calculate dg as the derivative of g w.r.t. only the state, and omit its derivatives w.r.t. the parameters, p.
Below, I provide a minimal example. The dynamics is that of a simple pendulum, controlled at the upward equilibrium point.
An observation: if I make the function g independent of the parameters, then the derivatives do match.
using DifferentialEquations.OrdinaryDiffEq
using DiffEqSensitivity
using DiffEqBase
using Calculus
using ForwardDiff
using QuadGK
function pendulum_eom(dx, x, p, t)
dx[1] = x[2]
dx[2] = -sin(x[1]) + (-p[1]*sin(x[1]) + p[2]*x[2]) # Second term is a simple controller that stabilizes π
end
x0 = [0.1, 0.0]
tspan = (0.0, 10.0)
p = [-24.05, -19.137]
prob = ODEProblem(pendulum_eom, x0, tspan, p)
sol = solve(prob, Vern9(), abstol=1e-8, reltol=1e-8)
g(x, p, t) = 1.0*(x[1] - π)^2 + 1.0*x[2]^2 + 5.0*(-p[1]*sin(x[1]) + p[2]*x[2])^2
dg(out, y, p, t) = ForwardDiff.gradient!(out, x -> g(x, p, t), y)
res = adjoint_sensitivities(sol,Vern9(),g,nothing,dg,abstol=1e-8,
reltol=1e-8,iabstol=1e-8,ireltol=1e-8)
function G(p)
tmp_prob = remake(prob,p=p)
sol = solve(tmp_prob,Vern9(),abstol=1e-8,reltol=1e-8)
res,err = quadgk((t)-> g(sol(t), p, t), 0.0,10.0,atol=1e-8,rtol=1e-8)
res
end
res2 = Calculus.gradient(G,p)
res' ≈ res2
If the ODE function contains any terms involving matrix multiplications (dense or sparse), adjoint sensitivity analysis fails. It appears that the issue is somewhere in the automatic Jacobian calculations.
Consider the following minimal example for a discrete objective:
A = rand(10, 10)
function f(du, u, p, t)
du .= A*u
end
tspan = (0.0, 1.0)
u0 = rand(10)
p = []
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5());
function dg(out, u, p, t, i)
out .= 2(1 .- u)
end
tts = [0.1, 0.2, 0.5]
adjoint_sensitivities(sol, Vern9(), dg, tts)
which fails with the following error message:
MethodError: *(::Transpose{Float64,Array{Float64,2}}, ::TrackedArray{…,Array{Float64,1}}) is ambiguous. Candidates:
*(x::AbstractArray{T,2} where T, y::Tracker.TrackedArray{T,1,A} where A where T) in Tracker at /home/henrik/.julia/packages/Tracker/6wcYJ/src/lib/array.jl:368
*(transA::Transpose{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:84
Possible fix, define
*(::Transpose{#s623,#s622} where #s622<:AbstractArray{T,2} where #s623, ::Tracker.TrackedArray{S,1,A} where A)
Stacktrace:
[1] (::getfield(Tracker, Symbol("##487#488")){Array{Float64,2},TrackedArray{…,Array{Float64,1}}})(::TrackedArray{…,Array{Float64,1}}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/lib/array.jl:376
[2] back_(::Tracker.Grads, ::Tracker.Call{getfield(Tracker, Symbol("##487#488")){Array{Float64,2},TrackedArray{…,Array{Float64,1}}},Tuple{Nothing,Tracker.Tracked{Array{Float64,1}}}}, ::TrackedArray{…,Array{Float64,1}}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:110
[3] back(::Tracker.Grads, ::Tracker.Tracked{Array{Float64,1}}, ::TrackedArray{…,Array{Float64,1}}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:123
[4] #16 at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:113 [inlined]
[5] foreach at ./abstractarray.jl:1867 [inlined]
[6] back_(::Tracker.Grads, ::Tracker.Call{getfield(Tracker, Symbol("##364#366")){TrackedArray{…,Array{Float64,1}},Tuple{Int64}},Tuple{Tracker.Tracked{Array{Float64,1}},Nothing}}, ::Float64) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:113
[7] back(::Tracker.Grads, ::Tracker.Tracked{Float64}, ::Float64) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:125
[8] (::getfield(Tracker, Symbol("##362#363")){Tracker.Grads})(::Tracker.Tracked{Float64}, ::Float64) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/lib/real.jl:156
[9] foreach(::Function, ::Array{Tracker.Tracked{Float64},1}, ::Array{Float64,1}) at ./abstractarray.jl:1867
[10] back_(::Tracker.Grads, ::Tracker.Call{typeof(Tracker.collect),Tuple{Array{Tracker.Tracked{Float64},1}}}, ::Array{Float64,1}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/lib/real.jl:156
[11] back(::Tracker.Grads, ::Tracker.Tracked{Array{Float64,1}}, ::Array{Float64,1}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:125
[12] #18 at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:140 [inlined]
[13] (::getfield(Tracker, Symbol("##21#23")){getfield(Tracker, Symbol("##18#19")){Tracker.Params,TrackedArray{…,Array{Float64,1}}}})(::Array{Float64,1}) at /home/henrik/.julia/packages/Tracker/6wcYJ/src/back.jl:149
[14] (::DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing})(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Any,1}, ::Float64) at /home/henrik/.julia/packages/DiffEqSensitivity/DI6VG/src/adjoint_sensitivity.jl:120
[15] ODEFunction at /home/henrik/.julia/packages/DiffEqBase/EFqMn/src/diffeqfunction.jl:188 [inlined]
[16] ode_determine_initdt(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.ODEIntegrator{Vern9,true,Array{Float64,1},Float64,Array{Any,1},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,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem},Vern9,OrdinaryDiffEq.InterpolationData{ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}) at /home/henrik/.julia/packages/OrdinaryDiffEq/uJP6R/src/initdt.jl:23
[17] auto_dt_reset! at /home/henrik/.julia/packages/OrdinaryDiffEq/uJP6R/src/integrators/integrator_interface.jl:238 [inlined]
[18] handle_dt!(::OrdinaryDiffEq.ODEIntegrator{Vern9,true,Array{Float64,1},Float64,Array{Any,1},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,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem},Vern9,OrdinaryDiffEq.InterpolationData{ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}}}},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Vern9Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Vern9ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},DataStructures.BinaryHeap{Float64,DataStructures.GreaterThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1},Float64}) at /home/henrik/.julia/packages/OrdinaryDiffEq/uJP6R/src/solve.jl:370
[19] #__init#285(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Nothing, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Float64, ::Float64, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Float64, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::Vern9, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/henrik/.julia/packages/OrdinaryDiffEq/uJP6R/src/solve.jl:333
[20] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:abstol, :reltol, :save_everystep, :save_start),Tuple{Float64,Float64,Bool,Bool}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::Vern9, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0
[21] #__solve#284(::Base.Iterators.Pairs{Symbol,Real,NTuple{4,Symbol},NamedTuple{(:abstol, :reltol, :save_everystep, :save_start),Tuple{Float64,Float64,Bool,Bool}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::Vern9, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/henrik/.julia/packages/OrdinaryDiffEq/uJP6R/src/solve.jl:6
[22] #__solve at ./none:0 [inlined] (repeats 5 times)
[23] #solve#429(::Base.Iterators.Pairs{Symbol,Real,NTuple{4,Symbol},NamedTuple{(:abstol, :reltol, :save_everystep, :save_start),Tuple{Float64,Float64,Bool,Bool}}}, ::Function, ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::Vern9) at /home/henrik/.julia/packages/DiffEqBase/EFqMn/src/solve.jl:39
[24] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:abstol, :reltol, :save_everystep, :save_start),Tuple{Float64,Float64,Bool,Bool}}, ::typeof(solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Any,1},ODEFunction{true,DiffEqSensitivity.ODEAdjointSensitivityFunction{Array{Float64,1},Array{Float64,1},Array{Float64,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,SensitivityAlg{0,true,Val{:central}},Nothing,UniformScaling{Bool},Nothing,Nothing,Nothing,Array{Float64,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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}},Nothing},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},CallbackSet{Tuple{},Tuple{DiscreteCallback{getfield(DiffEqCallbacks, Symbol("##31#34")){Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}},getfield(DiffEqCallbacks, Symbol("##33#36")){typeof(DiffEqBase.INITIALIZE_DEFAULT),Bool,getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},Base.RefValue{Union{Nothing, Float64}},getfield(DiffEqCallbacks, Symbol("##32#35")){getfield(DiffEqSensitivity, Symbol("#time_choice#20")){Array{Float64,1},Base.RefValue{Int64}},getfield(DiffEqSensitivity, Symbol("##19#21")){typeof(dg),Bool,Array{Float64,1},Array{Float64,1},Base.RefValue{Int64},Int64},Base.RefValue{Union{Nothing, Float64}}}}}}},DiffEqBase.StandardODEProblem}, ::Vern9) at ./none:0
[25] #adjoint_sensitivities#25(::Float64, ::Float64, ::Float64, ::Float64, ::SensitivityAlg{0,true,Val{:central}}, ::Array{Float64,1}, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}}, ::Vern9, ::typeof(dg), ::Array{Float64,1}, ::Nothing) at /home/henrik/.julia/packages/DiffEqSensitivity/DI6VG/src/adjoint_sensitivity.jl:332
[26] adjoint_sensitivities(::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,Array{Any,1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),UniformScaling{Bool},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}}}}, ::Vern9, ::Function, ::Array{Float64,1}, ::Nothing) at /home/henrik/.julia/packages/DiffEqSensitivity/DI6VG/src/adjoint_sensitivity.jl:330 (repeats 2 times)
[27] top-level scope at In[8]:21
The error occurs whenever there is any sort of matrix multiplication in the ODE function, it is not restricted to linear ODEs.
@JuliaRegistrator register()
Just like the new Zygote @adjoint, we can define how to lower solve to different choices of forward mode and adjoint definitions. I might need help from @oxinabox here.
As exemplified in the following code, using ODELocalSensitivityProblem results in errors about missing fields that are present in ODEFunction. This code does run if switched to the ODEProblem counterpart.
using DifferentialEquations
function lorenz!(du, u, p, t)
du[1] = p[1]*(u[2]-u[1])
du[2] = u[1]*(p[2]-u[3]) - u[2]
du[3] = u[1]*u[2] - p[3]*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0, 100.0)
params = [10.0, 28.0, 8.0/3.0]
probS = ODELocalSensitivityProblem(lorenz!, u0, tspan, params)
prob = ODEProblem(lorenz!, u0, tspan, params)
sol = solve(probS, Rodas4())
print(sol)Gives:
ERROR: LoadError: type ODELocalSensitvityFunction has no field tgrad
Stacktrace:
[1] getproperty(::Any, ::Symbol) at ./sysimg.jl:18
[2] has_tgrad(::ODELocalSensitvityFunction{true,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Tuple{Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1},Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1}},ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1}}},SensitivityAlg{0,true,DataType,true},Array{Float64,1},Nothing,Array{Float64,2},LinearAlgebra.UniformScaling{Bool}}) at /Users/asm/.julia/packages/DiffEqBase/8usQ9/src/diffeqfunction.jl:473
[3] build_grad_config(::Rodas4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}, ::Function, ::Function, ::Array{Float64,1}, ::Float64) at /Users/asm/.julia/packages/OrdinaryDiffEq/7mKtS/src/derivative_wrappers.jl:70
[4] alg_cache(::Rodas4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}, ::Array{Float64,1}, ::Array{Float64,1}, ::Type, ::Type, ::Type, ::Array{Float64,1}, ::Array{Float64,1}, ::ODELocalSensitvityFunction{true,ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Nothing,DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Tuple{Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1},Array{ForwardDiff.Dual{:___jac_tag,Float64,1},1}},ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3,Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.ParamJacobianWrapper{ODEFunction{true,typeof(lorenz!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,3},1}}},SensitivityAlg{0,true,DataType,true},Array{Float64,1},Nothing,Array{Float64,2},LinearAlgebra.UniformScaling{Bool}}, ::Float64, ::Float64, ::Float64, ::Array{Float64,1}, ::Bool, ::Type{Val{true}}) at /Users/asm/.julia/packages/OrdinaryDiffEq/7mKtS/src/caches/rosenbrock_caches.jl:630
...
My actual problem is stiff and compatible with autodiff, but not ParameterizedFunctions, hence not using the macros available there. Thanks for putting together such an amazing set of tools!
@JuliaRegistrator register()
@JuliaRegistrator register()
See this https://royalsocietypublishing.org/doi/full/10.1098/rsif.2017.0871
A new class of functions, called the ‘information sensitivity functions’ (ISFs), which quantify the information gain about the parameters through the measurements/observables of a dynamical system are presented. These functions can be easily computed through classical sensitivity functions alone and are based on Bayesian and information-theoretic approaches. While marginal information gain is quantified by decrease in differential entropy, correlations between arbitrary sets of parameters are assessed through mutual information.
For individual parameters, these information gains are also presented as marginal posterior variances, and, to assess the effect of correlations, as conditional variances when other parameters are given. The easy to interpret ISFs can be used to (a) identify time intervals or regions in dynamical system behaviour where information about the parameters is concentrated; (b) assess the effect of measurement noise on the information gain for the parameters; (c) assess whether sufficient information in an experimental protocol (input, measurements and their frequency) is available to identify the parameters; (d) assess correlation in the posterior distribution of the parameters to identify the sets of parameters that are likely to be indistinguishable; and (e) assess identifiability problems for particular sets of parameters.
To be more efficient, the InterpolatingAdjoint should just be creating a new dense sol between checkpoints and interpolating that, instead of doing a full solve to each point. By only saving the solution object, this would also massively cut down on compile times.
@JuliaRegistrator register()
Using toy example from tutorial:
function f(du,u,p,t)
du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
du[2] = -3*u[2] + u[1]*u[2]
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0]
prob = ODEProblem(f,u0,tspan,p)
t = collect(range(0, stop=10, length=200))
regre_sensitivity = regression_sensitivity(prob,Tsit5(),t,[[1,5],[0.5,5]],[3,3],100;coeffs=:rank)MethodError: no method matching -(::Array{Float64,1}, ::Float64)
Closest candidates are:
-(!Matched::Float64, ::Float64) at float.jl:397
-(!Matched::Complex{Bool}, ::Real) at complex.jl:298
-(!Matched::Missing, ::Number) at missing.jl:93
...
Stacktrace:
[1] pcs_and_srcs(::VectorOfArray{Float64,3,Array{Array{Float64,N} where N,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64) at /opt/julia/packages/DiffEqSensitivity/bB5NZ/src/regression_sensitivity.jl:85
[2] #regression_sensitivity#48(::Symbol, ::Function, ::getfield(DiffEqSensitivity, Symbol("##56#57")){ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},Tsit5,Array{Float64,1}}, ::Array{Array{Float64,1},1}, ::Array{Int64,1}, ::Int64) at /opt/julia/packages/DiffEqSensitivity/bB5NZ/src/regression_sensitivity.jl:49
[3] #regression_sensitivity at ./none:0 [inlined]
[4] #regression_sensitivity#55 at /opt/julia/packages/DiffEqSensitivity/bB5NZ/src/regression_sensitivity.jl:156 [inlined]
[5] (::getfield(DiffEqSensitivity, Symbol("#kw##regression_sensitivity")))(::NamedTuple{(:coeffs,),Tuple{Symbol}}, ::typeof(regression_sensitivity), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Float64,1}, ::Array{Array{Float64,1},1}, ::Array{Int64,1}, ::Int64) at ./none:0
I expect I did this well.
It does not work when replacing t by tspan, or by a vector [1.0, 3.0] for example. It does not crash but returns empty result if t is replaced by a unique number vector ie [5.0].
@JuliaRegistrator register()
@JuliaRegistrator register()
In Sundials.jl there's a way to get parameter sensitivities for ODEs with discontinuities, but when I tried the same approach with ODELocalSensitivityProblem, there either doesn't seem to be a way to do it, or it's undocumented and I didn't find it. The obvious thing didn't work:
using DifferentialEquations
f = @ode_def FFF begin
dx = -a
dy = b
end a=>2.0 b=>1.0
cb = ContinuousCallback((t,u,i) -> u[1], (args...)->(println("Stopped.");f.b=0.0))
prob = ODELocalSensitivityProblem(f, [1.0, 0.0], (0.0, 1.0), callback=cb)
sol = solve(prob, DP8())
display(sol)
The solution is
[-1.0, 0.5, -1.0, 0.0, 0.0, 1.0]
and the derivative dy/da is reported as zero, but should be -0.25. This is an issue whenever the callback changes some state (here b) and the trigger function depends on the state.
Using extract_local_sensitivities in combination with save_idxs fails as it is based on sol.prob.f.numindvar i.e. it always wants to extract the full state vector. For further discussion and a minimal example, see here.
In the implementation for the Morris method, it looks like the absolute value of the elementary effects are being taken:
Then, the mean and variance of this distribution are taken.
According to Global Sensitivity Analysis by Saltelli, and the R package reference for the sensitivity analysis, this is the correct calculation for
I propose that the elementary effects are stored without having their absolute value taken, and the MorrisResult is restructured to:
struct MorrisResult{T1,T2}
means::T1
means_star::T1
variances::T1
elementary_effects::T2
end
where means_star is the current means, means is calculated without taking the absolute value, and elementary_effects are stored without taking the absolute value.
Currently, the representation is kind of hard to work with. Maybe we should have a solution type that stores sensitivities as a matrix, and stores the ODE solution separately, e.g.
struct ODELocalSensitivitySolution{S<:AbstractODESolution,D<:AbstractMatrix}
sol::T
sen::D
endI am solving the 1D KS equation using FDM. I transformed the PDE into a set of ODEs inside a function that I called DerSen.jl. I used RK4 to solve the ODEs and the solution I got is correct.
using DifferentialEquations
p=c # parameters
u0=U[:,1] #Float
tspan=(0.0,Tspan)
prob=ODEProblem(DerSen,u0,tspan,p)
sol = solve(prob,RK4(),saveat=Δt)
U=convert(Array,sol)
t_vec=sol.tI would like to go further and compute the adjoint sensitivity of this model. My code is as follows:
using Sundials, DiffEqBase
g(u,p,t) = (sum(u).^2) ./ 2
function dg(out,u,p,t)
for i=1:Nx
out[i]=sum(U[i,:])
end
end
res = adjoint_sensitivities(sol,Vern9(),g,nothing,dg,abstol=1e-8,
reltol=1e-8,iabstol=1e-8,ireltol=1e-8)and I get this error massage that I don't understand
LoadError: MethodError: no method matching similar(::Float64)
Closest candidates are:
similar(!Matched::JuliaInterpreter.Compiled, !Matched::Any) at C:\Users*.juliapro\JuliaPro_v1.2.0-1\packages\JuliaInterpreter\MXq3U\src\types.jl:7
similar(!Matched::DataFrames.StackedVector, !Matched::Type, !Matched::Union{Integer, AbstractUnitRange}...) at C:\Users*.juliapro\JuliaPro_v1.2.0-1\packages\DataFrames\yH0f6\src\abstractdataframe\reshape.jl:401
similar(!Matched::Sundials.NVector) at C:\Users**.juliapro\JuliaPro_v1.2.0-1\packages\Sundials\fVIue\src\nvector_wrapper.jl:69
...
If the ODEProblem has the field tgrad, and we construct a ODEAdjointProblem from it, it seems that the resulting problem does not have the field tgrad. Can you add support to it?
I'd like to prototype DiffEqFlux with Zygote, which will have a much nicer API and generally improve all the things. However, to do that we need to make DiffEqSensitivity a but more flexible in terms of what kinds of parameters we can integrate over.
Currently, we only support a fixed-length vector. The key generalisation we need is to make this vector growable (which should be mathematically sensible, just act as if new elements were 0 in all previous time steps). As long as we have this it's easy to linearise any complex nested data structure, even when the structure changes over time.
Growable vectors are the simplest way to support Zygote but may not be the best; e.g. perhaps DiffEq should support nested key-value pairs with key paths as indices directly. In any case, I'm happy to help out with the higher-level parts but will need some expertise on the DiffEq side for the lower-level changes.
Just simple wrappers over the reverse mode AD libraries to plug into the whole Zygote system
Hi,
Running sobol_sensitivity in DiffEqSensitivity gives some negative values - even while useing the equation ind documents.
function f(du,u,p,t)
du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
du[2] = -3*u[2] + u[1]*u[2]
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0]
prob = ODEProblem(f,u0,tspan,p)
t = collect(range(0, stop=10, length=200))
s0 = sobol_sensitivity(prob,Tsit5(),t,[[1,5],[0.5,5]],N,0)
Out[8]: 2-element Array{Array{Float64,2},1}:
[NaN 0.507831 … 1.00731 1.00436; NaN 1.92336 … 0.732384 0.730945]
[NaN 0.47214 … 0.676224 0.681525; NaN -1.68656 … 0.879557 0.877603]
s1 = sobol_sensitivity(prob,Tsit5(),t,[[1,5],[0.5,5]],N,1)
Out[9]: 2-element Array{Array{Float64,2},1}:
[NaN 0.39537 … 0.341697 0.343645; NaN -2.06101 … 0.10922 0.106976]
[NaN 0.652815 … 0.00910675 0.00815206; NaN 5.24832 … 0.296978 0.296639]
s2 = sobol_sensitivity(prob,Tsit5(),t,[[1,5],[0.5,5]],N,2)
Out[10]: 1-element Array{Array{Float64,2},1}:
[NaN -0.0596478 … 0.652303 0.657847; NaN -1.84504 … 0.645139 0.620036]
for example -0.0596478 in second element of last list.Is there anyway to check the CI or stdev of the index (not clear form documentation) to see if 0 is included in the index CI - This would be useful to determine if should assume the value of the index to be zero or if I have to rerun the sensitivity analysis with a higher N.
Thanks,
A
Currently it's all written for the expanded form, but it would be nice to auto-extract from the Dual numbers as well.
If t is not provided:
julia> adjoint_sensitivities(sol,Vern9(),dg)
ERROR: type ODEAdjointSensitivityFunction has no field g_gradient_configIt would be useful for me to understand how the continuous case is different from the discrete one; at present it doesn't appear to be tested and has bitrotted.
In the documentation there are a few examples where the parameter sensitivity of an ODE by dual numbers is compared to other methods, such as local sensitivity analysis or using Calculus.jl. But this sensitivity analysis is only done at the end point of the simulation. When calculating sensitivities at other time points these methods don't line up for me.
Comparing sensitivity of the first state to the first parameter by forward differentiation and local sensitivity analysis:
using ForwardDiff
using ForwardDiff: Partials, Dual
using DifferentialEquations
using Plots
using Calculus
#forward differentiation
function func(du,u,p,t)
du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
du[2] = -3 * u[2] + u[1]*u[2]
end
p1 = 1.5
p2 = 1.0
u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
p = [p1, p2]
p1dual = Dual{Float64}(p1, (1.0, 0.0))
p2dual = Dual{Float64}(p2, (0.0, 1.0))
pdual = [p1dual, p2dual]
prob_dual = ODEProblem(func,eltype(pdual).(u0),eltype(pdual).(tspan),pdual)
sol_dual = solve(prob_dual,Tsit5(),reltol = 1e-15)
timepoints = [i.value for i in sol_dual.t]
sensitivity_forward_diff = [i[1].partials.values[1] for i in sol_dual.u]
plot(timepoints,sensitivity_forward_diff)
# sensitivity ODE
f_ode_sen = @ode_def_nohes test_sensitivity begin
du1 = p1 * u1 - p2 * u1*u2
du2 = -3 * u2 + u1*u2
end p1 p2
prob_ode_sen = ODELocalSensitivityProblem(f_ode_sen,u0,tspan,p)
sol_ode_state_and_sen = solve(prob_ode_sen,Tsit5(),reltol = 1e-9)
timepoints2 = [i for i in sol_ode_state_and_sen.t]
sensitivity_ode_sol = [i[3] for i in sol_ode_state_and_sen.u]
state_ode_sol = [i[3] for i in sol_ode_state_and_sen.u]
plot!(timepoints2,sensitivity_ode_sol)
I'm not sure I did everything right, but the way the solution by dual numbers snaps to the correct position only at the end seems odd to me.
Sensitivity analysis (Morris) does not allows for parameter fixing -
could also be an issue in other global methods, but I have not tried yet
param_range=[[0,5],[2,4],[3,3]]
param_step=[10,10,1] # didn'r work with [10,10,0] either or [10,10,2] either
BoundsError: attempt to access 1-element StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} at index [2]
Thanks
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.