Add of a function in models that creates new model objects + tests.

Add of a major test based on Cosmos: tests/cosmos/distance_F/ER_1D.jl

The LHA distance F of this package is tested with Cosmos.
For several parameters in R1, R2, R3 experiments (see Bentriou,
Ballarini, Cournede 2018), Cosmos is called to
compute an average distance. Then we collect the number of simulations
and the estimated value computed by Cosmos.
We compute an average distance by the package methods with the same
number of simulations done by Cosmos

Then the two average distances are compared to each others with
regard to the width of the confidence interval set in Cosmos.
Test passed.

core/MarkovProcesses.jl

@@ -26,7 +26,7 @@ export isbounded, isaccepted, check_consistency
 export load_model, get_module_path
 
 # Utils
-export get_module_path
+export get_module_path, cosmos_get_values
 
 include("common.jl")
 

core/utils.jl

 
 get_module_path() = dirname(dirname(pathof(@__MODULE__)))
 
+function cosmos_get_values(name_file::String) 
+    output_file = open(name_file)
+    dict_values = Dict{String}{Float64}()
+    while !eof(output_file)
+        line = readline(output_file)
+        splitted_line = split(line, ':')
+        if (length(splitted_line) > 1) && tryparse(Float64, splitted_line[2]) !== nothing 
+            dict_values[splitted_line[1]] = parse(Float64, splitted_line[2])
+        end
+    end
+    close(output_file)
+    return dict_values
+end
+

models/ER.jl

@@ -46,5 +46,13 @@ isabsorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) =
 g_ER = ["P"]
 
 ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER)
-export ER
+
+function create_ER(new_p::Vector{Float64})
+    ER_new = deepcopy(ER)
+    @assert length(ER_new.p) == length(new_p)
+    set_param!(ER_new, new_p)
+    return ER_new
+end
+
+export ER, create_ER
 

models/SIR.jl

@@ -45,5 +45,13 @@ isabsorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p
 g_SIR = ["I"]
 
 SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR)
-export SIR
+
+function create_SIR(new_p::Vector{Float64})
+    SIR_new = deepcopy(SIR)
+    @assert length(SIR_new.p) == length(new_p)
+    set_param!(SIR_new, new_p)
+    return SIR_new
+end
+
+export SIR, create_SIR
 

tests/cosmos/distance_F/ER.gspn

+const int E_init=100;
+const int S_init=100;
+const int ES_init=0;
+const int P_init=0;
+
+const double k_1=1.0;
+const double k_2=1.0;
+const double k_3=1.0;
+
+
+NbPlaces = 4;
+NbTransitions = 3;
+
+PlacesList = { 
+E,S, ES, P
+} ;
+
+TransitionsList = { 
+TR1,   TR2,  TR3
+} ;
+
+Marking={
+(E , E_init); (S , S_init); (ES , ES_init);
+(P , P_init); 
+};
+
+Transitions={
+(TR1,EXPONENTIAL(k_1*E*S),1,1, SINGLE);
+(TR2,EXPONENTIAL(k_2*ES),1,1, SINGLE);
+(TR3,EXPONENTIAL(k_3*ES),1,1, SINGLE);
+};
+
+InArcs={
+(E, TR1,1); (S, TR1,1); 
+(ES, TR2,1); (ES, TR3,1); 
+};
+
+OutArcs={
+(TR1, ES,1); 
+(TR2, E,1);
+(TR2, S,1); 
+(TR3, E,1); 
+(TR3, P,1); 
+};
+

tests/cosmos/distance_F/ER_1D.jl

