Risk-Hazard Convolution Module

This script takes a fragility function as input and does the convolution with a hazard curve in order to estimate the annual probability of collapse.

Load hazard curve

In order to use this methodology, it is necessary to provide a hazard curve, defined according to the format established on the RMTK manual. Please provide the location of the file containing this input using the parameter capacity_curves_file.


In [30]:
from rmtk.vulnerability.common import utils
%matplotlib inline 

hazard_curve_file = '../../../../../../Google Drive 2/Otros GEM/RMTK/Convolution_module/Hazard_PGA.csv'
hazard_curve = utils.read_hazard(hazard_curve_file)
utils.plot_hazard_curve(hazard_curve)


Load fragility function

The fragility function representing the risk of the building class being analysed should be provided here. This function should follow the format established in the RMTK manual. Please provide the location of the file containing this input using the parameter fragility_function_file.


In [44]:
fragility_function_file = '../../../../../../Google Drive 2/Otros GEM/RMTK/Convolution_module/MUR_H1.csv'
fragility_model = utils.read_frag_model(fragility_function_file)
save = False
utils.plot_fragility_model(fragility_function_file,0.01,2,save)


Calculate annual probability of damage

In order to obtain the fragility model, it is necessary to input the return period for which the hazard curve was derived, as well as the damage state of interest (slight, collapse, etc.), as described in the RMTK manual.

def calculate_convolution(hazard_curve,fragility_model,return_period,damage_state)


In [45]:
from rmtk.vulnerability.Convolution_module import Convolution

return_period = 50
damage_state = 'Collapse'
save = True

APD = Convolution.calculate_convolution(hazard_curve,fragility_model,return_period,damage_state)


