gillepsie.jl 2.53 KB
Newer Older
Victor Boussange's avatar
Victor Boussange committed
1
cd(@__DIR__)
Victor's avatar
Victor committed
2
using Random
3
Random.seed!(0)
Victor's avatar
Victor committed
4
5
6
7
using LightGraphs
using Test
using Revise,ABMEv
using UnPack,JLD2
Victor Boussange's avatar
Victor Boussange committed
8

Victor's avatar
Victor committed
9
myspace = (RealSpace{1,Float64}(),)
10
11
sigma_K = .9;
sigma_a = .7;
12
K0 = 1000;
Victor's avatar
Victor committed
13
14
15
16
17
b(X) = gaussian(X[1],0.,sigma_K)
d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
Victor's avatar
Victor committed
18
tend = 1.5
Victor's avatar
Victor committed
19
20
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax

21
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
Victor's avatar
Victor committed
22
23
w0 = World(myagents,myspace,p,0.)
w1 = copy(w0)
24
25
26
@info "Running simulation with Gillepsie algorithm"
@time sim = run!(w1,Gillepsie(),tend)
@test typeof(sim) <: Simulation
Victor's avatar
Victor committed
27

Victor Boussange's avatar
Victor Boussange committed
28
# @save "gillepsie_test.jld2" world_alive
29
# @load "gillepsie_test.jld2" world_alive
Victor Boussange's avatar
Victor Boussange committed
30
## Testing
31
32
@testset "Gillepsie Algorithm" begin
        @testset "Testing global functioning" begin
33
                @test get_size(sim) > 1
Victor's avatar
Victor committed
34
                @test get_tend(sim) >= tend
35
36
        end
        ## Comparing simulation
37
38
39
40
        # @testset "Matching new vs old results " begin
        #         xarray = get_x(w1,1);xarray_test = get_x(world_alive,1);
        #         @test prod(xarray .≈ xarray_test)
        # end
41

42
        @testset "Testing update rates matrix" begin
43
44
45
46
47
48
49
50
51
52
                @testset "birth event" begin
                        myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie())
                        mum_idx = 1
                        updateBirthEvent!(w0,Gillepsie(),1)
                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
                        update_rates!(w1,Gillepsie());
                        bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1))
                        @test prod(bs_end . bs_recalculated)
                        @test prod(ds_end . ds_recalculated)
Victor Boussange's avatar
Victor Boussange committed
53

54
                end
Victor Boussange's avatar
Victor Boussange committed
55

56
57
58
59
60
61
62
63
64
65
66
                @testset "death event" begin
                        myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie())
                        mum_idx = 1
                        updateDeathEvent!(w0,Gillepsie(),1)
                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
                        update_rates!(w1,Gillepsie());
                        bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1))
                        @test prod(bs_end . bs_recalculated)
                        @test prod(ds_end . ds_recalculated)
                end
67
        end
68
end