#INSYDE 2.0 MAIN

rm(list = ls())
library(truncnorm)
library("MASS")

setwd(".. working directory ..")

tnorm <- function(minimum, maximum, d_mean, d_std, size = 1){
  return(rtruncnorm(size, a = minimum, b = maximum, mean = d_mean, sd = d_std))
}

change_dist <- function(old, new){
  indices <- order(old)
  new_s <- sort(new)
  for (i in 1:length(old)){
    old[indices[i]] <- new_s[i]
  }
  return(old)
}

he_dist=rweibull(250000, shape = 1.25, scale = 1)
he_dist[he_dist < 0.01] <- runif(sum(he_dist < 0.01), 0.01, 0.03)
d_func <- 10 + 6 * sqrt(he_dist) * rnorm(250000, mean = 1, sd = 0.2)
d_dist <- rweibull(250000, shape = 1.25, scale = 36)
d_dist <- change_dist(d_func, d_dist)
d_dist[d_dist < 1] <- runif(length(d_dist[d_dist < 1]), min = 1, max = 2)
tt=d_dist/max(d_dist)
v_func <- (1.1-tt) * rnorm(250000, mean = 1, sd = 0.2*1/(1.1-tt))
v_dist <- rweibull(250000, shape = 1.15, scale = 0.35)
v_dist <- change_dist(v_func, v_dist)
v_dist[v_dist < 0.05] <- runif(length(v_dist[v_dist < 0.05]), min = 0.05, max = 0.1)
s_func <- 10 + 6 * sqrt(v_dist) * rnorm(250000, mean = 1, sd = 0.2)
s_dist <- runif(250000, min = 0.05, max = 0.2)
s_dist <- change_dist(s_func, s_dist)


default_std=0.2
copy = read.table(".. input file ..", sep=";", header=TRUE)  

# Load INSYDE
source("insyde_2_function.R")
survey <- read.csv('SurveyInput.csv', sep=';', header=TRUE)
IH_s <- survey$IH[!is.na(survey$IH)]
uncert <- 0

BH_dist <- read.csv("BH.csv",sep=',', header = TRUE)

up <- read.table("UnitPrice.txt")

cat('absoluteDamage', 'basementDamage', 'storeyDamage', 'floorDamage', 'boilerDamage', 'relativeDamage', 'Id', 'he', 'BT', 'Danno', "\n", file=".. output file ..", append=FALSE)

repetitions=1000
feats=c('he', 'v', 'd', 's', 'q', 'FA', 'BE', 'BA', 'EP', 'IH', 'BH', 'GL', 'NF', 'BT', 'BS', 'PD', 'PT', 'FL', 'YY', 'LM','Danno21') #lista parametri richiesti IA 'BP'
not_provided=!is.element(feats,colnames(copy)) 

if (any(not_provided)) {copy[,feats[not_provided]]=NA
	dati=as.data.frame(sapply(copy,rep.int,times=repetitions))
	cat('mancanti: ', feats[not_provided])
} else {dati=copy[rowSums(is.na(copy[feats]))==0,]
	dati=rbind(dati,sapply(copy[rowSums(is.na(copy[feats]))>0,],rep.int,times=repetitions))}

if (sum(is.na(dati$he))>0) {dati$he[is.na(dati$he)]=pmax(rweibull(sum(is.na(dati$he)), shape = 1.25, scale = 1),runif(sum(is.na(dati$he)), min = 0.01, max = 0.03))}
if (sum(is.na(dati$d))>0) dati$d[is.na(dati$d)]=quantile(d_dist,ecdf(d_func)(10 + 6 * sqrt(dati$he[is.na(dati$d)]) * rnorm(sum(is.na(dati$d)), mean = 1, sd = 0.2)))
if (sum(is.na(dati$v))>0) dati$v[is.na(dati$v)]=quantile(v_dist,ecdf(v_func)((1.1 - dati$d[is.na(dati$v)]/max(d_dist)) * rnorm(sum(is.na(dati$v)), mean = 1, sd = 0.2*1/(1.1-dati$d[is.na(dati$v)]/max(d_dist)))))
if (sum(is.na(dati$s))>0) dati$s[is.na(dati$s)]=quantile(s_dist,ecdf(s_func)(10 + 6 * sqrt(dati$v[is.na(dati$s)]) * rnorm(sum(is.na(dati$s)), mean = 1, sd = 0.2)))
if (sum(is.na(dati$q))>0) dati$q[is.na(dati$q)]=sample(c(0,1), size = sum(is.na(dati$q)), replace = TRUE, prob = c(0.5, 0.5))
if (sum(is.na(dati$BE))>0) dati$BE[is.na(dati$BE)]=sample(c(0,1), size = sum(is.na(dati$BE)), replace = TRUE, prob = c(0.2, 0.8)) 
if (sum(is.na(dati$BT))>0) dati$BT[is.na(dati$BT)]=sample(c(1, 2, 3, 4, 5), replace = TRUE, prob = c(0.535712, 0.10268, 0.129464, 0.10268, 0.129464), size = sum(is.na(dati$BT)))
if (sum(is.na(dati$FA))>0) {dati$FA[is.na(dati$FA) & dati$BT==1] = tnorm(min = 50, max = 320, d_mean = 160, d_std = 60, size = sum(is.na(dati$FA) & dati$BT==1))
	dati$FA[is.na(dati$FA) & dati$BT==2] = tnorm(min = 50, max = 160, d_mean = 110, d_std = 20, size = sum(is.na(dati$FA) & dati$BT==2))
	dati$FA[is.na(dati$FA) & dati$BT==3] = tnorm(min = 60, max = 160, d_mean = 95, d_std = 20, size = sum(is.na(dati$FA) & dati$BT==3))
	dati$FA[is.na(dati$FA) & dati$BT==4] = tnorm(min = 45, max = 140, d_mean = 85, d_std = 15, size = sum(is.na(dati$FA) & dati$BT==4))
	dati$FA[is.na(dati$FA) & dati$BT==5] = tnorm(min = 45, max = 120, d_mean = 85, d_std = 15, size = sum(is.na(dati$FA) & dati$BT==5))}