1%
[0.021815, 0.034322500000000006, 0.046830000000000004, 0.071845]
2%
[0.020418, 0.0304645, 0.040511, 0.060604]
3%
[0.022193, 0.035367499999999996, 0.048542, 0.074891]
4%
[0.021648, 0.033862, 0.046076, 0.070504]
5%
[0.021948, 0.03469175, 0.0474355, 0.072923]
6%
[0.018117, 0.02410725, 0.0300975, 0.042078]
7%
[0.01898, 0.02649125, 0.0340025, 0.049025]
8%
[0.020176, 0.02979575, 0.0394155, 0.058655]
9%
[0.019343, 0.027493999999999998, 0.035644999999999996, 0.051947]
10%
[0.02104, 0.0321815, 0.043323, 0.065606]
11%
[0.019523, 0.0279905, 0.036458000000000004, 0.053393]
12%
[0.020571, 0.030885999999999997, 0.041201, 0.061831]
13%
[0.020137, 0.02968775, 0.0392385, 0.05834]
14%
[0.022588, 0.036458500000000005, 0.050329, 0.07807]
15%
[0.020935, 0.03189225, 0.0428495, 0.064764]
16%
[0.021835, 0.03437900000000001, 0.046923000000000006, 0.072011]
17%
[0.018722, 0.025779499999999997, 0.032837, 0.046952]
18%
[0.019444, 0.02777325, 0.0361025, 0.052761]
19%
[0.021865, 0.03446175, 0.047058499999999996, 0.072252]
20%
[0.017204, 0.02158425, 0.0259645, 0.034725]
21%
[0.022039, 0.034941749999999994, 0.0478445, 0.07365]
22%
[0.021893, 0.0345395, 0.047186, 0.072479]
23%
[0.02042, 0.0304705, 0.040521, 0.060622]
24%
[0.023139, 0.03798125, 0.052823499999999995, 0.082508]
25%
[0.021608, 0.03375125, 0.0458945, 0.070181]
26%
[0.021744, 0.03412675, 0.0465095, 0.071275]
27%
[0.022487, 0.03618025, 0.0498735, 0.07726]
28%
[0.019985, 0.029268, 0.038551, 0.057117]
29%
[0.018196, 0.0243255, 0.030455000000000003, 0.042714]
30%
[0.021775, 0.03421325, 0.0466515, 0.071528]
31%
[0.019797, 0.02874825, 0.0376995, 0.055602]
32%
[0.01976, 0.02864525, 0.0375305, 0.055301]
33%
[0.021717, 0.0340525, 0.046388, 0.071059]
34%
[0.018348, 0.024746249999999997, 0.0311445, 0.043941]
35%
[0.017144, 0.02141875, 0.0256935, 0.034243]
36%
[0.021134, 0.03244275, 0.0437515, 0.066369]
37%
[0.022675, 0.03669875, 0.050722500000000004, 0.07877]
38%
[0.020946, 0.03192225, 0.042898500000000006, 0.064851]
39%
[0.019559, 0.028090749999999998, 0.0366225, 0.053686]
40%
[0.018639, 0.0255495, 0.03246, 0.046281]
41%
[0.021205, 0.03263725, 0.0440695, 0.066934]
42%
[0.023368, 0.03861325, 0.0538585, 0.084349]
43%
[0.018866, 0.026176, 0.033486, 0.048106]
44%
[0.017858, 0.023391, 0.028924, 0.03999]
45%
[0.020648, 0.03110025, 0.0415525, 0.062457]
46%
[0.02253, 0.036297750000000004, 0.0500655, 0.077601]
47%
[0.020385, 0.030373749999999998, 0.040362499999999996, 0.06034]
48%
[0.018641, 0.0255545, 0.032468000000000004, 0.046295]
49%
[0.01831, 0.024641, 0.030972, 0.043634]
50%
[0.02239, 0.035912, 0.049434000000000006, 0.076478]
51%
[0.021735, 0.034103, 0.046471000000000005, 0.071207]
52%
[0.01777, 0.023148250000000002, 0.0285265, 0.039283]
53%
[0.022939, 0.03742925, 0.0519195, 0.0809]
54%
[0.020525, 0.030759750000000002, 0.0409945, 0.061464]
55%
[0.020339, 0.030245, 0.040151, 0.059963]
56%
[0.019546, 0.028054000000000003, 0.036562, 0.053578]
57%
[0.021572, 0.033651749999999994, 0.045731499999999994, 0.069891]
58%
[0.016987, 0.020985749999999997, 0.0249845, 0.032982]
59%
[0.023639, 0.0393615, 0.055083999999999994, 0.086529]
60%
[0.02227, 0.03557975000000001, 0.0488895, 0.075509]
61%
[0.021705, 0.034018999999999994, 0.046333, 0.070961]
62%
[0.021734, 0.034098500000000004, 0.046463000000000004, 0.071192]
63%
[0.019831, 0.028843, 0.037855, 0.055879]
64%
[0.020033, 0.0294, 0.038766999999999996, 0.057501]
65%
[0.018699, 0.0257145, 0.032729999999999995, 0.046761]
66%
[0.021515, 0.033494499999999996, 0.045474, 0.069433]
67%
[0.018726, 0.02579025, 0.032854499999999995, 0.046983]
68%
[0.020316, 0.03018175, 0.0400475, 0.059779]
69%
[0.021254, 0.032774250000000005, 0.0442945, 0.067335]
70%
[0.019646, 0.0283305, 0.037015, 0.054384]
71%
[0.020172, 0.029784249999999998, 0.0393965, 0.058621]
72%
[0.022119, 0.0351625, 0.048206, 0.074293]
73%
[0.019557, 0.028086, 0.036615, 0.053673]
74%
[0.022231, 0.035472500000000004, 0.048714, 0.075197]
75%
[0.020292, 0.030116749999999998, 0.0399415, 0.059591]
76%
[0.01788, 0.02345275, 0.0290255, 0.040171]
77%
[0.018533, 0.025257, 0.031980999999999996, 0.045429]
78%
[0.01752, 0.022457249999999998, 0.0273945, 0.037269]
79%
[0.01876, 0.025883749999999997, 0.033007499999999995, 0.047255]
80%
[0.018354, 0.02476225, 0.031170499999999997, 0.043987]
81%
[0.020515, 0.030731249999999998, 0.0409475, 0.06138]
82%
[0.018885, 0.0262295, 0.033574, 0.048263]
83%
[0.018237, 0.024438, 0.030639, 0.043041]
84%
[0.021364, 0.03307675, 0.044789499999999996, 0.068215]
85%
[0.018692, 0.025695250000000003, 0.032698500000000005, 0.046705]
86%
[0.017855, 0.023383750000000002, 0.0289125, 0.03997]
87%
[0.019354, 0.02752525, 0.0356965, 0.052039]
88%
[0.022364, 0.03584, 0.049316, 0.076268]
89%
[0.01912, 0.0268785, 0.034637, 0.050154]
90%
[0.01937, 0.027568999999999996, 0.035767999999999994, 0.052166]
91%
[0.019867, 0.028942000000000002, 0.038017, 0.056167]
92%
[0.018653, 0.02558825, 0.0325235, 0.046394]
93%
[0.02146, 0.03334275, 0.0452255, 0.068991]
94%
[0.021108, 0.032369499999999995, 0.043631, 0.066154]
95%
[0.022084, 0.03506625, 0.048048499999999994, 0.074013]
96%
[0.01828, 0.024557000000000002, 0.030834, 0.043388]
97%
[0.022221, 0.03544525, 0.048669500000000004, 0.075118]
98%
[0.021839, 0.034390000000000004, 0.046941, 0.072043]
99%
[0.020541, 0.030802999999999997, 0.041065, 0.061589]
100%
[0.017478, 0.0223425, 0.027207, 0.036936]

