To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

ABMEv_Gillepsie.jl 3.16 KB
Newer Older
1
# In this file lie the function for Gillepsie algorithm
2
3
4
"""
$(TYPEDEF)
"""
Victor's avatar
Victor committed
5
struct Gillepsie <: AbstractAlg end
6
export Gillepsie
7
"""
Victor's avatar
Victor committed
8
    function give_birth(a::Agent,t,p::Dict)
9
10
Used for Gillepsie setting
"""
11
function give_birth(mum_idx::Int,w::World)
Victor's avatar
Victor committed
12
13
    new_a = copy(w[mum_idx])
    increment_x!(new_a,space(w),parameters(w),time(w))
14
15
16
    return new_a
end

Victor's avatar
Victor committed
17
function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int)
18
    # updating competition only the two columns corresponding to agent idx
Victor's avatar
Victor committed
19
20
    @unpack d,b = parameters(w)
    offspring = give_birth(mum_idx,w)
Victor's avatar
Victor committed
21
    x_offspring = get_x(offspring)
Victor's avatar
Victor committed
22
    _agents = agents(w)
Victor's avatar
Victor committed
23
    for a in _agents
24
        a.d += d(get_x(a),x_offspring,w.t)
25
26
    end
    # Now updating new agent
27
28
    offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring),w.t)) #- d(x_offspring,x_offspring)
    offspring.b = b(x_offspring,w.t)
Victor's avatar
Victor committed
29
    addAgent!(w,offspring)
30
31
end

Victor's avatar
Victor committed
32
function updateDeathEvent!(world::World,::Gillepsie,i_event::Int)
33
    @unpack d = parameters(world)
Victor's avatar
Victor committed
34
    x_death = get_x(world[i_event])
35
    # updating death rate only the two columns corresponding to agent idx
Victor's avatar
Victor committed
36
37
    removeAgent!(world,i_event)
    for a in agents(world)
38
        a.d -= d(get_x(a),x_death,world.t)
39
40
41
    end
end

42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
"""
    update_rates_std!(world,p::Dict,t::Float64)
This standard updates takes
    - competition kernels of the form α(x,y) and
    - carrying capacity of the form K(x)
"""
function  update_rates!(w::World,::Gillepsie)
    @unpack b,d = parameters(w)
    _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
58
            C = d(traits[i],traits[j],w.t)
59
60
61
62
63
64
65
            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]
66
        a.b = b(traits[i],w.t)
67
68
69
    end
end

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

"""
    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))