Merge pull request #30 from dslosky-usgs/updates

Handling intersections at 0,0 better

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,

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,

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:

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