search

SciML / ExponentialUtilities.jl

Fast and differentiable implementations of matrix exponentials, Krylov exponential matrix-vector multiplications ("expmv"), KIOPS, ExpoKit functions, and more. All your exponential needs in SciML form.
https://docs.sciml.ai/ExponentialUtilities/stable/
Other
95 stars 30 forks source link

Error using expv with Turing.jl #81

Closed aelmokadem closed 2 years ago

aelmokadem commented 2 years ago

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
ChrisRackauckas commented 2 years ago

gebal! is a LAPACK function which is not compatible with Dual numbers which are used by Turing's automatic differentiation. @sethaxen do you have a fallback?

sethaxen commented 2 years ago

gebal! is a LAPACK function which is not compatible with Dual numbers which are used by Turing's automatic differentiation. @sethaxen do you have a fallback?

No, ExponentialAction's algorithm doesn't use matrix balancing. One would need to implement matrix balancing in ExponentialUtilities to replace its LAPACK.gebal! calls.

Using ExponentialAction.expv works but very slow so I wanted to compare the performance to ExponentialUtilities!

@aelmokadem I'm still happy to look into this if you can provide a case that is much slower with ExponentialAction. A benchmark not using automatic differentiation is fine.

Ideally, I would be able to even use ExponentialUtilities.expv_timestep to get the full matrix without looping through timesteps!

Are your time steps evenly spaced? If so, there's a variant of the algorithm used by ExponentialAction that I haven't implemented yet that computes all expv(t, A, v) more efficiently. If that might be of use to you, feel free to open an issue there.

aelmokadem commented 2 years ago

@ChrisRackauckas and @sethaxen .

@sethaxen I am still trying to put together that minimal example that would reproduce the difference in performance. Unfortunately, I need the timesteps to match the observed data timepoints and those are not evenly spaced. That expv implementation still sounds very useful.

sethaxen commented 2 years ago

@ChrisRackauckas:

gebal! is a LAPACK function which is not compatible with Dual numbers which are used by Turing's automatic differentiation. @sethaxen do you have a fallback?

You might try GenericSchur.balance!, which is a pure Julia implementation of gebal!

ChrisRackauckas commented 2 years ago

It needs a little work first.

https://github.com/RalphAS/GenericSchur.jl/issues/6 https://github.com/RalphAS/GenericSchur.jl/pull/7

RalphAS commented 2 years ago

@ChrisRackauckas I've registered a new version of GenericSchur including your PR.

ChrisRackauckas commented 2 years ago

Hmm, it still doesn't seem to swap right in:

Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.002503127605795 0.0; 0.005025067838738853 1.007528195444534] ≈ [1.007528195444534 0.005025067838738853; 0.0 1.002503127605795]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.0451565331287087 0.0; 0.09652248296047598 1.1416790160891848] ≈ [1.1416790160891848 0.09652248296047598; 0.0 1.0451565331287087]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.2214027581601699 0.0; 0.6007160422303388 1.8221188003905084] ≈ [1.8221188003905084 0.6007160422303389; 0.0 1.2214027581601699]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.6625180212410025 0.0; 2.9326255480656873 4.595143569306688] ≈ [4.595143569306687 2.9326255480656864; 0.0 1.6625180212410025]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [3.4903429574618414 0.0; 39.03073904260098 42.52108200006281] ≈ [42.52108200006297 39.030739042601105; 0.0 3.4903429574618405]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [14.879731724872842 0.0; 3279.5883435589662 3294.4680752838403] ≈ [3294.4680752838403 3279.588343558967; 0.0 14.879731724872828]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [221.40641620418734 0.0; 1.0853298492648203e7 1.0853519899064412e7] ≈ [1.0853519899064412e7 1.0853298492648207e7; 0.0 221.4064162041869]     
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [49020.80113638183 0.0; 1.1779889415036631e14 1.1779889419938717e14] ≈ [1.1779889419938717e14 1.1779889415036636e14; 0.0 49020.80113638164]  
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [2.403038944052694e9 0.0; 1.3876579474598405e28 1.3876579474598412e28] ≈ [1.3876579474598412e28 1.387657947459841e28; 0.0 2.4030389440526757e9]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [5.774596166633886e18 0.0; 1.9255945791484585e56 1.9255945791484593e56] ≈ [1.9255945791484593e56 1.925594579148459e56; 0.0 5.774596166633799e18]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [3.3345960887702767e37 0.0; 3.7079144832459305e112 3.707914483245932e112] ≈ [3.707914483245932e112 3.7079144832459317e112; 0.0 3.334596088770176e37]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.267963432611071e50 0.0; 2.0385444553099337e150 2.0385444553099333e150] ≈ [2.0385444553096833e150 2.0385444553096833e150; 0.0 1.267963432611071e50]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Exp generated: Test Failed at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56
  Expression: A ≈ expA
   Evaluated: [1.3609151355757742e100 0.0; 2.5205373219390364e300 2.520537321939036e300] ≈ [2.520537321939602e300 2.5205373219396014e300; 0.0 1.36091513557564e100]
Stacktrace:
 [1] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:445 [inlined]
 [2] macro expansion
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:56 [inlined]
 [3] macro expansion
   @ C:\Users\accou\.julia\juliaup\julia-1.7.2+0~x64\share\julia\stdlib\v1.7\Test\src\Test.jl:1283 [inlined]
 [4] top-level scope
   @ C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:39
Test Summary: | Fail  Total
Exp generated |   13     13
ERROR: LoadError: Some tests did not pass: 0 passed, 13 failed, 0 errored, 0 broken.
in expression starting at C:\Users\accou\.julia\dev\ExponentialUtilities\test\runtests.jl:37
ERROR: Package ExponentialUtilities errored during testing