How to bring differential equations to wanted value with DiffEqFlux?

Something weird happened. Until the day before yesterday, when I was introducing 200 million phages, the system went into a cyclic phase and then reached an equilibrium:


Yesterday I ran this code:

parms = [mu, kappa, phi, omega, eta, beta, Vp]
# modification
condition1(u, t, integrator) = t==Tϕ
affect1!(integrator) = integrator.u[3] += integrator.p[7]
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u[1]-1
function affect2!(integrator)
  integrator.u[1] = 0
  integrator.u[2] = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u[3]-1
function affect3!(integrator)
  integrator.u[3] = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(sir!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp], sensealg=ForwardDiffSensitivity())

and, for the same amount I got the extinction:


Since the amount inserted is the same, I thought the solver itself has changed. So, I closed Julia’s core and started a new session but the results are the same.
What happened?
Can the solver be reversed as it was?
I think the question I have (found the amount of phage that extinguishes the bacterium) cannot be found mathematically, but from the other post (obtained with the solver as it was) it can be seen that there must be an amount that triggers such extinction.
So, I’ll rephrase the question: is there a computational savvy way to run the sir function in a loop providing different Vp?
Thanks

You can change versions at any time. Are you sure it’s a version thing?

https://diffeq.sciml.ai/stable/features/ensemble/

1 Like

I don’t think is a solver’s version problem, for I did not install a new version of DifferentialEquations: I just added BoundaryValueDiffEq, Flux, Optim, DiffEqFlux, DiffEqSensitivity. But how can it be that now the solver brings the system to zero even when using a fresh session but before it did not? I assumed that the code I used changed the solver in a stable manner (if possible)…

repro?

Hello, I don’t know what “repro” means but I tried again. If I give the value Vp = 200 million phages, I now get extinction, whereas before that was not the case. Moreover, If I use up to Vp/21 I still get extinction but I don’t with amounts below Vp/22. Standing that I don’t understand why I was getting a different result before running DiffEqFlux, it is still evident that the minimal amount of phages is not Vp but another value that DiffEqFlux did not help finding. I’ll try the ensemble simulations you’ve suggested. Thanks

Hello, how do i set the code for the parallel ensemble simulation? I can see that:
EnsembleProblem(prob::DEProblem; output_func = (sol, I),-> (sol, false), prob_func = (prob, i, repeat) -> (prob), reduction = (u, data, I) -> (append!(u, data), false) where prob is the one I have set in the first post. Is this verbatim OK? and how do I tell what parameter to check and in what range? The example in the link gives the code function prob_func(prob, I, repeat) @. [...] but I don’t understand what it does.

You just say how to change the base problem to new ones. So, remake it with new parameters.

Hello, could you please write an example of how should I set the function? I tried to focalize the problem with this post but I got no answer… Thank you

https://diffeq.sciml.ai/stable/features/ensemble/#Random-Initial-Conditions

I have seen that but I don’t know how to use it…

What’s the exact question here? There’s plenty of examples already, so without a simple exact question it’s hard to do much more. It’s literally just a function that spits out the problem to solve for trajectory i.

I have been looking at example 2 of the vignette since it covers a Lotka-Volterra system. This example has a main function that resembles the one I called growth in my second post. The problem is that the example than sets another function:

function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
end

What is the use of this function in my case?
I just want to test over u[3] in growth. How do I write that down?
shall I do:

function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
  du[3] = p[5]*u[3]
end

or simply

function g(du,u,p,t)
  du[3] = p[1]*u[3]
end

What about the indices?

What do you even mean by this? What are you trying to do mathematically? Do you mean you want to sample over the parameters of the noise? Then yes, make there be a parameter p that is used in the diffusion equation.

That’s completely dependent on the model. Is that the same parameter that is used in the drift equation f? Then use the same index. Is it a different one? Then use a different index. It’s literally just the same p that you pass in. There’s nothing fancy there.

I just need to pass different values to u[3] at the selected step. I got:

condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += i

I just want to test several combinations of Tp and i.

condition1(u, t, integrator) = t==integrator.p[4]
affect1!(integrator) = integrator.u[3] += integrator.p[5]

and follow the examples.

I modified as you suggested, but how do I write the solver?

julia> condition1(u, t, integrator) = t==integrator.p[4]         ###
condition1 (generic function with 1 method)
julia> affect1!(integrator) = integrator.u[3] += integrator.p[5] ###
affect1! (generic function with 1 method)
julia> cb1 = DiscreteCallback(condition1, affect1!)
DiscreteCallback{typeof(condition1),typeof(affect1!),typeof(DiffEqBase.INITIALIZE_DEFAULT),typeof(DiffEqBase.FINALIZE_DEFAULT)}(condition1, affect1!, DiffEqBase.INITIALIZE_DEFAULT, DiffEqBase.FINALIZE_DEFAULT, Bool[1, 1])
...
julia> prob = ODEProblem(growth!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 4000.0)
u0: [1.0e7, 0.0, 0.0]

julia> EnsembleProblem(prob;
                       output_func = (sol,i) -> (sol,false),
                       prob_func= (prob,i,repeat)->(prob),
                       reduction = (u,data,I)->(append!(u,data),false),
                       u_init = [], safetycopy = prob_func !== DEFAULT_PROB_FUNC)
ERROR: UndefVarError: prob_func not defined
Stacktrace:
 [1] top-level scope at none:1

julia> solve(ensembleprob,alg,EnsembleThreads();trajectories=1000)
ERROR: UndefVarError: ensembleprob not defined
Stacktrace:
 [1] top-level scope at none:1

Here I think the problem is that I do not need prob_func but simply pass an array with the values to test for both p[4] and p[5]… How would I do that? Thanks

Don’t copy the kwargs.

I tried with:

julia> prob = ODEProblem(growth!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 4000.0)
u0: [1.0e7, 0.0, 0.0]

julia> EnsembleProblem(prob;
                       output_func = (sol,i) -> (sol,false),
                       reduction = (u,data,I)->(append!(u,data),false))
EnsembleProblem with problem ODEProblem

julia> soln = solve(prob,EnsembleThreads();trajectories=1000)
retcode: Success
Interpolation: automatic order switching interpolation
t: 33-element Array{Float64,1}:
    0.0
    0.15303314059544662
...

Looks like it starts to work… Thank you

1 Like

An update: the process did not crash, but how should it be interpreted? I soln is a dataframe of 3 columns and 33 rows. Where do all the combinations come into play? If I do EnsembleSummary(soln) I get the error:

┌ Warning: `MonteCarloSummary(args...)` is deprecated, use `EnsembleSummary(args...)` instead.
│   caller = ip:0x0
└ @ Core :-1

and then the core stops the execution after few hours. The plot of soln columns gives:


How shall I read the results from EnsembleProblem? Thanks

I have no idea what the question is. Start a new thread with an MWE.