SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
139 stars 35 forks source link

Tracking issue for Documenter and design of docstrings and new tutorials. #238

Closed 00krishna closed 2 weeks ago

00krishna commented 2 years ago

I am creating a new issue to help track the work on migrating some of the DiffEqJump.jl documentation to Documenter.jl. In talking with @isaacsas , he indicated that he wanted to revamp the introductory tutorials for DiffEqJump.

Originally I thought I would just migrate the tutorials to Documenter and incorporate the docstrings into the tutorials. But instead of that, perhaps it is better to review how others have implemented the docstrings for some other larger packages in the combined SciML docs.

This issue can just help us to coordinate efforts and agree on a design.

isaacsas commented 2 years ago

Thanks! Yes, my plan, which hopefully I can get to this weekend or next week was to split the big tutorial into several pieces and clean them up. Something along the lines of pure jump models with ConstantRate and MassActionJumps (perhaps split into two pieces itself), followed by a tutorial involving VariableRateJumps.

The FAQ would be split out of the tutorial to a separate section. Basically I'm envisioning something like the Catalyst doc structure: https://catalyst.sciml.ai

We could also add and cleanup the jump diffusion tutorial -- I'd leave that decision to @ChrisRackauckas on where he feels that should live. (It would also make sense in StochasticDiffEq.)

Then we need to figure out how to reconcile the current JumpProblem page, and to a lesser extent the jump solve page, with the API doc strings to make an API page. They don't completely contain the same information, so it may be the doc strings need to be cleaned up and updated so the API docs are equivalent.

isaacsas commented 2 years ago

I don't think docstrings should be in the tutorials -- see the Catalyst style. I'd prefer the tutorials to introduce how to use DiffEqJump in the most common use cases, and then leave the docstrings for more detailed info about the specific components.

isaacsas commented 2 years ago

@00krishna if you want to take a first pass the simplest thing to do would be to initially split the tutorial in two parts, and have the second start with a note that one should be familiar with the first before looking at it. We could also convert it to use @example blocks so the code runs dynamically. Then put the FAQ on its own page/section, and just migrate the JumpProblem and jump solving pages (perhaps under a "Reference and Guidance" section). Finally, there can be an API section with docstrings. Initially that should probably be something like

JumpProblem
SSAStepper

# Jump Types
ConstantRateJump
VariableRateJump
MassActionJump
RegularJump

# Aggregators (see src/aggregators/aggregators.jl for a list)
Direct
FRM
NRM
...

That would be a good starting point to then build out from I think.

00krishna commented 2 years ago

@isaacsas Yeah, I can run with this. I will look at the Catalyst.jl docs as a model then, and try to emulate. I will post my design here so that you can make sure I am on the right track @isaacsas .

isaacsas commented 2 years ago

TODO still:

00krishna commented 2 years ago

@isaacsas say, I was looking through the tutorials as a basis for the docstrings. Seems like a lot of the tutorial code is focused on using Catalyst.jl and those examples--the main one being the SIR model. For the docstrings, should I put so much emphasis on Catalyst, or should I setup the problems as regular DifferentialEquations.jl problems, like ODEProblem or I believe DiscreteProblem would also work in some cases.

isaacsas commented 2 years ago

I'm not sure what you mean? The tutorial starts with Catalyst, but everything after it is pure DiffEqJump. Also, most of the core functions have docstrings showing now at https://jump.sciml.ai/dev/api/

00krishna commented 2 years ago

Ahh, I see. I was looking at the tutorial on https://jump.sciml.ai/dev/tutorials/discrete_stochastic_example/ . Right, so the original model is Catalyst for the SIR example, but then subsequent examples are regular DiffEq. Yeah, I think the DiffEq examples are one line functions, like rate1(u,p,t) = p[1]*u[1]*u[2], so I just read over them too quickly.

Thanks for setting me straight.

isaacsas commented 2 years ago

But maybe I should update that tutorial to move the Catalyst stuff to the end. Then it can start with DiffEqJump but present Catalyst at the end for people that are interested only in reactions systems and want an easier interface.

00krishna commented 2 years ago

