Extended testing fixes, simplified np.array reading, citation file

CITATION.cff

@@ -0,0 +1,49 @@
+cff-version: 1.0.0
+message: "If you use this software, please cite it as below."
+authors:
+- family-names: "Zaidi"
+  given-names: "Faraz"
+- family-names: "Pham"
+  given-names: "Huy G."
+  orcid: "https://orcid.org/0000-0002-8043-7784"
+- family-names: "Wang"
+  given-names: "Ziyi"
+- family-names: "Cook"
+  given-names: "Dustin T."
+  orcid: "https://orcid.org/0000-0002-4841-4515"
+- family-names: "Zsarnoczay"
+  given-names: "Adam"
+  orcid: "https://orcid.org/0000-0001-6895-1417"
+- family-names: "Molina Hutt"
+  given-names: "Carlos"
+  orcid: "https://orcid.org/0000-0003-2116-1201"
+- family-names: "Hulsey"
+  given-names: "Anne M."
+  orcid: "https://orcid.org/0000-0002-4265-055X"
+title: "OpenPBEE/Functional-Recovery-Python"
+version: 1.3.0
+doi: https://doi.org/10.5281/zenodo.20045328
+date-released: 2026-03-31
+url: "https://github.com/OpenPBEE/Functional-Recovery-Python"
+preferred-citation:
+  type: article
+  authors:
+  - family-names: "Cook"
+    given-names: "Dustin T."
+    orcid: "https://orcid.org/0000-0002-4841-4515"
+  - family-names: "Liel"
+    given-names: "Abbie B."
+    orcid: "https://orcid.org/0000-0002-9241-5144"
+  - family-names: "Haselton"
+    given-names: "Curt B."
+  - family-names: "Koliou"
+    given-names: "Maria"
+    orcid: "https://orcid.org/0000-0002-0686-493X"
+  doi: https://doi.org/10.1177/87552930221081538
+  journal: "Earthquake Spectra"
+  start: 1972 # First page number
+  end: 2007 # Last page number
+  title: "A framework for operationalizing the assessment of post-earthquake functional recovery of buildings"
+  issue: 3
+  volume: 38
+  year: 2022

LICENSE

@@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
+SOFTWARE.

README.md

@@ -258,4 +258,4 @@ python tests/compare_runs.py <reference_dir> <python_output_dir>
 
 ### Validation Testing Against MATLAB
 
-For detailed guidance on validating the Python implementation against the original MATLAB framework, see [tests/detailed_inspection_guide.md](tests/detailed_inspection_guide.md).
+For detailed guidance on validating the Python implementation against the original MATLAB framework, see [tests/README.md](tests/README.md).

plotters/other_plot_functions.py

@@ -249,9 +249,9 @@ def plt_recovery_trajectory(recovery, full_repair_time, plot_dir):
     # Plot Recovery Trajectory
     plt.rcParams["font.family"] = "Times New Roman"
     plt.figure(figsize=(6,4)) 
-    plt.plot(reoc, level_of_repair,'r-', linewidth = 1.5, label = 'Re-Occupancy') 
-    plt.plot(func, level_of_repair,'b-', linewidth = 1.5, label= 'Functional') 
-    plt.plot([full, full], [0, 1],'k-', linewidth = 1.5, label= 'Fully Repaired') 
+    plt.plot(reoc, level_of_repair,'r:', linewidth = 1.5, label = 'Re-Occupancy', zorder=2) 
+    plt.plot(func, level_of_repair,'b--', linewidth = 1.5, label= 'Functional',zorder=1) 
+    plt.plot([full, full], [0, 1],'k-', linewidth = 1.5, label= 'Fully Repaired', zorder=0) 
     plt.xlim([0,np.ceil((full+1)/10)*10])
     plt.xlabel('Days After Earthquake')
     plt.ylabel('Fraction of Floor Area')
