diff --git a/example_two_car_lane_switch.py b/example_two_car_lane_switch.py
index bd8f4798a1e57b00365033daff6d1abcf15958cd..69b9b9250725ea088b49dd9df42f3763f6e2f498 100644
--- a/example_two_car_lane_switch.py
+++ b/example_two_car_lane_switch.py
@@ -1,30 +1,39 @@
-from enum import Enum,auto
+import matplotlib.pyplot as plt
+from ourtool.agents.car_agent import CarAgent
+import numpy as np
+from user.simple_map import SimpleMap, SimpleMap2
+from ourtool.scenario.scenario import Scenario
+from enum import Enum, auto
from ourtool.map.lane_map import LaneMap
+
class VehicleMode(Enum):
Normal = auto()
SwitchLeft = auto()
SwitchRight = auto()
Brake = auto()
+
class LaneMode(Enum):
Lane0 = auto()
Lane1 = auto()
Lane2 = auto()
+
class State:
x = 0.0
y = 0.0
theta = 0.0
v = 0.0
- vehicle_mode:VehicleMode = VehicleMode.Normal
- lane_mode:LaneMode = LaneMode.Lane0
+ vehicle_mode: VehicleMode = VehicleMode.Normal
+ lane_mode: LaneMode = LaneMode.Lane0
- def __init__(self,x,y,theta,v,vehicle_mode:VehicleMode, lane_mode:LaneMode):
+ def __init__(self, x, y, theta, v, vehicle_mode: VehicleMode, lane_mode: LaneMode):
self.data = []
-def controller(ego:State, other:State, lane_map):
+
+def controller(ego: State, other: State, lane_map):
output = ego
if ego.vehicle_mode == VehicleMode.Normal:
if other.x - ego.x > 3 and other.x - ego.x < 5 and ego.lane_mode == other.lane_mode:
@@ -34,58 +43,53 @@ def controller(ego:State, other:State, lane_map):
if lane_map.has_right(ego.lane_mode):
output.vehicle_mode = VehicleMode.SwitchRight
if ego.vehicle_mode == VehicleMode.SwitchLeft:
- if ego.y >= lane_map.lane_geometry(ego.lane_mode)-0.5:
+ if ego.y >= lane_map.lane_geometry(ego.lane_mode)-2.5:
output.vehicle_mode = VehicleMode.Normal
output.lane_mode = lane_map.left_lane(ego.lane_mode)
if ego.vehicle_mode == VehicleMode.SwitchRight:
- if ego.y <= lane_map.lane_geometry(ego.lane_mode)+0.5:
+ if ego.y <= lane_map.lane_geometry(ego.lane_mode)+2.5:
output.vehicle_mode = VehicleMode.Normal
output.lane_mode = lane_map.right_lane(ego.lane_mode)
-
+
return output
-
-from ourtool.agents.car_agent import CarAgent
-from ourtool.scenario.scenario import Scenario
-from user.simple_map import SimpleMap, SimpleMap2
-import matplotlib.pyplot as plt
-import numpy as np
+
if __name__ == "__main__":
input_code_name = 'example_two_car_lane_switch.py'
scenario = Scenario()
-
+
car = CarAgent('car1', file_name=input_code_name)
scenario.add_agent(car)
car = CarAgent('car2', file_name=input_code_name)
scenario.add_agent(car)
scenario.add_map(SimpleMap2())
scenario.set_init(
- [[0,-3,0,1.0], [10,-3,0,0.5]],
+ [[10, 3, 0, 0.5], [0, 3, 0, 1.0]],
[
- (VehicleMode.Normal, LaneMode.Lane2),
- (VehicleMode.Normal, LaneMode.Lane2)
+ (VehicleMode.Normal, LaneMode.Lane0),
+ (VehicleMode.Normal, LaneMode.Lane0)
]
)
# simulator = Simulator()
- traces = scenario.simulate(40)
+ traces = scenario.verify(40)
- plt.plot([0,40],[3,3],'g')
- plt.plot([0,40],[0,0],'g')
- plt.plot([0,40],[-3,-3],'g')
+ plt.plot([0, 40], [3, 3], 'g')
+ plt.plot([0, 40], [0, 0], 'g')
+ plt.plot([0, 40], [-3, -3], 'g')
queue = [traces]
- while queue!=[]:
+ while queue != []:
node = queue.pop(0)
traces = node.trace
# for agent_id in traces:
agent_id = 'car2'
trace = np.array(traces[agent_id])
- plt.plot(trace[:,1], trace[:,2], 'r')
+ plt.plot(trace[:, 1], trace[:, 2], 'r')
agent_id = 'car1'
trace = np.array(traces[agent_id])
- plt.plot(trace[:,1], trace[:,2], 'b')
+ plt.plot(trace[:, 1], trace[:, 2], 'b')
# if node.child != []:
- queue += node.child
+ queue += node.child
plt.show()
diff --git a/ourtool/simulator/__init__.py b/ourtool/analysis/__init__.py
similarity index 100%
rename from ourtool/simulator/__init__.py
rename to ourtool/analysis/__init__.py
diff --git a/ourtool/analysis/analysis_tree_node.py b/ourtool/analysis/analysis_tree_node.py
new file mode 100644
index 0000000000000000000000000000000000000000..55606657c6979470576165eb2d5701203249042b
--- /dev/null
+++ b/ourtool/analysis/analysis_tree_node.py
@@ -0,0 +1,25 @@
+from typing import List, Dict
+
+class AnalysisTreeNode:
+ """AnalysisTreeNode class
+ A AnalysisTreeNode stores the continous execution of the system without transition happening"""
+ trace: Dict
+ """The trace for each agent.
+ The key of the dict is the agent id and the value of the dict is simulated traces for each agent"""
+ init: Dict
+
+ def __init__(
+ self,
+ trace={},
+ init={},
+ mode={},
+ agent={},
+ child=[],
+ start_time = 0
+ ):
+ self.trace:Dict = trace
+ self.init:Dict = init
+ self.mode:Dict = mode
+ self.agent:Dict = agent
+ self.child:List[AnalysisTreeNode] = child
+ self.start_time:float = start_time
diff --git a/ourtool/simulator/simulator.py b/ourtool/analysis/simulator.py
similarity index 82%
rename from ourtool/simulator/simulator.py
rename to ourtool/analysis/simulator.py
index ea80c9afa95483c2422d276ae279042b058f0312..10a1b23f77a4d51b4b84c0581310773df007bdda 100644
--- a/ourtool/simulator/simulator.py
+++ b/ourtool/analysis/simulator.py
@@ -1,31 +1,10 @@
from typing import List, Dict
-from ourtool.agents.base_agent import BaseAgent
-import numpy as np
import copy
-class SimulationTreeNode:
- """SimulationTreeNode class
- A SimulationTreeNode stores the continous execution of the system without transition happening"""
- trace: Dict
- """The trace for each agent.
- The key of the dict is the agent id and the value of the dict is simulated traces for each agent"""
-
- init: Dict
- def __init__(
- self,
- trace={},
- init={},
- mode={},
- agent={},
- child=[],
- start_time = 0
- ):
- self.trace:Dict = trace
- self.init:Dict = init
- self.mode:Dict = mode
- self.agent:Dict = agent
- self.child:List[SimulationTreeNode] = child
- self.start_time:float = start_time
+import numpy as np
+
+from ourtool.agents.base_agent import BaseAgent
+from ourtool.analysis.analysis_tree_node import AnalysisTreeNode
class Simulator:
def __init__(self):
@@ -33,7 +12,7 @@ class Simulator:
def simulate(self, init_list, init_mode_list, agent_list:List[BaseAgent], transition_graph, time_horizon, lane_map):
# Setup the root of the simulation tree
- root = SimulationTreeNode()
+ root = AnalysisTreeNode()
for i, agent in enumerate(agent_list):
root.init[agent.id] = init_list[i]
init_mode = init_mode_list[i][0].name
@@ -46,7 +25,7 @@ class Simulator:
simulation_queue.append(root)
# Perform BFS through the simulation tree to loop through all possible transitions
while simulation_queue != []:
- node:SimulationTreeNode = simulation_queue.pop(0)
+ node:AnalysisTreeNode = simulation_queue.pop(0)
print(node.mode)
remain_time = time_horizon - node.start_time
if remain_time <= 0:
@@ -88,7 +67,7 @@ class Simulator:
# copy the traces that are not under transition
for transition in transitions:
transit_agent_idx, src_mode, dest_mode, next_init, idx = transition
- # next_node = SimulationTreeNode(trace = {},init={},mode={},agent={}, child = [], start_time = 0)
+ # next_node = AnalysisTreeNode(trace = {},init={},mode={},agent={}, child = [], start_time = 0)
next_node_mode = copy.deepcopy(node.mode)
next_node_mode[transit_agent_idx] = dest_mode
next_node_agent = node.agent
@@ -101,7 +80,7 @@ class Simulator:
else:
next_node_trace[agent_idx] = truncated_trace[agent_idx]
- tmp = SimulationTreeNode(
+ tmp = AnalysisTreeNode(
trace = next_node_trace,
init = next_node_init,
mode = next_node_mode,
@@ -112,7 +91,7 @@ class Simulator:
node.child.append(tmp)
simulation_queue.append(tmp)
# Put the node in the child of current node. Put the new node in the queue
- # node.child.append(SimulationTreeNode(
+ # node.child.append(AnalysisTreeNode(
# trace = next_node_trace,
# init = next_node_init,
# mode = next_node_mode,
diff --git a/ourtool/analysis/verifier.py b/ourtool/analysis/verifier.py
new file mode 100644
index 0000000000000000000000000000000000000000..528a83e4eea4c4077b1f72d8cf117a6c45c49242
--- /dev/null
+++ b/ourtool/analysis/verifier.py
@@ -0,0 +1,85 @@
+from typing import List, Dict
+import copy
+
+import numpy as np
+
+from ourtool.agents.base_agent import BaseAgent
+from ourtool.analysis.analysis_tree_node import AnalysisTreeNode
+from ourtool.dryvr.core.dryvrcore import calc_bloated_tube
+import ourtool.dryvr.common.config as userConfig
+
+class Verifier:
+ def __init__(self):
+ self.reachtube_tree_root = None
+ self.unsafe_set = None
+ self.verification_result = None
+
+ def compute_full_reachtube(
+ self,
+ init_list,
+ init_mode_list,
+ agent_list:List[BaseAgent],
+ transition_graph,
+ time_horizon,
+ lane_map
+ ):
+ root = AnalysisTreeNode()
+ for i, agent in enumerate(agent_list):
+ root.init[agent.id] = init_list[i]
+ init_mode = [elem.name for elem in init_mode_list[i]]
+ init_mode =','.join(init_mode)
+ root.mode[agent.id] = init_mode
+ root.agent[agent.id] = agent
+ self.reachtube_tree_root = root
+ verification_queue = []
+ verification_queue.append(root)
+ while verification_queue != []:
+ node:AnalysisTreeNode = verification_queue.pop(0)
+ print(node.mode)
+ remain_time = time_horizon - node.start_time
+ if remain_time <= 0:
+ continue
+ # For reachtubes not already computed
+ for agent_id in node.agent:
+ if agent_id not in node.trace:
+ # Compute the trace starting from initial condition
+ mode = node.mode[agent_id]
+ init = node.init[agent_id]
+ # trace = node.agent[agent_id].TC_simulate(mode, init, remain_time,lane_map)
+ # trace[:,0] += node.start_time
+ # node.trace[agent_id] = trace.tolist()
+
+ cur_bloated_tube = calc_bloated_tube(mode,
+ init,
+ remain_time,
+ node.agent[agent_id].TC_simulate,
+ 'GLOBAL',
+ None,
+ userConfig.SIMTRACENUM,
+ lane_map = lane_map
+ )
+ trace = np.array(cur_bloated_tube)
+ trace[:,0] += node.start_time
+ node.trace[agent_id] = trace.tolist()
+ print("here")
+
+ trace_length = len(list(node.trace.values())[0])/2
+ transitions = []
+ for idx in range(0,trace_length,2):
+ # For each trace, check with the guard to see if there's any possible transition
+ # Store all possible transition in a list
+ # A transition is defined by (agent, src_mode, dest_mode, corresponding reset, transit idx)
+ # Here we enforce that only one agent transit at a time
+ all_agent_state = {}
+ for agent_id in node.agent:
+ all_agent_state[agent_id] = (node.trace[agent_id][idx:idx+2], node.mode[agent_id])
+ possible_transitions = transition_graph.get_all_transitions_set(all_agent_state)
+ if possible_transitions != []:
+ any_contained = False
+ for agent_idx, src_mode, dest_mode, next_init, contained in possible_transitions:
+ transitions.append((agent_idx, src_mode, dest_mode, next_init, idx))
+ any_contained = any_contained or contained
+ if any_contained:
+ break
+ pass
+
diff --git a/ourtool/dryvr/__init__.py b/ourtool/dryvr/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ourtool/dryvr/common/__init__.py b/ourtool/dryvr/common/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..7dbea944864d1e4d5337014013f9b49ae5ce1bb1
--- /dev/null
+++ b/ourtool/dryvr/common/__init__.py
@@ -0,0 +1,22 @@
+# _oo0oo_
+# o8888888o
+# 88" . "88
+# (| -_- |)
+# 0\ = /0
+# ___/`---'\___
+# .' \\| |// '.
+# / \\||| : |||// \
+# / _||||| -:- |||||- \
+# | | \\\ - /// | |
+# | \_| ''\---/'' |_/ |
+# \ .-\__ '-' ___/-. /
+# ___'. .' /--.--\ `. .'___
+# ."" '< `.___\_<|>_/___.' >' "".
+# | | : `- \`.;`\ _ /`;.`/ - ` : | |
+# \ \ `_. \_ __\ /__ _/ .-` / /
+# =====`-.____`.___ \_____/___.-`___.-'=====
+# `=---='
+#
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Codes are far away from bugs with the protection
diff --git a/ourtool/dryvr/common/config.py b/ourtool/dryvr/common/config.py
new file mode 100644
index 0000000000000000000000000000000000000000..6b23f518d74bd9d178fe2f56ad613ad65c0ea1a5
--- /dev/null
+++ b/ourtool/dryvr/common/config.py
@@ -0,0 +1,11 @@
+""" Global Configurations """
+
+# Verification config
+SIMUTESTNUM = 10
+SIMTRACENUM = 10
+REFINETHRES = 10
+CHILDREFINETHRES = 2
+
+# Synthesis config
+RANDMODENUM = 3
+RANDSECTIONNUM = 3
diff --git a/ourtool/dryvr/common/constant.py b/ourtool/dryvr/common/constant.py
new file mode 100644
index 0000000000000000000000000000000000000000..0eaaa69353ee78fe8b7023ec6d73a12550213335
--- /dev/null
+++ b/ourtool/dryvr/common/constant.py
@@ -0,0 +1,28 @@
+"""
+This file contains constant for DryVR
+"""
+
+# Logic constant
+DEBUG = False
+PLOTGRAPH = True
+
+# Flag value
+SAFE = 1
+UNSAFE = -1
+UNKNOWN = 0
+NOSTATUS = 99
+BLOATDEBUG = False
+PLOTDIM = 2
+
+# File Paths
+GRAPHOUTPUT = 'output/curGraph.png'
+RRTGRAPHPOUTPUT = 'output/rrtGraph.png'
+SIMRESULTOUTPUT = 'output/Traj.txt'
+REACHTUBEOUTPUT = 'output/reachtube.txt'
+RRTOUTPUT = 'output/rrtTube.txt'
+UNSAFEFILENAME = 'output/unsafeTube.txt'
+
+# Useful constant string
+NEWLINE = '----------------------------------------------'
+PW = "PW"
+GLOBAL = "GLOBAL"
diff --git a/ourtool/dryvr/common/io.py b/ourtool/dryvr/common/io.py
new file mode 100644
index 0000000000000000000000000000000000000000..c5b55e206966a7b8a5adb4fc016f5da359ec2b3d
--- /dev/null
+++ b/ourtool/dryvr/common/io.py
@@ -0,0 +1,156 @@
+"""
+This file contains IO functions for DryVR
+"""
+
+import six
+
+from ourtool.dryvr.common.utils import DryVRInput, RrtInput, checkVerificationInput, checkSynthesisInput
+
+
+def writeReachTubeFile(result, path):
+ """
+ Write reach tube to a file
+
+ reach tube example:
+ mode1
+ [0.0, 1, 2]
+ [0.1, 2, 3]
+ [0.1, 2, 3]
+ ....
+ mode1->mode2
+ [1.0, 3, 4]
+ ....
+
+ Args:
+ result (list): list of reachable state.
+ path (str): file name.
+
+ Returns:
+ None
+
+ """
+ with open(path, 'w') as f:
+ for line in result:
+ if isinstance(line, six.string_types):
+ f.write(line + '\n')
+ elif isinstance(line, list):
+ f.write(' '.join(map(str, line)) + '\n')
+
+
+def writeRrtResultFile(modes, traces, path):
+ """
+ Write control synthesis result to a file
+
+ Args:
+ modes (list): list of mode.
+ traces (list): list of traces corresponding to modes
+ path (str): file name.
+
+ Returns:
+ None
+
+ """
+ with open(path, 'w') as f:
+ for mode, trace in zip(modes, traces):
+ f.write(mode + '\n')
+ for line in trace:
+ f.write(" ".join(map(str, line)) + '\n')
+
+
+def parseVerificationInputFile(data):
+ """
+ Parse the json input for DryVR verification
+
+ Args:
+ data (dict): dictionary contains parameters
+
+ Returns:
+ DryVR verification input object
+
+ """
+
+ # If resets is missing, fill with empty resets
+ if not 'resets' in data:
+ data['resets'] = ["" for _ in range(len(data["edge"]))]
+
+ # If initialMode is missing, fill with empty initial mode
+ if not 'initialVertex' in data:
+ data['initialVertex'] = -1
+
+ # If deterministic is missing, default to non-deterministic
+ if not 'deterministic' in data:
+ data['deterministic'] = False
+
+ # If bloating method is missing, default global descrepancy
+ if not 'bloatingMethod' in data:
+ data['bloatingMethod'] = 'GLOBAL'
+
+ # Set a fake kvalue since kvalue is not used in this case
+
+ if data['bloatingMethod'] == "GLOBAL":
+ data['kvalue'] = [1.0 for i in range(len(data['variables']))]
+
+ # Set a fake directory if the directory is not provided, this means the user provides
+ # simulation function to DryVR directly
+ if not 'directory' in data:
+ data['directory'] = ""
+
+ checkVerificationInput(data)
+ return DryVRInput(
+ vertex=data["vertex"],
+ edge=data["edge"],
+ guards=data["guards"],
+ variables=data["variables"],
+ initialSet=data["initialSet"],
+ unsafeSet=data["unsafeSet"],
+ timeHorizon=data["timeHorizon"],
+ path=data["directory"],
+ resets=data["resets"],
+ initialVertex=data["initialVertex"],
+ deterministic=data["deterministic"],
+ bloatingMethod=data['bloatingMethod'],
+ kvalue=data['kvalue'],
+ )
+
+
+def parseRrtInputFile(data):
+ """
+ Parse the json input for DryVR controller synthesis
+
+ Args:
+ data (dict): dictionary contains parameters
+
+ Returns:
+ DryVR controller synthesis input object
+
+ """
+
+ # If bloating method is missing, default global descrepancy
+ if not 'bloatingMethod' in data:
+ data['bloatingMethod'] = 'GLOBAL'
+
+ # set a fake kvalue since kvalue is not used in this case
+
+ if data['bloatingMethod'] == "GLOBAL":
+ data['kvalue'] = [1.0 for i in range(len(data['variables']))]
+
+ # Set a fake directory if the directory is not provided, this means the user provides
+ # simulation function to DryVR directly
+ if not 'directory' in data:
+ data['directory'] = ""
+
+ checkSynthesisInput(data)
+
+ return RrtInput(
+ modes=data["modes"],
+ variables=data["variables"],
+ initialSet=data["initialSet"],
+ unsafeSet=data["unsafeSet"],
+ goalSet=data["goalSet"],
+ timeHorizon=data["timeHorizon"],
+ minTimeThres=data["minTimeThres"],
+ path=data["directory"],
+ goal=data["goal"],
+ bloatingMethod=data['bloatingMethod'],
+ kvalue=data['kvalue'],
+ )
diff --git a/ourtool/dryvr/common/utils.py b/ourtool/dryvr/common/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..21d64ff954b9f63fff37f7595b33a37808234b39
--- /dev/null
+++ b/ourtool/dryvr/common/utils.py
@@ -0,0 +1,304 @@
+"""
+This file contains common utils for DryVR
+"""
+
+import importlib
+import random
+
+from collections import namedtuple
+
+# This is the tuple for input file parsed by DryVR
+DryVRInput = namedtuple(
+ 'DryVRInput',
+ 'vertex edge guards variables initialSet unsafeSet timeHorizon path resets initialVertex deterministic '
+ 'bloatingMethod kvalue '
+)
+
+# This is the tuple for rtt input file parsed by DryVR
+RrtInput = namedtuple(
+ 'RttInput',
+ 'modes variables initialSet unsafeSet goalSet timeHorizon minTimeThres path goal bloatingMethod kvalue'
+)
+
+
+def importSimFunction(path):
+ """
+ Load simulation function from given file path
+ Note the folder in the examples directory must have __init__.py
+ And the simulation function must be named TC_Simulate
+ The function should looks like following:
+ TC_Simulate(Mode, initialCondition, time_bound)
+
+ Args:
+ path (str): Simulator directory.
+
+ Returns:
+ simulation function
+
+ """
+ try:
+ # FIXME TC_Simulate should just be a simple call back function
+ import os
+ import sys
+ sys.path.append(os.path.abspath(path))
+ mod_name = path.replace('/', '.')
+ module = importlib.import_module(mod_name)
+ sys.path.pop()
+
+ return module.TC_Simulate
+ except ImportError as e:
+ print("Import simulation function failed!")
+ print(e)
+ exit() # TODO Proper return
+
+
+def randomPoint(lower, upper):
+ """
+ Pick a random point between lower and upper bound
+ This function supports both int or list
+
+ Args:
+ lower (list or int or float): lower bound.
+ upper (list or int or float): upper bound.
+
+ Returns:
+ random point (either float or list of float)
+
+ """
+ # random.seed(4)
+
+ if isinstance(lower, int) or isinstance(lower, float):
+ return random.uniform(lower, upper)
+
+ if isinstance(lower, list):
+ assert len(lower) == len(upper), "Random Point List Range Error"
+
+ return [random.uniform(lower[i], upper[i]) for i in range(len(lower))]
+
+
+def calcDelta(lower, upper):
+ """
+ Calculate the delta value between the lower and upper bound
+ The function only supports list since we assue initial set is always list
+
+ Args:
+ lower (list): lowerbound.
+ upper (list): upperbound.
+
+ Returns:
+ delta (list of float)
+
+ """
+ # Convert list into float in case they are int
+ lower = [float(val) for val in lower]
+ upper = [float(val) for val in upper]
+
+ assert len(lower) == len(upper), "Delta calc List Range Error"
+ return [(upper[i] - lower[i]) / 2 for i in range(len(upper))]
+
+
+def calcCenterPoint(lower, upper):
+ """
+ Calculate the center point between the lower and upper bound
+ The function only supports list since we assue initial set is always list
+
+ Args:
+ lower (list): lowerbound.
+ upper (list): upperbound.
+
+ Returns:
+ delta (list of float)
+
+ """
+
+ # Convert list into float in case they are int
+ lower = [float(val) for val in lower]
+ upper = [float(val) for val in upper]
+ assert len(lower) == len(upper), "Center Point List Range Error"
+ return [(upper[i] + lower[i]) / 2 for i in range(len(upper))]
+
+
+def buildModeStr(g, vertex):
+ """
+ Build a unique string to represent a mode
+ This should be something like "modeName,modeNum"
+
+ Args:
+ g (igraph.Graph): Graph object.
+ vertex (int or str): vertex number.
+
+ Returns:
+ str: a string to represent a mode
+
+ """
+ return g.vs[vertex]['label'] + ',' + str(vertex)
+
+
+def handleReplace(input_str, keys):
+ """
+ Replace variable in inputStr to self.varDic["variable"]
+ For example:
+ input
+ And(y<=0,t>=0.2,v>=-0.1)
+ output:
+ And(self.varDic["y"]<=0,self.varDic["t"]>=0.2,self.varDic["v"]>=-0.1)
+
+ Args:
+ input_str (str): original string need to be replaced
+ keys (list): list of variable strings
+
+ Returns:
+ str: a string that all variables have been replaced into a desire form
+
+ """
+ idxes = []
+ i = 0
+ original = input_str
+
+ keys.sort(key=lambda s: len(s))
+ for key in keys[::-1]:
+ for i in range(len(input_str)):
+ if input_str[i:].startswith(key):
+ idxes.append((i, i + len(key)))
+ input_str = input_str[:i] + "@" * len(key) + input_str[i + len(key):]
+
+ idxes = sorted(idxes)
+
+ input_str = original
+ for idx in idxes[::-1]:
+ key = input_str[idx[0]:idx[1]]
+ target = 'self.varDic["' + key + '"]'
+ input_str = input_str[:idx[0]] + target + input_str[idx[1]:]
+ return input_str
+
+
+def neg(orig):
+ """
+ Neg the original condition
+ For example:
+ input
+ And(y<=0,t>=0.2,v>=-0.1)
+ output:
+ Not(And(y<=0,t>=0.2,v>=-0.1))
+
+ Args:
+ orig (str): original string need to be neg
+
+ Returns:
+ a neg condition string
+
+ """
+ return 'Not(' + orig + ')'
+
+
+def trimTraces(traces):
+ """
+ trim all traces to the same length
+
+ Args:
+ traces (list): list of traces generated by simulator
+ Returns:
+ traces (list) after trim to the same length
+
+ """
+
+ ret_traces = []
+ trace_lengths = []
+ for trace in traces:
+ trace_lengths.append(len(trace))
+ trace_len = min(trace_lengths)
+
+ for trace in traces:
+ ret_traces.append(trace[:trace_len])
+
+ return ret_traces
+
+
+def checkVerificationInput(data):
+ """
+ Check verification input to make sure it is valid
+
+ Args:
+ data (obj): json data object
+ Returns:
+ None
+
+ """
+ assert len(data['variables']) == len(data['initialSet'][0]), "Initial set dimension mismatch"
+
+ assert len(data['variables']) == len(data['initialSet'][1]), "Initial set dimension mismatch"
+
+ assert len(data['edge']) == len(data["guards"]), "guard number mismatch"
+
+ assert len(data['edge']) == len(data["resets"]), "reset number mismatch"
+
+ if data["bloatingMethod"] == "PW":
+ assert 'kvalue' in data, "kvalue need to be provided when bloating method set to PW"
+
+ for i in range(len(data['variables'])):
+ assert data['initialSet'][0][i] <= data['initialSet'][1][i], "initial set lowerbound is larger than upperbound"
+
+
+def checkSynthesisInput(data):
+ """
+ Check Synthesis input to make sure it is valid
+
+ Args:
+ data (obj): json data object
+ Returns:
+ None
+
+ """
+ assert len(data['variables']) == len(data['initialSet'][0]), "Initial set dimension mismatch"
+
+ assert len(data['variables']) == len(data['initialSet'][1]), "Initial set dimension mismatch"
+
+ for i in range(len(data['variables'])):
+ assert data['initialSet'][0][i] <= data['initialSet'][1][i], "initial set lowerbound is larger than upperbound"
+
+ assert data["minTimeThres"] < data["timeHorizon"], "min time threshold is too large!"
+
+ if data["bloatingMethod"] == "PW":
+ assert 'kvalue' in data, "kvalue need to be provided when bloating method set to PW"
+
+
+def isIpynb():
+ """
+ Check if the code is running on Ipython notebook
+ """
+ try:
+ cfg = get_ipython().config
+ if "IPKernelApp" in cfg:
+ return True
+ else:
+ return False
+ except NameError:
+ return False
+
+
+def overloadConfig(configObj, userConfig):
+ """
+ Overload user config to config module
+
+ Args:
+ configObj (module): config module
+ userConfig (dict): user specified config
+ """
+
+ if "SIMUTESTNUM" in userConfig:
+ configObj.SIMUTESTNUM = userConfig["SIMUTESTNUM"]
+
+ if "SIMTRACENUM" in userConfig:
+ configObj.SIMTRACENUM = userConfig["SIMTRACENUM"]
+
+ if "REFINETHRES" in userConfig:
+ configObj.REFINETHRES = userConfig["REFINETHRES"]
+
+ if "CHILDREFINETHRES" in userConfig:
+ configObj.CHILDREFINETHRES = userConfig["CHILDREFINETHRES"]
+
+ if "RANDMODENUM" in userConfig:
+ configObj.RANDMODENUM = userConfig["RANDMODENUM"]
+
+ if "RANDSECTIONNUM" in userConfig:
+ configObj.RANDSECTIONNUM = userConfig["RANDSECTIONNUM"]
diff --git a/ourtool/dryvr/core/__init__.py b/ourtool/dryvr/core/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ourtool/dryvr/core/distance.py b/ourtool/dryvr/core/distance.py
new file mode 100644
index 0000000000000000000000000000000000000000..8da15a12091e3501880ca1411bb5fc4e4156099d
--- /dev/null
+++ b/ourtool/dryvr/core/distance.py
@@ -0,0 +1,52 @@
+"""
+This file contains a distance checker class for controller synthesis
+"""
+from math import sqrt
+
+
+class DistChecker:
+ """
+ This is the class to calculate the distance between
+ current set and the goal set for DryVR controller synthesis.
+ The distance it calculate is the euclidean distance.
+ """
+
+ def __init__(self, goal_set, variables):
+ """
+ Distance checker class initialization function.
+
+ Args:
+ goal_set (list): a list describing the goal set.
+ [["x_1","x_2"],[19.5,-1],[20.5,1]]
+ which defines a goal set
+ 19.5<=x_1<=20.5 && -1<=x_2<=1
+ variables (list): list of variable name
+ """
+ var, self.lower, self.upper = goal_set
+ self.idx = []
+ for v in var:
+ self.idx.append(variables.index(v) + 1)
+
+ def calc_distance(self, lower_bound, upper_bound):
+ """
+ Calculate the euclidean distance between the
+ current set and goal set.
+
+ Args:
+ lower_bound (list): the lower bound of the current set.
+ upper_bound (list): the upper bound of the current set.
+
+ Returns:
+ the euclidean distance between current set and goal set
+
+ """
+ dist = 0.0
+ for i in range(len(self.idx)):
+ max_val = max(
+ (self.lower[i] - lower_bound[self.idx[i]]) ** 2,
+ (self.upper[i] - lower_bound[self.idx[i]]) ** 2,
+ (self.lower[i] - upper_bound[self.idx[i]]) ** 2,
+ (self.upper[i] - upper_bound[self.idx[i]]) ** 2
+ )
+ dist += max_val
+ return sqrt(dist)
diff --git a/ourtool/dryvr/core/dryvrcore.py b/ourtool/dryvr/core/dryvrcore.py
new file mode 100644
index 0000000000000000000000000000000000000000..29ff8cd3152d260cd3d85eadc64aaf6dcd57e595
--- /dev/null
+++ b/ourtool/dryvr/core/dryvrcore.py
@@ -0,0 +1,309 @@
+"""
+This file contains core functions used by DryVR
+"""
+from __future__ import print_function
+import random
+
+import networkx as nx
+import numpy as np
+import igraph
+
+
+from ourtool.dryvr.common.constant import *
+from ourtool.dryvr.common.io import writeReachTubeFile
+from ourtool.dryvr.common.utils import randomPoint, calcDelta, calcCenterPoint, trimTraces
+from ourtool.dryvr.discrepancy.Global_Disc import get_reachtube_segment
+# from ourtool.dryvr.tube_computer.backend.reachabilityengine import ReachabilityEngine
+# from ourtool.dryvr.tube_computer.backend.initialset import InitialSet
+
+def build_graph(vertex, edge, guards, resets):
+ """
+ Build graph object using given parameters
+
+ Args:
+ vertex (list): list of vertex with mode name
+ edge (list): list of edge that connects vertex
+ guards (list): list of guard corresponding to each edge
+ resets (list): list of reset corresponding to each edge
+
+ Returns:
+ graph object
+
+ """
+ g = igraph.Graph(directed=True)
+ g.add_vertices(len(vertex))
+ g.add_edges(edge)
+
+ g.vs['label'] = vertex
+ g.vs['name'] = vertex
+ labels = []
+ for i in range(len(guards)):
+ cur_guard = guards[i]
+ cur_reset = resets[i]
+ if not cur_reset:
+ labels.append(cur_guard)
+ else:
+ labels.append(cur_guard + '|' + cur_reset)
+
+ g.es['label'] = labels
+ g.es['guards'] = guards
+ g.es['resets'] = resets
+
+ # if PLOTGRAPH:
+ # graph = igraph.plot(g, GRAPHOUTPUT, margin=40)
+ # graph.save()
+ return g
+
+
+def build_rrt_graph(modes, traces, is_ipynb):
+ """
+ Build controller synthesis graph object using given modes and traces.
+ Note this function is very different from buildGraph function.
+ This is white-box transition graph learned from controller synthesis algorithm
+ The reason to build it is to output the transition graph to file
+
+ Args:
+ modes (list): list of mode name
+ traces (list): list of trace corresponding to each mode
+ is_ipynb (bool): check if it's in Ipython notebook environment
+
+ Returns:
+ None
+
+ """
+ if is_ipynb:
+ vertex = []
+ # Build unique identifier for a vertex and mode name
+ for idx, v in enumerate(modes):
+ vertex.append(v + "," + str(idx))
+
+ edge_list = []
+ edge_label = {}
+ for i in range(1, len(modes)):
+ edge_list.append((vertex[i - 1], vertex[i]))
+ lower = traces[i - 1][-2][0]
+ upper = traces[i - 1][-1][0]
+ edge_label[(vertex[i - 1], vertex[i])] = "[" + str(lower) + "," + str(upper) + "]"
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+ nx_graph = nx.DiGraph()
+ nx_graph.add_edges_from(edge_list)
+ pos = nx.spring_layout(nx_graph)
+ colors = ['green'] * len(nx_graph.nodes())
+ fig.suptitle('transition graph', fontsize=10)
+ nx.draw_networkx_labels(nx_graph, pos)
+ options = {
+ 'node_color': colors,
+ 'node_size': 1000,
+ 'cmap': plt.get_cmap('jet'),
+ 'arrowstyle': '-|>',
+ 'arrowsize': 50,
+ }
+ nx.draw_networkx(nx_graph, pos, arrows=True, **options)
+ nx.draw_networkx_edge_labels(nx_graph, pos, edge_labels=edge_label)
+ fig.canvas.draw()
+
+ else:
+ g = igraph.Graph(directed=True)
+ g.add_vertices(len(modes))
+ edges = []
+ for i in range(1, len(modes)):
+ edges.append([i - 1, i])
+ g.add_edges(edges)
+
+ g.vs['label'] = modes
+ g.vs['name'] = modes
+
+ # Build guard
+ guard = []
+ for i in range(len(traces) - 1):
+ lower = traces[i][-2][0]
+ upper = traces[i][-1][0]
+ guard.append("And(t>" + str(lower) + ", t<=" + str(upper) + ")")
+ g.es['label'] = guard
+ graph = igraph.plot(g, RRTGRAPHPOUTPUT, margin=40)
+ graph.save()
+
+
+def simulate(g, init_condition, time_horizon, guard, sim_func, reset, init_vertex, deterministic):
+ """
+ This function does a full hybrid simulation
+
+ Args:
+ g (obj): graph object
+ init_condition (list): initial point
+ time_horizon (float): time horizon to simulate
+ guard (src.core.guard.Guard): list of guard string corresponding to each transition
+ sim_func (function): simulation function
+ reset (src.core.reset.Reset): list of reset corresponding to each transition
+ init_vertex (int): initial vertex that simulation starts
+ deterministic (bool) : enable or disable must transition
+
+ Returns:
+ A dictionary obj contains simulation result.
+ Key is mode name and value is the simulation trace.
+
+ """
+
+ ret_val = igraph.defaultdict(list)
+ # If you do not declare initialMode, then we will just use topological sort to find starting point
+ if init_vertex == -1:
+ computer_order = g.topological_sorting(mode=igraph.OUT)
+ cur_vertex = computer_order[0]
+ else:
+ cur_vertex = init_vertex
+ remain_time = time_horizon
+ cur_time = 0
+
+ # Plus 1 because we need to consider about time
+ dimensions = len(init_condition) + 1
+
+ sim_result = []
+ # Avoid numeric error
+ while remain_time > 0.01:
+
+ if DEBUG:
+ print(NEWLINE)
+ print((cur_vertex, remain_time))
+ print('Current State', g.vs[cur_vertex]['label'], remain_time)
+
+ if init_condition is None:
+ # Ideally this should not happen
+ break
+
+ cur_successors = g.successors(cur_vertex)
+ transit_time = remain_time
+ cur_label = g.vs[cur_vertex]['label']
+
+ cur_sim_result = sim_func(cur_label, init_condition, transit_time)
+ if isinstance(cur_sim_result, np.ndarray):
+ cur_sim_result = cur_sim_result.tolist()
+
+ if len(cur_successors) == 0:
+ # Some model return numpy array, convert to list
+ init_condition, truncated_result = guard.guard_sim_trace(
+ cur_sim_result,
+ ""
+ )
+ cur_successor = None
+
+ else:
+ # First find all possible transition
+ # Second randomly pick a path and time to transit
+ next_modes = []
+ for cur_successor in cur_successors:
+ edge_id = g.get_eid(cur_vertex, cur_successor)
+ cur_guard_str = g.es[edge_id]['guards']
+ cur_reset_str = g.es[edge_id]['resets']
+
+ next_init, truncated_result = guard.guard_sim_trace(
+ cur_sim_result,
+ cur_guard_str
+ )
+
+ next_init = reset.reset_point(cur_reset_str, next_init)
+ # If there is a transition
+ if next_init:
+ next_modes.append((cur_successor, next_init, truncated_result))
+ if next_modes:
+ # It is a non-deterministic system, randomly choose next state to transit
+ if not deterministic:
+ cur_successor, init_condition, truncated_result = random.choice(next_modes)
+ # This is deterministic system, choose earliest transition
+ else:
+ shortest_time = float('inf')
+ for s, i, t in next_modes:
+ cur_tube_time = t[-1][0]
+ if cur_tube_time < shortest_time:
+ cur_successor = s
+ init_condition = i
+ truncated_result = t
+ shortest_time = cur_tube_time
+ else:
+ cur_successor = None
+ init_condition = None
+
+ # Get real transit time from truncated result
+ transit_time = truncated_result[-1][0]
+ ret_val[cur_label] += truncated_result
+ sim_result.append(cur_label)
+ for simRow in truncated_result:
+ simRow[0] += cur_time
+ sim_result.append(simRow)
+
+ remain_time -= transit_time
+ print("transit time", transit_time, "remain time", remain_time)
+ cur_time += transit_time
+ cur_vertex = cur_successor
+
+ writeReachTubeFile(sim_result, SIMRESULTOUTPUT)
+ return ret_val
+
+
+def calc_bloated_tube(
+ mode_label,
+ initial_set,
+ time_horizon,
+ sim_func,
+ bloating_method,
+ kvalue,
+ sim_trace_num,
+ guard_checker=None,
+ guard_str="",
+ lane_map = None
+ ):
+ """
+ This function calculate the reach tube for single given mode
+
+ Args:
+ mode_label (str): mode name
+ initial_set (list): a list contains upper and lower bound of the initial set
+ time_horizon (float): time horizon to simulate
+ sim_func (function): simulation function
+ bloating_method (str): determine the bloating method for reach tube, either GLOBAL or PW
+ sim_trace_num (int): number of simulations used to calculate the discrepancy
+ kvalue (list): list of float used when bloating method set to PW
+ guard_checker (src.core.guard.Guard or None): guard check object
+ guard_str (str): guard string
+
+ Returns:
+ Bloated reach tube
+
+ """
+ random.seed(4)
+ cur_center = calcCenterPoint(initial_set[0], initial_set[1])
+ cur_delta = calcDelta(initial_set[0], initial_set[1])
+ traces = [sim_func(mode_label, cur_center, time_horizon, lane_map)]
+ # Simulate SIMTRACENUM times to learn the sensitivity
+ for _ in range(sim_trace_num):
+ new_init_point = randomPoint(initial_set[0], initial_set[1])
+ traces.append(sim_func(mode_label, new_init_point, time_horizon, lane_map))
+
+ # Trim the trace to the same length
+ traces = trimTraces(traces)
+ if guard_checker is not None:
+ # pre truncated traces to get better bloat result
+ max_idx = -1
+ for trace in traces:
+ ret_idx = guard_checker.guard_sim_trace_time(trace, guard_str)
+ max_idx = max(max_idx, ret_idx + 1)
+ for i in range(len(traces)):
+ traces[i] = traces[i][:max_idx]
+
+ # The major
+ if bloating_method == GLOBAL:
+ # TODO: Replace this with ReachabilityEngine.get_reachtube_segment
+ cur_reach_tube: np.ndarray = get_reachtube_segment(np.array(traces), np.array(cur_delta), "PWGlobal")
+ # cur_reach_tube: np.ndarray = ReachabilityEngine.get_reachtube_segment_wrapper(np.array(traces), np.array(cur_delta))
+ elif bloating_method == PW:
+ # TODO: Replace this with ReachabilityEngine.get_reachtube_segment
+ cur_reach_tube: np.ndarray = get_reachtube_segment(np.array(traces), np.array(cur_delta), "PW")
+ # cur_reach_tube: np.ndarray = ReachabilityEngine.get_reachtube_segment_wrapper(np.array(traces), np.array(cur_delta))
+ else:
+ raise ValueError("Unsupported bloating method '" + bloating_method + "'")
+ final_tube = np.zeros((cur_reach_tube.shape[0]*2, cur_reach_tube.shape[2]))
+ final_tube[0::2, :] = cur_reach_tube[:, 0, :]
+ final_tube[1::2, :] = cur_reach_tube[:, 1, :]
+ print(final_tube.tolist()[-2], final_tube.tolist()[-1])
+ return final_tube.tolist()
diff --git a/ourtool/dryvr/core/dryvrmain.py b/ourtool/dryvr/core/dryvrmain.py
new file mode 100644
index 0000000000000000000000000000000000000000..af6ddd87302b6811d056b6ccffb4930da690bfea
--- /dev/null
+++ b/ourtool/dryvr/core/dryvrmain.py
@@ -0,0 +1,544 @@
+"""
+This file contains a single function that verifies model
+"""
+from __future__ import print_function
+import time
+
+import src.common.config as userConfig
+from src.common.io import parseVerificationInputFile, parseRrtInputFile, writeRrtResultFile
+from src.common.utils import buildModeStr, isIpynb, overloadConfig
+from src.core.distance import DistChecker
+from src.core.dryvrcore import *
+from src.core.goalchecker import GoalChecker
+from src.core.graph import Graph
+from src.core.guard import Guard
+from src.core.initialset import InitialSet
+from src.core.initialsetstack import InitialSetStack, GraphSearchNode
+from src.core.reachtube import ReachTube
+from src.core.reset import Reset
+from src.core.uniformchecker import UniformChecker
+# from src.tube_computer.backend.reachabilityengine import ReachabilityEngine
+# from src.tube_computer.backend.initialset import InitialSet
+
+def verify(data, sim_function, param_config=None):
+ """
+ DryVR verification algorithm.
+ It does the verification and print out the verify result.
+
+ Args:
+ data (dict): dictionary that contains params for the input file
+ sim_function (function): black-box simulation function
+ param_config (dict or None): user-specified configuration
+
+ Returns:
+ Safety (str): safety of the system
+ Reach (obj): reach tube object
+
+ """
+ if param_config is None:
+ param_config = {}
+ # There are some fields can be config by user,
+ # If user specified these fields in paramConfig,
+ # overload these parameters to userConfig
+ overloadConfig(userConfig, param_config)
+
+ refine_counter = 0
+
+ params = parseVerificationInputFile(data)
+ # Build the graph object
+ graph = build_graph(
+ params.vertex,
+ params.edge,
+ params.guards,
+ params.resets
+ )
+
+ # Build the progress graph for jupyter notebook
+ # isIpynb is used to detect if the code is running
+ # on notebook or terminal, the graph will only be shown
+ # in notebook mode
+ progress_graph = Graph(params, isIpynb())
+
+ # Make sure the initial mode is specified if the graph is dag
+ # FIXME should move this part to input check
+ # Bolun 02/12/2018
+ assert graph.is_dag() or params.initialVertex != -1, "Graph is not DAG and you do not have initial mode!"
+
+ checker = UniformChecker(params.unsafeSet, params.variables)
+ guard = Guard(params.variables)
+ reset = Reset(params.variables)
+ t_start = time.time()
+
+ # Step 1) Simulation Test
+ # Random generate points, then simulate and check the result
+ for _ in range(userConfig.SIMUTESTNUM):
+ rand_init = randomPoint(params.initialSet[0], params.initialSet[1])
+
+ if DEBUG:
+ print('Random checking round ', _, 'at point ', rand_init)
+
+ # Do a full hybrid simulation
+ sim_result = simulate(
+ graph,
+ rand_init,
+ params.timeHorizon,
+ guard,
+ sim_function,
+ reset,
+ params.initialVertex,
+ params.deterministic
+ )
+
+ # Check the traces for each mode
+ for mode in sim_result:
+ safety = checker.check_sim_trace(sim_result[mode], mode)
+ if safety == -1:
+ print('Current simulation is not safe. Program halt')
+ print('simulation time', time.time() - t_start)
+ return "UNSAFE", None
+ sim_end_time = time.time()
+
+ # Step 2) Check Reach Tube
+ # Calculate the over approximation of the reach tube and check the result
+ print("Verification Begin")
+
+ # Get the initial mode
+ if params.initialVertex == -1:
+ compute_order = graph.topological_sorting(mode=igraph.OUT)
+ initial_vertex = compute_order[0]
+ else:
+ initial_vertex = params.initialVertex
+
+ # Build the initial set stack
+ cur_mode_stack = InitialSetStack(initial_vertex, userConfig.REFINETHRES, params.timeHorizon,0)
+ cur_mode_stack.stack.append(InitialSet(params.initialSet[0], params.initialSet[1]))
+ cur_mode_stack.bloated_tube.append(buildModeStr(graph, initial_vertex))
+ while True:
+ # backward_flag can be SAFE, UNSAFE or UNKNOWN
+ # If the backward_flag is SAFE/UNSAFE, means that the children nodes
+ # of current nodes are all SAFE/UNSAFE. If one of the child node is
+ # UNKNOWN, then the backward_flag is UNKNOWN.
+ backward_flag = SAFE
+
+ while cur_mode_stack.stack:
+ print(str(cur_mode_stack))
+ print(cur_mode_stack.stack[-1])
+
+ if not cur_mode_stack.is_valid():
+ # A stack will be invalid if number of initial sets
+ # is more than refine threshold we set for each stack.
+ # Thus we declare this stack is UNKNOWN
+ print(cur_mode_stack.mode, "is not valid anymore")
+ backward_flag = UNKNOWN
+ break
+
+ # This is condition check to make sure the reach tube output file
+ # will be readable. Let me try to explain this.
+ # A reachtube output will be something like following
+ # MODEA->MODEB
+ # [0.0, 1.0, 1.1]
+ # [0.1, 1.1, 1.2]
+ # .....
+ # Once we have refinement, we will add multiple reach tube to
+ # this cur_mode_stack.bloatedTube
+ # However, we want to copy MODEA->MODEB so we know that two different
+ # reach tube from two different refined initial set
+ # The result will be look like following
+ # MODEA->MODEB
+ # [0.0, 1.0, 1.1]
+ # [0.1, 1.1, 1.2]
+ # .....
+ # MODEA->MODEB (this one gets copied!)
+ # [0.0, 1.5, 1.6]
+ # [0.1, 1.6, 1.7]
+ # .....
+ if isinstance(cur_mode_stack.bloated_tube[-1], list):
+ cur_mode_stack.bloated_tube.append(cur_mode_stack.bloated_tube[0])
+
+ cur_stack = cur_mode_stack.stack
+ cur_vertex = cur_mode_stack.mode
+ cur_remain_time = cur_mode_stack.remain_time
+ cur_label = graph.vs[cur_vertex]['label']
+ cur_successors = graph.successors(cur_vertex)
+ cur_initial = [cur_stack[-1].lower_bound, cur_stack[-1].upper_bound]
+ # Update the progress graph
+ progress_graph.update(buildModeStr(graph, cur_vertex), cur_mode_stack.bloated_tube[0],
+ cur_mode_stack.remain_time)
+
+ if len(cur_successors) == 0:
+ # If there is not successor
+ # Calculate the current bloated tube without considering the guard
+ cur_bloated_tube = calc_bloated_tube(cur_label,
+ cur_initial,
+ cur_remain_time,
+ sim_function,
+ params.bloatingMethod,
+ params.kvalue,
+ userConfig.SIMTRACENUM,
+ )
+
+ candidate_tube = []
+ shortest_time = float("inf")
+ shortest_tube = None
+
+ for cur_successor in cur_successors:
+ edge_id = graph.get_eid(cur_vertex, cur_successor)
+ cur_guard_str = graph.es[edge_id]['guards']
+ cur_reset_str = graph.es[edge_id]['resets']
+ # Calculate the current bloated tube with guard involved
+ # Pre-check the simulation trace so we can get better bloated result
+ cur_bloated_tube = calc_bloated_tube(cur_label,
+ cur_initial,
+ cur_remain_time,
+ sim_function,
+ params.bloatingMethod,
+ params.kvalue,
+ userConfig.SIMTRACENUM,
+ guard_checker=guard,
+ guard_str=cur_guard_str,
+ )
+
+ # Use the guard to calculate the next initial set
+ # TODO: Made this function return multiple next_init
+ next_init, truncated_result, transit_time = guard.guard_reachtube(
+ cur_bloated_tube,
+ cur_guard_str,
+ )
+
+ if next_init is None:
+ continue
+
+ # Reset the next initial set
+ # TODO: Made this function handle multiple next_init
+ next_init = reset.reset_set(cur_reset_str, next_init[0], next_init[1])
+
+ # Build next mode stack
+ next_mode_stack = InitialSetStack(
+ cur_successor,
+ userConfig.CHILDREFINETHRES,
+ cur_remain_time - transit_time,
+ start_time=transit_time+cur_mode_stack.start_time
+ )
+ next_mode_stack.parent = cur_mode_stack
+ # TODO: Append all next_init into the next_mode_stack
+ next_mode_stack.stack.append(InitialSet(next_init[0], next_init[1]))
+ next_mode_stack.bloated_tube.append(
+ cur_mode_stack.bloated_tube[0] + '->' + buildModeStr(graph, cur_successor))
+ cur_stack[-1].child[cur_successor] = next_mode_stack
+ if len(truncated_result) > len(candidate_tube):
+ candidate_tube = truncated_result
+
+ # In case of must transition
+ # We need to record shortest tube
+ # As shortest tube is the tube invoke transition
+ if truncated_result[-1][0] < shortest_time:
+ shortest_time = truncated_result[-1][0]
+ shortest_tube = truncated_result
+
+ # Handle must transition
+ if params.deterministic and len(cur_stack[-1].child) > 0:
+ next_modes_info = []
+ for next_mode in cur_stack[-1].child:
+ next_modes_info.append((cur_stack[-1].child[next_mode].remain_time, next_mode))
+ # This mode gets transit first, only keep this mode
+ max_remain_time, max_time_mode = max(next_modes_info)
+ # Pop other modes because of deterministic system
+ for _, next_mode in next_modes_info:
+ if next_mode == max_time_mode:
+ continue
+ cur_stack[-1].child.pop(next_mode)
+ candidate_tube = shortest_tube
+ print("Handle deterministic system, next mode", graph.vs[list(cur_stack[-1].child.keys())[0]]['label'])
+
+ if not candidate_tube:
+ candidate_tube = cur_bloated_tube
+
+ for i in range(len(candidate_tube)):
+ candidate_tube[i][0] += cur_mode_stack.start_time
+
+ # Check the safety for current bloated tube
+ safety = checker.check_reachtube(candidate_tube, cur_label)
+ if safety == UNSAFE:
+ print("System is not safe in Mode ", cur_label)
+ # Start back Tracking from this point and print tube to a file
+ # push current unsafe_tube to unsafe tube holder
+ unsafe_tube = [cur_mode_stack.bloated_tube[0]] + candidate_tube
+ while cur_mode_stack.parent is not None:
+ prev_mode_stack = cur_mode_stack.parent
+ unsafe_tube = [prev_mode_stack.bloated_tube[0]] + prev_mode_stack.stack[-1].bloated_tube \
+ + unsafe_tube
+ cur_mode_stack = prev_mode_stack
+ print('simulation time', sim_end_time - t_start)
+ print('verification time', time.time() - sim_end_time)
+ print('refine time', refine_counter)
+ writeReachTubeFile(unsafe_tube, UNSAFEFILENAME)
+ ret_reach = ReachTube(cur_mode_stack.bloated_tube, params.variables, params.vertex)
+ return "UNSAFE", ret_reach
+
+ elif safety == UNKNOWN:
+ # Refine the current initial set
+ print(cur_mode_stack.mode, "check bloated tube unknown")
+ discard_initial = cur_mode_stack.stack.pop()
+ init_one, init_two = discard_initial.refine()
+ cur_mode_stack.stack.append(init_one)
+ cur_mode_stack.stack.append(init_two)
+ refine_counter += 1
+
+ elif safety == SAFE:
+ print("Mode", cur_mode_stack.mode, "check bloated tube safe")
+ if cur_mode_stack.stack[-1].child:
+ cur_mode_stack.stack[-1].bloated_tube += candidate_tube
+ next_mode, next_mode_stack = cur_mode_stack.stack[-1].child.popitem()
+ cur_mode_stack = next_mode_stack
+ print("Child exist in cur mode inital", cur_mode_stack.mode, "is cur_mode_stack Now")
+ else:
+ cur_mode_stack.bloated_tube += candidate_tube
+ cur_mode_stack.stack.pop()
+ print("No child exist in current initial, pop")
+
+ if cur_mode_stack.parent is None:
+ # We are at head now
+ if backward_flag == SAFE:
+ # All the nodes are safe
+ print("System is Safe!")
+ print("refine time", refine_counter)
+ writeReachTubeFile(cur_mode_stack.bloated_tube, REACHTUBEOUTPUT)
+ ret_reach = ReachTube(cur_mode_stack.bloated_tube, params.variables, params.vertex)
+ print('simulation time', sim_end_time - t_start)
+ print('verification time', time.time() - sim_end_time)
+ return "SAFE", ret_reach
+ elif backward_flag == UNKNOWN:
+ print("Hit refine threshold, system halt, result unknown")
+ print('simulation time', sim_end_time - t_start)
+ print('verification time', time.time() - sim_end_time)
+ return "UNKNOWN", None
+ else:
+ if backward_flag == SAFE:
+ prev_mode_stack = cur_mode_stack.parent
+ prev_mode_stack.stack[-1].bloated_tube += cur_mode_stack.bloated_tube
+ print('back flag safe from', cur_mode_stack.mode, 'to', prev_mode_stack.mode)
+ if len(prev_mode_stack.stack[-1].child) == 0:
+ # There is no next mode from this initial set
+ prev_mode_stack.bloated_tube += prev_mode_stack.stack[-1].bloated_tube
+ prev_mode_stack.stack.pop()
+ cur_mode_stack = prev_mode_stack
+ print("No child in prev mode initial, pop,", prev_mode_stack.mode, "is cur_mode_stack Now")
+ else:
+ # There is another mode transition from this initial set
+ next_mode, next_mode_stack = prev_mode_stack.stack[-1].child.popitem()
+ cur_mode_stack = next_mode_stack
+ print("Child exist in prev mode inital", next_mode_stack.mode, "is cur_mode_stack Now")
+ elif backward_flag == UNKNOWN:
+ prev_mode_stack = cur_mode_stack.parent
+ print('back flag unknown from', cur_mode_stack.mode, 'to', prev_mode_stack.mode)
+ discard_initial = prev_mode_stack.stack.pop()
+ init_one, init_two = discard_initial.refine()
+ prev_mode_stack.stack.append(init_one)
+ prev_mode_stack.stack.append(init_two)
+ cur_mode_stack = prev_mode_stack
+ refine_counter += 1
+
+
+def graph_search(data, sim_function, param_config=None):
+ """
+ DryVR controller synthesis algorithm.
+ It does the controller synthesis and print out the search result.
+ tube and transition graph will be stored in output folder if algorithm finds one
+
+ Args:
+ data (dict): dictionary that contains params for the input file
+ sim_function (function): black-box simulation function
+ param_config (dict or None): user-specified configuration
+
+ Returns:
+ None
+
+ """
+ if param_config is None:
+ param_config = {}
+ # There are some fields can be config by user,
+ # If user specified these fields in paramConfig,
+ # overload these parameters to userConfig
+ overloadConfig(userConfig, param_config)
+ # Parse the input json file and read out the parameters
+ params = parseRrtInputFile(data)
+ # Construct objects
+ checker = UniformChecker(params.unsafeSet, params.variables)
+ goal_set_checker = GoalChecker(params.goalSet, params.variables)
+ distance_checker = DistChecker(params.goal, params.variables)
+ # Read the important param
+ available_modes = params.modes
+ start_modes = params.modes
+ remain_time = params.timeHorizon
+ min_time_thres = params.minTimeThres
+
+ # Set goal reach flag to False
+ # Once the flag is set to True, It means we find a transition Graph
+ goal_reached = False
+
+ # Build the initial mode stack
+ # Current Method is ugly, we need to get rid of the initial Mode for GraphSearch
+ # It helps us to achieve the full automate search
+ # TODO Get rid of the initial Mode thing
+ random.shuffle(start_modes)
+ dummy_node = GraphSearchNode("start", remain_time, min_time_thres, 0)
+ for mode in start_modes:
+ dummy_node.children[mode] = GraphSearchNode(mode, remain_time, min_time_thres, dummy_node.level + 1)
+ dummy_node.children[mode].parent = dummy_node
+ dummy_node.children[mode].initial = (params.initialSet[0], params.initialSet[1])
+
+ cur_mode_stack = dummy_node.children[start_modes[0]]
+ dummy_node.visited.add(start_modes[0])
+
+ t_start = time.time()
+ while True:
+
+ if not cur_mode_stack:
+ break
+
+ if cur_mode_stack == dummy_node:
+ start_modes.pop(0)
+ if len(start_modes) == 0:
+ break
+
+ cur_mode_stack = dummy_node.children[start_modes[0]]
+ dummy_node.visited.add(start_modes[0])
+ continue
+
+ print(str(cur_mode_stack))
+
+ # Keep check the remain time, if the remain time is less than minTime
+ # It means it is impossible to stay in one mode more than minTime
+ # Therefore, we have to go back to parents
+ if cur_mode_stack.remain_time < min_time_thres:
+ print("Back to previous mode because we cannot stay longer than the min time thres")
+ cur_mode_stack = cur_mode_stack.parent
+ continue
+
+ # If we have visited all available modes
+ # We should select a new candidate point to proceed
+ # If there is no candidates available,
+ # Then we can say current node is not valid and go back to parent
+ if len(cur_mode_stack.visited) == len(available_modes):
+ if len(cur_mode_stack.candidates) < 2:
+ print("Back to previous mode because we do not have any other modes to pick")
+ cur_mode_stack = cur_mode_stack.parent
+ # If the tried all possible cases with no luck to find path
+ if not cur_mode_stack:
+ break
+ continue
+ else:
+ print("Pick a new point from candidates")
+ cur_mode_stack.candidates.pop(0)
+ cur_mode_stack.visited = set()
+ cur_mode_stack.children = {}
+ continue
+
+ # Generate bloated tube if we haven't done so
+ if not cur_mode_stack.bloated_tube:
+ print("no bloated tube find in this mode, generate one")
+ cur_bloated_tube = calc_bloated_tube(
+ cur_mode_stack.mode,
+ cur_mode_stack.initial,
+ cur_mode_stack.remain_time,
+ sim_function,
+ params.bloatingMethod,
+ params.kvalue,
+ userConfig.SIMTRACENUM
+ )
+
+ # Cut the bloated tube once it intersect with the unsafe set
+ cur_bloated_tube = checker.cut_tube_till_unsafe(cur_bloated_tube)
+
+ # If the tube time horizon is less than minTime, it means
+ # we cannot stay in this mode for min thres time, back to the parent node
+ if not cur_bloated_tube or cur_bloated_tube[-1][0] < min_time_thres:
+ print("bloated tube is not long enough, discard the mode")
+ cur_mode_stack = cur_mode_stack.parent
+ continue
+ cur_mode_stack.bloated_tube = cur_bloated_tube
+
+ # Generate candidates points for next node
+ random_sections = cur_mode_stack.random_picker(userConfig.RANDSECTIONNUM)
+
+ if not random_sections:
+ print("bloated tube is not long enough, discard the mode")
+ cur_mode_stack = cur_mode_stack.parent
+ continue
+
+ # Sort random points based on the distance to the goal set
+ random_sections.sort(key=lambda x: distance_checker.calc_distance(x[0], x[1]))
+ cur_mode_stack.candidates = random_sections
+ print("Generate new bloated tube and candidate, with candidates length", len(cur_mode_stack.candidates))
+
+ # Check if the current tube reaches goal
+ result, tube = goal_set_checker.goal_reachtube(cur_bloated_tube)
+ if result:
+ cur_mode_stack.bloated_tube = tube
+ goal_reached = True
+ break
+
+ # We have visited all next mode we have, generate some thing new
+ # This is actually not necessary, just shuffle all modes would be enough
+ # There should not be RANDMODENUM things since it does not make any difference
+ # Anyway, for each candidate point, we will try to visit all modes eventually
+ # Therefore, using RANDMODENUM to get some random modes visit first is useless
+ # TODO, fix this part
+ if len(cur_mode_stack.visited) == len(cur_mode_stack.children):
+ # leftMode = set(available_modes) - set(cur_mode_stack.children.keys())
+ # random_modes = random.sample(leftMode, min(len(leftMode), RANDMODENUM))
+ random_modes = available_modes
+ random.shuffle(random_modes)
+
+ random_sections = cur_mode_stack.random_picker(userConfig.RANDSECTIONNUM)
+ for mode in random_modes:
+ candidate = cur_mode_stack.candidates[0]
+ cur_mode_stack.children[mode] = GraphSearchNode(mode, cur_mode_stack.remain_time - candidate[1][0],
+ min_time_thres, cur_mode_stack.level + 1)
+ cur_mode_stack.children[mode].initial = (candidate[0][1:], candidate[1][1:])
+ cur_mode_stack.children[mode].parent = cur_mode_stack
+
+ # Random visit a candidate that is not visited before
+ for key in cur_mode_stack.children:
+ if key not in cur_mode_stack.visited:
+ break
+
+ print("transit point is", cur_mode_stack.candidates[0])
+ cur_mode_stack.visited.add(key)
+ cur_mode_stack = cur_mode_stack.children[key]
+
+ # Back track to print out trace
+ print("RRT run time", time.time() - t_start)
+ if goal_reached:
+ print("goal reached")
+ traces = []
+ modes = []
+ while cur_mode_stack:
+ modes.append(cur_mode_stack.mode)
+ if not cur_mode_stack.candidates:
+ traces.append([t for t in cur_mode_stack.bloated_tube])
+ else:
+ # Cut the trace till candidate
+ temp = []
+ for t in cur_mode_stack.bloated_tube:
+ if t == cur_mode_stack.candidates[0][0]:
+ temp.append(cur_mode_stack.candidates[0][0])
+ temp.append(cur_mode_stack.candidates[0][1])
+ break
+ else:
+ temp.append(t)
+ traces.append(temp)
+ if cur_mode_stack.parent != dummy_node:
+ cur_mode_stack = cur_mode_stack.parent
+ else:
+ break
+ # Reorganize the content in modes list for plotter use
+ modes = modes[::-1]
+ traces = traces[::-1]
+ build_rrt_graph(modes, traces, isIpynb())
+ for i in range(1, len(modes)):
+ modes[i] = modes[i - 1] + '->' + modes[i]
+
+ writeRrtResultFile(modes, traces, RRTOUTPUT)
+ else:
+ print("could not find graph")
diff --git a/ourtool/dryvr/core/goalchecker.py b/ourtool/dryvr/core/goalchecker.py
new file mode 100644
index 0000000000000000000000000000000000000000..ba84296c159c2d102123179eb75a5cb640e71dcc
--- /dev/null
+++ b/ourtool/dryvr/core/goalchecker.py
@@ -0,0 +1,107 @@
+"""
+This file contains uniform checker class for DryVR
+"""
+
+from src.common.utils import handleReplace, neg
+from z3 import *
+
+
+class GoalChecker:
+ """
+ This is class to check if the goal set is reached
+ by reach tube
+ """
+
+ def __init__(self, goal, variables):
+ """
+ Goal checker class initialization function.
+
+ Args:
+ goal (str): a str describing the goal set.
+ For example: "And(x_1>=19.5, x_1<=20.5, x_2>=-1.0, x_2<=1.0)"
+ variables (list): list of variable name
+ """
+ self.varDic = {'t': Real('t')}
+ self.variables = variables
+ for var in variables:
+ self.varDic[var] = Real(var)
+
+ goal = handleReplace(goal, list(self.varDic.keys()))
+
+ self.intersectChecker = Solver()
+ self.containChecker = Solver()
+
+ self.intersectChecker.add(eval(goal))
+ self.containChecker.add(eval(neg(goal)))
+
+ def goal_reachtube(self, tube):
+ """
+ Check if the reach tube satisfied the goal
+
+ Args:
+ tube (list): the reach tube.
+
+ Returns:
+ A bool indicates if the goal is reached
+ The truncated tube if the goal is reached, otherwise the whole tube
+ """
+ for i in range(0, len(tube), 2):
+ lower = tube[i]
+ upper = tube[i + 1]
+ if self._check_intersection(lower, upper):
+ if self._check_containment(lower, upper):
+ return True, tube[:i + 2]
+ return False, tube
+
+ def _check_intersection(self, lower, upper):
+ """
+ Check if the goal set intersect with the current set
+ #FIXME Maybe this is not necessary since we only want to check
+ the fully contained case
+ Bolun 02/13/2018
+
+ Args:
+ lower (list): the list represent the set's lowerbound.
+ upper (list): the list represent the set's upperbound.
+
+ Returns:
+ A bool indicates if the set intersect with the goal set
+ """
+ cur_solver = self.intersectChecker
+ cur_solver.push()
+ cur_solver.add(self.varDic["t"] >= lower[0])
+ cur_solver.add(self.varDic["t"] <= upper[0])
+ for i in range(1, len(lower)):
+ cur_solver.add(self.varDic[self.variables[i - 1]] >= lower[i])
+ cur_solver.add(self.varDic[self.variables[i - 1]] <= upper[i])
+ if cur_solver.check() == sat:
+ cur_solver.pop()
+ return True
+ else:
+ cur_solver.pop()
+ return False
+
+ def _check_containment(self, lower, upper):
+ """
+ Check if the current set contained in goal set.
+
+ Args:
+ lower (list): the list represent the set's lowerbound.
+ upper (list): the list represent the set's upperbound.
+
+ Returns:
+ A bool indicates if the set if contained in the goal set
+ """
+ cur_solver = self.containChecker
+ cur_solver.push()
+ cur_solver.add(self.varDic["t"] >= lower[0])
+ cur_solver.add(self.varDic["t"] <= upper[0])
+ for i in range(1, len(lower)):
+ cur_solver.add(self.varDic[self.variables[i - 1]] >= lower[i])
+ cur_solver.add(self.varDic[self.variables[i - 1]] <= upper[i])
+ if cur_solver.check() == sat:
+ cur_solver.pop()
+ return False
+ else:
+ cur_solver.pop()
+ return True
diff --git a/ourtool/dryvr/core/graph.py b/ourtool/dryvr/core/graph.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4770ff2001d436b7c4e0c9408812cea809c2fa1
--- /dev/null
+++ b/ourtool/dryvr/core/graph.py
@@ -0,0 +1,80 @@
+"""
+This file contains graph class for DryVR
+"""
+# import matplotlib.pyplot as plt
+import networkx as nx
+
+
+class Graph:
+ """
+ This is class to plot the progress graph for dryvr's
+ function, it is supposed to display the graph in jupyter
+ notebook and update the graph as dryvr is running
+ """
+
+ def __init__(self, params, is_ipynb):
+ """
+ Guard checker class initialization function.
+
+ Args:
+ params (obj): An object contains the parameter
+ is_ipynb (bool): check if the code is running on ipython or not
+ """
+ self._is_ipynb = is_ipynb
+ if not is_ipynb:
+ return
+ vertex = []
+ # Build unique identifier for a vertex and mode name
+ for idx, v in enumerate(params.vertex):
+ vertex.append(v + "," + str(idx))
+
+ edges = params.edge
+ self.edgeList = []
+ for e in edges:
+ self.edgeList.append((vertex[e[0]], vertex[e[1]]))
+ # Initialize the plot
+ # self.fig = plt.figure()
+ # self.ax = self.fig.add_subplot(111)
+
+ # Initialize the graph
+ self.G = nx.DiGraph()
+ self.G.add_edges_from(self.edgeList)
+ self.pos = nx.spring_layout(self.G)
+ self.colors = ['green'] * len(self.G.nodes())
+ # self.fig.suptitle('', fontsize=10)
+ # Draw the graph when initialize
+ # self.draw()
+ # plt.show()
+
+ def draw(self):
+ """
+ Draw the white-box transition graph.
+ """
+
+ nx.draw_networkx_labels(self.G, self.pos)
+ options = {
+ 'node_color': self.colors,
+ 'node_size': 1000,
+ 'cmap': plt.get_cmap('jet'),
+ 'arrowstyle': '-|>',
+ 'arrowsize': 50,
+ }
+ nx.draw_networkx(self.G, self.pos, arrows=True, **options)
+ self.fig.canvas.draw()
+
+ def update(self, cur_node, title, remain_time):
+ """
+ update the graph
+ Args:
+ cur_node (str): current vertex dryvr is verifying
+ title (str): current transition path as the title
+ remain_time (float): remaining time
+ """
+ if not self._is_ipynb:
+ return
+ self.ax.clear()
+ self.colors = ['green'] * len(self.G.nodes())
+ self.colors[list(self.G.nodes()).index(cur_node)] = 'red'
+ self.fig.suptitle(title, fontsize=10)
+ self.ax.set_title("remain:" + str(remain_time))
+ self.draw()
diff --git a/ourtool/dryvr/core/guard.py b/ourtool/dryvr/core/guard.py
new file mode 100644
index 0000000000000000000000000000000000000000..76e7e8238defec4aaa34af58a41ce1609455306f
--- /dev/null
+++ b/ourtool/dryvr/core/guard.py
@@ -0,0 +1,213 @@
+"""
+This file contains guard class for DryVR
+"""
+
+import random
+
+import sympy
+from z3 import *
+
+from src.common.utils import handleReplace
+
+
+class Guard:
+ """
+ This is class to calculate the set in the
+ reach tube that intersect with the guard
+ """
+
+ def __init__(self, variables):
+ """
+ Guard checker class initialization function.
+
+ Args:
+ variables (list): list of variable name
+ """
+ self.varDic = {'t': Real('t')}
+ self.variables = variables
+ for var in variables:
+ self.varDic[var] = Real(var)
+
+ def _build_guard(self, guard_str):
+ """
+ Build solver for current guard based on guard string
+
+ Args:
+ guard_str (str): the guard string.
+ For example:"And(v>=40-0.1*u, v-40+0.1*u<=0)"
+
+ Returns:
+ A Z3 Solver obj that check for guard.
+ A symbol index dic obj that indicates the index
+ of variables that involved in the guard.
+ """
+ cur_solver = Solver()
+ # This magic line here is because SymPy will evaluate == to be False
+ # Therefore we are not be able to get free symbols from it
+ # Thus we need to replace "==" to something else
+ sympy_guard_str = guard_str.replace("==", ">=")
+
+ symbols = list(sympy.sympify(sympy_guard_str, evaluate=False).free_symbols)
+ symbols = [str(s) for s in symbols]
+ symbols_idx = {s: self.variables.index(s) + 1 for s in symbols if s in self.variables}
+ if 't' in symbols:
+ symbols_idx['t'] = 0
+
+ guard_str = handleReplace(guard_str, list(self.varDic.keys()))
+ cur_solver.add(eval(guard_str)) # TODO use an object instead of `eval` a string
+ return cur_solver, symbols_idx
+
+ def guard_sim_trace(self, trace, guard_str):
+ """
+ Check the guard for simulation trace.
+ Note we treat the simulation trace as the set as well.
+ Consider we have a simulation trace as following
+ [0.0, 1.0, 1.1]
+ [0.1, 1.02, 1.13]
+ [0.2, 1.05, 1.14]
+ ...
+ We can build set like
+ lower_bound: [0.0, 1.0, 1.1]
+ upper_bound: [0.1, 1.02, 1.13]
+
+ lower_bound: [0.1, 1.02, 1.13]
+ upper_bound: [0.2, 1.05, 1.14]
+ And check guard for these set. This is because if the guard
+ is too small, it is likely for simulation point ignored the guard.
+ For example:
+ . . . . |guard| . . .
+ In this case, the guard gets ignored
+
+ Args:
+ trace (list): the simulation trace
+ guard_str (str): the guard string.
+ For example:"And(v>=40-0.1*u, v-40+0.1*u<=0)"
+
+ Returns:
+ A initial point for next mode,
+ The truncated simulation trace
+ """
+ if not guard_str:
+ return None, trace
+ cur_solver, symbols = self._build_guard(guard_str)
+ guard_set = {}
+
+ for idx in range(len(trace) - 1):
+ lower = trace[idx]
+ upper = trace[idx + 1]
+ cur_solver.push()
+ for symbol in symbols:
+ cur_solver.add(self.varDic[symbol] >= min(lower[symbols[symbol]], upper[symbols[symbol]]))
+ cur_solver.add(self.varDic[symbol] <= max(lower[symbols[symbol]], upper[symbols[symbol]]))
+ if cur_solver.check() == sat:
+ cur_solver.pop()
+ guard_set[idx] = upper
+ else:
+ cur_solver.pop()
+ if guard_set:
+ # Guard set is not empty, randomly pick one and return
+ # idx, point = random.choice(list(guard_set.items()))
+ idx, point = list(guard_set.items())[0]
+ # Return the initial point for next mode, and truncated trace
+ return point[1:], trace[:idx + 1]
+
+ if guard_set:
+ # Guard set is not empty, randomly pick one and return
+ # idx, point = random.choice(list(guard_set.items()))
+ idx, point = list(guard_set.items())[0]
+ # Return the initial point for next mode, and truncated trace
+ return point[1:], trace[:idx + 1]
+
+ # No guard hits for current tube
+ return None, trace
+
+ def guard_sim_trace_time(self, trace, guard_str):
+ """
+ Return the length of the truncated traces
+
+ Args:
+ trace (list): the simulation trace
+ guard_str (str): the guard string.
+ For example:"And(v>=40-0.1*u, v-40+0.1*u<=0)"
+
+ Returns:
+ the length of the truncated traces.
+ """
+ next_init, trace = self.guard_sim_trace(trace, guard_str)
+ return len(trace)
+
+ def guard_reachtube(self, tube, guard_str):
+ """
+ Check the guard intersection of the reach tube
+
+
+ Args:
+ tube (list): the reach tube
+ guard_str (str): the guard string.
+ For example:"And(v>=40-0.1*u, v-40+0.1*u<=0)"
+
+ Returns:
+ Next mode initial set represent as [upper_bound, lower_bound],
+ Truncated tube before the guard,
+ The time when elapsed in current mode.
+
+ """
+ if not guard_str:
+ return None, tube
+
+ cur_solver, symbols = self._build_guard(guard_str)
+ guard_set_lower = []
+ guard_set_upper = []
+ for i in range(0, len(tube), 2):
+ cur_solver.push()
+ lower_bound = tube[i]
+ upper_bound = tube[i + 1]
+ for symbol in symbols:
+ cur_solver.add(self.varDic[symbol] >= lower_bound[symbols[symbol]])
+ cur_solver.add(self.varDic[symbol] <= upper_bound[symbols[symbol]])
+ if cur_solver.check() == sat:
+ # The reachtube hits the guard
+ cur_solver.pop()
+ guard_set_lower.append(lower_bound)
+ guard_set_upper.append(upper_bound)
+
+ # TODO: If the reachtube completely fall inside guard, break
+ tmp_solver = Solver()
+ tmp_solver.add(Not(cur_solver.assertions()[0]))
+ for symbol in symbols:
+ tmp_solver.add(self.varDic[symbol] >= lower_bound[symbols[symbol]])
+ tmp_solver.add(self.varDic[symbol] <= upper_bound[symbols[symbol]])
+ if tmp_solver.check() == unsat:
+ print("Full intersect, break")
+ break
+ else:
+ cur_solver.pop()
+ if guard_set_lower:
+ # Guard set is not empty, build the next initial set and return
+ # At some point we might further reduce the initial set for next mode
+ init_lower = guard_set_lower[0][1:]
+ init_upper = guard_set_upper[0][1:]
+ for j in range(1, len(guard_set_lower)):
+ for k in range(1, len(guard_set_lower[0])):
+ init_lower[k - 1] = min(init_lower[k - 1], guard_set_lower[j][k])
+ init_upper[k - 1] = max(init_upper[k - 1], guard_set_upper[j][k])
+ # Return next initial Set, the result tube, and the true transit time
+ return [init_lower, init_upper], tube[:i], guard_set_lower[0][0]
+
+ # Construct the guard if all later trace sat the guard condition
+ if guard_set_lower:
+ # Guard set is not empty, build the next initial set and return
+ # At some point we might further reduce the initial set for next mode
+ init_lower = guard_set_lower[0][1:]
+ init_upper = guard_set_upper[0][1:]
+ for j in range(1, len(guard_set_lower)):
+ for k in range(1, len(guard_set_lower[0])):
+ init_lower[k - 1] = min(init_lower[k - 1], guard_set_lower[j][k])
+ init_upper[k - 1] = max(init_upper[k - 1], guard_set_upper[j][k])
+ # TODO: Treat tau as a special clock variable that don't have variation
+ # init_upper[0] = init_lower[0]
+
+ # Return next initial Set, the result tube, and the true transit time
+ return [init_lower, init_upper], tube[:i], guard_set_lower[0][0]
+
+ return None, tube, tube[-1][0]
diff --git a/ourtool/dryvr/core/initialset.py b/ourtool/dryvr/core/initialset.py
new file mode 100644
index 0000000000000000000000000000000000000000..0bc2a0a4f00648f93e859d83b6d8c0e172b5d69b
--- /dev/null
+++ b/ourtool/dryvr/core/initialset.py
@@ -0,0 +1,69 @@
+"""
+This file contains initial set class for DryVR
+"""
+
+
+class InitialSet:
+ """
+ This is class to represent the initial set
+ """
+
+ def __init__(self, lower, upper):
+ """
+ Initial set class initialization function.
+
+ Args:
+ lower (list): lower bound of the initial set
+ upper (list): upper bound of the initial set
+ """
+
+ self.lower_bound = lower
+ self.upper_bound = upper
+ self.delta = [(upper[i] - lower[i]) / 2.0 for i in range(len(upper))]
+ # Child point points to children InitialSetStack obj
+ # This it how it works
+ # Since a initial set can generate a reach tube that intersect
+ # with different guards
+ # So there can be multiple children pointers
+ # Therefore this is a dictionary obj
+ # self.child["MODEA"] = InitialSetStack for MODEA
+ self.child = {}
+ self.bloated_tube = []
+
+ def refine(self):
+ """
+ This function refine the current initial set into two smaller set
+
+ Returns:
+ two refined initial set
+
+ """
+ # Refine the initial set into two smaller set
+ # based on index with largest delta
+ idx = self.delta.index(max(self.delta))
+ # Construct first smaller initial set
+ init_set_one_ub = list(self.upper_bound)
+ init_set_one_lb = list(self.lower_bound)
+ init_set_one_lb[idx] += self.delta[idx]
+ # Construct second smaller initial set
+ init_set_two_ub = list(self.upper_bound)
+ init_set_two_lb = list(self.lower_bound)
+ init_set_two_ub[idx] -= self.delta[idx]
+
+ return (
+ InitialSet(init_set_one_lb, init_set_one_ub),
+ InitialSet(init_set_two_lb, init_set_two_ub),
+ )
+
+ def __str__(self):
+ """
+ Build string representation for the initial set
+
+ Returns:
+ A string describes the initial set
+ """
+ ret = ""
+ ret += "Lower Bound: " + str(self.lower_bound) + "\n"
+ ret += "Upper Bound: " + str(self.upper_bound) + "\n"
+ ret += "Delta: " + str(self.delta) + "\n"
+ return ret
diff --git a/ourtool/dryvr/core/initialsetstack.py b/ourtool/dryvr/core/initialsetstack.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3474a07eb00a97e1fe734eaccc323ea410d2500
--- /dev/null
+++ b/ourtool/dryvr/core/initialsetstack.py
@@ -0,0 +1,132 @@
+"""
+This file contains initial set class stack for DryVR
+"""
+import random
+
+
+class InitialSetStack:
+ """
+ This is class is for list of initial sets for a single node
+ """
+
+ def __init__(self, mode, threshold, remain_time, start_time):
+ """
+ Initial set class initialization function.
+
+ Args:
+ mode (str): the mode name for current node
+ threshold (int): threshold for refinement
+ remain_time (float): remaining time of the current node
+ """
+
+ # Mode name for the current mode
+ self.mode = mode
+ # Threshold (int) for the refinement
+ self.threshold = threshold
+ # Pointer point to the parent node
+ # If you wonder where child pointers are
+ # It's actually in the InitialSet class.........
+ self.parent = None
+ # A stack for InitialSet object
+ self.stack = []
+ self.remain_time = remain_time
+ # Stores bloated tube result for all tubes start from initial sets in stack
+ self.bloated_tube = []
+ self.start_time = start_time
+
+ def is_valid(self):
+ """
+ Check if current node is still valid or not.
+ (if number of initial sets in stack are more than threshold)
+
+ Returns:
+ A bool to check if the stack is valid or not valid
+ """
+ return len(self.stack) <= self.threshold
+
+ def __str__(self):
+ """
+ Build string representation for the initial set stack
+
+ Args:
+ None
+ Returns:
+ A string describes the stack
+ """
+ ret = "===============================================\n"
+ ret += "Mode: " + str(self.mode) + "\n"
+ ret += "stack size: " + str(len(self.stack)) + "\n"
+ ret += "remainTime: " + str(self.remain_time) + "\n"
+ return ret
+
+
+class GraphSearchNode:
+ """
+ This is class for graph search node
+ Contains some helpful stuff in graph search tree
+ """
+
+ def __init__(self, mode, remain_time, min_time_thres, level):
+ """
+ GraphSearchNode class initialization function.
+
+ Args:
+ mode (str): the mode name for current node
+ remain_time (float): remaining time of the current node
+ min_time_thres (float): minimal time should stay in this mode
+ level (int): tree depth
+ """
+ self.mode = mode
+ # Pointer point to parent node
+ self.parent = None
+ # Initial set
+ self.initial = None
+ self.bloated_tube = []
+ self.remain_time = remain_time
+ self._min_time_thres = min_time_thres
+ # Pointer point to children
+ self.children = {}
+ # keep track which child is visited
+ self.visited = set()
+ # select number of candidates initial set for children
+ self.candidates = []
+ self.level = level
+
+ def random_picker(self, k):
+ """
+ Randomly pick k candidates initial set
+
+ Args:
+ k (int): number of candidates
+ Returns:
+ A list of candidate initial sets
+ """
+ i = 0
+ ret = []
+ while i < len(self.bloated_tube):
+ if self.bloated_tube[i][0] >= self._min_time_thres:
+ ret.append((self.bloated_tube[i], self.bloated_tube[i + 1]))
+ i += 2
+
+ random.shuffle(ret)
+ return ret[:k]
+
+ def __str__(self):
+ """
+ Build string representation for the Graph Search Node
+
+ Returns:
+ A string describes the stack
+ """
+ ret = ""
+ ret += "Level: " + str(self.level) + "\n"
+ ret += "Mode: " + str(self.mode) + "\n"
+ ret += "Initial: " + str(self.initial) + "\n"
+ ret += "visited: " + str(self.visited) + "\n"
+ ret += "num candidates: " + str(len(self.candidates)) + "\n"
+ ret += "remaining time: " + str(self.remain_time) + "\n"
+ if self.bloated_tube:
+ ret += "bloatedTube: True" + "\n"
+ else:
+ ret += "bloatedTube: False" + "\n"
+ return ret
diff --git a/ourtool/dryvr/core/reachtube.py b/ourtool/dryvr/core/reachtube.py
new file mode 100644
index 0000000000000000000000000000000000000000..f985f43778cc4fe2b8ef7f165025a6394e32a17b
--- /dev/null
+++ b/ourtool/dryvr/core/reachtube.py
@@ -0,0 +1,88 @@
+"""
+This file contains reach tube class for DryVR
+"""
+
+import six
+
+
+class ReachTube:
+ """
+ This is class is an object for reach tube
+ Ideally it should support to fetch reachtube by mode and variable name
+ And it should allow users to plot the reach tube in different ways
+ """
+
+ def __init__(self, tube, variables, modes):
+ """
+ ReachTube class initialization function.
+
+ Args:
+ tube (list): raw reach tube (that used to print to file)
+ variables (list): list of variables in the reach tube
+ modes (list): list of modes in the reach ReachTube
+ """
+ self._tube = tube
+ self._variables = variables
+ self._modes = modes
+
+ # Build the raw string representation so user can print it
+
+ self.raw = ""
+ for line in tube:
+ if isinstance(line, str):
+ self.raw += line + "\n"
+ else:
+ self.raw += " ".join(map(str, line)) + '\n'
+
+ # Build dictionary object so you can easily pull out part of the list
+ self._tube_dict = {}
+ for mode in modes:
+ self._tube_dict[mode] = {}
+ for var in variables + ["t"]:
+ self._tube_dict[mode][var] = []
+
+ cur_mode = ""
+ for line in tube:
+ if isinstance(line, six.string_types):
+ cur_mode = line.split('->')[-1].split(',')[0] # Get current mode name
+ for var in ['t'] + self._variables:
+ self._tube_dict[cur_mode][var].append(line)
+ else:
+ for var, val in zip(['t'] + self._variables, line):
+ self._tube_dict[cur_mode][var].append(val)
+
+ def __str__(self):
+ """
+ print the raw tube
+ """
+ return str(self.raw)
+
+ def filter(self, mode=None, variable=None, contain_label=True):
+ """
+ This is a filter function that allows you to select
+ Args:
+ mode (str, list): single mode name or list of mode name
+ variable (str, list): single variable or list of variables
+ """
+ if mode is None:
+ mode = self._modes
+ if variable is None:
+ variable = ["t"] + self._variables
+
+ if isinstance(mode, str):
+ mode = [mode]
+ if isinstance(variable, str):
+ variable = [variable]
+
+ res = []
+ for m in mode:
+ temp = []
+ for i in range(len(self._tube_dict[m]["t"])):
+ if isinstance(self._tube_dict[m]["t"][i], str):
+ if contain_label:
+ temp.append([self._tube_dict[m]["t"][i]] + variable)
+ continue
+
+ temp.append([self._tube_dict[m][v][i] for v in variable])
+ res.append(temp)
+ return res
diff --git a/ourtool/dryvr/core/reset.py b/ourtool/dryvr/core/reset.py
new file mode 100644
index 0000000000000000000000000000000000000000..684b9c5e2b3af8de2ba8231d51e8d3dafdd32875
--- /dev/null
+++ b/ourtool/dryvr/core/reset.py
@@ -0,0 +1,209 @@
+"""
+This file contains reset class for DryVR
+"""
+
+import sympy
+
+from src.common.utils import randomPoint
+
+
+class Reset:
+ """
+ This is class for resetting the initial set
+ """
+
+ def __init__(self, variables):
+ """
+ Reset class initialization function.
+
+ Args:
+ variables (list): list of varibale name
+ """
+ self.variables = variables
+
+ def reset_set(self, raw_eqs_str, lower_bound, upper_bound):
+ """
+ Reset the initial set based on reset expressions
+
+ Args:
+ raw_eqs_str (str): reset expressions separated by ';'
+ lower_bound (list): lower bound of the initial set
+ upper_bound (list): upper bound of the initial set
+
+ Returns:
+ lower bound and upper bound of the initial set after reset
+ """
+ if not raw_eqs_str:
+ return lower_bound, upper_bound
+
+ raw_eqs = raw_eqs_str.split(';')
+ lb_list = []
+ ub_list = []
+ for rawEqu in raw_eqs:
+ lb, ub = self._handle_reset(rawEqu, lower_bound, upper_bound)
+ lb_list.append(lb)
+ ub_list.append(ub)
+
+ return self._merge_result(lb_list, ub_list, lower_bound, upper_bound)
+
+ def reset_point(self, raw_eqs_str, point):
+ """
+ Reset the initial point based on reset expressions
+
+ Args:
+ raw_eqs_str (str): list of reset expression
+ point (list): the initial point need to be reset
+
+ Returns:
+ a point after reset
+ """
+ if point == [] or not point:
+ return point
+ lower, upper = self.reset_set(raw_eqs_str, point, point)
+ return randomPoint(lower, upper)
+
+ @staticmethod
+ def _merge_result(lb_list, ub_list, lower_bound, upper_bound):
+ """
+ Merge the a list of reset result
+ Since we allow multiple reset per transition,
+ we get list of reset result, each result corresponding to one reset expression
+ We need to merge all reset result together
+
+ Args:
+ lb_list (list): list of reset lower bound results
+ ub_list (list): list of reset upper bound results
+ lower_bound(list): original lower bound
+ upper_bound(list): original upper bound
+
+ Returns:
+ Upper bound and lower bound after merge the reset result
+ """
+ ret_lb = list(lower_bound)
+ ret_ub = list(upper_bound)
+
+ for i in range(len(lb_list)):
+ cur_lb = lb_list[i]
+ cur_ub = ub_list[i]
+ for j in range(len(cur_lb)):
+ if cur_lb[j] != lower_bound[j]:
+ ret_lb[j] = cur_lb[j]
+ if cur_ub[j] != upper_bound[j]:
+ ret_ub[j] = cur_ub[j]
+ return ret_lb, ret_ub
+
+ def _build_all_combo(self, symbols, lower_bound, upper_bound):
+ """
+ This function allows us to build all combination given symbols
+ For example, if we have a 2-dimension set for dim A and B.
+ symbols = [A,B]
+ lowerBound = [1.0, 2.0]
+ upperBound = [3.0, 4.0]
+ Then the result should be all possible combination of the value of A and B
+ result:
+ [[1.0, 2.0], [3.0, 4.0], [3.0, 2.0], [1.0, 4.0]]
+
+ Args:
+ symbols (list): symbols we use to create combo
+ lower_bound (list): lower bound of the set
+ upper_bound (list): upper bound of the set
+
+ Returns:
+ List of combination value
+ """
+ if not symbols:
+ return []
+
+ cur_symbol = str(symbols[0])
+ idx = self.variables.index(cur_symbol)
+ lo = lower_bound[idx]
+ up = upper_bound[idx]
+ ret = []
+ next_level = self._build_all_combo(symbols[1:], lower_bound, upper_bound)
+ if next_level:
+ for n in next_level:
+ ret.append(n + [(cur_symbol, lo)])
+ ret.append(n + [(cur_symbol, up)])
+ else:
+ ret.append([cur_symbol, lo])
+ ret.append([cur_symbol, up])
+ return ret
+
+ def _handle_wrapped_reset(self, raw_eq, lower_bound, upper_bound):
+ """
+ This is a function to handle reset such as V = [0, V+1]
+
+ Args:
+ raw_eq (str): reset equation
+ lower_bound (list): lower bound of the set
+ upper_bound (list): upper bound of the set
+
+ Returns:
+ Upper bound and lower bound after the reset
+ """
+ final_equ = sympy.sympify(raw_eq)
+ rhs_symbols = list(final_equ.free_symbols)
+ combos = self._build_all_combo(rhs_symbols, lower_bound, upper_bound)
+ min_reset = float('inf')
+ max_reset = float('-inf')
+ if combos:
+ for combo in combos:
+ if len(combo) == 2:
+ result = float(final_equ.subs(combo[0], combo[1]))
+ else:
+ result = float(final_equ.subs(combo))
+ min_reset = min(min_reset, float(result))
+ max_reset = max(max_reset, float(result))
+ else:
+ min_reset = float(final_equ)
+ max_reset = float(final_equ)
+ return (min_reset, max_reset)
+
+ def _handle_reset(self, raw_equ, lower_bound, upper_bound):
+ """
+ Handle the reset with single reset expression
+
+ Args:
+ raw_equ (str): reset equation
+ lower_bound (list): lower bound of the set
+ upper_bound (list): upper bound of the set
+
+ Returns:
+ Upper bound and lower bound after the reset
+ """
+ equ_split = raw_equ.split('=')
+ lhs, rhs = equ_split[0], equ_split[1]
+ target = sympy.sympify(lhs)
+ # Construct the equation
+ final_equ = sympy.sympify(rhs)
+ if not isinstance(final_equ, list):
+ rhs_symbols = list(sympy.sympify(rhs).free_symbols)
+ else:
+ rhs_symbols = None
+ # print target, rhs_symbols
+ combos = self._build_all_combo(rhs_symbols, lower_bound, upper_bound)
+ # final_equ = solve(equ, target)[0]
+
+ min_reset = float('inf')
+ max_reset = float('-inf')
+ if combos:
+ for combo in combos:
+ if len(combo) == 2:
+ result = float(final_equ.subs(combo[0], combo[1]))
+ else:
+ result = float(final_equ.subs(combo))
+ min_reset = min(min_reset, float(result))
+ max_reset = max(max_reset, float(result))
+ elif isinstance(final_equ, list):
+ min_reset = min(self._handle_wrapped_reset(final_equ[0], lower_bound, upper_bound))
+ max_reset = max(self._handle_wrapped_reset(final_equ[1], lower_bound, upper_bound))
+ else:
+ min_reset = float(final_equ)
+ max_reset = float(final_equ)
+
+ ret_lb = list(lower_bound)
+ ret_ub = list(upper_bound)
+ target_idx = self.variables.index(str(target))
+ ret_lb[target_idx] = min_reset
+ ret_ub[target_idx] = max_reset
+ return ret_lb, ret_ub
diff --git a/ourtool/dryvr/core/uniformchecker.py b/ourtool/dryvr/core/uniformchecker.py
new file mode 100644
index 0000000000000000000000000000000000000000..57641ad630bed35fed0199d015a90b7af6ba9210
--- /dev/null
+++ b/ourtool/dryvr/core/uniformchecker.py
@@ -0,0 +1,213 @@
+"""
+This file contains uniform checker class for DryVR
+"""
+import sympy
+from z3 import *
+
+from src.common.constant import *
+from src.common.utils import handleReplace, neg
+
+
+class UniformChecker:
+ """
+ This is class for check unsafe checking
+ """
+
+ def __init__(self, unsafe, variables):
+ """
+ Reset class initialization function.
+
+ Args:
+ unsafe (str): unsafe constraint
+ variables (list): list of variable name
+ """
+ self.varDic = {'t': Real('t')}
+ self._variables = variables
+ self._solver_dict = {}
+ for var in variables:
+ self.varDic[var] = Real(var)
+
+ if not unsafe:
+ return
+
+ original = unsafe
+
+ unsafe = handleReplace(unsafe, list(self.varDic.keys()))
+ unsafe_list = unsafe[1:].split('@')
+ for unsafe in unsafe_list:
+ mode, cond = unsafe.split(':')
+ self._solver_dict[mode] = [Solver(), Solver()]
+ self._solver_dict[mode][0].add(eval(cond))
+ self._solver_dict[mode][1].add(eval(neg(cond)))
+
+ unsafe_list = original[1:].split('@')
+ for unsafe in unsafe_list:
+ mode, cond = unsafe.split(':')
+ # This magic line here is because SymPy will evaluate == to be False
+ # Therefore we are not be able to get free symbols from it
+ # Thus we need to replace "==" to something else, which is >=
+ cond = cond.replace("==", ">=")
+ symbols = list(sympy.sympify(cond).free_symbols)
+ symbols = [str(s) for s in symbols]
+ symbols_idx = {s: self._variables.index(s) + 1 for s in symbols if s in self._variables}
+ if 't' in symbols:
+ symbols_idx['t'] = 0
+ self._solver_dict[mode].append(symbols_idx) # TODO Fix typing
+
+ def check_sim_trace(self, traces, mode):
+ """
+ Check the simulation trace
+
+ Args:
+ traces (list): simulation traces
+ mode (str): mode need to be checked
+
+ Returns:
+ An int for checking result SAFE = 1, UNSAFE = -1
+ """
+ if mode in self._solver_dict:
+ cur_solver = self._solver_dict[mode][0]
+ symbols = self._solver_dict[mode][2]
+ elif 'Allmode' in self._solver_dict:
+ cur_solver = self._solver_dict['Allmode'][0]
+ symbols = self._solver_dict['Allmode'][2]
+ else:
+ # Return True if we do not check this mode
+ return SAFE
+
+ for t in traces:
+ cur_solver.push()
+ for symbol in symbols:
+ cur_solver.add(self.varDic[symbol] == t[symbols[symbol]])
+
+ if cur_solver.check() == sat:
+ cur_solver.pop()
+ return UNSAFE
+ else:
+ cur_solver.pop()
+ return SAFE
+
+ def check_reachtube(self, tube, mode):
+ """
+ Check the bloated reach tube
+
+ Args:
+ tube (list): reach tube
+ mode (str): mode need to be checked
+
+ Returns:
+ An int for checking result SAFE = 1, UNSAFE = -1, UNKNOWN = 0
+ """
+ if mode not in self._solver_dict and 'Allmode' not in self._solver_dict:
+ # Return True if we do not check this mode
+ return SAFE
+
+ safe = SAFE
+ for i in range(0, len(tube), 2):
+ lower = tube[i]
+ upper = tube[i + 1]
+ if self._check_intersection(lower, upper, mode):
+ if self._check_containment(lower, upper, mode):
+ # The unsafe region is fully contained
+ return UNSAFE
+ else:
+ # We do not know if it is truly unsafe or not
+ safe = UNKNOWN
+ return safe
+
+ def cut_tube_till_unsafe(self, tube):
+ """
+ Truncated the tube before it intersect with unsafe set
+
+ Args:
+ tube (list): reach tube
+
+ Returns:
+ truncated tube
+ """
+ if not self._solver_dict:
+ return tube
+ # Cut the reach tube till it intersect with unsafe
+ for i in range(0, len(tube), 2):
+ lower = tube[i]
+ upper = tube[i + 1]
+ if self._check_intersection(lower, upper, 'Allmode'):
+ # we need to cut here
+ return tube[:i]
+
+ return tube
+
+ def _check_intersection(self, lower, upper, mode):
+ """
+ Check if current set intersect with the unsafe set
+
+ Args:
+ lower (list): lower bound of the current set
+ upper (list): upper bound of the current set
+ mode (str): the mode need to be checked
+
+ Returns:
+ Return a bool to indicate if the set intersect with the unsafe set
+ """
+ if mode in self._solver_dict:
+ cur_solver = self._solver_dict[mode][0]
+ symbols = self._solver_dict[mode][2]
+ elif 'Allmode' in self._solver_dict:
+ cur_solver = self._solver_dict['Allmode'][0]
+ symbols = self._solver_dict['Allmode'][2]
+ else:
+ raise ValueError("Unknown mode '" + mode + "'")
+
+ cur_solver.push()
+ for symbol in symbols:
+ cur_solver.add(self.varDic[symbol] >= lower[symbols[symbol]])
+ cur_solver.add(self.varDic[symbol] <= upper[symbols[symbol]])
+
+ check_result = cur_solver.check()
+
+ if check_result == sat:
+ cur_solver.pop()
+ return True
+ if check_result == unknown:
+ print("Z3 return unknown result")
+ exit() # TODO Proper return instead of exit
+ else:
+ cur_solver.pop()
+ return False
+
+ def _check_containment(self, lower, upper, mode):
+ """
+ Check if the current set is fully contained in unsafe region
+
+ Args:
+ lower (list): lower bound of the current set
+ upper (list): upper bound of the current set
+ mode (str): the mode need to be checked
+
+ Returns:
+ Return a bool to indicate if the set is fully contained in unsafe region
+ """
+ if mode in self._solver_dict:
+ cur_solver = self._solver_dict[mode][1]
+ symbols = self._solver_dict[mode][2]
+ elif 'Allmode' in self._solver_dict:
+ cur_solver = self._solver_dict['Allmode'][1]
+ symbols = self._solver_dict['Allmode'][2]
+ else:
+ raise ValueError("Unknown mode '" + mode + "'")
+
+ cur_solver.push()
+ for symbol in symbols:
+ cur_solver.add(self.varDic[symbol] >= lower[symbols[symbol]])
+ cur_solver.add(self.varDic[symbol] <= upper[symbols[symbol]])
+ check_result = cur_solver.check()
+
+ if check_result == sat:
+ cur_solver.pop()
+ return False
+ if check_result == unknown:
+ print("Z3 return unknown result")
+ exit() # TODO Proper return instead of exit
+ else:
+ cur_solver.pop()
+ return True
diff --git a/ourtool/dryvr/discrepancy/Global_Disc.py b/ourtool/dryvr/discrepancy/Global_Disc.py
new file mode 100644
index 0000000000000000000000000000000000000000..d84b5176e7c3ae8aca917af0ea533e5f6aab137d
--- /dev/null
+++ b/ourtool/dryvr/discrepancy/Global_Disc.py
@@ -0,0 +1,125 @@
+"""
+This file contains core bloating algorithm for dryvr
+"""
+
+from typing import List, Tuple
+import numpy as np
+import scipy as sp
+import scipy.spatial
+
+_TRUE_MIN_CONST = -10
+_EPSILON = 1.0e-6
+_SMALL_EPSILON = 1e-10
+
+
+def get_reachtube_segment(training_traces: np.ndarray, initial_radii: np.ndarray, method='PWGlobal') -> np.array:
+ num_traces: int = training_traces.shape[0]
+ ndims: int = training_traces.shape[2] # This includes time
+ trace_len: int = training_traces.shape[1]
+ center_trace: np.ndarray = training_traces[0, :, :]
+ trace_initial_time = center_trace[0, 0]
+ x_points: np.ndarray = center_trace[:, 0] - trace_initial_time
+ assert np.all(training_traces[0, :, 0] == training_traces[1:, :, 0])
+ y_points: np.ndarray = all_sensitivities_calc(training_traces, initial_radii)
+ points: np.ndarray = np.zeros((ndims - 1, trace_len, 2))
+ points[np.where(initial_radii != 0), 0, 1] = 1.0
+ points[:, :, 0] = np.reshape(x_points, (1, x_points.shape[0]))
+ points[:, 1:, 1] = y_points
+ normalizing_initial_set_radii: np.ndarray = initial_radii.copy()
+ normalizing_initial_set_radii[np.where(normalizing_initial_set_radii == 0)] = 1.0
+ df: np.ndarray = np.zeros((trace_len, ndims))
+ if method == 'PW':
+ df[:, 1:] = np.transpose(
+ points[:, :, 1] * np.reshape(normalizing_initial_set_radii, (normalizing_initial_set_radii.size, 1)))
+ elif method == 'PWGlobal':
+ # replace zeros with epsilons
+ # points[np.where(points[:, 0, 1] == 0), 0, 1] = 1.0e-100
+ # to fit exponentials make y axis log of sensitivity
+ points[:, :, 1] = np.maximum(points[:, :, 1], _EPSILON)
+ points[:, :, 1] = np.log(points[:, :, 1])
+ for dim_ind in range(1, ndims):
+ new_min = min(np.min(points[dim_ind - 1, 1:, 1]) + _TRUE_MIN_CONST, -10)
+ if initial_radii[dim_ind - 1] == 0:
+ # exclude initial set, then add true minimum points
+ new_points: np.ndarray = np.row_stack(
+ (np.array((points[dim_ind - 1, 1, 0], new_min)), np.array((points[dim_ind - 1, -1, 0], new_min))))
+ else:
+ # start from zero, then add true minimum points
+ new_points: np.ndarray = np.row_stack((points[dim_ind - 1, 0, :],
+ np.array((points[dim_ind - 1, 0, 0], new_min)),
+ np.array((points[dim_ind - 1, -1, 0], new_min))))
+ df[0, dim_ind] = initial_radii[dim_ind - 1]
+ # Tuple order is start_time, end_time, slope, y-intercept
+ cur_dim_points = np.concatenate((points[dim_ind - 1, 1:, :], new_points), axis=0)
+ cur_hull: sp.spatial.ConvexHull = sp.spatial.ConvexHull(cur_dim_points)
+ linear_separators: List[Tuple[float, float, float, float, int, int]] = []
+ vert_inds = list(zip(cur_hull.vertices[:-1], cur_hull.vertices[1:]))
+ vert_inds.append((cur_hull.vertices[-1], cur_hull.vertices[0]))
+ for end_ind, start_ind in vert_inds:
+ if cur_dim_points[start_ind, 1] != new_min and cur_dim_points[end_ind, 1] != new_min:
+ slope = (cur_dim_points[end_ind, 1] - cur_dim_points[start_ind, 1]) / (
+ cur_dim_points[end_ind, 0] - cur_dim_points[start_ind, 0])
+ y_intercept = cur_dim_points[start_ind, 1] - cur_dim_points[start_ind, 0] * slope
+ start_time = cur_dim_points[start_ind, 0]
+ end_time = cur_dim_points[end_ind, 0]
+ assert start_time < end_time
+ if start_time == 0:
+ linear_separators.append((start_time, end_time, slope, y_intercept, 0, end_ind + 1))
+ else:
+ linear_separators.append((start_time, end_time, slope, y_intercept, start_ind + 1, end_ind + 1))
+ linear_separators.sort()
+ prev_val = 0
+ prev_ind = 1 if initial_radii[dim_ind - 1] == 0 else 0
+ for linear_separator in linear_separators:
+ _, _, slope, y_intercept, start_ind, end_ind = linear_separator
+ assert prev_ind == start_ind
+ assert start_ind < end_ind
+ segment_t = center_trace[start_ind:end_ind + 1, 0]
+ segment_df = normalizing_initial_set_radii[dim_ind - 1] * np.exp(y_intercept) * np.exp(
+ slope * segment_t)
+ segment_df[0] = max(segment_df[0], prev_val)
+ df[start_ind:end_ind + 1, dim_ind] = segment_df
+ prev_val = segment_df[-1]
+ prev_ind = end_ind
+ else:
+ print('Discrepancy computation method,', method, ', is not supported!')
+ raise ValueError
+ assert (np.all(df >= 0))
+ reachtube_segment: np.ndarray = np.zeros((trace_len - 1, 2, ndims))
+ reachtube_segment[:, 0, :] = np.minimum(center_trace[1:, :] - df[1:, :], center_trace[:-1, :] - df[:-1, :])
+ reachtube_segment[:, 1, :] = np.maximum(center_trace[1:, :] + df[1:, :], center_trace[:-1, :] + df[:-1, :])
+ # assert 100% training accuracy (all trajectories are contained)
+ for trace_ind in range(training_traces.shape[0]):
+ if not (np.all(reachtube_segment[:, 0, :] <= training_traces[trace_ind, 1:, :]) and np.all(reachtube_segment[:, 1, :] >= training_traces[trace_ind, 1:, :])):
+ assert np.any(np.abs(training_traces[trace_ind, 0, 1:]-center_trace[0, 1:]) > initial_radii)
+ print(f"Warning: Trace #{trace_ind}", "of this initial set is sampled outside of the initial set because of floating point error and is not contained in the initial set")
+ return reachtube_segment
+
+
+def all_sensitivities_calc(training_traces: np.ndarray, initial_radii: np.ndarray):
+ num_traces: int
+ trace_len: int
+ ndims: int
+ num_traces, trace_len, ndims = training_traces.shape
+ normalizing_initial_set_radii: np.array = initial_radii.copy()
+ y_points: np.array = np.zeros((normalizing_initial_set_radii.shape[0], trace_len - 1))
+ normalizing_initial_set_radii[np.where(normalizing_initial_set_radii == 0)] = 1.0
+ for cur_dim_ind in range(1, ndims):
+ normalized_initial_points: np.array = training_traces[:, 0, 1:] / normalizing_initial_set_radii
+ initial_distances = sp.spatial.distance.pdist(normalized_initial_points, 'chebyshev') + _SMALL_EPSILON
+ for cur_time_ind in range(1, trace_len):
+ y_points[cur_dim_ind - 1, cur_time_ind - 1] = np.max((sp.spatial.distance.pdist(
+ np.reshape(training_traces[:, cur_time_ind, cur_dim_ind],
+ (training_traces.shape[0], 1)), 'chebychev'
+ )/ normalizing_initial_set_radii[cur_dim_ind - 1]) / initial_distances)
+ return y_points
+
+if __name__=="__main__":
+ with open("test.npy", "rb") as f:
+ training_traces = np.load(f)
+ initial_radii = np.array([1.96620653e-06, 2.99999995e+00, 3.07000514e-07, 8.84958773e-13, 1.05625786e-16, 3.72500000e+00, 0.00000000e+00, 0.00000000e+00])
+ result = get_reachtube_segment(training_traces, initial_radii, method='PWGlobal')
+ print(training_traces.dtype)
+ # plot_rtsegment_and_traces(result, training_traces[np.array((0, 6))])
+
+
diff --git a/ourtool/dryvr/discrepancy/PW_Discrepancy.py b/ourtool/dryvr/discrepancy/PW_Discrepancy.py
new file mode 100644
index 0000000000000000000000000000000000000000..53cac39418bfb66277fc7a37b25387648e5c076e
--- /dev/null
+++ b/ourtool/dryvr/discrepancy/PW_Discrepancy.py
@@ -0,0 +1,341 @@
+from __future__ import division, print_function
+
+import math
+
+
+def read_data(traces):
+ """ Read in all the traces """
+ error_thred_time = 1e-3
+
+ trace = traces[0]
+ delta_time = trace[1][0] - trace[0][0]
+
+ # Calculate variables
+ dimensions = len(trace[0])
+ dimensions_nt = dimensions - 1
+ end_time = trace[-1][0]
+
+ # Align all the traces
+ for i in range(len(traces)):
+ initial_time = traces[i][0][0]
+ for j in range(len(traces[i])):
+ traces[i][j][0] = traces[i][j][0] - initial_time
+
+ # reassign the start time and end time
+ start_time = 0
+ for i in range(len(traces)):
+ end_time = min(end_time, traces[i][-1][0])
+
+ # trim all the points after the end time
+ traces_trim = traces
+
+ # reassign trace_len
+ trace_len = len(traces_trim[0])
+
+ return (traces_trim, dimensions, dimensions_nt, trace_len, end_time, delta_time, start_time)
+
+
+def compute_diff(traces):
+ """ Compute difference between traces """
+ # Iterate over all combinations
+ traces_diff = []
+ for i in range(0, len(traces)):
+ for j in range(i + 1, len(traces)):
+
+ trace_diff = []
+
+ # Iterate over the data of the trace
+ for t in range(0, len(traces[i])):
+ diff = [abs(x_i - y_i) for x_i, y_i in zip(traces[i][t],
+ traces[j][t])]
+ trace_diff.append(diff[1:])
+
+ # Append to traces diff minus time difference
+ traces_diff.append(trace_diff)
+
+ # TODO Eliminate hardcoded file name
+ with open('output/tracediff2.txt', 'w') as write_file:
+ for i in range(len(traces_diff)):
+ for j in range(len(traces_diff[0])):
+ write_file.write(str(traces_diff[i][j]) + '\n')
+ write_file.write('\n')
+ return traces_diff
+
+
+def find_time_intervals(traces_diff, dimensions_nt, end_time, trace_len, delta_time, K_value):
+ """ Compute the time intervals """
+ # FIXME just do 1 dimension for now
+ # Iterate through all dimensions
+ num_ti = []
+ time_intervals = []
+
+ for i in range(0, dimensions_nt):
+
+ time_dim = []
+
+ # Obtain difference at start of interval
+ diff_0 = []
+ t_0 = 0.0
+ time_dim.append(t_0)
+ for k in range(0, len(traces_diff)):
+ diff_0.append(traces_diff[k][0][i])
+ # Iterate through all points in trace
+
+ for j in range(1, trace_len):
+ # Obtain difference at ith time of interval
+ diff_i = []
+ try:
+ for k in range(0, len(traces_diff)):
+ diff_i.append(traces_diff[k][j][i])
+ except IndexError:
+ print(trace_len)
+ print(k, j, i)
+ print(len(traces_diff[k]))
+ print(len(traces_diff[k][j]))
+
+ # Check time
+ t_i = j * delta_time
+ t = t_i - t_0
+ if t <= 0:
+ continue
+
+ # Compute ratios
+ ratio = []
+ for d_0, d_i in zip(diff_0, diff_i):
+ if d_i < 1E-3:
+ continue
+ elif d_0 < 1E-3:
+ continue
+
+ # NOTE not sure if this is right?
+ # ratio.append((1 / t) * math.log(d_i / d_0))
+ ratio.append(d_i / d_0)
+
+ # Check ratios if less than constant
+ # new_int = all(r <= 2.0*K_value[i] for r in ratio)
+ # new_int = all(r <= 2**(2*t)*K_value[i] for r in ratio)
+ new_int = all(r <= 1 for r in ratio)
+ if new_int == False:
+ if t_i != end_time:
+ time_dim.append(t_i)
+ diff_0 = diff_i
+ t_0 = t_i
+
+ # Append the time intervals
+ time_dim.append(end_time)
+ time_intervals.append(time_dim)
+ # record the number of time intervals
+ num_ti.append(len(time_intervals[i]) - 1)
+
+ return (time_intervals, num_ti)
+
+
+# Compute discrepancies
+def calculate_discrepancies(time_intervals, traces_diff, dimensions_nt, delta_time, K_value):
+ # FIXME
+ # Iterate over all dimensions
+ discrepancies = []
+ for nd in range(0, dimensions_nt):
+ # for nd in xrange(0, P_DIM):
+ disc = []
+
+ # Iterate over all time intervals
+ for ni in range(0, len(time_intervals[nd]) - 1):
+ t_0 = time_intervals[nd][ni]
+ t_e = time_intervals[nd][ni + 1]
+ # t_i = t_0 + delta_time
+
+ # FIXME (???)
+ # print "note",delta_time
+ points = int((t_e - t_0) / delta_time + 0.5) + 1
+ idx = int(t_0 / delta_time)
+
+ # try to find the best K and gamma
+ tmp_K_value = K_value[nd]
+ # Iterate over all trace difference
+ glpk_rows = []
+ close_flag = 0
+ for k in range(0, len(traces_diff)):
+
+ # Compute initial
+ diff_0 = traces_diff[k][0][nd]
+ if diff_0 <= 1E-3:
+ # print('Find two traces to be too closed!')
+ # print('use the default value!')
+ close_flag = 1
+ break
+ ln_0 = math.log(diff_0)
+
+ # FIXME need to reset the delta_time here
+ t_i = t_0 + delta_time
+ # print(disc)
+ # Obtain rows for GLPK
+ for r in range(1, points):
+ t_d = t_i - t_0
+ t_i += delta_time
+ diff_i = traces_diff[k][idx + r][nd]
+
+ if diff_i < 1E-3:
+ continue
+
+ ln_i = math.log(diff_i)
+
+ # compute the existing previous time interval discrepancy
+ discrepancy_now = 0
+ if len(disc) != 0:
+ for time_prev in range(0, len(disc)):
+ discrepancy_now = discrepancy_now + disc[time_prev] * (
+ time_intervals[nd][time_prev + 1] - time_intervals[nd][time_prev])
+
+ ln_d = ln_i - ln_0 - math.log(tmp_K_value) - discrepancy_now
+ glpk_rows.append([t_d, ln_d])
+
+ # Debugging algebraic solution
+ if close_flag == 0:
+ alg = [d / t for t, d in glpk_rows]
+ if len(alg) != 0:
+ alg_max = max(alg)
+ else:
+ alg_max = 0
+ else:
+ alg_max = 0
+
+ disc.append(alg_max)
+
+ # Append discrepancies
+ discrepancies.append(disc)
+
+ return discrepancies
+
+
+# Obtain bloated tube
+def generate_bloat_tube(traces, time_intervals, discrepancies, Initial_Delta, end_time, trace_len, dimensions_nt,
+ delta_time, K_value):
+
+ # Iterate over all dimensions
+ # FIXME
+ bloat_tube = []
+ for i in range(trace_len):
+ bloat_tube.append([])
+ bloat_tube.append([])
+
+ for nd in range(0, dimensions_nt):
+ # for nd in xrange(P_DIM - 1, P_DIM):
+
+ time_bloat = []
+ low_bloat = []
+ up_bloat = []
+
+ # To construct the reach tube
+ time_tube = []
+ tube = []
+
+ prev_delta = Initial_Delta[nd]
+
+ # Iterate over all intervals
+ previous_idx = -1
+
+ for ni in range(0, len(time_intervals[nd]) - 1):
+ t_0 = time_intervals[nd][ni]
+ t_e = time_intervals[nd][ni + 1]
+
+ if t_e == end_time:
+ points = int((t_e - t_0) / delta_time + 0.5) + 1
+ else:
+ points = int((t_e - t_0) / delta_time + 0.5)
+ idx = int(t_0 / delta_time)
+
+ gamma = discrepancies[nd][ni]
+
+ # Iterate over all points in center trace
+ for r in range(0, points):
+
+ current_idx = idx + r
+
+ if current_idx != previous_idx + 1:
+ # print('Index mismatch found!')
+ if current_idx == previous_idx:
+ idx += 1
+ elif current_idx == previous_idx + 2:
+ idx -= 1
+
+ pnt = traces[0][idx + r]
+ pnt_time = pnt[0]
+ pnt_data = pnt[nd + 1]
+
+ cur_delta = prev_delta * math.exp(gamma * delta_time)
+ max_delta = max(prev_delta, cur_delta)
+
+ time_bloat.append(pnt_time)
+ low_bloat.append(pnt_data - max_delta * K_value[nd])
+ up_bloat.append(pnt_data + max_delta * K_value[nd])
+
+ if nd == 0:
+ bloat_tube[2 * (idx + r)].append(pnt_time)
+ bloat_tube[2 * (idx + r)].append(pnt_data - max_delta * K_value[nd])
+ bloat_tube[2 * (idx + r) + 1].append(pnt_time + delta_time)
+ bloat_tube[2 * (idx + r) + 1].append(pnt_data + max_delta * K_value[nd])
+ else:
+ bloat_tube[2 * (idx + r)].append(pnt_data - max_delta * K_value[nd])
+ bloat_tube[2 * (idx + r) + 1].append(pnt_data + max_delta * K_value[nd])
+
+ prev_delta = cur_delta
+
+ previous_idx = idx + r
+
+ return bloat_tube
+
+# Print out the intervals and discrepancies
+def print_int_disc(discrepancies, time_intervals):
+ for nd in range(0, len(discrepancies)):
+ for p in range(0, len(discrepancies[nd])):
+ print('idx: ' + str(p) + ' int: ' + str(time_intervals[nd][p])
+ + ' to ' + str(time_intervals[nd][p + 1]) + ', disc: ' +
+ str(discrepancies[nd][p]))
+ print('')
+
+def PW_Bloat_to_tube(Initial_Delta, plot_flag, plot_dim, traces, K_value):
+ # Read data in
+ # if Mode == 'Const':
+ # K_value = [1.0,1.0,2.0]
+ # elif Mode == 'Brake':
+ # K_value = [1.0,1.0,7.0]
+
+ # if Mode == 'Const;Const':
+ # K_value = [1.0,1.0,2.0,1.0,1.0,2.0]
+ # elif Mode == 'Brake;Const':
+ # K_value = [1.0,1.0,2.0,1.0,1.0,2.0]
+
+ # elif Mode == 'Brake;Brake':
+ # K_value = [1.0,1.0,5.0,1.0,1.0,2.0]
+
+ traces, dimensions, dimensions_nt, trace_len, end_time, delta_time, start_time = read_data(traces)
+ # Compute difference between traces
+ traces_diff = compute_diff(traces)
+ # print traces_diff
+
+ # Find time intervals for discrepancy calculations
+ time_intervals, num_ti = find_time_intervals(traces_diff, dimensions_nt, end_time, trace_len, delta_time, K_value)
+ # print('number of time intervals:')
+ # print num_ti
+ # Discrepancy calculation
+ discrepancies = calculate_discrepancies(time_intervals, traces_diff, dimensions_nt, delta_time, K_value)
+ # print('The K values')
+ # print K_value,
+ # system.exit('test')
+
+ # Write discrepancies to file
+ # write_to_file(time_intervals,discrepancies,write_path +' disp.txt', 'disc')
+
+ # Nicely print the intervals and discrepancies
+ # print_int_disc(discrepancies,time_intervals)
+
+ # Bloat the tube using time intervals
+ reach_tube = generate_bloat_tube(traces, time_intervals, discrepancies, Initial_Delta, end_time, trace_len,
+ dimensions_nt, delta_time, K_value)
+
+ # if plot_flag:
+ # plot_traces(traces, plot_dim, reach_tube)
+ # plt.show()
+
+ return reach_tube
diff --git a/ourtool/dryvr/discrepancy/__init__.py b/ourtool/dryvr/discrepancy/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ourtool/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb b/ourtool/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2fd64429bf421126b7000c94ce0f6fd186fbd01f
--- /dev/null
+++ b/ourtool/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/ourtool/dryvr/plotter/README.md b/ourtool/dryvr/plotter/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..2ea9534697e8de9a5fa0ab9831ce4102499684ab
--- /dev/null
+++ b/ourtool/dryvr/plotter/README.md
@@ -0,0 +1 @@
+This folder consist plotter code for DryVR reachtube output
diff --git a/ourtool/dryvr/plotter/__init__.py b/ourtool/dryvr/plotter/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..5dbc287d927065e02e00ed2cccd0389194764818
--- /dev/null
+++ b/ourtool/dryvr/plotter/__init__.py
@@ -0,0 +1,22 @@
+# _oo0oo_
+# o8888888o
+# 88" . "88
+# (| -_- |)
+# 0\ = /0
+# ___/`---'\___
+# .' \\| |// '.
+# / \\||| : |||// \
+# / _||||| -:- |||||- \
+# | | \\\ - /// | |
+# | \_| ''\---/'' |_/ |
+# \ .-\__ '-' ___/-. /
+# ___'. .' /--.--\ `. .'___
+# ."" '< `.___\_<|>_/___.' >' "".
+# | | : `- \`.;`\ _ /`;.`/ - ` : | |
+# \ \ `_. \_ __\ /__ _/ .-` / /
+# =====`-.____`.___ \_____/___.-`___.-'=====
+# `=---='
+#
+#
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Codes are far away from bugs with the protection
\ No newline at end of file
diff --git a/ourtool/dryvr/plotter/linkednode.py b/ourtool/dryvr/plotter/linkednode.py
new file mode 100644
index 0000000000000000000000000000000000000000..21ee9f462f066b926a39cea050e034db2ab3c078
--- /dev/null
+++ b/ourtool/dryvr/plotter/linkednode.py
@@ -0,0 +1,18 @@
+"""
+This is a data structure to hold reachtube data per node
+"""
+
+
+class LinkedNode:
+ def __init__(self, node_id, file_name):
+ self.file_name = file_name.strip()
+ self.nodeId = node_id
+ self.lower_bound = {}
+ self.upper_bound = {}
+ self.child = {}
+
+ def print_tube(self):
+ for key in sorted(self.lower_bound):
+ if key in self.upper_bound:
+ print(self.upper_bound[key])
+ print(self.lower_bound[key])
diff --git a/ourtool/dryvr/plotter/parser.py b/ourtool/dryvr/plotter/parser.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d84ecd8a55fac0e7ab966cbfb3774d1fd30b4bc
--- /dev/null
+++ b/ourtool/dryvr/plotter/parser.py
@@ -0,0 +1,84 @@
+"""
+This file consist parser code for DryVR reachtube output
+"""
+
+from src.plotter.linkednode import LinkedNode
+
+
+def parse(data):
+ init_node = None
+ prev_node = None
+ cur_node = None
+ lower_bound = {}
+ upper_bound = {}
+ y_min = [float("inf") for _ in range(len(data[2].strip().split()))]
+ y_max = [float("-inf") for _ in range(len(data[2].strip().split()))]
+
+ for line in data:
+ # This is a mode indicator
+ if ',' in line or '->' in line or line.strip().isalpha() or len(line.strip()) == 1:
+ insert_data(cur_node, lower_bound, upper_bound)
+ # There is new a transition
+ if '->' in line:
+ mode_list = line.strip().split('->')
+ prev_node = init_node
+ for i in range(1, len(mode_list) - 1):
+ prev_node = prev_node.child[mode_list[i]]
+ cur_node = prev_node.child.setdefault(mode_list[-1], LinkedNode(mode_list[-1], line))
+ else:
+ cur_node = LinkedNode(line.strip(), line)
+ if not init_node:
+ init_node = cur_node
+ else:
+ cur_node = init_node
+ # Using dictionary because we want to concat data
+ lower_bound = {}
+ upper_bound = {}
+ LOWER = True
+
+ else:
+ line = list(map(float, line.strip().split()))
+ if len(line) <= 1:
+ continue
+ if LOWER:
+ LOWER = False
+ # This data appeared in lower_bound before, concat the data
+ if line[0] in lower_bound:
+ for i in range(1, len(line)):
+ lower_bound[line[0]][i] = min(lower_bound[line[0]][i], line[i])
+ else:
+ lower_bound[line[0]] = line
+
+ for i in range(len(line)):
+ y_min[i] = min(y_min[i], line[i])
+ else:
+ LOWER = True
+ if line[0] in upper_bound:
+ for i in range(1, len(line)):
+ upper_bound[line[0]][i] = max(upper_bound[line[0]][i], line[i])
+ else:
+ upper_bound[line[0]] = line
+
+ for i in range(len(line)):
+ y_max[i] = max(y_max[i], line[i])
+ insert_data(cur_node, lower_bound, upper_bound)
+ return init_node, y_min, y_max
+
+
+def insert_data(node, lower_bound, upper_bound):
+ if not node or len(lower_bound) == 0:
+ return
+
+ for key in lower_bound:
+ if key in node.lower_bound:
+ for i in range(1, len(lower_bound[key])):
+ node.lower_bound[key][i] = min(node.lower_bound[key][i], lower_bound[key][i])
+ else:
+ node.lower_bound[key] = lower_bound[key]
+
+ for key in sorted(upper_bound):
+ if key in node.upper_bound:
+ for i in range(1, len(upper_bound[key])):
+ node.upper_bound[key][i] = max(node.upper_bound[key][i], upper_bound[key][i])
+ else:
+ node.upper_bound[key] = upper_bound[key]
diff --git a/ourtool/dryvr/plotter/plot.py b/ourtool/dryvr/plotter/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..b07e575197d200745f9a83a5fbaa866bc6229fb1
--- /dev/null
+++ b/ourtool/dryvr/plotter/plot.py
@@ -0,0 +1,78 @@
+"""
+This file consist main plotter code for DryVR reachtube output
+"""
+
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+
+colors = ['red', 'green', 'blue', 'yellow', 'black']
+
+
+def plot(node, dim, y_min, y_max, x_dim):
+ fig1 = plt.figure()
+ ax1 = fig1.add_subplot('111')
+ lower_bound = []
+ upper_bound = []
+ for key in sorted(node.lower_bound):
+ lower_bound.append(node.lower_bound[key])
+ for key in sorted(node.upper_bound):
+ upper_bound.append(node.upper_bound[key])
+
+ for i in range(min(len(lower_bound), len(upper_bound))):
+ lb = list(map(float, lower_bound[i]))
+ ub = list(map(float, upper_bound[i]))
+
+ for ci, d in enumerate(dim):
+ rect = patches.Rectangle((lb[x_dim], lb[d]), ub[x_dim] - lb[x_dim], ub[d] - lb[d],
+ color=colors[ci % len(colors)], alpha=0.7)
+ ax1.add_patch(rect)
+
+ y_axis_min = min([y_min[i] for i in dim])
+ y_axis_max = max([y_max[i] for i in dim])
+ ax1.set_title(node.nodeId, fontsize=12)
+ ax1.set_ylim(bottom=y_axis_min, top=y_axis_max)
+ ax1.plot()
+ # TODO Eliminate hardcoded folders
+ fig1.savefig('output/' + node.file_name + '.png', format='png', dpi=200)
+
+
+def rrt_plot(lower_bound, upper_bound, x_dim, y_dim, goal, unsafe_list, region, initial):
+ fig1 = plt.figure()
+ ax1 = fig1.add_subplot('111')
+ x_min = region[0][0]
+ y_min = region[0][1]
+ x_max = region[1][0]
+ y_max = region[1][1]
+
+ # Draw the path
+ for i in range(min(len(lower_bound), len(upper_bound))):
+ lb = list(map(float, lower_bound[i]))
+ ub = list(map(float, upper_bound[i]))
+
+ rect = patches.Rectangle((lb[x_dim], lb[y_dim]), ub[x_dim] - lb[x_dim], ub[y_dim] - lb[y_dim], color='blue',
+ alpha=0.7)
+ ax1.add_patch(rect)
+
+ # Draw the goal
+ if goal:
+ lb, ub = goal
+ rect = patches.Rectangle((lb[0], lb[1]), ub[0] - lb[0], ub[1] - lb[1], color='green', alpha=0.7)
+ ax1.add_patch(rect)
+
+ if initial:
+ lb, ub = initial
+ rect = patches.Rectangle((lb[0], lb[1]), ub[0] - lb[0], ub[1] - lb[1], color='yellow', alpha=0.7)
+ ax1.add_patch(rect)
+
+ # Draw the unsafe
+ if unsafe_list:
+ for unsafe in unsafe_list:
+ lb, ub = unsafe
+ rect = patches.Rectangle((lb[0], lb[1]), ub[0] - lb[0], ub[1] - lb[1], color='red', alpha=0.7)
+ ax1.add_patch(rect)
+
+ ax1.set_title("RRT", fontsize=12)
+ ax1.set_ylim(bottom=y_min, top=y_max)
+ ax1.set_xlim(left=x_min, right=x_max)
+ ax1.plot()
+ fig1.savefig('output/rrt.png', format='png', dpi=200)
diff --git a/ourtool/dryvr/plotter/plot_util.py b/ourtool/dryvr/plotter/plot_util.py
new file mode 100644
index 0000000000000000000000000000000000000000..6718c923b4dd9148d8f6b55464a662daeab2a30f
--- /dev/null
+++ b/ourtool/dryvr/plotter/plot_util.py
@@ -0,0 +1,46 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.patches import Rectangle
+
+
+def plot_rtsegment_and_traces(rtsegment: np.ndarray, traces: np.ndarray):
+ for dim_ind in range(1, traces.shape[2]):
+ fig, ax = plt.subplots(1)
+ facecolor = 'r'
+ for trace_ind in range(traces.shape[0]):
+ ax.plot(traces[trace_ind, :, 0], traces[trace_ind, :, dim_ind])
+ for hrect_ind in range(rtsegment.shape[0]):
+ ax.add_patch(Rectangle((rtsegment[hrect_ind, 0, 0], rtsegment[hrect_ind, 0, dim_ind]), rtsegment[hrect_ind, 1, 0]-rtsegment[hrect_ind, 0, 0],
+ rtsegment[hrect_ind, 1, dim_ind] - rtsegment[hrect_ind, 0, dim_ind], alpha=0.1, facecolor='r'))
+ ax.set_title(f'dim #{dim_ind}')
+ fig.canvas.draw()
+ plt.show()
+
+
+def plot_traces(traces, dim, bloat_tube):
+ """ Plot the traces """
+ # Iterate over all individual traces
+ for i in range(0, len(traces)):
+ trace = traces[i]
+
+ # Obtain desired dimension
+ time = []
+ data = []
+ for j in range(0, len(trace)):
+ # for j in xrange(0,2):
+ time.append(trace[j][0])
+ data.append(trace[j][dim])
+
+ # Plot data
+ if i == 0:
+ plt.plot(time, data, 'b')
+ else:
+ plt.plot(time, data, 'r')
+
+ time = [row[0] for row in bloat_tube]
+ value = [row[dim] for row in bloat_tube]
+ time_bloat = [time[i] for i in range(0, len(value), 2)]
+ lower_bound = [value[i] for i in range(0, len(value), 2)]
+ upper_bound = [value[i + 1] for i in range(0, len(value), 2)]
+ plt.plot(time_bloat, lower_bound, 'g')
+ plt.plot(time_bloat, upper_bound, 'g')
\ No newline at end of file
diff --git a/ourtool/dryvr/plotter/plotter_interface.ipynb b/ourtool/dryvr/plotter/plotter_interface.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..de6c0d1e2aa8693652246f3ea36fa58d135ed351
--- /dev/null
+++ b/ourtool/dryvr/plotter/plotter_interface.ipynb
@@ -0,0 +1,132 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window.Plotly) {{require(['plotly'],function(plotly) {window.Plotly=plotly;});}}</script>"
+ ],
+ "text/vnd.plotly.v1+html": [
+ "<script>requirejs.config({paths: { 'plotly': ['https://cdn.plot.ly/plotly-latest.min']},});if(!window.Plotly) {{require(['plotly'],function(plotly) {window.Plotly=plotly;});}}</script>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "data": [
+ {
+ "type": "scatter",
+ "x": [
+ 1,
+ 2,
+ 3,
+ 4
+ ],
+ "y": [
+ 10,
+ 15,
+ 13,
+ 17
+ ]
+ },
+ {
+ "type": "scatter",
+ "x": [
+ 1,
+ 2,
+ 3,
+ 4
+ ],
+ "y": [
+ 16,
+ 5,
+ 11,
+ 9
+ ]
+ }
+ ],
+ "layout": {}
+ },
+ "text/html": [
+ "<div id=\"2839cfe2-5e9b-4166-bdb3-f660cb538d3e\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"2839cfe2-5e9b-4166-bdb3-f660cb538d3e\", [{\"y\": [10, 15, 13, 17], \"x\": [1, 2, 3, 4], \"type\": \"scatter\"}, {\"y\": [16, 5, 11, 9], \"x\": [1, 2, 3, 4], \"type\": \"scatter\"}], {}, {\"linkText\": \"Export to plot.ly\", \"showLink\": true})});</script>"
+ ],
+ "text/vnd.plotly.v1+html": [
+ "<div id=\"2839cfe2-5e9b-4166-bdb3-f660cb538d3e\" style=\"height: 525px; width: 100%;\" class=\"plotly-graph-div\"></div><script type=\"text/javascript\">require([\"plotly\"], function(Plotly) { window.PLOTLYENV=window.PLOTLYENV || {};window.PLOTLYENV.BASE_URL=\"https://plot.ly\";Plotly.newPlot(\"2839cfe2-5e9b-4166-bdb3-f660cb538d3e\", [{\"y\": [10, 15, 13, 17], \"x\": [1, 2, 3, 4], \"type\": \"scatter\"}, {\"y\": [16, 5, 11, 9], \"x\": [1, 2, 3, 4], \"type\": \"scatter\"}], {}, {\"linkText\": \"Export to plot.ly\", \"showLink\": true})});</script>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "from plotly.offline import init_notebook_mode, iplot\n",
+ "from plotly.graph_objects import Data, Scatter\n",
+ "\n",
+ "init_notebook_mode(connected=True) # initiate notebook for offline plot\n",
+ "\n",
+ "trace0 = Scatter(\n",
+ " x=[1, 2, 3, 4],\n",
+ " y=[10, 15, 13, 17]\n",
+ ")\n",
+ "trace1 = Scatter(\n",
+ " x=[1, 2, 3, 4],\n",
+ " y=[16, 5, 11, 9]\n",
+ ")\n",
+ "data = Data([trace0, trace1])\n",
+ "\n",
+ "iplot(data) # use plotly.offline.iplot for offline plot"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import networkx as nx\n",
+ "\n",
+ "G=nx.random_geometric_graph(200,0.125)\n",
+ "pos=nx.get_node_attributes(G,'pos')\n",
+ "\n",
+ "dmin=1\n",
+ "ncenter=0\n",
+ "for n in pos:\n",
+ " x,y=pos[n]\n",
+ " d=(x-0.5)**2+(y-0.5)**2\n",
+ " if d<dmin:\n",
+ " ncenter=n\n",
+ " dmin=d\n",
+ "\n",
+ "p=nx.single_source_shortest_path_length(G,ncenter)"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 2",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.10"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
\ No newline at end of file
diff --git a/ourtool/map/lane_map.py b/ourtool/map/lane_map.py
index c62e2f48993e6e478dd6cb888b1879e15ecf8c8f..31cd4b7ab4b0d3e924d7115412eb71b0e681c53d 100644
--- a/ourtool/map/lane_map.py
+++ b/ourtool/map/lane_map.py
@@ -1,5 +1,6 @@
from typing import Dict, List
import copy
+from enum import Enum,auto
from ourtool.map.lane_segment import LaneSegment
@@ -20,6 +21,8 @@ class LaneMap:
self.right_lane_dict[lane_seg.id] = []
def has_left(self, lane_idx):
+ if isinstance(lane_idx, Enum):
+ lane_idx = lane_idx.name
if lane_idx not in self.lane_segment_dict:
Warning(f'lane {lane_idx} not available')
return False
@@ -28,12 +31,16 @@ class LaneMap:
def left_lane(self, lane_idx):
assert all((elem in self.left_lane_dict) for elem in self.lane_segment_dict)
+ if isinstance(lane_idx, Enum):
+ lane_idx = lane_idx.name
if lane_idx not in self.left_lane_dict:
raise ValueError(f"lane_idx {lane_idx} not in lane_segment_dict")
left_lane_list = self.left_lane_dict[lane_idx]
return copy.deepcopy(left_lane_list)
def has_right(self, lane_idx):
+ if isinstance(lane_idx, Enum):
+ lane_idx = lane_idx.name
if lane_idx not in self.lane_segment_dict:
Warning(f'lane {lane_idx} not available')
return False
@@ -42,10 +49,23 @@ class LaneMap:
def right_lane(self, lane_idx):
assert all((elem in self.right_lane_dict) for elem in self.lane_segment_dict)
+ if isinstance(lane_idx, Enum):
+ lane_idx = lane_idx.name
if lane_idx not in self.right_lane_dict:
raise ValueError(f"lane_idx {lane_idx} not in lane_segment_dict")
right_lane_list = self.right_lane_dict[lane_idx]
return copy.deepcopy(right_lane_list)
def lane_geometry(self, lane_idx):
- return self.lane_segment_dict[lane_idx].get_geometry()
\ No newline at end of file
+ if isinstance(lane_idx, Enum):
+ lane_idx = lane_idx.name
+ return self.lane_segment_dict[lane_idx].get_geometry()
+
+ def get_longitudinal_error(self, lane_idx, agent_state):
+ raise NotImplementedError
+
+ def get_lateral_error(self, lane_idx, agent_state):
+ raise NotImplementedError
+
+ def get_altitude_error(self, lane_idx, agent_state):
+ raise NotImplementedError
\ No newline at end of file
diff --git a/ourtool/scenario/scenario.py b/ourtool/scenario/scenario.py
index e5cbe70efc5372021b55d08b4f22733fa5e16e0a..fac09725138356405c73021ecdda8587050849ac 100644
--- a/ourtool/scenario/scenario.py
+++ b/ourtool/scenario/scenario.py
@@ -3,11 +3,15 @@ import copy
import itertools
import ast
+import numpy as np
+from rsa import verify
+
from ourtool.agents.base_agent import BaseAgent
from ourtool.automaton.guard import GuardExpressionAst
from pythonparser import Guard
from pythonparser import Reset
-from ourtool.simulator.simulator import Simulator
+from ourtool.analysis.simulator import Simulator
+from ourtool.analysis.verifier import Verifier
from ourtool.map.lane_map import LaneMap
class FakeSensor:
@@ -56,6 +60,7 @@ class Scenario:
def __init__(self):
self.agent_dict = {}
self.simulator = Simulator()
+ self.verifier = Verifier()
self.init_dict = {}
self.init_mode_dict = {}
self.map = None
@@ -84,6 +89,20 @@ class Scenario:
agent_list.append(self.agent_dict[agent_id])
return self.simulator.simulate(init_list, init_mode_list, agent_list, self, time_horizon, self.map)
+ def verify(self, time_horizon):
+ init_list = []
+ init_mode_list = []
+ agent_list = []
+ for agent_id in self.agent_dict:
+ init = self.init_dict[agent_id]
+ tmp = np.array(init)
+ if tmp.ndim < 2:
+ init = [init, init]
+ init_list.append(init)
+ init_mode_list.append(self.init_mode_dict[agent_id])
+ agent_list.append(self.agent_dict[agent_id])
+ return self.verifier.compute_full_reachtube(init_list, init_mode_list, agent_list, self, time_horizon, self.map)
+
def get_all_transition(self, state_dict):
lane_map = self.map
satisfied_guard = []
diff --git a/user/simple_map.py b/user/simple_map.py
index 3c777ad471c8ce7e326d3c1f7f6f2cc6d2ea60fd..848171f5e7748577bfefcd39e5cd9a0d294719c7 100644
--- a/user/simple_map.py
+++ b/user/simple_map.py
@@ -22,6 +22,27 @@ class SimpleMap2(LaneMap):
self.right_lane_dict[segment1.id].append(segment2.id)
self.right_lane_dict[segment2.id].append(segment3.id)
+class SimpleMap3(LaneMap):
+ def __init__(self):
+ super().__init__()
+ segment1 = LaneSegment('Lane0', ((0,-30),33))
+ segment2 = LaneSegment('Lane1', ((0,-30),30))
+ segment3 = LaneSegment('Lane3', ((0,-30),27))
+ self.add_lanes([segment1,segment2,segment3])
+ self.left_lane_dict[segment2.id].append(segment1.id)
+ self.left_lane_dict[segment3.id].append(segment2.id)
+ self.right_lane_dict[segment1.id].append(segment2.id)
+ self.right_lane_dict[segment2.id].append(segment3.id)
+
+ def get_longitudinal_error(self, lane_idx, agent_state):
+ return super().get_longitudinal_error(lane_idx, agent_state)
+
+ def get_lateral_error(self, lane_idx, agent_state):
+ return super().get_lateral_error(lane_idx, agent_state)
+
+ def get_altitude_error(self, lane_idx, agent_state):
+ return super().get_altitude_error(lane_idx, agent_state)
+
if __name__ == "__main__":
test_map = SimpleMap()
print(test_map.left_lane_dict)