Problem running simulations using a function-based interface model definition

Hi! I am having trouble using the “function-based interface” to define a Pumas model and then run simulations. I hope you could help me figure it out.
I am using the example code from the function-based interface documentation (https://docs.pumas.ai/basics/models/#The-PumasModel-Function-Based-Interface) and defining the parameter values, observation times, dose regimes and subjects so then I can run simobs in the same way as I would have done when defining the model through the macro @model. At the end I get an error that I can’t understand what I should change…
This is the code:

    paramset = ParamSet((θ = VectorDomain(4, lower=zeros(4)), # parameters
                  Ω = PSDDomain(2),
                  Σ = RealDomain(lower=0.0)))

    function random(p)
      ParamSet((η=MvNormal(p.Ω),))
    end

    function pre(param,randeffs,subject)
      function pre_t(t)
        cov_t = subject.covariates(t)
        CL = param.θ[2] * ((cov_t.wt/70)^0.75) * (param.θ[4]^cov_t.sex) * exp(randeffs.η[1])
        return (Σ=param.Σ, Ka=param.θ[1], CL=CL, V=param.θ[3] * exp(randeffs.η[1]))
      end
    end

    function init(col,t0)
       [0.0,0.0]
    end

    function onecompartment_f(du,u,p,t)
        du[1] = -p.Ka*u[1]
        du[2] =  p.Ka*u[1] - (p.CL/p.V)*u[2]
    end

    prob = ODEProblem(onecompartment_f,nothing,nothing,nothing)

    function derived(col, sol, obstimes, subject)
        central = sol(obstimes)
        conc = @. central / col.V
        dv = @. Normal(conc, conc*col.Σ)
        (dv = dv, conc = conc, )
    end

    function observed(col,sol,obstimes,samples,subject)
        (obs_cmax = maximum(samples.dv),
         T_max = maximum(obstimes),
         dv = samples.dv)
    end

    param = (θ = Diagonal([0.05,0.05,0.05,0.05]),
                       Ω = Diagonal([0.05]),
                       Σ = 0.02)

    regimen = DosageRegimen(150, rate = 10, cmt = 1)

    sub = Subject(id=1, events=regimen, covariates=(sex=0,wt=70))

    obstimes = [0.25, 0.5, 0.75, 1, 2, 4, 8,
                    12, 16, 20, 24, 36, 48, 60, 71.9] # single dose observation times

    p_model = PumasModel(paramset,
                          random,
                          pre,
                          init,
                          prob,
                          derived,
                          observed)

    Random.seed!(123)

    sims = simobs(p_model, sub, param, obstimes= obstimes)

This is the error:

    <strong>MethodError: no method matching derived(::var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}, ::OrdinaryDiffEq.ODECompositeSolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}},ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:callback, :saveat, :tstops, :d_discontinuities, :save_first),Tuple{CallbackSet{Tuple{},Tuple{DiscreteCallback{Pumas.var"#783#786"{Array{Float64,1},Base.RefValue{Bool},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Bool},Base.RefValue{Int64}},Pumas.var"#affect!#788"{DataType,Array{Float64,1},Bool,Symbol,Float64,Float64,Int64,Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Base.RefValue{Bool},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Float64},Array{Float64,1},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Float64},Base.RefValue{Int64},Base.RefValue{Bool},Base.RefValue{Int64},Base.RefValue{Int64},Base.RefValue{Int64},Base.RefValue{Int64},Array{Float64,1},Array{Float64,1},Base.RefValue{Float64},Base.RefValue{Bool}},typeof(Pumas.subject_cb_initialize!)}}},Array{Float64,1},Array{Float64,1},Float64,Bool}}},DiffEqBase.StandardODEProblem},CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType},Rational{Int64},Int64}},OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23Tableau{Float64},DiffEqBase.TimeGradientWrapper{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}},DiffEqBase.UJacobianWrapper{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}},Float64},Float64,2},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.UJacobianWrapper{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}},Float64},Float64,2},1},Array{Float64,1},Array{Array{Tuple{Bool,Bool},1},1},UnitRange{Int64},Nothing},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.TimeGradientWrapper{ODEFunction{true,Pumas.DiffEqWrapper{typeof(onecompartment_f),Array{Float64,1}},UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},var"#pre_t#26"{NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}},NamedTuple{(:η,),Tuple{Array{Float64,1}}},Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}}},Float64},Float64,1},1}}},AutoSwitch{Tsit5,Rosenbrock23{0,true,DefaultLinSolve,DataType},Rational{Int64},Int64}}},DiffEqBase.DEStats}, ::Array{Float64,1}, ::Subject{Nothing,Pumas.ConstantCovar{NamedTuple{(:sex, :wt),Tuple{Int64,Int64}}},Array{Pumas.Event,1},Nothing}, ::NamedTuple{(:θ, :Ω, :Σ),Tuple{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}},Float64}}, ::NamedTuple{(:η,),Tuple{Array{Float64,1}}})</strong>

    <strong>Closest candidates are:
    derived(::Any, ::Any, ::Any, ::Any) at C:\Users\workspace\test.jl:46</strong>

    in top-level scope at [universal_ODE_PKPD.jl:79](#)

    in at [Pumas\RVlQK\src\models\model_api.jl:423](#)

    in at [Pumas\RVlQK\src\models\model_api.jl:423](#)

    in #simobs#177 at [Pumas\RVlQK\src\models\model_api.jl:439](#)

Thank you for posting this issue. There were indeed some parts there that had not been updated to Pumas v1.0 syntax. I will update the documentation, but before that is going to be visible, here is a version that works. I suspect this might have come from you trying to debug the issue, but you seem to have made the θ in your input param a Diagonal which I removed, and you had also removed the idxs=2 keyword to the solution in the derived block. The point of this keyword is to grab the second compartment’s values and store it in a variable central.

The other changes I made to make this work was the handling of the covariates in the derived block and the fact that the positional arguments to derived has changed from previous versions of Pumas. I will update the documentation. Thanks.

using Pumas, Random
paramset = ParamSet((θ = VectorDomain(4, lower=zeros(4)), # parameters
                  Ω = PSDDomain(2),
                  Σ = RealDomain(lower=0.0)))

function random(p)
  ParamSet((η=MvNormal(p.Ω),))
end

function pre(param,randeffs,subject)
  function pre_t(t)
    cov_t = subject.covariates(t)
    CL = param.θ[2] * ((cov_t.wt/70)^0.75) * (param.θ[4]^cov_t.sex) * exp(randeffs.η[1])
    return (Ka=param.θ[1], CL=CL, Vc=param.θ[3] * exp(randeffs.η[1]))
  end
end

function init(col,t0)
   [0.0,0.0]
end

function onecompartment_f(du,u,p,t)
    du[1] = -p.Ka*u[1]
    du[2] =  p.Ka*u[1] - (p.CL/p.Vc)*u[2]
end

prob = ODEProblem(onecompartment_f,nothing,nothing,nothing)

function derived(col, sol, obstimes, subject, param, random)
    Vc = map(t->col(t).Vc, obstimes)
    central = sol(obstimes; idxs=2)
    conc = @. central / Vc
    dv = @. Normal(conc, conc*param.Σ)
    (dv = dv, conc = conc, )
end

function observed(col,sol,obstimes,samples,subject)
    (obs_cmax = maximum(samples.dv),
     T_max = maximum(obstimes),
     dv = samples.dv)
end

param = (θ = [0.05,0.05,0.05,0.05],
                   Ω = Diagonal([0.05]),
                   Σ = 0.02)

regimen = DosageRegimen(150, rate = 10, cmt = 1)

sub = Subject(id=1, events=regimen, covariates=(sex=0,wt=70))

obstimes = [0.25, 0.5, 0.75, 1, 2, 4, 8,
                12, 16, 20, 24, 36, 48, 60, 71.9] # single dose observation times

p_model = PumasModel(paramset,
                      random,
                      pre,
                      init,
                      prob,
                      derived,
                      observed)



sims = simobs(p_model, sub, param, obstimes= obstimes)
1 Like

Thanks for the working code and the explanation! I am looking forward to the documentation update.

I updated the code. I just removed the sigma parameter from pre and just took it from param directly. This will be the code that is in the docs as soon as the website updates to reflect my changes. Thanks again for reporting this issue!