SIR.jl 1.14 KB
Newer Older
1
2
3
4
5
6
7
8

import StaticArrays: SVector, SMatrix, @SMatrix

State = SVector{3, Int}
Parameters = SVector{2, Real}

d=3
k=2
9
10
11
dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
dict_p = Dict("ki" => 1, "kr" => 2)
l_tr = ["R1","R2"]
12
13
p = Parameters(0.0012, 0.05)
x0 = State(95, 5, 0)
14
15
t0 = 0.0
function f(xn::State, tn::Real, p::Parameters)
16
17
18
19
20
    a1 = p[1] * xn[1] * xn[2]
    a2 = p[2] * xn[2]
    l_a = SVector(a1, a2)
    asum = sum(l_a)
    # column-major order
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
21
22
23
    l_nu = @SMatrix [-1 0;
                     1 -1;
                     0 1]
24
25
26
27
28
29
    
    u1, u2 = rand(), rand()
    tau = - log(u1) / asum
    b_inf = 0.0
    b_sup = a1
    reaction = 0
30
    for i = 1:2 
31
32
33
34
35
36
37
38
39
40
41
42
43
        if b_inf < asum*u2 < b_sup
            reaction = i
            break
        end
        b_inf += l_a[i]
        b_sup += l_a[i+1]
    end
 
    nu = l_nu[:,reaction]
    xnplus1 = State(xn[1]+nu[1], xn[2]+nu[2], xn[3]+nu[3])
    tnplus1 = tn + tau
    transition = "R$(reaction)"

44
    return xnplus1, tnplus1, transition
45
end
46
47
is_absorbing_sir(p::Parameters,xn::State) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) == 0.0
g = SVector("I")
48

49
50
SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f,is_absorbing_sir; g=g)
export SIR
51