## Figures for manuscript "Representing storylines with causal networks 
## to support decision making: framework and example"

library(tidyverse, warn.conflicts = FALSE) # ggplot
library(bnlearn, warn.conflicts = FALSE) # bn
library(gRain, warn.conflicts = FALSE) # for exact inference 
library(Rgraphviz, warn.conflicts = FALSE) # additional options for figures
library(readxl)
library(ggh4x)
library(cowplot)
library(magick)


### Setting up model

# data # pre-processed into rds files
capital <- readRDS("capital_under_combs.rds") # main data 
comb <- readRDS("combinations.rds")

# combining the data into one dataframe
n_scenario <- sum(capital$comb_number==0) # number of scenarios
comb1 <- comb[rep(1:nrow(comb), each = n_scenario), ]
names(comb1)[1] <- "scenario"
capital <- cbind(capital, comb1)

# cyclones to consider
cycl_year_1 <- c("Carlos", "Enawo", "Harvey", "Irma", "Maria", "Ophelia")
cycl_year_2 <- c("Berguitta", "Fakir", "Isaac", 
                 "Helene", "Kirk", "Leslie", "Ava")

# cleaning dataframe
capital <- capital[c("capy1", "capy2", "haz_incr", "exp_incr", "cap_incr", 
                     cycl_year_1, cycl_year_2)]
capital$haz_incr <- as.factor(capital$haz_incr)
capital$exp_incr <- as.factor(capital$exp_incr)
capital$cap_incr <- as.factor(capital$cap_incr)



# create DAG for inference
nodes_year_1 <- paste0("[", paste0(cycl_year_1, collapse="]["), "]")
nodes_year_2 <- paste0("[", paste0(cycl_year_2, collapse="]["), "]")
parents_year_1 <- paste0(cycl_year_1, collapse=":")
parents_year_2 <- paste0(cycl_year_2, collapse=":")

dag <- model2network(paste0("[haz_incr][exp_incr][cap_incr]",
                            nodes_year_1, nodes_year_2,
                            "[capy1|haz_incr:exp_incr:cap_incr:",
                            parents_year_1, "]",
                            "[capy2|haz_incr:exp_incr:cap_incr:capy1:",
                            parents_year_2, "]"))


# learn parameters - uses MLE
bn_model <- bn.fit(dag, data = capital)









### Figure 3: Show distributions

data <- 
  tibble("dist" = rep("increase", 5),
         "Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
         "Probability" = c(0, 0.1, 0.2, 0.3, 0.4)) 
data <- rbind(data, 
              tibble("dist" = rep("flat", 5), 
                     "Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
                     "Probability" = c(0.2, 0.2, 0.2, 0.2, 0.2)))
data <- rbind(data, 
              tibble("dist" = rep("decrease", 5), 
                     "Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
                     "Probability" = c(0.4, 0.3, 0.2, 0.1, 0)))

ggplot(data = data, mapping = aes(x = `Hazard increase`, y = `Probability`)) + 
  geom_bar(stat = "identity") +
  facet_wrap(~ dist) +
  theme(text = element_text(size = 20), strip.text.x = element_blank())      





### Figure 5: Matrix of possible distributions (Figure 4 uses part of this)

## prior inputs
# hazard
low_increase <- array(c(0.4, 0.3, 0.2, 0.1, 0), 
                      dimnames = list(c(1, 3, 5, 7, 9)))
uniform_increase <- array(c(0.2, 0.2, 0.2, 0.2, 0.2), 
                          dimnames = list(c(1, 3, 5, 7, 9)))
high_increase <- array(c(0, 0.1, 0.2, 0.3, 0.4), 
                       dimnames = list(c(1, 3, 5, 7, 9)))
haz_dist <- list(low_increase, uniform_increase, high_increase)
haz_name <- c("low", "uniform", "high") # names for each pattern

# gdp
low_increase <- array(c(0.4, 0.3, 0.2, 0.1, 0), 
                      dimnames = list(c(1, 6, 11, 16, 21)))
uniform_increase <- array(c(0.2, 0.2, 0.2, 0.2, 0.2), 
                          dimnames = list(c(1, 6, 11, 16, 21)))
high_increase <- array(c(0, 0.1, 0.2, 0.3, 0.4), 
                       dimnames = list(c(1, 6, 11, 16, 21)))
