From 2388fcc6d59df1936aa02b547b0b5ae2107e145d Mon Sep 17 00:00:00 2001 From: jinyan1214 Date: Tue, 19 Nov 2024 11:09:10 -0800 Subject: [PATCH 1/2] jz - Modify Capacity Spectrum Method to copy PGD IMs to response.csv --- .../capacitySpectrum/runCMS.py | 28 +++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/modules/performSIMULATION/capacitySpectrum/runCMS.py b/modules/performSIMULATION/capacitySpectrum/runCMS.py index 085658b46..3533c920f 100644 --- a/modules/performSIMULATION/capacitySpectrum/runCMS.py +++ b/modules/performSIMULATION/capacitySpectrum/runCMS.py @@ -396,11 +396,35 @@ def write_RV(AIM_input_path, EVENT_input_path): # noqa: C901, N802, N803, D103 EDP_output = np.concatenate([index, EDP_output], axis=1) # noqa: N806 + # Concatenate the original IM to CMS in case some (e.g., PGD) are needed + EDP_output = np.concatenate([EDP_output, IM_samples], axis = 1) + # prepare the header + header_out = ['1-PFA-0-0', '1-PRD-1-1'] + for h_label in header: + # remove leading and trailing whitespace + h_label = h_label.strip() # noqa: PLW2901 + + # convert suffixes to the loc-dir format used by the SimCenter + if h_label.endswith('_h'): # horizontal + header_out.append(f'1-{h_label[:-2]}-1-1') + + elif h_label.endswith('_v'): # vertical + header_out.append(f'1-{h_label[:-2]}-1-3') + + elif h_label.endswith('_x'): # x direction + header_out.append(f'1-{h_label[:-2]}-1-1') + + elif h_label.endswith('_y'): # y direction + header_out.append(f'1-{h_label[:-2]}-1-2') + + else: # if none of the above is given, default to 1-1 + header_out.append(f'1-{h_label.strip()}-1-1') + + working_dir = Path(PurePath(EVENT_input_path).parent) # working_dir = posixpath.dirname(EVENT_input_path) - # prepare the header - header_out = ['1-PFA-0-0', '1-PRD-1-1'] + np.savetxt( working_dir / 'response.csv', From b3c93779ff399542cc4f94a97b059abd43f6bb8a Mon Sep 17 00:00:00 2001 From: jinyan1214 Date: Sat, 23 Nov 2024 17:17:54 -0800 Subject: [PATCH 2/2] jz - residual demand mark travel time of incomplete trips as inf --- .../ResidualDemand/run_residual_demand.py | 39 +- .../ResidualDemand/transportation.py | 1159 +++++++++-------- 2 files changed, 643 insertions(+), 555 deletions(-) diff --git a/modules/systemPerformance/ResidualDemand/run_residual_demand.py b/modules/systemPerformance/ResidualDemand/run_residual_demand.py index c2dc17df4..e4261764a 100644 --- a/modules/systemPerformance/ResidualDemand/run_residual_demand.py +++ b/modules/systemPerformance/ResidualDemand/run_residual_demand.py @@ -357,15 +357,26 @@ def aggregate_delay_results(undamaged_time, damaged_time, od_file_pre, od_file_p compare_df.loc[undamaged_time['agent_id'], 'mean_time_used_undamaged'] = ( undamaged_time['data'].mean(axis=1) ) - compare_df.loc[undamaged_time['agent_id'], 'std_time_used_undamaged'] = ( - undamaged_time['data'].std(axis=1) - ) + compare_df.loc[damaged_time['agent_id'], 'mean_time_used_damaged'] = ( damaged_time['data'].mean(axis=1) ) - compare_df.loc[damaged_time['agent_id'], 'std_time_used_damaged'] = damaged_time[ - 'data' - ].std(axis=1) + + std_time_used_undamaged = np.zeros(len(undamaged_time['agent_id'])) + rows_with_inf = np.any(np.isinf(undamaged_time['data']), axis=1) + std_time_used_undamaged[rows_with_inf] = np.nan + std_time_used_undamaged[~rows_with_inf] = undamaged_time['data'][~rows_with_inf].std( + axis=1 + ) + compare_df.loc[undamaged_time['agent_id'], 'std_time_used_undamaged'] = std_time_used_undamaged + + std_time_used_damaged = np.zeros(len(damaged_time['agent_id'])) + rows_with_inf = np.any(np.isinf(damaged_time['data']), axis=1) + std_time_used_damaged[rows_with_inf] = np.nan + std_time_used_damaged[~rows_with_inf] = damaged_time['data'][~rows_with_inf].std( + axis=1 + ) + compare_df.loc[damaged_time['agent_id'], 'std_time_used_damaged'] = std_time_used_damaged inner_agents = od_df_pre.merge(od_df_post, on='agent_id', how='inner')[ 'agent_id' @@ -385,13 +396,24 @@ def aggregate_delay_results(undamaged_time, damaged_time, od_file_pre, od_file_p - undamaged_time['data'][indices_in_undamaged, :] ) delay_ratio = delay_duration / undamaged_time['data'][indices_in_undamaged, :] + + std_delay_duration = np.zeros(len(inner_agents)) + rows_with_inf = np.any(np.isinf(delay_duration), axis=1) + std_delay_duration[rows_with_inf] = np.nan + std_delay_duration[~rows_with_inf] = delay_duration[~rows_with_inf].std(axis=1) + + std_delay_ratio = np.zeros(len(inner_agents)) + rows_with_inf = np.any(np.isinf(delay_ratio), axis=1) + std_delay_ratio[rows_with_inf] = np.nan + std_delay_ratio[~rows_with_inf] = delay_ratio[~rows_with_inf].std(axis=1) + delay_df = pd.DataFrame( data={ 'agent_id': inner_agents, 'mean_delay_duration': delay_duration.mean(axis=1), 'mean_delay_ratio': delay_ratio.mean(axis=1), - 'std_delay_duration': delay_duration.std(axis=1), - 'std_delay_ratio': delay_ratio.std(axis=1), + 'std_delay_duration': std_delay_duration, + 'std_delay_ratio': std_delay_ratio, } ) @@ -770,6 +792,7 @@ def run_one_realization( trip_info_compare['delay_duration'] / trip_info_compare['travel_time_used_undamaged'] ) + trip_info_compare = trip_info_compare.replace([np.inf, -np.inf], 'inf') trip_info_compare.to_csv('trip_info_compare.csv', index=False) return True diff --git a/modules/systemPerformance/ResidualDemand/transportation.py b/modules/systemPerformance/ResidualDemand/transportation.py index 8f54702ca..933e7c15f 100644 --- a/modules/systemPerformance/ResidualDemand/transportation.py +++ b/modules/systemPerformance/ResidualDemand/transportation.py @@ -393,577 +393,538 @@ def get_graph_network(self, csv_file_dir) -> None: # noqa: D102 nodes_file = os.path.join(csv_file_dir, self.csv_files['network_nodes']) # noqa: PTH118 nodes_gpd.to_csv(nodes_file, index=False) return # noqa: PLR1711 + def substep_assignment( # noqa: PLR0913 + self, + nodes_df=None, + weighted_edges_df=None, + od_ss=None, + quarter_demand=None, + assigned_demand=None, + quarter_counts=4, + trip_info=None, + agent_time_limit=0, + sample_interval=1, + agents_path=None, + hour=None, + quarter=None, + ss_id=None, + alpha_f=0.3, + beta_f=3, + two_way_edges=False, # noqa: FBT002 + ): + # open_edges_df = weighted_edges_df.loc[weighted_edges_df['fft'] < 36000] # noqa: PLR2004 + open_edges_df = weighted_edges_df + + net = pdna.Network( + nodes_df['x'], + nodes_df['y'], + open_edges_df['start_nid'], + open_edges_df['end_nid'], + open_edges_df[['weight']], + twoway=two_way_edges, + ) - # @abstractmethod - def system_performance( # noqa: C901, D102, PLR0915 - self, state # noqa: ARG002 - ) -> None: # Move the CSV creation here - def substep_assignment( # noqa: PLR0913 - nodes_df=None, - weighted_edges_df=None, - od_ss=None, - quarter_demand=None, - assigned_demand=None, - quarter_counts=4, - trip_info=None, - agent_time_limit=0, - sample_interval=1, - agents_path=None, - hour=None, - quarter=None, - ss_id=None, - alpha_f=0.3, - beta_f=3, - two_way_edges=False, # noqa: FBT002 - ): - open_edges_df = weighted_edges_df.loc[weighted_edges_df['fft'] < 36000] # noqa: PLR2004 - - net = pdna.Network( - nodes_df['x'], - nodes_df['y'], - open_edges_df['start_nid'], - open_edges_df['end_nid'], - open_edges_df[['weight']], - twoway=two_way_edges, - ) - - # print('network') # noqa: RUF100, T201 - net.set(pd.Series(net.node_ids)) - # print('net') # noqa: RUF100, T201 - - nodes_origin = od_ss['origin_nid'].to_numpy() - nodes_destin = od_ss['destin_nid'].to_numpy() - nodes_current = od_ss['current_nid'].to_numpy() - agent_ids = od_ss['agent_id'].to_numpy() - agent_current_links = od_ss['current_link'].to_numpy() - agent_current_link_times = od_ss['current_link_time'].to_numpy() - paths = net.shortest_paths(nodes_current, nodes_destin) - - # check agent time limit - path_lengths = net.shortest_path_lengths(nodes_current, nodes_destin) - remove_agent_list = [] - if agent_time_limit is None: - pass - else: - for agent_idx, agent_id in enumerate(agent_ids): - planned_trip_length = path_lengths[agent_idx] - # agent_time_limit[agent_id] - trip_length_limit = agent_time_limit - if planned_trip_length > trip_length_limit + 0: - remove_agent_list.append(agent_id) - - edge_travel_time_dict = weighted_edges_df['t_avg'].T.to_dict() - edge_current_vehicles = weighted_edges_df['veh_current'].T.to_dict() - edge_quarter_vol = weighted_edges_df['vol_true'].T.to_dict() - # edge_length_dict = weighted_edges_df['length'].T.to_dict() - od_residual_ss_list = [] - # all_paths = [] - path_i = 0 - for path in paths: - trip_origin = nodes_origin[path_i] - trip_destin = nodes_destin[path_i] - agent_id = agent_ids[path_i] - # remove some agent (path too long) - if agent_id in remove_agent_list: - path_i += 1 - # no need to update trip info - continue - remaining_time = ( - 3600 / quarter_counts + agent_current_link_times[path_i] - ) - used_time = 0 - for edge_s, edge_e in zip(path, path[1:]): - edge_str = f'{edge_s}-{edge_e}' - edge_travel_time = edge_travel_time_dict[edge_str] - - if (remaining_time > edge_travel_time) and ( - edge_travel_time < 36000 # noqa: PLR2004 - ): - # all_paths.append(edge_str) - # p_dist += edge_travel_time - remaining_time -= edge_travel_time - used_time += edge_travel_time - edge_quarter_vol[edge_str] += 1 * sample_interval - trip_stop = edge_e - - if edge_str == agent_current_links[path_i]: - edge_current_vehicles[edge_str] -= 1 * sample_interval - else: - if edge_str != agent_current_links[path_i]: - edge_current_vehicles[edge_str] += 1 * sample_interval - new_current_link = edge_str - new_current_link_time = remaining_time - trip_stop = edge_s - od_residual_ss_list.append( - [ - agent_id, - trip_origin, - trip_destin, - trip_stop, - new_current_link, - new_current_link_time, - ] - ) - break - trip_info[(agent_id, trip_origin, trip_destin)][0] += ( - 3600 / quarter_counts - ) - trip_info[(agent_id, trip_origin, trip_destin)][1] += used_time - trip_info[(agent_id, trip_origin, trip_destin)][2] = trip_stop - trip_info[(agent_id, trip_origin, trip_destin)][3] = hour - trip_info[(agent_id, trip_origin, trip_destin)][4] = quarter - trip_info[(agent_id, trip_origin, trip_destin)][5] = ss_id + # print('network') # noqa: RUF100, T201 + net.set(pd.Series(net.node_ids)) + # print('net') # noqa: RUF100, T201 + + nodes_origin = od_ss['origin_nid'].to_numpy() + nodes_destin = od_ss['destin_nid'].to_numpy() + nodes_current = od_ss['current_nid'].to_numpy() + agent_ids = od_ss['agent_id'].to_numpy() + agent_current_links = od_ss['current_link'].to_numpy() + agent_current_link_times = od_ss['current_link_time'].to_numpy() + paths = net.shortest_paths(nodes_current, nodes_destin) + + # check agent time limit + path_lengths = net.shortest_path_lengths(nodes_current, nodes_destin) + remove_agent_list = [] + if agent_time_limit is None: + pass + else: + for agent_idx, agent_id in enumerate(agent_ids): + planned_trip_length = path_lengths[agent_idx] + # agent_time_limit[agent_id] + trip_length_limit = agent_time_limit + if planned_trip_length > trip_length_limit + 0: + remove_agent_list.append(agent_id) + + edge_travel_time_dict = weighted_edges_df['t_avg'].T.to_dict() + edge_current_vehicles = weighted_edges_df['veh_current'].T.to_dict() + edge_quarter_vol = weighted_edges_df['vol_true'].T.to_dict() + # edge_length_dict = weighted_edges_df['length'].T.to_dict() + od_residual_ss_list = [] + # all_paths = [] + path_i = 0 + for path in paths: + trip_origin = nodes_origin[path_i] + trip_destin = nodes_destin[path_i] + agent_id = agent_ids[path_i] + # remove some agent (path too long) + if agent_id in remove_agent_list: path_i += 1 - - new_edges_df = weighted_edges_df[ - [ - 'uniqueid', - 'u', - 'v', - 'start_nid', - 'end_nid', - 'fft', - 'capacity', - 'normal_fft', - 'normal_capacity', - 'length', - 'vol_true', - 'vol_tot', - 'veh_current', - 'geometry', - ] - ].copy() - # new_edges_df = new_edges_df.join(edge_volume, how='left') - # new_edges_df['vol_ss'] = new_edges_df['vol_ss'].fillna(0) - # new_edges_df['vol_true'] += new_edges_df['vol_ss'] - new_edges_df['vol_true'] = new_edges_df.index.map(edge_quarter_vol) - new_edges_df['veh_current'] = new_edges_df.index.map( - edge_current_vehicles - ) - # new_edges_df['vol_tot'] += new_edges_df['vol_ss'] - new_edges_df['flow'] = ( - new_edges_df['vol_true'] * quarter_demand / assigned_demand - ) * quarter_counts - new_edges_df['t_avg'] = new_edges_df['fft'] * ( - 1 - + alpha_f - * (new_edges_df['flow'] / new_edges_df['capacity']) ** beta_f + # no need to update trip info + continue + remaining_time = ( + 3600 / quarter_counts + agent_current_link_times[path_i] ) - new_edges_df['t_avg'] = np.where( - new_edges_df['t_avg'] > 36000, 36000, new_edges_df['t_avg'] # noqa: PLR2004 - ) - new_edges_df['t_avg'] = new_edges_df['t_avg'].round(2) - - return new_edges_df, od_residual_ss_list, trip_info, agents_path - - def write_edge_vol( - edges_df=None, - simulation_outputs=None, - quarter=None, - hour=None, - scen_nm=None, - ): - if 'flow' in edges_df.columns: - if edges_df.shape[0] < 10: # noqa: PLR2004 - edges_df[ + used_time = 0 + for edge_s, edge_e in zip(path, path[1:]): + edge_str = f'{edge_s}-{edge_e}' + edge_travel_time = edge_travel_time_dict[edge_str] + if (remaining_time > edge_travel_time) and ( + edge_travel_time < 36000 # noqa: PLR2004 + ): + # all_paths.append(edge_str) + # p_dist += edge_travel_time + remaining_time -= edge_travel_time + used_time += edge_travel_time + edge_quarter_vol[edge_str] += 1 * sample_interval + trip_stop = edge_e + if edge_str == agent_current_links[path_i]: + edge_current_vehicles[edge_str] -= 1 * sample_interval + else: + if edge_str != agent_current_links[path_i]: + edge_current_vehicles[edge_str] += 1 * sample_interval + new_current_link = edge_str + new_current_link_time = remaining_time + trip_stop = edge_s + od_residual_ss_list.append( [ - 'uniqueid', - 'start_nid', - 'end_nid', - 'capacity', - 'veh_current', - 'vol_true', - 'vol_tot', - 'flow', - 't_avg', - 'geometry', + agent_id, + trip_origin, + trip_destin, + trip_stop, + new_current_link, + new_current_link_time, ] - ].to_csv( - f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_' - f'qt{quarter}_{scen_nm}.csv', - index=False, ) + break + trip_info[(agent_id, trip_origin, trip_destin)][0] += ( + 3600 / quarter_counts + ) + trip_info[(agent_id, trip_origin, trip_destin)][1] += used_time + trip_info[(agent_id, trip_origin, trip_destin)][2] = trip_stop + trip_info[(agent_id, trip_origin, trip_destin)][3] = hour + trip_info[(agent_id, trip_origin, trip_destin)][4] = quarter + trip_info[(agent_id, trip_origin, trip_destin)][5] = ss_id + path_i += 1 + + new_edges_df = weighted_edges_df[ + [ + 'uniqueid', + 'u', + 'v', + 'start_nid', + 'end_nid', + 'fft', + 'capacity', + 'normal_fft', + 'normal_capacity', + 'length', + 'vol_true', + 'vol_tot', + 'veh_current', + 'geometry', + ] + ].copy() + # new_edges_df = new_edges_df.join(edge_volume, how='left') + # new_edges_df['vol_ss'] = new_edges_df['vol_ss'].fillna(0) + # new_edges_df['vol_true'] += new_edges_df['vol_ss'] + new_edges_df['vol_true'] = new_edges_df.index.map(edge_quarter_vol) + new_edges_df['veh_current'] = new_edges_df.index.map( + edge_current_vehicles + ) + # new_edges_df['vol_tot'] += new_edges_df['vol_ss'] + new_edges_df['flow'] = ( + new_edges_df['vol_true'] * quarter_demand / assigned_demand + ) * quarter_counts + new_edges_df['t_avg'] = new_edges_df['fft'] * ( + 1 + + alpha_f + * (new_edges_df['flow'] / new_edges_df['capacity']) ** beta_f + ) + new_edges_df['t_avg'] = np.where( + new_edges_df['t_avg'] > 36000, 36000, new_edges_df['t_avg'] # noqa: PLR2004 + ) + new_edges_df['t_avg'] = new_edges_df['t_avg'].round(2) + return new_edges_df, od_residual_ss_list, trip_info, agents_path + def write_edge_vol( + self, + edges_df=None, + simulation_outputs=None, + quarter=None, + hour=None, + scen_nm=None, + ): + if 'flow' in edges_df.columns: + if edges_df.shape[0] < 10: # noqa: PLR2004 + edges_df[ + [ + 'uniqueid', + 'start_nid', + 'end_nid', + 'capacity', + 'veh_current', + 'vol_true', + 'vol_tot', + 'flow', + 't_avg', + 'geometry', + ] + ].to_csv( + f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_' + f'qt{quarter}_{scen_nm}.csv', + index=False, + ) + else: + edges_df.loc[ + edges_df['vol_true'] > 0, + [ + 'uniqueid', + 'start_nid', + 'end_nid', + 'capacity', + 'veh_current', + 'vol_true', + 'vol_tot', + 'flow', + 't_avg', + 'geometry', + ], + ].to_csv( + f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_' + f'qt{quarter}_{scen_nm}.csv', + index=False, + ) + def write_final_vol( + self, + edges_df=None, + simulation_outputs=None, + quarter=None, + hour=None, + scen_nm=None, + ): + edges_df.loc[ + edges_df['vol_tot'] > 0, + ['uniqueid', 'start_nid', 'end_nid', 'vol_tot', 'geometry'], + ].to_csv( + simulation_outputs / 'edge_vol' / f'final_edge_vol_hr{hour}_qt' + f'{quarter}_{scen_nm}.csv', + index=False, + ) + def assignment( # noqa: C901, PLR0913 + self, + quarter_counts=6, + substep_counts=15, + substep_size=30000, + edges_df=None, + nodes_df=None, + od_all=None, + simulation_outputs=None, + scen_nm=None, + hour_list=None, + quarter_list=None, + cost_factor=None, # noqa: ARG001 + closure_hours=None, + closed_links=None, + agent_time_limit=None, + sample_interval=1, + agents_path=None, + alpha_f=0.3, + beta_f=4, + two_way_edges=False, # noqa: FBT002 + save_edge_vol=True, + ): + if closure_hours is None: + closure_hours = [] + + # Check if all od pairs has path. If not, + orig = od_all['origin_nid'].to_numpy() + dest = od_all['destin_nid'].to_numpy() + open_edges_df = edges_df[ + ~edges_df['uniqueid'].isin(closed_links['uniqueid'].to_numpy()) + ] + net = pdna.Network( + nodes_df['x'], + nodes_df['y'], + open_edges_df['start_nid'], + open_edges_df['end_nid'], + open_edges_df[['fft']], + twoway=two_way_edges, + ) + net.set(pd.Series(net.node_ids)) + paths = net.shortest_paths(orig, dest) + no_path_ind = [i for i in range(len(paths)) if len(paths[i]) == 0] + od_no_path = od_all.iloc[no_path_ind].copy() + od_all = od_all.drop(no_path_ind) + + od_all['current_nid'] = od_all['origin_nid'] + trip_info = { + (od.agent_id, od.origin_nid, od.destin_nid): [ + 0, + 0, + od.origin_nid, + 0, + od.hour, + od.quarter, + 0, + 0, + ] + for od in od_all.itertuples() + } + # Quarters and substeps + # probability of being in each division of hour + if quarter_list is None: + quarter_counts = 4 + else: + quarter_counts = len(quarter_list) + quarter_ps = [1 / quarter_counts for i in range(quarter_counts)] + quarter_ids = list(range(quarter_counts)) + + # initial setup + od_residual_list = [] + # accumulator + edges_df['vol_tot'] = 0 + edges_df['veh_current'] = 0 + + # Loop through days and hours + for _day in ['na']: + for hour in hour_list: + gc.collect() + if hour in closure_hours: + for row in closed_links.itertuples(): + edges_df.loc[ + (edges_df['uniqueid'] == row.uniqueid), 'capacity' + ] = 1 + edges_df.loc[ + (edges_df['uniqueid'] == row.uniqueid), 'fft' + ] = 36000 + else: + edges_df['capacity'] = edges_df['normal_capacity'] + edges_df['fft'] = edges_df['normal_fft'] + + # Read OD + od_hour = od_all[od_all['hour'] == hour].copy() + if od_hour.shape[0] == 0: + od_hour = pd.DataFrame([], columns=od_all.columns) + od_hour['current_link'] = None + od_hour['current_link_time'] = 0 + + # Divide into quarters + if 'quarter' in od_all.columns: + pass else: - edges_df.loc[ - edges_df['vol_true'] > 0, + od_quarter_msk = np.random.choice( + quarter_ids, size=od_hour.shape[0], p=quarter_ps + ) + od_hour['quarter'] = od_quarter_msk + + if quarter_list is None: + quarter_list = quarter_ids + for quarter in quarter_list: + # New OD in assignment period + od_quarter = od_hour.loc[ + od_hour['quarter'] == quarter, [ - 'uniqueid', - 'start_nid', - 'end_nid', - 'capacity', - 'veh_current', - 'vol_true', - 'vol_tot', - 'flow', - 't_avg', - 'geometry', + 'agent_id', + 'origin_nid', + 'destin_nid', + 'current_nid', + 'current_link', + 'current_link_time', + ], + ] + # Add resudal OD + od_residual = pd.DataFrame( + od_residual_list, + columns=[ + 'agent_id', + 'origin_nid', + 'destin_nid', + 'current_nid', + 'current_link', + 'current_link_time', ], - ].to_csv( - f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_' - f'qt{quarter}_{scen_nm}.csv', - index=False, ) + od_residual['quarter'] = quarter + # Total OD in each assignment period is the combined + # of new and residual OD: + od_quarter = pd.concat( + [od_quarter, od_residual], sort=False, ignore_index=True + ) + # Residual OD is no longer residual after it has been + # merged to the quarterly OD: + od_residual_list = [] + od_quarter = od_quarter[ + od_quarter['current_nid'] != od_quarter['destin_nid'] + ] + # total demand for this quarter, including total and + # residual demand: + quarter_demand = od_quarter.shape[0] + # how many among the OD pairs to be assigned in this + # quarter are actually residual from previous quarters + residual_demand = od_residual.shape[0] + assigned_demand = 0 + + substep_counts = max(1, (quarter_demand // substep_size) + 1) + logging.info( + f'HR {hour} QT {quarter} has ' + f'{quarter_demand}/{residual_demand} od' + f's/residuals {substep_counts} substeps' + ) + substep_ps = [ + 1 / substep_counts for i in range(substep_counts) + ] + substep_ids = list(range(substep_counts)) + od_substep_msk = np.random.choice( + substep_ids, size=quarter_demand, p=substep_ps + ) + od_quarter['ss_id'] = od_substep_msk + # reset volume at each quarter + edges_df['vol_true'] = 0 - def write_final_vol( - edges_df=None, - simulation_outputs=None, - quarter=None, - hour=None, - scen_nm=None, - ): - edges_df.loc[ - edges_df['vol_tot'] > 0, - ['uniqueid', 'start_nid', 'end_nid', 'vol_tot', 'geometry'], - ].to_csv( - simulation_outputs / 'edge_vol' / f'final_edge_vol_hr{hour}_qt' - f'{quarter}_{scen_nm}.csv', - index=False, - ) - - def assignment( # noqa: C901, PLR0913 - quarter_counts=6, - substep_counts=15, - substep_size=30000, - edges_df=None, - nodes_df=None, - od_all=None, - simulation_outputs=None, - scen_nm=None, - hour_list=None, - quarter_list=None, - cost_factor=None, # noqa: ARG001 - closure_hours=None, - closed_links=None, - agent_time_limit=None, - sample_interval=1, - agents_path=None, - alpha_f=0.3, - beta_f=4, - two_way_edges=False, # noqa: FBT002 - ): - if closure_hours is None: - closure_hours = [] - - # Check if all od pairs has path. If not, - orig = od_all['origin_nid'].to_numpy() - dest = od_all['destin_nid'].to_numpy() - open_edges_df = edges_df[ - ~edges_df['uniqueid'].isin(closed_links['uniqueid'].to_numpy()) - ] - net = pdna.Network( - nodes_df['x'], - nodes_df['y'], - open_edges_df['start_nid'], - open_edges_df['end_nid'], - open_edges_df[['fft']], - twoway=two_way_edges, - ) - paths = net.shortest_paths(orig, dest) - no_path_ind = [i for i in range(len(paths)) if len(paths[i]) == 0] - od_no_path = od_all.iloc[no_path_ind].copy() - od_all = od_all.drop(no_path_ind) - - od_all['current_nid'] = od_all['origin_nid'] - trip_info = { - (od.agent_id, od.origin_nid, od.destin_nid): [ - 0, - 0, - od.origin_nid, - 0, - od.hour, - od.quarter, - 0, - 0, - ] - for od in od_all.itertuples() - } - - # Quarters and substeps - # probability of being in each division of hour - if quarter_list is None: - quarter_counts = 4 - else: - quarter_counts = len(quarter_list) - quarter_ps = [1 / quarter_counts for i in range(quarter_counts)] - quarter_ids = list(range(quarter_counts)) - - # initial setup - od_residual_list = [] - # accumulator - edges_df['vol_tot'] = 0 - edges_df['veh_current'] = 0 - - # Loop through days and hours - for _day in ['na']: - for hour in hour_list: - gc.collect() - if hour in closure_hours: - for row in closed_links.itertuples(): - edges_df.loc[ - (edges_df['uniqueid'] == row.uniqueid), 'capacity' - ] = 1 - edges_df.loc[ - (edges_df['uniqueid'] == row.uniqueid), 'fft' - ] = 36000 - else: - edges_df['capacity'] = edges_df['normal_capacity'] - edges_df['fft'] = edges_df['normal_fft'] - - # Read OD - od_hour = od_all[od_all['hour'] == hour].copy() - if od_hour.shape[0] == 0: - od_hour = pd.DataFrame([], columns=od_all.columns) - od_hour['current_link'] = None - od_hour['current_link_time'] = 0 - - # Divide into quarters - if 'quarter' in od_all.columns: - pass - else: - od_quarter_msk = np.random.choice( - quarter_ids, size=od_hour.shape[0], p=quarter_ps - ) - od_hour['quarter'] = od_quarter_msk - - if quarter_list is None: - quarter_list = quarter_ids - for quarter in quarter_list: - # New OD in assignment period - od_quarter = od_hour.loc[ - od_hour['quarter'] == quarter, - [ - 'agent_id', - 'origin_nid', - 'destin_nid', - 'current_nid', - 'current_link', - 'current_link_time', - ], - ] - # Add resudal OD - od_residual = pd.DataFrame( - od_residual_list, - columns=[ - 'agent_id', - 'origin_nid', - 'destin_nid', - 'current_nid', - 'current_link', - 'current_link_time', - ], - ) - od_residual['quarter'] = quarter - # Total OD in each assignment period is the combined - # of new and residual OD: - od_quarter = pd.concat( - [od_quarter, od_residual], sort=False, ignore_index=True - ) - # Residual OD is no longer residual after it has been - # merged to the quarterly OD: - od_residual_list = [] - od_quarter = od_quarter[ - od_quarter['current_nid'] != od_quarter['destin_nid'] - ] - - # total demand for this quarter, including total and - # residual demand: - quarter_demand = od_quarter.shape[0] - # how many among the OD pairs to be assigned in this - # quarter are actually residual from previous quarters - residual_demand = od_residual.shape[0] - assigned_demand = 0 - - substep_counts = max(1, (quarter_demand // substep_size) + 1) + for ss_id in substep_ids: + gc.collect() + time_ss_0 = time.time() logging.info( - f'HR {hour} QT {quarter} has ' - f'{quarter_demand}/{residual_demand} od' - f's/residuals {substep_counts} substeps' + f'Hour: {hour}, Quarter: {quarter}, ' + 'SS ID: {ss_id}' ) - substep_ps = [ - 1 / substep_counts for i in range(substep_counts) - ] - substep_ids = list(range(substep_counts)) - od_substep_msk = np.random.choice( - substep_ids, size=quarter_demand, p=substep_ps + od_ss = od_quarter[od_quarter['ss_id'] == ss_id] + assigned_demand += od_ss.shape[0] + if assigned_demand == 0: + continue + # calculate weight + weighted_edges_df = edges_df.copy() + # weight by travel distance + # weighted_edges_df['weight'] = edges_df['length'] + # weight by travel time + # weighted_edges_df['weight'] = edges_df['t_avg'] + # weight by travel time + # weighted_edges_df['weight'] = (edges_df['t_avg'] + # - edges_df['fft']) * 0.5 + edges_df['length']*0.1 + # + cost_factor*edges_df['length']*0.1*( + # edges_df['is_highway']) + # 10 yen per 100 m --> 0.1 yen per m + weighted_edges_df['weight'] = edges_df['t_avg'] + # weighted_edges_df['weight'] = np.where( + # weighted_edges_df['weight']<0.1, 0.1, + # weighted_edges_df['weight']) + # traffic assignment with truncated path + ( + edges_df, + od_residual_ss_list, + trip_info, + agents_path, + ) = self.substep_assignment( + nodes_df=nodes_df, + weighted_edges_df=weighted_edges_df, + od_ss=od_ss, + quarter_demand=quarter_demand, + assigned_demand=assigned_demand, + quarter_counts=quarter_counts, + trip_info=trip_info, + agent_time_limit=agent_time_limit, + sample_interval=sample_interval, + agents_path=agents_path, + hour=hour, + quarter=quarter, + ss_id=ss_id, + alpha_f=alpha_f, + beta_f=beta_f, + two_way_edges=two_way_edges, ) - od_quarter['ss_id'] = od_substep_msk - - # reset volume at each quarter - edges_df['vol_true'] = 0 - - for ss_id in substep_ids: - gc.collect() - - time_ss_0 = time.time() - logging.info( - f'Hour: {hour}, Quarter: {quarter}, ' - 'SS ID: {ss_id}' - ) - od_ss = od_quarter[od_quarter['ss_id'] == ss_id] - assigned_demand += od_ss.shape[0] - if assigned_demand == 0: - continue - # calculate weight - weighted_edges_df = edges_df.copy() - # weight by travel distance - # weighted_edges_df['weight'] = edges_df['length'] - # weight by travel time - # weighted_edges_df['weight'] = edges_df['t_avg'] - # weight by travel time - # weighted_edges_df['weight'] = (edges_df['t_avg'] - # - edges_df['fft']) * 0.5 + edges_df['length']*0.1 - # + cost_factor*edges_df['length']*0.1*( - # edges_df['is_highway']) - # 10 yen per 100 m --> 0.1 yen per m - weighted_edges_df['weight'] = edges_df['t_avg'] - # weighted_edges_df['weight'] = np.where( - # weighted_edges_df['weight']<0.1, 0.1, - # weighted_edges_df['weight']) - - # traffic assignment with truncated path - ( - edges_df, - od_residual_ss_list, - trip_info, - agents_path, - ) = substep_assignment( - nodes_df=nodes_df, - weighted_edges_df=weighted_edges_df, - od_ss=od_ss, - quarter_demand=quarter_demand, - assigned_demand=assigned_demand, - quarter_counts=quarter_counts, - trip_info=trip_info, - agent_time_limit=agent_time_limit, - sample_interval=sample_interval, - agents_path=agents_path, - hour=hour, - quarter=quarter, - ss_id=ss_id, - alpha_f=alpha_f, - beta_f=beta_f, - two_way_edges=two_way_edges, - ) - - od_residual_list += od_residual_ss_list - # write_edge_vol(edges_df=edges_df, - # simulation_outputs=simulation_outputs, - # quarter=quarter, - # hour=hour, - # scen_nm='ss{}_{}'.format(ss_id, scen_nm)) - logging.info( - f'HR {hour} QT {quarter} SS {ss_id}' - ' finished, max vol ' - f'{np.max(edges_df["vol_true"])}, ' - f'time {time.time() - time_ss_0}' - ) - - # write quarterly results - edges_df['vol_tot'] += edges_df['vol_true'] - if True: # hour >=16 or (hour==15 and quarter==3): - write_edge_vol( - edges_df=edges_df, - simulation_outputs=simulation_outputs, - quarter=quarter, - hour=hour, - scen_nm=scen_nm, - ) - - if hour % 3 == 0: - trip_info_df = pd.DataFrame( - [ - [ - trip_key[0], - trip_key[1], - trip_key[2], - trip_value[0], - trip_value[1], - trip_value[2], - trip_value[3], - trip_value[4], - trip_value[5], - ] - for trip_key, trip_value in trip_info.items() - ], - columns=[ - 'agent_id', - 'origin_nid', - 'destin_nid', - 'travel_time', - 'travel_time_used', - 'stop_nid', - 'stop_hour', - 'stop_quarter', - 'stop_ssid', - ], - ) - trip_info_df.to_csv( - simulation_outputs / 'trip_info' - f'/trip_info_{scen_nm}_hr{hour}' - '.csv', - index=False, + od_residual_list += od_residual_ss_list + # write_edge_vol(edges_df=edges_df, + # simulation_outputs=simulation_outputs, + # quarter=quarter, + # hour=hour, + # scen_nm='ss{}_{}'.format(ss_id, scen_nm)) + logging.info( + f'HR {hour} QT {quarter} SS {ss_id}' + ' finished, max vol ' + f'{np.max(edges_df["vol_true"])}, ' + f'time {time.time() - time_ss_0}' ) - # output individual trip travel time and stop location + # write quarterly results + edges_df['vol_tot'] += edges_df['vol_true'] + if save_edge_vol: # hour >=16 or (hour==15 and quarter==3): + self.write_edge_vol( + edges_df=edges_df, + simulation_outputs=simulation_outputs, + quarter=quarter, + hour=hour, + scen_nm=scen_nm, + ) - trip_info_df = pd.DataFrame( + trip_info_df = pd.DataFrame( + [ [ - [ - trip_key[0], - trip_key[1], - trip_key[2], - trip_value[0], - trip_value[1], - trip_value[2], - trip_value[3], - trip_value[4], - trip_value[5], - ] - for trip_key, trip_value in trip_info.items() - ], - columns=[ - 'agent_id', - 'origin_nid', - 'destin_nid', - 'travel_time', - 'travel_time_used', - 'stop_nid', - 'stop_hour', - 'stop_quarter', - 'stop_ssid', - ], - ) - # Add the no path OD to the trip info - trip_info_no_path = od_no_path.drop( - columns=[ - col - for col in od_no_path.columns - if col not in ['agent_id', 'origin_nid', 'destin_nid'] + trip_key[0], + trip_key[1], + trip_key[2], + trip_value[0], + trip_value[1], + trip_value[2], + trip_value[3], + trip_value[4], + trip_value[5], ] - ) - trip_info_no_path['travel_time'] = 360000 - trip_info_no_path['travel_time_used'] = np.nan - trip_info_no_path['stop_nid'] = np.nan - trip_info_no_path['stop_hour'] = np.nan - trip_info_no_path['stop_quarter'] = np.nan - trip_info_no_path['stop_ssid'] = np.nan - trip_info_df = pd.concat( - [trip_info_df, trip_info_no_path], ignore_index=True - ) - + for trip_key, trip_value in trip_info.items() + ], + columns=[ + 'agent_id', + 'origin_nid', + 'destin_nid', + 'travel_time', + 'travel_time_used', + 'stop_nid', + 'stop_hour', + 'stop_quarter', + 'stop_ssid', + ], + ) + # If there are trips incompleted mark them as np.nan + incomplete_trips_agent_id = [x[0] for x in od_residual_list] + trip_info_df.loc[trip_info_df.agent_id.isin(incomplete_trips_agent_id),'travel_time_used'] = 'inf' + # Add the no path OD to the trip info + trip_info_no_path = od_no_path.drop( + columns=[ + col + for col in od_no_path.columns + if col not in ['agent_id', 'origin_nid', 'destin_nid'] + ] + ) + trip_info_no_path['travel_time'] = 360000 + trip_info_no_path['travel_time_used'] = 'inf' + trip_info_no_path['stop_nid'] = np.nan + trip_info_no_path['stop_hour'] = np.nan + trip_info_no_path['stop_quarter'] = np.nan + trip_info_no_path['stop_ssid'] = np.nan + trip_info_df = pd.concat( + [trip_info_df, trip_info_no_path], ignore_index=True + ) + if save_edge_vol: trip_info_df.to_csv( simulation_outputs / 'trip_info' / f'trip_info_{scen_nm}.csv', index=False, ) - write_final_vol( - edges_df=edges_df, - simulation_outputs=simulation_outputs, - quarter=quarter, - hour=hour, - scen_nm=scen_nm, - ) + self.write_final_vol( + edges_df=edges_df, + simulation_outputs=simulation_outputs, + quarter=quarter, + hour=hour, + scen_nm=scen_nm, + ) + return trip_info_df + # @abstractmethod + def system_performance( # noqa: C901, D102, PLR0915 + self, state # noqa: ARG002 + ) -> None: # Move the CSV creation here network_edges = self.csv_files['network_edges'] network_nodes = self.csv_files['network_nodes'] @@ -1051,7 +1012,7 @@ def assignment( # noqa: C901, PLR0913 t_od_1 = time.time() logging.info('%d sec to read %d OD pairs', t_od_1 - t_od_0, od_all.shape[0]) # run residual_demand_assignment - assignment( + self.assignment( edges_df=edges_df, nodes_df=nodes_df, od_all=od_all, @@ -1251,6 +1212,110 @@ def update_nodes_file(nodes, det): ) od_df.to_csv('updated_od.csv') return od_df + + +class pyrecodes_residual_demand(TransportationPerformance): + def __init__( + self, + edges_file, + nodes_file, + od_pre_file, + hour_list, + results_dir, + two_way_edges=False + ): + # Prepare edges and nodes files + edges_gdf = gpd.read_file(edges_file).to_crs(epsg=6500) + if 'length' not in edges_gdf.columns: + edges_gdf['length'] = edges_gdf['geometry'].apply(lambda x: x.length) + edges_gdf = edges_gdf.to_crs(epsg=4326) + if two_way_edges: + edges_gdf_copy = edges_gdf.copy() + edges_gdf_copy['StartNode'] = edges_gdf['EndNode'] + edges_gdf_copy['EndNode'] = edges_gdf['StartNode'] + edges_gdf = pd.concat([edges_gdf, edges_gdf_copy], ignore_index=True) + edges_gdf = edges_gdf.reset_index() + edges_gdf = edges_gdf.rename( + columns={ + 'index': 'uniqueid', + 'StartNode': 'start_nid', + 'EndNode': 'end_nid', + 'NumOfLanes': 'lanes', + 'MaxMPH': 'maxspeed', + } + ) + # Assume that the capacity for each lane is 1800 + edges_gdf['capacity'] = edges_gdf['lanes'] * 1800 + edges_gdf['normal_capacity'] = edges_gdf['capacity'] + edges_gdf['normal_maxspeed'] = edges_gdf['maxspeed'] + edges_gdf['start_nid'] = edges_gdf['start_nid'].astype(int) + edges_gdf['end_nid'] = edges_gdf['end_nid'].astype(int) + + nodes_gdf = gpd.read_file(nodes_file) + nodes_gdf = nodes_gdf.rename(columns={'nodeID': 'node_id'}) + nodes_gdf['y'] = nodes_gdf['geometry'].apply(lambda x: x.y) + nodes_gdf['x'] = nodes_gdf['geometry'].apply(lambda x: x.x) + + edges_gdf['fft'] = edges_gdf['length'] / edges_gdf['maxspeed'] * 2.23694 + edges_gdf['normal_fft'] = ( + edges_gdf['length'] / edges_gdf['normal_maxspeed'] * 2.23694 + ) + edges_gdf = edges_gdf.sort_values(by='fft', ascending=False).drop_duplicates( + subset=['start_nid', 'end_nid'], keep='first' + ) + edges_gdf['edge_str'] = ( + edges_gdf['start_nid'].astype('str') + + '-' + + edges_gdf['end_nid'].astype('str') + ) + edges_gdf['t_avg'] = edges_gdf['fft'] + edges_gdf['u'] = edges_gdf['start_nid'] + edges_gdf['v'] = edges_gdf['end_nid'] + edges_gdf = edges_gdf.set_index('edge_str') + + # delete geometry columns to save memory however, traffic assignment need to be updated accordingly + # edges_gdf = edges_gdf.drop(columns=['geometry']) + # nodes_gdf = nodes_gdf.drop(columns=['geometry']) + + + self.edges_df = edges_gdf + self.nodes_df = nodes_gdf + + self.od_pre_file = od_pre_file + self.hour_list = hour_list + self.results_dir = results_dir + self.two_way_edges = two_way_edges + # If needed run on undamaged network here + def simulate( + self, + r2d_dict + ): + def update_edges(edges, r2d_dict): + return edges + + od_matrix = self.compute_od(r2d_dict) + + edges_df = update_edges(self.edges_df, r2d_dict) + + # closed_links is simulated as normal links with very high fft + closed_links = pd.DataFrame([], columns=['uniqueid']) + trip_info_df = self.assignment( + edges_df = edges_df, + nodes_df=self.nodes_df, + od_all=od_matrix, + simulation_outputs=self.results_dir, + hour_list=self.hour_list, + quarter_list=[0, 1, 2, 3, 4, 5], + closure_hours=self.hour_list, + closed_links=closed_links, + two_way_edges=self.two_way_edges, + save_edge_vol = False + ) + return trip_info_df + + def compute_od(self, r2d_dict): + od_matrix = pd.read_csv(self.od_pre_file) + return od_matrix # def get_graph_network(self, # det_file_path: str,