OpenPBEE · dustin-cook · Dec 6, 2025 · Dec 8, 2025 · Dec 10, 2025 · Dec 11, 2025
diff --git a/README.md b/README.md
@@ -1,5 +1,35 @@
 # Functional-Recovery---Python
 This is translation of Matlab codebase into Python for quantifying building-specific functional recovery and reoccupancy based on a probabilistic performance-based earthquake engineering framework.
+
+## Requirements
+- The `requirements.txt` file defines the Python package dependencies required to run this codebase. Follow the instructions below to install all required depenedancies listed in the 'requirements.txt' file.
+- Recommended Python version: `3.9` (the codebase was developed and tested with Python 3.9).
+
+Installation (using a virtual environment is recommended):
+
+```powershell
+# create a virtual environment
+python -m venv .venv
+
+# activate the virtual environment (PowerShell)
+.\.venv\Scripts\Activate.ps1
+
+# upgrade pip (optional but recommended)
+python -m pip install --upgrade pip
+
+# install dependencies from requirements.txt
+pip install -r requirements.txt
+```
+
+If you prefer conda:
+
+```bash
+conda create -n frec python=3.9
+conda activate frec
+pip install -r requirements.txt
+```
+
+If you run into platform-specific dependency issues, please refer to the package error messages and install any missing system libraries before re-running `pip install -r requirements.txt`.
 Original Matlab code is from Dr. Dustin Cook's Github directory https://github.com/OpenPBEE/PBEE-Recovery.
 
 ### Method Description

