diff --git a/.gitignore b/.gitignore
index 6f9cf12591ea2083c6185bf0713d653f21547341..28affb7420b24d8cb24638534a99a43d89b56720 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,4 @@
 __pycache__/
-**/__pycache__/
\ No newline at end of file
+**/__pycache__/
+.vscode/
+.idea/
\ No newline at end of file
diff --git a/README.md b/README.md
index 915a12648b8a7c10c9ebed8c507ab8344fe26fe5..411bc03f7d93fa8011c33f31912bb4d029708edd 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,43 @@
-# Verification/Simulation Tool
+# GraphGeneration
+## Installation
+The package requires python 3.8+. All the required packages can be installed through
 
-ourtool folder contains the verification/simulation work
-pythonparser has the python parsing classes (should be moved)
-Example files are found in first level directory such as example_car_lane_switch.py. These types of files are what you can run to use the tool
+```
+python3 -m pip install -r requirements.txt
+```
+
+## Examples
+The package comes with two controller examples
+- The first example consists a scenario with a single vehicle, which can perform lane switch based on its location. The 
+first example can be run by using command
+
+```
+python3 example_car_lane_switch.py
+```
+
+- The second example consists a scenario with two vehicles, which can perform lane switch based on their relative position.
+The second example can be run using command
+
+```
+python3 example_two_car_lane_switch.py
+```
+
+## Package Structure
+
+The source code of the package is contained in the src folder, which contains the following sub-directories.
+
+- **scene_verifier**, which contains building blocks for creating and analyzing scenarios.
+  
+  - **scene_verifier/scenario** contains code for the scenario base class. A scenario is constructed by several **agents** with continuous dynamics and controller, a **map** and a **sensor** defining how different agents interact with each other.
+  - **scene_verifier/agents** contains code for the agent base class in the scenario. 
+  - **scene_verifier/map** contains code for the lane map base class and corresponding utilities in the scenario.
+  - **scene_verifier/code_parser** contains code for converting the controller code to ASTs. 
+  - **scene_verifier/automaton** contains code implementing components in hybrid-automaton
+  - **scene_verifier/analysis** contains the **Simulator** and **Verifier** and related utilities for doing analysis of the scenario
+  - **scene_verifier/dryvr** dryvr for computing reachable sets
+
+
+- **example** contains example map, sensor and agents that we provided
+
+
+- **plotter** contains code for visualizing the computed results
diff --git a/billiard_input.json b/billiard_input.json
deleted file mode 100644
index ce0b5c076e89d76ee15184c75a36937f6c748d3b..0000000000000000000000000000000000000000
--- a/billiard_input.json
+++ /dev/null
@@ -1,8 +0,0 @@
-{
-	"initialVertex":0,
-	"initialSet":[[1.5, 0.0, 1.0, 1.0],[1.51, 0.0, 1.0, 1.0]],
-	"unsafeSet":"@Allmode:posy>=12.0",
-	"timeHorizon":16,
-	"directory":"examples/billiard",
-	"deterministic":true
-}
diff --git a/billiard_out.json b/billiard_out.json
deleted file mode 100644
index 7ec65a276cc452a9a8fa220e80465a8d47c795af..0000000000000000000000000000000000000000
--- a/billiard_out.json
+++ /dev/null
@@ -1,135 +0,0 @@
-{
-    "initialVertex": 0,
-    "initialSet": [
-        [
-            1.5,
-            0.0,
-            1.0,
-            1.0
-        ],
-        [
-            1.51,
-            0.0,
-            1.0,
-            1.0
-        ]
-    ],
-    "unsafeSet": "@Allmode:posy>=12.0",
-    "timeHorizon": 16,
-    "directory": "examples/billiard",
-    "deterministic": true,
-    "variables": [
-        "posx",
-        "posy",
-        "vx",
-        "vy"
-    ],
-    "vertex": [
-        "NormalA",
-        "NormalB",
-        "NormalC",
-        "NormalD"
-    ],
-    "edge": [
-        [
-            0,
-            0
-        ],
-        [
-            0,
-            1
-        ],
-        [
-            0,
-            2
-        ],
-        [
-            0,
-            3
-        ],
-        [
-            1,
-            1
-        ],
-        [
-            1,
-            0
-        ],
-        [
-            1,
-            2
-        ],
-        [
-            1,
-            3
-        ],
-        [
-            2,
-            2
-        ],
-        [
-            2,
-            0
-        ],
-        [
-            2,
-            1
-        ],
-        [
-            2,
-            3
-        ],
-        [
-            3,
-            3
-        ],
-        [
-            3,
-            0
-        ],
-        [
-            3,
-            1
-        ],
-        [
-            3,
-            2
-        ]
-    ],
-    "guards": [
-        "And(posy<0,posy>=-0.01)",
-        "And(posx<0,posx>=-0.01)",
-        "And(posx<=5.01,posx>5)",
-        "And(posy>5,posy<=5.01)",
-        "And(posy<0,posy>=-0.01)",
-        "And(posx<0,posx>=-0.01)",
-        "And(posx<=5.01,posx>5)",
-        "And(posy>5,posy<=5.01)",
-        "And(posy<0,posy>=-0.01)",
-        "And(posx<0,posx>=-0.01)",
-        "And(posx<=5.01,posx>5)",
-        "And(posy>5,posy<=5.01)",
-        "And(posy<0,posy>=-0.01)",
-        "And(posx<0,posx>=-0.01)",
-        "And(posx<=5.01,posx>5)",
-        "And(posy>5,posy<=5.01)"
-    ],
-    "resets": [
-        "vy=-vy;posy=0",
-        "vx=-vx;posx=0",
-        "vx=-vx;posx=5",
-        "vy=-vy;posy=5",
-        "vy=-vy;posy=0",
-        "vx=-vx;posx=0",
-        "vx=-vx;posx=5",
-        "vy=-vy;posy=5",
-        "vy=-vy;posy=0",
-        "vx=-vx;posx=0",
-        "vx=-vx;posx=5",
-        "vy=-vy;posy=5",
-        "vy=-vy;posy=0",
-        "vx=-vx;posx=0",
-        "vx=-vx;posx=5",
-        "vy=-vy;posy=5"
-    ]
-}
\ No newline at end of file
diff --git a/billiardsrc.c b/billiardsrc.c
deleted file mode 100644
index 5e56e4814156785b5c85ec88f248c3251a89bd8d..0000000000000000000000000000000000000000
--- a/billiardsrc.c
+++ /dev/null
@@ -1,119 +0,0 @@
-enum modes {NormalA, NormalB, NormalC, NormalD};
-
-struct State {
-int posx;
-int posy;
-int vx;
-int vy;
-enum modes mode;
-};
-
-	
-	
-struct State P(struct State s) {
-	int posx = s.posx;
-	int posy = s.posy;
-	int vx = s.vx;
-	int vy = s.vy;
-	enum modes state = s.mode;
-	if (state ==NormalA) {
-		if (posy<0 && posy>=-0.01) {
-			vy=-vy;
-			posy=0;
-			state=NormalA;
-		}
-		if (posx<0 && posx>=-0.01) {
-			vx=-vx;
-			posx=0;
-			state=NormalB;
-		}
-		if (posx<=5.01 && posx>5) {
-			vx=-vx;
-			posx=5;
-			state=NormalC;
-		}
-		if (posy>5 && posy<=5.01) {
-			vy=-vy;
-			posy=5;
-			state = NormalD;
-		}
-	
-	}
-	if (state ==NormalB) {
-		if (posy<0 && posy>=-0.01) {
-			vy=-vy;
-			posy=0;
-			state=NormalB;
-		}
-		if (posx<0 && posx>=-0.01) {
-			vx=-vx;
-			posx=0;
-			state=NormalA;
-		}
-		if (posx<=5.01 && posx>5) {
-			vx=-vx;
-			posx=5;
-			state=NormalC;
-		}
-		if (posy>5 && posy<=5.01) {
-			vy=-vy;
-			posy=5;
-			state = NormalD;
-		}
-	
-	}
-	if (state ==NormalC) {
-		if (posy<0 && posy>=-0.01) {
-			vy=-vy;
-			posy=0;
-			state=NormalC;
-		}
-		if (posx<0 && posx>=-0.01) {
-			vx=-vx;
-			posx=0;
-			state=NormalA;
-		}
-		if (posx<=5.01 && posx>5) {
-			vx=-vx;
-			posx=5;
-			state=NormalB;
-		}
-		if (posy>5 && posy<=5.01) {
-			vy=-vy;
-			posy=5;
-			state = NormalD;
-		}
-	
-	}
-	
-	if (state ==NormalD) {
-		if (posy<0 && posy>=-0.01) {
-			vy=-vy;
-			posy=0;
-			state=NormalD;
-		}
-		if (posx<0 && posx>=-0.01) {
-			vx=-vx;
-			posx=0;
-			state=NormalA;
-		}
-		if (posx<=5.01 && posx>5) {
-			vx=-vx;
-			posx=5;
-			state=NormalB;
-		}
-		if (posy>5 && posy<=5.01) {
-			vy=-vy;
-			posy=5;
-			state = NormalC;
-		}
-	
-	}
-	s.posx = posx;
-	s.posy= posy;
-	s.vx = vx;
-	s.vy = vy;
-	s.mode = state;
-	return s;
-
-}
diff --git a/cartoy-out.json b/cartoy-out.json
deleted file mode 100644
index d12a408bf00d1e227a818344a87bee93bd3e6e09..0000000000000000000000000000000000000000
--- a/cartoy-out.json
+++ /dev/null
@@ -1,89 +0,0 @@
-{
-    "unsafeSet": "",
-    "initialSet": [
-        [
-            0,
-            -0.2,
-            0
-        ],
-        [
-            0,
-            0.2,
-            0
-        ]
-    ],
-    "timeHorizon": 10,
-    "directory": "examples/curve_controller_hybrid",
-    "initialVertex": 0,
-    "deterministic": true,
-    "variables": [
-        "vx",
-        "vy",
-        "dist"
-    ],
-    "vertex": [
-        "0",
-        "1",
-        "2"
-    ],
-    "edge": [
-        [
-            "0",
-            "0"
-        ],
-        [
-            "1",
-            "0"
-        ],
-        [
-            "2",
-            "0"
-        ],
-        [
-            "0",
-            "1"
-        ],
-        [
-            "1",
-            "1"
-        ],
-        [
-            "2",
-            "1"
-        ],
-        [
-            "0",
-            "2"
-        ],
-        [
-            "1",
-            "2"
-        ],
-        [
-            "2",
-            "2"
-        ]
-    ],
-    "guards": [
-        "And(dist < 10,(And(dist < 2,(state != move_over))))",
-        "And(dist < 10,(And(dist < 2,(state != move_over))))",
-        "And(dist < 10,(And(dist < 2,(state != move_over))))",
-        "And(dist < 10,(And(dist < 2,(dist > 2))))",
-        "And(dist < 10,(And(dist < 2,(dist > 2))))",
-        "And(dist < 10,(And(dist < 2,(dist > 2))))",
-        "And(dist < 10,(state == move_over))",
-        "And(dist < 10,(state == move_over))",
-        "And(dist < 10,(state == move_over))"
-    ],
-    "resets": [
-        "vy=vy -5;vx= vx + 5",
-        "vy=vy -5;vx= vx + 5",
-        "vy=vy -5;vx= vx + 5",
-        "vx=vx - 5",
-        "vx=vx - 5",
-        "vx=vx - 5",
-        "vx = vx + 5",
-        "vx = vx + 5",
-        "vx = vx + 5"
-    ]
-}
\ No newline at end of file
diff --git a/cartoy.c b/cartoy.c
deleted file mode 100644
index 9c9bcee6b61fd0c73392df5b7f196348d91674d4..0000000000000000000000000000000000000000
--- a/cartoy.c
+++ /dev/null
@@ -1,41 +0,0 @@
-enum modes {go_forward, move_over};
-
-struct State {
-int vx;
-int vy;
-int dist;
-//enum modes mode;
-};
-
-	
-	
-struct State P(struct State s) {
-	int posx = s.posx;
-	int posy = s.posy;
-	int vx = s.vx;
-	int vy = s.vy;
-	int dist = s.dist
-	enum modes state = s.mode;
-	if (dist < 10) {
-		if (dist < 2) {
-			if (state != move_over) {
-			//if super close, slow down and begin move over if we aren't already
-                        	vy=vy -5;
-                		vx= vx + 5;
-			}
-		}
-		else if (dist > 2) {
-			vx=vx - 5;
-		}
-	
-	}
-	else if (state == move_over) {
-		vx = vx + 5;
-	
-	}
-	s.vx = vx;
-	s.vy = vy;
-	s.mode = state;
-	return s;
-
-}
diff --git a/example_car_lane_switch.py b/example_car_lane_switch.py
index 4d82cba6834770d9474ae28b1063d49d57e2960c..5861a91a770dbdd88bd17a8476fe275c2e507734 100644
--- a/example_car_lane_switch.py
+++ b/example_car_lane_switch.py
@@ -9,44 +9,60 @@ class VehicleMode(Enum):
 class LaneMode(Enum):
     Lane0 = auto()
 
-def controller(x,y,theta,v,vehicle_mode, lane_mode):
-    output_vehicle_mode = vehicle_mode
-    output_lane_mode = lane_mode
-    if vehicle_mode == VehicleMode.Normal:
-        if lane_mode == LaneMode.Lane0:
-            if x > 3 and x < 5:
+class State:
+    x = 0.0
+    y = 0.0
+    theta = 0.0
+    v = 0.0
+    vehicle_mode: VehicleMode = VehicleMode.Normal
+    lane_mode: LaneMode = LaneMode.Lane0
+
+    def __init__(self, x, y, theta, v, vehicle_mode: VehicleMode, lane_mode: LaneMode):
+        self.data = []
+
+def controller(ego:State):
+    output_vehicle_mode = ego.vehicle_mode
+    output_lane_mode = ego.lane_mode
+    if ego.vehicle_mode == VehicleMode.Normal:
+        if ego.lane_mode == LaneMode.Lane0:
+            if ego.x > 3 and ego.x < 5:
                 output_vehicle_mode = VehicleMode.SwitchLeft
-            if x > 3 and x < 5:
+            if ego.x > 3 and ego.x < 5:
                 output_vehicle_mode = VehicleMode.SwitchRight
-    if vehicle_mode == VehicleMode.SwitchLeft:
-        if lane_mode == LaneMode.Lane0:
-            if x > 10:
+    if ego.vehicle_mode == VehicleMode.SwitchLeft:
+        if ego.lane_mode == LaneMode.Lane0:
+            if ego.x > 10:
                 output_vehicle_mode = VehicleMode.Normal
-    if vehicle_mode == VehicleMode.SwitchRight:
-        if lane_mode == LaneMode.Lane0:
-            if x > 10:
+    if ego.vehicle_mode == VehicleMode.SwitchRight:
+        if ego.lane_mode == LaneMode.Lane0:
+            if ego.x > 10:
                 output_vehicle_mode = VehicleMode.Normal
 
     return output_vehicle_mode, output_lane_mode
         
-from ourtool.agents.car_agent import CarAgent
-from ourtool.scenario.scenario import Scenario
-import matplotlib.pyplot as plt 
+from src.example.example_agent.car_agent import CarAgent
+from src.scene_verifier.scenario.scenario import Scenario
+from src.example.example_map.simple_map import SimpleMap2
+from src.example.example_sensor.fake_sensor import FakeSensor1
+import matplotlib.pyplot as plt
 import numpy as np
 
 if __name__ == "__main__":
-    input_code_name = 'car_lane_switch.py'
+    input_code_name = 'example_car_lane_switch.py'
     scenario = Scenario()
     
     car = CarAgent('ego', file_name=input_code_name)
     scenario.add_agent(car)
+    scenario.add_map(SimpleMap2())
+    scenario.set_sensor(FakeSensor1())
     
     # simulator = Simulator()
+    scenario.set_init(
+        [[0,3,0,0.5]],
+        [(VehicleMode.Normal, LaneMode.Lane0)]
+    )
+
     traces = scenario.simulate(
-        [[0,0,0,0.5]],
-        [(VehicleMode.Normal, LaneMode.Lane0)],
-        [car],
-        scenario,
         40
     )
 
diff --git a/example_two_car_lane_switch.py b/example_two_car_lane_switch.py
index 14b3d14dba6048d3e6cccc756669d6ba8b3403da..c270eb0ca16cbc11f68a301d4c0fa4dd2280fc44 100644
--- a/example_two_car_lane_switch.py
+++ b/example_two_car_lane_switch.py
@@ -1,4 +1,5 @@
-from enum import Enum,auto
+from enum import Enum, auto
+
 
 class VehicleMode(Enum):
     Normal = auto()
@@ -6,76 +7,105 @@ class VehicleMode(Enum):
     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):
+    def __init__(self, x, y, theta, v, vehicle_mode: VehicleMode, lane_mode: LaneMode):
         self.data = []
 
-def controller(ego:State, other:State, map):
+
+def controller(ego: State, other: State, lane_map):
     output = ego
     if ego.vehicle_mode == VehicleMode.Normal:
-        if ego.lane_mode == LaneMode.Lane0:
-            if other.x - ego.x > 3 and other.x - ego.x < 5 and map.has_left(ego.lane_mode):
+        # A simple example to demonstrate how our tool can handle change in controller
+        # if ego.x > 30 and ego.lane_mode == LaneMode.Lane0:
+        #     output.vehicle_mode = VehicleMode.SwitchRight
+        
+        if other.x - ego.x > 3 and other.x - ego.x < 5 and ego.lane_mode == other.lane_mode:
+            if lane_map.has_left(ego.lane_mode):
                 output.vehicle_mode = VehicleMode.SwitchLeft
