Error using expv with Turing.jl #81
Description
Hi,
I am running a Bayesian inference on an ODE model using Turing.jl. To try and speed up things, I represented the model as a linear system and expressed the coefficients in the form of a matrix to solve the system using matrix exponential (ExponentialUtilities.expv()). When I run the inference, I keep getting the error below (this is only the head of the error but I can try to share the rest if necessary). The error message is associated to the following line but I am guessing all other calls of ExponetialUtilities.expv would produce it as well:
u0_eoi = inv(K)*ExponentialUtilities.expv(durs[i], K, tmp_rate) - inv(K)*tmp_rate
Using ExponentialAction.expv works but very slow so I wanted to compare the performance to ExponentialUtilities! Ideally, I would be able to even use ExponentialUtilities.expv_timestep to get the full matrix without looping through timesteps!
Any ideas on how to solve this! Also, I would greatly appreciate it if you could share any alternative more efficient method to carry out the analysis with a linear solver!
Thank you.
Ahmed.
PS: the full model structure is also below. It is a PBPK model where I am modeling a drug's pharmacokinetics in patients. The drug is introduced via IV infusion so the model is set to calculate the state variables using the relation u(t)=A^-1 * exp(At) * rate - A^-1 * rate. Once infusion duration is over I switch to u(t)=exp(At) * u0 The CLint parameter includes random effects relevant to interindividual variability. The function PBPK() would take in the parameter set and generate the full model coefficients matrix. Within this function there is a line:
K = zeros(Base.promote_eltype(p[1]), ncmt, ncmt);
This line promotes the input parameter type when creating the K matrix.
I haven't attached the full function description but I can if necessary!
The error message:
ERROR: MethodError: no method matching gebal!(::Char, ::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Essential.var"#f#4"{DynamicPPL.TypedVarInfo{NamedTuple{(:σ, :ĈLint, :KbBR, :KbMU, :KbAD, :KbBO, :KbRB, :ω, :ηᵢ), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:σ, Setfield.IdentityLens}, Int64}, Vector{Truncated{Cauchy{Float64}, Continuous, Float64}}, Vector{AbstractPPL.VarName{:σ, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:ĈLint, Setfield.IdentityLens}, Int64}, Vector{LogNormal{Float64}}, Vector{AbstractPPL.VarName{:ĈLint, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}},
The model structure:
@model function fitPBPK(data, nSubject, rates, times, durs, wts, VVBs, BP, starts, ends, ncmt) # data should be a Vector
# priors
σ ~ truncated(Cauchy(0, 0.5), 0.0, Inf) # residual error
ĈLint ~ LogNormal(7.1,0.25)
KbBR ~ LogNormal(1.1,0.25)
KbMU ~ LogNormal(0.3,0.25)
KbAD ~ LogNormal(2.0,0.25)
KbBO ~ LogNormal(0.03, 0.25)
KbRB ~ LogNormal(0.3, 0.25)
ω ~ truncated(Cauchy(0, 0.5), 0.0, 1.0)
ηᵢ ~ filldist(Normal(0.0, 1.0), nSubject)
# individual params
CLintᵢ = ĈLint .* exp.(ω .* ηᵢ) # non-centered parameterization
# solve
tmp_rate = zeros(ncmt)
predicted = []
ps_tmp = [CLintᵢ[1], KbBR, KbMU, KbAD, KbBO, KbRB, wts[1]]
K_tmp = PBPK(ps_tmp)
u0_eoi_tmp = inv(K_tmp)*ExponentialAction.expv(durs[1], K_tmp, tmp_rate) - inv(K_tmp)*tmp_rate
ut = zeros(Base.promote_eltype(times[1], K_tmp, u0_eoi_tmp), length(u0_eoi_tmp), length(data))
for i in 1:nSubject
tmp_rate[15] = rates[i]
ps = [CLintᵢ[i], KbBR, KbMU, KbAD, KbBO, KbRB, wts[i]]
K = PBPK(ps)
#u0_eoi = inv(K)*ExponentialAction.expv(durs[i], K, tmp_rate) - inv(K)*tmp_rate
u0_eoi = inv(K)*ExponentialUtilities.expv(durs[i], K, tmp_rate) - inv(K)*tmp_rate
n = 1
for j in starts[i]:ends[i]
tmp_time = times[i][n]
if tmp_time <= durs[i]
#ut[:,j] = inv(K)*ExponentialAction.expv(tmp_time, K, tmp_rate) - inv(K)*tmp_rate
ut[:,j] = inv(K)*ExponentialUtilities.expv(tmp_time,K,tmp_rate) - inv(K)*tmp_rate
else
#ut[:,j] = ExponentialAction.expv(tmp_time, K, u0_eoi)
ut[:,j] = ExponentialUtilities.expv(tmp_time,K,u0_eoi)
end
n = n + 1
end
tmp_sol = ut[15,starts[i]:ends[i]] ./ (VVBs[i]*BP/1000.0)
append!(predicted, tmp_sol)
end
# likelihood
for i = 1:length(predicted)
data[i] ~ LogNormal(log.(max(predicted[i], 1e-12)), σ)
end
end