dati$M=1
dati$M[dati$BT==3]=as.integer(round(dati$FA[dati$BT==3]/160, digits = 0))
dati$M[dati$M==0]=1
dati$FA = dati$FA/dati$M
dati$EP = dati$EP/dati$M
dati$BA = dati$BA/dati$M


if (sum(is.na(dati$EP))>0) {dati$EP[is.na(dati$EP) & dati$BT==1] = 4.1*dati$FA[is.na(dati$EP) & dati$BT==1]^0.5 * rnorm(sum(is.na(dati$EP) & dati$BT==1), mean = 1, sd = default_std)
	dati$EP[is.na(dati$EP) & dati$BT==2] = 3*dati$FA[is.na(dati$EP) & dati$BT==2]^0.5 * rnorm(sum(is.na(dati$EP) & dati$BT==2), mean = 1, sd = default_std)
	dati$EP[is.na(dati$EP) & dati$BT==3] = 0.2885*dati$FA[is.na(dati$EP) & dati$BT==3] * rnorm(sum(is.na(dati$EP) & dati$BT==3), mean = 1, sd = default_std)
	dati$EP[is.na(dati$EP) & dati$BT==4] = 3*dati$FA[is.na(dati$EP) & dati$BT==4]^0.5 * rnorm(sum(is.na(dati$EP) & dati$BT==4), mean = 1, sd = default_std)
	dati$EP[is.na(dati$EP) & dati$BT==5] = 2*dati$FA[is.na(dati$EP) & dati$BT==5]^0.5 * rnorm(sum(is.na(dati$EP) & dati$BT==5), mean = 1, sd = default_std)
}
if (sum(is.na(dati$NF))>0) dati$NF[is.na(dati$NF)]= sample(c(1, 2, 3, 4), size = sum(is.na(dati$NF)), replace = TRUE, prob = c(0.092760181, 0.561085973, 0.248868778, 0.097285068))
dati$BA[dati$BE==0]=0
dati$BP[dati$BE==0]=0
dati$BH[dati$BE==0]=0
if (sum(is.na(dati$BA))>0) dati$BA[is.na(dati$BA)& dati$BE==1]= 0.5 * dati$FA[is.na(dati$BA)& dati$BE==1] * rnorm(sum(is.na(dati$BA)& dati$BE==1), mean = 1, sd = default_std)
if (sum(is.na(dati$BH))>0) dati$BH[is.na(dati$BH)& dati$BE==1]= sample(BH_dist$X0, size = sum(is.na(dati$BH)& dati$BE==1), replace = TRUE)
if (sum(is.na(dati$BP))>0) dati$BP[is.na(dati$BP)& dati$BE==1]= 4.2 * dati$FA[is.na(dati$BP)& dati$BE==1]^0.5 * rnorm(sum(is.na(dati$BP)& dati$BE==1), mean = 1, sd = default_std)
if (sum(is.na(dati$BS))>0) dati$BS[is.na(dati$BS)]= sample(c(1,2), size = sum(is.na(dati$BS)), replace = TRUE, prob = c(0.33, 0.67))

if (sum(is.na(dati$YY))>0) dati$YY[is.na(dati$YY)]= sample(c(1900,1930,1955,1965,1975,1985,1995,2003,2010), size= sum(is.na(dati$YY)), replace = TRUE , prob = c(0.16801854,0.106604867,0.139049826,0.17612978,0.163383546,0.095017381,0.070683662,0.046349942,0.034762457))