gdp_dist <- list(low_increase, uniform_increase, high_increase)
gdp_name <- c("low", "uniform", "high") # names for each pattern


cap_lab <- c("0%", "30%", 
             "60%", "90%", 
             "120%", "150%")
names(cap_lab) <- c("0", "30", "60", "90", "120", "150")
haz_lab <- c("Low hazard increase more likely",
             "Uniform distribution for hazard increase",
             "High hazard increase more likely")
names(haz_lab) <- c("low", "uniform", "high")
gdp_lab <- c("GDP low",
             "GDP uniform",
             "GDP high")
names(gdp_lab) <- c("low", "uniform", "high")

# years to consider 
years <- c("capy1", "capy2")
names(years) <-  c(2017, 2018)

# plot function for figures, using ggplot

dist_all <- data.frame(val = NULL, haz = NULL, gdp = NULL, cap = NULL)

for (haz in 1:3){
  for (gdp in 1:3){
    for (cap in names(cap_lab)){
      
        
      bn_model$haz_incr <- haz_dist[[haz]]
      bn_model$exp_incr <- gdp_dist[[gdp]]
      
      # Realized/unrealized cyclones and capital value input
      evidence <- vector(mode = "list")
      for (cyc in c(cycl_year_1, cycl_year_2)) {
          evidence[[cyc]] <- "1"
      }
      
      evidence["cap_incr"] <- cap
      
      # Monte Carlo
      val <- cpdist(bn_model, 
                    nodes = "capy2",
                    evidence = evidence,
                    method = "lw")[["capy2"]]
      # add results to dataframe
      dist_all <-
        rbind(dist_all,
              cbind(val = val,
                    cap = cap,
                    haz = haz_name[haz],
                    gdp = gdp_name[gdp]))
    }
  }
}


# fix data type
dist_all$val <- as.numeric(dist_all$val)
dist_all$haz <- factor(dist_all$haz,
                       levels = c("low", "uniform", "high"))
dist_all$gdp <- factor(dist_all$gdp,
                       levels = c("low", "uniform", "high"))
dist_all$cap <- factor(dist_all$cap,
                       levels = c("0", "30", "60", "90", "120", "150"))

# total number of simulations for each case
n_each <- nrow(dist_all[dist_all$haz == "high" & 
                          dist_all$gdp == "high" & 
                          dist_all$cap == "0", ])



dist_all$xaxis <- "Hazard increase and GDP increase"
dist_all$yaxis <- "Capital increase"
dist_all$xaxis <- factor(dist_all$xaxis,
                         levels = c("Hazard increase and GDP increase"))
dist_all$yaxis <- factor(dist_all$yaxis,
                         levels = c("Capital increase"))
dist_all$empty <- ""



# User number 
dist_all[dist_all$haz == "high" & dist_all$gdp == "low", ]$empty <- "(User 1)"
dist_all[dist_all$haz == "low" & dist_all$gdp == "high", ]$empty <- "(User 2)"
dist_all[dist_all$haz == "uniform" & 
           dist_all$gdp == "uniform", ]$empty <- "(User 3)"


numbers <- dist_all %>% 
  group_by(haz, gdp, cap, xaxis, yaxis, empty) %>% 
  summarize(prop_p = sum(val >= 0) / n_each,
            prop_n = 1 - sum(val >= 0) / n_each)
numbers$haz <- factor(numbers$haz,
                  levels = c("low", "uniform", "high"))
numbers$gdp <- factor(numbers$gdp,
                      levels = c("low", "uniform", "high"))
numbers$cap <- factor(numbers$cap,
                      levels = c("0", "30", "60", "90", "120", "150"))
numbers$xaxis <- factor(numbers$xaxis,
                        levels = c("Hazard increase and GDP increase"))
numbers$yaxis <- factor(numbers$yaxis,
                        levels = c("Capital increase"))
numbers$empty <- factor(numbers$empty,
                        levels = c("", "(User 1)", "(User 2)", "(User 3)"))


# set plotting region

if (min(dist_all$val) > 1800) {
  lim <- 3000
} else {
  lim <- 1800
}



