SIR.jl 1.24 KB
Newer Older
1
2
3
4
5

import StaticArrays: SVector, SMatrix, @SMatrix

d=3
k=2
6
7
dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
dict_p = Dict("ki" => 1, "kr" => 2)
8
9
10
11
l_tr_SIR = ["R1","R2"]
p_SIR = [0.0012, 0.05]
x0_SIR = [95, 5, 0]
t0_SIR = 0.0
12
function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
13
           xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
14
15
16
17
18
    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
19
20
21
    l_nu = @SMatrix [-1 0;
                     1 -1;
                     0 1]
22
23
24
25
26
27
    
    u1, u2 = rand(), rand()
    tau = - log(u1) / asum
    b_inf = 0.0
    b_sup = a1
    reaction = 0
28
    for i = 1:2 
29
30
31
32
33
34
35
36
        if b_inf < asum*u2 < b_sup
            reaction = i
            break
        end
        b_inf += l_a[i]
        b_sup += l_a[i+1]
    end
 
37
    nu = view(l_nu, :, reaction) # macro for avoiding a copy
38
39
40
41
42
    for i = 1:3
        mat_x[idx,i] = xn[i]+nu[i]
    end
    l_t[idx] = tn + tau
    l_tr[idx] = "R$(reaction)"
43
end
44
is_absorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
45
g_SIR = ["I"]
46

47
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,is_absorbing_SIR; g=g_SIR)
48
export SIR
49