-                output.lane_mode = map.left_lane(ego.lane_mode)
-            if other.x - ego.x > 3 and other.x - ego.x < 5:
+        if other.x - ego.x > 3 and other.x - ego.x < 5 and ego.lane_mode == other.lane_mode:
+            if lane_map.has_right(ego.lane_mode):
                 output.vehicle_mode = VehicleMode.SwitchRight
     if ego.vehicle_mode == VehicleMode.SwitchLeft:
-        if ego.lane_mode == LaneMode.Lane0:
-            if ego.x - other.x > 10:
-                output.vehicle_mode = VehicleMode.Normal
+        if  lane_map.lane_geometry(ego.lane_mode) - ego.y <= -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.lane_mode == LaneMode.Lane0:
-            if ego.x - other.x > 10:
-                output.vehicle_mode = VehicleMode.Normal
+        if lane_map.lane_geometry(ego.lane_mode)-ego.y >= 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.sensor import SimpleSensor
-from user.map import SimpleMap
-import matplotlib.pyplot as plt 
-import numpy as np
+
+
+from src.example.example_agent.car_agent import CarAgent
+from src.scene_verifier.scenario.scenario import Scenario
+from src.example.example_map.simple_map import SimpleMap2
+from src.plotter.plotter2D import plot_tree
+from src.example.example_sensor.fake_sensor import FakeSensor2
+
+import matplotlib.pyplot as plt
 
 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(SimpleMap())
-    # scenario.set_sensor(SimpleSensor())
+    scenario.add_map(SimpleMap2())
+    scenario.set_sensor(FakeSensor2())
     scenario.set_init(
-        [[0,0,0,1.0], [10,0,0,0.5]],
         [
-            (VehicleMode.Normal, LaneMode.Lane0),
-            (VehicleMode.Normal, LaneMode.Lane0)
+            [[10, 0, 0, 0.5],[10, 0, 0, 0.5]], 
+            [[-0.2, -0.2, 0, 1.0],[0.2, 0.2, 0, 1.0]],
+        ],
+        [
+            (VehicleMode.Normal, LaneMode.Lane1),
+            (VehicleMode.Normal, LaneMode.Lane1)
         ]
     )
     # simulator = Simulator()
-    traces = scenario.simulate(40)
-
-    queue = [traces]
-    while queue!=[]:
-        node = queue.pop(0)
-        traces = node.trace
-        agent_id = 'ego'
-        # for agent_id in traces:
-        trace = np.array(traces[agent_id])
-        plt.plot(trace[:,0], trace[:,2], 'b')
-        # if node.child != []:
-        queue += node.child 
+    # traces = scenario.simulate(40)
+    traces = scenario.verify(40)
+
+    fig = plt.figure()
+    fig = plot_tree(traces, 'car1', 1, [2], 'b', fig)
+    fig = plot_tree(traces, 'car2', 1, [2], 'r', fig)
+
     plt.show()
