Gillepsie.jl 3.36 KB
Newer Older
1
# In this file lie the function for Gillepsie algorithm
Victor's avatar
Victor committed
2

3
4
"""
$(TYPEDEF)
Victor's avatar
Victor committed
5
6
7
8
Gillespie algorithm.

Denote by ``b_i`` and ``d_i`` the birth and death rates of 
agent ``i``. The total rate is given by the sum of all individual rates
Victor's avatar
Victor committed
9
10
11
``R(t) = \\left[ \\sum_i b_i(t) + d_i(t) \\right]``.
A particular event, birth or death, 
is chosen at random with a probability equal to the rate of this event divided by the total rate ``R``.
Victor's avatar
Victor committed
12
13
14


# The original article by Gillsepie:
Victor's avatar
Victor committed
15
16
[**A general method for numerically simulating the stochastic 
time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub)
Victor's avatar
Victor committed
17

18
"""
Victor's avatar
Victor committed
19
struct Gillepsie <: AbstractAlg end
20
export Gillepsie
21

Victor's avatar
Victor committed
22
function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int,b,d)
23
    # updating competition only the two columns corresponding to agent idx
Victor's avatar
Victor committed
24
    offspring = give_birth(mum_idx,w)
Victor's avatar
Victor committed
25
    x_offspring = get_x(offspring)
Victor's avatar
Victor committed
26
    _agents = agents(w)
Victor's avatar
Victor committed
27
    for a in _agents
Victor's avatar
Victor committed
28
29
30
31
        # Warning : symmetric competition
        tmp = d(get_x(a),x_offspring,w.t)
        a.d += tmp
        offspring.d += tmp
32
33
    end
    # Now updating new agent
34
    offspring.b = b(x_offspring,w.t)
Victor's avatar
Victor committed
35
    addAgent!(w,offspring)
36
37
end

Victor's avatar
Victor committed
38
function updateDeathEvent!(world::World,::Gillepsie,i_event::Int,d)
Victor's avatar
Victor committed
39
    x_death = get_x(world[i_event])
40
    # updating death rate only the two columns corresponding to agent idx
Victor's avatar
Victor committed
41
42
    removeAgent!(world,i_event)
    for a in agents(world)
43
        a.d -= d(get_x(a),x_death,world.t)
44
45
46
    end
end

47
"""
Victor's avatar
Victor committed
48
    update_rates!(w::World,::Gillepsie,b,d)
49
50
51
52
This standard updates takes
    - competition kernels of the form α(x,y) and
    - carrying capacity of the form K(x)
"""
Victor's avatar
Victor committed
53
function  update_rates!(w::World,::Gillepsie,b,d)
54
55
56
57
58
59
60
61
    _agents = agents(w)
    traits = get_x.(_agents)
    # traits = get_xhist.(world)
    n = size(w)
    D = zeros(n)
    # Here you should do a shared array to compute in parallel
    for i in 1:(n-1)
        for j in i+1:n
62
            C = d(traits[i],traits[j],w.t)
63
64
65
66
67
68
69
            D[i] += C
            D[j] += C
        end
    end
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(_agents)
        a.d = D[i]
70
        a.b = b(traits[i],w.t)
71
72
73
    end
end

74
"""
Victor's avatar
Victor committed
75
    function updateWorld!(w::World{A,S,T},g::G,b,d)
76
Updating rule for gillepsie setting.
Victor's avatar
Victor committed
77
Returning `dt` drawn from an exponential distribution with parameter the total rates of events.
78
"""
Victor's avatar
Victor committed
79
function updateWorld!(w::World{A,S,T},g::G,b,d) where {A,S,T,G <: Gillepsie}
80
    # total update world
81
82
    alive = agents(w)
    events_weights = ProbabilityWeights(vcat(get_d.(alive),get_b.(alive)))
83
    # Total rate of events
84
85
86
     = sum(events_weights)
    dt = - log(rand(T))/
    update_clock!(w,dt)
87
88
    if dt > 0.
        i_event = sample(events_weights)
Victor's avatar
Victor committed
89
        # This is ok since size is a multiple of 2
90
        n = size(w)
Victor's avatar
Victor committed
91
        if i_event <= n
92
            # DEATH EVENT
93
            # In this case i_event is also the index of the individual to die in the world_alive
Victor's avatar
Victor committed
94
            updateDeathEvent!(w,g,i_event,d)
95
96
        else
            # birth event
Victor's avatar
Victor committed
97
            # i_event - n is also the index of the individual to give birth in the world_alive
Victor's avatar
Victor committed
98
            updateBirthEvent!(w,g,i_event-n,b,d)
99
100
101
        end
        return dt
    else
Victor's avatar
Victor committed
102
        return -1
103
104
    end
end
105
106
107
108
109
110
111

"""
    clean_world(world::Array{T}) where {T <: Agent}
Get rid of missing value in `world`
"""
#TODO : put some type specs here
clean_world(world) = collect(skipmissing(world))