Merge pull request #22 from dslosky-usgs/updates

Triple intersection

setup.py

@@ -29,7 +29,7 @@
     # 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.0a6',
+    version='0.0a7',
 
     description='Potential impact calculation for well defined facilities due to an input earthquake hazard',
     long_description=long_description,

shakecastaebm/capacity.py

@@ -523,13 +523,13 @@ def get_capacity(mbt, sdl, bid, height, stories, year, performance_rating='basel
 
     Args:
         mbt: HAZUS model building type
-        sdl: seasmic design level (special, high, moderate, low)
+        sdl: seasmic design level (special_high, high, moderate, low)
         bid: basis ID, 1-7, describes structural deficiencies from FEMA 155
         height: height of the structure in feet
         stories: count of stories above ground.
         year: year the building was constructured or retrofit
         performance_rating: DEFAULT 'baseline', structural performance rating (baseline, poor, or very_poor)
-        quality_rating: DEFAULT 'poor', The quality of this structural data (high, moderate, poor, very_poor) elastic_period=None, 
+        quality_rating: DEFAULT 'poor', The quality of this structural data (best, very_good, good, poor, very_poor) elastic_period=None, 
         elastic_damping=None, design_period=None, ultimate_period=None, design_coefficient=None, modal_weight=None,
         modal_height=None, modal_response=None, pre_yield=None, post_yield=None,
         max_strength=None, ductility=None, default_damage_state_beta=None

shakecastaebm/core.py

@@ -57,10 +57,10 @@ def run(capacity, hazard, hazard_beta,
     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'],
+        lower_intersections['disp'],
+        lower_intersections['acc'],
+        upper_intersections['disp'],
+        upper_intersections['acc'],
         hazard_beta,
         capacity['quality_rating'],
         capacity['performance_rating'],
@@ -71,7 +71,7 @@ def run(capacity, hazard, hazard_beta,
     damage_probabilities = get_damage_probabilities(
         capacity['damage_state_medians'],
         capacity['calcucated_beta'],
-        med_intersections[0]['disp']
+        med_intersections['disp']
     )
 
     return damage_probabilities, capacity, demand, lower_demand, upper_demand, med_intersections, lower_intersections, upper_intersections

shakecastaebm/performance_point.py

@@ -1,4 +1,110 @@
 import math
+from .spectrum import linear_interpolate, build_spectrum
+
+def average_intersections(intersections, capacity, demand):
+    first_point = intersections[0]
+    second_point = intersections[1]
+    third_point = intersections[2]
+
+    capacity_seg1 = chop_curve(capacity,
+            first_point['disp'], second_point['disp'], 'disp', 'acc')
+    capacity_seg2 = chop_curve(capacity,
+            second_point['disp'], third_point['disp'], 'disp', 'acc')
+
+    demand_seg1 = chop_curve(demand,
+            first_point['disp'], second_point['disp'], 'disp', 'acc')
+    demand_seg2 = chop_curve(demand,
+            second_point['disp'], third_point['disp'], 'disp', 'acc')
+
+    capacity_area1 = get_area_under_curve(capacity_seg1, 'disp', 'acc')
+    capacity_area2 = get_area_under_curve(capacity_seg2, 'disp', 'acc')
+    demand_area1 = get_area_under_curve(demand_seg1, 'disp', 'acc')
+    demand_area2 = get_area_under_curve(demand_seg2, 'disp', 'acc')
+
+    area_diff1 = abs(capacity_area1 - demand_area1)
+    area_diff2 = abs(capacity_area2 - demand_area2)
+
+    x_difference = third_point['disp'] - first_point['disp']
+
+    if area_diff1 > area_diff2:
+        adjust_percentage = area_diff2 / (area_diff1 + area_diff2)
+        performance_point_disp = first_point['disp'] + (x_difference * adjust_percentage)
+    else:
+        adjust_percentage = area_diff1 / (area_diff1 + area_diff2)
+        performance_point_disp = third_point['disp'] - (x_difference * adjust_percentage)
+
+    performance_point = find_point_on_curve(performance_point_disp,
+            capacity, 'disp', 'acc')
+
+    return performance_point
+
+def chop_curve(curve, start_val, stop_val, x='x', y='y'):
+    idx = 0
+    chopped_curve = []
+    for point in curve:
+        if point[x] == start_val:
+            chopped_curve += [point]
+        elif point[x] > start_val and point[x] < stop_val:
+            if idx > 0 and len(chopped_curve) == 0:
+                # interpolate to get the right value
+                chopped_curve += [
+                    {
+                        'acc': linear_interpolate(start_val, curve[idx - 1], curve[idx], x, y),
+                        'disp': start_val
+                    },
+                    point
+                ]
+            else:
+                chopped_curve += [point]
+
+        elif point[x] > stop_val:
+            # interpolate final point
+            chopped_curve += [
+                {
+                    'acc': linear_interpolate(stop_val, curve[idx - 1], curve[idx], x, y),
+                    'disp': stop_val
+                }
+            ]
+            return chopped_curve
+
+        idx += 1
+    return chopped_curve
+
+def find_point_on_curve(x_val, curve, x = 'x', y = 'y'):
+    point = {
+        x: x_val,
+        y: None
+    }
+
+    for idx in range(len(curve) - 1):
+        if curve[idx][x] > x_val:
+            if idx > 0:
+                point[y] = linear_interpolate(x_val, curve[idx-1], curve[idx], x, y)
+            else:
+                point[y] = linear_interpolate(x_val, curve[idx], curve[idx + 1], x, y)
+
+            return point
+
+    return None
+
+def get_area_under_curve(curve, x='x', y='y'):
+    '''
+    Integrates a curve of discrete values
+    '''
+    x_vals = [curve[0][x] + x_int * .001 for x_int in 
+            range(0, int((curve[-1][x] - curve[0][x]) * 1000))]
+    expanded_curve = build_spectrum(curve, x_vals, x='disp', y='acc')
+
+    total_area = 0
+    for idx in range(len(expanded_curve) - 2):
+        p1 = expanded_curve[idx]
+        p2 = expanded_curve[idx + 1]
+
+        area = (p2[x] - p1[x]) * p2[y]
+        total_area += area
+
+    return total_area
+
 
 def get_performance_point(capacity, demand):
     intersections = find_intersections(capacity, demand, 'disp', 'acc')
@@ -9,7 +115,11 @@ def get_performance_point(capacity, demand):
         intersection['period'] = round(period*100) / 100
 
     if len(intersections) == 1:
-        return intersections
+        performance_point = intersections[0]
+    else:
+        performance_point = average_intersections(intersections, capacity, demand)
+
+    return performance_point
 
     # determine performance point from multiple intersections
 

shakecastaebm/spectrum.py

@@ -82,6 +82,21 @@ def interpolate(interpX, p1, p2, x='x', y='y'):
 
     return val
 
+def linear_interpolate(interpX, p1, p2, x='x', y='y'):
+    # swap input points if they're in the wrong order
+    if p1[x] > p2[x]:
+      p1, p2 = p2, p1
+
+    x1 = p1[x] if p1[x] != 0 else .000000000001
+    y1 = p1[y] if p1[y] != 0 else .000000000001
+    x2 = p2[x] if p2[x] != 0 else .000000000001
+    y2 = p2[y] if p2[y] != 0 else .000000000001
+
+    val = (y1 * (1 - (interpX - x1) / (x2 - x1)) + 
+            y2 * ((interpX - x1) / (x2 - x1)))
+
+    return val
+
 def insert_vals(lst, vals):
     lst = sorted(list(set(lst)))
     vals = sorted(list(set(vals)))

shakecastaebm/tests/performance_point.py

@@ -14,14 +14,12 @@ def test_performancePoint(self):
         ]
 
         pp = get_performance_point(capacity, demand)
-        self.assertTrue(isinstance(pp, list))
-
-        pp = pp[0]
+        self.assertTrue(isinstance(pp, dict))
         self.assertTrue('disp' in pp.keys())
         self.assertTrue('acc' in pp.keys())
         self.assertTrue('period' in pp.keys())
 
-    def test_intersection(self):
+    def test_Intersection(self):
         capacity = [
             {'disp': 1, 'acc': .1},
             {'disp': 2, 'acc': 1}
@@ -31,7 +29,7 @@ def test_intersection(self):
             {'disp': 2, 'acc': .1}
         ]
 
-        pp = get_performance_point(capacity, demand)[0]
+        pp = get_performance_point(capacity, demand)
 
         # slopes for the intersection with the demand
         dem_slope1 = ((demand[0]['acc'] - pp['acc']) /
@@ -48,6 +46,74 @@ def test_intersection(self):
         self.assertAlmostEqual(dem_slope1, dem_slope2)
         self.assertAlmostEqual(cap_slope1, cap_slope2)
 
+    def test_TripleIntersection(self):
+        capacity = [
+            {'disp': .1, 'acc': 1},
+            {'disp': .5, 'acc': 1},
+            {'disp': 1, 'acc': 1},
+            {'disp': 2, 'acc': 1},
+            {'disp': 3, 'acc': 1},
+            {'disp': 4, 'acc': 1}
+        ]
+
+        demand = [
+            {'disp': .1, 'acc': 2},
+            {'disp': .5, 'acc': .5},
+            {'disp': 2, 'acc': 1.5},
+            {'disp': 3, 'acc': .1}
+        ]
+
+        perf_point = get_performance_point(capacity, demand)
+        intersections = find_intersections(capacity, demand, 'disp', 'acc')
+
+        self.assertTrue(perf_point is not None)
+        self.assertTrue(perf_point['disp'] < intersections[-1]['disp'])
+
+    def test_ChopCurve(self):
+        curve = [
+            {'disp': .1, 'acc': 1},
+            {'disp': .5, 'acc': 1},
+            {'disp': 1, 'acc': 1},
+            {'disp': 2, 'acc': 1},
+            {'disp': 3, 'acc': 1},
+            {'disp': 4, 'acc': 1}
+        ]
+
+        start = .3
+        end = 2.2
+
+        chopped = chop_curve(curve, start, end, 'disp', 'acc')
+
+        self.assertAlmostEqual(start, chopped[0]['disp'])
+        self.assertAlmostEqual(end, chopped[-1]['disp'])
+
+    def test_FindPointOnCurve(self):
+        curve = [
+            {'disp': .1, 'acc': 1},
+            {'disp': .5, 'acc': 1},
+            {'disp': 1, 'acc': 1},
+            {'disp': 2, 'acc': 1},
+            {'disp': 3, 'acc': 1},
+            {'disp': 4, 'acc': 1}
+        ]
+
+        point_disp = 1.7
+        point = find_point_on_curve(point_disp, curve, 'disp', 'acc')
+
+        self.assertAlmostEqual(point_disp, point['disp'])
+        self.assertAlmostEqual(1, point['acc'])
+
+    def test_GetAreaUnderCurve(self):
+        curve = [
+            {'disp': 1, 'acc': 1},
+            {'disp': 2, 'acc': 1},
+            {'disp': 3, 'acc': 1},
+            {'disp': 4, 'acc': 1}
+        ]
+
+        area = get_area_under_curve(curve, 'disp', 'acc')
+
+        self.assertAlmostEqual(area, 3, 2)
 
 if __name__ == '__main__':
     unittest.main()