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

Added crows 15,7 benchmark

parent 6b69f02e
......@@ -157,14 +157,17 @@ def BCDF(p, n, m):
def CDF(x):
return (1.0 + math.erf(x/math.sqrt(2.0)))/2.0
def determine_discarded(N=1000, beta=1e-06, eta=0.86):
def determine_discarded(N=1000, beta=1e-6, eta=0.86):
k = 0
res = -1
while res < beta and k < N:
res = sum([math.comb(N, i)*(eta)**(N-i)*eta**(i) for i in range(k+1)])
#res = math.comb(N+1, N) * sum([math.comb(N, i)*(eta)**(N-i)*(1-eta)**(i) for i in range(k+1)])
res = sum([math.comb(N, i)*(eta)**(N-i)*(1-eta)**(i) for i in range(k+1)])
k+=1
print(k,':',res)
\ No newline at end of file
# print(k,':',res)
print('Number of discarded constraints:',k)
\ No newline at end of file
import itertools
import math
import numpy as np
template_header = """
mdp
const int MAXX = 10;
const int MAXY = 10;
const int MAXZ = 10;
formula valid = (1 < x & x < MAXX-1 & 1 < y & y < MAXY-1 & 1 < z & z < MAXZ-1) & ok=1;
formula attarget = x > MAXX - 3;
"""
control_module = """
module control
c : [0..3] 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);
endmodule
"""
mod_head = """
module mod
ok : [0..1] init 1;
cond : [0..3] init 0;
wdir : [0..7] init 0;
x : [0..MAXX] init 2;
y : [0..MAXY] init floor(MAXY/2);
z : [0..MAXZ] init floor(MAXZ/2);
phase : [0..1] init 0;
[move] x < MAXX -> 1: (x'=x+1);
[move] x > 0 -> 1: (x'=x-1);
[move] y < MAXY -> 1: (y'=y+1);
[move] y > 0 -> 1: (y'=y-1);
[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);
"""
mod_foot="""
endmodule
"""
labels="""
"""
#import itertools
#for element in itertools.product(*somelists):
# print(element)
xsplits = [4,6]
ysplits = [7]
zsplits = [5]
x_area_splits = [0] + xsplits + ["MAXX+1"]
y_area_splits = [0] + ysplits + ["MAXY+1"]
z_area_splits = [0] + zsplits + ["MAXZ+1"]
def _areadefinitions(dir, splits):
return "\n".join(["formula {}area{} = {} <= {} & {} < {};".format(dir, i, l, dir, dir, h) for i, (l, h) in
enumerate(zip(splits[:-1], splits[1:]))])
x_area_definitions = _areadefinitions("x", x_area_splits)
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]
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:
# straight
affected_var, other_var = ("x","y") if dir in [0,4] else ("y","x")
main_var_sym = "+" if dir < 4 else "-"
sec_var_op = (["","-1","+1"][eff%3]) if dir in [0,6] else (["","+1","-1"][eff%3])
return "({}' = {}{}{})".format(affected_var, affected_var, main_var_sym, math.ceil(float(eff)/3)) + " & " + "({}' = {}{})".format(other_var,other_var,sec_var_op)
else:
# diagionals
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,int(eff/3),ysym,int(eff/3))
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+2) & (y'=y+1)")
if dir == 3:
return ("(y'=y+2) & (x'=x-1)")
if dir == 5:
return ("(x'=x-2) & (y'=y-1)")
if dir == 7:
return ("(y'=y-2) & (x'=x+1)")
if eff == 5:
if dir == 1:
return ("(y'=y+2) & (x'=x+1)")
if dir == 3:
return ("(x'=x-2) & (y'=y+1)")
if dir == 5:
return ("(y'=y-2) & (x'=x-1)")
if dir == 7:
return ("(x'=x+2) & (y'=y-1)")
assert False
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)]
weighted_updates = ["{}: {}".format(p, u) for p, u in zip(probs, updates)]
update_str = " + ".join(weighted_updates)
return "\t[env] {} -> {};".format(guard, update_str)
commands = [create_command(c,w,x,y,z) for (c,w,x,y,z) in itertools.product(conditions, winddirs, xzones, yzones, zzones)]
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 (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)])
return groups
def parameter_definitions(groups,init):
if init:
return "\n".join(["const double {} = 1/{};".format(par, len(group)) for group in groups for par in group])
else:
return "\n".join(["const double {};".format(par) for group in groups for par in group])
def parameter_dirichlet_instantiations(groups):
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):
array = np.random.random_sample((len(group),))
array[0]=2*array[0]
array[1]=1*array[1]
array[2]=4*array[2]
array[3]=10*array[3]
array[4]=5*array[4]
array[5]=3*array[5]
array[6]=3*array[6]
#create uniformly
else:
array = np.random.random_sample((len(group),))
instantiations[tuple(group)]=array
return instantiations
#for each group, gather the dirichlet weights and sample from the weights
def parameter_dirichlet_samples(parameter_instantiations):
instantiations_sample=dict()
for group in parameter_instantiations:
#print(group)
array = parameter_instantiations[tuple(group)]
sample=np.random.dirichlet((array), 1)
#print(sample,array)
instantiations_sample[tuple(group)]=sample
#print(array)
return instantiations_sample
#groups = parameter_groups()
#parameter_defs = parameter_definitions(groups,False)
#parameter_instantiations=parameter_dirichlet_instantiations(groups)
#parameter_samples=parameter_dirichlet_samples(parameter_instantiations)
#for element in parameter_samples:
#print(np.sum(parameter_samples[tuple(element)]))
#print(parameter_instantiations)
#print(parameter_defs)
with open("drone_par.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)
// Exported by storm
// Original model type: DTMC
@type: DTMC
@parameters
p1_0 p2_0 p3_0 p4_0 p5_0 p8_0 p9_0 p11_0 r2_0 r1_0 p3_1 p3_2 p3_3
@reward_models
@nr_states
30
@model
state 0 [0] goal
action 0 [0]
0 : 1
state 1 [r1_0]
action 0 [0]
4 : p3_0+p3_3
12 : p3_1+p3_2
state 2 [0]
action 0 [0]
14 : p5_0
28 : (-1)*p5_0+1
state 3 [1]
action 0 [0]
2 : 1
state 4 [r2_0]
action 0 [0]
16 : (-1)*p2_0+1
25 : p2_0
state 5 [0] init
action 0 [0]
1 : 1/13
2 : 1/13
3 : 1/13
4 : 1/13
6 : 1/13
8 : 1/13
13 : 1/13
18 : 1/13
20 : 1/13
23 : 1/13
24 : 1/13
25 : 1/13
27 : 1/13
state 6 [r2_0]
action 0 [0]
1 : p2_0
17 : (-1)*p2_0+1
state 7 [1-r2_0]
action 0 [0]
1 : p11_0
13 : (-1)*p11_0+1
state 8 [0]
action 0 [0]
7 : (-1)*p5_0+1
10 : p5_0
state 9 [1]
action 0 [0]
0 : (-1)*p11_0+1
8 : p11_0
state 10 [1]
action 0 [0]
8 : 1
state 11 [1]
action 0 [0]
6 : p8_0
8 : (-1)*p8_0+1
state 12 [1-r2_0]
action 0 [0]
6 : p9_0
8 : (-1)*p9_0+1
state 13 [0]
action 0 [0]
9 : (-1)*p5_0+1
14 : p5_0
state 14 [1]
action 0 [0]
13 : 1
state 15 [r2_0]
action 0 [0]
4 : p2_0
11 : (-1)*p2_0+1
state 16 [1-r1_0]
action 0 [0]
4 : (-1)*p8_0+1
15 : p8_0
state 17 [1]
action 0 [0]
6 : (-1)*p8_0+1
18 : p8_0
state 18 [1]
action 0 [0]
6 : p1_0
20 : (-1)*p1_0+1
state 19 [1]
action 0 [0]
18 : p11_0
23 : (-1)*p11_0+1
state 20 [0]
action 0 [0]
19 : (-1)*p5_0+1
22 : p5_0
state 21 [1]
action 0 [0]
20 : p11_0
24 : (-1)*p11_0+1
state 22 [1]
action 0 [0]
20 : 1
state 23 [0]
action 0 [0]
21 : (-1)*p5_0+1
24 : p5_0
state 24 [1]
action 0 [0]
23 : 1
state 25 [1]
action 0 [0]
4 : p4_0
27 : (-1)*p4_0+1
state 26 [1]
action 0 [0]
2 : (-1)*p11_0+1
25 : p11_0
state 27 [0]
action 0 [0]
26 : (-1)*p5_0+1
29 : p5_0
state 28 [1]
action 0 [0]
3 : (-1)*p11_0+1
27 : p11_0
state 29 [1]
action 0 [0]
27 : 1
import itertools
import math
import numpy as np
template_header = """
mdp
const int MAXX = 10;
const int MAXY = 10;
const int MAXZ = 10;
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;
"""
control_module = """
module control
c : [0..3] 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);
endmodule
"""
mod_head = """
module mod
ok : [0..1] init 1;
cond : [0..3] init 0;
wdir : [0..7] init 0;
x : [0..MAXX] init 2;
y : [0..MAXY] init floor(MAXY/2);
z : [0..MAXZ] init floor(MAXZ/2);
phase : [0..1] init 0;
[move] x < MAXX -> 1: (x'=x+1);
[move] x > 0 -> 1: (x'=x-1);
[move] y < MAXY -> 1: (y'=y+1);
[move] y > 0 -> 1: (y'=y-1);
[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);
"""
mod_foot="""
endmodule
"""
reward_head="""
rewards "time" // flight time
"""
reward_foot="""
endrewards
"""
labels="""
"""
#import itertools
#for element in itertools.product(*somelists):
# print(element)
xsplits = [4,6]
ysplits = [7]
zsplits = [5]
x_area_splits = [0] + xsplits + ["MAXX"]
y_area_splits = [0] + ysplits + ["MAXY"]
z_area_splits = [0] + zsplits + ["MAXZ"]
def _areadefinitions(dir, splits):
return "\n".join(["formula {}area{} = {} <= {} & {} < {};".format(dir, i, l, dir, dir, h) for i, (l, h) in
enumerate(zip(splits[:-1], splits[1:]))])
x_area_definitions = _areadefinitions("x", x_area_splits)
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]
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:
# straight
affected_var, other_var = ("x","y") if dir in [0,4] else ("y","x")
main_var_sym = "+" if dir < 4 else "-"
sec_var_op = (["","-1","+1"][eff%3]) if dir in [0,6] else (["","+1","-1"][eff%3])
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]
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 == 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)")
assert False
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)]
weighted_updates = ["{}: {}".format(p, u) for p, u in zip(probs, updates)]
update_str = " + ".join(weighted_updates)
return "\t[env] {} -> {};".format(guard, update_str)
commands = [create_command(c,w,x,y,z) for (c,w,x,y,z) in itertools.product(conditions, winddirs, xzones, yzones, zzones)]
def create_reward(w,x,y,z):
guard = "wdir = {} & xarea{} & yarea{} & zarea{} & valid".format(w, x, y, z)
#print(guard)
#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)]
rewards = "RewardEnvw{}x{}y{}z{}".format(w, x, y, z,)
#print(rewards)
#weighted_updates = ["{}".format(p,) for p in zip(rewards)]
#update_str = " + ".join(weighted_updates)
return "\t[env] true & {}: {};".format(guard, rewards)
rewards = [create_reward(w,x,y,z) for (w,x,y,z) in itertools.product(winddirs, xzones, yzones, zzones)]
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 (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)])
#for r in range(1):
groups.append(["RewardEnv","RewardChd"])
#for (w) in (winddirs):
# groups.append(["RewardEnvw{}x{}y{}z{}".format(w,x,y,z) for x,y,z in itertools.product(xzones,yzones,zzones)])
#print(groups)
return groups