In [208]:
# -*- coding: utf-8 -*-
import os
import numpy
import math
import csv
import matplotlib.pyplot as plt
from rmtk.vulnerability.common import utils
from scipy.stats import lognorm

def read_hazard(input_file):    
#This function reads a hazard curve and stores it in a dictionary

    file = open(input_file)
    data = csv.reader(file)
    IMLs = []
    prob_exceedance = []
    IM_type = []
    for line in data:
        if line[0] == 'PGA':
            IM_type.append(line[0])
            for value in line[1:]:
                if isNumber(value):    
                    IMLs.append(float(value))
        if line[0] == 'Sa':
            IM_type.append(line[0])
            for value in line[1:]:
                if isNumber(value):
                    IMLs.append(float(value))
        if line[0] == 'PoE':
            for value in line[1:]:
                if isNumber(value):
                    prob_exceedance.append(float(value))
                
    #Store all the data in the dictionary
    hazard_curve = {'IMLs':None,'PoE':None,'IM_Type':None}
    hazard_curve['IMLs'] = IMLs
    hazard_curve['PoE'] = prob_exceedance
    hazard_curve['IM_Type'] = IM_type
      
    return hazard_curve

def isNumber(s):
    try:
        float(s)
        return True
    except ValueError:
        return False
    
def plot_hazard_curve(hazard_curve):
#This function plots the hazard curve

    IMLs = hazard_curve['IMLs']
    PoE = hazard_curve['PoE']
    plt.plot(IMLs,PoE,color='g',linewidth=2)
    plt.xlabel(hazard_curve['IM_Type'][0] + ' [g]',fontsize = 10)
    plt.ylabel('Probability of Exceedance',fontsize = 10)

In [209]:
input_file = 'Hazard_PGA.csv'
%matplotlib inline 
hazard_curve = read_hazard(input_file)
print(hazard_curve)
plot_hazard_curve(hazard_curve)


{'IM_Type': ['PGA'], 'PoE': [0.8417, 0.4357, 0.1229, 0.0459, 0.0205, 0.0103, 0.0056, 0.0032, 0.0019, 0.0012, 0.0008, 0.0003], 'IMLs': [0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2]}

In [210]:
# -*- coding: utf-8 -*-
import os
import numpy
import math
import csv
import matplotlib.pyplot as plt
from rmtk.vulnerability.common import utils
from scipy.stats import lognorm

