ABMEv_Gillepsie.jl 2.38 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# In this file lie the function for Gillepsie algorithm


"""
    give_birth(a::Agent,p)
Used for Gillepsie setting
"""
function give_birth(a::Agent,p::Dict,reflected)
    new_a = copy(a)
    increment_x!(new_a,p,reflected=reflected)
    return new_a
end

function update_afterbirth_std!(world,C,idx::Int,p::Dict) where T
    traits = get_x.(world)
    # updating competition only the two columns corresponding to agent idx
Victor Boussange's avatar
Victor Boussange committed
17
    α = p["alpha"];K=p["K"]
18
    for i in 1:length(traits)
Victor Boussange's avatar
Victor Boussange committed
19
        C[i,idx] = α(traits[i],traits[idx])
20
21
22
23
        C[idx,i] = C[i,idx]
    end
    # updating death rate only the two columns corresponding to agent idx
    for (i,a) in enumerate(world)
Victor Boussange's avatar
Victor Boussange committed
24
        a.d += C[i,idx]
25
26
    end
    # Now updating new agent
Victor Boussange's avatar
Victor Boussange committed
27
    world[idx].d = sum(C[idx,:])
28
    world[idx].b = K(traits[idx])
29
30
31
32
33
34
end

function update_afterdeath_std!(world,C,idx::Int,p::Dict) where T
    traits = get_x.(world)
    # updating death rate only the two columns corresponding to agent idx
    for (i,a) in enumerate(world)
Victor Boussange's avatar
Victor Boussange committed
35
        a.d -= C[i,idx]
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    end
    # updating competition only the two columns corresponding to agent idx
    for i in 1:length(traits)
        C[i,idx] = .0
        C[idx,i] = .0
    end
end

"""
    function updateWorld_G!(world,C,tspan,p)
Updating rule for gillepsie setting.
Returning dt drawn from an exponential distribution with parameter the total rates of events.
"""
function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
    # 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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
    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
        I = Int(length(events_weights) / 2)
        if i_event <= I
            # DEATH EVENT
            idx_offspring = idx_world[i_event]
            world[idx_offspring] = missing
            update_afterdeath_std!(world_alive,C,idx_offspring,p)
        else
            # birth event
            idx_offspring = findfirst(ismissing,world)
            world[idx_offspring] = give_birth(world[idx_world[i_event-I]],p,reflected)
            update_afterbirth_std!(world_alive,C,idx_offspring,p)
71

72
73
74
75
        end
        return dt
    else
        return -1.
76
77
    end
end