Skip to content
Snippets Groups Projects
Commit 7b6c1f44 authored by Vincent Grande's avatar Vincent Grande
Browse files

minor updates to the markdown

parent 988f5605
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ExperimentsBlackBoxStorm.ipynb
+ 4
− 2
View file @ 7b6c1f44
%% Cell type:code id: tags:
``` python
import stormpy
import stormpy.examples
import stormpy.examples.files
import stormpy.simulator
import numpy as np
import scipy
import scipy.sparse
import multiprocess
#import time
import os
```
%% Cell type:markdown id: tags:
## Synthetic Markov Chains
Here we create the two synthetic Markov Chains used in the paper in the Storm format.
%% Cell type:code id: tags:
``` python
def createSimpleMC(R, p, Pgood=0.5):
num_states = R+3
rows = []
cols = []
probs = []
state_labeling = stormpy.storage.StateLabeling(np.int64(num_states))
state_labels = {'golden', 'init'}
for label in state_labels:
state_labeling.add_label(label)
#initial state:
state_labeling.add_label_to_state('init', 0)
rows.append(0)
cols.append(1)
probs.append(1-Pgood)
rows.append(0)
cols.append(2)
probs.append(Pgood)
#Bad state
rows.append(1)
cols.append(1)
probs.append(1)
#normal states
for i in range(2, R):
rows.append(i)
cols.append(i+1)
probs.append(1)
#Decision state
rows.append(R)
cols.append(R+1)
probs.append(1-p)
rows.append(R)
cols.append(R+2)
probs.append(p)
#Last normal state
rows.append(R+1)
cols.append(2)
probs.append(1)
#Golden state
state_labeling.add_label_to_state('golden', np.int64(R)+2)
rows.append(R+2)
cols.append(2)
probs.append(1)
transition_matrix = scipy.sparse.csc_matrix(
(probs, (rows, cols)), shape=(num_states, num_states)).toarray()
transition_matrix = stormpy.build_sparse_matrix(transition_matrix)
components = stormpy.SparseModelComponents(
transition_matrix=transition_matrix, state_labeling=state_labeling)
dtmc = stormpy.storage.SparseDtmc(components)
return dtmc
def createPGoodMC(k, pgamma):
num_states = k+2
rows = []
cols = []
probs = []
state_labeling = stormpy.storage.StateLabeling(np.int64(num_states))
state_labels = {'golden', 'init'}
for label in state_labels:
state_labeling.add_label(label)
#initial state:
state_labeling.add_label_to_state('init', 1)
#Bad state
rows.append(0)
cols.append(0)
probs.append(1)
#normal states
for i in range(1, k+1):
state_labeling.add_label_to_state('golden', i)
rows.append(i)
cols.append(i+1)
probs.append(pgamma)
rows.append(i)
cols.append(0)
probs.append(1-pgamma)
#goal state
rows.append(k+1)
cols.append(k+1)
probs.append(1)
state_labeling.add_label_to_state('golden', k+1)
transition_matrix = scipy.sparse.csc_matrix(
(probs, (rows, cols)), shape=(num_states, num_states)).toarray()
transition_matrix = stormpy.build_sparse_matrix(transition_matrix)
components = stormpy.SparseModelComponents(
transition_matrix=transition_matrix, state_labeling=state_labeling)
dtmc = stormpy.storage.SparseDtmc(components)
return dtmc
```
%% Cell type:markdown id: tags:
## The Code
We call the *test_multiple_strategies_multiprocess* function, to initiate the experiments. It then handles the test_parameters, initiating *do_one_testrun* for each repetition of each set of parameters, and then collects the data in the end.
*do_one_testrun* initiates the Storm Simulator based on the test_parameters, and then calls *c_strategy_whole_run*. *c_strategy_whole_run* keeps track of the restarts and steps, and calls repeatedly *c_strategy_fragment*, until a call of *c_strategy_fragment* results in reaching the step threshold.
*c_strategy_fragment* implements one fragment of the strategy, e.g. what happens between two restarts, or after the last restart. It checks for the occurence of the Rabin labels, and executes entire blocks of steps with *simulator_step_block*.
%% Cell type:code id: tags:
``` python
#Does 2*num_restart^c steps and checks for the occurence of Rabin Labels
def simulator_step_block(simulator,num_restart,rabin_pairs,c):
num_pairs=len(rabin_pairs)
golden_state = [False for i in range(num_pairs)]
bad_state = [False for i in range(num_pairs)]
for i in range(np.int32(num_restart**c)):
cur_labels=simulator.step()[2]
for i in range(num_pairs):
if (rabin_pairs[i][0][0] in cur_labels) == rabin_pairs[i][0][1]:
golden_state[i]=True
if (rabin_pairs[i][1][0] in cur_labels) == rabin_pairs[i][1][1]:
bad_state[i]=True
return golden_state, bad_state
#Implements a Rabin automaton for the condition needed in the not_gridworld example
def simulator_step_block_not_gridworld(simulator, num_restart, rabin_pairs, rabin_state, c):
num_pairs = len(rabin_pairs)
golden_state = [False for i in range(num_pairs)]
bad_state = [False for i in range(num_pairs)]
for i in range(num_restart**c):
cur_labels = simulator.step()[2]
for i in range(num_pairs):
if (rabin_pairs[i][0][0] in cur_labels) == rabin_pairs[i][0][1] and rabin_state != i:
golden_state[i] = True
rabin_state = i
if (rabin_pairs[i][1][0] in cur_labels) == rabin_pairs[i][1][1]:
bad_state[i] = True
return golden_state, bad_state, rabin_state
#Checks whether there is a Rabin pair such that the second half of the run is good
def check_last_half(last_golden_block, last_bad_block, num_blocks, num_pairs):
half_time = num_blocks//2
for i in range(num_pairs):
if last_golden_block[i]>half_time and last_bad_block[i]<=half_time:
return True
return False
def c_strategy_fragment(simulator,num_restart,c,rabin_pairs,threshold,mode):
simulator.restart()
num_blocks=0
rabin_state = 10 #Anything will work. Keeps track of the state of the Rabin automaton for the not_gridworld case
num_pairs=len(rabin_pairs)
last_golden_block = [1 for i in range(num_pairs)] #Initialised to 1 so that the while loop condition is fulfilled for num_blocks = 0
last_bad_block = [0 for i in range(num_pairs)]
if threshold<10*num_restart**c: #If the step-threshold is too low, manually set it to a higher value
threshold=np.int32(10*num_restart**c)
while check_last_half(last_golden_block, last_bad_block, num_blocks, num_pairs):
for i in range(2):
num_blocks += 1
if mode == 'not_gridworld' or mode == 'not_gridworld_Big':
golden_state, bad_state, rabin_state = simulator_step_block_not_gridworld(simulator, num_restart, rabin_pairs, rabin_state, c)
else:
golden_state, bad_state = simulator_step_block(
simulator, num_restart, rabin_pairs,c)
for i in range(num_pairs):
if golden_state[i]:
last_golden_block[i] = num_blocks
if bad_state[i]:
last_bad_block[i] = num_blocks
if num_blocks*np.int32(num_restart**c) > threshold: #Threshold is reached and run is assumed to continue in a good run
return num_blocks*np.int32(num_restart**c), 'threshold'
if num_blocks==2:
return num_blocks*np.int32(num_restart **c), 'after initialisation'
return num_blocks*np.int32(num_restart**c), 'No valid Rabin pair'
#Repeatedly calls c_fragment... until the threshold is reached for a run.
def c_strategy_whole_run(simulator,c,rabin_pairs, threshold_fragment, threshold_restarts, mode):#=1000000,threshold_restarts=10000, mode):
num_restarts=1
total_steps=0
for i in range(threshold_restarts):
cur_steps,outcome = c_strategy_fragment(simulator, num_restarts,c, rabin_pairs, threshold_fragment, mode)
if outcome == "threshold":
return total_steps,num_restarts,True
total_steps=total_steps+cur_steps
num_restarts=num_restarts+1
return total_steps,num_restarts-1,False
def do_one_testrun(test_parameters):
os.nice(10)
[test_parameters, mode] = test_parameters
print('Began with:'+str(test_parameters))
#Initialise the simulator in the correct way depending on the mode
if mode == 'syntheticPgood':
[c, k, Pgamma, threshold_fragment, threshold_restarts] = test_parameters
k = np.int64(k)
testdtmc = createPGoodMC(k, Pgamma)
if mode == 'syntheticPmRm':
[c, P, R, Pgood, threshold_fragment, threshold_restarts] = test_parameters
R = np.int64(R)
testdtmc = createSimpleMC(R, P, Pgood)
if mode == 'syntheticPmRm' or mode == 'syntheticPgood':
simulator = stormpy.simulator.create_simulator(testdtmc)
rabin_pairs = [[['golden', True], ['bad', True]]]
else:
[c, path, program, property_file, rabin_pairs,
threshold_fragment, threshold_restarts] = test_parameters
prism_program = stormpy.parse_prism_program(
path+program)
if mode == 'not_gridworld_Big': #In this case, we will switch to a Program Based Simulator, sacrificing runtime to memory efficiency
options = stormpy.BuilderOptions()
options.set_build_all_labels()
simulator = stormpy.simulator.create_simulator(prism_program, options=options)
else:
properties = stormpy.parse_properties(
path+property_file, prism_program)
model = stormpy.build_model(prism_program, properties)
simulator = stormpy.simulator.create_simulator(model)
cur_steps, cur_restarts, success = c_strategy_whole_run(
simulator, c, rabin_pairs, threshold_fragment, threshold_restarts, mode)
print('Done with: mode:'+str(mode)+', '+str(test_parameters)+': steps='+str(cur_steps)+', restarts='+str(cur_restarts)+', success='+str(success))
if not success:
print('ERROR! For c='+str(c)+' and parameters ' +
test_parameters+' Restart threshold was reached.')
raise ValueError('For c='+str(c)+' and program '+test_parameters+' Restart threshold was reached.')
return (cur_steps, cur_restarts)
def test_multiple_strategies_multiprocess(test_parameters, n):
mean_steps = []
mean_restarts = []
long_test_parameters =[]
results_list=[]
k = 0
kmax = len(test_parameters)
print("kmax: "+str(kmax))
for params in test_parameters:
for j in range(n):
long_test_parameters.append((params,))
k+= 1
PROCESSES = 10 #Adjust to desired number of cores
with multiprocess.Pool(PROCESSES) as pool:
results = [pool.map_async(do_one_testrun, p) for p in long_test_parameters]
for r in results:
results_list.append(r.get()[0])
for i in range(len(test_parameters)):
print(i)
cur_total_steps=0
cur_total_restarts=0
for j in range(n):
cur_total_steps += results_list[i*n+j][0]
cur_total_restarts += results_list[i*n+j][1]
print('cur mean steps:'+str(cur_total_steps/n))
mean_steps.append(cur_total_steps/n)
mean_restarts.append(cur_total_restarts/n)
return mean_steps, mean_restarts
```
%% Cell type:code id: tags:
``` python
modes = ['syntheticPgood','syntheticPmRm', 'normal', 'not_gridworld', 'not_gridworld_Big']
```
%% Cell type:markdown id: tags:
## Test parameters
In the following, we will define the test parameters that were used to conduct the experiments. The first entry denotes the used c. In the PRISM examples, this is followed by the path and the name of the program and property file. Then comes a list of the Rabin pairs. Every rabin pair consists of a good label and a bad label, followed by True if the presence of the label indicates the good/bad label, or False if the absence of the label indicates the good/bad label. Then we have the threshold for the number of steps. If 10num_restart^c is larger than this threshold, the program uses 10num_restart^c instead. Then we set the maximum number of restarts. Ideally, if Pgood>0, this number will never be reached. Finally, we have the type of example ['syntheticPgood','syntheticPmRm', 'normal', 'not_gridworld', 'not_gridworld_Big'].
In the following, we will define the test parameters that were used to conduct the experiments. The first entry denotes the used c. In the PRISM examples, this is followed by the path and the name of the program and property file. Then comes a list of the Rabin pairs. Every rabin pair consists of a good label and a bad label, followed by True if the presence of the label indicates the good/bad label, or False if the absence of the label indicates the good/bad label. If there is no good or no bad labels, just write down any string which does not appear as a label. Next comes the threshold for the number of steps. If 10num_restart^c is larger than this threshold, the program uses 10num_restart^c instead. Then we set the maximum number of restarts. Ideally, if Pgood>0, this number will never be reached. Finally, we have the type of example ['syntheticPgood','syntheticPmRm', 'normal', 'not_gridworld', 'not_gridworld_Big'].
For the synthetic examples, the parameter sets looks slightly different because the Markov Chains will be generated directly by the Python code.
%% Cell type:code id: tags:
``` python
testparamsNandMultipleC = [[[i/5+1, 'PRISM/', 'nand_VG.pm', 'propNand.pctl', [
[['(z = 4)', True], ['x', True]]], 1000000, 100000], 'normal'] for i in range(11)]
testparamsCrowdsMultipleC = [[[i/5+1, 'PRISM/',
'crowds_nodl_VG.pm', 'propCrowds.pctl', [[['((observe0 = 0) & (observe1 = 0))', True], ['((observe0 = 0) & (observe1 = 0))', False]]],1000000,100000],'normal'] for i in range(11)]
```
%% Cell type:code id: tags:
``` python
testparamsNandStorm = [[[i+1, 'PRISM/', 'nand_VG.pm',
'propNand.pctl', [[['(z = 4)', True], ['x', True]]],1000000,100000], 'normal'] for i in range(3)]
testparamsScaleStorm = [[[i+1, 'PRISM/',
'scale_VG.pm', 'propScale.pctl', [[['goal', True], ['x', True]]],1000000,100000], 'normal'] for i in range(3)]
testparamsCrowdsStorm = [[[i+1, 'PRISM/',
'crowds_nodl_VG.pm', 'propCrowds.pctl', [[['((observe0 = 0) & (observe1 = 0))', True], ['((observe0 = 0) & (observe1 = 0))', False]]],1000000,100000], 'normal'] for i in range(3)]
testparamsHermanStorm = [[[i+1, 'PRISM/',
'herman17_oneinit_VG.pm', 'propHerman.pctl', [[['(x1 = 0)', True], ['x', True]]],1000000,100000], 'normal'] for i in range(3)]
testparamsGridWorldStorm = [[[1+i, 'PRISM/','gridworld_VG.pm', 'propGridworld.pctl', [[['(robot_row = 1)', True], ['(robot_row = 1)', False]],[['(janitor_row = 1)', False], ['(janitor_row = 1)', True]]],1000000,100000], 'normal'] for i in range(3)]
testparamsBluetoothStorm = [[[i+1, 'PRISM/',
'bluetooth_oneinit_VG.pm', 'propBluetooth.pctl', [[['(y1 < 50)', True], ['ee', True]]],1000000,100000], 'normal'] for i in range(3)]
```
%% Cell type:code id: tags:
``` python
testparamsPmRm = np.concatenate([[[[c+1, 0.5*0.5**(i/4), 2*np.sqrt(np.sqrt(2))**i, 0.5, 1000000, 100000], 'syntheticPmRm']
for i in range(27)] for c in range(3)], axis=0)
testparamsPgood = np.concatenate([[[[c+1, i, 0.75, 1000000, 100000], 'syntheticPgood']
for i in range(12)] for c in range(3)], axis=0)
```
%% Cell type:code id: tags:
``` python
testparamsNotGridWorldBig = [[[i+1, 'PRISM/',
'not_gridworld_VG_Big.pm', 'propNotGridworld.pctl', [[['janitorrow', True], ['x', True]], [['robotrow', False], ['x', True]]], 30000000, 100000], 'not_gridworld_Big'] for i in range(3)]
testparamsNotGridWorld = [[[i+1, 'PRISM/',
'not_gridworld_VG.pm', 'propNotGridworld.pctl', [[['(janitor_row = 1)', True], ['x', True]], [['(robot_row = 1)', False], ['x', True]]], 30000000,100000], 'not_gridworld'] for i in range(3)]
```
%% Cell type:markdown id: tags:
## Experiments
Now we just run the experiments. The method test_multiple_strategies_multiprocess takes two arguments: The test parameters, and the number of repetitions. You can alter the maximum number of cores used in the the function itself (Under PROCESSES). You can decide how to increment the nice-values of the individual processes in the do_one_testrun function.
*test_multiple_strategies_multiprocess* returns two lists: A list of the average number of steps before the final restart, and a list of the average number of restarts in total.
%% Cell type:code id: tags:
``` python
BluetoothResults = test_multiple_strategies_multiprocess(testparamsBluetoothStorm, 300)
NandResults = test_multiple_strategies_multiprocess(testparamsNandStorm, 300)
ScaleResults = test_multiple_strategies_multiprocess(testparamsScaleStorm, 300)
CrowdsResults = test_multiple_strategies_multiprocess(testparamsCrowdsStorm, 300)
HermanResults = test_multiple_strategies_multiprocess(testparamsHermanStorm, 100)
GridWorldResults = test_multiple_strategies_multiprocess(testparamsGridWorldStorm, 300)
NandMultipleCResults = test_multiple_strategies_multiprocess(testparamsNandMultipleC,1000)
CrowdsMultipleCResults = test_multiple_strategies_multiprocess(testparamsCrowdsMultipleC,1000)
```
%% Cell type:code id: tags:
``` python
NotGridWorldResults = test_multiple_strategies_multiprocess(testparamsNotGridWorld, 300)
NotGridWorldBigResults = test_multiple_strategies_multiprocess(testparamsNotGridWorldBig, 100)
```
%% Cell type:code id: tags:
``` python
SyntheticPmRmResults = test_multiple_strategies_multiprocess(testparamsPmRm,500)
SyntheticPGoodResults = test_multiple_strategies_multiprocess(testparamsPgood,500)
```
%% Cell type:markdown id: tags:
## Plots
Now we plot the data for the synthetic experiments and the Mulitple-C experiments.
%% Cell type:code id: tags:
``` python
import pylab
import matplotlib.pyplot as plt
import scipy.special
import seaborn as sns
sns.set_theme()
testparamsCNandMultipleC = [parameter[0][0] for parameter in testparamsNandMultipleC]
testparamsCCrowdsMultipleC = [parameter[0][0]
for parameter in testparamsCrowdsMultipleC]
c2Nand = np.array(
NandMultipleCResults[0])[0]
c2Crowds = np.array(
CrowdsMultipleCResults[0])[0]
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot()
line, = ax.plot(testparamsCNandMultipleC, np.array(
NandMultipleCResults[0])/c2Nand)
line.set_label('Nand dataset')
line2, = ax.plot(testparamsCCrowdsMultipleC, np.array(
CrowdsMultipleCResults[0])/c2Crowds)
line2.set_label('Crowds dataset')
```
%% Cell type:code id: tags:
``` python
nc=3
kmax = len(testparamsPgood)
testparamsPgoodX = np.array([parameters[0][1] for parameters in testparamsPgood])
ResultsPGood = SyntheticPGoodResults
k_ind=np.int64(kmax/(nc))
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot()
ofs=0
Pgood=0.75
ofsb =7
for i in range(nc):
c=i+1
line, = ax.plot(np.log(0.75**(-testparamsPgoodX[i*k_ind:(i+1)*k_ind-ofs])),
np.log(ResultsPGood[0][i*k_ind:(i+1)*k_ind-ofs]), color=plt.cm.tab10(i))
line.set_label('c='+str(c))
coef = np.polyfit(np.log(Pgood**-testparamsPgoodX[np.int64((i)*k_ind)+ofsb:(i+1)*k_ind]),
np.log(ResultsPGood[0][np.int64((i)*k_ind)+ofsb:(i+1)*k_ind]), 1)
poly1d_fn = np.poly1d(coef)
print(coef)
line2, = ax.plot(np.log(Pgood**-testparamsPgoodX[i*k_ind+ofsb:(i+1)*k_ind-ofs]),
poly1d_fn(np.log(Pgood**-testparamsPgoodX[i*k_ind+ofsb:(i+1)*k_ind-ofs])), '--k')
line2.set_label('y=x^'+str(np.round(coef[0],2)))
ax.legend()
ax.set_xlabel("log(1/P_good)")
ax.set_ylabel("log(E[#steps])")
#plt.savefig("SyntheticExperimentPGood.pdf")
plt.show()
```
%% Cell type:code id: tags:
``` python
import pylab
import matplotlib.pyplot as plt
import scipy.special
nc=3
resultsX = SyntheticPmRmResults
testparamsX = np.array([parameters[0] for parameters in testparamsPmRm])
kmax=len(testparamsX)
k_ind=np.int64(kmax/(nc))
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot()
ofs=5
ofsb =10
for i in range(nc):
c=i+1
line, = ax.plot(np.log(testparamsX[i*k_ind:(i+1)*k_ind-ofs, 2])-np.log(testparamsX[i*k_ind:(i+1)*k_ind-ofs, 1]),
np.log(resultsX[0][i*k_ind:(i+1)*k_ind-ofs]), color=plt.cm.tab10(i))
line.set_label('c='+str(i+1))
coef = np.polyfit(np.log(testparamsX[np.int64((i+0.5)*k_ind):(i+1)*k_ind, 2])-np.log(testparamsX[np.int64((i+0.5)*k_ind):(i+1)*k_ind, 1]),
np.log(resultsX[0][np.int64((i+0.5)*k_ind):(i+1)*k_ind]), 1)
poly1d_fn = np.poly1d(coef)
print(coef)
line2, = ax.plot(np.log(testparamsX[i*k_ind+ofsb:(i+1)*k_ind-ofs, 2])-np.log(testparamsX[i*k_ind+ofsb:(i+1)*k_ind-ofs, 1]),
poly1d_fn(np.log(testparamsX[i*k_ind+ofsb:(i+1)*k_ind-ofs, 2])-np.log(testparamsX[i*k_ind+ofsb:(i+1)*k_ind-ofs, 1])), '--k')
line2.set_label('y=x^'+str(np.round(coef[0],2)))
ax.legend()
ax.set_xlabel("log(R_m/P_m)")
ax.set_ylabel("log(E[#steps])")
#plt.savefig("SyntheticExperimentPmRm.pdf")
plt.show()
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment