diff --git a/model/spatial_queue_array_updated.py b/model/spatial_queue_array_updated.py
new file mode 100755
index 0000000..5f7556c
--- /dev/null
+++ b/model/spatial_queue_array_updated.py
@@ -0,0 +1,547 @@
+#!/usr/bin/env python
+# coding: utf-8
+import os, sys
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+from shapely.wkt import loads
+from shapely.geometry import Point, LineString
+from networkx import Graph, DiGraph, simple_cycles
+
+### usr
+abs_path = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, abs_path)
+sys.path.insert(0, "/home/bingyu/Documents")
+from sp import interface
+
+np.random.seed(1)
+
+class Network:
+    def __init__(self):
+        self.nodes = None
+        self.links = None
+        ### no conflict movement directions by node
+        self.no_conflict_moves = dict()
+        ### simulation/shortest path network
+        self.g = None
+        ### some frequently used columns
+        self.link_nid_dict = dict()
+        self.link_length_dict = dict()
+        self.link_fft_dict = dict()
+        self.link_lanes_dict = dict()
+        self.link_capacity_dict = dict()
+    
+    def prepare_nodes(self, nodes_df, ods_df):
+        ### nodes_df must contain the following columns
+        self.nodes = nodes_df[['node_id', 'lon', 'lat']].copy()
+        ### if virtual node is hard coded in the input data
+        self.nodes['node_type'] = 'real'
+
+        ### add virtual nodes
+        virtual_nodes_df = self.nodes.loc[self.nodes['node_id'].isin(np.unique(ods_df['origin_nid']))].copy()
+        virtual_nodes_df['node_id'] = virtual_nodes_df['node_id'].apply(lambda x: 'vn_{}'.format(x))
+        virtual_nodes_df['node_type'] = 'virtual'
+        self.nodes = pd.concat([self.nodes, virtual_nodes_df])
+        return virtual_nodes_df
+    
+    def prepare_links(self, links_df, virtual_nodes_df):
+        ### links must contain the following columns; geometry must be in epsg:4326
+        self.links = links_df[['link_id', 'start_node_id', 'end_node_id', 'geometry']].copy()
+        self.links['link_type'] = 'real'
+        
+        ### add attributes
+        self.links = gpd.GeoDataFrame(self.links, crs='epsg:4326', geometry=links_df['geometry'].map(loads)).to_crs('epsg:3857')
+        try: self.links['length'] = links_df['length'] 
+        except KeyError:
+            self.links['length'] = self.links['geometry'].length
+            print('link length not specified; derive from geometry instead.')
+        try: self.links['maxmph'] = links_df['maxmph']
+        except KeyError:
+            self.links['maxmph'] = 25
+            print('link speed limit (mph); use 25mph.')
+        try: self.links['lanes'] = links_df['lanes']
+        except KeyError:
+            self.links['lanes'] = 1
+            print('link length not specified; assume one lane.')
+        
+        ### add virtual links
+        virtual_links_df = virtual_nodes_df[['node_id', 'lon', 'lat']].copy()
+        virtual_links_df['start_node_id'] = virtual_links_df['node_id']
+        virtual_links_df['end_node_id'] = virtual_links_df['node_id'].apply(lambda x: int(x.replace('vn_', ''))) ### remove the 'vn_' from node id
+        virtual_links_df['link_id'] = virtual_links_df['node_id'].apply(lambda x: x.replace('vn_', 'vl_'))
+        virtual_links_df['link_type'] = 'virtual'
+        virtual_links_df['geometry'] = None
+        virtual_links_df['length'] = 1
+        virtual_links_df['maxmph'] = 1000
+        virtual_links_df['lanes'] = 1e8
+
+        ### add other attributes
+        self.links = pd.concat([self.links, virtual_links_df[self.links.columns]])
+        self.links['fft'] = self.links['length']/(self.links['maxmph']*1609/3600)
+        self.links['capacity'] = self.links['lanes'] * 1800/3600
+        self.links['storage'] = self.links['length'] * self.links['lanes']
+        self.links['storage'] = np.where(self.links['storage']<18, 18, self.links['storage'])
+        self.links['storage_add_times'] = 0 # default is 0
+        self.links['rtt'] = 0
+        self.links['rtt'] = np.where(self.links['rtt'] < self.links['fft'], self.links['fft'], self.links['rtt']) # real average travel time, updated every 10 mins, if equals 0, let it be the fft
+
+        self.links = self.links[['link_id', 'start_node_id', 'end_node_id', 'link_type', 'length', 'lanes', 'maxmph', 'fft', 'capacity', 'storage', 'geometry','storage_add_times']]
+
+    def node_movements(self):
+
+        link_out_coord = dict()
+        link_in_coord = dict()
+        node_coords = dict()
+        incoming_links = {n: [] for n in self.nodes['node_id'].values.tolist()}
+        outgoing_links = {n: [] for n in self.nodes['node_id'].values.tolist()}
+        
+        for link in self.links[self.links['link_type']=='real'].itertuples():
+            ### virtual incoming links not considered here
+            link_id = getattr(link, 'link_id')
+            link_geometry = getattr(link, 'geometry')
+            link_snid = getattr(link, 'start_node_id')
+            link_enid = getattr(link, 'end_node_id')
+            
+            ### near end coordinates
+            ### first upstream node except from the start node
+            link_in_coord[link_id] = link_geometry.coords[1]
+            ### last downstream node except from the end node
+            link_out_coord[link_id] = link_geometry.coords[-2]
+            
+            ### node coordinates
+            if link_snid not in node_coords: node_coords[link_snid] = link_geometry.coords[0]
+            if link_enid not in node_coords: node_coords[link_enid] = link_geometry.coords[-1] 
+            
+            ### node to all links
+            incoming_links[link_enid].append(link_id)
+            outgoing_links[link_snid].append(link_id)
+        
+        #display(self.nodes[self.nodes['node_id']==194563])
+        for node in self.nodes[self.nodes['node_type']=='real'].itertuples():
+            node_id = getattr(node, 'node_id')
+            nc = node_coords[node_id]
+            node_movements_info = dict()
+            for il in incoming_links[node_id]:
+                for ol in outgoing_links[node_id]:
+                    ### virtual incoming links not considered here
+                    ilc, olc = link_out_coord[il], link_in_coord[ol]
+                    move_1 = np.arctan2(nc[1]-ilc[1], nc[0]-ilc[0])*180/np.pi
+                    move_2 = np.arctan2(olc[1]-nc[1], olc[0]-nc[0])*180/np.pi
+                    if move_1 < move_2: move_2 -= 360
+                    if (move_1 - move_2)<=45 or (move_1 - move_2)>315: move_dir='sa'
+                    elif (move_1 - move_2)>45 and (move_1 - move_2)<170: move_dir='r'
+                    elif (move_1 - move_2)>=170 and (move_1 - move_2)<=315: move_dir='left' ### left or u turn
+                    else: move_dir = None
+                    #if node_id==29: print(il, ol, ilc, nc, olc, move_1, move_2, move_dir)
+                    node_movements_info[(il, ol)] = [ilc, olc, move_dir]
+            ### pair-wise evaluation of movements in this intersection
+            for (il, ol), [ilc, olc, move_dir] in node_movements_info.items():
+                for (il_, ol_), [ilc_, olc_, move_dir_] in node_movements_info.items():
+                    if (il==il_) and (ol==ol_):
+                        ### same direction can go together
+                        self.no_conflict_moves[(il, ol, il_, ol_)] = 1
+                        continue
+                    if (il==il_) and (ol!=ol_):
+                        ### same in direction but different out direction
+                        if (move_dir in ['left']) or (move_dir_ in ['left']):
+                            continue
+                        else: self.no_conflict_moves[(il, ol, il_, ol_)] = 1
+                    ### find opposite direction
+                    move_1 = np.arctan2(nc[1]-ilc[1], nc[0]-ilc[0])*180/np.pi
+                    move_2 = np.arctan2(nc[1]-ilc_[1], nc[0]-ilc_[0])*180/np.pi
+                    if move_1 < move_2: move_2 -= 360
+                    if (move_1-move_2)>=170 and (move_1-move_2)<=190:
+                        ### opposite direction
+                        if (move_dir in ['left']) and (move_dir_ in ['sa', 'r']):
+                            continue
+                        if (move_dir in ['sa', 'r']) and (move_dir_ in ['left']):
+                            continue
+                        #if (move_dir in ['left']) and (move_dir_ in ['left']):
+                        #    continue
+                        self.no_conflict_moves[(il, ol, il_, ol_)] = 1
+                    else:
+                        continue
+    
+    def prepare_network(self, nodes_df, links_df, ods_df):
+        
+        ### add virtual nodes and virtual links
+        ### only return virtual nodes as it will be used in the creation of virtual links in the next function
+        virtual_nodes_df = self.prepare_nodes(nodes_df, ods_df) 
+        self.prepare_links(links_df, virtual_nodes_df)
+        
+        ### find no conflict movements
+        self.node_movements()
+
+        ### some lookup or calculation variables
+        self.g = interface.from_dataframe(self.links[self.links['link_type']=='real'], 'start_node_id', 'end_node_id', 'fft')
+        self.link_nid_dict = {(getattr(link, 'start_node_id'), getattr(link, 'end_node_id')): getattr(link, 'link_id') for link in self.links.itertuples()}
+        self.link_length_dict = {getattr(link, 'link_id'): getattr(link, 'length') for link in self.links.itertuples()}
+        self.link_lanes_dict = {getattr(link, 'link_id'): getattr(link, 'lanes') for link in self.links.itertuples()}
+        self.link_fft_dict = {getattr(link, 'link_id'): getattr(link, 'fft') for link in self.links.itertuples()}
+        self.link_capacity_dict = {getattr(link, 'link_id'): getattr(link, 'capacity') for link in self.links.itertuples()}
+
+    def find_next_link(self, agent, agent_routes, just_loaded=False):
+        agent_id = int(getattr(agent, 'agent_id'))
+        agent_cl_enid = int(getattr(agent, 'cl_enid'))
+        ### just loaded
+        if just_loaded:
+            just_loaded_nl_enid = agent_routes[agent_id][agent_cl_enid]
+            just_loaded_nl = self.link_nid_dict[(agent_cl_enid, just_loaded_nl_enid)]
+            just_loaded_nl_lanes = self.link_lanes_dict[just_loaded_nl]
+            just_loaded_nl_capacity = self.link_capacity_dict[just_loaded_nl]
+            return ['vl_{}'.format(agent_cl_enid), 'vn_{}'.format(agent_cl_enid), agent_cl_enid, 0.1, 1e8, 1e8/3.6, just_loaded_nl, just_loaded_nl_enid, just_loaded_nl_lanes, just_loaded_nl_capacity]
+        ### current link end is destination
+        if agent_cl_enid == int(getattr(agent, 'destin_nid')):
+            return ['vl_sink', agent_cl_enid, 'sink_node', 0.1, 1e8, 1e8/3.6, None, None, None, None]
+        #print('current_agent_id:{}, current_cl_enid:{}, current_agent_route:{}'.format(agent_id, agent_cl_enid, agent_routes[agent_id]))
+        agent_nl_enid = agent_routes[agent_id][agent_cl_enid] #agent_routes[agent_id][agent_cl_enid]
+        agent_nl = self.link_nid_dict[(agent_cl_enid, agent_nl_enid)] #self.link_nid_dict[(agent_cl_enid, agent_nl_enid)]
+        agent_nl_fft = self.link_fft_dict[agent_nl]
+        agent_nl_lanes = self.link_lanes_dict[agent_nl]
+        agent_nl_capacity = self.link_capacity_dict[agent_nl]
+        ### next link end is destination, next next link does not exist
+        if agent_nl_enid == int(getattr(agent, 'destin_nid')):
+            #return [agent_nl, agent_cl_enid, agent_nl_enid, agent_nl_fft, agent_nl_lanes, agent_nl_capacity, 'vl_sink', 'sink_node', 1e8, 1e8/3.6]
+            return [agent_nl, agent_cl_enid, agent_nl_enid, agent_nl_fft, 1000, 1e7/3.6, 'vl_sink', 'sink_node', 1e8, 1e8/3.6]
+        #print('agent_id:{}, cl_enid: {}, nl_enid: {}, agent_routes:{}'.format(agent_id, agent_cl_enid, agent_nl_enid, agent_routes[agent_id]))
+        agent_nnl_enid = agent_routes[agent_id][agent_nl_enid]
+        agent_nnl = self.link_nid_dict[(agent_nl_enid, agent_nnl_enid)]
+        agent_nnl_lanes = self.link_lanes_dict[agent_nnl]
+        agent_nnl_capacity = self.link_capacity_dict[agent_nnl]
+        return [agent_nl, agent_cl_enid, agent_nl_enid, agent_nl_fft, agent_nl_lanes, agent_nl_capacity, agent_nnl, agent_nnl_enid, agent_nnl_lanes, agent_nnl_capacity]
+
+    def run_link_model(self, agents=None, t=None):
+
+        agents_df = agents.agents
+        
+        ### running vehicles joins queue
+        agents_df['agent_status'] = np.where(
+            (agents_df['agent_status']==1) &
+            (agents_df['cl_fft']<(t-agents_df['current_link_enter_time'])),
+            2, agents_df['agent_status'])
+
+        ### summarize from agents_df
+        link_undeparted_dict = agents_df.loc[agents_df['agent_status']==0].groupby('current_link').size().to_dict()
+        link_queue_dict = agents_df.loc[agents_df['agent_status']==2].groupby('current_link').size().to_dict()
+        link_run_dict = agents_df.loc[agents_df['agent_status']==1].groupby('current_link').size().to_dict()
+ 
+        ### update links
+        self.links['undeparted'] = self.links['link_id'].map(link_undeparted_dict).fillna(0)
+        self.links['queue'] = self.links['link_id'].map(link_queue_dict).fillna(0)
+        self.links['run'] = self.links['link_id'].map(link_run_dict).fillna(0)
+
+        if t % 900 == 0:
+            link_travel_time_dict = ((t * agents_df.loc[agents_df['agent_status'].isin([2])].groupby('current_link').size() - 
+                                    agents_df.loc[agents_df['agent_status'].isin([2])].groupby('current_link')['current_link_enter_time'].sum())/
+                                    agents_df.loc[agents_df['agent_status'].isin([2])].groupby('current_link').size()).to_dict()
+            self.links['rtt'] = self.links['link_id'].map(link_travel_time_dict).fillna(0)
+            self.links['rtt'] = np.where(self.links['rtt'] < self.links['fft'], self.links['fft'], self.links['rtt'])
+            self.g = interface.from_dataframe(self.links[self.links['link_type']=='real'], 'start_node_id', 'end_node_id', 'rtt') # update g
+
+        current_time_storage_remain = (self.links['storage'] - (self.links['run'] + self.links['queue'])*8).values
+        self.links['storage_remain'] = current_time_storage_remain
+
+
+        return agents_df
+    
+    def run_node_model(self, agents=None, t=None, special_nodes=None):
+        
+        agents_df = agents.agents
+        agent_routes = agents.agent_routes
+        
+        ### queue vehicles
+        go_agents_df = agents_df.loc[
+            agents_df['agent_status']==2].copy()
+        if go_agents_df.shape[0] == 0:
+            ### no queuing agents 
+            return agents_df
+        
+        ### choose primary movement direction
+        ### 1. first come first serve
+        #go_agents_df = go_agents_df.sort_values(by='current_link_enter_time', ascending=True)
+        ### 2. random (e.g., a random traffic light)
+        go_agents_df = go_agents_df.sample(frac=1)
+        ### 3.highest control capacity can go
+        #go_agents_df = go_agents_df.sort_values(by='control_capacity', ascending=True)
+        ### user-specified movement priority at nodes (e.g., roundabouts)
+        if len(special_nodes) != 0:
+            go_agents_midx = pd.MultiIndex.from_frame(go_agents_df[['current_link', 'next_link']])
+            go_agents_df['default_priority'] = 3
+            go_agents_df['default_priority'] = np.where(go_agents_midx.isin(
+                special_nodes), 1, go_agents_df['default_priority'])
+
+            go_agents_df = go_agents_df.sort_values(by='default_priority', ascending=True)
+        ### the results: directions of the first vehicle in the queue after sorting
+        go_agents_df['primary_cl'] = go_agents_df.groupby('cl_enid', sort=False)['current_link'].transform('first')
+        go_agents_df['primary_nl'] = go_agents_df.groupby('cl_enid', sort=False)['next_link'].transform('first')
+        
+        ### limit by capacity
+        go_agents_df['control_capacity'] = np.floor(
+            go_agents_df[['cl_capacity', 'nl_capacity']].min(axis=1) ### mostly zero (except highways)
+            + np.random.random(size=go_agents_df.shape[0]))
+        #print(go_agents_df.head(1))
+        go_agents_df['control_capacity'] = go_agents_df.groupby(
+            ['current_link', 'next_link'])['control_capacity'].transform('first')
+        ### filter for vehicles that can move according to the control capacity
+        #print(go_agents_df.head(1))
+        go_agents_df['current_link_queue_time'] = t - go_agents_df['cl_fft'] - go_agents_df['current_link_enter_time']
+        #print(go_agents_df.head(1))
+        go_agents_df = go_agents_df.sort_values(by='current_link_queue_time', ascending=False)
+        #print(go_agents_df.head(1))
+        go_agents_df['clnl_position'] = go_agents_df.groupby(['current_link', 'next_link'], sort=False).cumcount()
+        go_agents_df = go_agents_df.loc[
+            go_agents_df['clnl_position']<go_agents_df['control_capacity']]
+        if go_agents_df.shape[0] == 0:
+            return agents_df
+        
+        ### all nodes where the primary movement can go
+        selected_nodes = go_agents_df.loc[
+            (go_agents_df['current_link']==go_agents_df['primary_cl']) &
+            (go_agents_df['next_link']==go_agents_df['primary_nl']), 'cl_enid'
+        ].values
+        go_agents_df = go_agents_df.loc[go_agents_df['cl_enid'].isin(selected_nodes)].copy()
+        if go_agents_df.shape[0] == 0:
+            return agents_df
+        
+        ### label non-conflict movements
+        go_agents_df['no_conflict'] = list(map(self.no_conflict_moves.get, go_agents_df[['primary_cl', 'primary_nl', 'current_link', 'next_link']].to_records(index=False).tolist()))
+        ### try if it is faster?
+        #go_agents_midx = pd.MultiIndex.from_frame(go_agents_df[['primary_cl', 'primary_nl', 'current_link', 'next_link']])
+        #go_agents_df['no_conflict'] = np.where(go_agents_midx.isin(self.no_conflict_moves.keys()), 1, 0)
+        ### Coming from virtual source or going to virtual sink. Assume no conflict with any direction.
+        go_agents_df['no_conflict'] = np.where(
+            (go_agents_df['cl_lanes']>=100) | 
+            (go_agents_df['nl_lanes']>=100), 1, go_agents_df['no_conflict'])
+        go_agents_df = go_agents_df.loc[go_agents_df['no_conflict']==1].copy()
+        if go_agents_df.shape[0] == 0:
+            return agents_df
+
+        ### find next link storage constraints, changes each time step
+        links_to_use_df = self.links.loc[(self.links['link_type']=='real') & (self.links['link_id'].isin(go_agents_df['current_link']) | self.links['link_id'].isin(go_agents_df['next_link']))]
+        links_storage_remain_dict = dict(zip(links_to_use_df['link_id'], links_to_use_df['storage_remain']))
+        ### no constraints if not in dictionary; meaning the links are either virtual source or sink
+        go_agents_df['nl_storage_remain'] = go_agents_df['next_link'].map(links_storage_remain_dict).fillna(1e5)
+        
+        ### keep only vehicles that the link outflow capacity allows (next link inflow is not controlling as two movements goign to the same downstream link is usually considered conflicting)
+        ### no sorting again, keep the same with the order random or higher control capacity first, etc.
+        ### MAYBE PROBLEMATIC FOR HIGHWAYS
+        go_agents_df['cl_position'] = go_agents_df.groupby('current_link', sort=False).cumcount()
+        go_agents_df = go_agents_df.loc[
+            (go_agents_df['cl_position']<go_agents_df['control_capacity']) |
+            (go_agents_df['nl_lanes']>=100)
+        ].copy()
+        go_agents_df['nl_position'] = go_agents_df.groupby('next_link', sort=False).cumcount()
+        go_agents_df = go_agents_df.loc[
+            (go_agents_df['nl_position']<go_agents_df['control_capacity']) &
+            (go_agents_df['nl_position']*8<go_agents_df['nl_storage_remain'])
+        ].copy()
+        
+        if go_agents_df.shape[0] == 0:
+            return agents_df
+        #test = go_agents_df.loc[go_agents_df['cl_enid']==17]
+        #if test.shape[0]>0:
+        #    print(t, test[['current_link', 'next_link', 'control_capacity', 'cl_capacity', 'nl_capacity']].values.tolist())
+        
+        ### update agent position
+        go_ids = go_agents_df['agent_id']
+        agent_clnl_update = {getattr(agent, 'agent_id'): self.find_next_link(agent, agent_routes) for agent in go_agents_df.itertuples()}
+        #if len(keep)>0:
+        #    print(agent_clnl_update[keep[0]])
+        agents_df.loc[agents_df['agent_id'].isin(go_ids), ['current_link', 'cl_snid', 'cl_enid', 'cl_fft', 'cl_lanes', 'cl_capacity', 'next_link', 'nl_enid', 'nl_lanes', 'nl_capacity']] = pd.DataFrame(agents_df.loc[agents_df['agent_id'].isin(go_ids), 'agent_id'].map(agent_clnl_update).to_list()).values
+        ### set agent status to run
+        agents_df.loc[agents_df['agent_id'].isin(go_ids), 'agent_status'] = 1
+        agents_df.loc[agents_df['agent_id'].isin(go_ids), 'current_link_enter_time'] = t
+        agents_df['agent_status'] = np.where(agents_df['current_link']=='vl_sink', -1, agents_df['agent_status'])
+
+        return agents_df
+        
+
+class Agents:
+    def __init__(self):
+        ### hold all agents in a class
+        self.agents = None ### a pandas dataframe
+        self.agent_routes = dict() ### key: agent_id, value: route
+
+    def prepare_agents(self, ods_df):
+        ### ods_df should have at least the following columns
+        ### `nid` correspond to node_id
+        self.agents = ods_df[['origin_nid', 'destin_nid']].copy()
+        ### other columns to initialize
+        self.agents['agent_status'] = 0 ### 0: unloaded; 1: run; 2:queue, -1: arrive; -2:no_route
+        try: self.agents['agent_id'] = ods_df['agent_id']
+        except KeyError:
+            self.agents['agent_id'] = np.arange(self.agents.shape[0])
+            print('Agent IDs not specified; use sequential number.')
+        try: self.agents['departure_time'] = ods_df['departure_time']
+        except KeyError:
+            self.agents['departure_time'] = 0
+            print('agents departure times not specified; assume leave immediately.')
+        for agent_attribute in ['current_link', 'current_link_enter_time', 'cl_snid', 'cl_enid', 'cl_fft', 'cl_lanes', 'next_link', 'nl_enid', 'nl_lanes']:
+            self.agents[agent_attribute] = np.nan
+
+    def get_agent_routes(self, origin, destin, g):
+        ### note this applies to individual agents
+        sp = g.dijkstra(origin, destin)
+        sp_dist = sp.distance(destin)
+        if sp_dist > 1e8:
+            sp.clear()
+            route = {}
+        else:
+            route = {start_node_id: end_node_id for (start_node_id, end_node_id) in sp.route(destin)}
+            sp.clear()
+        return route
+
+    def update_route_info(self, reroute_agents, network, just_loaded=False):
+
+        unarrived_list = []
+        noneedreroute_list = []
+        agent_clnl_update = dict()
+        for agent in reroute_agents.itertuples():
+            ### agent id
+            agent_id = getattr(agent, 'agent_id')
+            if (int(getattr(agent, 'cl_enid')) == int(getattr(agent, 'destin_nid'))) & (just_loaded == False):
+                noneedreroute_list.append(agent_id)
+                #print('noneedreroute_list include {}'.format(agent_id))
+                continue
+            ### get initial or new routes
+            route = self.get_agent_routes(int(getattr(agent, 'cl_enid')), int(getattr(agent, 'destin_nid')), network.g)
+            ### update route
+            self.agent_routes[agent_id] = route
+            ### no path
+            if len(route)==0: 
+                unarrived_list.append(agent_id)
+                print('no route for {}, cl_enid {}, destin_nid {}'.format(agent_id, int(getattr(agent, 'cl_enid')), int(getattr(agent, 'destin_nid'))))
+            else:
+                agent_clnl_update[agent_id] = network.find_next_link(agent, self.agent_routes, just_loaded)
+        if len(unarrived_list)>0:
+            self.agents.loc[self.agents['agent_id'].isin(unarrived_list), 'agent_status'] = -2
+            reroute_agents = reroute_agents.loc[~reroute_agents['agent_id'].isin(unarrived_list), :].copy()
+        if len(noneedreroute_list) > 0:
+            reroute_agents = reroute_agents.loc[~reroute_agents['agent_id'].isin(noneedreroute_list), :].copy()
+        if len(agent_clnl_update) == 0:
+            return
+
+        # because the function of reroutes messes up the order, we cannot directly use pd.DataFrame(reroute_agents['agent_id'].map(xxx)
+        self.agents.loc[(self.agents['agent_id'].isin(reroute_agents['agent_id'])), 
+                            ['current_link', 'cl_snid', 'cl_enid', 'cl_fft', 'cl_lanes', 'cl_capacity',
+                            'next_link', 'nl_enid', 'nl_lanes', 'nl_capacity'
+                            ]] = pd.DataFrame(self.agents.loc[(self.agents['agent_id'].isin(reroute_agents['agent_id'])),
+                                'agent_id'].map(agent_clnl_update).to_list()).values
+
+    def load_agents(self, network=None, t=None):
+        ### load agent
+        load_ids = t == self.agents['departure_time']
+        if np.sum(load_ids) == 0:
+            return
+        self.agents['cl_enid'] = np.where(load_ids, self.agents['origin_nid'], self.agents['cl_enid'])
+        ### loaded to the queue on virtual source links
+        self.agents['agent_status'] = np.where(load_ids, 2, self.agents['agent_status'])
+        self.agents['current_link_enter_time'] = np.where(load_ids, t, self.agents['current_link_enter_time'])
+        #print(self.agents[load_ids])
+        self.update_route_info(self.agents[load_ids], network, just_loaded=True)
+
+    def reroutes(self, network=None, t=None, reroute_frequency=None):
+        
+        if reroute_frequency is None:
+            ### no reroute
+            return
+        else:
+            if (t >=1) & (t%reroute_frequency == 0):
+                ### staggered rerouting
+                reroute_agents = self.agents[(self.agents['agent_status'].isin([1, 2]))].sample(frac=0.2)
+                if len(reroute_agents) == 0:
+                    return
+                self.update_route_info(reroute_agents, network)
+            else:
+                return
+
+
+class Simulation():
+    def __init__(self):
+        self.agents = None
+        self.network = None
+        self.special_nodes = [] #set special nodes to give priority to the agents in the roundabouts
+        self.loops = {} # record the time when each loop exists
+
+    def initialize_simulation(self, nodes_df, links_df, ods_df):
+        
+        ### initialize network
+        self.network = Network()
+        self.network.prepare_network(nodes_df, links_df, ods_df)
+        
+        ### initialize demand
+        self.agents = Agents()
+        self.agents.prepare_agents(ods_df)
+    
+    def find_loops(self, t, agents_df, links_df):
+        self.special_nodes = []
+        loop_agents_df = agents_df.copy()
+        loop_agents_df = loop_agents_df[loop_agents_df['agent_status'].isin([1,2])]
+        loop_agents_df['current_link'] = loop_agents_df['current_link'].astype('str')
+        loop_agents_df['current_link_id_first_letter'] = loop_agents_df['current_link'].apply(lambda x: x[0])
+        loop_agents_df = loop_agents_df.loc[(loop_agents_df['current_link_id_first_letter'] != 'v'), :]
+        loop_agents_df['current_link'] = loop_agents_df['current_link'].astype('float')
+
+        loop_agents_df['next_link'] = loop_agents_df['next_link'].astype('str')
+        loop_agents_df['next_link_id_first_letter'] = loop_agents_df['next_link'].apply(lambda x: x[0])
+        loop_agents_df = loop_agents_df.loc[loop_agents_df['next_link_id_first_letter'] != 'v',  :]
+        loop_agents_df['next_link'] = loop_agents_df['next_link'].astype('float')
+
+        loop_agents_df = loop_agents_df[['agent_id','current_link', 'next_link']]
+        loop_agents_df.reset_index(drop = True, inplace = True)
+
+        link_total = np.unique(loop_agents_df['current_link'].values)
+
+        next_link_dict = [{link_tmp: np.unique(loop_agents_df.loc[loop_agents_df['current_link'] == 
+                        link_tmp, 'next_link'].values).tolist()} for link_tmp in link_total]
+        DG = DiGraph()
+        DG.add_edges_from([(int(getattr(agent_tmp, 'current_link')), int(getattr(agent_tmp, 'next_link'))) for agent_tmp in loop_agents_df.itertuples()])
+        loop_list = list(simple_cycles(DG))
+        if len(loop_list) == 0:
+            print('check: currently, there is no loop in the network')
+            return
+        else:
+            loop_link_total= []
+            for loop_tmp in loop_list:
+                loop_tmp_storage_remain = links_df.loc[links_df['link_id'].isin(loop_tmp),'storage_remain']
+                if len(loop_tmp_storage_remain[loop_tmp_storage_remain <= 0]) == len(loop_tmp_storage_remain):
+                    loop_link_total =  loop_link_total + loop_tmp
+                    if t in self.loops.keys():
+                        self.loops[t].append(loop_tmp)
+                    else:
+                        self.loops[t] = [loop_tmp]
+                    print('check: there exists one deadlock loop')
+                [self.special_nodes.append(tuple((loop_tmp[i], loop_tmp[i+1]))) for i in range(len(loop_tmp)-1)]
+                self.special_nodes.append(tuple((loop_tmp[-1],loop_tmp[0])))
+
+            loop_link_total = np.unique(loop_link_total).tolist()
+            if len(loop_link_total) == 0:
+                print('check: currently, there is no loop in the network')
+                return
+            else:
+                self.network.links['storage_remain'] = np.where(self.network.links['link_id'].isin(loop_link_total), self.network.links['storage_remain'] + 8, self.network.links['storage_remain'])
+                self.network.links['storage_add_times'] = np.where(self.network.links['link_id'].isin(loop_link_total), self.network.links['storage_add_times'] + 1, 
+                                                            self.network.links['storage_add_times'])
+
+
+
+    def run_one_step(self, t, reroute_frequency=None):
+        
+        ### load agents that are ready to departure; find routes and next links for them
+        self.agents.load_agents(network=self.network, t=t)
+
+        ### reroute agents
+        self.agents.reroutes(network=self.network, t=t, reroute_frequency=reroute_frequency)
+
+
+        self.agents.agents = self.network.run_link_model(agents=self.agents, t=t)
+
+        if (t > 1) & (t % 600 == 0):
+            self.find_loops(t=t, agents_df = self.agents.agents, links_df = self.network.links)
+
+
+        self.agents.agents = self.network.run_node_model(agents=self.agents, t=t, special_nodes=self.special_nodes)
+
+
+
+