ERROR: PosDefException

githubmanojkumar · October 14, 2023, 4:41am

Code: obs1= simobs(Meropenam,pop,params1,obstimes = [0,0.5,2,8])

Error:PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] checkpositivedefinite
    @ /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/factorization.jl:18 [inlined]
  [2] cholesky!(A::LinearAlgebra.Hermitian{Float64, Matrix{Float64}}, ::LinearAlgebra.NoPivot; check::Bool)
    @ LinearAlgebra /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:268
  [3] cholesky!
    @ /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:266 [inlined]
  [4] cholesky!(A::Matrix{Float64}, ::LinearAlgebra.NoPivot; check::Bool)
    @ LinearAlgebra /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:300
  [5] cholesky! (repeats 2 times)
    @ /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:294 [inlined]
  [6] #cholesky#162
    @ /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:400 [inlined]
  [7] cholesky (repeats 2 times)
    @ /opt/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/cholesky.jl:400 [inlined]
  [8] PDMats.PDMat(mat::Matrix{Float64})
    @ PDMats /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/PDMats/CbBv1/src/pdmat.jl:19
  [9] MvNormal
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/uREy8/src/multivariate/mvnormal.jl:201 [inlined]
 [10] MvNormal
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Distributions/uREy8/src/multivariate/mvnormal.jl:218 [inlined]
 [11] #7
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/dsl/model_macro.jl:368 [inlined]
 [12] sample_randeffs(rng::TaskLocalRNG, model::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, egfronCL = 1, Ω = 4, σ²_prop = 1), 4, (:Central, :Peripheral), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :egfronCL, :Ω, :σ²_prop), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#7#13", Pumas.TimeDispatcher{var"#8#14", var"#9#15"}, Nothing, var"#10#16", Pumas.LinearODE, var"#11#17", var"#12#18", ModelingToolkit.ODESystem}, param::NamedTuple{(:tvcl, :tvvc, :tvvp, :tvq, :egfronCL, :Ω, :σ²_prop), Tuple{Float64, Float64, Float64, Float64, Float64, Matrix{Float64}, Float64}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:627
 [13] simobs(::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, egfronCL = 1, Ω = 4, σ²_prop = 1), 4, (:Central, :Peripheral), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :egfronCL, :Ω, :σ²_prop), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#7#13", Pumas.TimeDispatcher{var"#8#14", var"#9#15"}, Nothing, var"#10#16", Pumas.LinearODE, var"#11#17", var"#12#18", ModelingToolkit.ODESystem}, ::Subject{NamedTuple{(:DV,), Tuple{Vector{Missing}}}, Pumas.ConstantCovar{NamedTuple{(:EGFR, :group), Tuple{Int64, String}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Nothing}, ::NamedTuple{(:tvcl, :tvvc, :tvvp, :tvq, :egfronCL, :Ω, :σ²_prop), Tuple{Float64, Float64, Float64, Float64, Float64, Diagonal{Float64, Vector{Float64}}, Float64}}, ::Nothing; obstimes::Vector{Float64}, ensemblealg::EnsembleSerial, diffeq_options::NamedTuple{(:alg, :callback), Tuple{OrdinaryDiffEq.CompositeAlgorithm{Tuple{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitch{OrdinaryDiffEq.Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, true, LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, Nothing}}, rng::TaskLocalRNG, simulate_error::Val{true}, isfor_derived::Val{false}, return_model::Val{true})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:1527
 [14] simobs
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:1500 [inlined]
 [15] #_simobs!#358
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:1795 [inlined]
 [16] #_simobs#378
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:2112 [inlined]
 [17] _simobs
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:2023 [inlined]
 [18] simobs(::PumasModel{(tvcl = 1, tvvc = 1, tvq = 1, tvvp = 1, egfronCL = 1, Ω = 4, σ²_prop = 1), 4, (:Central, :Peripheral), ParamSet{NamedTuple{(:tvcl, :tvvc, :tvq, :tvvp, :egfronCL, :Ω, :σ²_prop), Tuple{RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, var"#7#13", Pumas.TimeDispatcher{var"#8#14", var"#9#15"}, Nothing, var"#10#16", Pumas.LinearODE, var"#11#17", var"#12#18", ModelingToolkit.ODESystem}, ::Vector{Subject{NamedTuple{(:DV,), Tuple{Vector{Missing}}}, Pumas.ConstantCovar{NamedTuple{(:EGFR, :group), Tuple{Int64, String}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Nothing}}, ::NamedTuple{(:tvcl, :tvvc, :tvvp, :tvq, :egfronCL, :Ω, :σ²_prop), Tuple{Float64, Float64, Float64, Float64, Float64, Diagonal{Float64, Vector{Float64}}, Float64}}, ::Nothing; rng::TaskLocalRNG, kwargs::Base.Pairs{Symbol, Vector{Float64}, Tuple{Symbol}, NamedTuple{(:obstimes,), Tuple{Vector{Float64}}}})
    @ Pumas /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:1991
 [19] simobs
    @ /build/_work/PumasSystemImages/PumasSystemImages/julia_depot/packages/Pumas/VyE8h/src/models/model_api.jl:1966 [inlined]
 [20] top-level scope
    @ ~/data/code/Jyothi Meropenam/Merpenam.jl:81

kindly give your inputs in solving this error.
Thanks in advance

andreasnoack · October 14, 2023, 1:10pm

Would you be able to show us the model?

vijay · October 17, 2023, 10:17pm

Looks like this had been resolved in an offline thread. @githubmanojkumar can you please post how you fixed this error

githubmanojkumar · October 18, 2023, 5:16am

[quote=“githubmanojkumar, post:4, topic:948, full:true”]
Dear all,

I would like to inform you that I have addressed the error by making the necessary correction, which has been incorporated within the code. I sincerely appreciate your support throughout this learning process, and I look forward to your continued assistance in the future.

The model:

Meropenam= @model begin
@param begin
tvcl ∈ RealDomain(lower=0.001)
tvvc ∈ RealDomain(lower=0.001)
tvq ∈ RealDomain(lower=0.001)
tvvp ∈ RealDomain(lower=0.001)
egfronCL ∈ RealDomain(lower=0.001)
Ω ∈ PDiagDomain(4)

before i Have kept 4 in defining randomeffects but actually the parameters were 5

correction Ω ∈ PDiagDomain(5)**

σ²_prop     ∈ RealDomain(lower=0.001)

end

@random begin
η ~ MvNormal(Ω)
end

@covariates EGFR

@pre begin
Cl = tvcl *(1+(egfr-47)*egfronCL)*exp(η[1])
Vc = tvvc * exp(η[2])
Vp = tvvp * exp(η[3])
Q = tvq * exp(η[4])
end

@dynamics begin
Central’ = -(Cl/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
Peripheral’ = (Q/Vc)*Central - (Q/Vp)*Peripheral
end

@derived begin
CONC = @. Central/Vc
dv ~ @. Normal(CONC, sqrt(CONC^2*σ²_prop))

end
end

params1 = (tvcl =17.26, tvvc = 32.61 ,tvvp =14.5,tvq=15.3, egfronCL= 0.01013, Ω = Diagonal([0.2,0.2,0.2,0.,0.2]), σ²_prop = 0.01)

correction Ω = Diagonal([0.2,0.2,0.2,0.2,0.2])

Random.seed!(123)

egfr1 = rand(90:120,18)
egfr2 = rand(60:89,15)
egfr3 = rand(30:59,33)
egfr4 = rand(15:29,39)
egfr5 = rand(10:15,5)

EGFR = vcat(egfr1,egfr2,egfr3,egfr4,egfr5)

ev1 = DosageRegimen(500, cmt=1, rate = 0.5)
ev2 = DosageRegimen(1000, cmt = 1, rate = 3)

pop1 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr1[i],group=“500mg”),observations = (DV = nothing,),),1:6))
pop2 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr2[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop3 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr3[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop4 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr4[i],group=“500mg”),observations = (DV = nothing,),),1:9))
pop5 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr5[i],group=“500mg”),observations = (DV = nothing,),),1:2))
pop1a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr1[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop2a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr2[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop3a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr3[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop4a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr4[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop5a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr5[i],group=“1000mg”),observations = (DV = nothing,),),1:3))

