The other question in the same general subject is if I have a DAE system with callbacks
would something in this general terms work?
# Example Callbacks
condition(u,t,integrator) = t in ([21, 42, 63, 84])
function affect!(integrator)
integrator.u[1] += Dose
integrator.u[11] = (4*Dose+integrator.u[11]*(integrator.u[1]+integrator.u[2]+(integrator.u[3]+
integrator.u[4]+integrator.u[5])*integrator.u[16]))/(Dose+(integrator.u[1]+
integrator.u[2]+(integrator.u[3]+integrator.u[4]+integrator.u[5])*integrator.u[16]))
end
condition1(u,t,integrator) = u[19] < 8
affect1!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition,affect!)
cb1 = DiscreteCallback(condition1, affect1!)
cbset = CallbackSet(cb1,cb)
# Example model
#p=is my parameters a vcetor of 36 parameters
MM = zeros(21,21)
[MM[i,i] = 1 for i in 1:15];
t = 0.0:1.0:100
tspan = (0.0,100)
ODESystem = ODEFunction(PKPD4;mass_matrix=MM)
model_ode(p_) = ODEProblem(ODESystem, IC, tspan,p_)
solve_model(mp_) = DifferentialEquations.solve(model_ode(mp_), Rodas5(),saveat=1, callback=cbset, tstops=[21,42,63,84])
mock_data = Array(solve_model(p))
### Building loss function (here is where I am most confused)
###Is it here where the callbacks would fit in the loss_objective function?
### I am having trouble figuring this out
loss_objective(mp_, dat) = build_loss_objective(model_ode(mp_), Rodas5(), callback=cbset, tstops=[21,42,63,84], L2Loss(t,dat),maxiters=100000)
juobj(args...) = loss_objective(args, mock_data)(args)
# JuMP model to get parameters
### What algorithm do you recommend for large non-linear problems? I am using :GN_ISRES, but maybe there are better options?
jumodel = Model()
JuMP.register(jumodel, :juobj, 36, juobj, autodiff=true)
@variables jumodel begin
#p1 to 36
end
@NLobjective(jumodel, Min, juobj(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_Tumor, 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, c, RCell))
setsolver(jumodel, NLoptSolver(algorithm=:GN_ISRES))
From what you see in this sort of pseudo code - Does it make sense to you? I’ve been trying to test it, but it takes long long time to run, so I am honestly not sure if my code is failing or its just slow.
I did not use the callbacks in my test to make things simpler on first pass.
The code structure seems to make sense (at least to me), but I am concerned I am not inputing the mass matrices or callbacks properly into the loss function.
Thanks,
A