# add label for right strip
p <- 
  ggplot(data = dist_all, mapping = aes(x = val, y = stat(count / n_each))) +
  geom_histogram(data = subset(dist_all, val < 0), fill = "red", col="grey",
                 boundary = 0) +
  geom_histogram(data = subset(dist_all, val >= 0), fill = "blue", col="grey", 
                 boundary = 0) +
  facet_nested(rows = c(vars(yaxis), vars(cap)), 
               cols = c(vars(xaxis), vars(haz), vars(gdp), vars(empty)), 
    # yaxis + cap ~ xaxis + haz + gdp,
             labeller = labeller(haz = haz_lab,
                                 gdp = gdp_lab,
                                 cap = cap_lab)) +
  geom_vline(xintercept=0, linetype = "longdash") +
  xlim(-lim, lim) +
  xlab("Capital Level after second year (Million EUR)") + 
  ylab("Probability") +
  theme(plot.title = element_text(hjust=0.5),
        strip.text = element_text(size = 15),
        panel.spacing = unit(4, "pt"),
        axis.text = element_text(size = 12),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 15))

p <- p +
  geom_text(data = numbers, 
            size = 5,
            mapping = aes(x = 1000, y = 0.3, 
                          label = paste0(round(prop_p * 100, 0), "%"))) + 
  geom_text(data = numbers, 
            size = 5,
            mapping = aes(x = -1000, y = 0.3, 
                          label = paste0(round(prop_n * 100, 0), "%"))) 

p

## Save figure
# ggsave("facet.png", plot = p)







### Figure 4: How policy affects the pdf


dist <- dist_all[dist_all$gdp == "uniform" & 
                   dist_all$haz == "uniform" & 
                   dist_all$cap %in% c("60", "120"), ]


cap_lab <- c("0%", "30%", 
             "60% capital increase", "90%", 
             "120% capital increase", "150%")
names(cap_lab) <- c("0", "30", "60", "90", "120", "150")


numbers <- dist%>% 
  group_by(haz, gdp, cap) %>% 
  summarize(prop_p = sum(val >= 0) / n_each,
            prop_n = 1 - sum(val >= 0) / n_each)
numbers$haz <- factor(numbers$haz,
                      levels = c("low", "uniform", "high"))
numbers$gdp <- factor(numbers$gdp,
                      levels = c("low", "uniform", "high"))
numbers$cap <- factor(numbers$cap,
                      levels = c("0", "30", "60", "90", "120", "150"))


g <- ggplot() +
  geom_histogram(data = subset(dist, val < 0),
                 mapping = aes(x = val, y = stat(count / n_each), 
                               color = "negative"),
                 fill = "red", col="grey", boundary = 0) +
  geom_histogram(data = subset(dist, val >= 0),
                 mapping = aes(x = val, y = stat(count / n_each), 
                               color = "positive"),
                 fill = "blue", col="grey", boundary = 0) +
  geom_vline(xintercept=0, linetype = "longdash") +
  facet_wrap(vars(cap), ncol = 1,
             labeller = labeller(cap = cap_lab)) +
  xlim(-1500, 1500) +
  xlab("Capital Level after second year (Million EUR)") + 
  ylab("Probability") +
  theme(plot.title = element_text(hjust=0.5),
        panel.spacing = unit(30, "pt"),
        strip.text.x = element_text(size = 15)) +
  geom_text(data = numbers, 
            size = 6,
            mapping = aes(x = 500, y = 0.12, 
                          label = paste0(round(prop_p * 100, 0), "%"))) + 
  geom_text(data = numbers, 
            size = 6,
            mapping = aes(x = -500, y = 0.12, 
                          label = paste0(round(prop_n * 100, 0), "%"))) +
  scale_color_manual(name = "Capital level", 
                     values = c("positive" = "blue", "negative" = "red")) +
  theme(legend.position = "bottom")



g


## Save figure
# 
# png("plot.png", width = 90, height = 60, units = 'cm', res = 500)
# plot(g)
# par(new = TRUE)
# legend("bottomleft", legend =c('specified (conditioned upon)', 'variable',
#                                'impact of concern'), 
#        pch=16, pt.cex=10, cex=5, bty='n',
#        col = c('steelblue3', 'lightblue', 'indianred2'))
# # mtext("Node colors", at=0.2, cex=2)
# # pdf("plot.pdf", width = 120, height = 80)
# dev.off()






