From 32d8b4ccef08b861643adb862e599651475153b0 Mon Sep 17 00:00:00 2001 From: dslosky-usgs Date: Fri, 15 Mar 2019 11:23:01 -0600 Subject: [PATCH] Handling intersections at 0,0 better --- setup.py | 3 ++- shakecastaebm/damage.py | 8 +++++++- shakecastaebm/damping.py | 10 +++++++++- shakecastaebm/performance_point.py | 24 ++++++++++++++++++++---- 4 files changed, 38 insertions(+), 7 deletions(-) diff --git a/setup.py b/setup.py index 5591c60..c52d964 100644 --- a/setup.py +++ b/setup.py @@ -29,7 +29,8 @@ # Versions should comply with PEP440. For a discussion on single-sourcing # the version across setup.py and the project code, see # https://packaging.python.org/en/latest/single_source_version.html - version='0.0a10', + + version='0.0b1', description='Potential impact calculation for well defined facilities due to an input earthquake hazard', long_description=long_description, diff --git a/shakecastaebm/damage.py b/shakecastaebm/damage.py index 52e3327..9dae94c 100644 --- a/shakecastaebm/damage.py +++ b/shakecastaebm/damage.py @@ -116,7 +116,13 @@ def get_damage_state_beta(default_beta, default_median, lower_bound_demand_disp, beta_c = uncertainty_lookup[quality_rating][performance_rating]['beta_c'] beta_t = uncertainty_lookup[quality_rating][performance_rating]['beta_t'] - return min(1.1 * default_beta, max(.9 * default_beta, math.sqrt(min(max(math.log(min(upper_bound_demand_disp, 1.2 * default_median) / min(lower_bound_demand_disp, 1.2 * default_median))/2, demand_uncertainty / 2), demand_uncertainty * 2)**2) + min(max(math.log(min(upper_bound_demand_acc, 1.2 * default_median) / min(lower_bound_demand_acc, 1.2 * default_median))/2, beta_c / 2), 2 * beta_c)**2 + beta_t**2)) + lower_bound_demand_acc = .0000000001 if lower_bound_demand_acc == 0 else lower_bound_demand_acc + upper_bound_demand_acc = .0000000001 if upper_bound_demand_acc == 0 else upper_bound_demand_acc + lower_bound_demand_disp = .0000000001 if lower_bound_demand_disp == 0 else lower_bound_demand_disp + upper_bound_demand_disp = .0000000001 if upper_bound_demand_disp == 0 else upper_bound_demand_disp + return min(1.1 * default_beta, + max(.9 * default_beta, + math.sqrt(min(max(math.log(min(upper_bound_demand_disp, 1.2 * default_median) / min(lower_bound_demand_disp, 1.2 * default_median))/2, demand_uncertainty / 2), demand_uncertainty * 2)**2) + min(max(math.log(min(upper_bound_demand_acc, 1.2 * default_median) / min(lower_bound_demand_acc, 1.2 * default_median))/2, beta_c / 2), 2 * beta_c)**2 + beta_t**2)) def get_damage_probabilities( damage_state_medians, diff --git a/shakecastaebm/damping.py b/shakecastaebm/damping.py index 7057707..016977d 100644 --- a/shakecastaebm/damping.py +++ b/shakecastaebm/damping.py @@ -98,7 +98,15 @@ def get_kappa(spr, year, mag, r_rup): lower_mag = 7 upper_mag = 7.25 - r_rup = int(round(r_rup)) if r_rup < 50 else 50 + if r_rup <= 15: + r_rup = int(round(r_rup)) + elif r_rup < 50: + r_rup = int(round(r_rup / 5) * 5) + elif r_rup < 1000: + r_rup = 50 + else: + r_rup = 1000 + spr = 'non-baseline' if spr != 'baseline' else spr if mag <= 6.25 or mag >= 7.25: diff --git a/shakecastaebm/performance_point.py b/shakecastaebm/performance_point.py index 6a94412..30e9ffe 100644 --- a/shakecastaebm/performance_point.py +++ b/shakecastaebm/performance_point.py @@ -2,6 +2,14 @@ from .spectrum import linear_interpolate, build_spectrum def average_intersections(intersections, capacity, demand): + if len(intersections) < 3: + if intersections[0]['acc'] == 0: + return intersections[1] + elif intersections[1]['acc'] == 0: + return intersections[0] + else: + return intersections[1] + first_point = intersections[0] second_point = intersections[1] third_point = intersections[2] @@ -111,14 +119,22 @@ def get_performance_point(capacity, demand): # calculate periods for intersections for intersection in intersections: - period = math.sqrt(intersection['disp'] / intersection['acc'] / 9.779738) - intersection['period'] = round(period*100) / 100 + if intersection['acc'] != 0: + period = math.sqrt(intersection['disp'] / intersection['acc'] / 9.779738) + intersection['period'] = round(period*100) / 100 + else: + intersection['period'] = 0 if len(intersections) == 1: performance_point = intersections[0] - else: + elif len(intersections) > 1: performance_point = average_intersections(intersections, capacity, demand) - + else: + performance_point = { + 'disp': 0, + 'acc': 0 + } + return performance_point # determine performance point from multiple intersections