EnsembleProblem change tstops for each trajectory

Hey,

I’am trying to use the EnsembleProblem to run an ODE with a big parameter set. Following the problem in the documentation (Event Handling and Callback Functions · DifferentialEquations.jl example 1) I want to follow the patients drug concentration over time, except from single doses at specific time points I am giving doses over a period of time. I implemented for this two PresetTimeCallbacks, one for activating the dosage and one fore deactivating the dosage. My problem now is that I want to have different treatment schedules for each trajectory, i.e. for each trajectory the treatment starts and ends at different time points. I don’t manage to update the treatment schedules for each trajectory, it always stays the same. I tried to set up a small example:

using DifferentialEquations, DiffEqCallbacks,Plots
function f(du,u,p,t)
du[1] = -u[1]+p
end
u0 = [10.0]

function affect_treatment_monoA!(integrator)
integrator.p=a
end
function affect_notreat!(integrator)
integrator.p=0.0
end

p=0.0

dosage=[10.0, 5.0]
rndsched=Array{Any}(undef,2)
rndsched[1]=[[1.0 ,8.0], [3.0 ,9.0]] #treatment start and end of first trajectory
rndsched[2]=[[5.0] ,[7.0] ]
#treatment start and end of second trajectory
treatStart=rndsched[1][1]
treatEnd=rndsched[1][2]
treatMonoSchedA=PresetTimeCallback(treatStart,affect_treatment_monoA!)
NOtreatSched=PresetTimeCallback(treatEnd,affect_notreat!)

function prob_func_Mono1(prob,i,repeat)
global treatStart=rndsched[i][1]
global treatEnd=rndsched[i][2]
global a=dosage[i]
remake(prob,u0=u0,p=p,callback=CallbackSet(treatMonoSchedA,NOtreatSched))
end

prob = ODEProblem(f,u0,(0.0,10.0),0.0,callback=CallbackSet(treatMonoSchedA,NOtreatSched))
emble_prob_mono1=EnsembleProblem(prob,prob_func=prob_func_Mono1)
solMonoSched=solve(emble_prob_mono1,EnsembleThreads(),trajectories=2)
plot(solMonoSched[1],label=“trajectory1”)
plot(solMonoSched[2], label=“trajectory2”)

(Also the update of the dosage parameter a does not seem to work in my small example)
For any hints for what I am doing wrong I would be grateful.

1 Like

Thanks i figured it out! The the whole PresetTimeCallback should be redefined within the prob_func like:

treatMonoSchedA=PresetTimeCallback(rndsched[1][1],affect_treatment_monoA!)
NOtreatSched=PresetTimeCallback(rndsched[1][2],affect_notreat!)

function prob_func_Mono1(prob,i,repeat)
treatMonoSchedA=PresetTimeCallback(rndsched[i][1],affect_treatment_monoA!)
NOtreatSched=PresetTimeCallback(rndsched[i][2],affect_notreat!)
global a=dosage[i]
remake(prob,u0=u0,p=p,callback=CallbackSet(treatMonoSchedA,NOtreatSched))
end
(I only still need to figure out how to update the dosage for each trajectory)