Error with trunc!

[aelmokadem]

Hi,

I am running a Bayesian inference on an ODE model using Turing.jl. To 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 (ExponentialAction.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).

Any ideas on how to solve this!

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(A*t) * rate - A^-1 * rate. Once infusion duration is over I switch to u(t)=exp(A*t) * 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. I haven't attached the full function description but I can if necessary!

ERROR: InexactError: trunc(Int64, 1.0868108130524078e19)
Stacktrace:

trunc at [./float.jl]()

ceil at [./float.jl]()

parameters(t::Float64, A::Matrix{Float64}, n0::Int64, m_max::Int64, p_max::Int64, tol::Float64) at [/Users/ahmede/.julia/packages/ExponentialAction/GSSyv/src/parameters.jl]()

expv(t::Float64, A::Matrix{Float64}, B::Vector{Float64}; shift::Bool, tol::Float64) at [/Users/ahmede/.julia/packages/ExponentialAction/GSSyv/src/expv.jl]()

expv(t::Float64, A::Matrix{Float64}, B::Vector{Float64}) at [/Users/ahmede/.julia/packages/ExponentialAction/GSSyv/src/expv.jl]()

Full 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)
    #ut = zeros(ncmt,length(data))
    predicted = []

    ps_tmp = [CLintᵢ[1], KbBR, KbMU, KbAD, KbBO, KbRB, wts[1]]
    K_tmp = PBPK(ps_tmp)
    #tmp_rate[15] = rates[1]
    #u0_eoi_tmp = inv(K_tmp)*exp(K_tmp*durs[1])*tmp_rate - inv(K_tmp)*tmp_rate
    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)*exp(K*durs[i])*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.exp_generic(K*tmp_time)*tmp_rate - inv(K)*tmp_rate
            else
                ut[:,j] = ExponentialAction.expv(tmp_time, K, u0_eoi)
                #ut[:,j] = ExponentialUtilities.exp_generic(K*tmp_time)*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.(predicted[i]), σ)
    end
end

[sethaxen]

The line where the error is thrown is https://github.com/sethaxen/ExponentialAction.jl/blob/6ee4452473f53423f5a8c7e34fae1e0d67701d2e/src/parameters.jl#L24. The truncation error is raised because the ratio is >1e19, which is too large to be encoded by the integers. This is likely because t or the operator norm of A is very large. I'll need to look into if there's a way for the algorithm to be more robust in this case (e.g. maybe not using integers).

But if you could provide a combination of inputs (durs[1], K_tmp, tmp_rate) that raise this error, that would be very helpful for debugging. You could e.g. call something like _expv(args...) = (print(args); expv(args...)) and provide the output if your arrays are small; otherwise you could write them to a file.

[aelmokadem]

Thanks a lot @sethaxen for responding. I really appreciate your help with this.

Found a combination that would give the error. Seems that the error comes when tmp_time exceeds some value (maybe > 5)

The problematic line:

ut[:,j] = ExponentialAction.expv(tmp_time, K, u0_eoi)

variables:

tmp_time = 12.2

K = [-284.23904664215075 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 87.52510874621326 0.0; 0.0 -13.520975517897782 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.002008699697061 0.0 0.0; 0.0 0.0 -11.496159804401913 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 21.00602609909118 0.0 0.0; 0.0 0.0 0.0 -0.9334895320387491 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 29.758536973712506 0.0 0.0; 0.0 0.0 0.0 0.0 -0.20544544006353113 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8.752510874621327 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -12.011005732426296 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8.752510874621327 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 -41.66847554687845 0.0 0.0 0.0 0.0 0.0 0.0 5.251506524772795 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -25.5745116248461 0.0 0.0 0.0 0.0 0.0 1.7505021749242653 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 41.66847554687845 25.5745116248461 -12.17145662311623 17.13125180611701 14.407974904677381 0.0 0.0 11.37826413700772 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -17.13125180611701 0.0 0.0 0.0 1.7505021749242653 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.407974904677381 0.0 0.0 24.507030448939712 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.187611178419207 0.0 8.752510874621327 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -57.45351404671386 33.25954132356104 0.0 0.0; 284.23904664215075 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -175.05021749242653 0.0 0.0; 0.0 13.520975517897782 11.496159804401913 0.9334895320387491 0.20544544006353113 12.011005732426296 0.0 0.0 8.393517993108697 0.0 0.0 4.187611178419207 57.45351404671386 0.0 -87.52510874621326 6.276140811583681; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 13.128766311931969 0.0 -6.276140811583681]

u0_eoi = [0.6494829319865283, 0.4763956456413778, 1.6303307144127004, 6.310148634770343, 2.0414672101868234, 0.656029914050817, 0.12736772023833332, 0.06749413530797493, 2.9703026680965188, 0.09719953864088937, 1.5807160635494348, 1.269714633246279, 0.5909292617482294, 1.0464468508457871, 2.1191493487017077, 1.54796784340739]

The full error message when I isolated that for loop:

julia> ut[:,j] = ExponentialAction.expv(tmp_time, K, tmp_rate)
ERROR: InexactError: trunc(Int64, 1.0549418223667896e19)
Stacktrace:
 [1] trunc
   @ ./float.jl:805 [inlined]
 [2] ceil
   @ ./float.jl:368 [inlined]
 [3] parameters(t::Int64, A::Matrix{Float64}, n0::Int64, m_max::Int64, p_max::Int64, tol::Float64)
   @ ExponentialAction ~/.julia/packages/ExponentialAction/GSSyv/src/parameters.jl:46
 [4] expv(t::Int64, A::Matrix{Float64}, B::Vector{Float64}; shift::Bool, tol::Float64)
   @ ExponentialAction ~/.julia/packages/ExponentialAction/GSSyv/src/expv.jl:35
 [5] expv(t::Int64, A::Matrix{Float64}, B::Vector{Float64})
   @ ExponentialAction ~/.julia/packages/ExponentialAction/GSSyv/src/expv.jl:24
 [6] top-level scope
   @ REPL[123]:1

[sethaxen]

Thanks! I was able to confirm the error, and it looks like an easy fix. I'll open a PR.