if (sum(is.na(dati$PD))>0) dati$PD[is.na(dati$PD) & dati$YY>=1990]= sample(c(1,2), size = sum(is.na(dati$PD) & dati$YY>=1990), replace = TRUE, prob = c(0.106, 0.894))
if (sum(is.na(dati$PD))>0) dati$PD[is.na(dati$PD) & dati$YY<1990]= sample(c(1,2), size = sum(is.na(dati$PD) & dati$YY<1990), replace = TRUE, prob = c(0.6, 0.4))
if (sum(is.na(dati$FL))>0) dati$FL[is.na(dati$FL)]= sample(c(0.8, 1, 1.2), size = sum(is.na(dati$FL)), replace = TRUE, prob = c(0.0505, 0.4245, 0.525))

if (sum(is.na(dati$PT))>0) dati$PT[is.na(dati$PT) & dati$YY>2000 & dati$FL> 1]= sample(c(1,2), size = sum(is.na(dati$PT) & dati$YY>2000 & dati$FL> 1), replace = TRUE, prob = c(0.2, 0.8))
if (sum(is.na(dati$PT))>0) dati$PT[is.na(dati$PT)]= sample(c(1,2), size = sum(is.na(dati$PT)), replace = TRUE, prob = c(0.8, 0.2))

if (sum(is.na(dati$LM))>0) dati$LM[is.na(dati$LM)]= sample(c(0.9, 1, 1.1), size = sum(is.na(dati$LM)), replace = TRUE, prob = c(0.13, 0.47, 0.4))
if (sum(is.na(dati$GL))>0) dati$GL[is.na(dati$GL)]= rnorm(sum(is.na(dati$GL)), mean = 0.1, sd = 0.09)
if (sum(is.na(dati$IH))>0) dati$IH[is.na(dati$IH)]= sample(IH_s, size=sum(is.na(dati$IH)), replace=TRUE) + tnorm(-0.15, 0.15, 0, 0.5, size=1) 

N <- nrow(dati)

repValData <- read.table("replacement_values.txt", header = TRUE)
ad  <- matrix(, N, 1)
s_d  <- matrix(, N, 1)
bo_d <- matrix(, N, 1)
f_d  <- matrix(, N, 1)
ba_d <- matrix(, N, 1)
rd   <- matrix(, N, 1)
RVN  <- matrix(, N, 1)
C    <- matrix(, N, 1)
R    <- matrix(, N, 1)
NS   <- matrix(, N, 1)
S    <- matrix(, N, 1)
F    <- matrix(, N, 1)
SY   <- matrix(, N, 1)

inp <- matrix(,N, 23)
C1 <- matrix(, N, 1)
C2 <- matrix(, N, 1)
C3 <- matrix(, N, 1)
C4 <- matrix(, N, 1)
R1 <- matrix(, N, 1)
R2 <- matrix(, N, 1)
R3 <- matrix(, N, 1)
R4 <- matrix(, N, 1)
R5 <- matrix(, N, 1)
R6 <- matrix(, N, 1)
R7 <- matrix(, N, 1)
R8 <- matrix(, N, 1)
R9 <- matrix(, N, 1)
R10<- matrix(, N, 1)
N1 <- matrix(, N, 1)
N2 <- matrix(, N, 1)
N3 <- matrix(, N, 1)
S1 <- matrix(, N, 1)
S2 <- matrix(, N, 1)
S3 <- matrix(, N, 1)
F1 <- matrix(, N, 1)
F2 <- matrix(, N, 1)
F3 <- matrix(, N, 1)
F4 <- matrix(, N, 1)
F5 <- matrix(, N, 1)
F6 <- matrix(, N, 1)
W1 <- matrix(, N, 1)
W2 <- matrix(, N, 1)
P1 <- matrix(, N, 1)
P2 <- matrix(, N, 1)
P3 <- matrix(, N, 1)
P4 <- matrix(, N, 1)
P5 <- matrix(, N, 1)

for (i in 1:N){

  M = dati$M[i]
  he  = dati$he[i]
  v   = dati$v[i]
  d   = dati$d[i]
  s   = dati$s[i]
  q   = dati$q[i]
  FA <- dati$FA[i]
  BE <- dati$BE[i]
  BA <- dati$BA[i]
  EP <- dati$EP[i]                            
  IH <- dati$IH[i]
  BH <- dati$BH[i]	  
  GL <- dati$GL[i]
  NF <- dati$NF[i]
  BT <- dati$BT[i]
  BS <- dati$BS[i]
  YY <- dati$YY[i]
  PD <- dati$PD[i]
  PT <- dati$PT[i]
  FL <- dati$FL[i]
  LM <- dati$LM[i]
  IP <- dati$IP[i]  
  BP <- dati$BP[i]
  IA <- FA * 0.9
  if (is.na(IP)||is.null(IP)) { 
	if (BT==1)        {IP <- 20.151 + 0.6254 * (FA * rnorm(1, mean = 1, sd = default_std)) # Internal perimeter (m)
	} else if (BT==2) {IP <- 20.119 + 0.6105 * (FA * rnorm(1, mean = 1, sd = default_std))
	} else if (BT==3) {IP <- 20.336 + 0.6576 * (FA * rnorm(1, mean = 1, sd = default_std))
	} else if (BT==4) {IP <- 9.709 + 0.6902 * (FA * rnorm(1, mean = 1, sd = default_std))
	} else if (BT==5) {IP <- 16.801 + 0.559 * (FA * rnorm(1, mean = 1, sd = default_std))}}
    
  repVal <- repValData[BS,BT]
  
  RVN[i,1] <- repVal * FA * NF
  inp[i,]=c(he, v, d, s, q, FA, M, IA, BE, BA, EP, IH, BH, GL, NF, BT, BS, PD, PT, FL, YY, LM, BP)
  modelOutput <- ComputeDamage(he, v, d, s, q, 
			  FA, IA, BE, BA, EP, IH, BH, GL, NF, BT, BS, PD, PT, FL, YY, LM, BP, IP, repVal, up, uncert) #
  ad[i,1] <- modelOutput$absDamage * M
  rd[i, 1]   <- modelOutput$relDamage * M
  s_d[i, 1]  <- modelOutput$storeyDamage * M
  bo_d[i, 1] <- modelOutput$boilerDamage * M
  f_d[i, 1]  <- modelOutput$floorDamage * M
  ba_d[i, 1] <- modelOutput$basementDamage * M
  C[i, 1]  <- modelOutput$cleanDamage * M
  R[i, 1]  <- modelOutput$removalDamage * M
  NS[i, 1] <- modelOutput$nonstructuralDamage * M
  S[i, 1]  <- modelOutput$structuralDamage * M
  F[i, 1]  <- modelOutput$finishingDamage * M
  SY[i, 1] <- modelOutput$systemsDamage * M
  
  C1[i, 1] <- modelOutput$C1 * M
  C2[i, 1] <- modelOutput$C2 * M
  C3[i, 1] <- modelOutput$C3 * M
  C4[i, 1] <- modelOutput$C4 * M
  R1[i, 1] <- modelOutput$R1 * M
  R2[i, 1] <- modelOutput$R2 * M
  R3[i, 1] <- modelOutput$R3 * M
  R4[i, 1] <- modelOutput$R4 * M
  R5[i, 1] <- modelOutput$R5 * M
  R6[i, 1] <- modelOutput$R6 * M
  R7[i, 1] <- modelOutput$R7 * M
  R8[i, 1] <- modelOutput$R8 * M
  R9[i, 1] <- modelOutput$R9 * M
  R10[i, 1] <- modelOutput$R10 * M
  N1[i, 1] <- modelOutput$N1 * M
  N2[i, 1] <- modelOutput$N2 * M
  N3[i, 1] <- modelOutput$N3 * M
  S1[i, 1] <- modelOutput$S1 * M
  S2[i, 1] <- modelOutput$S2 * M
  S3[i, 1] <- modelOutput$S3 * M
  F1[i, 1] <- modelOutput$F1 * M
  F2[i, 1] <- modelOutput$F2 * M
  F3[i, 1] <- modelOutput$F3 * M 
  F4[i, 1] <- modelOutput$F4 * M
  F5[i, 1] <- modelOutput$F5 * M
  F6[i, 1] <- modelOutput$F6 * M
  W1[i, 1] <- modelOutput$W1 * M
  W2[i, 1] <- modelOutput$W2 * M
  P1[i, 1] <- modelOutput$P1 * M
  P2[i, 1] <- modelOutput$P2 * M
  P3[i, 1] <- modelOutput$P3 * M
  P4[i, 1] <- modelOutput$P4 * M
  P5[i, 1] <- modelOutput$P5 * M
    cat(ad[i,1], ba_d[i,1], s_d[i,1], f_d[i,1], bo_d[i,1], rd[i,1], dati$ID[i], he, BT, dati$Danno21[i], "\n", file="damage_Lodi_1000M.txt", append=TRUE)
}

damage <- cbind(ad, ba_d, s_d, f_d, bo_d, rd, dati$ID, dati$he, dati$BT, dati$Danno21)
colnames(damage) <- c('absoluteDamage', 'basementDamage', 'storeyDamage', 'floorDamage', 'boilerDamage', 'relativeDamage', 'Id', 'he', 'BT','Danno21')

write.csv(damage,paste('.. output file 2 ..', sep=''))