def read_frag_model(input_file):    
#This function reads a fragility model and stores it in a dictionary

    damage_states = []
    model_type = []
    cent_value = []
    dispersion = []
    log_mean = []
    log_stdev = []
        
    file = open(input_file)
    data = file.readlines()
    line = data[0]
    line = line.strip().split(',')
    model_type = line[1]+' - '+line[2]
    
    file = open(input_file)
    data = csv.reader(file)
    data = [row for row in data]
    for iline in range(len(data)-1):
        line = data[iline+1]
        damage_states.append(line[0])
        cent_value.append(float(line[1]))
        dispersion.append(float(line[2]))
    
    if model_type == ' median -  dispersion':
        for ids in range(len(cent_value)):
            mean = math.log(cent_value[ids])
            stdev = dispersion[ids]
            log_mean.append(mean)
            log_stdev.append(stdev)
    elif model_type == ' mean of x -  cov of x':
        for ids in range(len(cent_value)):
            mu = cent.value[ids]
            cov = dispersion[ids]
            median = (mu**2)/math.sqrt(mu**2+cov**2)
            mean = math.log(median)
            stdev = math.sqrt(math.log(1+(cov**2)/mu**2))
            log_mean.append(mean)
            log_stdev.append(stdev)
    elif model_type == 'mean of ln(x) -  st. dev. of ln(x)':
        for ids in range(len(cent_value)):
            mean = cent_value[ids]
            stdev = line[ids]
            log_mean.append(mean)
            log_stdev.append(stdev)
                
    #Store all the data in the dictionary
    fragility_model =  {'Damage_states':None,'log_mean':None,'log_stdev':None}
    fragility_model['Damage_states'] = damage_states
    fragility_model['log_mean'] = log_mean
    fragility_model['log_stdev'] = log_stdev
       
    return fragility_model

In [211]:
fragility_function_file = 'MUR_H1.csv'
fragility_model = read_frag_model(fragility_function_file)
save = False
print(fragility_model)


{'Damage_states': ['Slight', 'Moderate', 'Extensive', 'Collapse'], 'log_stdev': [0.393208, 0.39895, 0.401788, 0.427387], 'log_mean': [-1.5282312670143712, -0.994517173291625, -0.7509097768813768, -0.47061501610459816]}

In [227]:
print(hazard_curve['IMLs'])
print(len(hazard_curve['IMLs']))


[0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2]
12

In [345]:
# -*- coding: utf-8 -*-
import os
import numpy
import math
import csv
import matplotlib.pyplot as plt
from rmtk.vulnerability.common import utils
from scipy.stats import lognorm

def calculate_convolution(hazard_curve,fragility_model,return_period,damage_state):
#This function calculates the annual probability of collapse given a hazard 
#curve and a fragility model

    rate_exceed = []
    segment_limit = []
    segment_width = []
    imls = []
    prob_exceed = []
    frag_curve = []
    damage_rate_dist = []
   
    for iIML in range(len(hazard_curve['IMLs'])):
        #1. Calculate the annual rate of exceedance
        imls.append(hazard_curve['IMLs'][iIML])
        prob_exceed.append(hazard_curve['PoE'][iIML])
        rate_exceed.append(-math.log(1-prob_exceed[iIML])/return_period)  
            
    #2. Divide the curve into segments and derive the rate of occurrence 
    #   of the associated central IML value
    for iIML in range(len(imls)-1):
        segment_limit.append((imls[iIML]+imls[iIML+1])/2)
        
    segment_width.append(segment_limit[0])
    for iDelta in range(len(imls)-2):
        k = iDelta + 1
        segment_width.append(segment_limit[k]-segment_limit[(k-1)])
    segment_width.append(imls[-1]-segment_limit[-1])
    
    #3. Read fragility curve of interest and estimate the probability of 
    #   reaching the DS of interest, given the values of IML
    ds = fragility_model['Damage_states'].index(damage_state)
    mean = fragility_model['log_mean'][ds]
    stdev = fragility_model['log_stdev'][ds]
    for iIML in range(len(imls)):
        frag_curve.append(lognorm.cdf(imls[iIML], stdev, loc=0, scale=math.exp(mean)))
    
    #4. Calculate the damage rate distribution for the given IML
    for iIML in range(len(imls)):
        damage_rate_dist.append(rate_exceed[iIML]*frag_curve[iIML]*segment_width[iIML])
    
    #5. Calculate the annual damage rate (ADR) and annual probability of damage (APD)
    ADR = sum(damage_rate_dist)
    APD = (math.exp(ADR*1)-1)/math.exp(ADR*1)
      
    return APD

    print(ADR)

    #print(rate_exceed)
    #print(frag_curve)
    #print(damage_rate_dist)
   
    
    #print(rate_exceed)

In [346]:
return_period = 50
damage_state = 'Collapse'

APD = calculate_convolution(hazard_curve,fragility_model,return_period,damage_state)
print(APD)


3.38339719991e-05