Inference stops at 10+ jumps

[EddaSchulz]

Hi,

I have some problems to understand why sometimes my simulations become very slow. I set up an example below. There I'm simulating a Gillespie system with 10 reactions. You can find the code below. When I simulate only nine of them (Jumpset all_r), one simulation takes about 0.04sec. As soon as I add the reaction 10 (Jumpset all_r2), which is essentially identical to reaction 9, the simulation takes around 0.6 sec. I don't understand why it becomes about 10 time slower all of a sudden. Where does this come from and what can I do about it?

Thanks a lot for your help!
Edda

Define reactions

rate1(u,p,t) = p[21](0.5(u[5]+u[6]))^p[11]/((0.5*(u[5]+u[6]))^p[11]+(p[22]p[12])^p[11])(1-u[7]^p[13]/(u[7]^p[13]+(p[23]*p[14])^p[13]))
function affect1!(integrator)
integrator.u[1] += 1
end
r1 = ConstantRateJump(rate1,affect1!)

rate2(u,p,t) = p[21](0.5(u[5]+u[6]))^p[11]/((0.5*(u[5]+u[6]))^p[11]+(p[22]p[12])^p[11])(1-u[8]^p[13]/(u[8]^p[13]+(p[23]*p[14])^p[13]))
function affect2!(integrator)
integrator.u[2] += 1
end
r2= ConstantRateJump(rate2,affect2!)

rate3(u,p,t) = p[23]*(1-u[1]^p[5]/(u[1]^p[5]+(p[21]*p[6])^p[5]))
function affect3!(integrator)
integrator.u[7] += 1
end
r3= ConstantRateJump(rate3,affect3!)

rate4(u,p,t) = p[23]*(1-u[2]^p[5]/(u[2]^p[5]+(p[21]*p[6])^p[5]))
function affect4!(integrator)
integrator.u[8] += 1
end
r4= ConstantRateJump(rate4,affect4!)

rate5(u,p,t) = p[22]*(1-u[1]^p[3]/(u[1]^p[3]+(p[21]*p[4])^p[3]))
function affect5!(integrator)
integrator.u[5] += 1
end
r5= ConstantRateJump(rate5,affect5!)

rate6(u,p,t) = p[22]*(1-u[2]^p[3]/(u[2]^p[3]+(p[21]*p[4])^p[3]))
function affect6!(integrator)
integrator.u[6] += 1
end
r6= ConstantRateJump(rate6,affect6!)

rate7(u,p,t) = 1*u[1]
function affect7!(integrator)
integrator.u[1] -= 1
end
r7= ConstantRateJump(rate7,affect7!)

rate8(u,p,t) = 1*u[2]
function affect8!(integrator)
integrator.u[2] -= 1
end
r8= ConstantRateJump(rate8,affect8!)

rate9(u,p,t) = 1*u[7]
function affect9!(integrator)
integrator.u[7] -= 1
end
r9= ConstantRateJump(rate9,affect9!)

rate10(u,p,t) = 1*u[8]
function affect10!(integrator)
integrator.u[8] -= 1
end
r10= ConstantRateJump(rate10,affect10!)

all_r=JumpSet(r1,r2,r3,r4,r5,r6,r7,r8,r9)
all_r2=JumpSet(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10)

Function to simulate one cell

function one_cell(p::Array{Float64,1},u0::Array{Int64,1},tspan::Tuple{Float64,Float64},r::DiffEqJump.JumpSet)
prob = DiscreteProblem(u0,tspan,p)
jump_prob = JumpProblem(prob,Direct(),r)
sol = solve(jump_prob,FunctionMap())
return sol
end

Perform simulation

p=[0,0,3,0.11,3,0.11,1,1,0,0,3,0.6,3,0.46,1,1,1,1,1,1,100,100,100]
tspan=(0.0,100.0)
x0 = zeros(10)
x0[1:8] = [100,1,0,0,p[22],p[22],p[23],p[23]]
x0_int = convert(Array{Int64},x0)

simulate 9 reactions

@time s=one_cell(p,x0_int,tspan,all_r)

simulate 10 reactions

@time s=one_cell(p,x0_int,tspan,all_r2)

[Armavica]

Hi, I don't know about your issue but you might be interested to know that surrounding your code with three backticks ( ``` ) in a GitHub issue gives it a nice formatting, to which you can even add syntactic coloring if you specify the language:

``` Julia
# simulate 9 reactions
@time s = one_cell(p, x0_int, tspan, all_r)
```

will render as:

# simulate 9 reactions
@time s = one_cell(p, x0_int, tspan, all_r)