I think that the SIR model example is a good one. You could define the SIR model in both DiffEq, as well as in Catayst, and present them together? I think that threw me off was that the tutorials starts with a detailed description of the SIR model in Catalyst form. So once that expectation was set, I did not notice the DiffEq examples below. So using an extended example is good, but perhaps showing the SIR example in both Catalyst and DiffEq forms? Like can I define a standard Diffeq model f(du, u, p, t) and then use that in the DiscreteProblem just like you used the Catalyst model setup for sir_model. I will try it out, but that was my thinking on what I want to check.

isaacsas commented 2 years ago

Yeah, I think it makes sense to reverse the order. i.e. show SIR in DiffEqJump first, and then at the end point out it is easier to do these types of chemical systems with Catalyst (which has some other advantages). That way the information is there, but people who just want to learn DiffEqJump and don't care about chemical systems won't get hung up on it. I'm working on such an update now. Thanks for the suggestion!

00krishna commented 2 years ago

Say, @isaacsas I had an idea to adjust the formatting of the code blocks in the docstrings, but I wanted to check with you first. I thought an adjustment might make the DiffEqJump code blocks a bit more consistent with the DifferentialEquations codeblocks.

The current code blocks are formatted like this:

using DiffEqJump
β = 0.1 / 1000.0; ν = .01;
p = (β,ν)
rate1(u,p,t) = p[1]*u[1]*u[2]  # β*S*I
function affect1!(integrator)
  integrator.u[1] -= 1         # S -> S - 1
  integrator.u[2] += 1         # I -> I + 1
end
jump = ConstantRateJump(rate1,affect1!)

rate2(u,p,t) = p[2]*u[2]      # ν*I
function affect2!(integrator)
  integrator.u[2] -= 1        # I -> I - 1
  integrator.u[3] += 1        # R -> R + 1
end
jump2 = ConstantRateJump(rate2,affect2!)
u₀    = [999,1,0]
tspan = (0.0,250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump, jump2)
sol = solve(jump_prob, SSAStepper())

But most of the original DifferentialEquations.jl codeblocks have the function first and the parameters afterwards. This adjustment also gives a bit more prominence to the model function itself, so that a reader does not miss it. I was thinking to reformat the block like this.

using DiffEqJump

function rate1(u,p,t) 
    p[1]*u[1]*u[2]  # β*S*I
end

function affect1!(integrator)
  integrator.u[1] -= 1         # S -> S - 1
  integrator.u[2] += 1         # I -> I + 1
end

jump = ConstantRateJump(rate1,affect1!)

function rate2(u,p,t)
    p[2]*u[2]      # ν*I
end

function affect2!(integrator)
  integrator.u[2] -= 1        # I -> I - 1
  integrator.u[3] += 1        # R -> R + 1
end

jump2 = ConstantRateJump(rate2,affect2!)

β = 0.1 / 1000.0; ν = .01;
p = (β,ν)
u₀    = [999,1,0]
tspan = (0.0,250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump, jump2)
sol = solve(jump_prob, SSAStepper())

Let me know either way, and I will use that as a model. Like I said, it is just an idea.

isaacsas commented 2 years ago

I don't think putting the parameters after is really needed, and maybe it is a bit confusing (since the earlier declaration makes clear the parameters are beta and nu, which are then referenced in subsequent comments -- moving it after means the parameter names are referenced before they are actually defined).

An alternative would be to move the p declaration after, but make p a named tuple and then access it via the named fields.

00krishna commented 2 years ago

Okay. So then something like this below? It adds a bit of space and makes sure that the rate and affect functions themselves are visible on quick glance.


β = 0.1 / 1000.0; ν = .01;
p = (β,ν)

function rate1(u,p,t) 
    p[1]*u[1]*u[2]  # β*S*I
end

function affect1!(integrator)
  integrator.u[1] -= 1         # S -> S - 1
  integrator.u[2] += 1         # I -> I + 1
end

jump = ConstantRateJump(rate1,affect1!)

function rate2(u,p,t)
    p[2]*u[2]      # ν*I
end

function affect2!(integrator)
  integrator.u[2] -= 1        # I -> I - 1
  integrator.u[3] += 1        # R -> R + 1
end

jump2 = ConstantRateJump(rate2,affect2!)

u₀    = [999,1,0]
tspan = (0.0,250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump, jump2)
sol = solve(jump_prob, SSAStepper())
isaacsas commented 2 years ago

