Non-diagonal and commutative SDE benchmarks

[ChrisRackauckas]

@onoderat has moved us forward finally in the Non-diagonal SDE department, and so we should think about creating some benchmark notebooks in this area.

[onoderat]

Agreed. Where are the ipython notebook stored in the repo?

[ChrisRackauckas]

https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl/tree/master/NonStiffSDE and https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl/tree/master/StiffSDE . README is just links to the nbviewer versions.

[onoderat]

These notebook are very nice. I would like to add the PCEuler to the following notebooks,

1. QuadraticStiffness.ipynb
2. BasicSDEWorkPrecision.ipynb

I probably wont have any time to work on this for the next 2-3 weeks. Would it be ok for me to finish this around the beginning of March?

On a slightly unrelated note, I personally think that PCEuler should not be released until the documentation plus the benchmarks have been finished - just to be extra careful. What do you think?

[onoderat]

I forgot to add that I can work on the non-diagonal and commutative test as well. That should more or less be a copy and paste of the testing code in nondiagonal_test.jl .

The first thing to do would probably be to refactor the problem code into the DiffEqProblemLibrary.

[ChrisRackauckas]

On a slightly unrelated note, I personally think that PCEuler should not be released until the documentation plus the benchmarks have been finished - just to be extra careful. What do you think?

If it's passing tests I think it's safe to release. People won't actually make use of it until it's in the documentation so that is some extra safety.

The first thing to do would probably be to refactor the problem code into the DiffEqProblemLibrary.

That would be heroic 👍 . I always think about doing it but never get around to it...

[ChrisRackauckas]

@Gaussia do you have some problems which could make good benchmarks? They don't need analytical solutions, though if you know of any that's even better. Other noise types are welcome too of course.

[TorkelE]

Not sure exactly what you are looking for, but here are some ideas:

f(du,u,p,t) = du[1] = 5 + 18*u[1]^4/(u[1]^4+11.9^4) - u[1]
function g(du,u,p,t)
  du[1,1] = .1*sqrt(5 + 18*u[1]^4/(u[1]^4+11.9^4))
  du[1,2] = -.1*sqrt(u[1])
end
prob = SDEProblem(f,g,[3.],(0.,300.),noise_rate_prototype=zeros(1,2))
sol = solve(prob,ISSEM())

This will generate stuff like (several solutions plotted) this:

Firs there is a small jump to an initial equilibrium. At some random point in time there will be a sudden jump to a second "more stable" equilibrium were the solution will noise around for eternity.

The Brusselator is always fun...

function f(du,u,p,t)
    du[1] = 1 + u[1]*u[1]*u[2] - 2.5*u[1]-u[1]
    du[2] = 2.5*u[1] - u[1]*u[1]*u[2]
end
function g(du,u,p,t)
  du[1,1] = 0.15*sqrt(1)
  du[1,2] = 0.15*sqrt(u[1]*u[1]*u[2])
  du[1,3] = -0.15*sqrt(2.5*u[1])
  du[1,4] = -0.15*sqrt(u[1])
  du[2,1] = 0
  du[2,2] = -0.15*sqrt(u[1]*u[1]*u[2])
  du[2,3] = 0.15*sqrt(2.5*u[1])
  du[2,4] = 0
end
prob = SDEProblem(f,g,[3.,2.],(0.,100.),noise_rate_prototype=zeros(2,4))
 sol = solve(prob,ISSEM())

resulting in something like:

If you want something big and messy you can always try

network = @reaction_network rnType  begin
    0.01, (X,Y,Z) --> 0
    hill(X,3.,100.,-4), 0 --> Y
    hill(Y,3.,100.,-4), 0 --> Z
    hill(Z,4.5,100.,-4), 0 --> X
    hill(X,2.,100.,6), 0 --> R
    hill(Y,15.,100.,4)*0.002, R --> 0
    20, 0 --> S
    R*0.005, S --> SP
    0.01, SP + SP --> SP2
    0.05, SP2 --> 0
end;
prob = SDEProblem(network, [200.,60.,120.,100.,50.,50.,50.], (0.,4000.))
sol = solve(prob,ISSEM())

the solution to this will look something like

[ChrisRackauckas]

Those look great! Any background information on them? I know the Brusselator but not the others. Usually for DiffEqProblemLibrary we make sure to cite the source model (and then in papers of course the citation is necessary)

[ChrisRackauckas]

For future reference, here's some stochastic Brusselators and they utilize a much simpler noise term probably to make it easier to compute:

https://www.sciencedirect.com/science/article/pii/S0377042704000731
https://arxiv.org/abs/1405.0390

[TorkelE]

The first one is just a very simple bistable system (X activating its own production and getting degraded). I guess people come across it quite often but I don't think there is any special resources on it. The last one is just a random system I made myself while experimenting with biochemical modelling, it does not have any deeper meaning.

[ChrisRackauckas]

Got it, thanks!

[ChrisRackauckas]

@Gaussia shouldn't the reaction rates in the noise term not be square rooted?

https://aip.scitation.org/doi/pdf/10.1063/1.481811 (Equation 23)

[TorkelE]

You mean that it should rather be

function g(du,u,p,t)
...
  du[1,3] = -0.15*2.5*sqrt(u[1])
...
  du[2,3] = 0.15*2.5*sqrt(u[1])
...
end

I might be wrong but I think that is not the case? According to Gillespie 2000, Equation 23, the noise term is
v_i,j * sqrt(a_i,j(X))
Here a_i,j should be the probability of a certain reaction happening (including both reaction rates and concentrations of reactants, Equation 4) and v_i,j the part not being square rooted is only the stoichiometry (Equation 3).

If I am wrong then I should probably have this cleared up asap, but I am fairly confident that this is the case.

(Also: If it is possible to use Markdown or something to make e.g. v_i,j look pretty I would be happy to know, have searched around for it but found nothing)

[ChrisRackauckas]

v_i,j the part not being square rooted is only the stoichiometry (Equation 3).

Oh yeah, that is stoich, in which case you're correct. I thought it was reaction rate when re-reading.