jz - Modify Capacity Spectrum Method to copy PGD IMs to response.csv

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',

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