Plot_ag_trade_loss_ExpRes_dev.R

rm(list = ls())

# Libraries
library(tidyverse)
library(viridis)
library(patchwork)
library(hrbrthemes)
library(circlize)
library(chorddiag)  #devtools::install_github("mattflor/chorddiag")
library(reticulate)  # Load the reticulate package for working with Python objects
library(dplyr)
library(readr)

np <- import("numpy")

# Import label data
label <- read_csv("./country_name_merged.csv",
                    col_types = cols(
                      ISO = col_character(),
                      Country = col_character(),
                      region = col_character(),
                      subregion = col_character(),
                      continent = col_character(),
                      Country_Code = col_character(),
                      perIncome_percentiles = col_character()
                    ))

RR <- 189

years <- 2016:2050
TT <- length(years)


# Load no Export restriction order

PATH1 <- "./Results/TechProExtreme/temp/ag/"

bs_cty <- np$load(paste0(PATH1, 'order_0_cty.npy'))
pExtreme <- np$load(paste0(PATH1, 'pExtreme_126_order.npy'))
nExtreme <- np$load(paste0(PATH1, 'nExtreme_126_order.npy'))
cExtreme <- np$load(paste0(PATH1, 'cExtreme_126_order.npy'))

# Load with Export restriction order

expResSce <- c('Combined_5y', 'USA_5y', 'CHN_5y', 'IDN_5y', 'IND_5y')
expRes <- expResSce[1]

PATH2 <- "./Results/TechProExpResExtreme/temp/ag/"

pExtreme_res <- np$load(paste0(PATH2, expRes, '_pExtreme_126_order.npy'))
nExtreme_res <- np$load(paste0(PATH2, expRes, '_nExtreme_126_order.npy'))
cExtreme_res <- np$load(paste0(PATH2, expRes, '_cExtreme_126_order.npy'))


# Total intermediate and final demand
combine_and_sum <- function(data) {
  aa = data[, 1:RR, ] + data[, (RR + 1):ncol(data), ]
  return(apply(aa, c(1, 2), sum))
}

# Intermediate order
inter <- function(data) {
  return(data[, 1:RR, ])
}

# calculate loss of flow
calculate_dev <- function(data_res, data) {
  return(data_res - data)
}

calculate_dev_rate <- function(data_res, data, base){
  loss_dev <- data_res - data
  loss_rate_dev <- loss_dev / base * 100
  return (loss_rate_dev)
}

base_tot <- (bs_cty[, 1:RR] + bs_cty[, (RR+1):378]) * TT
pExtreme_tot <- combine_and_sum(pExtreme)
nExtreme_tot <- combine_and_sum(nExtreme)
cExtreme_tot <- combine_and_sum(cExtreme)
pExtreme_res_tot <- combine_and_sum(pExtreme_res)
nExtreme_res_tot <-  combine_and_sum(nExtreme_res)
cExtreme_res_tot <- combine_and_sum(cExtreme_res)

# calculate the deviation between with-restriction and non-restriction
dev_pExtreme <- calculate_dev(pExtreme_res_tot, pExtreme_tot) / 1e6 # convert thousand to billion
dev_nExtreme <- calculate_dev(nExtreme_res_tot, nExtreme_tot) / 1e6
dev_cExtreme <- calculate_dev(cExtreme_res_tot, cExtreme_tot) / 1e6


dev_pExtreme <- as.data.frame(-dev_pExtreme)
dev_nExtreme <- as.data.frame(-dev_nExtreme)
dev_cExtreme <- as.data.frame(-dev_cExtreme)

rownames(dev_pExtreme) <- colnames(dev_pExtreme) <- label$Country_Code
rownames(dev_nExtreme) <- colnames(dev_nExtreme) <- label$Country_Code
rownames(dev_cExtreme) <- colnames(dev_cExtreme) <- label$Country_Code


# 找到 'Dominican Republic' 在rownames和colnames中的索引
# row_index <- which(rownames(dev_pExtreme) == 'Dominican Republic')
# col_index <- which(colnames(dev_pExtreme) == 'Dominican Republic')
# rownames(dev_pExtreme)[row_index] <- colnames(dev_pExtreme)[col_index] <- 'Dominican'

# row_index <- which(rownames(dev_nExtreme) == 'Dominican Republic')
# col_index <- which(colnames(dev_nExtreme) == 'Dominican Republic')
# rownames(dev_nExtreme)[row_index] <- colnames(dev_nExtreme)[col_index] <- 'Dominican'

#row_index <- which(rownames(dev_cExtreme) == 'Dominican Republic')
#col_index <- which(colnames(dev_cExtreme) == 'Dominican Republic')
#rownames(dev_cExtreme)[row_index] <- colnames(dev_cExtreme)[col_index] <- 'Dominican'

# Calculate the sum for the row named "United States"
sum_dev_pExtreme_US <- sum(dev_pExtreme["USA", ])
sum_dev_nExtreme_US <- sum(dev_nExtreme["USA", ])
sum_dev_cExtreme_US <- sum(dev_cExtreme["USA", ])

# Print the results
print(paste("Sum of deviations in positive extreme for United States:", sum_dev_pExtreme_US))
print(paste("Sum of deviations in negative extreme for United States:", sum_dev_nExtreme_US))
print(paste("Sum of deviations in combined extreme for United States:", sum_dev_cExtreme_US))