## Figure 6: Bar chart of policy recommendations

### Probability distributions without axes

# 1
data <- tibble("Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
               "Probability" = c(0, 0.1, 0.2, 0.3, 0.4)) 
ggplot(data = data, mapping = aes(x = `Hazard increase`, y = `Probability`)) + 
  geom_bar(stat = "identity") +
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank())

# 2
data <- tibble("Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
               "Probability" = c(0.2, 0.2, 0.2, 0.2, 0.2)) 
ggplot(data = data, mapping = aes(x = `Hazard increase`, y = `Probability`)) + 
  geom_bar(stat = "identity") +
  ylim(0, 0.4) +
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank())

# 3
data <- tibble("Hazard increase" = c("1%", "3%", "5%", "7%", "9%"), 
               "Probability" = c(0.4, 0.3, 0.2, 0.1, 0)) 
ggplot(data = data, mapping = aes(x = `Hazard increase`, y = `Probability`)) + 
  geom_bar(stat = "identity") +
  theme(axis.title = element_blank(),
        axis.ticks = element_blank(),
        axis.text = element_blank())


# Bar chart


# simulate for each capital level to determine sufficient level

### Choose input  priors ###

# GDP
gdp_incr <- c(0.2, 0.2, 0.2, 0.2, 0.2) # neutral
# gdp_incr <- c(0, 0.1, 0.2, 0.3, 0.4) # optimistic
# gdp_incr <- c(0.4, 0.3, 0.2, 0.1, 0) # pessimistic

# Hazard intensity
haz_incr <- c(0.2, 0.2, 0.2, 0.2, 0.2) # neutral
# haz_incr <- c(0, 0.1, 0.2, 0.3, 0.4) # pessimistic
# haz_incr <- c(0.4, 0.3, 0.2, 0.1, 0) # optimistic

bn_model$exp_incr <- 
  array(gdp_incr,
        dimnames = list(c("1", "6", "11", "16", "21")))

bn_model$haz_incr <- 
  array(haz_incr,
        dimnames = list(c("1", "3", "5", "7", "9")))

result <- data.frame(NA, NA, NA, NA, NA)
evidence <- vector(mode = "list")
for (cyc in c(cycl_year_1, cycl_year_2)) {
  evidence[[cyc]] <- "1"
}
cap_lab <- c("0", "30", "60", "90", "120", "150")

for (cap in cap_lab) {
  
  evidence["cap_incr"] <- cap
  
  val <- cpdist(bn_model, nodes = "capy2",
                evidence = evidence,
                method = "lw")[["capy2"]]
  
  prop <- sum(val >= 0) / length(val)
  
  if (prop >= 0.1 & is.na(result[1,1])) {
    result[1, 1] <- paste0(cap, "% increase needed")
  } 
  if (prop >= 0.33 & is.na(result[1,2])) {
    result[1, 2] <-  paste0(cap, "% increase needed")  
  } 
  if (prop >= 0.5 & is.na(result[1,3])) {
    result[1, 3] <-  paste0(cap, "% increase needed")
  } 
  if (prop >= 0.67 & is.na(result[1,4])) {
    result[1, 4] <-  paste0(cap, "% increase needed")
  } 
  if (prop >= 0.9 & is.na(result[1,5])) {
    result[1, 5] <-  paste0(cap, "% increase needed") 
  }
  
}



data <- tibble(neutral = result[1, 3] |> parse_number(), 
               averse = result[1, 5] |> parse_number()) %>% 
  pivot_longer(cols = 1:2) %>% 
  mutate(name = factor(name, levels = c("neutral", "averse")))
g <- ggplot(data = data) + 
  geom_bar(mapping = aes(x = name, y = value), 
           fill = c("#00AFBB", "#FC4E07"),
           color = "grey",
           lwd = 1,
           stat = "identity", show.legend = F
  ) +
  geom_text(aes(x = name, y = value, label = paste0(value, "%")), 
            vjust = -0.3, size = 15) + 
  ylim(0, 170) +
  xlab("") + ylab("") +
  theme(axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank())


g


## Save figure
# png("bar6090.png", width = 12, height = 14, units = 'cm', res = 50)
# g
# dev.off()

