I sent you one by direct mail some time ago, in private through this channel -
If you don’t have it - I can make it public again - Please let me know 
A
Here is the code - thanks
I don’t have a toy example so I’ll share the full code with you. Hope you don’t mind .
could the issue be that my server runf julia 1.0 and not 1.1?
I am happy to share it with you but I’ll have to take it down soon as this is not public yet - hope that’s not an issue.
thanks,
A
using Distributed
using SymPy
#using Sundials
using DifferentialEquations
@everywhere using DiffEqSensitivity
#using ParameterizedFunctions # not necesary
#using LinearAlgebra
using Plots
using Unitful
using LaTeXStrings
#using DelimitedFiles
using JLD2
using IterableTables
using DataFrames
using CSV
# system of equations
function PKPD3(du,u,p,t)
#=
X1 = X1_ADC
X2 = X2_ADC
X3 = ADC_Tumor_Free_Extracellular
X4 = ADC_Tumor_Bound_Extracellular
X5 = ADC_Tumor_Intracellular
X6 = PL_Tumor_Intracellular
X7 = PL_Bound_Tubulin
X8 = PL_Tumor_Extracellular
X9 = C1_PL
X10 = C2_PL
X11 = DAR
X12 = V1
X13 = V2
X14 = V3
X15 = V4
X16 = TV = V1 + V2 + V3 + V4
X17 = Growth
X18 = Kill
X19 = R_tum
X20 = Diff_adc
X21 = Diff_drug
=#
(CL_adc, CLD_adc, V1_adc, V2_adc, Kdec, e_adc, P_adc, D_adc,
Kon_adc, Koff_adc, Agtot, Kint_ag, Kdeg, Kdiff_drug, Kout_drug,
Kon_tub, Tub, Koff_tub, e_drug, P_drug, D_drug,
CL_drug, CLD_drug, V1_drug, V2_drug,
Vmax, τ, Kgex, Kglin, ψ, KKillmax, KC50, R_cap,R_kro) = p
# Define varaibles:
X1 = u[1] ; dX1 = du[1]
X2 = u[2] ; dX2 = du[2]
X3 = u[3] ; dX3 = du[3]
X4 = u[4] ; dX4 = du[4]
X5 = u[5] ; dX5 = du[5]
X6 = u[6] ; dX6 = du[6]
X7 = u[7] ; dX7 = du[7]
X8 = u[8] ; dX8 = du[8]
X9 = u[9] ; dX9 = du[9]
X10 = u[10] ; dX10 = du[10]
X11 = u[11] ; dX11 = du[11]
DAR = X11
X12 = u[12] ; dX12 = du[12]
X13 = u[13] ; dX13 = du[13]
X14 = u[14] ; dX14 = du[14]
X15 = u[15] ; dX15 = du[15]
X16 = u[16]
X17 = u[17]
X18 = u[18]
X19 = u[19]
X20 = u[20]
X21 = u[21]
TV = X16
Growth = X17
Kill = X18
R_tum = X19
Diff_adc= X20
Diff_drug = X21
# ADC 2 Compartment Equations
du[1] =-((CL_adc/V1_adc) + (CLD_adc/V1_adc))*X1 +
(CLD_adc/V2_adc)*X2 - ((e_adc*X1/V1_adc -X3 )*TV*((2*P_adc*R_cap/(R_kro^2))+Diff_adc))#- Kdec*X1 - X3
du[2] = (CLD_adc/V1_adc)*X1 - (CLD_adc/V2_adc)*X2
# ADC Cellular Equations
du[3] = (e_adc*(X1/V1_adc)- X3)*((2*P_adc*R_cap/(R_kro^2))+Diff_adc) +
- Kon_adc*X3*((Agtot - X4)/e_adc) + Koff_adc*X4#- Kdec*X3
du[4] = Kon_adc*X3*((Agtot - X4)/e_adc) +
-(Koff_adc+Kint_ag)*X4 # -Kdec*X4 -dX4
du[5] = Kint_ag*X4 - Kdeg*X5
# Drug Cellular Equations
du[6] = Kdeg*X5*DAR + Kdiff_drug*(X8-X6) - Kout_drug*X6 +
- Kon_tub*X6*(Tub - X7) + Koff_tub*X7
du[7] = Kon_tub*X6*(Tub - X7) - Koff_tub*X7
du[8] = (e_drug*X9 - X8)*((2*P_drug*R_cap/(R_kro^2)) + Diff_drug) +
- Kdiff_drug*(X8 - X6) + Kout_drug*X6 + DAR*Kdec*(X3 + X4)
# Drug 2 Compartment Equations
du[9] = -(CL_drug/V1_drug)*X9 - (CLD_drug/V1_drug)*X9 +
+ (CLD_drug/V1_drug)*X10 + (X1*DAR*Kdec)/V1_drug +
+ (CL_adc*DAR*(X1/V1_adc))/V1_drug +
- (e_drug*X9 - X8)*((2*P_drug*R_cap/(R_kro^2)) + Diff_drug) #*TV
du[10] = (CLD_drug/V2_drug)*X9 - (CLD_drug/V2_drug)*X10
# ADC deacy Equation
du[11] = -Kdec*X11
# PD Equations
du[12] = (Growth-Kill)*X12
du[13] = Kill*X12 - (1/τ)*X13
du[14] = (1/τ)*(X13 - X14)
du[15] = (1/τ)*(X14 - X15)
# Equations 16 and on are algebraic
du[16] = (X12 + X13 + X14 + X15)-TV
du[17] = Kgex*(1-(TV/Vmax))/(1+abs((Kgex*TV)/Kglin)^ψ)^(1/ψ)-Growth
du[18] = KKillmax*(X6+X7)/(KC50 + X6+X7)-Kill
# Dynamic Parameters
du[19] = cbrt((3/(4π))*abs(TV))*100000 - R_tum
du[20] = 6*D_adc/(abs(R_tum)^2) -Diff_adc
du[21] = 6*D_drug/(abs(R_tum)^2) - Diff_drug
end
# parameters
weight = 20.0u"g"
MW_adc = uconvert(u"g/nmol",153000u"g/mol")
MW_drug = uconvert(u"g/nmol",717.993u"g/mol")
Dose = uconvert(u"g/g",2u"mg/kg")*weight/MW_adc
CL_adc = 0.025u"L/(d*kg)"*uconvert(u"kg",weight)
CLD_adc = 0.036u"L/(d*kg)"*uconvert(u"kg",weight)
V1_adc = 0.054u"L/kg"*uconvert(u"kg",weight)
V2_adc = 0.1u"L/kg"*uconvert(u"kg",weight)
Kdec = 0.072u"1/d"
e_adc = 0.241
P_adc = 334u"μm/d"
D_adc = 0.022u"cm^2/d"
D_adc = uconvert(u"μm^2/d",D_adc)
Kon_adc = 6.22u"L/(nmol*d)"
Koff_adc = 12.5u"1/d"
Agtot = 166u"nmol/L"
Kint_ag = 0.61u"1/d"
Kdeg = 0.72u"1/d"
Kdiff_drug = 9.66u"1/d"
Kout_drug = 1.1u"1/d"
Kon_tub = 0.44u"L/(nmol*d)"#6.22u"L/(nmol*d)"
Tub = 65u"nmol/L"
Koff_tub = 13.1u"1/(d)"
e_drug = 0.44
P_drug = 2.1e4u"μm/d"
D_drug = 0.25u"cm^2/d"
D_drug = uconvert(u"μm^2/d",D_drug)
CL_drug = 0.93u"L/(d*kg)"*uconvert(u"kg",weight)
CLD_drug = 3.03u"L/(d*kg)"*uconvert(u"kg",weight)
V1_drug = 1.26u"L/kg"*uconvert(u"kg",weight)
V2_drug = 0.799u"L/kg"*uconvert(u"kg",weight)
Vmax = 500u"mm^3"
Vmax = uconvert(u"L",Vmax)
τ = 0.128u"d"
Kgex = 0.245u"1/(d)"
Kglin = 43.3u"mm^3/(d)"
Kglin = uconvert(u"L/d",Kglin)
ψ = 20
KKillmax = 0.761u"1/(d)"
KC50 = 100u"nmol/L"
R_cap = 8u"μm"
R_kro = 75u"μm"
#R_tum = 200u"μm"
ITV = 100u"mm^3"
ITV = uconvert(u"L",ITV)
IR_tum = cbrt((3/(4π))*(uconvert(u"μm^3",ITV)))
IDiff_adc = 6*D_adc/(IR_tum^2)
IDiff_drug = 6*D_drug/(IR_tum)^2
Perm_adc = (2*P_adc*R_cap/(R_kro^2))
Perm_drug = (2*P_drug*R_cap/(R_kro^2))
p=ustrip.([CL_adc, CLD_adc, V1_adc, V2_adc,
Kdec, e_adc, P_adc, D_adc, Kon_adc, Koff_adc,
Agtot, Kint_ag, Kdeg, Kdiff_drug, Kout_drug, Kon_tub, Tub, Koff_tub,
e_drug, P_drug, D_drug, CL_drug, CLD_drug, V1_drug, V2_drug, #, TV])
Vmax, τ, Kgex, Kglin, ψ, KKillmax, KC50, R_cap, R_kro])
;
#IC
IX1 = Dose
IX2 = 0u"nmol"
IX3 = 0u"nmol/L"
IX4 = 0u"nmol/L"
IX5 = 0u"nmol/L"
IX6 = 0u"nmol/L"
IX7 = 0u"nmol/L"
IX8 = 0u"nmol/L"
IX9 = 0u"nmol/L"
IX10 = 0u"nmol/L"
IX11 = 4
IX12 = ITV
IX13 = 0u"L"
IX14 = 0u"L"
IX15 = 0u"L"
IX16 = ITV
IX17 = Kgex*(1-(ITV/Vmax))/(1+((Kgex*ITV)/Kglin)^ψ)^(1/ψ)
IX18 = KKillmax*(IX6+IX7)/(KC50 + IX6+IX7)
IX19 = IR_tum
IX20 = IDiff_adc
IX21 = IDiff_drug
IC = ustrip.([IX1, IX2, IX3, IX4, IX5, IX6, IX7, IX8, IX9,
IX10, IX11, IX12, IX13, IX14, IX15, IX16, IX17, IX18, IX19, IX20, IX21])
IC2 = copy(IC)
IC2[1]=IC[1]/4
# tspan and others
tspan = ustrip.((0.0u"d",100u"d"))
tspan2 = ustrip.((0.0u"d",28u"d"))
# diff_vars = [true,true,true,true,true, true, true, true,true,true, true, true, true, true, true,
# false, false, false, false, false, false]
varnum=21
algnum=6
dv=varnum - algnum
diff_vars=[if i≤dv true else false end for i in 1:varnum ]
MM = zeros(varnum,varnum)
[MM[i,i] = 1 for i in 1:dv]
;
# Solve equations
PKPD3_test = ODEFunction(PKPD3;mass_matrix=MM)
PKPD3_Prob = ODEProblem(PKPD3_test,IC,tspan,p)
PKPD3_sol = DifferentialEquations.solve(PKPD3_Prob,Rodas5())#;dt=1//2^(5))#, alghint=:stiff)
# convert solution to dataframe
DataFrame(PKPD3_sol)
I know it is a lot of text, but I am not sure why it doesn’t work so I rather share the actual code other than make a toy sample -
thanks
a-
That’s not a good excuse. The best way to find out why it doesn’t work is to make a toy example.
using OrdinaryDiffEq
f_2dlinear = (du,u,p,t) -> du.=1.01u
prob = ODEProblem(f_2dlinear,rand(2,2),(0.0,1.0))
sol1 =solve(prob,Tsit5())
using IterableTables, DataFrames
df = DataFrame(sol1)
works while
using OrdinaryDiffEq
f_2dlinear = (du,u,p,t) -> du.=1.01u
prob = ODEProblem(f_2dlinear,rand(2),(0.0,1.0))
sol1 =solve(prob,Tsit5())
using IterableTables, DataFrames
df = DataFrame(sol1)
fails. Probably an issue with IterableTables.jl. I opened an issue: Forcing the most specialized iterator · Issue #92 · queryverse/IterableTables.jl · GitHub
Ok-
My apologies.
I was really not sure if the issue was related to my code- that’s why I wanted to share the the real problem.
A
Next time I’ll make sure to simplify.
As always thanks for the help!
A
Hi Chirs-
I hope you are doing well!
Sorry to bother you again-
Is measurement type compatible with DAEs?
I wrote the following test code:
# Equation
function test(du,u,p,t)
du[1]=-p[1]*u[1]
du[2]=-(p[2])*u[2]
du[3]=2*u[2] - u[1]
end
# test = @ode_def begin
# dx=-k*x
# end k
x0 = [10.0 ± 1, 5.0 ± 1, 1±0]
k=1 ± 0.1
k2=5 ± 0.05
# Solve Equation
# Mass=generate_MassMatrix(3,1)
# Mass=[1±0 0.0 0.0; 0.0 1±0 0.0; 0.0 0.0 0.0]
Mass=[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0]
System = ODEFunction(test;mass_matrix=Mass)
prob = ODEProblem(System,x0,(0.0,10.0),[k,k2])
condition(u,t,integrator) = t in [2,4,6,8] # && u[1]/V<1
affect!(integrator) = integrator.u[1] += 10±1
affect2!(integrator) = integrator.u[2] += 5±1
cb = DiscreteCallback(condition,affect!;save_positions=(false,false))
cb2 = DiscreteCallback(condition,affect2!;save_positions=(false,false))
cbset=CallbackSet(cb,cb2)
sol = solve(prob,Rodas5())#,callback=cbset,tstops=[2,4,6,8])
If I use Tsit() and no mass_matrix- the problems solves.
However once I change the solver to Rodas5() or other Rosenbrock methods to make the problem into a DAE - the code fails. This happens even if I use Rodas to solve the system as ODE (no mass matrix).
Is the issue with my code or just DAE/ Rodas solvers do not work with measurement type, Is there any solution or solver I could use to actually solev this problem?
The error message seems to indicate that the issue is with forward differentiation, but that I don;t know how to fix.
Thanks,
A
Truncated Error:
MethodError: *(::Measurement{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Measurement{Float64},1}},Measurement{Float64}},Measurement{Float64},3}) is ambiguous. Candidates:
*(a::Measurement, b::Real) in Measurements at /opt/julia/packages/Measurements/VjRMd/src/math.jl:197
*(x::AbstractFloat, y::ForwardDiff.Dual{Ty,V,N} where N where V) where Ty in ForwardDiff at /opt/julia/packages/ForwardDiff/N0wMF/src/dual.jl:140
Possible fix, define
*(::Measurement, ::ForwardDiff.Dual{Ty,V,N} where N where V)
Stacktrace:
[1] test(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Measurement{Float64},1}},Measurement{Float64}},Measurement{Float64},3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Measurement{Float64},1}},Measurement{Float64}},Measurement{Float64},3},1}, ::Array{Measurement{Float64},1}, ::Float64) at ./In[39]:3
[2] (::ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Measurement{Float64},1}},Measurement{Float64}},Measurement{Float64},3},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(test),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Measurement{Float64},1}},Measurement{Float64}},Measurement{Float64},3},1}, ::Vararg{Any,N} where N) at /opt/julia/packages/DiffEqBase/6ewdP/src/diffeqfunction.jl:188
[3]
Try sol = solve(prob,Rodas5(autodiff=false))
That worked - once I chabged the toy example to one that wasn’t singular
Thanks,
A
so this worked well for the toy system - once I tried to apply to a bigger system it doesn’t.
Do you think it might be an issue with error porpagation - The sstem solves perfectly with Float64.
Thanks,
Ale
Here is the error
MethodError: no method matching gemm!(::Char, ::Char, ::Measurement{Float64}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::Measurement{Float64}, ::SubArray{Measurement{Float64},2,Array{Measurement{Float64},2},Tuple{UnitRange{Int64},UnitRange{Int64}},false})
Closest candidates are:
gemm!(::AbstractChar, ::AbstractChar, !Matched::Float64, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, !Matched::Float64, !Matched::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109
gemm!(::AbstractChar, ::AbstractChar, !Matched::Float32, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}, !Matched::Float32, !Matched::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109
gemm!(::AbstractChar, ::AbstractChar, !Matched::Complex{Float64}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}, !Matched::Complex{Float64}, !Matched::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}) at /opt/julia-1.0.0/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:1109
It’s a missing method gemm! for measurements. Please do not continue one thread for so many different questions. It is not inviting other people to answer.
OK -
Last question here!
Is this something I need to install? update?
Is there a version of Measurements that you use?
Or just not available?
Thanks,
Ale
PS-I’ll start new thread after this.
It’s probably a method missing in Measurements.jl. You should either open an issue in the repo mentioning that gemm! is missing for matrices of measurements, or open a PR to Measurements.jl that fixes it.
Thanks
Next question will be in new thread