+
+    # 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 != []:
+    #     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')
+
+    #     agent_id = 'car1'
+    #     trace = np.array(traces[agent_id])
+    #     plt.plot(trace[:, 1], trace[:, 2], 'b')
+
+    #     # if node.child != []:
+    #     queue += node.child
+    # plt.show()
diff --git a/full_run b/full_run
deleted file mode 100644
index d95c8c760cdf68d3769218a5d627b426853ed023..0000000000000000000000000000000000000000
--- a/full_run
+++ /dev/null
@@ -1,27 +0,0 @@
-#!/bin/bash
-
-START="$(date +%s.%N)"
-
-
-input_code=$1
-input_info=$2
-#output should have path to correct place in dryvr directory
-output=$3
-
-GSTART="$(date +%s.%N)"
-
-python generateGraph-new.py $input_code $input_info $output
-
-GEND="$(date +%s.%N)"
-
-python main.py $output
-
-END="$(date +%s.%N)"
-#DURATION=$[ $'expr (date +%s.%N) - ${START}' ]
-DURATION=$( echo "$END - $START" | bc -l  )
-GDURATION=$( echo "$GEND - $GSTART" | bc -l  )
-
-echo "Running time of pipeline is"
-echo ${DURATION}
-echo "Running time of model generation is"
-echo ${GDURATION}
diff --git a/generateGraph-new.py b/generateGraph-new.py
deleted file mode 100644
index 849f1464f73a056d863fa1631e7942622866936c..0000000000000000000000000000000000000000
--- a/generateGraph-new.py
+++ /dev/null
@@ -1,483 +0,0 @@
-#NEW VERSION OF CODE WHICH ANALYZES PATHS AND COLLECTS VARAIBLES ALONG THE WAY
-
-import clang.cindex
-import typing
-import json
-import sys
-#from cfg import CFG
-#import clang.cfg
-
-
-#index = clang.cindex.Index.create()
-#translation_unit = index.parse('my_source.c', args=['-std=c99'])
-
-#print all the cursor types
-#for i in translation_unit.get_tokens(extent=translation_unit.cursor.extent):
-#    print (i.kind)
-
-#get the structs
-def filter_node_list_by_node_kind(
-    nodes: typing.Iterable[clang.cindex.Cursor],
-    kinds: list
-) -> typing.Iterable[clang.cindex.Cursor]:
-    result = []
-    for i in nodes:
-        if i.kind in kinds:
-            result.append(i)
-    return result
-
-#exclude includes
-def filter_node_list_by_file(
-    nodes: typing.Iterable[clang.cindex.Cursor],
-    file_name: str
-) -> typing.Iterable[clang.cindex.Cursor]:
-    result = []
-    for i in nodes:
-        if i.location.file.name == file_name:
-            result.append(i)
-    return result
-
-
-#all_classes = filter_node_list_by_node_kind(translation_unit.cursor.get_children(), [clang.cindex.CursorKind.ENUM_DECL, clang.cindex.CursorKind.STRUCT_DECL])
-#for i in all_classes:
-#    print (i.spelling)
-
-#captured fields in classes
-def is_exposed_field(node):return node.access_specifier == clang.cindex.AccessSpecifier.PUBLIC
-def find_all_exposed_fields(
-    cursor: clang.cindex.Cursor
-):
-    result = []
-    field_declarations = filter_node_list_by_node_kind(cursor.get_children(), [clang.cindex.CursorKind.FIELD_DECL])
-    for i in field_declarations:
-        if not is_exposed_field(i):
-            continue
-        result.append(i.displayname)
-    return result
-
-
-def populate_field_list_recursively(class_name: str,class_field_map,class_inheritance_map):
-    field_list = class_field_map.get(class_name)
-    if field_list is None:
-        return []
-    baseClasses = class_inheritance_map[class_name]
-    for i in baseClasses:
-        field_list = populate_field_list_recursively(i.spelling) + field_list
-    return field_list
-#rtti_map = {}
-
-
-
-def get_state_variables(translation_unit):
-    rtti_map = {}
-    source_nodes = filter_node_list_by_file(translation_unit.cursor.get_children(), translation_unit.spelling)
-    all_classes = filter_node_list_by_node_kind(source_nodes, [clang.cindex.CursorKind.ENUM_DECL, clang.cindex.CursorKind.STRUCT_DECL])
-    class_inheritance_map = {}
-    class_field_map = {}
-    for i in all_classes:
-        bases = []
-        for node in i.get_children():
-            if node.kind == clang.cindex.CursorKind.CXX_BASE_SPECIFIER:
-                referenceNode = node.referenced
-                bases.append(node.referenced)
-        class_inheritance_map[i.spelling] = bases
-    for i in all_classes:
-        fields = find_all_exposed_fields(i)
-        class_field_map[i.spelling] = fields
-    #print("foo")
-    for class_name, class_list in class_inheritance_map.items():
-        rtti_map[class_name] = populate_field_list_recursively(class_name,class_field_map,class_inheritance_map)
-    for class_name, field_list in rtti_map.items():
-        rendered_fields = []
-        for f in field_list:
-            rendered_fields.append(f)
-    return rendered_fields
-
-
-
-#def populate_field_list_recursively(class_name: str):
-#    field_list = class_field_map.get(class_name)
-#    if field_list is None:
-#        return []
-#    baseClasses = class_inheritance_map[class_name]
-#    for i in baseClasses:
-#        field_list = populate_field_list_recursively(i.spelling) + field_list
-#    return field_list
-#rtti_map = {}
-
-#def get_state_variables():
-#    print("foo")
-#    for class_name, class_list in class_inheritance_map.items():
-#        rtti_map[class_name] = populate_field_list_recursively(class_name)
-#    for class_name, field_list in rtti_map.items():
-#        rendered_fields = []
-#        for f in field_list:
-#            rendered_fields.append(f)
-#    return rendered_fields
-
-
-
-#how can we get the structure with the if statements
-#need to get ahold of if statement, then using if statement methods, we can get the inner statement
-
-
-#all_ifs = filter_node_list_by_node_kind(translation_unit.cursor.get_children(), [clang.cindex.CursorKind.IF_STMT])
-#for i in all_ifs:
-#    print (i.spelling)
-
-#print(translation_unit.spelling)
-
-def find_ifstmt(node):
-    if node.kind.is_statement():
-        #ref_node = clang.cindex.Cursor_ref(node)
-        #print ("found %s"" % node.location.line)
-        print(node.displayname)
-        print(node.location.line)
-    for c in node.get_children():
-        find_ifstmt(c)
-
-#finds the if statement, but cant get anything out of it
-#find_ifstmt(translation_unit.cursor)
-
-
-def rectree(node, indent):
-    print("%s item %s of kind %s" % (indent, node.spelling, node.kind))
-    for c in node.get_children():
-        rectree(c, indent + "  ")
-#goes through every part of code but doesn't give info about non variable names
-#rectree(translation_unit.cursor, "")
-
-#CursorKind.IF_STMT
-
-#given the if statement as node, visit each of the else if statements
-def visit_elif(node):
-    print("visiting inner elseif statement")
-    for i in node.get_children():
-        #print(i.kind)
-        #print(i.spelling)
-        if i.kind == clang.cindex.CursorKind.IF_STMT or i.kind == clang.cindex.CursorKind.COMPOUND_STMT:
-            print("found el-if of compound")
-            print(i.spelling)
-            visit_elif(i)
-    print("finished elseif statement")    
-#if (node.hasElseStorage()):
-        #print("foo")
-
-
-#TODO, but a method to get us the info from an if statemen
-def visitif(node):
-    print(node.spelling)
-
-def visit_inner_if(node):
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            print("inner if %s" % i.kind)
-            visit_elif(node)
-        else:
-            visit_inner_if(i)
-
-def get_cursor_id(cursor, cursor_list = []):
-    if cursor is None:
-        return None
-
-    # FIXME: This is really slow. It would be nice if the index API exposed
-    # something that let us hash cursors.
-    for i,c in enumerate(cursor_list):
-        if cursor == c:
-            return i
-    cursor_list.append(cursor)
-    return len(cursor_list) - 1
-
-def get_info(node, depth=0):
-    children = [get_info(c, depth+1)
-                    for c in node.get_children()]
-    return { 'id' : get_cursor_id(node),
-             'kind' : node.kind,
-             'usr' : node.get_usr(),
-             'spelling' : node.spelling,
-             'location' : node.location,
-             'extent.start' : node.extent.start,
-             'extent.end' : node.extent.end,
-             'is_definition' : node.is_definition(),
-             'definition id' : get_cursor_id(node.get_definition())}
-             #'children' : children }
-
-def code_from_cursor(cursor):
-    code = []
-    line = ""
-    prev_token = None
-    for tok in cursor.get_tokens():
-        if prev_token is None:
-            prev_token = tok
-        prev_location = prev_token.location
-        prev_token_end_col = prev_location.column + len(prev_token.spelling)
-        cur_location = tok.location
-        if cur_location.line > prev_location.line:
-            code.append(line)
-            line = " " * (cur_location.column - 1)
-        else:
-            if cur_location.column > (prev_token_end_col):
-                line += " "
-        line += tok.spelling
-        prev_token = tok
-    if len(line.strip()) > 0:
-        code.append(line)
-    return code
-
-
-def visit_outter_if(node):
-    if_statements = []
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            #print("Outter if statement %s" % i.kind)
-            #visit_inner_if(i)
-            if_statements.append(code_from_cursor(i)[0])
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = visit_outter_if(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-#visit_outter_if(translation_unit.cursor)
-
-
-def get_outter_if_cursors(node):
-    if_statements = []
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            #print("Outter if statement %s" % i.kind)
-            #visit_inner_if(i)
-            if_statements.append(i)
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = get_outter_if_cursors(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-
-def get_inner_if_cursors(node):
-    if_statements = []
-    for i in node.get_children():
-        #print(i.kind)
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            if_statements.append(i)
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = get_inner_if_cursors(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-
-
-#input: current node, the set of guards/conditionals up to this point, the set of resets
-#return: the sets of guards/resets for all the children
-def traverse_tree(node, guards, resets, indent,hasIfParent):
-    #print(guards)
-    childrens_guards=[]
-    childrens_resets=[]
-    results = []
-    found = []
-    #print(resets)
-    #if (indent== '------'):
-    #    print(guards)
-    #    return [guards,resets]
-    #if (node.kind == clang.cindex.CursorKind.RETURN_STMT):
-    #    return [guards, resets]
-    #elseifstatement = False
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT: #TODO: add else statements
-            temp_results = traverse_tree(i, guards, resets, indent+ "-", True)
-            for item in temp_results:
-                found.append(item)
-            #childrens_guards.append(result_guards)
-            #childrens_resets.append(result_resets)
-
-        else:
-            reset_copy = resets
-            guard_copy = guards
-            if hasIfParent:
-                if node.kind == clang.cindex.CursorKind.IF_STMT and i.kind != clang.cindex.CursorKind.COMPOUND_STMT:
-                    childrens_guards.append(code_from_cursor(i))
-                    #print(indent + "Guard statement")
-                    #print(code_from_cursor(i))
-                    #print(get_info(i))
-                    #found.append(code_from_cursor(i))
-                    #temp_results = traverse_tree(i, guard_copy, reset_copy, indent+ "-", hasIfParent)
-                    #for result in temp_results:
-                    #    result.append(code_from_cursor(i))
-                    #    results.append(result)
-                    #if len(temp_results)== 0:
-                    #    result.append([code_from_cursor(i)])
-                    #    print("foo")
-                if node.kind == clang.cindex.CursorKind.COMPOUND_STMT and i.kind != clang.cindex.CursorKind.DECL_STMT:
-                    #print(indent + "Reset statement")
-                    #print(code_from_cursor(i))
-                    childrens_resets.append(code_from_cursor(i))
-                    #found.append(code_from_cursor(i))
-                    #temp_results = traverse_tree(i, guard_copy, reset_copy, indent+ "-", hasIfParent)
-                    #for result in temp_results:
-                    #    result.append(code_from_cursor(i))
-                    #    results.append(result)
-                    #if len(temp_results) == 0:
-                    #    results.append([code_from_cursor(i)])
-                    #    print("foo")
-                #else:
-            temp_results = traverse_tree(i, guard_copy, reset_copy, indent+ "-", hasIfParent)
-            for item in temp_results:
-                 found.append(item)
-            #childrens_guards.append(results[0])
-            #childrens_resets.append(results[1])
-    #if my parent is an if statement and i'm and else if statement, I want to add the negation of my parent, not myself
-    if len(found) == 0 and len(childrens_resets) > 0:
-        found.append([])
-    for item in found:
-        for reset in childrens_resets:
-        #for item in found:
-            item.append([reset, 'r'])
-        results.append(item)
-        for guard in childrens_guards:
-            item.append([guard, 'g'])
-
-    return results
-
-#assumption, word mode is not in states
-def pretty_vertices(vertices):
-    output = []
-    for vertex_code in vertices:
-        parts = vertex_code.split("==")
-        nonwhitespace = parts[-1].split()
-        if "mode" not in nonwhitespace:
-            output.append(nonwhitespace[0].strip(')'))
-    return output
- 
-def get_next_state(code):
-    line = code[-2]
-    parts = line.split("=")
-    state = parts[-1].strip(';')
-    return state.strip()
-
-
-#TODO: And( , )
-def pretty_guards(guards):
-    output = ""
-    first = True
-    for condition in guards: 
-        #print(type(condition))
-        if first:
-            output+= condition[0]
-        else:
-            output = "And(" + condition[0] + ",(" + output + "))"
-        first = False
-    return output
-
-#assumption: reset;reset;reset
-def pretty_resets(resets):
-    outstring = ""
-    for reset in resets:
-       outstring+=reset[0] + ';'
-    return outstring.strip(';')
-
-
-##main code###
-#print(sys.argv)
-if len(sys.argv) < 4:
-    print("incorrect usage. call createGraph.py program inputfile outputfilename")
-    quit()
-
-input_code_name = sys.argv[1]
-input_file_name = sys.argv[2] 
-output_file_name = sys.argv[3]
-
-
-
-with open(input_file_name) as in_json_file:
-    input_json = json.load(in_json_file)
-
-output_dict = {
-}
-
-output_dict.update(input_json)
-
-#file parsing set up
-index = clang.cindex.Index.create()
-translation_unit = index.parse(input_code_name, args=['-std=c99'])
-print("testing cfg")
-
-#FunctionDecl* func = /* ... */ ;
-#ASTContext* context = /* ... */ ;
-#Stmt* funcBody = func->getBody();
-#CFG* sourceCFG = CFG::buildCFG(func, funcBody, context, clang::CFG::BuildOptions());
-
-#add each variable in state struct as variable
-variables = get_state_variables(translation_unit)
-output_dict['variables'] = variables
-
-#traverse the tree for all the paths
-paths = traverse_tree(translation_unit.cursor, [],[], "", False)
-
-vertices = []
-counter = 0
-for path in paths:
-    vertices.append(str(counter))
-    counter += 1
-    print(path)
-output_dict['vertex'] = vertices
-
-
-
-
-#traverse outter if statements and find inner statements
-edges = []
-guards = []
-resets = []
-
-counter = 0
-for path in paths:
-    guard = []
-    reset = []
-    for item in path:
-        if item[1] == 'g':
-            guard.append(item[0])
-        elif item[1] == 'r':
-            reset.append(item[0])
-    #create an edge from all other nodes to me
-    for vertex in vertices:
-        edges.append([vertex, str(counter)])
-        guards.append(pretty_guards(guard))
-        resets.append(pretty_resets(reset))
-    counter+= 1
-
-#add edge, transition(guards) and resets
-output_dict['edge'] = edges
-output_dict['guards'] = guards
-output_dict['resets'] = resets
-
-output_json = json.dumps(output_dict, indent=4)
-#print(output_json)
-outfile = open(output_file_name, "w")
-outfile.write(output_json)
-outfile.close()
-
-print("wrote json to " + output_file_name)
-
-
-
-
-
-
-
-
-
-
diff --git a/generateGraph.py b/generateGraph.py
deleted file mode 100644
index 5fb5919f1d2f62dab80d86fb82dd0b6704f5b0aa..0000000000000000000000000000000000000000
--- a/generateGraph.py
+++ /dev/null
@@ -1,401 +0,0 @@
-#OLD CODE WHICH ONLY ALLOWS TWO LEVELS OF IF STATEMENTS AND EXPLICT MODES 
-#NOT THE UPDATED CODE
-
-import clang.cindex
-import typing
-import json
-import sys
-
-#index = clang.cindex.Index.create()
-#translation_unit = index.parse('my_source.c', args=['-std=c99'])
-
-#print all the cursor types
-#for i in translation_unit.get_tokens(extent=translation_unit.cursor.extent):
-#    print (i.kind)
-
-#get the structs
-def filter_node_list_by_node_kind(
-    nodes: typing.Iterable[clang.cindex.Cursor],
-    kinds: list
-) -> typing.Iterable[clang.cindex.Cursor]:
-    result = []
-    for i in nodes:
-        if i.kind in kinds:
-            result.append(i)
-    return result
-
-#exclude includes
-def filter_node_list_by_file(
-    nodes: typing.Iterable[clang.cindex.Cursor],
-    file_name: str
-) -> typing.Iterable[clang.cindex.Cursor]:
-    result = []
-    for i in nodes:
-        if i.location.file.name == file_name:
-            result.append(i)
-    return result
-
-
-#all_classes = filter_node_list_by_node_kind(translation_unit.cursor.get_children(), [clang.cindex.CursorKind.ENUM_DECL, clang.cindex.CursorKind.STRUCT_DECL])
-#for i in all_classes:
-#    print (i.spelling)
-
-#captured fields in classes
-def is_exposed_field(node):return node.access_specifier == clang.cindex.AccessSpecifier.PUBLIC
-def find_all_exposed_fields(
-    cursor: clang.cindex.Cursor
-):
-    result = []
-    field_declarations = filter_node_list_by_node_kind(cursor.get_children(), [clang.cindex.CursorKind.FIELD_DECL])
-    for i in field_declarations:
-        if not is_exposed_field(i):
-            continue
-        result.append(i.displayname)
-    return result
-
-
-def populate_field_list_recursively(class_name: str,class_field_map,class_inheritance_map):
-    field_list = class_field_map.get(class_name)
-    if field_list is None:
-        return []
-    baseClasses = class_inheritance_map[class_name]
-    for i in baseClasses:
-        field_list = populate_field_list_recursively(i.spelling) + field_list
-    return field_list
-#rtti_map = {}
-
-
-
-def get_state_variables(translation_unit):
-    rtti_map = {}
-    source_nodes = filter_node_list_by_file(translation_unit.cursor.get_children(), translation_unit.spelling)
-    all_classes = filter_node_list_by_node_kind(source_nodes, [clang.cindex.CursorKind.ENUM_DECL, clang.cindex.CursorKind.STRUCT_DECL])
-    class_inheritance_map = {}
-    class_field_map = {}
-    for i in all_classes:
-        bases = []
-        for node in i.get_children():
-            if node.kind == clang.cindex.CursorKind.CXX_BASE_SPECIFIER:
-                referenceNode = node.referenced
-                bases.append(node.referenced)
-        class_inheritance_map[i.spelling] = bases
-    for i in all_classes:
-        fields = find_all_exposed_fields(i)
-        class_field_map[i.spelling] = fields
-    #print("foo")
-    for class_name, class_list in class_inheritance_map.items():
-        rtti_map[class_name] = populate_field_list_recursively(class_name,class_field_map,class_inheritance_map)
-    for class_name, field_list in rtti_map.items():
-        rendered_fields = []
-        for f in field_list:
-            rendered_fields.append(f)
-    return rendered_fields
-
-
-
-#def populate_field_list_recursively(class_name: str):
-#    field_list = class_field_map.get(class_name)
-#    if field_list is None:
-#        return []
-#    baseClasses = class_inheritance_map[class_name]
-#    for i in baseClasses:
-#        field_list = populate_field_list_recursively(i.spelling) + field_list
-#    return field_list
-#rtti_map = {}
-
-#def get_state_variables():
-#    print("foo")
-#    for class_name, class_list in class_inheritance_map.items():
-#        rtti_map[class_name] = populate_field_list_recursively(class_name)
-#    for class_name, field_list in rtti_map.items():
-#        rendered_fields = []
-#        for f in field_list:
-#            rendered_fields.append(f)
-#    return rendered_fields
-
-
-
-#how can we get the structure with the if statements
-#need to get ahold of if statement, then using if statement methods, we can get the inner statement
-
-
-#all_ifs = filter_node_list_by_node_kind(translation_unit.cursor.get_children(), [clang.cindex.CursorKind.IF_STMT])
-#for i in all_ifs:
-#    print (i.spelling)
-
-#print(translation_unit.spelling)
-
-def find_ifstmt(node):
-    if node.kind.is_statement():
-        #ref_node = clang.cindex.Cursor_ref(node)
-        #print ("found %s"" % node.location.line)
-        print(node.displayname)
-        print(node.location.line)
-    for c in node.get_children():
-        find_ifstmt(c)
-
-#finds the if statement, but cant get anything out of it
-#find_ifstmt(translation_unit.cursor)
-
-
-def rectree(node, indent):
-    print("%s item %s of kind %s" % (indent, node.spelling, node.kind))
-    for c in node.get_children():
-        rectree(c, indent + "  ")
-#goes through every part of code but doesn't give info about non variable names
-#rectree(translation_unit.cursor, "")
-
-#CursorKind.IF_STMT
-
-#given the if statement as node, visit each of the else if statements
-def visit_elif(node):
-    print("visiting inner elseif statement")
-    for i in node.get_children():
-        #print(i.kind)
-        #print(i.spelling)
-        if i.kind == clang.cindex.CursorKind.IF_STMT or i.kind == clang.cindex.CursorKind.COMPOUND_STMT:
-            print("found el-if of compound")
-            print(i.spelling)
-            visit_elif(i)
-    print("finished elseif statement")    
-#if (node.hasElseStorage()):
-        #print("foo")
-
-
-#TODO, but a method to get us the info from an if statemen
-def visitif(node):
-    print(node.spelling)
-
-def visit_inner_if(node):
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            print("inner if %s" % i.kind)
-            visit_elif(node)
-        else:
-            visit_inner_if(i)
-
-def get_cursor_id(cursor, cursor_list = []):
-    if cursor is None:
-        return None
-
-    # FIXME: This is really slow. It would be nice if the index API exposed
-    # something that let us hash cursors.
-    for i,c in enumerate(cursor_list):
-        if cursor == c:
-            return i
-    cursor_list.append(cursor)
-    return len(cursor_list) - 1
-
-def get_info(node, depth=0):
-    children = [get_info(c, depth+1)
-                    for c in node.get_children()]
-    return { 'id' : get_cursor_id(node),
-             'kind' : node.kind,
-             'usr' : node.get_usr(),
-             'spelling' : node.spelling,
-             'location' : node.location,
-             'extent.start' : node.extent.start,
-             'extent.end' : node.extent.end,
-             'is_definition' : node.is_definition(),
-             'definition id' : get_cursor_id(node.get_definition())}
-             #'children' : children }
-
-def code_from_cursor(cursor):
-    code = []
-    line = ""
-    prev_token = None
-    for tok in cursor.get_tokens():
-        if prev_token is None:
-            prev_token = tok
-        prev_location = prev_token.location
-        prev_token_end_col = prev_location.column + len(prev_token.spelling)
-        cur_location = tok.location
-        if cur_location.line > prev_location.line:
-            code.append(line)
-            line = " " * (cur_location.column - 1)
-        else:
-            if cur_location.column > (prev_token_end_col):
-                line += " "
-        line += tok.spelling
-        prev_token = tok
-    if len(line.strip()) > 0:
-        code.append(line)
-    return code
-
-
-def visit_outter_if(node):
-    if_statements = []
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            #print("Outter if statement %s" % i.kind)
-            #visit_inner_if(i)
-            if_statements.append(code_from_cursor(i)[0])
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = visit_outter_if(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-#visit_outter_if(translation_unit.cursor)
-
-
-def get_outter_if_cursors(node):
-    if_statements = []
-    for i in node.get_children():
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            #print("Outter if statement %s" % i.kind)
-            #visit_inner_if(i)
-            if_statements.append(i)
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = get_outter_if_cursors(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-
-def get_inner_if_cursors(node):
-    if_statements = []
-    for i in node.get_children():
-        #print(i.kind)
-        if i.kind == clang.cindex.CursorKind.IF_STMT:
-            if_statements.append(i)
-            #print(if_statements)
-            #print(i.spelling)
-            #for child in i.get_children():
-            #    print(str(i.kind) + str(i.spelling))
-        else:
-            out_ifs = get_inner_if_cursors(i)
-            for statement in out_ifs:
-                if_statements.append(statement)
-    return if_statements
-
-#assumption, word mode is not in states
-def pretty_vertices(vertices):
-    output = []
-    for vertex_code in vertices:
-        parts = vertex_code.split("==")
-        nonwhitespace = parts[-1].split()
-        if "mode" not in nonwhitespace:
-            output.append(nonwhitespace[0].strip(')'))
-    return output
- 
-def get_next_state(code):
-    line = code[-2]
-    parts = line.split("=")
-    state = parts[-1].strip(';')
-    return state.strip()
-
-
-#TODO: consider or??
-def pretty_guards(code):
-    line = code[0]
-    conditional = line.strip('if').strip('{')
-    conditions = conditional.split('&&')
-    output = "And"
-    for condition in conditions: 
-        output += condition.strip() + ","
-    output = output.strip(",")
-    return output
-
-#assumption: last two lines of code are reset and bracket... not idea
-def pretty_resets(code):
-    outstring = ""
-    for index in range(0,len(code)):
-       if index != 0 and index != len(code)-1 and index != len(code)-2:
-        outstring += code[index].strip().strip('\n')
-    return outstring.strip(';')
-
-##main code###
-#print(sys.argv)
-if len(sys.argv) < 4:
-    print("incorrect usage. call createGraph.py program inputfile outputfilename")
-    quit()
-
-input_code_name = sys.argv[1]
-input_file_name = sys.argv[2] 
-output_file_name = sys.argv[3]
-
-with open(input_file_name) as in_json_file:
-    input_json = json.load(in_json_file)
-
-output_dict = {
-}
-
-output_dict.update(input_json)
-
-#file parsing set up
-index = clang.cindex.Index.create()
-translation_unit = index.parse(input_code_name, args=['-std=c99'])
-
-#add each variable in state struct as variable
-variables = get_state_variables(translation_unit)
-#assumption mode variable
-variables.remove('mode')
-output_dict['variables'] = variables
-
-#capture each outter if statement and create state
-vertices = pretty_vertices(visit_outter_if(translation_unit.cursor))
-#print(vertice
-output_dict['vertex'] = vertices
-
-
-#traverse outter if statements and find inner statements
-edges = []
-guards = []
-resets = []
-
-#print(type(translation_unit.cursor))
-
-index = 0
-outter_ifs = get_outter_if_cursors(translation_unit.cursor)
-#print(outter_ifs)
-for if_statement in outter_ifs:
-    inner_statements = get_inner_if_cursors(if_statement)
-    for in_statement in inner_statements:
-        code = code_from_cursor(in_statement)
-        guards.append(pretty_guards(code)) 
-        #reset = ""
-        #for line in range(1, len(code)):
-        #    reset += str(code[line])
-        resets.append(pretty_resets(code))
-        to_edge = 0
-        #key assumption, second to last line gives next state, fix!
-        for i in range(0,len(vertices)):
-            #print(vertices[i])
-            #print(get_next_state(code))
-            if vertices[i] == get_next_state(code):
-                to_edge = i
-        edges.append([index,to_edge]) 
-    index+=1    
-
-#each inner if
-#add edge, transition(guards) and resets
-output_dict['edge'] = edges
-output_dict['guards'] = guards
-output_dict['resets'] = resets
-
-output_json = json.dumps(output_dict, indent=4)
-#print(output_json)
-outfile = open(output_file_name, "w")
-outfile.write(output_json)
-outfile.close()
-
-print("wrote json to " + output_file_name)
-
-
-
-
-
-
-
-
-
-
diff --git a/ourtool/automaton/__init__.py b/ourtool/automaton/__init__.py
deleted file mode 100644
index 171bea6f2daef09e5dd888ad0e79d32db6695628..0000000000000000000000000000000000000000
--- a/ourtool/automaton/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-# try:
-#     from hybrid_automaton import *
-# except:
-#     from ourtool.automaton.hybrid_automaton import *
\ No newline at end of file
diff --git a/ourtool/automaton/guard.py b/ourtool/automaton/guard.py
deleted file mode 100644
index 3792687ea45119d8e11026726e656116ab73833c..0000000000000000000000000000000000000000
--- a/ourtool/automaton/guard.py
+++ /dev/null
@@ -1,255 +0,0 @@
-import re
-from typing import List, Dict
-from ourtool.automaton.hybrid_io_automaton import HybridIoAutomaton
-
-class LogicTreeNode:
-    def __init__(self, data, child = [], val = None, mode_guard = None):
-        self.data = data 
-        self.child = child
-        self.val = val
-        self.mode_guard = mode_guard
-
-class GuardExpression:
-    def __init__(self, root:LogicTreeNode=None, logic_str:str=None):
-        self.logic_tree_root = root
-        self.logic_string = logic_str
-        if self.logic_tree_root is None and logic_str is not None:
-            self.construct_tree_from_str(logic_str)
-
-    def logic_string_split(self, logic_string):
-        # Input:
-        #   logic_string: str, a python logic expression
-        # Output:
-        #   List[str], a list of string containing atomics, logic operator and brackets
-        # The function take a python logic expression and split the expression into brackets, atomics and logic operators
-        # logic_string = logic_string.replace(' ','')
-        res = re.split('( and )',logic_string)
-
-        tmp = []
-        for sub_str in res:
-            tmp += re.split('( or )',sub_str)
-        res = tmp
-
-        tmp = []
-        for sub_str in res:
-            tmp += re.split('(\()',sub_str)
-        res = tmp
-
-        tmp = []
-        for sub_str in res:
-            tmp += re.split('(\))',sub_str)
-        res = tmp
-
-        while("" in res) :
-            res.remove("")
-        while(" " in res):
-            res.remove(" ")
-        for i,sub_str in enumerate(res):
-            res[i]= sub_str.strip(' ')
-
-        # Handle spurious brackets in the splitted string
-        # Get all the index of brackets pairs in the splitted string
-        # Construct brackets tree
-        # class BracketTreeNode:
-        #     def __init__(self):
-        #         self.left_idx = None 
-        #         self.right_idx = None 
-        #         self.child = []
-        bracket_stack = []
-        for i in range(len(res)):
-            if res[i] == "(":
-                bracket_stack.append(i)
-            elif res[i] == ")":
-                left_idx = bracket_stack.pop()
-                sub_list = res[left_idx:i+1]
-                # Check for each brackets pairs if there's any logic operators in between
-                # If no, combine things in between and the brackets together, reconstruct the list
-                if "and" not in sub_list and "or" not in sub_list:   
-                    res[left_idx] = "".join(sub_list)
-                    for j in range(left_idx+1,i+1):
-                        res[j] = ""
-
-        # For each pair of logic operator
-        start_idx = 0
-        end_idx = 0
-        for i in range(len(res)):
-            if res[i]!="(":
-                start_idx = i
-                break
-        
-        for i in range(len(res)):
-            if res[i] == "and" or res[i] == "or":
-                end_idx = i 
-                sub_list = res[start_idx:end_idx]
-                # Check if there's any dangling brackents in between. 
-                # If no, combine things between logic operators
-                if "(" not in sub_list and ")" not in sub_list:
-                    res[start_idx] = "".join(sub_list)
-                    for j in range(start_idx+1, end_idx):
-                        res[j] = ""
-                start_idx = end_idx + 1
-        while("" in res) :
-            res.remove("")
-        return res 
-
-    def construct_tree_from_str(self, logic_string:str):
-        # Convert an infix expression notation to an expression tree
-        # https://www.geeksforgeeks.org/program-to-convert-infix-notation-to-expression-tree/
-
-        self.logic_string = logic_string
-        logic_string = "(" + logic_string + ")"
-        s = self.logic_string_split(logic_string)
-
-        stN = []
-        stC = []
-        p = {}
-        p["and"] = 1
-        p["or"] = 1
-        p[")"] = 0
-
-        for i in range(len(s)):
-            if s[i] == "(":
-                stC.append(s[i])
-            
-            elif s[i] not in p:
-                t = LogicTreeNode(s[i])
-                stN.append(t)
-            
-            elif(p[s[i]]>0):
-                while (len(stC) != 0 and stC[-1] != '(' and p[stC[-1]] >= p[s[i]]):
-                                    # Get and remove the top element
-                    # from the character stack
-                    t = LogicTreeNode(stC[-1])
-                    stC.pop()
-    
-                    # Get and remove the top element
-                    # from the node stack
-                    t1 = stN[-1]
-                    stN.pop()
-    
-                    # Get and remove the currently top
-                    # element from the node stack
-                    t2 = stN[-1]
-                    stN.pop()
-    
-                    # Update the tree
-                    t.child = [t1, t2]
-    
-                    # Push the node to the node stack
-                    stN.append(t)
-                stC.append(s[i])
-            elif (s[i] == ')'):
-                while (len(stC) != 0 and stC[-1] != '('):
-                    # from the character stack
-                    t = LogicTreeNode(stC[-1])
-                    stC.pop()
-    
-                    # Get and remove the top element
-                    # from the node stack
-                    t1 = stN[-1]
-                    stN.pop()
-    
-                    # Get and remove the currently top
-                    # element from the node stack
-                    t2 = stN[-1]
-                    stN.pop()
-    
-                    # Update the tree
-                    t.child = [t1, t2]
-    
-                    # Push the node to the node stack
-                    stN.append(t)
-                stC.pop()
-        t = stN[-1]
-        self.logic_tree_root = t
-
-    def generate_guard_string_python(self):
-        return self._generate_guard_string_python(self.logic_tree_root)
-
-    def _generate_guard_string_python(self, root: LogicTreeNode)->str:
-        if root.data!="and" and root.data!="or":
-            return root.data
-        else:
-            data1 = self._generate_guard_string_python(root.child[0])
-            data2 = self._generate_guard_string_python(root.child[1])
-            return f"({data1} {root.data} {data2})"
-
-    def generate_guard_string(self):
-        return self._generate_guard_string(self.logic_tree_root)
-
-    def _generate_guard_string(self, root: LogicTreeNode)->str:
-        if root.data!="and" and root.data!="or":
-            return root.data
-        else:
-            data1 = self._generate_guard_string(root.child[0])
-            data2 = self._generate_guard_string(root.child[1])
-            if root.data == "and":
-                return f"And({data1},{data2})"
-            elif root.data == "or":
-                return f"Or({data1},{data2})"
-
-    def execute_guard(self, discrete_variable_dict:Dict) -> bool:
-        # This function will execute guard, and remove guard related to mode from the tree
-        # We can do this recursively
-        res = self._execute_guard(self.logic_tree_root, discrete_variable_dict)
-        
-        return res
-
-    def _execute_guard(self, root:LogicTreeNode, discrete_variable_dict:Dict) -> bool:
-        # If is tree leaf
-        if root.child == []:
-            # Check if the expression involves mode
-            expr = root.data
-            is_mode_guard = False
-            for key in discrete_variable_dict:
-                if key in expr:
-                    is_mode_guard = True
-                    expr = expr.replace(key, discrete_variable_dict[key])
-            if is_mode_guard:
-                # Execute guard, assign type and  and return result
-                root.mode_guard = True
-                expr = expr.strip('(')
-                expr = expr.strip(')')
-                expr = expr.replace(' ','')
-                expr = expr.split('==')
-                res = expr[0] == expr[1]
-                # res = eval(expr)
-                root.val = res 
-                return res
-            # Otherwise, return True
-            else: 
-                root.mode_guard = False 
-                root.val = True 
-                return True
-        # For the two children, call _execute_guard and collect result
-        res1 = self._execute_guard(root.child[0],discrete_variable_dict)
-        res2 = self._execute_guard(root.child[1],discrete_variable_dict)
-        # Evaluate result for current node
-        if root.data == "and":
-            res = res1 and res2 
-        elif root.data == "or":
-            res = res1 or res2
-        else:
-            raise ValueError(f"Invalid root data {root.data}")
-        
-        # If the result is False, return False
-        if not res:
-            return False 
-        # Else if any child have false result, remove that child
-        else: 
-            if not res1 or root.child[0].mode_guard:
-                root.data = root.child[1].data
-                root.val = root.child[1].val 
-                root.mode_guard = root.child[1].mode_guard
-                root.child = root.child[1].child
-            elif not res2 or root.child[1].mode_guard:
-                root.data = root.child[0].data
-                root.val = root.child[0].val 
-                root.mode_guard = root.child[0].mode_guard
-                root.child = root.child[0].child
-            return True
-
-if __name__ == "__main__":
-    tmp = GuardExpression()
-    tmp.construct_tree_from_str('(other_x-ego_x<20) and other_x-ego_x>10 and other_vehicle_lane==ego_vehicle_lane')
-    print("stop")
\ No newline at end of file
diff --git a/ourtool/map/__init__.py b/ourtool/map/__init__.py
deleted file mode 100644
index e7b9ea180897f28e93a2da6b3d4aec642a08938b..0000000000000000000000000000000000000000
--- a/ourtool/map/__init__.py
+++ /dev/null
@@ -1,8 +0,0 @@
-# try:
-#     from lane_map import LaneMap
-#     from single_straight_lane import SingleStraightLaneMap
-#     from lane_segment import LaneSegment
-# except:
-#     from ourtool.map.lane_segment import LaneSegment
-#     # from ourtool.map.lane_map import LaneMap
-#     # from ourtool.map.single_straight_lane import SingleStraightLaneMap
\ No newline at end of file
diff --git a/ourtool/map/lane_map.py b/ourtool/map/lane_map.py
deleted file mode 100644
index 30e2fe8cb309183247126978222d5fcd7a3e0408..0000000000000000000000000000000000000000
--- a/ourtool/map/lane_map.py
+++ /dev/null
@@ -1,37 +0,0 @@
-from typing import Dict 
-
-from ourtool.map.lane_segment import LaneSegment
-
-class LaneMap:
-    def __init__(self):
-        self.lane_segment_dict:Dict[int:LaneSegment] = {}
-
-    def get_left_lane_idx(self, lane_idx):
-        if lane_idx not in self.lane_segment_dict:
-            raise ValueError(f"lane_idx {lane_idx} not in lane_segment_dict")
-        lane_segment:LaneSegment = self.lane_segment_dict[lane_idx]
-        return lane_segment.left_lane
-        
-    def get_left_lane_segment(self,lane_idx):
-        left_lane_idx = self.get_left_lane_idx(lane_idx)
-        return self.lane_segment_dict[left_lane_idx]
-
-    def get_right_lane_idx(self, lane_idx):
-        if lane_idx not in self.lane_segment_dict:
-            raise ValueError(f"lane_idx {lane_idx} not in lane_segment_dict")
-        lane_segment:LaneSegment = self.lane_segment_dict[lane_idx]
-        return lane_segment.right_lane
-        
-    def get_right_lane_segment(self,lane_idx):
-        right_lane_idx = self.get_right_lane_idx(lane_idx)
-        return self.lane_segment_dict[right_lane_idx]
-
-    def get_next_lane_idx(self, lane_idx):
-        if lane_idx not in self.lane_segment_dict:
-            raise ValueError(f"lane_idx {lane_idx} not in lane_segment_dict")
-        lane_segment:LaneSegment = self.lane_segment_dict[lane_idx]
-        return lane_segment.next_segment
-        
-    def get_next_lane_segment(self,lane_idx):
-        next_lane_idx = self.get_next_lane_idx(lane_idx)
-        return self.lane_segment_dict[next_lane_idx]
diff --git a/ourtool/map/lane_segment.py b/ourtool/map/lane_segment.py
deleted file mode 100644
index ab541e5798aed68e296ccf74f78c2461d7018cfa..0000000000000000000000000000000000000000
--- a/ourtool/map/lane_segment.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from typing import List
-
-class LaneSegment:
-    def __init__(self, id, left_lane, right_lane, next_segment, lane_parameter = None):
-        self.id:int = id
-        self.left_lane:int = left_lane
-        self.right_lane:int = right_lane 
-        self.next_segment:int = next_segment
-
-        self.lane_parameter = None 
-        if lane_parameter is not None:
-            self.lane_parameter = lane_parameter
-
-        
\ No newline at end of file
diff --git a/ourtool/map/single_straight_lane.py b/ourtool/map/single_straight_lane.py
deleted file mode 100644
index e968f17fd0ca9fd5e929858cfc4d2b38bcbda8de..0000000000000000000000000000000000000000
--- a/ourtool/map/single_straight_lane.py
+++ /dev/null
@@ -1,8 +0,0 @@
-from ourtool.map.lane_map import LaneMap
-from ourtool.map.lane_segment import LaneSegment
-
-class SingleStraightLaneMap(LaneMap):
-    def __init__(self):
-        super().__init__()
-        segment = LaneSegment(0, None, None, None, None)
-        self.lane_segment_dict[segment.id] = segment
\ No newline at end of file
diff --git a/ourtool/scenario/__init__.py b/ourtool/scenario/__init__.py
deleted file mode 100644
index 76c3abd5dff7acdee36f5d23104f86e6ee455f19..0000000000000000000000000000000000000000
--- a/ourtool/scenario/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-# try:
-#     from Scenario import *
-# except:
-#     from ourtool.scenario.Scenario import *
\ No newline at end of file
diff --git a/ourtool/scenario/scenario.py b/ourtool/scenario/scenario.py
deleted file mode 100644
index d5ca98718ac65447b333f2eb9f0f27d9c1d87e60..0000000000000000000000000000000000000000
--- a/ourtool/scenario/scenario.py
+++ /dev/null
@@ -1,100 +0,0 @@
-from typing import Dict, List
-import copy
-
-from ourtool.agents.base_agent import BaseAgent
-from ourtool.automaton.guard import GuardExpression
-from pythonparser import Guard
-from pythonparser import Reset
-from ourtool.simulator.simulator import Simulator
-
-class Scenario:
-    def __init__(self):
-        self.agent_dict = {}
-        self.simulator = Simulator()
-        self.init_dict = {}
-        self.init_mode_dict = {}
-
-    def add_agent(self, agent:BaseAgent):
-        self.agent_dict[agent.id] = agent
-
-    def set_init(self, init_list, init_mode_list):
-        assert len(init_list) == len(self.agent_dict)
-        assert len(init_mode_list) == len(self.agent_dict)
-        for i,agent_id in enumerate(self.agent_dict.keys()):
-            self.init_dict[agent_id] = copy.deepcopy(init_list[i])
-            self.init_mode_dict[agent_id] = copy.deepcopy(init_mode_list[i])
-
-    def simulate(self, time_horizon):
-        init_list = []
-        init_mode_list = []
-        agent_list = []
-        for agent_id in self.agent_dict:
-            init_list.append(self.init_dict[agent_id])
-            init_mode_list.append(self.init_mode_dict[agent_id])
-            agent_list.append(self.agent_dict[agent_id])
-        return self.simulator.simulate(init_list, init_mode_list, agent_list, self, time_horizon)
-
-    def get_all_transition(self, state_dict):
-        guard_hit = False
-        satisfied_guard = []
-        for agent_id in state_dict:
-            agent:BaseAgent = self.agent_dict[agent_id]
-            agent_state, agent_mode = state_dict[agent_id]
-            t = agent_state[0]
-            agent_state = agent_state[1:]
-            paths = agent.controller.getNextModes(agent_mode)
-            for path in paths:
-                guard_list = []
-                reset_list = []
-                for item in path:
-                    if isinstance(item, Guard):
-                        guard_list.append("(" + item.code + ")")
-                    elif isinstance(item, Reset):
-                        reset_list.append(item.code)
-                guard_str = ' and '.join(guard_list)
-                guard_expression = GuardExpression(logic_str = guard_str)
-                # print(guard_expression.generate_guard_string_python())
-                discrete_variable_dict = {}
-                agent_mode_split = agent_mode.split(',')
-                assert len(agent_mode_split)==len(agent.controller.discrete_variables)
-                for dis_var,dis_val in zip(agent.controller.discrete_variables, agent_mode_split):
-                    for key in agent.controller.modes:
-                        if dis_val in agent.controller.modes[key]:
-                            tmp = key+'.'+dis_val
-                            break
-                    discrete_variable_dict[dis_var] = tmp
-                guard_can_satisfy = guard_expression.execute_guard(discrete_variable_dict)
-                if guard_can_satisfy:
-                    dryvr_guard_string = guard_expression.generate_guard_string_python()
-                    # Substitute value into dryvr guard string
-                    for i, variable in enumerate(agent.controller.variables):
-                        dryvr_guard_string = dryvr_guard_string.replace(variable, str(agent_state[i]))
-                    # Evaluate the guard strings 
-                    res = eval(dryvr_guard_string)
-                    # if result is true, check reset and construct next mode             
-                    if res:
-                        next_init = agent_state
-                        dest = agent_mode.split(',')
-                        for reset in reset_list:
-                            # Specify the destination mode
-                            if "mode" in reset:
-                                for i, discrete_variable in enumerate(agent.controller.discrete_variables):
-                                    if discrete_variable in reset:
-                                        break
-                                tmp = reset.split('=')
-                                tmp = tmp[1].split('.')
-                                if tmp[0].strip(' ') in agent.controller.modes:
-                                    dest[i] = tmp[1]
-                                
-                            else: 
-                                for i, cts_variable in enumerate(agent.controller.variables):
-                                    if cts_variable in reset:
-                                        break 
-                                tmp = reset.split('=')
-                                next_init[i] = float(tmp[1])
-                        dest = ','.join(dest)
-                        next_transition = (
-                            agent_id, agent_mode, dest, next_init, 
-                        )
-                        satisfied_guard.append(next_transition)
-        return satisfied_guard
\ No newline at end of file
diff --git a/out.json b/out.json
deleted file mode 100644
index 0984a68184e7819e8d0ce7dea15d944487635cac..0000000000000000000000000000000000000000
--- a/out.json
+++ /dev/null
@@ -1,21 +0,0 @@
-{
-    "initialVertex": 0,
-    "initialSet": [
-        [
-            1.5,
-            0.0,
-            1.0,
-            1.0
-        ],
-        [
-            1.51,
-            0.0,
-            1.0,
-            1.0
-        ]
-    ],
-    "unsafeSet": "@Allmode:posy>=12.0",
-    "timeHorizon": 16,
-    "directory": "examples/billiard",
-    "deterministic": true
-}
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2af2f8061800760758da46e37292b86ee2a26373
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,10 @@
+networkx~=2.2
+sympy~=1.6.2
+numpy~=1.22.1
+six~=1.14.0
+matplotlib~=3.4.2
+scipy~=1.8.0
+astunparse~=1.6.3
+treelib~=1.6.1
+polytope~=0.2.3
+pyvista~=0.32.1
\ No newline at end of file
diff --git a/singlevehiclesat-out.json b/singlevehiclesat-out.json
deleted file mode 100644
index 28aa3e5796fad7e36ee8b1bc1252975d2cccfa35..0000000000000000000000000000000000000000
--- a/singlevehiclesat-out.json
+++ /dev/null
@@ -1,89 +0,0 @@
-{
-    "unsafeSet": "",
-    "initialSet": [
-        [
-            0,
-            -0.2,
-            0
-        ],
-        [
-            0,
-            0.2,
-            0
-        ]
-    ],
-    "timeHorizon": 10,
-    "directory": "examples/curve_controller_hybrid",
-    "initialVertex": 0,
-    "deterministic": true,
-    "variables": [
-        "vx",
-        "vy",
-        "dist"
-    ],
-    "vertex": [
-        "0",
-        "1",
-        "2"
-    ],
-    "edge": [
-        [
-            "0",
-            "0"
-        ],
-        [
-            "1",
-            "0"
-        ],
-        [
-            "2",
-            "0"
-        ],
-        [
-            "0",
-            "1"
-        ],
-        [
-            "1",
-            "1"
-        ],
-        [
-            "2",
-            "1"
-        ],
-        [
-            "0",
-            "2"
-        ],
-        [
-            "1",
-            "2"
-        ],
-        [
-            "2",
-            "2"
-        ]
-    ],
-    "guards": [
-        "And(dist < 10,And((dist > 2),And(dist < 2,state != move_over)))",
-        "And(dist < 10,And((dist > 2),And(dist < 2,state != move_over)))",
-        "And(dist < 10,And((dist > 2),And(dist < 2,state != move_over)))",
-        "dist < 10",
-        "dist < 10",
-        "dist < 10",
-        "And(dist < 10,state == move_over)",
-        "And(dist < 10,state == move_over)",
-        "And(dist < 10,state == move_over)"
-    ],
-    "resets": [
-        "vy=vy - 5;vx= vx + 5;{",
-        "vy=vy - 5;vx= vx + 5;{",
-        "vy=vy - 5;vx= vx + 5;{",
-        "vx=vx - 5;{",
-        "vx=vx - 5;{",
-        "vx=vx - 5;{",
-        "vx = vx + 5",
-        "vx = vx + 5",
-        "vx = vx + 5"
-    ]
-}
\ No newline at end of file
diff --git a/singlevehiclesat.c b/singlevehiclesat.c
deleted file mode 100644
index a5d84ff81926110b2f81f7d918539349d46b2b30..0000000000000000000000000000000000000000
--- a/singlevehiclesat.c
+++ /dev/null
@@ -1,42 +0,0 @@
-enum modes {Normal,Sat_low,Sat_high};
-
-struct State {
-double tau;
-double yf;
-double thetaf;
-enum modes mode;
-};
-
-	
-	
-struct State P(struct State s) {
-	double tau = s.tau;
-	double yf = s.yf;
-	double thetaf = s.theta;
-	enum modes state = s.mode;
-	if (s.mode==Normal) {
-		if (-0.155914*yf-thetaf <= -0.60871) {
-			state=Sat_low;
-		}
-		if (-0.155914*yf-thetaf >= 0.60871) {
-                        state=Sat_high;
-                }
-	
-	}
-	if (s.mode ==Sat_low) {
-		if (-0.155914*yf-thetaf >= -0.60871) {
-                        state=Normal;
-                }
-	
-	}
-	if (s.mode == Sat_high) {
-		if (-0.155914*yf-thetaf <= 0.60871) {
-                        state=Normal;
-                }
-	}
-	
-	}
-	s.mode = state;
-	return s;
-
-}
diff --git a/singlevehiclesat.json b/singlevehiclesat.json
deleted file mode 100644
index de98f9d2eecf699a060e692af9506d920116aa4b..0000000000000000000000000000000000000000
--- a/singlevehiclesat.json
+++ /dev/null
@@ -1,8 +0,0 @@
-{
-"unsafeSet": "",
-"initialSet": [[0,-0.2,0, "Normal"],[0,0.2,0, "Normal"]],
-"timeHorizon": 10,
-"directory": "examples/curve_controller_hybrid",
-"initialVertex":0,
-"deterministic":true
-}
diff --git a/ourtool/__init__.py b/src/__init__.py
similarity index 100%
rename from ourtool/__init__.py
rename to src/__init__.py
diff --git a/ourtool/agents/__init__.py b/src/example/__init__.py
similarity index 100%
rename from ourtool/agents/__init__.py
rename to src/example/__init__.py
diff --git a/ourtool/simulator/__init__.py b/src/example/example_agent/__init__.py
similarity index 100%
rename from ourtool/simulator/__init__.py
rename to src/example/example_agent/__init__.py
diff --git a/ourtool/agents/car_agent.py b/src/example/example_agent/car_agent.py
similarity index 86%
rename from ourtool/agents/car_agent.py
rename to src/example/example_agent/car_agent.py
index e698ac6d3cc3fb24a00a941f5e993d80e71dba2a..f8adc8878dd8e36b533e1c5bb420cc867bce7510 100644
--- a/ourtool/agents/car_agent.py
+++ b/src/example/example_agent/car_agent.py
@@ -1,6 +1,7 @@
-from ourtool.agents.base_agent import BaseAgent
+from src.scene_verifier.agents.base_agent import BaseAgent
 import numpy as np 
 from scipy.integrate import ode
+from src.scene_verifier.map.lane_map import LaneMap
 
 class CarAgent(BaseAgent):
     def __init__(self, id, code = None, file_name = None):
@@ -16,7 +17,7 @@ class CarAgent(BaseAgent):
         v_dot = a 
         return [x_dot, y_dot, theta_dot, v_dot]
 
-    def TC_simulate(self, mode, initialCondition, time_bound, map=None):
+    def TC_simulate(self, mode, initialCondition, time_bound, lane_map:LaneMap=None):
         mode = mode.split(',')
         vehicle_mode = mode[0]
         vehicle_lane = mode[1]
@@ -27,10 +28,11 @@ class CarAgent(BaseAgent):
 
         init = initialCondition
         trace = [[0]+init]
+        lane_parameter = lane_map.lane_geometry(vehicle_lane)
         if vehicle_mode == "Normal":
             for i in range(len(t)):
                 x,y,theta,v = init
-                d = -y
+                d = -y+lane_parameter
                 psi = -theta
                 steering = psi + np.arctan2(0.45*d, v)
                 steering = np.clip(steering, -0.61, 0.61)
@@ -43,7 +45,7 @@ class CarAgent(BaseAgent):
         elif vehicle_mode == "SwitchLeft":
             for i in range(len(t)):
                 x,y,theta,v = init
-                d = -y+1
+                d = -y+3+lane_parameter
                 psi = -theta
                 steering = psi + np.arctan2(0.45*d, v)
                 steering = np.clip(steering, -0.61, 0.61)
@@ -56,7 +58,7 @@ class CarAgent(BaseAgent):
         elif vehicle_mode == "SwitchRight":
             for i in range(len(t)):
                 x,y,theta,v = init
-                d = -y-1
+                d = -y-3+lane_parameter
                 psi = -theta
                 steering = psi + np.arctan2(0.45*d, v)
                 steering = np.clip(steering, -0.61, 0.61)
diff --git a/src/example/example_map/__init__.py b/src/example/example_map/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/example/example_map/simple_map.py b/src/example/example_map/simple_map.py
new file mode 100644
index 0000000000000000000000000000000000000000..2665aa9e49a049f3aba09dccdb3ad1c9014a3b35
--- /dev/null
+++ b/src/example/example_map/simple_map.py
@@ -0,0 +1,50 @@
+from src.scene_verifier.map.lane_map import LaneMap
+from src.scene_verifier.map.lane_segment import LaneSegment
+
+class SimpleMap(LaneMap):
+    def __init__(self):
+        super().__init__()
+        segment1 = LaneSegment('Lane0', 0)
+        segment2 = LaneSegment('Lane1', 3)
+        self.add_lanes([segment1,segment2])
+        self.left_lane_dict[segment1.id].append(segment2.id)
+        self.right_lane_dict[segment2.id].append(segment1.id)
+
+class SimpleMap2(LaneMap):
+    def __init__(self):
+        super().__init__()
+        segment1 = LaneSegment('Lane0', 3)
+        segment2 = LaneSegment('Lane1', 0)
+        segment3 = LaneSegment('Lane2', -3)
+        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)
+
+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)
+    print(test_map.right_lane_dict)
+    print(test_map.lane_segment_dict)
diff --git a/src/example/example_sensor/__init__.py b/src/example/example_sensor/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/example/example_sensor/fake_sensor.py b/src/example/example_sensor/fake_sensor.py
new file mode 100644
index 0000000000000000000000000000000000000000..d415d1d555d8fab301810937ae89a227c8788f3d
--- /dev/null
+++ b/src/example/example_sensor/fake_sensor.py
@@ -0,0 +1,98 @@
+import numpy as np
+
+class FakeSensor1:
+    def sense(self, scenario, agent, state_dict, lane_map):
+        cnts = {}
+        disc = {}
+        state = state_dict['ego'][0]
+        mode = state_dict['ego'][1].split(',')
+        cnts['ego.x'] = state[1]
+        cnts['ego.y'] = state[2]
+        cnts['ego.theta'] = state[3]
+        cnts['ego.v'] = state[4]
+        disc['ego.vehicle_mode'] = mode[0]
+        disc['ego.lane_mode'] = mode[1]
+        return cnts, disc
+
+
+class FakeSensor2:
+    def sense(self, scenario, agent, state_dict, lane_map):
+        cnts = {}
+        disc = {}
+        tmp = np.array(state_dict['car1'][0])
+        if tmp.ndim < 2:
+            if agent.id == 'car1':
+                state = state_dict['car1'][0]
+                mode = state_dict['car1'][1].split(',')
+                cnts['ego.x'] = state[1]
+                cnts['ego.y'] = state[2]
+                cnts['ego.theta'] = state[3]
+                cnts['ego.v'] = state[4]
+                disc['ego.vehicle_mode'] = mode[0]
+                disc['ego.lane_mode'] = mode[1]
+
+                state = state_dict['car2'][0]
+                mode = state_dict['car2'][1].split(',')
+                cnts['other.x'] = state[1]
+                cnts['other.y'] = state[2]
+                cnts['other.theta'] = state[3]
+                cnts['other.v'] = state[4]
+                disc['other.vehicle_mode'] = mode[0]
+                disc['other.lane_mode'] = mode[1]
+            elif agent.id == 'car2':
+                state = state_dict['car2'][0]
+                mode = state_dict['car2'][1].split(',')
+                cnts['ego.x'] = state[1]
+                cnts['ego.y'] = state[2]
+                cnts['ego.theta'] = state[3]
+                cnts['ego.v'] = state[4]
+                disc['ego.vehicle_mode'] = mode[0]
+                disc['ego.lane_mode'] = mode[1]
+
+                state = state_dict['car1'][0]
+                mode = state_dict['car1'][1].split(',')
+                cnts['other.x'] = state[1]
+                cnts['other.y'] = state[2]
+                cnts['other.theta'] = state[3]
+                cnts['other.v'] = state[4]
+                disc['other.vehicle_mode'] = mode[0]
+                disc['other.lane_mode'] = mode[1]
+            return cnts, disc
+        else:
+            if agent.id == 'car1':
+                state = state_dict['car1'][0]
+                mode = state_dict['car1'][1].split(',')
+                cnts['ego.x'] = [state[0][1], state[1][1]]
+                cnts['ego.y'] = [state[0][2], state[1][2]]
+                cnts['ego.theta'] = [state[0][3], state[1][3]]
+                cnts['ego.v'] = [state[0][4], state[1][4]]
+                disc['ego.vehicle_mode'] = mode[0]
+                disc['ego.lane_mode'] = mode[1]
+
+                state = state_dict['car2'][0]
+                mode = state_dict['car2'][1].split(',')
+                cnts['other.x'] = [state[0][1], state[1][1]]
+                cnts['other.y'] = [state[0][2], state[1][2]]
+                cnts['other.theta'] = [state[0][3], state[1][3]]
+                cnts['other.v'] = [state[0][4], state[1][4]]
+                disc['other.vehicle_mode'] = mode[0]
+                disc['other.lane_mode'] = mode[1]
+            elif agent.id == 'car2':
+                state = state_dict['car2'][0]
+                mode = state_dict['car2'][1].split(',')
+                cnts['ego.x'] = [state[0][1], state[1][1]]
+                cnts['ego.y'] = [state[0][2], state[1][2]]
+                cnts['ego.theta'] = [state[0][3], state[1][3]]
+                cnts['ego.v'] = [state[0][4], state[1][4]]
+                disc['ego.vehicle_mode'] = mode[0]
+                disc['ego.lane_mode'] = mode[1]
+
+                state = state_dict['car1'][0]
+                mode = state_dict['car1'][1].split(',')
+                cnts['other.x'] = [state[0][1], state[1][1]]
+                cnts['other.y'] = [state[0][2], state[1][2]]
+                cnts['other.theta'] = [state[0][3], state[1][3]]
+                cnts['other.v'] = [state[0][4], state[1][4]]
+                disc['other.vehicle_mode'] = mode[0]
+                disc['other.lane_mode'] = mode[1]
+            return cnts, disc
diff --git a/src/plotter/__init__.py b/src/plotter/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..5dbc287d927065e02e00ed2cccd0389194764818
--- /dev/null
+++ b/src/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/src/plotter/parser.py b/src/plotter/parser.py
new file mode 100644
index 0000000000000000000000000000000000000000..e5567fbd046e05a2895aeef5a2e7a7078c74a2ca
--- /dev/null
+++ b/src/plotter/parser.py
@@ -0,0 +1,38 @@
+"""
+This file consist parser code for DryVR reachtube output
+"""
+
+from typing import TextIO
+import re 
+
+class Parser:
+    def __init__(self, f: TextIO):
+        data = f.readlines()
+        curr_key = ""
+        self.data_dict = {}
+        i = 0
+        while i < len(data):
+            line = data[i]
+            if not re.match('^[-+0-9+.+0-9+e+0-9 ]+$', line):
+                self.data_dict[line] = []
+                curr_key = line
+                i += 1
+            else:
+                line_lower = data[i]
+                line_lower_list = line_lower.split(' ')
+                line_lower_list = [float(elem) for elem in line_lower_list]
+                line_upper = data[i+1]
+                line_upper_list = line_upper.split(' ')
+                line_upper_list = [float(elem) for elem in line_upper_list]
+                rect = [line_lower_list, line_upper_list]
+                self.data_dict[curr_key].append(rect)
+                i += 2
+                
+    
+    def get_all_data(self):
+        res = []
+        for key in self.data_dict:
+            res += self.data_dict[key]
+        return res
+
+    
\ No newline at end of file
diff --git a/src/plotter/plotter2D.py b/src/plotter/plotter2D.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c79b5fc7a7f6571ec93b6f517567130c1aa44c3
--- /dev/null
+++ b/src/plotter/plotter2D.py
@@ -0,0 +1,59 @@
+"""
+This file consist main plotter code for DryVR reachtube output
+"""
+
+import matplotlib.patches as patches
+import matplotlib.pyplot as plt
+import numpy as np 
+from typing import List 
+
+colors = ['red', 'green', 'blue', 'yellow', 'black']
+
+def plot(
+    data, 
+    x_dim: int = 0, 
+    y_dim_list: List[int] = [1], 
+    color = 'b', 
+    fig = None, 
+    x_lim = (float('inf'), -float('inf')), 
+    y_lim = (float('inf'), -float('inf'))
+):
+    if fig is None:
+        fig = plt.figure()
+    ax = plt.gca()
+    x_min, x_max = x_lim
+    y_min, y_max = y_lim
+    for rect in data:
+        lb = rect[0]
+        ub = rect[1]
+        for y_dim in y_dim_list:
+            rect_patch = patches.Rectangle((lb[x_dim], lb[y_dim]), ub[x_dim]-lb[x_dim], ub[y_dim]-lb[y_dim], color = color)
+            ax.add_patch(rect_patch)
+            x_min = min(lb[x_dim], x_min)
+            y_min = min(lb[y_dim], y_min)
+            x_max = max(ub[x_dim], x_max)
+            y_max = max(ub[y_dim], y_max)
+
+    ax.set_xlim([x_min-1, x_max+1])
+    ax.set_ylim([y_min-1, y_max+1])
+    return fig, (x_min, x_max), (y_min, y_max)
+
+def plot_tree(root, agent_id, x_dim: int=0, y_dim_list: List[int]=[1], color='b', fig = None):
+    if fig is None:
+        fig = plt.figure()
+    
+    queue = [root]
+    x_lim = (float('inf'), -float('inf')) 
+    y_lim = (float('inf'), -float('inf'))
+    while queue != []:
+        node = queue.pop(0)
+        traces = node.trace
+        trace = traces[agent_id]
+        data = []
+        for i in range(0,len(trace),2):
+            data.append([trace[i], trace[i+1]])
+        fig, x_lim, y_lim = plot(data, x_dim, y_dim_list, color, fig, x_lim, y_lim)
+
+        queue += node.child
+
+    return fig
\ No newline at end of file
diff --git a/src/plotter/plotter3D.py b/src/plotter/plotter3D.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ede38e0f2e077b92f048da1d9a9c88c9dedbdb2
--- /dev/null
+++ b/src/plotter/plotter3D.py
@@ -0,0 +1,68 @@
+# Plot polytope in 3d
+# Written by: Kristina Miller
+
+import numpy as np
+import matplotlib.pyplot as plt
+import mpl_toolkits.mplot3d as a3
+
+from scipy.spatial import ConvexHull
+import polytope as pc
+import pyvista as pv
+
+def plot3d(node, x_dim, y_dim, z_dim, ax = None):
+    if ax is None:
+        ax = pv.Plotter()
+    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]))
+
+        box = [[lb[x_dim], lb[y_dim], lb[z_dim]],[ub[x_dim], ub[y_dim], ub[z_dim]]]
+        poly = pc.box2poly(np.array(box).T)
+        plot_polytope_3d(poly.A, poly.b, ax = ax, color = '#b3de69')
+
+def plot_polytope_3d(A, b, ax = None, color = 'red', trans = 0.2, edge = True):
+	if ax is None:
+		ax = pv.Plotter() 
+
+	poly = pc.Polytope(A = A, b = b)
+	vertices = pc.extreme(poly)
+	cloud = pv.PolyData(vertices)
+	volume = cloud.delaunay_3d()
+	shell = volume.extract_geometry()
+	ax.add_mesh(shell, opacity=trans, color=color)
+	if edge:
+		edges = shell.extract_feature_edges(20)
+		ax.add_mesh(edges, color="k", line_width=1)
+
+def plot_line_3d(start, end, ax = None, color = 'blue', line_width = 1):
+	if ax is None:
+		ax = pv.Plotter() 
+
+	a = start
+	b = end
+
+	# Preview how this line intersects this mesh
+	line = pv.Line(a, b)
+
+	ax.add_mesh(line, color=color, line_width=line_width)
+
+if __name__ == '__main__':
+	A = np.array([[-1, 0, 0],
+				  [1, 0, 0],
+				  [0, -1, 0],
+				  [0, 1, 0],
+				  [0, 0, -1],
+				  [0, 0, 1]])
+	b = np.array([[1], [1], [1], [1], [1], [1]])
+	b2 = np.array([[-1], [2], [-1], [2], [-1], [2]])
+	ax1 = a3.Axes3D(plt.figure())
+	plot_polytope_3d(A, b, ax = ax1, color = 'red')
+	plot_polytope_3d(A, b2, ax = ax1, color = 'green')
+	plt.show()
\ No newline at end of file
diff --git a/src/scene_verifier/__init__.py b/src/scene_verifier/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/scene_verifier/agents/__init__.py b/src/scene_verifier/agents/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ourtool/agents/base_agent.py b/src/scene_verifier/agents/base_agent.py
similarity index 78%
rename from ourtool/agents/base_agent.py
rename to src/scene_verifier/agents/base_agent.py
index a533b3579b3e110806a09538e43b24b16f1ca2cc..e0044528e0c7ce5c0d0119b932d049c8d65b4802 100644
--- a/ourtool/agents/base_agent.py
+++ b/src/scene_verifier/agents/base_agent.py
@@ -1,4 +1,4 @@
-from pythonparser import ControllerAst
+from src.scene_verifier.code_parser.pythonparser import ControllerAst
 
 class BaseAgent:
     def __init__(self, id, code = None, file_name = None):  