pop = vcat(pop1,pop2,pop3,pop4, pop5, pop1a, pop2a, pop3a, pop4a, pop5a)

obs1= simobs(Meropenam,pop,params1,obstimes = [0,0.5,2,8])

With regards,
Dr. Manojkumar V, Pharm D.,
Research Scholar,
JSS College of Pharmacy, Ooty
Nilgiris - 643001

andreasnoack · October 18, 2023, 7:46am

Thanks for sharing the solution. However, I’m not sure I follow. It would be helpful if you could share the complete @model the didn’t work and the complete @model that resolved the issue.

githubmanojkumar · October 18, 2023, 8:17am

Certainly,
Here is the model that didnt work

Meropenam= @model begin
@param begin
tvcl ∈ RealDomain(lower=0.001)
tvvc ∈ RealDomain(lower=0.001)
tvq ∈ RealDomain(lower=0.001)
tvvp ∈ RealDomain(lower=0.001)
egfronCL ∈ RealDomain(lower=0.001)
Ω ∈ PDiagDomain(4)
σ²_prop ∈ RealDomain(lower=0.001)

end

@random begin
η ~ MvNormal(Ω)
end

@covariates EGFR

@pre begin
Cl = tvcl *(1+(egfr-47)*egfronCL)*exp(η[1])
Vc = tvvc * exp(η[2])
Vp = tvvp * exp(η[3])
Q = tvq * exp(η[4])
end