Sure

00krishna commented 2 years ago

@isaacsas I created the first set of doctests for the jump types. One thing to note is that the doctests are building, but they are indicating a failure a warning. I created the docstring for RegularJumps and included the code below. The warning I am getting is about

The RegularJump interface has changed to be matrix-free. See the documentation for more details.
│    └ @ DiffEqJump /media/hayagriva/Dropbox/sandbox/DiffEqJump.jl/src/jumps.jl:181

I used the code below.

using DiffEqJump

function rate(out,u,p,t)
    out[1] = (0.1/1000.0)*u[1]*u[2]
    out[2] = 0.01u[2]
end

function c(dc,u,p,t,mark)
    dc[1,1] = -1
    dc[2,1] = 1
    dc[2,2] = -1
    dc[3,2] = 1
end

dc = zeros(3,2)
rj = RegularJump(rate,c,dc;constant_c=true)

So the problem seems to be that dc is an AbstractMatrix and that is dispatching to the function with the warning. But I think it is okay to publish this code, because it will just produce a warning, but the code will otherwise work. If you tell me how to adjust this code to avoid the warning, then I can incorporate that as well.

isaacsas commented 2 years ago

Someone changed the interface a while back but I don't think it has gotten much use since the change. You can look at the regular jump tests

https://github.com/SciML/DiffEqJump.jl/blob/master/test/regular_jumps.jl

and the definitions to see how it works now (notice it takes a counts argument now too)

https://github.com/SciML/DiffEqJump.jl/blob/277cb15fcd3be6cf80b6a3f1fc9fcd69219af831/src/jumps.jl#L94-L118

It is also used in StochasticDiffEq at

https://github.com/SciML/StochasticDiffEq.jl/blob/master/test/tau_leaping.jl

isaacsas commented 2 years ago

Getting a doc string on RegularJumps with an example would be good.

00krishna commented 2 years ago

@isaacsas Okay good. I got a working example--no warnings. Note that in order to include this example, I will need to include StableRNGs, and LinearAlgebra in the Project.toml. Technically I could probably get away without StableRNGs, but keeping it seems helpful if people want to make replicable examples.

I will go ahead and put it in the PR with this code.

using DiffEqJump
using StableRNGs
using LinearAlgebra
rng = StableRNG(12345)

function regular_rate(out, u, p, t)
  out[1] = (0.1 / 1000.0) * u[1] * u[2]
  out[2] = 0.01u[2]
end

function regular_c(du, u, p, t, counts, mark)
    mul!(du, dc, counts)
end

dc = zeros(3, 2)
rj = RegularJump(regular_rate, regular_c, 2)
jumps = JumpSet(rj)
prob = DiscreteProblem([999, 1, 0], (0.0, 250.0))
jump_prob = JumpProblem(prob, Direct(), rj; rng = rng)
sol = solve(jump_prob, SimpleTauLeaping(); dt = 1.0)
00krishna commented 2 years ago

@isaacsas i am traveling the yesterday and today, but will get back to the PR in the next day or so. Just wanted to let you know.

isaacsas commented 2 years ago

No worries. Have a good trip!

00krishna commented 2 years ago

@isaacsas Gradually getting back on the horse :). So I did more research and thinking and it seems like doctests in the docstrings only support test cases where some result is returned and compared to a desired result. Those @example blocks only work in markdown documents, and are not run when pulling docstrings from source files.

So we could just go with julia blocks for now, and then add tests to the tutorials. There are a couple of ways to do this as well. We could either see if the existing tutorials can be setup with @example blocks that are sufficient tests for CI purposes. Or, we can also create a set of hidden documents where we can write more comprehensive tests that are run at build time, but that users never see. This hidden option works if say the tutorials have code blocks that are not full tests or code runs, etc.

I will read over the new tutorials to figure out the options.

isaacsas commented 2 years ago

I think julia blocks for docstrings with simple examples aimed at helping users see how to use the functions are fine. I don't really see the point of silent tests via the documentation -- at that point we might as well add the associated test as a normal package test instead. We have testing pretty well covered with the normal tests, and I've setup the new documentation to run all the tutorials as example blocks.

isaacsas commented 2 weeks ago

Closing as I think this is out of date now. Happy to start doc discussions in a new thread.