#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jul 17 19:04:01 2021

@author: helencrowley
"""
import glob
import os
import pandas as pd
import numpy as np
from scipy import stats


#%% read list of typologies from capacity curves

current_dir=os.getcwd()
cap_in_files=glob.glob(os.path.join(os.path.join(current_dir,'../ESRM20_capacity_curves','*.csv')))
cap_in_files=sorted(cap_in_files)

#%% read fragility data

fragility_PGA = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='PGA')
fragility_SA03 = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='SA(0.3s)')
fragility_SA06 = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='SA(0.6s)')
fragility_SA10 = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='SA(1.0s)')
fragility_final = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='Final')
fragility_ShakeMap = pd.read_excel('../esrm20_fragility_various_IM_lognormal.xlsx',sheet_name='ShakeMap')

#%%  damage-loss model for loss of life

dam2loss= pd.read_excel('fatality_damage_model_ESRM20.xlsx')
p_fatcol_com = dam2loss['P_lethal-building | DS4']
coll_factor = dam2loss['Collapse_factor']
p_entrap_fatcol_day = dam2loss['P_entrapment_day']
p_entrap_fatcol_night = dam2loss['P_entrapment_night']
p_entrap_fatcol = (p_entrap_fatcol_day+p_entrap_fatcol_night)/2;
fat_entrap = dam2loss['P_loss-life | entrapment']
fat_ratio_day = p_fatcol_com * coll_factor * p_entrap_fatcol_day * fat_entrap
fat_ratio_night = p_fatcol_com * coll_factor * p_entrap_fatcol_night * fat_entrap
fat_ratio_ave = p_fatcol_com * coll_factor * p_entrap_fatcol * fat_entrap # average models over 24 hour period
#fat_ratio_ave = fat_ratio_day # uncomment when running models for day
#fat_ratio_ave = fat_ratio_night # uncomment when running models for night

#%% calculate vulnerability for all IMs

vuln_PGA = np.zeros([len(cap_in_files),50])
vuln_SA03 = np.zeros([len(cap_in_files),50])
vuln_SA06 = np.zeros([len(cap_in_files),50])
vuln_SA10 = np.zeros([len(cap_in_files),50])
Class = np.empty([len(cap_in_files)],dtype='<U25')

for k in range(len(cap_in_files)):
    cap_file=cap_in_files[k]
    building_class=cap_file.split('/')[-1].split('\\')[-1].split('.')[0][10:]
    Class[k] = building_class
    x_vec = np.linspace(np.log(0.05),np.log(5),endpoint=True)
    #PGA
    med4 =  fragility_PGA[fragility_PGA['Typology']==building_class]['Median_DS4']
    beta =  fragility_PGA[fragility_PGA['Typology']==building_class]['Beta']
    
    probability_damage_state_PGA4=np.zeros(len(x_vec))
    for i in range(len(x_vec)):       
        probability_damage_state_PGA4[i]=stats.norm.cdf(x_vec[i],loc=np.log(med4),scale=beta) 
        
    vuln_PGA[k,:]= fat_ratio_ave[k]*probability_damage_state_PGA4
    
    
    #SA03
    med4 =  fragility_SA03[fragility_SA03['Typology']==building_class]['Median_DS4']
    beta =  fragility_SA03[fragility_SA03['Typology']==building_class]['Beta']
    
    probability_damage_state_SA034=np.zeros(len(x_vec))
    for i in range(len(x_vec)):       
        probability_damage_state_SA034[i]=stats.norm.cdf(x_vec[i],loc=np.log(med4),scale=beta) 
        
    vuln_SA03[k,:]= fat_ratio_ave[k]*probability_damage_state_SA034
    
    #SA06
    med4 =  fragility_SA06[fragility_SA06['Typology']==building_class]['Median_DS4']
    beta =  fragility_SA06[fragility_SA06['Typology']==building_class]['Beta']
    
    probability_damage_state_SA064=np.zeros(len(x_vec))
    for i in range(len(x_vec)):       
        probability_damage_state_SA064[i]=stats.norm.cdf(x_vec[i],loc=np.log(med4),scale=beta) 
        
    vuln_SA06[k,:]= fat_ratio_ave[k]*probability_damage_state_SA064
            
    #SA10
    med4 =  fragility_SA10[fragility_SA10['Typology']==building_class]['Median_DS4']
    beta =  fragility_SA10[fragility_SA10['Typology']==building_class]['Beta']
    
    probability_damage_state_SA104=np.zeros(len(x_vec))
    for i in range(len(x_vec)):       
        probability_damage_state_SA104[i]=stats.norm.cdf(x_vec[i],loc=np.log(med4),scale=beta) 
        
    vuln_SA10[k,:]= fat_ratio_ave[k]*probability_damage_state_SA104
 
    
#%% write Excel sheets for viewer

out_PGA = pd.DataFrame(vuln_PGA)
out_SA03 = pd.DataFrame(vuln_SA03)
out_SA06 = pd.DataFrame(vuln_SA06)
out_SA10 = pd.DataFrame(vuln_SA10)
imls_PGA = pd.DataFrame(np.exp(x_vec), columns=['PGA'])
imls_SA03 = pd.DataFrame(np.exp(x_vec), columns=['SA(0.3)'])
imls_SA06 = pd.DataFrame(np.exp(x_vec), columns=['SA(0.6)'])
imls_SA10 = pd.DataFrame(np.exp(x_vec), columns=['SA(1.0)'])

typ = pd.DataFrame(Class)


writer=pd.ExcelWriter('../esrm20_vulnerability_loss-of-life.xlsx', engine='xlsxwriter')    
imls_PGA.to_excel(writer,sheet_name='PGA', header=True, index=False, startrow=0, startcol=0)
typ.T.to_excel(writer,sheet_name='PGA', header=False, index=False, startrow=0, startcol=1)
out_PGA.T.to_excel(writer,sheet_name='PGA', header=False, index=False, startrow=1, startcol=1)
imls_SA03.to_excel(writer,sheet_name='SA(0.3s)', header=True, index=False, startrow=0, startcol=0)
typ.T.to_excel(writer,sheet_name='SA(0.3s)', header=False, index=False, startrow=0, startcol=1)
out_SA03.T.to_excel(writer,sheet_name='SA(0.3s)', header=False, index=False, startrow=1, startcol=1)
imls_SA06.to_excel(writer,sheet_name='SA(0.6s)', header=True, index=False, startrow=0, startcol=0)
typ.T.to_excel(writer,sheet_name='SA(0.6s)', header=False, index=False, startrow=0, startcol=1)
out_SA06.T.to_excel(writer,sheet_name='SA(0.6s)', header=False, index=False, startrow=1, startcol=1)
imls_SA10.to_excel(writer,sheet_name='SA(1.0s)', header=True, index=False, startrow=0, startcol=0)
typ.T.to_excel(writer,sheet_name='SA(1.0s)', header=False, index=False, startrow=0, startcol=1)
out_SA10.T.to_excel(writer,sheet_name='SA(1.0s)', header=False, index=False, startrow=1, startcol=1)

writer.save()

#%% write csv files for viewer

(pd.read_excel('../esrm20_vulnerability_loss-of-life.xlsx',sheet_name='PGA')).to_csv('../viewer/ESRM20_vulnerability_loss-of-life_PGA.csv', index=False)
(pd.read_excel('../esrm20_vulnerability_loss-of-life.xlsx',sheet_name='SA(0.3s)')).to_csv('../viewer/ESRM20_vulnerability_loss-of-life_SA03.csv', index=False)
(pd.read_excel('../esrm20_vulnerability_loss-of-life.xlsx',sheet_name='SA(0.6s)')).to_csv('../viewer/ESRM20_vulnerability_loss-of-life_SA06.csv', index=False)
(pd.read_excel('../esrm20_vulnerability_loss-of-life.xlsx',sheet_name='SA(1.0s)')).to_csv('../viewer/ESRM20_vulnerability_loss-of-life_SA10.csv', index=False)


#%% xml code generation for loss of life, most efficient intensity measures

from xml.dom import minidom

doc = minidom.Document()

nrml=doc.createElement('nrml')
nrml.setAttribute('xmlns','http://openquake.org/xmlns/nrml/0.5')
doc.appendChild(nrml)

vulnerabilityModel = doc.createElement('vulnerabilityModel')
vulnerabilityModel.setAttribute('id','vm-global')
vulnerabilityModel.setAttribute('assetCategory','buildings')
vulnerabilityModel.setAttribute('lossCategory','occupants')
nrml.appendChild(vulnerabilityModel)

description=doc.createElement('description')
text_d=doc.createTextNode('ESRM20 loss of life vulnerability model for Europe in terms of various SA')
description.appendChild(text_d)
vulnerabilityModel.appendChild(description)


for x in range(len(cap_in_files)):
    
    imt = fragility_final.IMT.iloc[x]
          
    imls2=np.exp(x_vec)
    if imt == 'PGA':
        vul_mean= vuln_PGA[x]
    elif imt == 'SA(0.3)':
        vul_mean = vuln_SA03[x]
    elif imt == 'SA(0.6)':
        vul_mean = vuln_SA06[x]
    else:
        vul_mean = vuln_SA10[x]
    
    s=" "
    imls_chain=s.join([str(i) for i in np.ravel(np.around(imls2,9))])
    vul_mean_chain=s.join([str(i) for i in np.ravel(np.around(vul_mean,9))])
           
    vulnerabilityFunction = doc.createElement('vulnerabilityFunction')
    vulnerabilityFunction.setAttribute('id',fragility_final['Typology'].iloc[x])
    
    vulnerabilityFunction.setAttribute('dist','LN')
    vulnerabilityModel.appendChild(vulnerabilityFunction)
    
    imls = doc.createElement('imls')
    text = doc.createTextNode(imls_chain)
    imls.appendChild(text)
    imls.setAttribute('imt', imt)
    vulnerabilityFunction.appendChild(imls)
    
    meanLRs = doc.createElement('meanLRs')
    text = doc.createTextNode(vul_mean_chain)
    meanLRs.appendChild(text)
    vulnerabilityFunction.appendChild(meanLRs)
       
                
    prob=np.zeros((1,50))
    s=" "
    prob_chain=s.join([str(i) for i in np.ravel(np.around(prob,6))])
    covLRs = doc.createElement('covLRs')
    text = doc.createTextNode(prob_chain)
    covLRs.appendChild(text)
    vulnerabilityFunction.appendChild(covLRs)

        
    xml_str = doc.toprettyxml(indent="")

    
    with open("../../esrm20_openquake/Vulnerability/vulnerability_loss-of-life_ESRM20_VariousIM.xml", "w") as f:
        f.write(xml_str)
  #  with open("../../esrm20_scenario_tests/ruptures/vulnerability/vulnerability_loss-of-life_ESRM20_VariousIM.xml", "w") as f:
   #     f.write(xml_str)
    
           
        
# #%% xml code generation for loss of life, ShakeMap intensity measures

# from xml.dom import minidom

# doc = minidom.Document()

# nrml=doc.createElement('nrml')
# nrml.setAttribute('xmlns','http://openquake.org/xmlns/nrml/0.5')
# doc.appendChild(nrml)

# vulnerabilityModel = doc.createElement('vulnerabilityModel')
# vulnerabilityModel.setAttribute('id','vm-global')
# vulnerabilityModel.setAttribute('assetCategory','buildings')
# vulnerabilityModel.setAttribute('lossCategory','occupants')
# nrml.appendChild(vulnerabilityModel)

# description=doc.createElement('description')
# text_d=doc.createTextNode('ESHM20 loss of life vulnerability model for Europe in terms of ShakeMap intensity measures')
# description.appendChild(text_d)
# vulnerabilityModel.appendChild(description)


# for x in range(len(cap_in_files)):
    
#     imt = fragility_ShakeMap.IMT.iloc[x]
          
#     imls2=np.exp(x_vec)
#     if imt == 'PGA':
#         vul_mean= vuln_PGA[x]
#     elif imt == 'SA(0.3)':
#         vul_mean = vuln_SA03[x]
#     elif imt == 'SA(0.6)':
#         vul_mean = vuln_SA06[x]
#     else:
#         vul_mean = vuln_SA10[x]
    
#     s=" "
#     imls_chain=s.join([str(i) for i in imls2])
#     vul_mean_chain=s.join([str(i) for i in vul_mean])
           
#     vulnerabilityFunction = doc.createElement('vulnerabilityFunction')
#     vulnerabilityFunction.setAttribute('id',fragility_final['Typology'].iloc[x])
    
#     vulnerabilityFunction.setAttribute('dist','LN')
#     vulnerabilityModel.appendChild(vulnerabilityFunction)
    
#     imls = doc.createElement('imls')
#     text = doc.createTextNode(imls_chain)
#     imls.appendChild(text)
#     imls.setAttribute('imt', imt)
#     vulnerabilityFunction.appendChild(imls)
    
#     meanLRs = doc.createElement('meanLRs')
#     text = doc.createTextNode(vul_mean_chain)
#     meanLRs.appendChild(text)
#     vulnerabilityFunction.appendChild(meanLRs)
       
                
#     prob=np.zeros((1,50))
#     s=" "
#     prob_chain=s.join([str(i) for i in np.ravel(np.around(prob,6))])
#     covLRs = doc.createElement('covLRs')
#     text = doc.createTextNode(prob_chain)
#     covLRs.appendChild(text)
#     vulnerabilityFunction.appendChild(covLRs)

        
#     xml_str = doc.toprettyxml(indent="")

    
#     with open("../../esrm20_scenario_tests/shakemaps/shakemap_vulnerability/vulnerability_loss-of-life_ESRM20_ShakeMapIM.xml", "w") as f:
#         f.write(xml_str)
        
        