I’m using Julia since today so I might miss something fundamental. I recently discovered the UDE for SciML paper and the super nice SciML library. I try to rebuild something similar to the Lotka-Voltera example from the paper, but with a Stiff System (Oregonator), but fail to do so.

When running

using DifferentialEquations, Plots, DiffEqFlux, Flux, Optim, OrdinaryDiffEq, DiffEqSensitivity

#build the true solution

tspan = (0.0, 360.0)

tsteps = 0.0:1.0:360.0

u0 = [1.0,2.0,3.0]function orego(du,u,p,t)

s,q,w = p

y1,y2,y3 = u

du[1] = s*(y2+y1*(1-q*y1-y2))

du[2] = (y3-(1+y1)y2)/s(y1-y3)

du[3] = w

endp_true = [77.27,8.375e-6,0.161]

prob_true = ODEProblem(orego,u0,tspan,p_true)sol = solve(prob_true,Rodas5(),saveat = tsteps)

true_sol = Array(sol)## plot(sol)

#build the UDE

ann = FastChain(FastDense(3,16,tanh), FastDense(16,16,tanh), FastDense(16,1))

p1 = initial_params(ann)

p2 = [0.5,1, 0.]

p_init = [p1;p2]

theta = Flux.params(p_init)function orego_ude!(du,u,p,t)

s,q,w = p[end-2:end]

y1,y2,y3 = u

du[1] = s*(y2+y1*(1-qy1-y2))(y1-y3)

du[2] = ann(u,p[1:length(p1)])[1]

du[3] = w

endprob_mod = ODEProblem(orego_ude!,u0,tspan,p_init)

sol_mod_test = Array(solve(prob_mod,Rodas4P(), u0=u0, p=p_init,saveat = tsteps, abstol=1e-12, reltol=1e-12,sensealg = InterpolatingAdjoint(checkpointing=true)))#ML part

function predict_mod(theta)

return Array(solve(prob_mod,Rodas4P(), u0=u0, p=p_init,saveat = tsteps, abstol=1e-12, reltol=1e-12,sensealg = InterpolatingAdjoint(checkpointing=true)))

endloss(theta) = sum(abs2, predict_mod(theta).-true_sol)

l = loss(p_init)cb = function (theta,l)

println(l)

return false

endresult_mod = DiffEqFlux.sciml_train(loss, p_init, ADAM(0.01), cb = cb,maxiters = 200)

Julia throws a huge stack on the last line, beginning with

LoadError: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,12}

in expression starting at /somepath/oregonator_ude.jl:57

Which confuses me quite a bit in many ways, the first being me not understanding Julia’s approach to types: I was assuming that one does not have to take care of explicitly defining types.

Also I do not understand why sol_mod_test gets assigned without problem, the same line later however fails as long as it contains a stiff solver and or the sensealg argument (With Rodas5() it works but reaches max_iter, which I kinda expect due to the nature of the System). Due to the behaviour I suspect it has to do something with the ODE solver evaluating Jacobians (or not) during the sensitivity calculations but I don’t see what I have to change, even after consulting the documentation.