#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 17 15:08:55 2020

@author: jcgomez
"""
#Script to calculate an inter-scheme compatibility matrix #as proposed by Gomez-Zapata et al.,2021
#the method is descibed in Section 3.2. Probabilistic inter-scheme compatibility matrix
#in the paper submitted to Bulletin of Earthquake Engineering in January 24, 2021.
#and available in preprint format in:  https://doi.org/10.21203/rs.3.rs-178120/v1
#In this example, SARA (source) and HAZUS (target) schemes and both represent the residential buildings classes of Valparaiso
import pandas as pd
import numpy as np
import itertools
import matplotlib.pyplot as plt
import seaborn as sns
import os

os.chdir('..')
from Exposure_fuzzy_scores.lib_fuzzy_scoring import scoring

def compute_tax_conv(tab_source,tab_target,input_tax,p_inputdata): 
    res1 = tab_target.transpose().mul(tab_source[input_tax]*p_inputdata).transpose()
    res2 = np.sum(res1/tab_source.mean()[input_tax])
    return(res2)

dict_file_SARA = "Exposure_fuzzy_scores/fuzzy_scored_schemes/Lima_SARA.json"
dict_file_Medina = "Exposure_fuzzy_scores/fuzzy_scored_schemes/Lima_MEDINA.json"   

SARA_cd = scoring.importClassDict(dict_file_SARA)
Medina_cd = scoring.importClassDict(dict_file_Medina)
from collections import OrderedDict
print(SARA_cd['weights'])
print(Medina_cd['weights'])

###########################################################################################################
###Modify weighting scheme (customizing if necesarry)
###########################################################################################################

Medina_cd['weights'] = OrderedDict([('mat_type', 0.5), 
                                      ('mat_tech', 0.1), 
                                      #('mat_prop', 0.0), 
                                      ('llrs', 0.1), 
                                      ('llrs_duct', 0.00), 
                                      #('floor_mat', 0.1),
                                      ('height_2', 0.3)])

SARA_cd['weights'] = OrderedDict([('mat_type', 0.5), 
                                      ('mat_tech', 0.1), 
                                      #('mat_prop', 0.05), 
                                      ('llrs', 0.1), 
                                      ('llrs_duct', 0.0), 
                                      #('floor_mat', 0.1),
                                      ('height_2', 0.3)])

print(sum(Medina_cd['weights'].values()))
print(sum(SARA_cd['weights'].values()))

###########################################################################################################
###Generate an exhaustive collection of attribute values
###########################################################################################################
#Please note that in this case, only the attributes types that comonly belong to both schemes are kept.
#Other attributes types and values that might belong to other building classes are provided in a commented (unused) manner.

def expand_grid(dct):
    rows = itertools.product(*dct.values())
    return pd.DataFrame.from_records(rows, columns=dct.keys())
    #return pd.DataFrame.from_records(rows, columns=dct.keys()).drop_duplicates()

df = expand_grid(
{'llrs' : ['L99', 'LDUAL', 'LFINF', 'LFM','LWAL', 'L99', 'LFLS', 'LFBR', 'PC_LWAL', 'LPB', 'LN'],

'llrs_duct' : ['DU99','DUC', 'DNO'],

 'mat_type' : ['M99','S', 'CR','ME', 'MUR','W','MR','MCF', 'SRC','EU' , 'C99', 'E99'],

 'mat_tech' : ['MATT99','ADO','CB99','CBH', 'CBS','CIP','CL99', 'CLBLH', 'CLBRH', 'CLBRS', 'ET99', 'ETC', 'ETO',                
                'ETR','MO','MUN99', 'S99','SL', 'SO', 'SR', 'ST99','STDRE','STRUB','WWD', 'WLI', 'WHE', 'WS','WO','WBB', 'PC'],
 
 #'mat_prop' : ['SC99','WEL','RIV','BOL','MO99','MON','MOM','MOL','MOC','MOCL','SP99','SPLI','SPSA','SPTU','SPSL',
 #             'SPGR','SPBA','SPO','MATP99'],

 #'str_irreg_type' : ['IR99','SOS', 'IRRE', 'SOS' ],

 'height_2': ['H99','1','2','3','4 - 7','8 - 19'],  

 #'floor_mat': ['F99','FC','FE','FN','FW'],

#'floor_type' :['FT99','FC1','FC2','FC3','FC4','FC99','FE99','FW1', 'FW2','FW3','FW4','FW99',
#'FM1', 'FM2', 'FM3', 'FM99', 'FME1','FME2','FME3','FME4']

#'floor_conn': ['FWC99', 'FWCN','FWCP'],

#'roofSystTyp': ['RST99','RC1','RC2','RC3','RC4','RC99','RME1','RME2','RME3','RME4','RME99','RWO','RWO1','RWO2','RWO3','RWO4','RWO99'] ,#only in hazus
#'roofSystMat': ['R99', 'RC', 'RC99','RME','RME99', 'RSM99','RWO','RWO99'], #only in hazus
#'roofCovMat': ['RMT99','RMN','RMT6', 'RMT7','RMTO'], #only in hazus
#'nonStrcExW' : ['EW99', 'EWC', 'EWG', 'EWMA', 'EWME', 'EWO', 'EWSL', 'EWW']
} )

df['height_1']='NO_DATA'
df['object_id']=df.index

input_set = df.sample(len(df))
##input_set

###########################################################################################################
#SARA Model WARNING: the following can take a loong time !
###########################################################################################################

###compute scores
expo_df_SARA,scores_SARA = scoring.compute_expo_model_fuzzy(input_set,SARA_cd['weights'], SARA_cd['definition'], SARA_cd['fuzzy_values'],bDefuz=True)
#normalization - only the scores > 0 are retained
sc_tmp = scores_SARA[SARA_cd['classes']]*(scores_SARA[SARA_cd['classes']]>0).astype(int)
proba_SARA= sc_tmp.div(sc_tmp.sum(axis=1),axis=0)

###########################################################################################################
#Medina Model ### WARNING: the following can take a loong time !
##########################################################################################################

###compute scores
expo_df_Medina,scores_Medina = scoring.compute_expo_model_fuzzy(input_set,Medina_cd['weights'], Medina_cd['definition'], Medina_cd['fuzzy_values'],bDefuz=True)
#normalization - only the scores > 0 are retained
sc_tmp = scores_Medina[Medina_cd['classes']]*(scores_Medina[Medina_cd['classes']]>0).astype(int)
proba_Medina = sc_tmp.div(sc_tmp.sum(axis=1),axis=0)

#compile a list of records for which at list a positive score was found
ok_SARA = ((proba_SARA.sum(1)>0.99) & (proba_SARA.sum(1)<1.01))
ok_Medina = ((proba_Medina.sum(1)>0.99) & (proba_Medina.sum(1)<1.01))

#filter the dataframes keeping only the good records
proba_Medina = proba_Medina[ok_Medina & ok_SARA]
proba_SARA = proba_SARA[ok_Medina & ok_SARA]

n_ok_records = len(proba_Medina)


####Create conversion matrix  WARNING: the following can take a loong time !
p_od = 1/n_ok_records
SARA_Medina_convmatrix = pd.DataFrame([compute_tax_conv(proba_SARA,proba_Medina,input_tax,p_inputdata=p_od) for input_tax in SARA_cd['classes']])
SARA_Medina_convmatrix.index = SARA_cd['classes']

pd.options.display.width=100
print(SARA_Medina_convmatrix.round(3))

#check the normalisation of the conversion matrix
SARA_Medina_convmatrix.sum(1)

#note= for the UNK (unknown) building type, there is no value > 0 obtained.
##the value is then, in a post-process= 1/6 for all of Medina classes

#SAVE
#SARA_Medina_convmatrix.to_csv('include_the_name_you prefer.csv')
SARA_Medina_convmatrix.to_csv("SARA-Medina_LIMA_conv_matrix_modweights_v1.csv")

###########################################################################################################
###PLOT
###short script used to generate figure 6b in the preprint Gomez-Zapata et al.,2021 (https://doi.org/10.21203/rs.3.rs-178120/v1)
###########################################################################################################
matrix = pd.read_csv('SARA-Medina_LIMA_conv_matrix_modweights_v1.csv')
tax_name= matrix.iloc[:,0]
matrix = matrix.drop(matrix.columns[0], axis=1)

cmap = 'coolwarm'
fontsize_pt = 10
dpi = 72.27
matrix_height_pt = fontsize_pt * matrix.shape[0]
matrix_height_in = matrix_height_pt / dpi
top_margin = 0.06  # in percentage of the figure height
bottom_margin = 0.06 # in percentage of the figure height
figure_height = matrix_height_in / (1 - top_margin - bottom_margin)

fig, ax = plt.subplots(
        figsize=(3,figure_height*1.3), 
        gridspec_kw=dict(top=1-top_margin, bottom=bottom_margin,))

ax = sns.heatmap(matrix, cmap=cmap, ax=ax, yticklabels=tax_name, vmin=0, vmax=1) 
# Set x-axis label
plt.ylabel('SARA')
# Set y-axis label
plt.xlabel('Medina scheme')
###########################################################################################################
#SAVE FIG
###########################################################################################################
#ax.figure.savefig('Fig_Matrix_Lima_SARA-Medina_matrix.png')


    
