ABMEv_Gillepsie.jl 2.38 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# 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)
    # ! carefull you also need to change get_xhist
    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
    for i in 1:length(traits)
        C[i,idx] = α(traits[i],traits[idx],p["n_alpha"],p["sigma_a"])
        C[idx,i] = C[i,idx]
    end
    # updating death rate only the two columns corresponding to agent idx
    for (i,a) in enumerate(world)
        a.d += C[i,idx] / p["K0"]
    end
    # Now updating new agent
    world[idx].d = sum(C[idx,:]) / p["K0"]
    world[idx].b = K(traits[idx][:,end],1,p["n_K"],p["sigma_K"])
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)
        a.d -= C[i,idx] / p["K0"]
    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
    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)

    end
    return dt
end