-
-
Notifications
You must be signed in to change notification settings - Fork 34
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Jump DDEs error #71
Comments
Is
and can be traced to https://github.com/JuliaDiffEq/DiffEqJump.jl/blob/master/src/problem.jl#L86 |
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 ( So the MWE may not be perfect, but it results in the same error that my real problem raises. |
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. |
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. |
Thanks! There is no rush on my part. I have already solved my problem using callbacks instead. I just wanted to report the issue. |
We were almost cool. Then we weren't. Now? 😎 |
Oh, but we didn't handle the ConstantRateJump 😢 |
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
Inf This means the |
Maybe DDEs aren't running the callback initialization step at the right point? |
See the offending LOCs. There are references to
The function call should be |
@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.
The text was updated successfully, but these errors were encountered: