Merge 02dac52 into 277cb15

.github/workflows/Documentation.yml

@@ -21,4 +21,5 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+          GKSwstype: "100" # https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
         run: julia --project=docs/ docs/make.jl

Project.toml

@@ -17,6 +17,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TreeViews = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
@@ -32,6 +33,7 @@ PoissonRandom = "0.4"
 RandomNumbers = "1.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1.0"
+SciMLBase = "1.35.1"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 TreeViews = "0.3"
 UnPack = "1.0.2"

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -2,13 +2,26 @@ using Documenter, DiffEqJump
 
 include("pages.jl")
 
+mathengine = MathJax3(Dict(:loader => Dict("load" => ["[tex]/require", "[tex]/mathtools"]),
+                           :tex => Dict("inlineMath" => [["\$", "\$"], ["\\(", "\\)"]],
+                                        "packages" => [
+                                            "base",
+                                            "ams",
+                                            "autoload",
+                                            "mathtools",
+                                            "require",
+                                        ])))
+
 makedocs(sitename = "DiffEqJump.jl",
          authors = "Chris Rackauckas",
          modules = [DiffEqJump],
-         clean = true, doctest = false,
-         format = Documenter.HTML(analytics = "UA-90474609-3",
+         clean = true,
+         doctest = false,
+         format = Documenter.HTML(; analytics = "UA-90474609-3",
                                   assets = ["assets/favicon.ico"],
-                                  canonical = "https://diffeqjump.sciml.ai/stable/"),
+                                  canonical = "https://jump.sciml.ai/stable/",
+                                  prettyurls = (get(ENV, "CI", nothing) == "true"),
+                                  mathengine),
          pages = pages)
 
 deploydocs(repo = "github.com/SciML/DiffEqJump.jl.git";

docs/old/cut.md

@@ -0,0 +1,64 @@
+##  Coupling Jump Problems
+THIS NEEDS UPDATING
+
+In many applications one is interested in coupling two stochastic processes.
+This has applications in Monte Carlo simulations and sensitivity analysis, for
+example. Currently, the coupling that is implemented for jump processes is known
+as the split coupling. The split coupling couples two jump processes by coupling
+the underlying Poisson processes driving the jump components.
+
+Suppose `prob` and `prob_control` are two problems we wish to couple. Then the
+coupled problem is obtained by
+```@example tut3
+prob_coupled =  SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map)
+```
+Here, `coupling_map` specifies which jumps to couple. If `(j,i)` is in
+`coupling_map`, then the `i`th jump in `prob` will be coupled to the `j`th jump
+in `prob_control`. Note that currently `SplitCoupledJumpProblem` is only
+implemented for `ConstantRateJump` problems.
+
+As an example, consider a doubly stochastic Poisson process, that is, a Poisson
+process whose rate is itself a stochastic process. In particular, we will take
+the rate to randomly switch between zero and `10` at unit rates:
+```@example tut3
+rate(u, p, t) = 10 * u[2]
+affect!(integrator) = integrator.u[1] += 1
+jump1 = ConstantRateJump(rate, affect!)
+
+rate(u, p, t) = u[2]
+affect!(integrator) = (integrator.u[2] -= 1; integrator.u[3] += 1)
+jump2 = ConstantRateJump(rate, affect!)
+
+rate(u,p,t) = u[3]
+affect!(integrator) = (integrator.u[2] += 1; integrator.u[3] -= 1)
+jump3 = ConstantRateJump(rate, affect!)
+
+prob = DiscreteProblem(u0, tspan)
+jump_prob = JumpProblem(prob, Direct(), jump1, jump2, jump3)
+```
+The doubly stochastic Poisson process has two sources of randomness: one due to
+the Poisson process, and another due to random evolution of the rate. This is
+typical of many multiscale stochastic processes appearing in applications, and
+it is often useful to compare such a process to one obtained by removing one
+source of randomness. In present context, this means looking at an ODE with
+constant jump rates, where the deterministic evolution between jumps is given by
+the expected value of the Poisson process:
+```@example tut3
+function f(du, u, p, t)
+  du[1] = 10 * u[2]
+  du[2] = 0
+  du[3] = 0
+end
+prob_control = ODEProblem(f, [0.0, 0.0, 1.0], tspan)
+jump_prob_control = JumpProblem(prob_control, Direct(), jump2, jump3)
+```
+Let's couple the two problems by coupling the jumps corresponding to the
+switching of the rate:
+```@example tut3
+coupling_map = [(2,1), (3,2)]
+prob_coupled =  SplitCoupledJumpProblem(jump_prob, jump_prob_control, Direct(), coupling_map)
+```
+Now `prob_coupled` can be solved like any other `JumpProblem`:
+```@example tut3
+sol = solve(prob_coupled, Tsit5())
+```

docs/pages.jl

@@ -2,7 +2,8 @@
 
 pages = [
     "Home" => "index.md",
-    "Tutorials" => Any["tutorials/discrete_stochastic_example.md",
+    "Tutorials" => Any["tutorials/simple_poisson_process.md",
+                       "tutorials/discrete_stochastic_example.md",
                        "tutorials/jump_diffusion.md"],
     "FAQ" => "faq.md",
     "Type Documentation" => Any["Jump types and JumpProblem" => "jump_types.md",

docs/src/api.md

@@ -11,10 +11,10 @@ SSAStepper
 
 ## Types of Jumps
 ```@docs
-MassActionJump
 ConstantRateJump
+MassActionJump
 VariableRateJump
-RegularJump
+JumpSet
 ```
 
 ## Aggregators

docs/src/faq.md

@@ -19,42 +19,6 @@ condition, and other components of a generated `JumpProblem`. This can be useful
 when trying to call `solve` many times while avoiding reallocations of the
 internal aggregators for each new parameter value or initial condition.
 
-## How do I use callbacks with `ConstantRateJump` or `MassActionJump` systems?
-
-Callbacks can be used with `ConstantRateJump`s and `MassActionJump`s. When
-solving a pure jump system with `SSAStepper`, only discrete callbacks can be
-used (otherwise a different time stepper is needed).
-
-*Note, when modifying `u` or `p` within a callback, you must call
-`reset_aggregated_jumps!(integrator)` after making updates.* This ensures that
-the underlying jump simulation algorithms know to reinitialize their internal
-data structures. Leaving out this call will lead to incorrect behavior!
-
-A simple example that uses a `MassActionJump` and changes the parameters at a
-specified time in the simulation using a `DiscreteCallback` is
-```julia
-using DiffEqJump
-rs = [[1 => 1], [2 => 1]]
-ns = [[1 => -1, 2 => 1], [1 => 1, 2 => -1]]
-p = [1.0, 0.0]
-maj = MassActionJump(rs, ns; param_idxs=[1, 2])
-u₀ = [100, 0]
-tspan = (0.0, 40.0)
-dprob = DiscreteProblem(u₀, tspan, p)
-jprob = JumpProblem(dprob, Direct(), maj)
-pcondit(u, t, integrator) = t == 20.0
-function paffect!(integrator)
-    integrator.p[1] = 0.0
-    integrator.p[2] = 1.0
-    reset_aggregated_jumps!(integrator)
-end
-sol = solve(jprob, SSAStepper(), tstops=[20.0], callback=DiscreteCallback(pcondit, paffect!))
-```
-Here at time `20.0` we turn off production of `u[2]` while activating production
-of `u[1]`, giving
-
-![callback_gillespie](assets/callback_gillespie.png)
-
 ## How can I define collections of many different jumps and pass them to `JumpProblem`?
 
 We can use `JumpSet`s to collect jumps together, and then pass them into
@@ -90,12 +54,14 @@ argument. Continuing the previous example:
 ```julia
 #] add RandomNumbers
 using RandomNumbers
-jprob = JumpProblem(dprob, Direct(), maj, rng=Xorshifts.Xoroshiro128Star(rand(UInt64)))
+jprob = JumpProblem(dprob, Direct(), maj,
+                    rng = Xorshifts.Xoroshiro128Star(rand(UInt64)))
 ```
 uses the `Xoroshiro128Star` generator from
 [RandomNumbers.jl](https://github.com/JuliaRandom/RandomNumbers.jl).
 
-On version 1.7 and up DiffEqJump uses Julia's builtin random number generator by default. On versions below 1.7 it uses `Xoroshiro128Star`.
+On version 1.7 and up DiffEqJump uses Julia's builtin random number generator by
+default. On versions below 1.7 it uses `Xoroshiro128Star`.
 
 ## What are these aggregators and aggregations in DiffEqJump?
 
@@ -134,4 +100,78 @@ The four jumps would be ordered by the first jump in `maj`, the second jump in
 then follow this ordering when assigning an integer id to each jump.
 
 See also [Constant Rate Jump Aggregators Requiring Dependency Graphs](@ref) for
-more on dependency graphs needed for the various SSAs.
+more on dependency graphs needed for the various SSAs.
+
+## How do I use callbacks with `ConstantRateJump` or `MassActionJump` systems?
+
+Callbacks can be used with `ConstantRateJump`s and `MassActionJump`s. When
+solving a pure jump system with `SSAStepper`, only discrete callbacks can be
+used (otherwise a different time stepper is needed). When using an ODE or SDE
+time stepper any callback should work.
+
+*Note, when modifying `u` or `p` within a callback, you must call
+`reset_aggregated_jumps!(integrator)` after making updates.* This ensures that
+the underlying jump simulation algorithms know to reinitialize their internal
+data structures. Leaving out this call will lead to incorrect behavior!
+
+A simple example that uses a `MassActionJump` and changes the parameters at a
+specified time in the simulation using a `DiscreteCallback` is
+```julia
+using DiffEqJump
+rs = [[1 => 1], [2 => 1]]
+ns = [[1 => -1, 2 => 1], [1 => 1, 2 => -1]]
+p = [1.0, 0.0]
+maj = MassActionJump(rs, ns; param_idxs=[1, 2])
+u₀ = [100, 0]
+tspan = (0.0, 40.0)
+dprob = DiscreteProblem(u₀, tspan, p)
+jprob = JumpProblem(dprob, Direct(), maj)
+pcondit(u, t, integrator) = t == 20.0
+function paffect!(integrator)
+    integrator.p[1] = 0.0
+    integrator.p[2] = 1.0
+    reset_aggregated_jumps!(integrator)
+    nothing
+end
+sol = solve(jprob, SSAStepper(), tstops=[20.0], callback=DiscreteCallback(pcondit, paffect!))
+```
+Here at time `20.0` we turn off production of `u[2]` while activating production
+of `u[1]`, giving
+
+![callback_gillespie](assets/callback_gillespie.png)
+
+
+## How can I access earlier solution values in callbacks?
+When using an ODE or SDE time-stepper that conforms to the [integrator
+interface](https://docs.sciml.ai/dev/modules/DiffEqDocs/basics/integrator/) one
+can simply use `integrator.uprev`. For efficiency reasons, the pure jump
+[`SSAStepper`](@ref) integrator does not have such a field. If one needs
+solution components at earlier times one can save them within the callback
+condition by making a functor:
+```julia
+# stores the previous value of u[2] and represents the callback functions
+mutable struct UprevCondition{T}
+     u2::T
+end
+
+# condition
+function (upc::UprevCondition)(u, t, integrator)
+    # condition for the callback is that the new value of u[2]
+    # is smaller than the previous value
+    condit = u[2] - upc.u2 < 0
+
+    # save the new value as the previous value
+    upc.u2 = u[2]
+
+    condit
+end
+
+# affect!
+function (upc::UprevCondition)(integrator)
+    integrator.u[4] -= 1
+    nothing
+end
+
+upc = UprevCondition(u0[2])
+cb = DiscreteCallback(upc, upc)
+```