@dynamics begin
Central’ = -(Cl/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
Peripheral’ = (Q/Vc)*Central - (Q/Vp)*Peripheral
end

@derived begin
CONC = @. Central/Vc
dv ~ @. Normal(CONC, sqrt(CONC^2*σ²_prop))

end
end

params1 = (tvcl =17.26, tvvc = 32.61 ,tvvp =14.5,tvq=15.3, egfronCL= 0.01013, Ω = Diagonal([0.2,0.2,0.2,0.,0.2]), σ²_prop = 0.01)

Random.seed!(123)

egfr1 = rand(90:120,18)
egfr2 = rand(60:89,15)
egfr3 = rand(30:59,33)
egfr4 = rand(15:29,39)
egfr5 = rand(10:15,5)

EGFR = vcat(egfr1,egfr2,egfr3,egfr4,egfr5)

ev1 = DosageRegimen(500, cmt=1, rate = 0.5)
ev2 = DosageRegimen(1000, cmt = 1, rate = 3)

pop1 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr1[i],group=“500mg”),observations = (DV = nothing,),),1:6))
pop2 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr2[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop3 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr3[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop4 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr4[i],group=“500mg”),observations = (DV = nothing,),),1:9))
pop5 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr5[i],group=“500mg”),observations = (DV = nothing,),),1:2))
pop1a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr1[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop2a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr2[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop3a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr3[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop4a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr4[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop5a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr5[i],group=“1000mg”),observations = (DV = nothing,),),1:3))

pop = vcat(pop1,pop2,pop3,pop4, pop5, pop1a, pop2a, pop3a, pop4a, pop5a)

obs1= simobs(Meropenam,pop,params1,obstimes = [0,0.5,2,8])

the model that actually worked
Meropenam= @model begin
@param begin
tvcl ∈ RealDomain(lower=0.001)
tvvc ∈ RealDomain(lower=0.001)
tvq ∈ RealDomain(lower=0.001)
tvvp ∈ RealDomain(lower=0.001)
egfronCL ∈ RealDomain(lower=0.001)
Ω ∈ PDiagDomain(5)
σ² ∈ RealDomain(lower=0.001)