diff --git a/fn_red_tag.py b/fn_red_tag.py
@@ -1,4 +1,4 @@
-def fn_red_tag( calculate_red_tag, damage, comps, simulated_replacement_time):
+def fn_red_tag( calculate_red_tag, damage, comps, simulated_replacement_time, red_tag_options ):
     '''Perform the ATC-138 functional recovery time assessement given similation
     of component damage for a single shaking intensity
 
@@ -27,18 +27,23 @@ def fn_red_tag( calculate_red_tag, damage, comps, simulated_replacement_time):
 
     import numpy as np
 
-    def simulate_tagging(damage, comps, sc_ids, sc_thresholds):
+    def simulate_tagging(damage, comps, sc_ids, sc_thresholds, red_tag_options ):
 
         red_tag_impact = np.zeros(np.shape(damage['tenant_units'][0]['qnt_damaged'])) # num reals by num comp_ds
+        red_tag_impact_cb = np.zeros(np.shape(damage['tenant_units'][0]['qnt_damaged'])) # num reals by num comp_ds, specific for counting red tags for coupling beams
 
         num_reals = len(red_tag_impact)
         # Go through each structural system and calc the realizations where the
         # building is red tagged
 
         sc_tag = np.zeros([num_reals,len(sc_ids)])
         for sc in range(len(sc_ids)):
-
             story_tag = np.zeros([num_reals,len(damage['story'])])
+
+            # resets the coupling beam counter for every safety class
+            coupling_beam_dmg = np.zeros([num_reals,3]) # three directions
+            coupling_beam_qty = np.zeros([1,3]) # three directions
+            red_tag_impact_cb_sc = np.zeros(np.shape(damage['tenant_units'][0]['qnt_damaged'])) # num reals by num comp_ds
 
             for s in range(len(damage['story'])):
                 sc_filt = np.array(damage['comp_ds_table']['safety_class']) >= sc_ids[sc]
@@ -54,14 +59,18 @@ def simulate_tagging(damage, comps, sc_ids, sc_thresholds):
 
                     sys_tag = np.zeros([num_reals, len(structural_systems)]).astype(bool)
 
+                    if red_tag_options['ignore_coupling_beam_for_red_tag']:
+                          # just ignore system 12 (coupling beams)
+                          structural_systems = structural_systems[structural_systems != 12]
+
                     for sys in range(len(structural_systems)):
                         ss_filt_ds = np.logical_or(np.array(damage['comp_ds_table']['structural_system']) == structural_systems[sys], np.array(damage['comp_ds_table']['structural_system_alt']) == structural_systems[sys])
                         ss_filt_comp = np.logical_or(np.array(comps['comp_table']['structural_system']) == structural_systems[sys], np.array(comps['comp_table']['structural_system_alt']) == structural_systems[sys])
 
                         # Check damage among each series within this structural system
                         series = np.unique(np.array(damage['comp_ds_table']['structural_series_id'])[ss_filt_ds])
                         ser_dmg = np.zeros([num_reals,len(series)])
-                        ser_qty = np.zeros([num_reals,len(series)])
+                        ser_qty = np.zeros([1,len(series)])
                         for ser in range(len(series)):
                             ser_filt_ds = np.array(damage['comp_ds_table']['structural_series_id']) == series[ser] 
                             ser_filt_comp = np.array(comps['comp_table']['structural_series_id']) == series[ser] 
@@ -73,17 +82,24 @@ def simulate_tagging(damage, comps, sc_ids, sc_thresholds):
                             ser_qty[:,ser] = np.sum(num_comps[ser_filt_comp & ss_filt_comp])
 
                         # Check if this system is causing a red tag
-                        sys_dmg = np.nanmax(ser_dmg, axis = 1)
-                        sys_qty = np.nanmax(ser_qty, axis = 1)
-                        sys_ratio = sys_dmg / sys_qty
-                        sys_tag[:,sys] = sys_ratio > sc_thresholds[sc]
-
-                        '''Calculate the impact that each component has on red tag
-                        (boolean, 1 = affects red tag, 0 = does not affect)
-                        Take all damage that is part of this system at this story
-                        in this direction that is damaged to this safety class
-                        level, only where damage exceeds tagging threshold'''
-                        red_tag_impact = np.fmax(red_tag_impact, 1*sys_tag[:,sys].reshape(len(sys_tag),1) * ss_filt_ds.reshape(1,len(ss_filt_ds)) * sc_filt.reshape(1,len(ss_filt_ds)) * (sc_dmg>0))
+                        if structural_systems[sys] == 12 and bool(red_tag_options.get('tag_coupling_beams_over_height', False)):
+                            # HARED CODED CHECK FOR COUPLING BEAMS
+                            # just do counts for now instead of adding to red tag per story check
+                            coupling_beam_dmg[:, direc-1] = coupling_beam_dmg[:,direc-1] + np.nanmax(ser_dmg, axis = 1)
+                            coupling_beam_qty[0, direc-1] = coupling_beam_qty[0,direc-1] + np.nanmax(ser_qty, axis = 1)
+                            red_tag_impact_cb_sc = np.fmax(red_tag_impact_cb_sc, ss_filt_ds.reshape(1,-1) * sc_filt.reshape(1,-1) * (sc_dmg>0))
+                        else:
+                            sys_dmg = np.nanmax(ser_dmg, axis = 1)
+                            sys_qty = np.nanmax(ser_qty, axis = 1)
+                            sys_ratio = sys_dmg / sys_qty
+                            sys_tag[:,sys] = sys_ratio > sc_thresholds[sc]
+
+                            '''Calculate the impact that each component has on red tag
+                            (boolean, 1 = affects red tag, 0 = does not affect)
+                            Take all damage that is part of this system at this story
+                            in this direction that is damaged to this safety class
+                            level, only where damage exceeds tagging threshold'''
+                            red_tag_impact = np.fmax(red_tag_impact, 1*sys_tag[:,sys].reshape(-1,1) * ss_filt_ds.reshape(1,-1) * sc_filt.reshape(1,-1) * (sc_dmg>0))
 
 
                     # Combine across all systems in this direction
@@ -96,7 +112,19 @@ def simulate_tagging(damage, comps, sc_ids, sc_thresholds):
 
             # Combine across all stories for this safety class
             sc_tag[:,sc] = np.nanmax(story_tag, axis = 1)
+
+            # Additional Hard Coded Check for Coupling Beams
+            if np.any(np.array(damage['comp_ds_table']['structural_system']) == 12):
+                CB_damage_ratio = coupling_beam_dmg / coupling_beam_qty
+                CB_tag_sc = np.nanmax(CB_damage_ratio, axis = 1) > sc_thresholds[sc]
+                sc_tag[:,sc] = np.fmax(sc_tag[:,sc], CB_tag_sc) # add it back to sum
+                red_tag_impact_cb = np.fmax(
+                    red_tag_impact_cb,
+                    red_tag_impact_cb_sc * sc_tag[:, sc].reshape(-1, 1)
+                )
 
+        red_tag_impact = np.fmax(red_tag_impact, red_tag_impact_cb)
+
         # Combine all safety class checks into one simulated red tag
         red_tag =  np.nanmax(sc_tag, axis = 1)
 
@@ -116,10 +144,11 @@ def simulate_tagging(damage, comps, sc_ids, sc_thresholds):
         # Simulate Red Tags
         sc_ids = np.array([1, 2, 3, 4])
         sc_thresholds = np.array([0.5, 0.25, 0.1, 0])
-        red_tag, red_tag_impact = simulate_tagging(damage, comps, sc_ids, sc_thresholds)
+        red_tag, red_tag_impact = simulate_tagging(damage, comps, sc_ids, sc_thresholds, red_tag_options )
 
         # Inspection is flagged for 50% of the red tag thresholds
-        inspection_tag = simulate_tagging(damage, comps, sc_ids, 0.5*sc_thresholds)[0]
+        inspection_thresholds = 0.5*sc_thresholds
+        inspection_tag, _ = simulate_tagging(damage, comps, sc_ids, inspection_thresholds, red_tag_options )
 
     else:
 

diff --git a/functionality/fn_calculate_reoccupancy.py b/functionality/fn_calculate_reoccupancy.py
@@ -67,18 +67,24 @@ def fn_calculate_reoccupancy(damage, damage_consequences, utilities,
     recovery_day['tenant_safety'], comp_breakdowns['tenant_safety'] = other_functionality_functions.fn_tenant_safety( damage, building_model, functionality_options, tenant_units)
 
     ## Combine Check to determine the day the each tenant unit is reoccupiable
-    # Check the day the building is safe
-    day_building_safe = np.fmax(recovery_day['building_safety']['red_tag'],
-                        np.fmax(recovery_day['building_safety']['shoring'],
-                        np.fmax(recovery_day['building_safety']['hazardous_material'],
-                        np.fmax(recovery_day['building_safety']['entry_door_access'],
-                        recovery_day['building_safety']['fire_suppression']))))
-    # Check the day each story is accessible
-    day_story_accessible = np.fmax(recovery_day['story_access']['stairs'], np.fmax(recovery_day['story_access']['stair_doors'], np.fmax(recovery_day['story_access']['flooding'], recovery_day['story_access']['horizontal_egress'])))
-
-    # Check the day each tenant unit is safe
-    day_tenant_unit_safe = np.fmax(recovery_day['tenant_safety']['interior'], np.fmax(recovery_day['tenant_safety']['exterior'], recovery_day['tenant_safety']['hazardous_material']))
+    # Go through each of the building safety checks and combine them to check the day the building is safe (max of all checks)
+    fault_tree_events_building_safety = list(recovery_day['building_safety'].keys())
+    day_building_safe = recovery_day['building_safety'][fault_tree_events_building_safety[0]]
+    for i in range(1, len(fault_tree_events_building_safety)):
+        day_building_safe = np.fmax(day_building_safe, recovery_day['building_safety'][fault_tree_events_building_safety[i]])
+
+    # Go through each of the story access checks and combine them to check the day each story is accessible (max of all checks)
+    fault_tree_events_story_access = list(recovery_day['story_access'].keys())
+    day_story_accessible = recovery_day['story_access'][fault_tree_events_story_access[0]]
+    for i in range(1, len(fault_tree_events_story_access)):
+        day_story_accessible = np.fmax(day_story_accessible, recovery_day['story_access'][fault_tree_events_story_access[i]])
 
+    # Go through each of the tenant unit safety checks and combine them to check the day each tenant unit is safe (max of all checks)
+    fault_tree_events_tenant_unit_safe = list(recovery_day['tenant_safety'].keys())
+    day_tenant_unit_safe = recovery_day['tenant_safety'][fault_tree_events_tenant_unit_safe[0]]
+    for i in range(1, len(fault_tree_events_tenant_unit_safe)):
+        day_tenant_unit_safe = np.fmax(day_tenant_unit_safe, recovery_day['tenant_safety'][fault_tree_events_tenant_unit_safe[i]])
+
     # Combine checks to determine when each tenant unit is re-occupiable
     day_tenant_unit_reoccupiable = np.fmax(np.fmax(day_building_safe.reshape(len(day_building_safe),1), day_story_accessible), day_tenant_unit_safe)
 

diff --git a/functionality/fn_check_habitability.py b/functionality/fn_check_habitability.py
@@ -41,9 +41,9 @@ def fn_check_habitability( damage, damage_consequences, reoc_meta, func_meta,
 
     # Go through each of the tenant function branches and combines checks
     day_tenant_unit_reoccupiable = np.zeros([num_reals,1])
-    fault_tree_events_LV1 = list(comp_breakdowns.keys())
+    fault_tree_events_LV1 = list(recovery_day.keys())
     for i in range(len(fault_tree_events_LV1)):
-        fault_tree_events_LV2 = list(comp_breakdowns[fault_tree_events_LV1[i]].keys())
+        fault_tree_events_LV2 = list(recovery_day[fault_tree_events_LV1[i]].keys())
         for j in range(len(fault_tree_events_LV2)):
             day_tenant_unit_reoccupiable = np.fmax(day_tenant_unit_reoccupiable, 
                                                recovery_day[fault_tree_events_LV1[i]][fault_tree_events_LV2[j]].reshape(num_reals, int(np.size(recovery_day[fault_tree_events_LV1[i]][fault_tree_events_LV2[j]])/num_reals)))

diff --git a/inputs/Inputs2Copy/optional_inputs.py b/inputs/Inputs2Copy/optional_inputs.py
@@ -71,6 +71,10 @@
 "calculate_red_tag" : 1,
 "red_tag_clear_time" : 7,
 "red_tag_clear_beta" : 0.6,
+"red_tag_options" : {
+    "tag_coupling_beams_over_height" : True,
+    "ignore_coupling_beam_for_red_tag" : False
+    },
 "include_local_stability_impact" : 1,
 "include_flooding_impact": 1,
 "egress_threshold" : 0.5,
@@ -80,14 +84,15 @@
 "exterior_safety_threshold" : 0.1,
 "interior_safety_threshold" : 0.25,
 "door_access_width_ft" : 9,
-"habitability_requirements": {"electrical" : 0,
-                             "water_potable" : 0,
-                             "water_sanitary" : 0,
-                             "hvac_ventilation" : 0,
-                             "hvac_heating" : 0,
-                             "hvac_cooling" : 0,
-                             "hvac_exhaust" : 0                       
-                             },
+"habitability_requirements": {
+    "electrical" : 0,
+    "water_potable" : 0,
+    "water_sanitary" : 0,
+    "hvac_ventilation" : 0,
+    "hvac_heating" : 0,
+    "hvac_cooling" : 0,
+    "hvac_exhaust" : 0                       
+    },
 "water_pressure_max_story" : 4,
 "heat_utility" : 'gas',
                         }

diff --git a/main_PBEE_recovery.py b/main_PBEE_recovery.py
@@ -71,7 +71,8 @@ def main_PBEE_recovery(damage, damage_consequences, building_model,
     ## Calculate Red Tags
     RT, RTI, IT = fn_red_tag(functionality_options['calculate_red_tag'], 
                                     damage, building_model['comps'],
-                                    np.array(damage_consequences['simulated_replacement_time']))
+                                    np.array(damage_consequences['simulated_replacement_time']),
+                                    functionality_options['red_tag_options'])
 
     damage_consequences['red_tag'] = RT 
     damage_consequences['red_tag_impact'] = RTI 

diff --git a/requirements.txt b/requirements.txt
@@ -0,0 +1,5 @@
+numpy
+pandas
+scipy
+matplotlib
+seaborn