docs: update README and examples (#8)

Author: wangl-cc

README.md

@@ -29,68 +29,60 @@ the state of host population can also be represented as a vector
 by concatenating the state of each component of host `u = vcat(S, vec(I), vec(R))`,
 where `S` is a vector of length `n` and `I`, `R` are matrixes of size `n × m`.
 However, in evolutionary biology,
-the "mutation" and "extinction" will change the types of "hosts" and "viruses",
+the "mutation" and "extinction" will change the types of hosts and viruses,
 which means the `n` and `m` changes during the evolution,
 and the length of the state vector `u` will also change.
 
 ## How to use?
 
-Similarly to `DifferentialEquations.jl`,
-you must define reactions firstly (jump in `DifferentialEquations.jl`),
-For example, the infect reaction of above SIR model is defined as:
+The package [`Catalyst.jl`](https://github.com/SciML/Catalyst.jl) provides a simple way
+to build biochemical reaction for `DifferentialEquations.jl`.
+Similarly, this package provides a macro `@reaction_eq` which generate reaction(s),
+with given equation.
 
+For example, the infect reaction of above SIR model with multi-type of viruses can be defined as:
 ```julia
-# reaction: S_i + I_ij -> 2I_ij
-@cfunc infect_c(β, S, I) = β * S' * I # function to calculate the rate of infection
-@ufunc function infect_u!(ind, S, I) # function to update the state when infection happens
-    # ind is the index of the reaction selected with weight calculated by `infect_c`
-    S[ind[2]] -= 1 # S_j -= 1
-    I[ind] += 1    # I_ij += 1
-    return nothing
-end
-infect = Reaction(infect_c, infect_u!)
+@reaction_eq infect β S + I[i] --> 2I[i]
 ```
+where `i` donates the virus type.  The equation `S + I[i] --> 2I[i]` means
+the an host infected by the virus `i` infect a susceptible host with rate `β`,
+then convert the the susceptible host to a infectious host.
+This expression will not only generate one reaction but a group of reactions trough the index `i`.
 
-where `@cfunc` is a macro to help you define a function to calculate the rate of reaction,
-`@ufunc` is a macro to help you define a function to update the state when reaction happens.
-Unlike `DifferentialEquations.jl`, the "calculate" function returns an array of rates,
-which allows variable length state and contains a lot sub-reactions.
-Additionally, the "update" function accepts a special argument `ind`,
-which tells you which sub-reaction is selected and should always be the first argument.
-All of these is designed to work with variable length state.
-
-Besides, there is another macro `@reaction` helps you define a reaction more easily.
-The below code is equivalent to the above code:
-
+However, the mutation and extinction can not be defined easily by the macro `@reaction_eq` currently.
+Thus the alternative macro `@reaction` provides a low-level way to build reaction(s)
 ```julia
-@reaction infect begin
-    β * S' * I
+@reaction mutation begin
+    @quickloop μ * I[i]
     begin
-        S[ind[2]] -= 1
-        I[ind] += 1
+        i = ind[1]
+        I[i] -= 1 # the host individual is converted to another type
+        push!(I, 1) # add a new type of infectious host to the system
+        push!(α, randn() + α[i]) # the virulence of the new host type is generated randomly with mean `α[i]`
     end
 end
 ```
+This expression defined mutation of virus which contains two parts:
+1. The `@quickloop μ * I[i]` defines how to calculate the rate of mutation,
+2. The `begin ... end` block defines what happens when the host is mutated,
+   where `ind` is a preserved variable which is used to store the index of the mutated host.
 
 Once you have defined all reactions, put them together as a tuple:
-
 ```julia
-reactions = (infect, ...)
+reactions = (infect, mutation, ...)
 ```
-
 and define the initial state and parameters of the system as a named tuple:
-
 ```julia
-params = (β = 0.1, ..., S = [100], I = [1;;], R = [0;;]) # the [1;;] syntax requires Julia >= 1.7
+params = (β = 0.1, μ = 0.001, ..., S = scalar(100), I = [1], R = [0], α = [0.5])
 ```
+where the `scalar` will create a type similar to `Number`,
+but it can be update in-place like `Ref` like `S[] += 1`.
 
-Then, you can use `gillespie` to simulate the system:
-
+Once reaction and parameters is defined, you can use `gillespie` to simulate the system:
 ```julia
 max_time = 100 # the maximum time of simulation
 params′, t, term = gillespie(max_time, params, reactions)
 ```
-
 where the `gillespie` function returns an tuple `t, ps′, term`
 where `t` is the time when the simulation ends, `params′` is an updated `params`,
 and `term` is a flag indicating whether the simulation is finished after the maximum time
@@ -102,15 +94,13 @@ but you can use my another package `RecordedArrays` to record them, like this:
 ```julia
 using RecordedArrays
 c = DiscreteClock(max_time) # clock store information of max_time
-S = recorded(DynamicEntry, c, [100]) # create a recorded vector as S with the clock c
-I = recorded(DynamicEntry, c, [1;;]) # create a recorded matrix as I with the clock c
-R = recorded(DynamicEntry, c, [0;;]) # create a recorded matrix as R with the clock c
-params = (; β = 0.1, ..., S, I, R) # create new params with recorded S, I, R
+S = recorded(DynamicEntry, c, 100) # create a recorded number as S with the clock c
+I = recorded(DynamicEntry, c, [1]) # create a recorded vector as I with the clock c
+R = recorded(DynamicEntry, c, [0]) # create a recorded vector as R with the clock c
+α = recorded(StaticEntry, c, [0.5]) # create a recorded vector as α with the clock c
+params = (; β = 0.1, μ = 0.001, ..., S, I, R, α) # create new params with recorded S, I, R, α
 gillespie(c, params, reactions) # run the simulation with the clock and new params
 ```
 
 More information about `RecordedArrays`, see its
 [documentation](https://wangl-cc.github.io/RecordedArrays.jl/dev).
-
-Examples and references of this package can be found in
-[documentation](https://wangl-cc.github.io/EvolutionaryModelingTools.jl/dev).

docs/src/example.md

@@ -47,8 +47,7 @@ ps = (; r, d, K, x) # define the parameters for the model
 rs = (growth, death, comp) # define the reactions for the model
 
 # simulate
-Random.seed!(1)
-ps_new, t, state = gillespie(c, ps, rs)
+ps_new, t, state = gillespie(MersenneTwister(1), c, ps, rs)
 
 # plot
 timeseries(ps_new.x;
@@ -134,12 +133,88 @@ ps = (; r, d, K, μ, x, m) # define the parameters for the model
 rs = (growth, mutation, competition) # define the reactions for the model
 
 # simulate
-Random.seed!(1)
-ps_new, t, state = gillespie(check_extinction, c, ps, rs)
+ps_new, t, state = gillespie(check_extinction, MersenneTwister(1), c, ps, rs)
 
 # plot
 timeseries(ps_new.x;
     grid=false, frame=:box, legend=false,
     title="Population Dynamics", xlabel="Time", ylabel="Population Size",
 )
 ```
+
+## SIR model with vital dynamics and logistic population
+
+A simple SIR model with vital dynamics and logistic population,
+where the host population follows the logistic growth model.
+
+This model may contains more reactions than the previous example,
+which can be generated easily by `@reaction_eq`.
+
+The deterministic differential equation of this model is:
+
+```math
+\begin{aligned}
+\dot{S} &= r (S + I + R) - S(d + c (S + I + R) + \beta I) \\
+\dot{I} &= I (\beta S - d - c (S + I + R) - \nu) \\
+\dot{R} &= R (\nu - d - c (S + I + R))
+\end{aligned}
+```
+
+```@example sir_with_vital_dynamics_and_logistic_population
+using RecordedArrays
+using EvolutionaryModelingTools
+using Random
+using Plots
+
+const T = 100.0
+const β = 0.005
+const ν = 0.01
+const α = 0.2
+const r = 0.5
+const d = 0.1
+const c = 0.001
+const S = 100
+const I = 10
+const R = 0
+
+# epidemic dynamics
+@reaction_eq infection β S + I --> I + I
+@reaction_eq recovery ν I --> R
+@reaction_eq virulence α I --> 0 # death of infection host caused by virus
+# demography dynamics, generate reactions with @eval
+for sym in (:S, :I, :R)
+    r_name = Symbol(:growth_, sym)
+    d_name = Symbol(:death_, sym)
+    @eval @reaction_eq $r_name r $sym --> S + $sym # growth
+    @eval @reaction_eq $d_name d $sym --> 0 # death
+    for sym′ in (:S, :I, :R)
+        c_name = Symbol(:competition_, sym, sym′)
+        @eval @reaction_eq $c_name c $sym + $sym′ --> $sym′ # competition
+    end
+end
+
+const REACTIONS = (
+    infection, recovery, virulence,
+    growth_S, growth_I, growth_R,
+    death_S, death_I, death_R,
+    competition_SS, competition_SI, competition_SR,
+    competition_IS, competition_II, competition_IR,
+    competition_RS, competition_RI, competition_RR
+)
+
+clock = ContinuousClock(T)
+ps = (;β, ν, α, r, d, c,
+    S=recorded(DynamicEntry, clock, S),
+    I=recorded(DynamicEntry, clock, I),
+    R=recorded(DynamicEntry, clock, R),
+)
+ps_new, t, term = gillespie(MersenneTwister(1), clock, ps, REACTIONS)
+
+plt = plot(grid=false, frame=:box,
+    title="Population Dynamics", xlabel="Time", ylabel="Population Size"
+)
+timeseries!(plt, ps_new.S; label="S")
+timeseries!(plt, ps_new.I; label="I")
+timeseries!(plt, ps_new.R; label="R")
+plt
+```

docs/src/index.md

@@ -0,0 +1 @@
+../../README.md

docs/src/index.md

@@ -4,9 +4,9 @@ using EvolutionaryModelingTools
 using EvolutionaryModelingTools.Scalar
 
 # parameters
-const T = 100.0
-const β = 0.001
-const ν = 0.01
+const T = 10.0
+const β = 0.005
+const ν = 0.1
 const α = 0.2
 const r = 0.5
 const d = 0.1