diff --git a/setup.py b/setup.py index d52dd19..2c65973 100644 --- a/setup.py +++ b/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, diff --git a/shakecastaebm/capacity.py b/shakecastaebm/capacity.py index e6430b6..e699bf8 100644 --- a/shakecastaebm/capacity.py +++ b/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 diff --git a/shakecastaebm/core.py b/shakecastaebm/core.py index 274c6e9..8544909 100644 --- a/shakecastaebm/core.py +++ b/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 diff --git a/shakecastaebm/performance_point.py b/shakecastaebm/performance_point.py index a7f563c..75e31d5 100644 --- a/shakecastaebm/performance_point.py +++ b/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 diff --git a/shakecastaebm/spectrum.py b/shakecastaebm/spectrum.py index 40975be..46b7816 100644 --- a/shakecastaebm/spectrum.py +++ b/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))) diff --git a/shakecastaebm/tests/performance_point.py b/shakecastaebm/tests/performance_point.py index aa02a43..cb1c12f 100644 --- a/shakecastaebm/tests/performance_point.py +++ b/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()