+
+
+using MarkovProcesses
+absolute_path = get_module_path() * "/tests/cosmos/"
+
+# Remarque: la valeur estimée par Cosmos varie plus que de 0.01
+
+# Values x1, x2  t1, t2
+dict_exp = Dict(
+                "R1" => [50, 75, 0.025, 0.05],
+                "R2" => [50, 75, 0.05, 0.075],
+                "R3" => [25, 50, 0.05, 0.075]
+               )
+
+str_model = "ER"
+if str_model == "ER"
+    load_model(str_model)
+    model = ER
+end
+load_automaton("automaton_F")
+
+test_all = true
+width = 0.1
+level = 0.95
+l_exp = ["R1", "R2", "R3"]
+#l_exp = ["R2"]
+for exp in l_exp
+    if !(exp in keys(dict_exp))
+        println("Unrecognized experiment: <<$exp>>")
+        continue
+    end
+    val_exp = dict_exp[exp]
+    x1, x2, t1, t2 = val_exp[1], val_exp[2], val_exp[3], val_exp[4]
+    A_F = create_automaton_F(model, x1, x2, t1, t2, "P")  
+    l_k3 = 0:10:100
+    nb_param = length(l_k3)
+    l_dist_cosmos = zeros(nb_param)
+    l_dist_pkg = zeros(nb_param)
+    for i in 1:nb_param
+        # Cosmos estimation
+        k3 = l_k3[i]
+        command = `Cosmos $(absolute_path * "distance_F/" * str_model * ".gspn") 
+        $(absolute_path * "distance_F/dist_F_"  * str_model * ".lha") --njob $(ENV["JULIA_NUM_THREADS"]) 
+        --const k_3=$(k3),x1=$x1,x2=$x2,t1=$t1,t2=$t2 
+        --level $(level) --width $(width) 
+        --verbose 0` 
+        run(pipeline(command, stderr=devnull))
+        dict_values = cosmos_get_values("Result_dist_F_$(str_model).res")
+        l_dist_cosmos[i] = dict_values["Estimated value"]
+        nb_sim = dict_values["Total paths"]
+        nb_accepted = dict_values["Accepted paths"]
+        nb_sim = convert(Int, nb_sim)
+        # MarkovProcesses estimation
+        set_param!(ER, "k3", convert(Float64, k3))
+        sync_ER = ER*A_F
+        l_dist_i_threads = zeros(Threads.nthreads()) 
+        Threads.@threads for j = 1:nb_sim
+            σ = simulate(sync_ER)
+            l_dist_i_threads[Threads.threadid()] += (σ.S)["d"]
+        end
+        l_dist_pkg[i] = sum(l_dist_i_threads) / nb_sim
+        global test_all = test_all && isapprox(l_dist_cosmos[i], l_dist_pkg[i]; atol = width)
+    end
+
+end
+
+rm("Result_dist_F_$(str_model).res")
+rm("Result.res")
+
+return test_all
+

tests/cosmos/distance_F/dist_F_ER.lha

+NbVariables = 3;
+NbLocations = 4;
+
+const int x1 = 50;
+const int x2 = 75;
+
+const double t1 = 0.025;
+const double t2 = 0.05;
+
+VariablesList = {t,d,n};
+
+LocationsList = {l0, l1, l2, l3};
+
+AVG(Last(d));
+
+InitialLocations={l0};
+FinalLocations={l2};
+
+Locations={
+(l0,   TRUE , (t:1));
+(l1,   TRUE , (t:1));
+(l2,   TRUE , (t:1));
+(l3,   TRUE , (t:1));
+};
+
+Edges={
+((l0,l1), #, t>=0, {n=P});
+((l1,l3), #, t<=t1 & n<=x1-1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
+((l1,l3), #, t<=t1 & n>=x2+1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
+((l1,l3), #, n>=x1 & n<=x2, {d=0});
+((l1,l3), #, n<=x1-1 & t >= t1 & t<=t2, {d=min(d ,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) )});
+((l1,l3), #, n>=x2+1 & t >= t1 & t<=t2, {d=min(d,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) ) });
+((l3,l1), ALL, t>=0, {n=P});
+
+((l3,l2), #, t>=t2 ,#);
+
+((l1,l2),#, n>=x1 & n<=x2 & t>=t1 & t<=t2 ,{d=0});
+((l1,l2),#, t>=t2 , #);
+((l1,l2),#, d=0 & t>=t1, #);
+};

tests/run_all.jl

@@ -3,4 +3,5 @@ include("run_unit.jl")
 include("run_simulation.jl")
 include("run_dist_lp.jl")
 include("run_automata.jl")
+include("run_cosmos.jl")
 

tests/run_cosmos.jl

+
+using Test
+
+@testset "Cosmos tests" begin
+    @test include("cosmos/distance_F/ER_1D.jl")
+end
+

tests/run_unit.jl

@@ -20,5 +20,6 @@ using Test
     @test include("unit/check_model_consistency.jl")
     @test include("unit/set_param.jl")
     @test include("unit/side_effects_1.jl")
+    @test include("unit/create_models.jl")
 end
 

tests/unit/create_models.jl

+
+using MarkovProcesses
+
+load_model("SIR")
+load_model("ER")
+
+test_all = true
+
+new_p = [0.01, 0.2]
+new_SIR = create_SIR(new_p)
+test_all = test_all && new_SIR.p == new_p && new_SIR !== SIR
+
+p_SIR = [0.0012, 0.05]
+new_SIR = create_SIR(p_SIR)
+test_all = test_all && new_SIR.p == p_SIR && new_SIR !== SIR
+
+new_p = [0.01, 0.2, 20.0]
+new_ER = create_ER(new_p)
+test_all = test_all && new_ER.p == new_p && new_ER !== ER
+
+p_ER = [1.0, 1.0, 1.0]
+new_ER = create_ER(p_ER)
+test_all = test_all && new_ER.p == p_ER && new_ER !== ER
+
+return true
+