Passing parameters to discrete problem.

[EddaSchulz]

I'm trying to solve a jump Problem and I would like to pass parameters to the rates. I however cannot find a way to do it. In the tutorial there is only and example for ODE with paramters...

Any help would be appreciated.

[ChrisRackauckas]

Hey, If you are defining the rates directly, then the rate equations can use parameters. For example:

rate(u,p,t) = (0.1/1000.0)*u[1]*u[2]
function affect!(integrator)
  integrator.u[1] -= 1
  integrator.u[2] += 1
end
jump = ConstantRateJump(rate,affect!)

is what's in the documentation, but you can make use of p in this equation. (Ref: http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html#Defining-the-Jumps-Directly-1 )

Note that currently the reaction network DSL is unable to handle parameters, but we are in the process of fixing that, and hopefully that will get released in the next few weeks. When that's released, we will update that tutorial to show it using parameters like the others.

[EddaSchulz]

Thanks for the advice, this is working now. I stumbled over another problem though. I can simulate my system (see code below), but when I try to plot the results with

plot(sol)

it takes forever and uses a large amount of memory (>8GB). What is wrong here?

Thanks for your help! All the best Edda

using DifferentialEquations using Plots

rate1(u,p,t) = p[7](0.5(u[5]+u[6]))^p[1]/((0.5(u[5]+u[6]))^p[1]+(p[8]p[2])^p[1])(1-u[3]^p[10]/(u[3]^p[10]+(p[9]p[11])^p[10])) function affect1!(integrator) integrator.u[1] += 1 end r1 = ConstantRateJump(rate1,affect1!)

rate2(u,p,t) = p[7](0.5(u[5]+u[6]))^p[1]/((0.5(u[5]+u[6]))^p[1]+(p[8]p[2])^p[1])(1-u[4]^p[10]/(u[4]^p[10]+(p[9]p[11])^p[10])) function affect2!(integrator) integrator.u[2] += 1 end r2= ConstantRateJump(rate2,affect2!)

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

rate4(u,p,t) = p[9](u[2]^p[5]/(u[2]^p[5]+(p[7]p[6])^p[5])) function affect4!(integrator) integrator.u[4] += 1 end r4= ConstantRateJump(rate4,affect4!)

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

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

rate7(u,p,t) = 1u[1] function affect7!(integrator) integrator.u[1] -= 1 end r7= ConstantRateJump(rate7,affect7!) rate8(u,p,t) = 1u[2] function affect8!(integrator) integrator.u[2] -= 1 end r8= ConstantRateJump(rate8,affect8!) rate9(u,p,t) = 1u[3] function affect9!(integrator) integrator.u[3] -= 1 end r9= ConstantRateJump(rate9,affect9!) rate10(u,p,t) = 1u[4] function affect10!(integrator) integrator.u[4] -= 1 end r10= ConstantRateJump(rate10,affect10!) rate11(u,p,t) = 1u[5] function affect11!(integrator) integrator.u[5] -= 1 end r11= ConstantRateJump(rate11,affect11!) rate12(u,p,t) = 1u[6] function affect12!(integrator) integrator.u[6] -= 1 end r12= ConstantRateJump(rate12,affect12!)

p = [2.8177, 0.134, 2.71, 0.130, 3.7, 0.41, 1000, 100, 100, 2.82,0.879, 3.44, 0.29] x0 = [999,1,1,99,1,99] tspan = (0.0,100.0)

prob = DiscreteProblem(x0,tspan,p); jump_prob = JumpProblem(prob,Direct(),JumpSet(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12)); sol = solve(jump_prob,FunctionMap(),dense=true);

plot(sol.t,sol.u)

Am 08.02.2018 um 15:42 schrieb Christopher Rackauckas notifications@github.com:

Hey, If you are defining the rates directly, then the rate equations can use parameters. For example:

rate(u,p,t) = (0.1/1000.0)u[1]u[2] function affect!(integrator) integrator.u[1] -= 1 integrator.u[2] += 1 end jump = ConstantRateJump(rate,affect!) is what's in the documentation, but you can make use of p in this equation. (Ref: http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html#Defining-the-Jumps-Directly-1 http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html#Defining-the-Jumps-Directly-1 )

Note that currently the reaction network DSL is unable to handle parameters, but we are in the process of fixing that, and hopefully that will get released in the next few weeks. When that's released, we will update that tutorial to show it using parameters like the others.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/246#issuecomment-364132163, or mute the thread https://github.com/notifications/unsubscribe-auth/Aigwuh2WvKKC0M1_qMS3WmDgS-FG2sC7ks5tSwffgaJpZM4R-YX8.

---

Edda Schulz Max Planck Research Group "Regulatory Networks in Stem Cells" Max Planck Institute for Molecular Genetics Ihnestraße 63-73 14195 Berlin Tel: +49-30-8413-1226 E-mail : Edda.Schulz@molgen.mpg.de mailto:Edda.Schulz@molgen.mpg.de Web: www.molgen.mpg.de/molgen/en/Forschung/Otto-Warburg-Laboratorium/networks_stem_cells http://www.molgen.mpg.de/molgen/en/Forschung/Otto-Warburg-Laboratorium/networks_stem_cells

[onoderat]

I suspect that it may be an issue with using the Plots.jl package. I once had to plot some data with over a million points and Plots.jl chocked on me. I would advice using PyPlot for large data sizes.

[ChrisRackauckas]

PyPlot is slow for large data. The GR backend handles it much better, or InspectDR. But the problem here may just be the defaults. The default number of points is here:

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/solutions/solution_interface.jl#L57

so for a discrete problem it will be:

max(1000,100*length(sol))

For a large number of points, that means it's plotting 100*length(sol) points. It was set high so that way it resolves the discontinuities well, but my guess is that needs to have a ceiling since it's just growing really huge. What happens if you set plotdensity to something like 100000?

[ChrisRackauckas]

I see what the issue is.

plot(sol.t,sol.u)

That's nonsensical since you're telling it to plot an array of arrays (so it plots 130,000 timeseries). The plot recipe works fine:

plot(sol,plotdensity=1000)

and we should cap the plotdensity here though. If you want to specifically plot the points, it notes in the documentation

http://docs.juliadiffeq.org/latest/basics/plot.html#Plotting-Without-the-Plot-Recipe-1

that you'd want to transpose if you want the timeseries, and

plot(sol.t,sol')

works fine.

[ChrisRackauckas]

The new reaction DSL can handle parameters:

http://docs.juliadiffeq.org/latest/models/biological.html

Let us know if you need anything else.