I am trying to estimate parameters for an ODE system using Turing, ForwardDiff, and DifferentialEquations. However, I consistently get the following error:
julia> TypeError: in setfield!, expected Float64, got ForwardDiff.Dual{Nothing,Float64,12}
Stacktrace:
[1] setproperty! at .\Base.jl:34 [inlined]
[2] set_W_γdt! at C:\Users\d.julia\packages\OrdinaryDiffEq\VPJBD\src\nlsolve\utils.jl:34 [inlined]
[3] set_W_γdt! at C:\Users\d.julia\packages\OrdinaryDiffEq\VPJBD\src\nlsolve\utils.jl:32 [inlined] …
I tried logging the full stack trace, which can be viewed here:
https://cl1p.net/dfasdkmjnf384h
I found some related Q&As here, including [1], [2], [3], [4], [5], and [6]. Following the advice there, I tried to make my code play well with Dual numbers, but haven’t found a solution yet. I’ve included (what I think are) the relevant parts of my code below. Any help would be greatly appreciated. Thanks!
using Glob, CSV, Plots, StatsPlots
using DifferentialEquations, OrdinaryDiffEq, DiffEqSensitivity
using DiffEqBase: get_tmp, dualcache
using Turing, MCMCChains
using ForwardDiff
using Random, Distributions
function tm(v::Float64, θ::AbstractVector{T}) where T
#elements of Q are functions of theta and `v`, e.g. a1 = theta[1]*(v/theta[2])
a0, sa, b0, sb, c0, sc, d0, sd, g1, g2, h1, h2 = θ
a1 = a0*exp(-v/sa)
b1 = b0*exp(v/sb)
a2 = c0*exp(-v/sc)
b2 = d0*exp(v/sd)
Q = T[a1 -(b1+g1+a1) -a1 h1-a1;
g2 -g2 -(a2+h2+g2) b2-g2;
0.0 g1 a2 -(b2 + h1)]
return Q
end
# get starting values
function get_inits(Qh)
return Qh[:,2:end] \ (-1 .* Qh[:,1])
end
# ODE
function rhs!(du, u, Q, t)
du .= Q[:,2:end]*u + Q[:,1]
end
# This gives a different error from the other things I've tried:
# ERROR: LoadError: MethodError: no method matching get_tmp(::Float64, ::Array{Float64,1})
# function rhs!(du, u, (Q, tmp), t)
# tmp = DiffEqBase.get_tmp(tmp, u)
# tmp .= Q[:,2:end]*u + Q[:,1]
# @. du = u + tmp
# nothing
# end
times = collect(1.0:1000.0)
Nt = length(times)
voltages = collect(-150.0:50.0:-50.0)
holds = [-20.0, -180.0]
data_to_fit = rand(Uniform(0,1), (length(voltages), length(holds)*Nt))
function do_sim(θ::AbstractVector{T};
ydata = data_to_fit, times=times,
Nt=Nt, voltages=voltages, holds=holds) where {T<:Real}
"""
θ = parameters; used to construct a matrix when solving ODEs
data = Array{Float64,2}
time, voltages, holds = Array{Float64,1}
Nt = Int64
"""
# array to hold simulation values
ysim = zeros(T, size(ydata))
ts = (times[1], times[end])
# find initial values
Qh1 = tm(holds[1], θ)
u1 = convert.(T, Qh1[:,2:end] \ (-1 .* Qh1[:,1]))
# tried this, too (doesn't work)
# u1 = get_inits(Qh1)
Qh2 = tm(holds[2], θ)
u2 = convert.(T, Qh2[:,2:end] \ (-1 .* Qh2[:,1]))
for (j, vt) in enumerate(voltages)
Qv = tm(vt, θ)
prob1 = ODEProblem(rhs!, u1, ts, Qv)
prob2 = ODEProblem(rhs!, u2, ts, Qv)
# tried this with `theta` and `Qv` to make u1/u2 the same type
# prob1 = ODEProblem(rhs!, eltype(θ).(u1), ts, Qv)
# solve ODEs
sol1 = solve(prob1, TRBDF2(), saveat=times)
sol2 = solve(prob2, TRBDF2(), saveat=times)
# tried TRBDF2(autodiff=false)
# add to array
ysim[j,1:Nt] = sum(sol1[2:end,:], dims=1)
ysim[j,1+Nt:end] = sum(sol2[2:end,:], dims=1)
end
return ysim
end
# some parameters
p0 = ([-10.783723199314082, 2.768277944084355, -0.9532448651472752,
3.806992452948674, -12.864645880101888, 2.262134729959288, -3.8200088837688195, 3.1458924991336534, -4.992066945500882, -7.758380598719609, -9.485975086743023, -4.298650310672511])
Np = length(p0)
Nv = length(voltages)
Turing.setadbackend(:forwarddiff)
@model function ChTurMod(y, ::Type{T} = Float64) where T
"""
y = data array
"""
# What I have right now, doesn't work.
θ = Vector{T}(undef, Np)
for i = 1:Np
θ[i] ~ Normal(0, 1)
end
# Tried this as well, didn't work.
# θ ~ MvNormal(zeros(T, 12),1)
ysim = do_sim( exp.(θ .+ p0) )
for j = 1:Nv
y[j,:] ~ MvNormal(ysim[j,:], 0.02)
end
end
# run a simulation
# ysim = do_sim(exp.(p0))
# Output of the above should look like:
# 3×2000 Array{Float64,2}:
# 0.030866 0.0319676 0.0342911 … 0.988392 0.988392
# 0.030866 0.0309165 0.0310629 0.945324 0.945307
# 0.030866 0.0308674 0.0308714 0.133803 0.133578
# sampling
to_fit = ChTurMod(data_to_fit)
chain = sample(to_fit, Turing.NUTS(0.65), 1000)
# plot simulation with random data and default parameters
#=
ysim = do_sim(exp.(p0))
plt = plot(size=(400,280))
plt2 = plot(size=(400,280))
for j = 1:size(ysim)[1]
plot!(plt, times, ysim[j,1:Nt], lw=2, label="")
plot!(plt2, times, ysim[j, Nt+1:end], lw=2, label="")
# plot random data
plot!(plt, times, data_to_fit[j,1:Nt],
lalpha=0.2, label="")
plot!(plt2, times, data_to_fit[j, Nt+1:end],
lalpha=0.2, label="")
end
display(plt)
display(plt2)
=#```