Commit 6b69f02e authored by Thom Badings's avatar Thom Badings
Browse files

Synch

parent c65386f8
### BRP (256,5)
python3 samplingprism.py --model-file brp/brp_256_5.pm --property-file brp/brp.prctl --threshold 0.5 --bisimulation weak --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.055 --outfile brp_256_5_sat.txt
python3 samplingprism.py --model-file brp/brp_256_5.pm --property-file brp/brp.prctl --threshold 0.5 --bisimulation weak --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.055 --outfile brp_256_5_sat.txt
python3 samplingprism.py --model-file brp/brp_256_5.pm --property-file brp/brp.prctl --threshold 0.5 --bisimulation weak --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.898 --outfile brp_256_5_unsat.txt
python3 samplingprism.py --model-file brp/brp_256_5.pm --property-file brp/brp.prctl --threshold 0.5 --bisimulation weak --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.898 --outfile brp_256_5_unsat.txt
### BRP (16,5)
python3 samplingprism.py --model-file brp/brp_rewards4_16_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation strong --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.275 --outfile brp_16-5_sat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_16_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation strong --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.275 --outfile brp_16-5_sat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_16_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation strong --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.676 --outfile brp_16-5_unsat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_16_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation strong --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.676 --outfile brp_16-5_unsat_output.txt
### BRP (32,5)
python3 samplingprism.py --model-file brp/brp_rewards4_32_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation weak --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.232 --outfile brp_32-5_sat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_32_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation weak --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.232 --outfile brp_32-5_sat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_32_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation weak --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 50 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.718 --outfile brp_32-5_unsat_output.txt
python3 samplingprism.py --model-file brp/brp_rewards4_32_5.pm --property-file brp/brp_rewards4.prctl --threshold 3 --bisimulation weak --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.718 --outfile brp_32-5_unsat_output.txt
### CROWDS (10,5)
......@@ -23,9 +22,9 @@ python3 samplingprism.py --model-file crowds/crowds_10_5.pm --property-file crow
### CROWDS (20,7)
python3 samplingprism.py --model-file crowds/crowds_20_7.pm --property-file crowds/crowds.prctl --threshold 0.9 --bisimulation weak --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.416 --outfile crowds_20-7_sat_output.txt
python3 samplingprism.py --model-file crowds/crowds_20_7.pm --property-file crowds/crowds.prctl --threshold 0.9 --bisimulation strong --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.416 --outfile crowds_20-7_sat_output.txt --verbose 1
python3 samplingprism.py --model-file crowds/crowds_20_7.pm --property-file crowds/crowds.prctl --threshold 0.9 --bisimulation weak --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.534 --outfile crowds_20-7_unsat_output.txt
python3 samplingprism.py --model-file crowds/crowds_20_7.pm --property-file crowds/crowds.prctl --threshold 0.9 --bisimulation strong --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.534 --outfile crowds_20-7_unsat_output.txt
### NAND (10,5)
......@@ -35,9 +34,9 @@ python3 samplingprism.py --model-file nand/nand_10_5.pm --property-file nand/nan
### NAND (25,5)
python3 samplingprism.py --model-file nand/nand_25_5.pm --property-file nand/nand.prctl --threshold 0.05 --bisimulation strong --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.206 --outfile nand_25-5_sat_output.txt
python3 samplingprism.py --model-file nand/nand_25_5.pm --property-file nand/nand.prctl --threshold 0.05 --bisimulation strong --direction geq --num_samples '(100,1000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.206 --outfile nand_25-5_sat_output.txt
python3 samplingprism.py --model-file nand/nand_25_5.pm --property-file nand/nand.prctl --threshold 0.05 --bisimulation strong --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.744 --outfile nand_25-5_unsat_output.txt
python3 samplingprism.py --model-file nand/nand_25_5.pm --property-file nand/nand.prctl --threshold 0.05 --bisimulation strong --direction leq --num_samples '(100,1000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.744 --outfile nand_25-5_unsat_output.txt
### CONSENSUS (2,2)
......@@ -47,6 +46,10 @@ python3 samplingprism.py --model-file consensus/coin2.pm --property-file consens
### CONSENSUS (4,2)
python3 samplingprism.py --model-file consensus/coin4.pm --property-file consensus/coin2.prctl --threshold 0.25 --bisimulation strong --direction geq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.063 --outfile consensus_4-2_sat_output.txt
python3 samplingprism.py --model-file consensus/coin4.pm --property-file consensus/coin2.prctl --threshold 0.25 --bisimulation strong --direction geq --num_samples '(100,1000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.063 --outfile consensus_4-2_sat_output.txt
python3 samplingprism.py --model-file consensus/coin4.pm --property-file consensus/coin2.prctl --threshold 0.25 --bisimulation strong --direction leq --num_samples '(100,1000,10000)' --num_iter 5 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.888 --outfile consensus_4-2_unsat_output.txt
### DRONE
python3 samplingprism.py --model-file consensus/coin4.pm --property-file consensus/coin2.prctl --threshold 0.25 --bisimulation strong --direction leq --num_samples '(100,500,1000,5000,10000)' --num_iter 10 --beta '(0.9,0.99,0.999,0.9999)' --eta 0.888 --outfile consensus_4-2_unsat_output.txt
python3 samplingprism.py --model-file drone/drone_par.nm --property-file drone/drone.prctl --threshold 0.9 --bisimulation none --direction geq --num_samples '1000' --num_iter 1 --beta '0.999999' --outfile drone_output.txt
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
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]=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
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:
array = parameter_instantiations[tuple(group)]
if any("RewardEnvw" in group_item for group_item in group):
#print(group)
sample=10*np.random.dirichlet((array), 1)
#print(sample)
else:
sample=np.random.dirichlet((array), 1)
#print(sample,array)
instantiations_sample[tuple(group)]=sample
#print(array)
#print(sample,group)
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_rewardparam.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)
Pmax=? [F attarget ]
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
from __future__ import division
import stormpy
import stormpy.core
import stormpy.logic
import stormpy.pars
import re
import stormpy.examples
import stormpy.examples.files
import stormpy.simulator
import numpy as np
import cvxpy
# import time
import numpy as np
# import time
# from tqdm import tqdm
# from gurobipy import *
# from collections import defaultdict
# import sobol_lib
import create_prism
def find_rew0_states(property, model):
formula = property.raw_formula
assert type(formula) == stormpy.logic.RewardOperator
path_formula = formula.subformula
if type(path_formula) == stormpy.logic.EventuallyFormula:
psi_formula = path_formula.subformula
else:
raise ValueError("Property type not supported")
psi_result = stormpy.model_checking(model, psi_formula)
psi_states = psi_result.get_truth_values()
return psi_states
def get_prob01States(model, formulas):
parameters = model.collect_probability_parameters()
instantiator = stormpy.pars.PMdpInstantiator(model)
print(parameters)
point = dict()
for p in parameters:
point[p] = stormpy.RationalRF(0.4)
print(p)
instantiated_model = instantiator.instantiate(point)
assert instantiated_model.nr_states == model.nr_states
assert not instantiated_model.has_parameters
pathform = formulas[0].raw_formula.subformula
assert type(pathform) == stormpy.logic.EventuallyFormula
labelform = pathform.subformula
labelprop = stormpy.core.Property("label-prop", labelform)
phiStates = stormpy.BitVector(instantiated_model.nr_states, True)
psiStates = stormpy.model_checking(instantiated_model, labelprop).get_truth_values()
(prob0, prob1) = stormpy.compute_prob01_states(model, phiStates, psiStates)
return prob0, prob1
def generate_sim_dict(path):
"""
:param path: path of the exported file
:return a dictionary that contains the label of each state:
"""
file1 = open(path, 'r')
Lines = file1.readlines()
count = 0
# Strips the newline character
stringflag = False
stringdict = dict()
#go through the file
for line in Lines:
count += 1
stringval = str(line.strip())
# print(stringval)
#this is a bit harcoded way to get the labels from th efile format
if len(stringval) > 1:
if stringflag:
stringdict[stateval] = str(line.strip())
stringflag = False
if stringval[0] == 's':
x = stringval.split(" ")
stateval = int(x[1])
# print("Line{}: {}".format(count, line.strip()),stateval)
stringflag = True
# print(stringdict)
return (stringdict)
def replaceAll(file,searchExp,replaceExp):
for line in fileinput.input(file, inplace=1):
if searchExp in line:
line = replaceExp
sys.stdout.write(line)
def simulator_stormpy(model,dict_labels,action_dict):
"""
:param model: stormpy model of the satellite problem
:param dict_labels: dictionary that contains the labels for each state
:param action_dict: dictionary for policy
:return:
path_drine: 3d positions of the drone path
"""
#set up the simulator
simulator = stormpy.simulator.create_simulator(model)
#list for trajectory for the states
paths= []
#list for the drone trajectory
path_drone = []
#start the simulator
state, reward = simulator.restart()
path = [f"{state}"]
#get the integer values in the label, hard coded
intvals2 = re.findall(r'\d+', dict_labels[int(state)])
#append the trajectory
path_drone.append([int(intvals2[4]), int(intvals2[5]), int(intvals2[6])])
#simulator loop
for n in range(2000):
#get the action
actions = simulator.available_actions()
# print(dict_labels[int(state)])
#get the label of the previous state
intvals2_prev = re.findall(r'\d+', dict_labels[int(state)])
#select action
select_action = action_dict[state]
#append to the state-action path list
path.append(f"--act={actions[select_action]}-->")
#step through the simulatorr
state, reward = simulator.step(actions[select_action])
# print(state)
path.append(f"{state}")
# print(dict_labels[int(state)])
#get the label of the next state
intvals2 = re.findall(r'\d+', dict_labels[int(state)])
#print(dict_labels[int(state)])
#appends the path of the drone
path_drone.append([int(intvals2[4]), int(intvals2[5]), int(intvals2[6])])
#break if simulator is done
if simulator.is_done():
break
return path_drone
def search(values, searchFor):
for k in values:
# for v in values[k]:
if searchFor in k:
return k
return None
if __name__ == "__main__":
threshold = 0.30
direction = "above" # can be "below" or "above"
# Reading file
path = "drone_par.nm"
print("Building model from {}".format(path))
program = stormpy.parse_prism_program(path)
# A spec that I made up in drone_par.nm file
formula_str = 'Pmax=? [F attarget ]'
print(formula_str)
properties = stormpy.parse_properties_for_prism_program(formula_str, program, None)
#build options
options = stormpy.BuilderOptions([p.raw_formula for p in properties])
options.set_build_state_valuations(True)
options.set_build_choice_labels(True)
options.set_build_with_choice_origins(True)
options.set_build_all_labels(True)
options.set_build_all_reward_models(True)
#build the model
model = stormpy.build_sparse_parametric_model_with_options(program, options)
numstate = model.nr_states
print("Number of states before bisim:", numstate)
path_pmc = path + str("_prism_export")
stormpy.export_to_drn(model, path_pmc)
numstate=model.nr_states
#print(numstate)
prob0E, prob1A = stormpy.prob01min_states(model, properties[0].raw_formula.subformula)
#print(prob0E,prob1A)
#print(model.initial_states)
parameters = model.collect_probability_parameters()
instantiator = stormpy.pars.PMdpInstantiator(model)
groups = create_prism.parameter_groups()
#print(groups)
#instantiate parameters
parameter_defs = create_prism.parameter_definitions(groups,False)
parameter_instantiations=create_prism.parameter_dirichlet_instantiations(groups)
#sample parameters from the dirichlet distribution
for i in range(10):
parameter_samples=create_prism.parameter_dirichlet_samples(parameter_instantiations)
#print(parameter_samples)
point = dict()
for x in parameters:
#instantiate parameters
parameter_group = search(parameter_instantiations, str(x.name))
element_int=0
for element in parameter_group:
parameter_sample_array = parameter_samples[tuple(parameter_group)]
#print(parameter_sample_array)
if (str(element)) == (str(x.name)):
point[x] = parameter_sample_array[0,element_int]
#print(point[x])
element_int=element_int+1
#model check the instantiated MDP
rational_parameter_assignments = dict(
[[x, stormpy.RationalRF(val)] for x, val in point.items()])
instantiated_model = instantiator.instantiate(rational_parameter_assignments)
mc_res = stormpy.model_checking(instantiated_model, properties[0]).at(model.initial_states[0])
print("Probability of satisfying the spec:",mc_res)
result = stormpy.model_checking(instantiated_model, properties[0], extract_scheduler=True)
scheduler = result.scheduler
assert scheduler.memoryless
assert scheduler.deterministic
# load the policy to a dict
action_dict = dict()
for state in model.states:
choice = scheduler.get_choice(state)
action = choice.get_deterministic_choice()
action_dict[state] = action
action_dict[int(state)] = action
# genetate the label
dict_labels = generate_sim_dict(path_pmc)