end

@random begin
η ~ MvNormal(Ω)
end

@covariates EGFR

@pre begin
Cl = tvcl *(1+(EGFR-47)*egfronCL)*exp(η[1])
Vc = tvvc * exp(η[2])
Vp = tvvp * exp(η[3])
Q = tvq * exp(η[4])
end

@dynamics begin
Central’ = -(Cl/Vc)*Central - (Q/Vc)*Central + (Q/Vp)*Peripheral
Peripheral’ = (Q/Vc)*Central - (Q/Vp)*Peripheral
end

@derived begin
CONC = @. Central/Vc
dv ~ @. Normal(CONC, sqrt(CONC^2*σ²))

end
end

params1 = (tvcl =17.26, tvvc = 32.61 ,tvvp =14.5,tvq=15.3, egfronCL= 0.01013, Ω = Diagonal([0.5,0.2,0.2,0.2,0.2]), σ² = 0.01)

Random.seed!(123)

egfr1 = rand(90:120,18)
egfr2 = rand(60:89,15)
egfr3 = rand(30:59,33)
egfr4 = rand(15:29,39)
egfr5 = rand(10:15,5)

EGFR = vcat(egfr1,egfr2,egfr3,egfr4,egfr5)

ev1 = DosageRegimen(500, cmt=1, rate = 0.5)
ev2 = DosageRegimen(1000, cmt = 1, rate = 3)

pop1 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr1[i],group=“500mg”),observations = (DV = nothing,),),1:6))
pop2 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr2[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop3 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr3[i],group=“500mg”),observations = (DV = nothing,),),1:3))
pop4 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr4[i],group=“500mg”),observations = (DV = nothing,),),1:9))
pop5 = Population(map(i → Subject(id = i,events = ev1,covariates = ( EGFR=egfr5[i],group=“500mg”),observations = (DV = nothing,),),1:2))
pop1a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr1[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop2a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr2[i],group=“1000mg”),observations = (DV = nothing,),),1:12))
pop3a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr3[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop4a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr4[i],group=“1000mg”),observations = (DV = nothing,),),1:30))
pop5a = Population(map(i → Subject(id = i,events = ev2,covariates = ( EGFR=egfr5[i],group=“1000mg”),observations = (DV = nothing,),),1:3))

pop = vcat(pop1,pop2,pop3,pop4, pop5, pop1a, pop2a, pop3a, pop4a, pop5a)

obs1= simobs(Meropenam,pop,params1,obstimes = [0,0.5,2,8])

If i made any basic silly mistakes in this also, kindly tell me, so that i can correct it next time

andreasnoack · October 18, 2023, 8:29am

Thanks. However, I don’t understand why you’d make it Ω ∈ PDiagDomain(5) since I don’t see η[5] in your code. I can see the in the first version, you wrote

Ω = Diagonal([0.2,0.2,0.2,0.,0.2])

Notice that the fourth element is zero. Should that instead have been

Ω = Diagonal([0.2,0.2,0.2,0.2])

githubmanojkumar · October 19, 2023, 4:16am

Thanks for the update, i will correct it next time.

Topic		Replies	Views
Cholesky factorization failed - Simulation for Population Not Working	2	240	November 16, 2022
Cholesky factorization failed - predict function	2	238	April 4, 2023
Error: LinearAlgebra.SingularException(4) Estimation	10	1991	February 26, 2020
Simulating with vcov Simulation	36	635	May 1, 2024
Problem running example	8	911	August 2, 2019

ERROR: PosDefException

before i Have kept 4 in defining randomeffects but actually the parameters were 5

correction Ω ∈ PDiagDomain(5)**

correction Ω = Diagonal([0.2,0.2,0.2,0.2,0.2])

Related topics