It makes it easier to copy & paste your code to reproduce your problem, as it prevents some characters from being interpreted as formatting instructions by GitHub (like '#' at the beginning of a line).

[EddaSchulz]

Hi, thank you!
So here is my code with better formatting...

# Define reactions
rate1(u,p,t) = p[21]*(0.5*(u[5]+u[6]))^p[11]/((0.5*(u[5]+u[6]))^p[11]+(p[22]*p[12])^p[11])*(1-u[7]^p[13]/(u[7]^p[13]+(p[23]*p[14])^p[13]))
function affect1!(integrator)
  integrator.u[1] += 1
end
r1 = ConstantRateJump(rate1,affect1!)

rate2(u,p,t) = p[21]*(0.5*(u[5]+u[6]))^p[11]/((0.5*(u[5]+u[6]))^p[11]+(p[22]*p[12])^p[11])*(1-u[8]^p[13]/(u[8]^p[13]+(p[23]*p[14])^p[13]))
function affect2!(integrator)
  integrator.u[2] += 1
end
r2= ConstantRateJump(rate2,affect2!)

rate3(u,p,t) = p[23]*(1-u[1]^p[5]/(u[1]^p[5]+(p[21]*p[6])^p[5]))
function affect3!(integrator)
  integrator.u[7] += 1
end
r3= ConstantRateJump(rate3,affect3!)

rate4(u,p,t) = p[23]*(1-u[2]^p[5]/(u[2]^p[5]+(p[21]*p[6])^p[5]))
function affect4!(integrator)
  integrator.u[8] += 1
end
r4= ConstantRateJump(rate4,affect4!)

rate5(u,p,t) = p[22]*(1-u[1]^p[3]/(u[1]^p[3]+(p[21]*p[4])^p[3]))
function affect5!(integrator)
  integrator.u[5] += 1
end
r5= ConstantRateJump(rate5,affect5!)

rate6(u,p,t) = p[22]*(1-u[2]^p[3]/(u[2]^p[3]+(p[21]*p[4])^p[3]))
function affect6!(integrator)
  integrator.u[6] += 1
end
r6= ConstantRateJump(rate6,affect6!)

rate7(u,p,t) = 1*u[1]
function affect7!(integrator)
  integrator.u[1] -= 1
end
r7= ConstantRateJump(rate7,affect7!)
rate8(u,p,t) = 1*u[2]
function affect8!(integrator)
  integrator.u[2] -= 1
end
r8= ConstantRateJump(rate8,affect8!)
rate9(u,p,t) = 1*u[7]
function affect9!(integrator)
  integrator.u[7] -= 1
end
r9= ConstantRateJump(rate9,affect9!)
rate10(u,p,t) = 1*u[8]
function affect10!(integrator)
  integrator.u[8] -= 1
end
r10= ConstantRateJump(rate10,affect10!)

all_r=JumpSet(r1,r2,r3,r4,r5,r6,r7,r8,r9)
all_r2=JumpSet(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10)

# Function to perform simulation
function one_cell(p::Array{Float64,1},u0::Array{Int64,1},tspan::Tuple{Float64,Float64},r::DiffEqJump.JumpSet)
  prob = DiscreteProblem(u0,tspan,p)
  jump_prob = JumpProblem(prob,Direct(),r)
  sol = solve(jump_prob,FunctionMap())
  return sol
end

# simulation
p=[0,0,3,0.11,3,0.11,1,1,0,0,3,0.6,3,0.46,1,1,1,1,1,1,100,100,100]
tspan=(0.0,100.0)
x0 = zeros(10)
x0[1:8] = [100,1,0,0,p[22],p[22],p[23],p[23]]
x0_int = convert(Array{Int64},x0)
@time s=one_cell(p,x0_int,tspan,all_r)
@time s=one_cell(p,x0_int,tspan,all_r2)

[ChrisRackauckas]

It's related to Julia's inference bug:

JuliaLang/julia#22255
JuliaLang/julia#22370

We have a separate jump problem design coming for this case: SciML/JumpProcesses.jl#20 . Basically, it'll use vectors of FunctionWrappers instead of fully-inlined tuple recursion, since the latter abuses inference but is only allowed to up to a small finite number of jumps.

[EddaSchulz]

Sorry, but I don't understand how to use this. Is there an example somewhere?
Thank you!
Edda

[ChrisRackauckas]

It's not released yet. It still needs improvements, but we merged some PRs yesterday so it's actively under development.

[ChrisRackauckas]

See the updated tutorial

http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html

along with the new jump aggregators

http://docs.juliadiffeq.org/latest/types/jump_types.html#Constant-Rate-Jump-Aggregators-1