@@ -414,9 +414,9 @@ def plt_gantt_chart(p_idx, recovery, full_repair_time, workers, schedule, impede
     G3.grid(which = 'major', axis = 'x', alpha=0.5) 
 
     # Plot Recovery Trajectory
-    G4.plot(recovery_trajectory['reoc'], recovery_trajectory['level_of_repair'],'r-', linewidth = 1.5, label = 'Re-Occupancy')
-    G4.plot(recovery_trajectory['func'], recovery_trajectory['level_of_repair'],'b-', linewidth = 1.5, label = 'Functional')
-    G4.plot(recovery_trajectory['ful_rep'], recovery_trajectory['level_of_repair'],'k-', linewidth = 1.5, label = 'Fully Repaired')
+    G4.plot(recovery_trajectory['func'], recovery_trajectory['level_of_repair'],'b--', linewidth = 1.5, label = 'Functional',zorder=1)
+    G4.plot(recovery_trajectory['reoc'], recovery_trajectory['level_of_repair'],'r:', linewidth = 1.5, label = 'Re-Occupancy',zorder=2)
+    G4.plot(recovery_trajectory['ful_rep'], recovery_trajectory['level_of_repair'],'k-', linewidth = 1.5, label = 'Fully Repaired', zorder=0)
     G4.set_ylim([0,1])
     G4.set_xlabel('Days After Earthquake')
     G4.set_ylabel('Fraction of Floor Area')

src/atc138/driver.py

@@ -1,3 +1,4 @@
+import numpy as np
 def run_analysis(input_dir, output_dir, 
                  output_file='recovery_outputs.json',
                  seed=None, force_rebuild=False):
@@ -163,36 +164,11 @@ def run_analysis(input_dir, output_dir,
         os.makedirs(output_dir)
 
     # Covert arrays to list for writing to json file   
-    fnc_keys_1 = list(functionality.keys())
-    for k_1 in fnc_keys_1:
-        if type(functionality[k_1]) == np.ndarray: 
-            functionality[k_1] = functionality[k_1].tolist() 
-        if type(functionality[k_1]) == dict:
-            fnc_keys_2 = list(functionality[k_1].keys())    
-
-            for k_2 in fnc_keys_2:
-                if type(functionality[k_1][k_2]) == np.ndarray: 
-                    functionality[k_1][k_2] = functionality[k_1][k_2].tolist()
-                if type(functionality[k_1][k_2]) == dict:
-                    fnc_keys_3 = list(functionality[k_1][k_2].keys())    
-
-                    for k_3 in fnc_keys_3:
-                        if type(functionality[k_1][k_2][k_3]) == np.ndarray: 
-                            functionality[k_1][k_2][k_3] = functionality[k_1][k_2][k_3].tolist()
-                        if type(functionality[k_1][k_2][k_3]) == dict:
-                            fnc_keys_4 = list(functionality[k_1][k_2][k_3].keys())
-
-                            for k_4 in fnc_keys_4:
-                                if type(functionality[k_1][k_2][k_3][k_4]) == np.ndarray: 
-                                    functionality[k_1][k_2][k_3][k_4] = functionality[k_1][k_2][k_3][k_4].tolist()
-                                if type(functionality[k_1][k_2][k_3][k_4]) == dict:
-                                    fnc_keys_5 = list(functionality[k_1][k_2][k_3][k_4].keys())
-
-                                    for k_5 in fnc_keys_5:
-                                        if type(functionality[k_1][k_2][k_3][k_4][k_5]) == np.ndarray: 
-                                            functionality[k_1][k_2][k_3][k_4][k_5] = functionality[k_1][k_2][k_3][k_4][k_5].tolist()
-                                        if type(functionality[k_1][k_2][k_3][k_4][k_5]) == dict:
-                                            fnc_keys_6 = list(functionality[k_1][k_2][k_3][k_4][k_5].keys())
+    functionality = recursive_convert_array(functionality)
+
+    # join consequences
+    consequences = recursive_convert_array(damage_consequences)
+    functionality = functionality | consequences
 
     output_json_object = json.dumps(functionality)
 
@@ -205,3 +181,32 @@ def run_analysis(input_dir, output_dir,
     print('Recovery assessment complete')
     print('time to run '+str(round(end_time - start_time,2))+'s')
 
+def recursive_convert_array(obj):
+    """
+    Recursively convert numpy arrays to lists,
+    traversing nested dicts/lists/tuples.
+    """
+
+    if isinstance(obj, np.ndarray):
+        return obj.tolist()
+
+    elif isinstance(obj, dict):
+        return {
+            key: recursive_convert_array(value)
+            for key, value in obj.items()
+        }
+
+    elif isinstance(obj, list):
+        return [recursive_convert_array(item) for item in obj]
+
+    elif isinstance(obj, tuple):
+        return tuple(recursive_convert_array(item) for item in obj)
+
+    elif isinstance(obj, np.integer):
+        return int(obj)
+
+    elif isinstance(obj, np.floating):
+        return float(obj)
+
+    else:
+        return obj

src/atc138/functionality/other_functionality_functions.py

@@ -159,7 +159,7 @@ def fn_building_safety(damage, building_model, damage_consequences, utilities,
 
     # Loop through component repair times to determine the day it stops affecting re-occupancy
     num_repair_time_increments = np.sum(damage['fnc_filters']['ext_fall_haz_all'])*num_units # possible unique number of loop increments
-    edge_lengths = np.row_stack((np.array(building_model['edge_lengths']), np.array(building_model['edge_lengths'])))
+    edge_lengths = np.vstack((np.array(building_model['edge_lengths']), np.array(building_model['edge_lengths'])))
 
     affected_ratio={}
     for side in range(4):

src/atc138/input_builder.py

@@ -136,7 +136,8 @@ def build_simulated_inputs(model_dir):
     comp_header = list(comp_population.columns)
 
     # component IDs from comp_population
-    comp_list = np.array(comp_header[2:])
+    comp_meta = ['story', 'dir']
+    comp_list = [cmp_label for cmp_label in comp_header if cmp_label not in comp_meta]
     comp_list = np.char.replace(comp_list, '_', '.')
     comp_list = comp_list.tolist()
 

tests/detailed_inspection_guide.md → tests/README.md

@@ -56,23 +56,31 @@ Threshold: ```≤ 4 %```
 
 **Metric 3 — Worst system mean absolute error (MAE)**
 
-System-level wise, the maximum MAE of the fraction of unresolved damage over time, across all functional-recovery and reoccupancy systems. 
+System-level wise, the highest MAE of the fraction of unresolved damage over time, across all functional-recovery and reoccupancy systems. 
 
-Threshold: ```≤ 0.02``` (maximum mean difference of 2% unresolved fraction)
+Threshold: ```≤ 0.02``` (no more than 2% unresolved fraction across all systems)
 
 </br>
 
 **Metric 4 — Worst system 95th percentile absolute error**
 
-System-level wise, the maximum 95th percentile aboslute error of the fraction of unresolved damage over time, across all functional-recovery and reoccupancy systems. 
+System-level wise, the maximum 95th percentile absolute error of the fraction of unresolved damage over time, across all functional-recovery and reoccupancy systems. 
 
-Threshold: ```≤ 0.04``` (maximum mean difference of 4% unresolved fraction)
-
-Both absolute and relative differences are reported to provide a more complete assessment —— absolute difference reflects the magnitude of discrepancies, while relative difference normalizes them across models with different recovery scales. Using both avoids misleading interpretations —— for example, a small absolute difference may appear large in relative terms for short recovery durations, while a small relative difference may still correspond to a large absolute difference for long recovery durations. 
+Threshold: ```≤ 0.04``` (no more than of 4% unresolved fraction across all systems)
 
 </br>
 
-#### 5. Recovery Trajectory and System-Level Recovery Comparison
+#### 5. Testing thresholds
+
+Both absolute and relative differences are reported to provide a more complete assessment: absolute difference reflects the magnitude of discrepancies, while relative difference normalizes them across models with different recovery scales. Using both avoids misleading interpretations. For example, a small absolute difference may appear large in relative terms for short recovery durations, while a small relative difference may still correspond to a large absolute difference for long recovery durations. 
+
+In the porting process from MATLAB to Python, the output Python results were compared against the reference MATLAB values. The test case is ideally averaged across multiple random seeds, to minimize variance between each realization. The test case is considered to adequately match if:
+
+1) At the high level, either Metric 1 _or_ Metric 2 pass their thresholds. That is, either an absolute difference in mean recovery time $\leq$ 2 days _or_ a relative difference in mean recovery time of $\leq$ 4% is considered passing.
+2) At the system level, both Metric 3 _and_ Metric 4 pass their thresholds. That is, for all systems at each assessed point in time, the average difference in unresolved fraction is no more than 2%, and the 95th-percentile worst difference in unresolved fraction is no more than 4%.
+
+
+#### 6. Recovery Trajectory and System-Level Recovery Comparison
 
 For each model, plots of recovery trajectories and system-level recovery curves are saved: 
 
@@ -81,7 +89,7 @@ For each model, plots of recovery trajectories and system-level recovery curves
 - `system_breakdowns_reoccupancy.jpg`
 
 
-#### 6. Validation Results
+#### 7. Validation Results
 
 The comparison results for the two example models are summarized below. Twenty batch runs were carried out for each model, and the mean values are compared against the defined metrics.
 

tests/compare_runs.py

@@ -359,11 +359,13 @@ def evaluate_tolerances(reference_dir, python_dir):
 
         if ref_stats['mean'] != 0:
             diff = 100.0 * abs(py_stats['mean'] - ref_stats['mean']) / ref_stats['mean']
+            abs_diff = abs(py_stats['mean'] - ref_stats['mean'])
         else:
             diff = float('nan')
+            abs_diff = float('nan')
 
-        # Tolerance: mean recovery days within 3% of reference
-        passed = diff <= 4.0 or (np.isnan(diff) and py_stats['mean'] == 0)
+        # Tolerance: mean recovery days within 4% of reference OR absolute difference within 2 days of reference
+        passed = diff <= 4.0 or abs_diff <= 2.0 or (np.isnan(diff) and py_stats['mean'] == 0)
         if not passed: results["all_passed"] = False
 
         results["high_level"][metric] = {

tests/test_integration.py

@@ -52,7 +52,8 @@ def test_reference_comparison(model_name, model_dir, tmp_path):
 
     reference_dir = model_dir / "reference"
 
-    run_analysis(str(example_dir), str(tmp_path))
+    default_seed = 123
+    run_analysis(str(example_dir), str(tmp_path), seed=default_seed)
     results = evaluate_tolerances(str(reference_dir), str(tmp_path))
 
     hl = results["high_level"]

tests/test_utils.py

@@ -0,0 +1,89 @@
+import atc138.driver
+from pathlib import Path
+import shutil
+import pandas as pd
+
+def copy_defaults(case_dir: Path, default_path: Path):
+    case_dir.mkdir(parents=True, exist_ok=True)
+
+    for f in ['optional_inputs.json',
+              'tenant_function_requirements.csv',
+              'tenant_unit_list.csv']:
+
+        shutil.copy(default_path / f, case_dir / f)
+
+def cleanup(case_dir: Path):
+    # remove intermediate files only
+    for f in case_dir.glob("simulated_inputs.json"):
+        f.unlink()
+
+def run_atc138(case_dir: Path, num_iter=1):
+    '''
+    Functionality to run the extended/canonical test suites. On the first run, rebuild the input file. On the last run,
+    cleanup the large input file.
+
+    case_dir: directory of folder containing the input files
+    num_iter: Number of iterations to be ran. Seeds will be created from 1-num_iter
+
+    '''
+    print(f"Running {case_dir}")
+
+    # path_to_defaults = case_dir.parent.parent / 'atcpy_input'
+    # copy_defaults(case_dir, default_path=path_to_defaults)
+
+    # define outputs explicitly
+    output_dir = case_dir / 'output'
+    output_dir.mkdir(parents=True, exist_ok=True)
+
+    comp_population_csv = case_dir / 'comp_population.csv'
+
+    if not comp_population_csv.is_file():
+        cmp_pop = pd.read_excel(case_dir / 'comp_population.xlsx')
+
+        # handle duplicate components in xlsx by summing
+        # strip all trailing duplicate naming convention of ".1" ".2"
+        # IMPORTANT: this can handle only up to 10 duplicates
+        cmp_pop.columns = cmp_pop.columns.str.replace(r'\.\d$', '', regex=True)
+
+        # group & sum
+        cmp_pop_sum = cmp_pop.T.groupby(level=0).sum().T
+        meta_cols = ['story', 'dir']
+        cmp_pop_sum = cmp_pop_sum.astype({'story': int, 'dir': int})
+        cmp_pop_sum = cmp_pop_sum[meta_cols + [col for col in cmp_pop_sum.columns if col not in meta_cols]]
+
+        cmp_pop_sum.to_csv(comp_population_csv, index=False)
+
+    # run analysis
+    if num_iter > 1:
+        for run_no in range(1, num_iter+1):
+            output_file = f'recovery_outputs_run{run_no}.json'
+            if run_no == 1:
+                atc138.driver.run_analysis(
+                    case_dir, output_dir, seed=run_no, output_file=output_file, force_rebuild=True)
+            else:
+                atc138.driver.run_analysis(
+                    case_dir, output_dir, seed=run_no, output_file=output_file, force_rebuild=False)
+
+    else:
+        atc138.driver.run_analysis(case_dir, output_dir, seed=985, force_rebuild=True)
+
+    # cleanup
+    cleanup(case_dir)
+
+def batch_run_folder(root_dir: Path, n_iter=1):
+    '''
+    go through each folder and run atc138
+
+    supports either a flat directory or a nested directory (i.e. of IM folders)
+    '''
+
+    # case 1: no IMs and building_model.json is directly in dir
+    building_model_file = root_dir / 'building_model.json'
+    if building_model_file.is_file():
+        run_atc138(root_dir, n_iter)
+    # case 2: models are nested within experiment directories
+    else:
+        for im_dir in root_dir.iterdir():
+            if not im_dir.is_dir():
+                continue
+            run_atc138(im_dir, n_iter)