gillepsie.jl 2.56 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;
13
14
b(X,t) = gaussian(X[1],0.,sigma_K)
d(X,Y,t) = gaussian(X[1],Y[1],sigma_a)/K0
15
D = [1e-2]
Victor's avatar
Victor committed
16
17
mu = [.1]
NMax = 10000
Victor's avatar
Victor committed
18
tend = 1.5
Victor's avatar
Victor committed
19
p = Dict{String,Any}();@pack! p = D,mu,NMax
Victor's avatar
Victor committed
20

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
@info "Running simulation with Gillepsie algorithm"
Victor's avatar
Victor committed
25
@time sim = run!(w1,Gillepsie(),tend,b,d)
26
@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
                @testset "birth event" begin
44
                        myagents = [Agent(myspace,[0],ancestors=true,rates=true) for i in 1:K0]
Victor's avatar
Victor committed
45
                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie(),b,d,)
46
                        mum_idx = 1
Victor's avatar
Victor committed
47
                        updateBirthEvent!(w0,Gillepsie(),1,b,d)
48
                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
Victor's avatar
Victor committed
49
                        update_rates!(w1,Gillepsie(),b,d);
50
51
52
                        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
                @testset "death event" begin
57
                        myagents = [Agent(myspace,[0],ancestors=true,rates=true) for i in 1:K0]
Victor's avatar
Victor committed
58
                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie(),b,d)
59
                        mum_idx = 1
Victor's avatar
Victor committed
60
                        updateDeathEvent!(w0,Gillepsie(),1,d)
61
                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
Victor's avatar
Victor committed
62
                        update_rates!(w1,Gillepsie(),b,d);
63
64
65
66
                        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