# define data
data_long <- dev_cExtreme %>%
  rownames_to_column() %>%
  gather(key = 'key', value = 'value', -rowname) %>%
  arrange(desc(value))

total_value <- sum(data_long$value)
data_long <- data_long %>%
  mutate(cumulative_sum = cumsum(value))

# Find the rows that constitute the first 80% of the cumulative distribution
threshold_value <- total_value * 0.90
data_long_select <- data_long %>%
  filter(cumulative_sum <= threshold_value) %>%
  select(rowname, key, value)

# data_long_select <- data_long %>%
#   arrange(desc(value)) %>%
#   slice(1:30)

# order of country suffer most inflow and outflow loss
sum_by_key <- aggregate(value ~ key, data = data_long_select, sum, na.rm = TRUE)
sum_by_rowname <- aggregate(value ~ rowname, data = data_long_select, sum, na.rm = TRUE)
merged_data <- merge(sum_by_key, sum_by_rowname, by.x = "key", by.y = "rowname", all = TRUE) # value.x: inflow; value_y: outflow
merged_data[is.na(merged_data)] <- 0
final_sum <- aggregate(cbind(value.x, value.y) ~ key, data = merged_data, sum)
final_sum$total_sum <- final_sum$value.x + final_sum$value.y
final_sum <- final_sum[order(-final_sum$total_sum), ]
#print(final_sum[, c("key", "total_sum")])


# color palette
num_colors <- nrow(final_sum)
mycolor <- colorRampPalette(c("#B4DE2CFF", "#FDE725FF"))(num_colors)

country_colors <- data.frame(
  country = unique(final_sum$key),
  color = mycolor[1:length(unique(final_sum$key))]
)

# 创建特定国家与颜色的映射
specific_color_map <- c(
  "CHN" = "#FF0000",
  "USA" = "#B22234",
  'IDN' = '#F3D449',
  'IND' = '#E88D46'

)

# 创建特定地区与颜色的映射
region_color_map <- c(
  "Asia" = "#8B0000",
  "Europe" = "#00008B",
  'Americas' = '#FF8C00',
  'Oceania' = '#006400',
  'Northern Africa' = '#B8860B',
  'Sub-Saharan Africa' = '#A9A9A9'
  
)

# Function to get color based on region
get_region_color <- function(country, exp_res_countries) {
  region <- label$region[label$Country == country][1]  # Selecting the first match
  subregion <- label$subregion[label$Country == country][1]  # Selecting the first match
  if (country %in% exp_res_countries) {
    return("#000000")
  } else if (region %in% names(region_color_map)) {
    return(region_color_map[region])
  } else if (subregion %in% names(region_color_map)) {
    return(region_color_map[subregion])
  } else {
    return("#000000") # Default to black if no match
  }
}


# Apply country-specific colors to data
country_colors$color[country_colors$country %in% names(specific_color_map)] <- specific_color_map[country_colors$country[country_colors$country %in% names(specific_color_map)]]

order_cty <- final_sum$key
exp_res_countries <- order_cty[order_cty %in% c("CHN", "USA", "IND", "IDN")]
no_res_countries <- setdiff(order_cty, exp_res_countries)


# Plot circos

circos.clear()
circos.par(start.degree = 180, # 180 or 270
           gap.degree = 2,
           gap.after = 5,
           track.margin = c(-0.15, 0.15),
           points.overflow.warning = FALSE,
           canvas.xlim = c(-1.15, 1.15),  # horizontal, default is c(-1,1)
           canvas.ylim = c(-1.15, 1.15)   # vertical, , default is c(-1,1)
           )
par(mar = c(0, 0, 0, 0))

PATH2 <- './FigExpRes_Extreme/ag_plot/'
png(file=paste0(PATH2, expRes, '_trade_loss_dev_cExtreme_126.png'), height=2200, width=2300, bg = "transparent",res = 600)


# Base plot
chordDiagram(
  x = data_long_select, 
  grid.col = setNames(country_colors$color, country_colors$country),
  order = c(exp_res_countries, no_res_countries),
  transparency = 0.25,
  directional = 1,
  direction.type = c("arrows", "diffHeight"), 
  diffHeight  = -0.04,
  link.arr.type = "big.arrow",
  annotationTrack = "grid",
  annotationTrackHeight = c(0.05, 0.1)

)
# Add text and axis
circos.trackPlotRegion(
  track.index = 1, 
  bg.border = NA, 
  panel.fun = function(x, y) {
    
    xlim = get.cell.meta.data("xlim")
    sector.index = get.cell.meta.data("sector.index")
    
    for (i in seq_along(sector.index)) {
      label = sector.index[i]
      text_color <- get_region_color(label, exp_res_countries)
      
      circos.text(
        x = mean(xlim), 
        y = 3,  # Adjusted y-coordinate to add spacing between label and plot
        labels = label, 
        facing = "clockwise",
        niceFacing = TRUE,
        adj = c(0, 0.5),
        cex = 0.8,
        col = text_color
      )
    }
    
    # Add graduation on axis, conditionally hide labels if xlim[2] < 25 or 100
    if (xlim[2] >= 10) {
      circos.axis(
        h = "top", 
        major.at = seq(from = 0, to = xlim[2], by = 10), 
        minor.ticks = 0, 
        major.tick.length = 0.1,
        labels.niceFacing = TRUE,
        labels.cex = 0.8
      )
    } 
    par(cex.axis = 0.2)
  }
)

dev.off()