Merge pull request #15 from dslosky-usgs/updates

Make line descriptions better and add some test coverage

dev-requirements.txt

@@ -14,6 +14,7 @@ matplotlib==2.2.2
 numpy==1.15.0
 packaging==17.1
 Pygments==2.2.0
+pylint==1.9.3
 pyparsing==2.2.0
 python-dateutil==2.7.3
 pytz==2018.5

shakecastaebm/capacity.py

@@ -422,13 +422,13 @@ def get_yield_point(design_coefficient, elastic_period, modal_weight, pre_yield)
   a_y = design_coefficient * pre_yield / modal_weight
   d_y = 386.08858 / (4 * math.pi**2) * a_y * elastic_period**2
 
-  return {'x': d_y, 'y': a_y}
+  return {'disp': d_y, 'acc': a_y}
 
 def get_ultimate_point(ductility, d_y, a_y, post_yield):
   a_u = a_y * post_yield
   d_u = ductility * d_y * post_yield
 
-  return {'x': d_u, 'y': a_u}
+  return {'disp': d_u, 'acc': a_u}
 
 def get_ultimate_period(d, a):
     return math.sqrt(d / (a * 9.779738))
@@ -461,48 +461,48 @@ def get_capacity_curve(d_y, a_y, d_u, a_u, elastic_points=5, elipse_points=15, u
     a = math.sqrt(d_y / a_y * b**2 * (d_u - d_y) / (a_y - k))
 
     # anchor at (0,0)
-    points = [{'x': 0, 'y': 0}]
+    points = [{'disp': 0, 'acc': 0}]
 
     # pick points for elestic section
     slope = a_y / d_y
     d = 0
     incr = d_y / elastic_points
     while d < d_y:
-        points += [{'x': d, 'y': d * slope}]
+        points += [{'disp': d, 'acc': d * slope}]
 
         d += incr
 
     # add elastic period
-    points += [{'x': d_y, 'y': a_y}]
+    points += [{'disp': d_y, 'acc': a_y}]
 
     # pick points between dy and du
     incr = (d_u - d_y) / elipse_points
     d = d_y + incr
     while d < d_u:
         point = {
-            'x': d,
-            'y': b * math.sqrt((1 - (d - d_u)**2 / a**2)) + k
+            'disp': d,
+            'acc': b * math.sqrt((1 - (d - d_u)**2 / a**2)) + k
         }
 
         points += [point]
         d += incr
 
     # add ultimate point
-    points += [{'x': d_u, 'y': a_u}]
+    points += [{'disp': d_u, 'acc': a_u}]
 
     # pick points between du and the end of the curve
     incr = (d_u * 10 - d_u) / (ultimate_points / 2)
     d = d_u + incr
     while d < d_u * 10:
-        points += [{'x': d, 'y': a_u}]
+        points += [{'disp': d, 'acc': a_u}]
 
         d += incr
 
     # extend the curve further, but less precisely
     incr = (d_u * 100 - d) / (ultimate_points / 2)
     d = d_u + incr
     while d < d_u * 100:
-        points += [{'x': d, 'y': a_u}]
+        points += [{'disp': d, 'acc': a_u}]
 
         d += incr
 
@@ -576,10 +576,10 @@ def get_capacity(mbt, sdl, bid, height, stories, year, performance_rating='basel
     )
 
     yield_point = get_yield_point(design_coefficient, elastic_period, modal_weight, pre_yield)
-    ultimate_point = get_ultimate_point(ductility, yield_point['x'], yield_point['y'], post_yield)
-    ultimate_period = ultimate_period if ultimate_period else get_ultimate_period(ultimate_point['x'], ultimate_point['y'])
+    ultimate_point = get_ultimate_point(ductility, yield_point['disp'], yield_point['acc'], post_yield)
+    ultimate_period = ultimate_period if ultimate_period else get_ultimate_period(ultimate_point['disp'], ultimate_point['acc'])
     damage_state_medians = get_damage_state_medians(mbt, sdl, performance_rating, height, modal_height, modal_response)
-    capacity_curve = get_capacity_curve(yield_point['x'], yield_point['y'], ultimate_point['x'], ultimate_point['y'])
+    capacity_curve = get_capacity_curve(yield_point['disp'], yield_point['acc'], ultimate_point['disp'], ultimate_point['acc'])
 
     return {
           'curve': capacity_curve,

shakecastaebm/core.py

@@ -1,7 +1,7 @@
 from . import damping
 from .demand import get_demand
 from .damage import get_damage_state_beta, get_damage_probabilities
-from .performance_point import performance_point
+from .performance_point import get_performance_point
 from .data_tables import pref_periods
 
 def run(capacity, hazard, hazard_beta,
@@ -50,12 +50,28 @@ def run(capacity, hazard, hazard_beta,
 
     '''
     demand, lower_demand, upper_demand = get_demand(hazard, hazard_beta, pref_periods, capacity, mag, r_rup)
-    med_intersections = performance_point(capacity['curve'], demand)
-    lower_intersections = performance_point(capacity['curve'], lower_demand)
-    upper_intersections = performance_point(capacity['curve'], upper_demand)
+    med_intersections = get_performance_point(capacity['curve'], demand)
+    lower_intersections = get_performance_point(capacity['curve'], lower_demand)
+    upper_intersections = get_performance_point(capacity['curve'], upper_demand)
 
-    capacity['calcucated_beta'] = get_damage_state_beta(capacity['default_damage_state_beta'], capacity['damage_state_medians']['complete'], lower_intersections[0]['x'], lower_intersections[0]['y'], upper_intersections[0]['x'], upper_intersections[0]['y'], hazard_beta, capacity['quality_rating'], capacity['performance_rating'], capacity['year'], capacity['stories'])
+    capacity['calcucated_beta'] = get_damage_state_beta(
+        capacity['default_damage_state_beta'],
+        capacity['damage_state_medians']['complete'],
+        lower_intersections[0]['disp'],
+        lower_intersections[0]['acc'],
+        upper_intersections[0]['disp'],
+        upper_intersections[0]['acc'],
+        hazard_beta,
+        capacity['quality_rating'],
+        capacity['performance_rating'],
+        capacity['year'],
+        capacity['stories']
+    )
 
-    damage_probabilities = get_damage_probabilities(capacity['damage_state_medians'], capacity['calcucated_beta'], med_intersections[0]['x'])
+    damage_probabilities = get_damage_probabilities(
+        capacity['damage_state_medians'],
+        capacity['calcucated_beta'],
+        med_intersections[0]['disp']
+    )
 
     return damage_probabilities, capacity, demand, lower_demand, upper_demand, med_intersections, lower_intersections, upper_intersections

shakecastaebm/damping.py

@@ -4,9 +4,9 @@
 from .data_tables import degradation_factor
 
 def get_b_eff(capacity, kappa):
-    last_b_h = 0;
-    b_eff = [0] * len(capacity['curve'])
-    b_eff[0] = [{'x': 0, 'y': capacity['elastic_damping'] * 100}]
+    last_b_h = 0
+    b_eff_lst = [0] * len(capacity['curve'])
+    b_eff_lst[0] = [{'x': 0, 'y': capacity['elastic_damping'] * 100}]
 
     for i in range(len(capacity['curve'])):
         point = capacity['curve'][i]
@@ -16,52 +16,52 @@ def get_b_eff(capacity, kappa):
         else:
             prev_point = None
 
-        t = 0;
-        if point['y'] > 0:
-            t = math.sqrt(point['x'] /
-                    (9.779738 * point['y']));
+        t = 0
+        if point['acc'] > 0:
+            t = math.sqrt(point['disp'] /
+                    (9.779738 * point['acc']))
 
         if t >= capacity['elastic_period'] and prev_point is not None:
-            b_h = 100 * (kappa * ( 2*(point['y'] + prev_point['y']) *
-                        (point['x']-(prev_point['x'] + (capacity['yield_point']['x']/capacity['yield_point']['y']) *
-                        (point['y']-prev_point['y'])))+(((last_b_h/100)/kappa)) *
-                        2 * math.pi * prev_point['x']*prev_point['y'])/(2*math.pi*point['x']*point['y']))
+            b_h = 100 * (kappa * ( 2*(point['acc'] + prev_point['acc']) *
+                        (point['disp']-(prev_point['disp'] + (capacity['yield_point']['disp']/capacity['yield_point']['acc']) *
+                        (point['acc']-prev_point['acc'])))+(((last_b_h/100)/kappa)) *
+                        2 * math.pi * prev_point['disp']*prev_point['acc'])/(2*math.pi*point['disp']*point['acc']))
 
-            last_b_h = b_h;
+            last_b_h = b_h
 
-            b_eff[i] = {'x': t, 'y': max(b_h, capacity['elastic_damping'] * 100)};
+            b_eff_lst[i] = {'x': t, 'y': max(b_h, capacity['elastic_damping'] * 100)}
         else:
-            b_eff[i] = {'x': t, 'y': (capacity['elastic_damping'] * 100)};
+            b_eff_lst[i] = {'x': t, 'y': (capacity['elastic_damping'] * 100)}
 
-    return b_eff
+    return b_eff_lst
 
 def get_dsf(beta, mag, rRup):
     dsf = []
     for i in range(len(beta)):
         lnDSF = ((sanaz.b0[i]['y'] + sanaz.b1[i]['y']*math.log(beta[i]['y']) + sanaz.b2[i]['y']*((math.log(beta[i]['y']))**2)) +
                         (sanaz.b3[i]['y'] + sanaz.b4[i]['y']*math.log(beta[i]['y']) + sanaz.b5[i]['y']*((math.log(beta[i]['y']))**2)) * mag +
-                        (sanaz.b6[i]['y'] + sanaz.b7[i]['y']*math.log(beta[i]['y']) + sanaz.b8[i]['y']*((math.log(beta[i]['y']))**2)) * math.log(rRup+1));
+                        (sanaz.b6[i]['y'] + sanaz.b7[i]['y']*math.log(beta[i]['y']) + sanaz.b8[i]['y']*((math.log(beta[i]['y']))**2)) * math.log(rRup+1))
 
-        dsf += [{'x': beta[i]['x'], 'y': round(math.exp(lnDSF), 3)}];
+        dsf += [{'x': beta[i]['x'], 'y': round(math.exp(lnDSF), 3)}]
 
     return dsf
 
 def damp(demand, capacity, mag, r_rup):
     kappa = get_kappa(capacity['performance_rating'], capacity['year'], mag, r_rup)
     b_eff_spectrum = get_b_eff(capacity, kappa)
 
-    beta = build_spectrum(b_eff_spectrum, sanaz.t);
+    beta = build_spectrum(b_eff_spectrum, sanaz.t)
     dsf = get_dsf(beta, mag, r_rup)
 
     # expand dsf to match demand spectrum periods
-    dsf = build_spectrum(dsf, [point['x'] for point in demand])
+    dsf = build_spectrum(dsf, [point['period'] for point in demand])
 
     damped_demand = [0] * len(demand)
     for i in range(len(demand)):
-        damp_disp = demand[i]['disp'] * dsf[i]['y'];
-        damp_acc = damp_disp/(9.779738 * demand[i]['x']**2);
+        damp_disp = demand[i]['disp'] * dsf[i]['y']
+        damp_acc = damp_disp/(9.779738 * demand[i]['period']**2)
 
-        damped_demand[i] = {'disp': damp_disp, 'y': damp_acc, 'x': demand[i]['x']};
+        damped_demand[i] = {'disp': damp_disp, 'acc': damp_acc, 'period': demand[i]['period']}
 
     return damped_demand
 

shakecastaebm/demand.py

@@ -3,14 +3,14 @@
 
 import math
 
-def make_demand_spectrum(input):
+def make_demand_spectrum(input, x='x', y='y'):
     demand = []
     for point in input:
-        disp = point['y'] * point['x']**2 * 9.779738
-        acc = disp/(9.779738 * (point['x']**2))
+        disp = point[y] * point[x]**2 * 9.779738
+        acc = disp/(9.779738 * (point[x]**2))
         demand += [{
-            'x': point['x'],
-            'y': acc,
+            'period': point[x],
+            'acc': acc,
             'disp': disp
         }]
 
@@ -27,13 +27,13 @@ def get_demand(hazard, hazard_beta, pref_periods, capacity, mag, rRup):
         lower_bound_demand += [
             {
                 'disp': point['disp'] / math.exp(hazard_beta),
-                'y': point['y'] / math.exp(hazard_beta)
+                'acc': point['acc'] / math.exp(hazard_beta)
             }
         ]
         upper_bound_demand += [
             {
                 'disp': point['disp'] * math.exp(hazard_beta),
-                'y': point['y'] * math.exp(hazard_beta)
+                'acc': point['acc'] * math.exp(hazard_beta)
             }
         ]
 

shakecastaebm/performance_point.py

@@ -1,58 +1,72 @@
-def performance_point(capacity, demand):
-    dem = [{'x': point['disp'], 'y': point['y']} for point in demand]
-    intersections = find_intersections(capacity, dem)
+import math
+
+def get_performance_point(capacity, demand):
+    intersections = find_intersections(capacity, demand, 'disp', 'acc')
+
+    # calculate periods for intersections
+    for intersection in intersections:
+        period = math.sqrt(intersection['disp'] / intersection['acc'] /  9.779738)
+        intersection['period'] = round(period*100) / 100
 
     if len(intersections) == 1:
         return intersections
 
     # determine performance point from multiple intersections
 
-def find_intersections(line1, line2):
-    intersections = [];
-    line1_idx = 0;
-    line2_idx = 0;
+def find_intersections(line1, line2, x='x', y='y'):
+    intersections = []
+    line1_idx = 0
+    line2_idx = 0
     while line1_idx < len(line1) - 1 and line2_idx < len(line2) - 1:
         seg1 = [line1[line1_idx], line1[line1_idx + 1]]
         seg2 = [line2[line2_idx], line2[line2_idx + 1]]
 
-        intersection = get_intersection(seg1, seg2)
+        intersection = get_intersection(seg1, seg2, x, y)
         if intersection is not False:
             intersections += [intersection]
 
-        if seg1[1]['x'] == seg2[1]['x']:
+        if seg1[1][x] == seg2[1][x]:
             line1_idx += 1
             line2_idx += 1
-        elif seg1[1]['x'] < seg2[1]['x']:
+        elif seg1[1][x] < seg2[1][x]:
             line1_idx += 1
         else:
             line2_idx += 1
 
     return intersections
 
-def get_intersection(seg1, seg2):
+def get_intersection(seg1, seg2, x, y):
     # seg1 and seg2 are 2d arrays 
     # [ {x: x1_val,  y: y1_val}, {x: x2_val, y: y2_val} ]
-    dx1 = seg1[1]['x'] - seg1[0]['x'];
-    dx2 = seg2[1]['x'] - seg2[0]['x'];
-    dy1 = seg1[1]['y'] - seg1[0]['y'];
-    dy2 = seg2[1]['y'] - seg2[0]['y'];
+    dx1 = seg1[1][x] - seg1[0][x]
+    dx2 = seg2[1][x] - seg2[0][x]
+    dy1 = seg1[1][y] - seg1[0][y]
+    dy2 = seg2[1][y] - seg2[0][y]
 
     # check for parallel segs
     if (dx2 * dy1 - dy2 * dx1) == 0:
         # The segments are parallel.
         return False
 
-    s = ((dx1 * (seg2[0]['y'] - seg1[0]['y']) + dy1 * (seg1[0]['x'] - seg2[0]['x'])) /
-                (dx2 * dy1 - dy2 * dx1));
-    t = ((dx2 * (seg1[0]['y'] - seg2[0]['y']) + dy2 * (seg2[0]['x'] - seg1[0]['x'])) /
-                (dy2 * dx1 - dx2 * dy1));
+    s = ((dx1 * (seg2[0][y] - seg1[0][y]) + dy1 * (seg1[0][x] - seg2[0][x])) /
+                (dx2 * dy1 - dy2 * dx1))
+    t = ((dx2 * (seg1[0][y] - seg2[0][y]) + dy2 * (seg2[0][x] - seg1[0][x])) /
+                (dy2 * dx1 - dx2 * dy1))
 
     if (s >= 0 and s <= 1 and t >= 0 and t <= 1):
-        return {'x': seg1[0]['x'] + t * dx1, 'y': seg1[0]['y'] + t * dy1};
+        x_val = abs(round(seg1[0][x] + t * dx1 * 1000) / 1000)
+        y_val = abs(round(seg1[0][y] + t * dy1 * 1000) / 1000)
+        return {x: x_val, y: y_val}
     else:
       return False
 
 def weight_intersections(intersections, line1, line2):
+    '''
+    TODO: WRITE THIS FUNCTION
+
+    In the event of multiple intersections, the performance point
+    should be a weighted average of the outer intersections
+    '''
     #on_top = []
     #inter_idx = 0
     #for idx in range(len(line1)):