diff --git a/src/scene_verifier/analysis/__init__.py b/src/scene_verifier/analysis/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/scene_verifier/analysis/analysis_tree_node.py b/src/scene_verifier/analysis/analysis_tree_node.py
new file mode 100644
index 0000000000000000000000000000000000000000..8f771685883c032ad9e06074cb869f92baaf58a3
--- /dev/null
+++ b/src/scene_verifier/analysis/analysis_tree_node.py
@@ -0,0 +1,26 @@
+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,
+        ndigits = 10
+    ):
+        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 = round(start_time,ndigits)
diff --git a/ourtool/simulator/simulator.py b/src/scene_verifier/analysis/simulator.py
similarity index 80%
rename from ourtool/simulator/simulator.py
rename to src/scene_verifier/analysis/simulator.py
index 0282a49e1f70541235f7d81b32569c21beb47e12..30f47ce5dd44ab46dc9f9bdfca34fc528fdcbcab 100644
--- a/ourtool/simulator/simulator.py
+++ b/src/scene_verifier/analysis/simulator.py
@@ -1,39 +1,18 @@
 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 src.scene_verifier.agents.base_agent import BaseAgent
+from src.scene_verifier.analysis.analysis_tree_node import AnalysisTreeNode
 
 class Simulator:
     def __init__(self):
         self.simulation_tree_root = None
 
