ABMEv_Gillepsie.jl 2.44 KB
Newer Older
1
2
3
4
5
6
7
# In this file lie the function for Gillepsie algorithm


"""
    give_birth(a::Agent,p)
Used for Gillepsie setting
"""
Victor's avatar
Victor committed
8
function give_birth(a::Agent,p::Dict)
9
    new_a = copy(a)
Victor's avatar
Victor committed
10
    increment_x!(new_a,p)
11
12
13
    return new_a
end

14
function update_afterbirth_std!(world,idx_offspring,p::Dict) where T
15
    # updating competition only the two columns corresponding to agent idx
16
17
18
19
    α = p["alpha"];K=p["K"];
    x_offspring = get_x(world[idx_offspring])
    for a in skipmissing(world)
        a.d += α(get_x(a),x_offspring)
20
21
    end
    # Now updating new agent
22
    world[idx_offspring].d = sum(α.(get_x.(skipmissing(world)),Ref(x_offspring))) - α(x_offspring,x_offspring)
23
    world[idx_offspring].b = K(x_offspring)
24
25
end

26
27
function update_afterdeath_std!(world,x_death,p::Dict) where T
    α = p["alpha"]
28
    # updating death rate only the two columns corresponding to agent idx
29
30
    for a in skipmissing(world)
        a.d -= α(get_x(a),x_death)
31
32
33
34
    end
end

"""
Victor's avatar
Victor committed
35
    function updateWorld_G!(world,tspan,p)
36
37
Updating rule for gillepsie setting.
Returning dt drawn from an exponential distribution with parameter the total rates of events.
38
39
 # Args
 tspan is for now not used but might be used for dynamic landscape
40
"""
Victor's avatar
Victor committed
41
function updateWorld_G!(world,p,update_rates!,tspan)
42
43
44
45
46
47
    # total update world
    world_alive = skipmissing(world)
    idx_world = collect(eachindex(world_alive))
    # Total rate of events
    U = sum(get_d.(world_alive) .+ get_b.(world_alive))
    dt = - log(rand(Float64))/U
48
49
50
51
    if dt > 0.
        events_weights = ProbabilityWeights(vcat(get_d.(world_alive),get_b.(world_alive)))
        i_event = sample(events_weights)
        # This is ok since length is a multiple of 2
52
53
        N = length(world) - count(ismissing,world)
        if i_event <= N
54
            # DEATH EVENT
55
            # In this case i_event is also the index of the individual to die in the world_alive
56
            idx_offspring = idx_world[i_event]
57
58
            x_death = get_x(world[idx_offspring])
            update_afterdeath_std!(world,x_death,p)
59
60
61
62
            world[idx_offspring] = missing
        else
            # birth event
            idx_offspring = findfirst(ismissing,world)
63
64
            # i_event - N is also the index of the individual to give birth in the world_alive
            mum = world[idx_world[i_event-N]]
Victor's avatar
Victor committed
65
            world[idx_offspring] = give_birth(mum,p)
66
            update_afterbirth_std!(world,idx_offspring,p)
67
68
69
70
        end
        return dt
    else
        return -1.
71
72
    end
end