Closed korsbo closed 3 years ago
alanderos91
commented
5 years ago Is ConstantRateJump appropriate here? The variable u[1] could change continuously with the DDE which means the jump rates do too. Using VariableRateJump throws a different error:
ERROR: type DDEProblem has no field lagsand can be traced to https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/problem.jl#L86
korsbo
commented
5 years ago You are right. However, that is a problem with my MWE and I don't think its an issue with my real problem. The system that I am really trying to solve has one variable with is exclusively governed by a Gillespie simulation (with no information taken from any other variable) and the other variables use the first one as input.
In my real problem, I also do not use a fixed time lag but rather calculate the integral of the entire history function. I thus use a different algorithm too (MethodOfSteps(RK4())).
So the MWE may not be perfect, but it results in the same error that my real problem raises.
korsbo
commented
5 years ago For the sake of bug-hunting I think that the MWE is fine. However, the problem I'm actually trying to simulate is much closer to:
using DifferentialEquations
using SpecialFunctions: gamma
using QuadGK
@inline gamma_model(t, p, n, r) = p*r^n*t^(n-1)*exp(-r*t)/gamma(n)
function dde(du,u,h,p,t)
du[1] = 0
integral, err = quadgk(
### Integrand
time -> h(p, time, idxs=1) * gamma_model(t - time, 1, p[2]-2, p[3]),
0., ## t_start
t; ## t_stop
rtol=1e-4
)
du[2] = p[3] * ( p[1] * integral - u[2])
end
h(p, t) = fill(0., 2)
p = [1., 8., 1.4]
tspan = (0., 10.)
u0 = [p[3], 0.]
prob = DDEProblem(dde, u0, h , tspan, p)
### Add some Gillespie action to the first variable
prod_rate(u,p,t) = (1 + 9 * u[1] / (5 + u[1]))/10
prod_affect!(i) = (i.u[1] += 1)
prod_jump = ConstantRateJump(prod_rate, prod_affect!)
deg_rate(u,p,t) = u[1]/10
deg_affect!(i) = (i.u[1] -= 1)
deg_jump = ConstantRateJump(deg_rate, deg_affect!)
jump_prob = JumpProblem(prob, Direct(), prod_jump, deg_jump)
sol = solve(jump_prob, MethodOfSteps(RK4()))
using Plots
plot(sol, vars=2, lw=5, label="DDE model with noisy input.")I just figured that this contained more complexity than the MWE needed.
Ps. for the computational biologist who cares: this model simulates the signal-transformation conferred by a multi-step linear pathway without having to explicitly include every pathway step. I'm hoping that the related paper will be submitted in a week or two.
alanderos91
commented
5 years ago Yeah your MWE is fine. My guess is that the constructor for DDEProblem changed and this was never updated: https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/problem.jl#L86
I'll give this a look but I may not get back with a solution until Monday.
korsbo
commented
5 years ago Thanks! There is no rush on my part. I have already solved my problem using callbacks instead. I just wanted to report the issue.
ChrisRackauckas
commented
5 years ago I tried to run a Jump DDE and, given the supreme coolness of DifferentialEquations.jl
We were almost cool. Then we weren't. Now? 😎
ChrisRackauckas
commented
5 years ago Oh, but we didn't handle the ConstantRateJump 😢
alanderos91
commented
5 years ago I assumed everything was working in the past and did not check for missing jumps.
Here's what I've found (based on the second example):
julia> integrator = init(jump_prob, MethodOfSteps(RK4()))
t: 0.0
u: 2-element Array{Float64,1}:
1.4
0.0
julia> integrator.opts.callback.discrete_callbacks[1].condition.next_jump
0
julia> integrator.opts.callback.discrete_callbacks[1].condition.next_jump_time
InfThis means the CallbackSet is never touched to initialize the jumps, in contrast to ODE/Jump problems. Seeing as how ODE/Jump problems work and that I've made my way down to the solve!(integrator) call, I suspect the error may be outside DiffEqJump. Maybe I missed something?
ChrisRackauckas
commented
5 years ago Maybe DDEs aren't running the callback initialization step at the right point?
alanderos91
commented
5 years ago See the offending LOCs. There are references to dde_int and the bootstrap integrator dde_int.integrator. I'm not sure what L394 is doing here:
u_modified = initialize!(integrator.opts.callback,dde_int.t,u,dde_int)The function call should be initialize!(callbacks, **u**, **t**, integrator). Changing this and replacing references to integrator with dde_int seems to yield the desired behavior:
I'll file a separate issue.
isaacsas
commented
3 years ago @korsbo both examples here seem to go through for me now, but if you are still having issues let me know and I can reopen.
Hi,
I tried to run a Jump DDE and, given the supreme coolness of DifferentialEquations.jl, I rather expected this to just work. However, it did not.
M(not)WE:
I can achieve what I want using callbacks, so I'm not desperate for a fix, but I thought I should report the problem.