-    def simulate(self, init_list, init_mode_list, agent_list:List[BaseAgent], transition_graph, time_horizon):
+    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:
@@ -57,7 +36,7 @@ class Simulator:
                     # Simulate 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)
+                    trace = node.agent[agent_id].TC_simulate(mode, init, remain_time,lane_map)
                     trace[:,0] += node.start_time
                     node.trace[agent_id] = trace.tolist()
 
@@ -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/src/scene_verifier/analysis/verifier.py b/src/scene_verifier/analysis/verifier.py
new file mode 100644
index 0000000000000000000000000000000000000000..53a3b501ba64c8b63d3feb33c1e5d4be0cbbf31b
--- /dev/null
+++ b/src/scene_verifier/analysis/verifier.py
@@ -0,0 +1,109 @@
+from typing import List, Dict
+import copy
+
+import numpy as np
+
+from src.scene_verifier.agents.base_agent import BaseAgent
+from src.scene_verifier.analysis.analysis_tree_node import AnalysisTreeNode
+from src.scene_verifier.dryvr.core.dryvrcore import calc_bloated_tube
+import src.scene_verifier.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
+            # TODO: can add parallalization for this loop
+            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")
+            
+            # Check safety conditions here
+
+            # Get all possible transitions to next mode
+            all_possible_transitions = transition_graph.get_all_transition_set(node)
+            max_end_idx = 0
+            for transition in all_possible_transitions:
+                transit_agent_idx, src_mode, dest_mode, next_init, idx = transition 
+                start_idx, end_idx = idx
+ 
+                truncated_trace = {}
+                for agent_idx in node.agent:
+                    truncated_trace[agent_idx] = node.trace[agent_idx][start_idx*2:]
+                if end_idx > max_end_idx:
+                    max_end_idx = end_idx
+                next_node_mode = copy.deepcopy(node.mode) 
+                next_node_mode[transit_agent_idx] = dest_mode 
+                next_node_agent = node.agent 
+                next_node_start_time = list(truncated_trace.values())[0][0][0]
+                next_node_init = {}
+                next_node_trace = {}
+                for agent_idx in next_node_agent:
+                    if agent_idx == transit_agent_idx:
+                        next_node_init[agent_idx] = next_init 
+                    else:
+                        next_node_trace[agent_idx] = truncated_trace[agent_idx]
+                
+                tmp = AnalysisTreeNode(
+                    trace = next_node_trace,
+                    init = next_node_init,
+                    mode = next_node_mode,
+                    agent = next_node_agent,
+                    child = [],
+                    start_time = next_node_start_time
+                )
+                node.child.append(tmp)
+                verification_queue.append(tmp)
+
+            """Truncate trace of current node based on max_end_idx"""
+            """Only truncate when there's transitions"""
+            if all_possible_transitions:
+                for agent_idx in node.agent:
+                    node.trace[agent_idx] = node.trace[agent_idx][:(max_end_idx+1)*2]
+        return root
\ No newline at end of file
diff --git a/src/scene_verifier/automaton/__init__.py b/src/scene_verifier/automaton/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..a1ae96f571beeb44895a0b155022d6c852604e6e
--- /dev/null
+++ b/src/scene_verifier/automaton/__init__.py
@@ -0,0 +1,4 @@
+# try:
+#     from hybrid_automaton import *
+# except:
+#     from scene_verifier.automaton.hybrid_automaton import *
\ No newline at end of file
diff --git a/src/scene_verifier/automaton/guard.py b/src/scene_verifier/automaton/guard.py
new file mode 100644
index 0000000000000000000000000000000000000000..3c65044756030181162ff058f0198b68070bc1fe
--- /dev/null
+++ b/src/scene_verifier/automaton/guard.py
@@ -0,0 +1,441 @@
+import enum
+import re
+from typing import List, Dict
+import pickle
+# from scene_verifier.automaton.hybrid_io_automaton import HybridIoAutomaton
+# from pythonparser import Guard
+import ast
+import copy
+
+from z3 import *
+import sympy
+
+import astunparse
+
+class LogicTreeNode:
+    def __init__(self, data, child = [], val = None, mode_guard = None):
+        self.data = data 
+        self.child = child
+        self.val = val
+        self.mode_guard = mode_guard
+
+class GuardExpressionAst:
+    def __init__(self, guard_list):
+        self.ast_list = []
+        for guard in guard_list:
+            self.ast_list.append(copy.deepcopy(guard.ast))
+        self.cont_variables = {}
+        self.varDict = {'t':Real('t')}
+
+    def _build_guard(self, guard_str, agent):
+        """
+        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("==", ">=")
+        for vars in self.cont_variables:
+            sympy_guard_str = sympy_guard_str.replace(vars, self.cont_variables[vars])
+
+        symbols = list(sympy.sympify(sympy_guard_str, evaluate=False).free_symbols)
+        symbols = [str(s) for s in symbols]
+        tmp = list(self.cont_variables.values())
+        symbols_map = {}
+        for s in symbols:
+            if s in tmp:
+                key = list(self.cont_variables.keys())[list(self.cont_variables.values()).index(s)]
+                symbols_map[s] = key
+
+        for vars in reversed(self.cont_variables):
+            guard_str = guard_str.replace(vars, self.cont_variables[vars])
+        guard_str = self._handleReplace(guard_str)
+        cur_solver.add(eval(guard_str))  # TODO use an object instead of `eval` a string
+        return cur_solver, symbols_map
+
+    def _handleReplace(self, input_str):
+        """
+        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 = list(self.varDict.keys())
+
+        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.varDict["' + key + '"]'
+            input_str = input_str[:idx[0]] + target + input_str[idx[1]:]
+        return input_str
+
+    def evaluate_guard_cont(self, agent, continuous_variable_dict, lane_map):
+        res = False
+        is_contained = False
+
+        for cont_vars in continuous_variable_dict:
+            self.cont_variables[cont_vars] = cont_vars.replace('.','_')
+            self.varDict[cont_vars.replace('.','_')] = Real(cont_vars.replace('.','_'))
+
+        z3_string = self.generate_z3_expression() 
+        if isinstance(z3_string, bool):
+            if z3_string:
+                return True, True 
+            else:
+                return False, False
+
+        cur_solver, symbols = self._build_guard(z3_string, agent)
+        cur_solver.push()
+        for symbol in symbols:
+            cur_solver.add(self.varDict[symbol] >= continuous_variable_dict[symbols[symbol]][0])
+            cur_solver.add(self.varDict[symbol] <= continuous_variable_dict[symbols[symbol]][1])
+        if cur_solver.check() == sat:
+            # The reachtube hits the guard
+            cur_solver.pop()
+            res = True
+            
+            # 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.varDict[symbol] >= continuous_variable_dict[symbols[symbol]][0])
+                tmp_solver.add(self.varDict[symbol] <= continuous_variable_dict[symbols[symbol]][1])
+            if tmp_solver.check() == unsat:
+                print("Full intersect, break")
+                is_contained = True
+
+        return res, is_contained
+
+    def generate_z3_expression(self):
+        """
+        The return value of this function will be a bool/str
+
+        If without evaluating the continuous variables the result is True, then
+        the guard will automatically be satisfied and is_contained will be True
+
+        If without evaluating the continuous variables the result is False, th-
+        en the guard will automatically be unsatisfied
+
+        If the result is a string, then continuous variables will be checked to
+        see if the guard can be satisfied 
+        """
+        res = []
+        for node in self.ast_list:
+            tmp = self._generate_z3_expression_node(node)
+            if isinstance(tmp, bool):
+                if not tmp:
+                    return False
+                else:
+                    continue
+            res.append(tmp)
+        if res == []:
+            return True
+        elif len(res) == 1:
+            return res[0]
+        res = "And("+",".join(res)+")"
+        return res
+
+    def _generate_z3_expression_node(self, node):
+        """
+        Perform a DFS over expression ast and generate the guard expression
+        The return value of this function can be a bool/str
+
+        If without evaluating the continuous variables the result is True, then
+        the guard condition will automatically be satisfied
+        
+        If without evaluating the continuous variables the result is False, then
+        the guard condition will not be satisfied
+
+        If the result is a string, then continuous variables will be checked to
+        see if the guard can be satisfied
+        """
+        if isinstance(node, ast.BoolOp):
+            # Check the operator
+            # For each value in the boolop, check results
+            if isinstance(node.op, ast.And):
+                z3_str = []
+                for i,val in enumerate(node.values):
+                    tmp = self._generate_z3_expression_node(val)
+                    if isinstance(tmp, bool):
+                        if tmp:
+                            continue 
+                        else:
+                            return False
+                    z3_str.append(tmp)
+                z3_str = 'And('+','.join(z3_str)+')'
+                return z3_str
+            elif isinstance(node.op, ast.Or):
+                z3_str = []
+                for val in node.values:
+                    tmp = self._generate_z3_expression_node(val)
+                    if isinstance(tmp, bool):
+                        if tmp:
+                            return True
+                        else:
+                            continue
+                    z3_str.append(tmp)
+                z3_str = 'Or('+','.join(z3_str)+')'
+                return z3_str
+            # If string, construct string
+            # If bool, check result and discard/evaluate result according to operator
+            pass 
+        elif isinstance(node, ast.Constant):
+            # If is bool, return boolean result
+            if isinstance(node.value, bool):
+                return node.value
+            # Else, return raw expression
+            else:
+                expr = astunparse.unparse(node)
+                expr = expr.strip('\n')
+                return expr
+        elif isinstance(node, ast.UnaryOp):
+            # If is UnaryOp, 
+            value = self._generate_z3_expression_node(node.operand)
+            if isinstance(node.op, ast.USub):
+                return -value
+        else:
+            # For other cases, we can return the expression directly
+            expr = astunparse.unparse(node)
+            expr = expr.strip('\n')
+            return expr
+
+    def evaluate_guard_disc(self, agent, discrete_variable_dict, lane_map):
+        """
+        Evaluate guard that involves only discrete variables. 
+        """
+        res = True
+        for i, node in enumerate(self.ast_list):
+            tmp, self.ast_list[i] = self._evaluate_guard_disc(node, agent, discrete_variable_dict, lane_map)
+            res = res and tmp 
+        return res
+            
+    def _evaluate_guard_disc(self, root, agent, disc_var_dict, lane_map):
+        if isinstance(root, ast.Compare):
+            expr = astunparse.unparse(root)
+            if any([var in expr for var in disc_var_dict]):
+                left, root.left = self._evaluate_guard_disc(root.left, agent, disc_var_dict, lane_map)
+                right, root.comparators[0] = self._evaluate_guard_disc(root.comparators[0], agent, disc_var_dict, lane_map)
+                if isinstance(left, bool) or isinstance(right, bool):
+                    return True, root
+                if isinstance(root.ops[0], ast.GtE):
+                    res = left>=right
+                elif isinstance(root.ops[0], ast.Gt):
+                    res = left>right 
+                elif isinstance(root.ops[0], ast.Lt):
+                    res = left<right
+                elif isinstance(root.ops[0], ast.LtE):
+                    res = left<=right
+                elif isinstance(root.ops[0], ast.Eq):
+                    res = left == right 
+                elif isinstance(root.ops[0], ast.NotEq):
+                    res = left != right 
+                else:
+                    raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+                if res:
+                    root = ast.parse('True').body[0].value
+                else:
+                    root = ast.parse('False').body[0].value    
+                return res, root
+            else:
+                return True, root
+        elif isinstance(root, ast.BoolOp):
+            if isinstance(root.op, ast.And):
+                res = True
+                for i,val in enumerate(root.values):
+                    tmp,root.values[i] = self._evaluate_guard_disc(val, agent, disc_var_dict, lane_map)
+                    res = res and tmp
+                    if not res:
+                        break
+                return res, root
+            elif isinstance(root.op, ast.Or):
+                res = False
+                for val in root.values:
+                    tmp,val = self._evaluate_guard_disc(val, agent, disc_var_dict, lane_map)
+                    res = res or tmp
+                    if res:
+                        break
+                return res, root     
+        elif isinstance(root, ast.BinOp):
+            # Check left and right in the binop and replace all attributes involving discrete variables
+            left, root.left = self._evaluate_guard_disc(root.left, agent, disc_var_dict, lane_map)
+            right, root.right = self._evaluate_guard_disc(root.right, agent, disc_var_dict, lane_map)
+            return True, root
+        elif isinstance(root, ast.Call):
+            expr = astunparse.unparse(root)
+            # Check if the root is a function
+            if any([var in expr for var in disc_var_dict]):
+                # tmp = re.split('\(|\)',expr)
+                # while "" in tmp:
+                #     tmp.remove("")
+                # for arg in tmp[1:]:
+                #     if arg in disc_var_dict:
+                #         expr = expr.replace(arg,f'"{disc_var_dict[arg]}"')
+                # res = eval(expr)
+                for arg in disc_var_dict:
+                    expr = expr.replace(arg, f'"{disc_var_dict[arg]}"')
+                res = eval(expr)
+                if isinstance(res, bool):
+                    if res:
+                        root = ast.parse('True').body[0].value
+                    else:
+                        root = ast.parse('False').body[0].value    
+                else:
+                    root = ast.parse(str(res)).body[0].value
+                return res, root
+            else:
+                return True, root
+        elif isinstance(root, ast.Attribute):
+            expr = astunparse.unparse(root)
+            expr = expr.strip('\n')
+            if expr in disc_var_dict:
+                val = disc_var_dict[expr]
+                for mode_name in agent.controller.modes:
+                    if val in agent.controller.modes[mode_name]:
+                        val = mode_name+'.'+val
+                        break
+                return val, root
+            elif root.value.id in agent.controller.modes:
+                return expr, root
+            else:
+                return True, root
+        elif isinstance(root, ast.Constant):
+            return root.value, root
+        elif isinstance(root, ast.UnaryOp):
+            if isinstance(root.op, ast.USub):
+                res, root.operand = self._evaluate_guard_disc(root.operand, agent, disc_var_dict, lane_map)
+            else:
+                raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+            return True, root
+        else:
+            raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+
+    def evaluate_guard(self, agent, continuous_variable_dict, discrete_variable_dict, lane_map):
+        res = True
+        for node in self.ast_list:
+            tmp = self._evaluate_guard(node, agent, continuous_variable_dict, discrete_variable_dict, lane_map)
+            res = tmp and res
+            if not res:
+                break
+        return res
+
+    def _evaluate_guard(self, root, agent, cnts_var_dict, disc_var_dict, lane_map):
+        if isinstance(root, ast.Compare):
+            left = self._evaluate_guard(root.left, agent, cnts_var_dict, disc_var_dict, lane_map)
+            right = self._evaluate_guard(root.comparators[0], agent, cnts_var_dict, disc_var_dict, lane_map)
+            if isinstance(root.ops[0], ast.GtE):
+                return left>=right
+            elif isinstance(root.ops[0], ast.Gt):
+                return left>right 
+            elif isinstance(root.ops[0], ast.Lt):
+                return left<right
+            elif isinstance(root.ops[0], ast.LtE):
+                return left<=right
+            elif isinstance(root.ops[0], ast.Eq):
+                return left == right 
+            elif isinstance(root.ops[0], ast.NotEq):
+                return left != right 
+            else:
+                raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+
+        elif isinstance(root, ast.BoolOp):
+            if isinstance(root.op, ast.And):
+                res = True
+                for val in root.values:
+                    tmp = self._evaluate_guard(val, agent, cnts_var_dict, disc_var_dict, lane_map)
+                    res = res and tmp
+                    if not res:
+                        break
+                return res
+            elif isinstance(root.op, ast.Or):
+                res = False
+                for val in root.values:
+                    tmp = self._evaluate_guard(val, agent, cnts_var_dict, disc_var_dict, lane_map)
+                    res = res or tmp
+                    if res:
+                        break
+                return res
+        elif isinstance(root, ast.BinOp):
+            left = self._evaluate_guard(root.left, agent, cnts_var_dict, disc_var_dict, lane_map)
+            right = self._evaluate_guard(root.right, agent, cnts_var_dict, disc_var_dict, lane_map)
+            if isinstance(root.op, ast.Sub):
+                return left - right
+            elif isinstance(root.op, ast.Add):
+                return left + right
+            else:
+                raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+        elif isinstance(root, ast.Call):
+            expr = astunparse.unparse(root)
+            # Check if the root is a function
+            if 'map' in expr:
+                # tmp = re.split('\(|\)',expr)
+                # while "" in tmp:
+                #     tmp.remove("")
+                # for arg in tmp[1:]:
+                #     if arg in disc_var_dict:
+                #         expr = expr.replace(arg,f'"{disc_var_dict[arg]}"')
+                # res = eval(expr)
+                for arg in disc_var_dict:
+                    expr = expr.replace(arg, f'"{disc_var_dict[arg]}"')
+                for arg in cnts_var_dict:
+                    expr = expr.replace(arg, str(cnts_var_dict[arg]))    
+                res = eval(expr)
+                return res
+        elif isinstance(root, ast.Attribute):
+            expr = astunparse.unparse(root)
+            expr = expr.strip('\n')
+            if expr in disc_var_dict:
+                val = disc_var_dict[expr]
+                for mode_name in agent.controller.modes:
+                    if val in agent.controller.modes[mode_name]:
+                        val = mode_name+'.'+val
+                        break
+                return val
+            elif expr in cnts_var_dict:
+                val = cnts_var_dict[expr]
+                return val
+            elif root.value.id in agent.controller.modes:
+                return expr
+        elif isinstance(root, ast.Constant):
+            return root.value
+        else:
+            raise ValueError(f'Node type {root} from {astunparse.unparse(root)} is not supported')
+
+if __name__ == "__main__":
+    with open('tmp.pickle','rb') as f:
+        guard_list = pickle.load(f)
+    tmp = GuardExpressionAst(guard_list)
+    # tmp.evaluate_guard()
+    # tmp.construct_tree_from_str('(other_x-ego_x<20) and other_x-ego_x>10 and other_vehicle_lane==ego_vehicle_lane')
+    print("stop")
\ No newline at end of file
diff --git a/ourtool/automaton/hybrid_automaton.py b/src/scene_verifier/automaton/hybrid_automaton.py
similarity index 100%
rename from ourtool/automaton/hybrid_automaton.py
rename to src/scene_verifier/automaton/hybrid_automaton.py
diff --git a/ourtool/automaton/hybrid_io_automaton.py b/src/scene_verifier/automaton/hybrid_io_automaton.py
similarity index 90%
rename from ourtool/automaton/hybrid_io_automaton.py
rename to src/scene_verifier/automaton/hybrid_io_automaton.py
index e0b936c5f24dd19a26b728172d6be13793702b4c..74f2c963c7ec16095d2ecbc846da58a21e336702 100644
--- a/ourtool/automaton/hybrid_io_automaton.py
+++ b/src/scene_verifier/automaton/hybrid_io_automaton.py
@@ -1,4 +1,4 @@
-from ourtool.automaton.hybrid_automaton import HybridAutomaton
+from src.scene_verifier.automaton.hybrid_automaton import HybridAutomaton
 
 class HybridIoAutomaton(HybridAutomaton):
     def __init__(
diff --git a/src/scene_verifier/automaton/reset.py b/src/scene_verifier/automaton/reset.py
new file mode 100644
index 0000000000000000000000000000000000000000..a0e9ddcb5bdeba9297aebfa86d29f929b2e76ccf
--- /dev/null
+++ b/src/scene_verifier/automaton/reset.py
@@ -0,0 +1,84 @@
+import itertools
+
+import numpy as np 
+
+class ResetExpression:
+    def __init__(self, reset_list):
+        self.ast_list = []
+        for reset in reset_list:
+            self.ast_list.append(reset.ast)
+        self.expr_list = []
+        for reset in reset_list:
+            self.expr_list.append(reset.code)
+
+    def apply_reset_continuous(self, agent, continuous_variable_dict, lane_map):
+        agent_state_lower = []
+        agent_state_upper = []
+        for var in agent.controller.vars_dict['ego']['cont']:
+            agent_state_lower.append(continuous_variable_dict['ego.'+var][0])
+            agent_state_upper.append(continuous_variable_dict['ego.'+var][1])
+        assert len(agent_state_lower) == len(agent_state_upper) == len(agent.controller.vars_dict['ego']['cont'])       
+        for expr in self.expr_list:
+            if 'mode' not in expr:
+                tmp = expr.split('=')
+                lhs, rhs = tmp[0], tmp[1]
+                for lhs_idx, cts_variable in enumerate(agent.controller.vars_dict['ego']['cont']):
+                    if "output."+cts_variable == lhs:
+                        break 
+
+                lower = float('inf')
+                upper = -float('inf')
+
+                symbols = []
+                for var in continuous_variable_dict:
+                    if var in expr:
+                        symbols.append(var)
+                
+                combinations = self._get_combinations(symbols, continuous_variable_dict)
+                # for cts_variable in continuous_variable_dict:
+                #     tmp = tmp.replace(cts_variable, str(continuous_variable_dict[cts_variable]))
+                # next_init[i] = eval(tmp)
+                for i in combinations.shape[0]:
+                    comb = combinations[i,:]
+                    for j in range(len(symbols)):
+                        tmp = rhs.replace(symbols[j], str(comb[i,j]))
+                        tmp = min(tmp, lower)
+                        tmp = max(tmp, upper)
+
+                agent_state_lower[lhs_idx] = lower 
+                agent_state_upper[lhs_idx] = upper
+        
+        return [agent_state_lower, agent_state_upper]
+
+    def _get_combinations(self, symbols, cont_var_dict):
+        all_vars = []
+        for symbol in symbols:
+            all_vars.append(cont_var_dict[symbol])
+        comb_array = np.array(np.meshgrid(*all_vars)).T.reshape(-1, len(symbols))    
+        return comb_array 
+
+    def get_dest(self, agent, agent_state, discrete_variable_dict, lane_map) -> str:
+        agent_mode = agent_state[1]
+        dest = agent_mode.split(',')
+        possible_dest = [[elem] for elem in dest]
+        for reset in self.expr_list:
+            if "mode" in reset:
+                for i, discrete_variable_ego in enumerate(agent.controller.vars_dict['ego']['disc']):
+                    if discrete_variable_ego in reset:
+                        break
+                tmp = reset.split('=')
+                if 'map' in tmp[1]:
+                    tmp = tmp[1]
+                    for var in discrete_variable_dict:
+                        tmp = tmp.replace(var, f"'{discrete_variable_dict[var]}'")
+                    possible_dest[i] = eval(tmp)
+                else:
+                    tmp = tmp[1].split('.')
+                    if tmp[0].strip(' ') in agent.controller.modes:
+                        possible_dest[i] = [tmp[1]]                            
+        all_dest = itertools.product(*possible_dest)
+        res = []
+        for dest in all_dest:
+            dest = ','.join(dest)
+            res.append(dest)
+        return res
\ No newline at end of file
diff --git a/src/scene_verifier/code_parser/__init__.py b/src/scene_verifier/code_parser/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/pythonparser.py b/src/scene_verifier/code_parser/pythonparser.py
similarity index 94%
rename from pythonparser.py
rename to src/scene_verifier/code_parser/pythonparser.py
index 02d940e6931c8e58b0792aef8875b71da750819c..885ab535ec24cd03ca8ae971655b861be327a7c5 100644
--- a/pythonparser.py
+++ b/src/scene_verifier/code_parser/pythonparser.py
@@ -28,11 +28,12 @@ Statement super class. Holds the code and mode information for a statement.
 If there is no mode information, mode and modeType are None.
 '''
 class Statement:
-    def __init__(self, code, mode, modeType):
+    def __init__(self, code, mode, modeType, func = None, args = None):
         self.code = code
         self.modeType = modeType
         self.mode = mode
-
+        self.func = func 
+        self.args = args
     
     def print(self):
         print(self.code)
@@ -42,8 +43,9 @@ class Statement:
 Guard class. Subclass of statement.
 '''
 class Guard(Statement):
-    def __init__(self, code, mode, modeType):
-        super().__init__(code, mode, modeType)
+    def __init__(self, code, mode, modeType, inp_ast, func=None, args=None):
+        super().__init__(code, mode, modeType, func, args)
+        self.ast = inp_ast
 
 
     '''
@@ -64,18 +66,27 @@ class Guard(Statement):
                 if ("Mode" in str(node.test.comparators[0].value.id)):
                     modeType = str(node.test.comparators[0].value.id)
                     mode = str(node.test.comparators[0].attr)
-                    return Guard(ast.get_source_segment(code, node.test), mode, modeType)
+                    return Guard(ast.get_source_segment(code, node.test), mode, modeType, node.test)
             else:
-                return Guard(ast.get_source_segment(code, node.test), None, None)
-        else:
-            return Guard(ast.get_source_segment(code, node.test), None, None)
-        
+                return Guard(ast.get_source_segment(code, node.test), None, None, node.test)
+        elif isinstance(node.test, ast.BoolOp):
+            return Guard(ast.get_source_segment(code, node.test), None, None, node.test)
+        elif isinstance(node.test, ast.Call):
+            source_segment = ast.get_source_segment(code, node.test)
+            if "map" in source_segment:
+                func = node.test.func.value.id + '.' + node.test.func.attr 
+                args = []
+                for arg in node.test.args:
+                    args.append(arg.value.id + '.' + arg.attr)
+                return Guard(source_segment, None, None, node.test, func, args)
+
 '''
 Reset class. Subclass of statement.
 '''
 class Reset(Statement):
-    def __init__(self, code, mode, modeType):
+    def __init__(self, code, mode, modeType, inp_ast):
         super().__init__(code, mode, modeType)
+        self.ast = inp_ast
 
     '''
     Returns true if a reset is updating our mode. 
@@ -95,8 +106,8 @@ class Reset(Statement):
             if ("Mode" in str(node.value.value.id)):
                 modeType = str(node.value.value.id)
                 mode = str(node.value.attr)
-            return Reset(ast.get_source_segment(code, node), mode, modeType)
-        return Reset(ast.get_source_segment(code, node), None, None)
+            return Reset(ast.get_source_segment(code, node), mode, modeType, node)
+        return Reset(ast.get_source_segment(code, node), None, None, node)
 
 
 '''
@@ -461,9 +472,9 @@ if __name__ == "__main__":
     #    print("incorrect usage. call createGraph.py program inputfile outputfilename")
     #    quit()
 
-    input_code_name = 'toythermomini.py' #sys.argv[1]
-    input_file_name = 'billiard_input.json' #sys.argv[2] 
-    output_file_name = 'out.json' #sys.argv[3]
+    input_code_name = sys.argv[1]
+    input_file_name = sys.argv[2]
+    output_file_name = sys.argv[3]
 
     with open(input_file_name) as in_json_file:
         input_json = json.load(in_json_file)
diff --git a/src/scene_verifier/dryvr/__init__.py b/src/scene_verifier/dryvr/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/scene_verifier/dryvr/common/__init__.py b/src/scene_verifier/dryvr/common/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..7dbea944864d1e4d5337014013f9b49ae5ce1bb1
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/common/config.py b/src/scene_verifier/dryvr/common/config.py
new file mode 100644
index 0000000000000000000000000000000000000000..6b23f518d74bd9d178fe2f56ad613ad65c0ea1a5
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/common/constant.py b/src/scene_verifier/dryvr/common/constant.py
new file mode 100644
index 0000000000000000000000000000000000000000..0eaaa69353ee78fe8b7023ec6d73a12550213335
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/common/io.py b/src/scene_verifier/dryvr/common/io.py
new file mode 100644
index 0000000000000000000000000000000000000000..bf864448b0ab522a1f5e17df03d67e2be0512390
--- /dev/null
+++ b/src/scene_verifier/dryvr/common/io.py
@@ -0,0 +1,156 @@
+"""
+This file contains IO functions for DryVR
+"""
+
+import six
+
+from src.scene_verifier.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 example 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 example 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/src/scene_verifier/dryvr/common/utils.py b/src/scene_verifier/dryvr/common/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..1923a84ea51badf888fc64852a0d89edd45d3575
--- /dev/null
+++ b/src/scene_verifier/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 example config to config module
+    
+    Args:
+        configObj (module): config module
+        userConfig (dict): example 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/src/scene_verifier/dryvr/core/__init__.py b/src/scene_verifier/dryvr/core/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/scene_verifier/dryvr/core/distance.py b/src/scene_verifier/dryvr/core/distance.py
new file mode 100644
index 0000000000000000000000000000000000000000..8da15a12091e3501880ca1411bb5fc4e4156099d
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/dryvrcore.py b/src/scene_verifier/dryvr/core/dryvrcore.py
new file mode 100644
index 0000000000000000000000000000000000000000..14df5292f4f25786be7c265288aaa585663d15c4
--- /dev/null
+++ b/src/scene_verifier/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 src.scene_verifier.dryvr.common.constant import *
+from src.scene_verifier.dryvr.common.io import writeReachTubeFile
+from src.scene_verifier.dryvr.common.utils import randomPoint, calcDelta, calcCenterPoint, trimTraces
+from src.scene_verifier.dryvr.discrepancy.Global_Disc import get_reachtube_segment
+# from scene_verifier.dryvr.tube_computer.backend.reachabilityengine import ReachabilityEngine
+# from scene_verifier.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/src/scene_verifier/dryvr/core/dryvrmain.py b/src/scene_verifier/dryvr/core/dryvrmain.py
new file mode 100644
index 0000000000000000000000000000000000000000..7b085aab083f01e9af8b9fe52bd82857493d544c
--- /dev/null
+++ b/src/scene_verifier/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): example-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 example,
+    # If example 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): example-specified configuration
+
+    Returns:
+        None
+
+    """
+    if param_config is None:
+        param_config = {}
+    # There are some fields can be config by example,
+    # If example 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/src/scene_verifier/dryvr/core/goalchecker.py b/src/scene_verifier/dryvr/core/goalchecker.py
new file mode 100644
index 0000000000000000000000000000000000000000..ba84296c159c2d102123179eb75a5cb640e71dcc
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/graph.py b/src/scene_verifier/dryvr/core/graph.py
new file mode 100644
index 0000000000000000000000000000000000000000..d4770ff2001d436b7c4e0c9408812cea809c2fa1
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/guard.py b/src/scene_verifier/dryvr/core/guard.py
new file mode 100644
index 0000000000000000000000000000000000000000..76e7e8238defec4aaa34af58a41ce1609455306f
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/initialset.py b/src/scene_verifier/dryvr/core/initialset.py
new file mode 100644
index 0000000000000000000000000000000000000000..0bc2a0a4f00648f93e859d83b6d8c0e172b5d69b
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/initialsetstack.py b/src/scene_verifier/dryvr/core/initialsetstack.py
new file mode 100644
index 0000000000000000000000000000000000000000..e3474a07eb00a97e1fe734eaccc323ea410d2500
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/reachtube.py b/src/scene_verifier/dryvr/core/reachtube.py
new file mode 100644
index 0000000000000000000000000000000000000000..850b8a5a1fd796541d4cade354e0ad3537b3ad6c
--- /dev/null
+++ b/src/scene_verifier/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 example 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/src/scene_verifier/dryvr/core/reset.py b/src/scene_verifier/dryvr/core/reset.py
new file mode 100644
index 0000000000000000000000000000000000000000..684b9c5e2b3af8de2ba8231d51e8d3dafdd32875
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/core/uniformchecker.py b/src/scene_verifier/dryvr/core/uniformchecker.py
new file mode 100644
index 0000000000000000000000000000000000000000..57641ad630bed35fed0199d015a90b7af6ba9210
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/discrepancy/Global_Disc.py b/src/scene_verifier/dryvr/discrepancy/Global_Disc.py
new file mode 100644
index 0000000000000000000000000000000000000000..d84b5176e7c3ae8aca917af0ea533e5f6aab137d
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/discrepancy/PW_Discrepancy.py b/src/scene_verifier/dryvr/discrepancy/PW_Discrepancy.py
new file mode 100644
index 0000000000000000000000000000000000000000..53cac39418bfb66277fc7a37b25387648e5c076e
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/discrepancy/__init__.py b/src/scene_verifier/dryvr/discrepancy/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/scene_verifier/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb b/src/scene_verifier/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2fd64429bf421126b7000c94ce0f6fd186fbd01f
--- /dev/null
+++ b/src/scene_verifier/dryvr/plotter/.ipynb_checkpoints/plotterInterface-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/scene_verifier/dryvr/plotter/README.md b/src/scene_verifier/dryvr/plotter/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..2ea9534697e8de9a5fa0ab9831ce4102499684ab
--- /dev/null
+++ b/src/scene_verifier/dryvr/plotter/README.md
@@ -0,0 +1 @@
+This folder consist plotter code for DryVR reachtube output
diff --git a/src/scene_verifier/dryvr/plotter/__init__.py b/src/scene_verifier/dryvr/plotter/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..5dbc287d927065e02e00ed2cccd0389194764818
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/plotter/linkednode.py b/src/scene_verifier/dryvr/plotter/linkednode.py
new file mode 100644
index 0000000000000000000000000000000000000000..21ee9f462f066b926a39cea050e034db2ab3c078
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/plotter/parser.py b/src/scene_verifier/dryvr/plotter/parser.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d84ecd8a55fac0e7ab966cbfb3774d1fd30b4bc
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/plotter/plot.py b/src/scene_verifier/dryvr/plotter/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..b07e575197d200745f9a83a5fbaa866bc6229fb1
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/plotter/plot_util.py b/src/scene_verifier/dryvr/plotter/plot_util.py
new file mode 100644
index 0000000000000000000000000000000000000000..6718c923b4dd9148d8f6b55464a662daeab2a30f
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/dryvr/plotter/plotter_interface.ipynb b/src/scene_verifier/dryvr/plotter/plotter_interface.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..de6c0d1e2aa8693652246f3ea36fa58d135ed351
--- /dev/null
+++ b/src/scene_verifier/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/src/scene_verifier/map/__init__.py b/src/scene_verifier/map/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c4fd24f6b4b1fa640294337d048c0fc8d7051183
--- /dev/null
+++ b/src/scene_verifier/map/__init__.py
@@ -0,0 +1,8 @@
+# try:
+#     from lane_map import LaneMap
+#     from single_straight_lane import SingleStraightLaneMap
+#     from lane_segment import LaneSegment
+# except:
+#     from scene_verifier.map.lane_segment import LaneSegment
+#     # from scene_verifier.map.lane_map import LaneMap
+#     # from scene_verifier.map.single_straight_lane import SingleStraightLaneMap
\ No newline at end of file
diff --git a/src/scene_verifier/map/lane_map.py b/src/scene_verifier/map/lane_map.py
new file mode 100644
index 0000000000000000000000000000000000000000..9e54323250cfc4cd834cae5ec083defbce9bb075
--- /dev/null
+++ b/src/scene_verifier/map/lane_map.py
@@ -0,0 +1,71 @@
+from typing import Dict, List
+import copy
+from enum import Enum,auto
+
+from src.scene_verifier.map.lane_segment import LaneSegment
+
+class LaneMap:
+    def __init__(self, lane_seg_list:List[LaneSegment] = []):
+        self.lane_segment_dict:Dict[str, LaneSegment] = {}
+        self.left_lane_dict:Dict[str, List[str]] = {}
+        self.right_lane_dict:Dict[str, List[str]] = {}
+        for lane_seg in lane_seg_list:
+            self.lane_segment_dict[lane_seg.id] = lane_seg
+            self.left_lane_dict[lane_seg.id] = []
+            self.right_lane_dict[lane_seg.id] = []
+
+    def add_lanes(self, lane_seg_list:List[LaneSegment]):
+        for lane_seg in lane_seg_list:
+            self.lane_segment_dict[lane_seg.id] = lane_seg
+            self.left_lane_dict[lane_seg.id] = []
+            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
+        left_lane_list = self.left_lane_dict[lane_idx]
+        return len(left_lane_list)>0
+
+    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
+        right_lane_list = self.right_lane_dict[lane_idx]
+        return len(right_lane_list)>0
+
+    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):
+        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/src/scene_verifier/map/lane_segment.py b/src/scene_verifier/map/lane_segment.py
new file mode 100644
index 0000000000000000000000000000000000000000..33143c69e68a73db84ce7a1ec851aa79e7382754
--- /dev/null
+++ b/src/scene_verifier/map/lane_segment.py
@@ -0,0 +1,15 @@
+from typing import List
+
+class LaneSegment:
+    def __init__(self, id, lane_parameter = None):
+        self.id = id
+        # self.left_lane:List[str] = left_lane
+        # self.right_lane:List[str] = right_lane 
+        # self.next_segment:int = next_segment
+
+        self.lane_parameter = None 
+        if lane_parameter is not None:
+            self.lane_parameter = lane_parameter
+
+    def get_geometry(self):
+        return self.lane_parameter
\ No newline at end of file
diff --git a/src/scene_verifier/scenario/__init__.py b/src/scene_verifier/scenario/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..559069149e378ddbcab6cecd08c7ec429a29962f
--- /dev/null
+++ b/src/scene_verifier/scenario/__init__.py
@@ -0,0 +1,4 @@
+# try:
+#     from Scenario import *
+# except:
+#     from scene_verifier.scenario.Scenario import *
\ No newline at end of file
diff --git a/src/scene_verifier/scenario/scenario.py b/src/scene_verifier/scenario/scenario.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ad1a5db8760aba1567b540ec309c3fb3b658b5f
--- /dev/null
+++ b/src/scene_verifier/scenario/scenario.py
@@ -0,0 +1,236 @@
+from typing import Tuple
+import copy
+import itertools
+
+import numpy as np
+
+from src.scene_verifier.agents.base_agent import BaseAgent
+from src.scene_verifier.automaton.guard import GuardExpressionAst
+from src.scene_verifier.automaton.reset import ResetExpression
+from src.scene_verifier.code_parser.pythonparser import Guard
+from src.scene_verifier.code_parser.pythonparser import Reset
+from src.scene_verifier.analysis.simulator import Simulator
+from src.scene_verifier.analysis.verifier import Verifier
+from src.scene_verifier.map.lane_map import LaneMap
+
+class Scenario:
+    def __init__(self):
+        self.agent_dict = {}
+        self.simulator = Simulator()
+        self.verifier = Verifier()
+        self.init_dict = {}
+        self.init_mode_dict = {}
+        self.map = None
+        self.sensor = None
+
+    def set_sensor(self, sensor):
+        self.sensor = sensor
+
+    def add_map(self, lane_map:LaneMap):
+        self.map = lane_map
+
+    def add_agent(self, agent:BaseAgent):
+        self.agent_dict[agent.id] = agent
+
+    def set_init(self, init_list, init_mode_list):
+        assert len(init_list) == len(self.agent_dict)
+        assert len(init_mode_list) == len(self.agent_dict)
+        for i,agent_id in enumerate(self.agent_dict.keys()):
+            self.init_dict[agent_id] = copy.deepcopy(init_list[i])
+            self.init_mode_dict[agent_id] = copy.deepcopy(init_mode_list[i])
+
+    def simulate(self, time_horizon):
+        init_list = []
+        init_mode_list = []
+        agent_list = []
+        for agent_id in self.agent_dict:
+            init_list.append(self.init_dict[agent_id])
+            init_mode_list.append(self.init_mode_dict[agent_id])
+            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 check_guard_hit(self, state_dict):
+        lane_map = self.map 
+        guard_hits = []
+        any_contained = False        # TODO: Handle this
+        for agent_id in state_dict:
+            agent:BaseAgent = self.agent_dict[agent_id]
+            agent_state, agent_mode = state_dict[agent_id]
+        
+            t = agent_state[0]
+            agent_state = agent_state[1:]
+            paths = agent.controller.getNextModes(agent_mode)
+            for path in paths:
+                # Construct the guard expression
+                guard_list = []
+                reset_list = []
+                for item in path:
+                    if isinstance(item, Guard):
+                        guard_list.append(item)
+                    elif isinstance(item, Reset):
+                        reset_list.append(item)
+                # guard_expression = GuardExpression(guard_list=guard_list)
+                guard_expression = GuardExpressionAst(guard_list)
+                # Map the values to variables using sensor
+                continuous_variable_dict, discrete_variable_dict = self.sensor.sense(self, agent, state_dict, self.map)
+                
+                '''Execute functions related to map to see if the guard can be satisfied'''
+                '''Check guards related to modes to see if the guards can be satisfied'''
+                '''Actually plug in the values to see if the guards can be satisfied'''
+                # Check if the guard can be satisfied
+                guard_can_satisfied = guard_expression.evaluate_guard_disc(agent, discrete_variable_dict, self.map)
+                if not guard_can_satisfied:
+                    continue
+                guard_satisfied, is_contained = guard_expression.evaluate_guard_cont(agent, continuous_variable_dict, self.map)
+                any_contained = any_contained or is_contained
+                if guard_satisfied:
+                    guard_hits.append((agent_id, guard_list, reset_list))
+        return guard_hits, any_contained
+
+    def get_all_transition_set(self, node):
+        possible_transitions = []
+        trace_length = int(len(list(node.trace.values())[0])/2)
+        guard_hits = []
+        guard_hit_bool = False
+
+        # TODO: can add parallalization for this loop
+        for idx in range(0,trace_length):
+            # 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*2:idx*2+2], node.mode[agent_id])
+            hits, is_contain = self.check_guard_hit(all_agent_state)
+            if hits != []:
+                guard_hits.append((hits, all_agent_state, idx))
+                guard_hit_bool = True
+            if hits == [] and guard_hit_bool:
+                break
+            if is_contain:
+                break
+
+        reset_dict = {}
+        reset_idx_dict = {}
+        for hits, all_agent_state, hit_idx in guard_hits:
+            for agent_id, guard_list, reset_list in hits:
+                dest_list,reset_rect = self.apply_reset(node.agent[agent_id], reset_list, all_agent_state)
+                if agent_id not in reset_dict:
+                    reset_dict[agent_id] = {}
+                    reset_idx_dict[agent_id] = {}
+                for dest in dest_list:
+                    if dest not in reset_dict[agent_id]:
+                        reset_dict[agent_id][dest] = []
+                        reset_idx_dict[agent_id][dest] = []
+                    reset_dict[agent_id][dest].append(reset_rect)
+                    reset_idx_dict[agent_id][dest].append(hit_idx)
+            
+        # Combine reset rects and construct transitions
+        for agent in reset_dict:
+            for dest in reset_dict[agent]:
+                combined_rect = None 
+                for rect in reset_dict[agent][dest]:
+                    rect = np.array(rect)
+                    if combined_rect is None:
+                        combined_rect = rect 
+                    else:
+                        combined_rect[0,:] = np.minimum(combined_rect[0,:], rect[0,:])
+                        combined_rect[1,:] = np.maximum(combined_rect[1,:], rect[1,:])
+                combined_rect = combined_rect.tolist()
+                min_idx = min(reset_idx_dict[agent][dest])
+                max_idx = max(reset_idx_dict[agent][dest])
+                transition = (agent, node.mode[agent], dest, combined_rect, (min_idx, max_idx))
+                possible_transitions.append(transition)
+        # Return result
+        return possible_transitions
+
+    def apply_reset(self, agent, reset_list, all_agent_state) -> Tuple[str, np.ndarray]:
+        reset_expr = ResetExpression(reset_list)
+        continuous_variable_dict, discrete_variable_dict = self.sensor.sense(self, agent, all_agent_state, self.map)
+        dest = reset_expr.get_dest(agent, all_agent_state[agent.id], discrete_variable_dict, self.map)
+        rect = reset_expr.apply_reset_continuous(agent, continuous_variable_dict, self.map)
+        return dest, rect
+
+    def get_all_transition(self, state_dict):
+        lane_map = self.map
+        satisfied_guard = []
+        for agent_id in state_dict:
+            agent:BaseAgent = self.agent_dict[agent_id]
+            agent_state, agent_mode = state_dict[agent_id]
+            t = agent_state[0]
+            agent_state = agent_state[1:]
+            paths = agent.controller.getNextModes(agent_mode)
+            for path in paths:
+                # Construct the guard expression
+                guard_list = []
+                reset_list = []
+                for item in path:
+                    if isinstance(item, Guard):
+                        guard_list.append(item)
+                    elif isinstance(item, Reset):
+                        reset_list.append(item)
+                # guard_expression = GuardExpression(guard_list=guard_list)
+                guard_expression = GuardExpressionAst(guard_list)
+                # Map the values to variables using sensor
+                continuous_variable_dict, discrete_variable_dict = self.sensor.sense(self, agent, state_dict, self.map)
+                
+                '''Execute functions related to map to see if the guard can be satisfied'''
+                '''Check guards related to modes to see if the guards can be satisfied'''
+                '''Actually plug in the values to see if the guards can be satisfied'''
+                # Check if the guard can be satisfied
+                guard_satisfied = guard_expression.evaluate_guard(agent, continuous_variable_dict, discrete_variable_dict, self.map)
+                if guard_satisfied:
+                    # If the guard can be satisfied, handle resets
+                    next_init = agent_state
+                    dest = agent_mode.split(',')
+                    possible_dest = [[elem] for elem in dest]
+                    for reset in reset_list:
+                        # Specify the destination mode
+                        reset = reset.code
+                        if "mode" in reset:
+                            for i, discrete_variable_ego in enumerate(agent.controller.vars_dict['ego']['disc']):
+                                if discrete_variable_ego in reset:
+                                    break
+                            tmp = reset.split('=')
+                            if 'map' in tmp[1]:
+                                tmp = tmp[1]
+                                for var in discrete_variable_dict:
+                                    tmp = tmp.replace(var, f"'{discrete_variable_dict[var]}'")
+                                possible_dest[i] = eval(tmp)
+                            else:
+                                tmp = tmp[1].split('.')
+                                if tmp[0].strip(' ') in agent.controller.modes:
+                                    possible_dest[i] = [tmp[1]]                            
+                        else: 
+                            for i, cts_variable in enumerate(agent.controller.vars_dict['ego']['cont']):
+                                if "output."+cts_variable in reset:
+                                    break 
+                            tmp = reset.split('=')
+                            tmp = tmp[1]
+                            for cts_variable in continuous_variable_dict:
+                                tmp = tmp.replace(cts_variable, str(continuous_variable_dict[cts_variable]))
+                            next_init[i] = eval(tmp)
+                    all_dest = itertools.product(*possible_dest)
+                    for dest in all_dest:
+                        dest = ','.join(dest)
+                        next_transition = (
+                            agent_id, agent_mode, dest, next_init, 
+                        )
+                        satisfied_guard.append(next_transition)
+
+        return satisfied_guard
\ No newline at end of file
diff --git a/toythermo.py b/toythermo.py
deleted file mode 100644
index 34a85817fe20c3ee18897c550a3c9e6fbed71508..0000000000000000000000000000000000000000
--- a/toythermo.py
+++ /dev/null
@@ -1,135 +0,0 @@
-#python logic toy code:
-
-from enum import Enum, auto
-from ssl import VERIFY_X509_PARTIAL_CHAIN
-from statistics import mode
-from termios import VSTART
-
-
-class Modes(Enum):
-    NormalA = auto()
-    NormalB = auto()
-    NormalC = auto()
-    NormalD = auto()
-
-class State:
-    posx = 0.0
-    posy = 0.0
-    vx = 1.0
-    vy = 1.0
-    mode = Modes.NormalA
-
-    def __init__(self):
-        self.data = []
-
-	
-if __name__ == "__main__":	
-#todo: how would this actually be given
-    s = State()
-
-    posx = s.posx
-    posy = s.posy
-    vx = s.vx
-    vy = s.vy
-    state = s.mode
-    if (state ==Modes.NormalA):
-        if posy<0 and posy>=-0.01: 
-            vy=-vy
-            posy=0
-            state=Modes.NormalA
-            
-        if (posx<0 and posx>=-0.01) :
-                vx=-vx
-                posx=0
-                state=Modes.NormalB
-            
-        if (posx<=5.01 and posx>5) :
-            vx=-vx
-            posx=5
-            state=Modes.NormalC
-            
-        if (posy>5 and posy<=5.01) :
-            vy=-vy
-            posy=5
-            state = Modes.NormalD
-            
-        
-        
-    if (state ==Modes.NormalB) :
-        if (posy<0 and posy>=-0.01) :
-            vy=-vy
-            posy=0
-            state=Modes.NormalB
-            
-        if (posx<0 and posx>=-0.01) :
-            vx=-vx
-            posx=0
-            state=Modes.NormalA
-            
-        if (posx<=5.01 and posx>5) :
-            vx=-vx
-            posx=5
-            state=Modes.NormalC
-            
-        if (posy>5 and posy<=5.01) :
-            vy=-vy
-            posy=5
-            state = Modes.NormalD
-            
-        
-        
-    if (state ==Modes.NormalC) :
-        if (posy<0 and posy>=-0.01) :
-            vy=-vy
-            posy=0
-            state=Modes.NormalC
-            
-        if (posx<0 and posx>=-0.01) :
-            vx=-vx
-            posx=0
-            state=Modes.NormalA
-            
-        if (posx<=5.01 and posx>5) :
-            vx=-vx
-            posx=5
-            state=Modes.NormalB
-            
-        if (posy>5 and posy<=5.01) :
-            vy=-vy
-            posy=5
-            state = Modes.NormalD
-            
-        
-        
-        
-    if (state ==State.Modes.NormalD) :
-        if (posy<0 and posy>=-0.01) :
-            vy=-vy
-            posy=0
-            state=Modes.NormalD
-            
-        if (posx<0 and posx>=-0.01) :
-            vx=-vx
-            posx=0
-            state=Modes.NormalA
-            
-        if (posx<=5.01 and posx>5) :
-            vx=-vx
-            posx=5
-            state=Modes.NormalB
-            
-        if (posy>5 and posy<=5.01) :
-            vy=-vy
-            posy=5
-            state = Modes.NormalC
-            
-        
-        
-    s.posx = posx
-    s.posy= posy
-    s.vx = vx
-    s.vy = vy
-    s.mode = state
-
-    #TODO: what is output?
-
diff --git a/toythermomini.py b/toythermomini.py
deleted file mode 100644
index 0423c8c84ddb00372658d640e9c30d39ea5e0904..0000000000000000000000000000000000000000
--- a/toythermomini.py
+++ /dev/null
@@ -1,45 +0,0 @@
-#python logic toy code:
-
-from enum import Enum, auto
-
-class Modes(Enum):
-    NormalA = auto()
-    NormalB = auto()
-    NormalC = auto()
-    NormalD = auto()
-
-class State:
-    posx = 0.0
-    posy = 0.0
-    mode = Modes.NormalA
-
-    def __init__(self):
-        self.data = []
-
-	
-def controller(posx, posy, mode):	
-    outmode = mode
-#todo: how would this actually be given
-    if (mode ==Modes.NormalA):
-        if posy<0 and posy>=-0.0: 
-            posy=0
-            outmode=Modes.NormalA
-        if posy>10 and posy==-10: 
-            posy=10
-            outmode=Modes.NormalB
-
-    if (mode ==Modes.NormalB):
-        if posy<0 and posy>=-0.0: 
-            posy=0
-            outmode=Modes.NormalA
-        if posy>10 and posy==-10: 
-            posy=10
-            outmode=Modes.NormalB
-            
-    
-
-
-    return outmode
-
-    #TODO: what is output?
-