Commit f0d20414 authored by Thom Badings's avatar Thom Badings
Browse files

Chance in coin probs

parent 18f66cdf
......@@ -4,35 +4,35 @@ import numpy as np
template_header = """
mdp
const int MAXX = 10;
const int MAXY = 10;
const int MAXZ = 10;
const int MAXX = 15;
const int MAXY = 15;
const int MAXZ = 15;
formula valid = (1 < x & x < MAXX-1 & 1 < y & y < MAXY-1 & 1 < z & z < MAXZ-1) & ok=1;
formula attarget = x > MAXX - 3;
formula attarget2 = x > MAXX - 1 & y > MAXY - 2 & z > MAXZ - 3;
formula attarget2 = x > MAXX - 3;
formula attarget = x > MAXX - 2 & y > MAXY - 2 & z > MAXZ - 2;
formula crash = (xarea2 & yarea4 & zarea1) | (xarea1 & yarea3 & zarea2) | (xarea4 & yarea2 & zarea2) | (xarea3 & yarea0 & zarea4);
formula valid = (1 <= x & x <= MAXX-1 & 1 <= y & y <= MAXY-1 & 1 <= z & z <= MAXZ-1) & !crash & ok=1;
"""
control_module = """
module control
c : [0..3] init 0;
c : [0..1] init 0;
[move] c = 0 -> 1: (c'=1);
[env] c = 1 -> 1: (c'=2);
[chc] c = 2 -> 1: (c'=3);
[chd] c = 3 -> 1: (c'=0);
[env] c = 1 -> 1: (c'=0);
endmodule
"""
mod_head = """
module mod
ok : [0..1] init 1;
cond : [0..3] init 0;
wdir : [0..7] init 0;
cond : [0..0] init 0;
wdir : [0..0] init 0;
x : [0..MAXX] init 2;
y : [0..MAXY] init floor(MAXY/2);
z : [0..MAXZ] init floor(MAXZ/2);
y : [0..MAXY] init 2;
z : [0..MAXZ] init 2;
phase : [0..1] init 0;
......@@ -43,17 +43,8 @@ module mod
[move] z < MAXX -> 1: (z'=z+1);
[move] z > 0 -> 1: (z'=z-1);
[move] attarget -> 1: (phase'=1);
[chc] cond = 0 -> p00: (cond' = 0) + p01: (cond' = 1) + p02: (cond' = 2) + p03: (cond' = 3);
[chc] cond = 1 -> p10: (cond' = 0) + p11: (cond' = 1) + p12: (cond' = 2) + p13: (cond' = 3);
[chc] cond = 2 -> p20: (cond' = 0) + p21: (cond' = 1) + p22: (cond' = 2) + p23: (cond' = 3);
[chc] cond = 3 -> p30: (cond' = 0) + p31: (cond' = 1) + p32: (cond' = 2) + p33: (cond' = 3);
[chd] cond = 0 -> c0pdircm2: (wdir' = mod(wdir + 6,8)) + c0pdircm1: (wdir' = mod(wdir + 7,8)) + c0pdirc0: (wdir' = wdir) + c0pdircp1: (wdir' = mod(wdir + 1,8)) + c0pdircp2: (wdir' = mod(wdir + 2,8));
[chd] cond = 1 -> c1pdircm2: (wdir' = mod(wdir + 6,8)) + c1pdircm1: (wdir' = mod(wdir + 7,8)) + c1pdirc0: (wdir' = wdir) + c1pdircp1: (wdir' = mod(wdir + 1,8)) + c1pdircp2: (wdir' = mod(wdir + 2,8));
[chd] cond = 2 -> c2pdircm2: (wdir' = mod(wdir + 6,8)) + c2pdircm1: (wdir' = mod(wdir + 7,8)) + c2pdirc0: (wdir' = wdir) + c2pdircp1: (wdir' = mod(wdir + 1,8)) + c2pdircp2: (wdir' = mod(wdir + 2,8));
[chd] cond = 3 -> c3pdircm2: (wdir' = mod(wdir + 6,8)) + c3pdircm1: (wdir' = mod(wdir + 7,8)) + c3pdirc0: (wdir' = wdir) + c3pdircp1: (wdir' = mod(wdir + 1,8)) + c3pdircp2: (wdir' = mod(wdir + 2,8));
[env] !valid -> 1: (ok'=0) & (x'=0) & (y'=0) & (z'=0) & (phase'=0);
"""
......@@ -78,9 +69,9 @@ labels="""
#for element in itertools.product(*somelists):
# print(element)
xsplits = [4,6]
ysplits = [7]
zsplits = [5]
xsplits = [2,4,7,9,11]
ysplits = [3,5,7,10,13]
zsplits = [2,4,8,10,12]
x_area_splits = [0] + xsplits + ["MAXX"]
y_area_splits = [0] + ysplits + ["MAXY"]
......@@ -95,14 +86,14 @@ y_area_definitions = _areadefinitions("y", y_area_splits)
z_area_definitions = _areadefinitions("z", z_area_splits)
area_definitions = "\n".join([x_area_definitions, y_area_definitions, z_area_definitions])
conditions = [0,1,2,3]
winddirs = [0,1,2,3,4,5,6,7]
conditions = [0]
winddirs = [0]
xzones = list(range(len(xsplits)+1))
yzones = list(range(len(ysplits)+1))
zzones = list(range(len(zsplits)+1))
def get_update(dir,eff):
if dir % 2 == 0:
if dir % 10 == 1:
# straight
affected_var, other_var = ("x","y") if dir in [0,4] else ("y","x")
main_var_sym = "+" if dir < 4 else "-"
......@@ -110,55 +101,81 @@ def get_update(dir,eff):
return "({}' = {}{}{})".format(affected_var, affected_var, main_var_sym, math.ceil(1)) + " & " + "({}' = {}{})".format(other_var,other_var,sec_var_op)
else:
# diagionals
assert dir in [1,3,5,7]
#assert dir in [1,3,5,7]
xsym = "+" if dir in [1,7] else "-"
ysym = "+" if dir < 4 else "-"
if eff == 0:
return "(x'=x) & (y'=y)"
if eff in [3,6]:
return "(x' = x{}{}) & (y' = y{}{})".format(xsym,1,ysym,1)
# if eff in [3,6]:
# return "(x' = x{}{}) & (y' = y{}{})".format(xsym,1,ysym,1)
# if eff == 1:
# if dir == 1:
# return ("(x'=x+1)")
# if dir == 3:
# return ("(y'=y+1)")
# if dir == 5:
# return ("(x'=x-1)")
# if dir == 7:
# return ("(y'=y-1)")
# if eff == 2:
# if dir == 1:
# return ("(y'=y+1)")
# if dir == 3:
# return ("(x'=x-1)")
# if dir == 5:
# return ("(y'=y-1)")
# if dir == 7:
# return ("(x'=x+1)")
# if eff == 4:
# if dir == 1:
# return ("(x'=x+1) & (y'=y+1)")
# if dir == 3:
# return ("(y'=y+1) & (x'=x-1)")
# if dir == 5:
# return ("(x'=x-1) & (y'=y-1)")
# if dir == 7:
# return ("(y'=y-1) & (x'=x+1)")
# if eff == 5:
# if dir == 1:
# return ("(y'=y+1) & (x'=x+1)")
# if dir == 3:
# return ("(x'=x-1) & (y'=y+1)")
# if dir == 5:
# return ("(y'=y-1) & (x'=x-1)")
# if dir == 7:
# return ("(x'=x+1) & (y'=y-1)")
if eff == 1:
if dir == 1:
return ("(x'=x+1)")
if dir == 3:
return ("(y'=y+1)")
if dir == 5:
return ("(x'=x-1)")
if dir == 7:
return ("(y'=y-1)")
return ("(x'=x+1)")
if eff == 2:
if dir == 1:
return ("(y'=y+1)")
if dir == 3:
return ("(x'=x-1)")
if dir == 5:
return ("(y'=y-1)")
if dir == 7:
return ("(x'=x+1)")
return ("(x'=x+2)")
if eff == 3:
return ("(x'=x-1)")
if eff == 4:
if dir == 1:
return ("(x'=x+1) & (y'=y+1)")
if dir == 3:
return ("(y'=y+1) & (x'=x-1)")
if dir == 5:
return ("(x'=x-1) & (y'=y-1)")
if dir == 7:
return ("(y'=y-1) & (x'=x+1)")
return ("(x'=x-2)")
if eff == 5:
if dir == 1:
return ("(y'=y+1) & (x'=x+1)")
if dir == 3:
return ("(x'=x-1) & (y'=y+1)")
if dir == 5:
return ("(y'=y-1) & (x'=x-1)")
if dir == 7:
return ("(x'=x+1) & (y'=y-1)")
return ("(y'=y+1)")
if eff == 6:
return ("(y'=y+2)")
if eff == 7:
return ("(y'=y-1)")
if eff == 8:
return ("(y'=y-2)")
if eff == 9:
return ("(y'=y+1) & (x'=x+1)")
if eff == 10:
return ("(y'=y-1) & (x'=x-1)")
if eff == 11:
return ("(y'=y+1) & (x'=x-1)")
if eff == 12:
return ("(y'=y-1) & (x'=x+1)")
assert False
num_effs=13
def create_command(c,w,x,y,z):
guard = "cond = {} & wdir = {} & xarea{} & yarea{} & zarea{} & valid".format(c, w, x, y, z)
updates = [get_update(w, i) for i in range(7)]
probs = ["pc{}w{}x{}y{}z{}eff{}".format(c, w, x, y, z, i) for i in range(7)]
updates = [get_update(w, i) for i in range(num_effs)]
probs = ["pc{}w{}x{}y{}z{}eff{}".format(c, w, x, y, z, i) for i in range(num_effs)]
weighted_updates = ["{}: {}".format(p, u) for p, u in zip(probs, updates)]
update_str = " + ".join(weighted_updates)
return "\t[env] {} -> {};".format(guard, update_str)
......@@ -180,12 +197,12 @@ commands = [create_command(c,w,x,y,z) for (c,w,x,y,z) in itertools.product(condi
def parameter_groups():
groups = []
for cfrom in range(4):
groups.append(["p{}{}".format(cfrom,cto) for cto in range(4)])
for cfrom in range(4):
groups.append(["c{}pdirc{}".format(cfrom,dirto) for dirto in ["m2","m1","0","p1","p2"]])
# for cfrom in range(1):
# groups.append(["p{}{}".format(cfrom,cto) for cto in range(1)])
# for cfrom in range(1):
# groups.append(["c{}pdirc{}".format(cfrom,dirto) for dirto in ["m2","m1","0","p1","p2"]])
for (c,w,x,y,z) in itertools.product(conditions, winddirs, xzones, yzones, zzones):
groups.append(["pc{}w{}x{}y{}z{}eff{}".format(c, w, x, y, z, i) for i in range(7)])
groups.append(["pc{}w{}x{}y{}z{}eff{}".format(c, w, x, y, z, i) for i in range(num_effs)])
#for r in range(1):
# groups.append(["RewardEnv","RewardChd"])
#for (w) in (winddirs):
......@@ -202,23 +219,55 @@ def parameter_definitions(groups,init):
return "\n".join(["const double {};".format(par) for group in groups for par in group])
def parameter_dirichlet_instantiations(groups):
def parameter_dirichlet_instantiations(groups,weather):
instantiations=dict()
for group in groups:
#print(group)
#if there is eff in the parameter group, create dirichlet weights non-uniformly
if any("eff" in group_item for group_item in group):
if weather == "uniform":
array = np.random.random_sample((len(group),))
array[0]=1*array[0]
array[1]=1*array[1]
array[2]=3*array[2]
array[3]=3*array[3]
array[4]=4*array[4]
array[5]=2*array[5]
array[6]=1*array[6]
#create uniformly
elif weather == "x-neg-bias":
if any("eff" in group_item for group_item in group):
array = np.random.random_sample((len(group),))
array[0]=1*array[0]
array[1]=1*array[1]
array[2]=1*array[2]
array[3]=2*array[3]
array[4]=2*array[4]
array[5]=1*array[5]
array[6]=1*array[6]
array[7]=1*array[7]
array[8]=1*array[8]
array[9]=1*array[9]
array[10]=2*array[10]
array[11]=2*array[11]
array[12]=1*array[12]
# print(array,len(array),group)
#create uniformly
else:
array = np.random.random_sample((len(group),))
elif weather == "y-pos-bias":
if any("eff" in group_item for group_item in group):
array = np.random.random_sample((len(group),))
array[0]=1*array[0]
array[1]=1*array[1]
array[2]=1*array[2]
array[3]=1*array[3]
array[4]=1*array[4]
array[5]=1*array[5]
array[6]=1*array[6]
array[7]=2*array[7]
array[8]=2*array[8]
array[9]=1*array[9]
array[10]=2*array[10]
array[11]=1*array[11]
array[12]=2*array[12]
# print(array,len(array),group)
#create uniformly
else:
array = np.random.random_sample((len(group),))
else:
array = np.random.random_sample((len(group),))
raise RuntimeError("Specificed weather is not compatible")
instantiations[tuple(group)]=array
return instantiations
......@@ -245,11 +294,12 @@ def parameter_dirichlet_samples(parameter_instantiations):
return instantiations_sample
# groups = parameter_groups()
groups = parameter_groups()
weather = "uniform"
# parameter_defs = parameter_definitions(groups,False)
# parameter_instantiations=parameter_dirichlet_instantiations(groups)
# parameter_samples=parameter_dirichlet_samples(parameter_instantiations)
parameter_defs = parameter_definitions(groups,False)
parameter_instantiations=parameter_dirichlet_instantiations(groups,weather)
parameter_samples=parameter_dirichlet_samples(parameter_instantiations)
# for element in parameter_samples:
# print(np.sum(parameter_samples[tuple(element)]))
......@@ -258,15 +308,15 @@ def parameter_dirichlet_samples(parameter_instantiations):
# with open("drone_par_param_example.nm",'w') as file:
# file.write(template_header)
# file.write(area_definitions)
# file.write("\n\n")
# file.write(parameter_defs)
# file.write(control_module)
# file.write(mod_head)
# file.write("\n".join(commands))
# file.write(mod_foot)
with open("drone_par_param_effect_example.nm",'w') as file:
file.write(template_header)
file.write(area_definitions)
file.write("\n\n")
file.write(parameter_defs)
file.write(control_module)
file.write(mod_head)
file.write("\n".join(commands))
file.write(mod_foot)
# file.write(reward_head)
# file.write("\n".join(rewards))
# file.write(reward_foot)
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par.nm --property-file drone_par.prctl --threshold 0.9 --bisimulation none --direction geq --num_samples 1000 --num_iter 5 --violation_prob 0.563
Building model from drone_par.nm
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ ^C
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 100 --num_iter 5 --weather "y-pos-bias" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 23513
Number of params before bisim: 932
Number of states after bisim: 23513
Number of params after bisim: 932
40%|██████████████████ | 2/5 [02:03<03:00, 60.29s/it]Average violation so far: 0.398
60%|███████████████████████████ | 3/5 [03:04<02:01, 60.72s/it]Average violation so far: 0.39666666666666667
80%|████████████████████████████████████ | 4/5 [04:12<01:02, 62.88s/it]Average violation so far: 0.39725
100%|█████████████████████████████████████████████| 5/5 [05:20<00:00, 64.10s/it]
Average solver time : 64.10429401397705
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|██████████████████ | 2/5 [00:28<00:42, 14.11s/it]Average violation so far: 0.225
60%|███████████████████████████ | 3/5 [00:42<00:28, 14.01s/it]Average violation so far: 0.22
80%|████████████████████████████████████ | 4/5 [00:57<00:14, 14.45s/it]Average violation so far: 0.225
100%|█████████████████████████████████████████████| 5/5 [01:12<00:00, 14.44s/it]
Average solver time : 14.442442750930786
printing specification violation array
[0.392, 0.404, 0.394, 0.399, 0.415]
[0.25, 0.2, 0.21, 0.24, 0.21]
probability of violating the spec less then the threshold for each iter:
[0.00227671 0.01913301 0.00336893 0.00840672 0.08523369]
[0.99999988 0.99976737 0.99993679 0.99999933 0.99993679]
average alpha_nu values of the array:
0.023683812925771354
0.9999280323951837
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 100 --num_iter 5 --weather "x-neg-bias" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|██████████████████ | 2/5 [00:28<00:43, 14.38s/it]Average violation so far: 0.25
60%|███████████████████████████ | 3/5 [00:43<00:28, 14.42s/it]Average violation so far: 0.24333333333333332
80%|████████████████████████████████████ | 4/5 [00:57<00:14, 14.46s/it]Average violation so far: 0.235
100%|█████████████████████████████████████████████| 5/5 [01:13<00:00, 14.62s/it]
Average solver time : 14.6180570602417
printing specification violation array
[0.24, 0.26, 0.23, 0.21, 0.16]
probability of violating the spec less then the threshold for each iter:
[0.99999933 0.99999998 0.9999966 0.99993679 0.98486986]
average alpha_nu values of the array:
0.9969605125533507
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 100 --num_iter 5 --weather "uniform" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|██████████████████ | 2/5 [00:26<00:38, 12.99s/it]Average violation so far: 0.03
60%|███████████████████████████ | 3/5 [00:40<00:26, 13.23s/it]Average violation so far: 0.04
80%|████████████████████████████████████ | 4/5 [00:53<00:13, 13.38s/it]Average violation so far: 0.04
100%|█████████████████████████████████████████████| 5/5 [01:07<00:00, 13.49s/it]
Average solver time : 13.488839817047118
printing specification violation array
[0.03, 0.03, 0.06, 0.04, 0.0]
probability of violating the spec less then the threshold for each iter:
[0.01513014 0.01513014 0.1216725 0.03337651 0.00077098]
average alpha_nu values of the array:
0.03721605539282802
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 1000 --num_iter 5 --weather "uniform" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|█████████████████▌ | 2/5 [05:12<07:44, 154.82s/it]Average violation so far: 0.03
60%|██████████████████████████▍ | 3/5 [07:47<05:09, 155.00s/it]Average violation so far: 0.03133333333333333
80%|███████████████████████████████████▏ | 4/5 [10:11<02:31, 151.70s/it]Average violation so far: 0.031
100%|████████████████████████████████████████████| 5/5 [12:44<00:00, 152.82s/it]
Average solver time : 152.8250288963318
printing specification violation array
[0.034, 0.026, 0.034, 0.03, 0.03]
probability of violating the spec less then the threshold for each iter:
[2.52242671e-12 4.66293670e-15 2.52242671e-12 1.18571819e-13
1.18571819e-13]
average alpha_nu values of the array:
1.0573319997320141e-12
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 1000 --num_iter 5 --weather "x-neg-bias" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|█████████████████▌ | 2/5 [05:07<07:36, 152.08s/it]Average violation so far: 0.20600000000000002
60%|██████████████████████████▍ | 3/5 [08:34<05:37, 168.62s/it]Average violation so far: 0.211
80%|███████████████████████████████████▏ | 4/5 [12:40<03:11, 191.86s/it]Average violation so far: 0.21225
100%|████████████████████████████████████████████| 5/5 [16:35<00:00, 199.19s/it]
Average solver time : 199.19611296653747
printing specification violation array
[0.202, 0.21, 0.221, 0.216, 0.22]
probability of violating the spec less then the threshold for each iter:
[1. 1. 1. 1. 1.]
average alpha_nu values of the array:
1.0
(stormpyenv-upt) mcubuktepe@asg-a32265:~/research/codes/scenario_code/drone$ python3 samplingprism_drone_scenario.py --model-file drone_par_param_effect_example.nm --property-file drone_par.prctl --threshold 0.95 --bisimulation none --direction geq --num_samples 1000 --num_iter 5 --weather "y-pos-bias" --violation_prob 0.90
Building model from drone_par_param_effect_example.nm
ModelType.MDP
Model supports parameters: True
Number of states before bisim: 7679
Number of params before bisim: 2756
Number of states after bisim: 7679
Number of params after bisim: 2756
40%|████████████████████████████████████████████████████████ | 2/5 [06:12<09:48, 196.22s/it]Average violation so far: 0.192
60%|████████████████████████████████████████████████████████████████████████████████████ | 3/5 [08:51<06:10, 185.10s/it]Average violation so far: 0.19299999999999998
80%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 4/5 [11:30<02:57, 177.17s/it]Average violation so far: 0.19649999999999998
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [14:06<00:00, 169.31s/it]
Average solver time : 169.31317949295044
printing specification violation array
[0.196, 0.188, 0.195, 0.207, 0.183]
probability of violating the spec less then the threshold for each iter:
[1. 1. 1. 1. 1.]
average alpha_nu values of the array:
1.0
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -15,6 +15,8 @@ import create_prism
from gurobipy import *
from ..functions import BCDF, CDF, etaLow, betaLow, compute_eta, compute_beta, \
ticDiff, tocDiff
import argparse
......@@ -32,7 +34,9 @@ parser.add_argument('--threshold', action="store", dest='threshold', default=0)
parser.add_argument('--num_samples', action="store", dest='num_samples', default=0)
parser.add_argument('--num_iter', action="store", dest='num_iter', default=0)
parser.add_argument('--direction', action="store", dest='direction', default=0)
parser.add_argument('--violation_prob', action="store", dest='violation_prob', default=0)
parser.add_argument('--eta', action="store", dest='eta', default=0)
parser.add_argument('--beta', action="store", dest='beta', default=0)
parser.add_argument('--weather', action="store", dest='weather', default=0)
# Now, parse the command line arguments and store the
# values in the `args` variable
......@@ -46,8 +50,16 @@ bisimulation_type=args.bisim
numiter=int(args.num_iter)
numsample=int(args.num_samples)
threshold=float(args.threshold)
violation_prob=float(args.violation_prob)
eta=float(args.eta)
beta=args.beta
direction_type=args.direction
weather = args.weather
if beta != 0:
try:
beta = [float(beta)]
except:
beta = list(literal_eval(beta))
if direction_type=="leq":
direction=True
......@@ -76,7 +88,7 @@ def search(values, searchFor):
# paraminit = [[0 for state in range(numstate)] for state in range(numstate)]
def run_sample(numiter,numsample,thres,direction,parameters,model,properties,model_file):
def run_sample(numiter,numsample,thres,direction,parameters,model,properties,model_file,weather):
'''
:param numiter: number of trials to compute the number of samples that satisfies the probability
:param numsample: number of sampled pMDPs/pMCs to check
......@@ -94,19 +106,16 @@ def run_sample(numiter,numsample,thres,direction,parameters,model,properties,mod
raise RuntimeError("Invalid model type, model should be pDTMC or pMDP.")
groups = create_prism.parameter_groups()
print(groups)
#print(groups)
#instantiate parameters
parameter_defs = create_prism.parameter_definitions(groups,False)
parameter_instantiations=create_prism.parameter_dirichlet_instantiations(groups)
start3 = time.time()
#for each iteration
for iter in tqdm(range(numiter)):
counter=0
# for each sample
for i in range(int(numsample)):
parameter_instantiations=create_prism.parameter_dirichlet_instantiations(groups,weather)
parameter_samples=create_prism.parameter_dirichlet_samples(parameter_instantiations)
# print(parameter_samples)
point = dict()
......@@ -130,7 +139,9 @@ def run_sample(numiter,numsample,thres,direction,parameters,model,properties,mod
[[x, stormpy.RationalRF(val)] for x, val in point.items()])
instantiated_model = instantiator.instantiate(rational_parameter_assignments)
result = stormpy.model_checking(instantiated_model, properties[0]).at(instantiated_model.initial_states[0])
#print(result)
# if result <1e-6 or result > 1-1e-6:
# print(result)
# print(result)
#append the counter according to the spec
if direction==True:
if float(result)>thres:
......@@ -139,7 +150,7 @@ def run_sample(numiter,numsample,thres,direction,parameters,model,properties,mod
if float(result)<thres:
counter=counter+1
counterarray[iter]=counter/((numsample)*1.0)
counterarray[iter]= counter / numsample
if iter>1:
#print(counterarray[0:iter])
......@@ -192,7 +203,7 @@ def load_problem(model_file, property_file,bisimulation_type):
return parameters,model,properties
def compute_avg_satprob(counterarray,N,eps,flag):
def compute_avg_beta(counterarray,N,eta):
'''
:param counterarray: approximate satisfaction probs
:param N: number of samples
......@@ -206,34 +217,110 @@ def compute_avg_satprob(counterarray,N,eps,flag):
for iters in range(len(counterarray)):
#compute number of constraints to remove in the LP
if flag:
removeconstraints = int(N * (counterarray[iters]))
else:
removeconstraints = int(N * (1 - counterarray[iters]))
removeconstraints = int(N * (counterarray[iters]))
# print(N,removeconstraints,eps)
#approximately compute the satisfaction prob
val3 = BCDF(eps, N, removeconstraints)
valuearray[iters] = val3
if flag:
print("probability of violating the spec less then the threshold for each iter:")
print(1-valuearray)
print("average alpha_nu values of the array:")
print(1-np.mean(valuearray))
else:
print("probability of violating the spec less then the threshold for each iter:")
print(valuearray)
print("average alpha_nu values of the array:")
print(np.mean(valuearray))