# ------------------------------------------------------------------
# Analyse anchor regression simulations for Fig. 06: 
# ------------------------------------------------------------------

# Sebastian Sippel
# 31.05.2021

library(raster)
library(ncdf4)
library(hydroGOF)
library(fields)
library(matrixStats)
library(bigmemory)
library(weights)



# ii. Load "maximum internal variability fraction" at "extremely likely" level:
load("../data/_final_processed/frac.internal.95percent.RData")


# 0. Function for plotting:
add.boxplot <- function(cur.data, at, outline = T, width=0.2, boxwex = 0.5, horizontal = T, whisker.quantiles = c(0.05, 0.95), border = "black", col = "white") {
  cur.boxplot.stats = boxplot(cur.data, outline=F, plot=F)
  cur.boxplot.stats$stats[1,] <- quantile(cur.data, whisker.quantiles[1], na.rm=T)
  cur.boxplot.stats$stats[5,] <- quantile(cur.data, whisker.quantiles[2], na.rm=T)
  if (cur.boxplot.stats$stats[1,] < 0) cur.boxplot.stats$stats[1,] = 0
  # bxp(cur.boxplot.stats, add = T, at = 5, frame=F, xaxt="n", yaxt = "n", )
  # bxp(cur.boxplot.stats, at = at, add=T, bty="n", boxwex=0.2, border = border, col = "black", frame=F, xaxt="n")
  bxp(cur.boxplot.stats, outline = outline, horizontal = horizontal, at = at, add=T, boxwex = boxwex, width=width, bty="n", notch = F, border = border, boxfill = col, col = col, frame=F, xaxt="n", yaxt = "n")
}

add.outlier <- function(x, at, whisker.quantiles = c(0.05, 0.95), pch = 16, col = "black", ...) {
  x.q = quantile(x, whisker.quantiles, na.rm=T)
  x. = x[which(x < x.q[1])]
  if (any(x. < 0)) x.[which(x.<0)] = 0;
  points(x = x., y = rep(at, length(x.)), pch = pch, col = col, ...)
  x. = x[which(x > x.q[2])]
  points(x = x., y = rep(at, length(x.)), pch = pch, col = col, ...)
}

add.points <- function(x, at, pch = 16, col = "black", ...) {
  x. = x
  if (any(x. < 0)) x.[which(x.<0)] = 0;
  points(x = x., y = rep(at, length(x.)), pch = pch, col = col, ...)
}




## ------------------------------------------------------------------------------------------
## Fig. 06: IPCC: Which fraction of observed trend *cannot* be explained by internal variability
## ------------------------------------------------------------------------------------------  
trend.lengths = c(30, 40, 50, 60); trend.names = c("y30", "y40", "y50", "y60")

## Plot ideas for signficance (no scaling!!):
## (1) Proportion of warming that could account for 2 or 3 sigma event (per model).
# 'Virtually certain': 99–100% probability, Very likely 90–100%, Likely 66–100
# 'Extremely likely': 95–100%
# 'Very likely': >90% chance


## FIG. 06:
# ----------------------------------------------------
set1 = brewer.pal(n = 8, name = "Set1")[c(1,4,5,2)]

## b) FRACTION THAT CAN BE EXPLAINED by internal variability  ('extremely likely' - 95% threshold):  
for(t.ix in 1:4) {
  
  cairo_pdf(filename = paste("_fig06/06_IPCC_", trend.lengths[t.ix], "_95perc.pdf", sep=""), width = 6, height = 5)
  {
    par(mar = c(5,1.5,3,1), mfrow = c(1,1))
    plot(c(1,1), type='n', xlim = c(0, 1), ylim = c(0.5, 5.5), main = "",
         # xlab = paste("Minimum external fraction of observed ", trend.lengths[t.ix],"-year trend at 'extremely likely' level [%]", sep=""), 
         xlab = paste("Fraction of attributable warming of observed ", trend.lengths[t.ix],"-year trend [%]", sep=""), 
         yaxs="i", ylab = "", yaxt="n", xaxt = "n", bty = "n", cex.lab = 1)
    mtext(text = "Detection metric", side = 2, line = 0, cex = 1)
    
    # Fraction of observed trend that internal variability could not account for at 'virtually certain' level [%]
    axis(side = 1, at = seq(0, 1, 0.1), labels = seq(0, 100, 10))
    axis(side = 1, las=2, at = seq(0, 1, 0.05), labels = F, tcl = 0.2)
    
    text(x = 0.02, y = 5.3, "Global mean temperature", pos = 4, col = "grey25", cex = 0.8)
    text(x = 0.02, y = 4.3, "Mean warming pattern", pos = 4, col = "grey25", cex = 0.8)
    text(x = 0.02, y = 3.3, "Ridge", pos = 4, col = "grey25", cex = 0.8)
    text(x = 0.02, y = 2.3, "Anchor (0.5/0.5)", pos = 4, col = "grey25", cex = 0.8)
    text(x = 0.02, y = 1.3, "Anchor (0.1/0.9)", pos = 4, col = "grey25", cex = 0.8)
    add.boxplot(cur.data = 1 - frac.internal2.AGMT[[2]][[t.ix]], at = 5, outline = F, horizontal = T, boxwex = 0.4, whisker.quantiles = c(0.05, 0.95), col = "darkgrey")
    add.outlier(x = 1 - frac.internal2.AGMT[[2]][[t.ix]], at = 5, pch = "|", col = "grey25", cex = 1.2)
    add.boxplot(cur.data = 1 - frac.internal2.EOF[[2]][[t.ix]], at = 4, outline = F, horizontal = T, boxwex = 0.4, whisker.quantiles = c(0.05, 0.95), col = "darkgrey")
    add.outlier(x = 1 - frac.internal2.EOF[[2]][[t.ix]], at = 4, pch = "|", col = "grey25", cex = 1.2)
    add.boxplot(cur.data = 1 - frac.internal2.ridge[[2]][[t.ix]], at = 3, outline = F, horizontal = T, boxwex = 0.4, whisker.quantiles = c(0.05, 0.95), col = "darkgrey")
    add.outlier(x = 1 - frac.internal2.ridge[[2]][[t.ix]], at = 3, pch = "|", col = "grey25", cex = 1.2)
    add.boxplot(cur.data = 1 - frac.internal2.opti_0.5[[2]][[t.ix]], at = 2, outline = F, horizontal = T, boxwex = 0.4, whisker.quantiles = c(0.05, 0.95), col = "darkgrey")
    add.outlier(x = 1 - frac.internal2.opti_0.5[[2]][[t.ix]], at = 2, pch = "|", col = "grey25", cex = 1.2)
    add.boxplot(cur.data = 1 - frac.internal2.opti_0.1[[2]][[t.ix]], at = 1, outline = F, horizontal = T, boxwex = 0.4, whisker.quantiles = c(0.05, 0.95), col = "darkgrey")
    add.outlier(x = 1 - frac.internal2.opti_0.1[[2]][[t.ix]], at = 1, pch = "|", col = "grey25", cex = 1.2)
  }
  dev.off()
  
  ### EOF10 scaling:
  col.bp = c(make.transparent.color("red", 120), "darkgrey", make.transparent.color("lightblue", alpha = 120))
  col.outlier = c("red", "black", "blue")
  
  cairo_pdf(filename = paste("_fig06/06_IPCC_", trend.lengths[t.ix], "-inv_EOF10_95perc.pdf", sep=""), width = 6, height = 6)
  {
    par(mar = c(5,1.5,1,1))
    plot(c(1,1), type='n', xlim = c(0, 1), ylim = c(0.5, 7), 
         # xlab = paste("Minimum external fraction of observed ", trend.lengths[t.ix],"-year trend at 'extremely likely' level [%]", sep=""), 
         xlab = paste("Fraction of attributable warming of observed ", trend.lengths[t.ix],"-year trend [%]", sep=""), 
         ylab = "",
         yaxs="i", yaxt="n", xaxt = "n", bty = "n", cex.lab = 1, main = "", cex = 1) 
    # Scaling of decadal-scale variability modes
    mtext(text = "Detection metric", side = 2, line = 0, cex = 1)
    
    # "scaling of decadal-scale variability modes (EOF10)"
    # Fraction of observed trend that internal variability could not account for at 'virtually certain' level [%]
    axis(side = 1, at = seq(0, 1, 0.1), labels = seq(0, 100, 10))
    axis(side = 1, las=2, at = seq(0, 1, 0.05), labels = F, tcl = 0.2)
    # axis(side = 2, las=2, at = 4 + c(-0.2, 0, +0.2), labels = c("x 1/2", "orig.", "x 2"), tcl = -0.2, cex.axis=0.7)
    # axis(side = 2, las=2, at = 3 + c(-0.2, 0, +0.2), labels = c("x 1/2", "orig.", "x 2"), tcl = -0.2, cex.axis=0.7)
    # axis(side = 2, las=2, at = 2 + c(-0.2, 0, +0.2), labels = c("x 1/2", "orig.", "x 2"), tcl = -0.2, cex.axis=0.7)
    # axis(side = 2, las=2, at = 1 + c(-0.2, 0, +0.2), labels = c("x 1/2", "orig.", "x 2"), tcl = -0.2, cex.axis=0.7)
    legend("top", c("Doubling of decadal varibility modes", "Original", "Halving of decadal variability modes"), lwd = 8, lty = 1, cex = 0.6, col = col.bp, inset = 0.02,
           title = "Scaling of internal decadal variability modes")
    
    text(x = 0.02, y = 5.4, "Global mean temperature", pos = 4, col = "grey25", cex = 0.8)
    sapply(1:3, FUN=function(i) {
      add.boxplot(cur.data = 1 - frac.internal2.AGMT_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 5 + c(0.2, 0, -0.2)[i], outline = F, horizontal = T, boxwex = 0.1, whisker.quantiles = c(0.05, 0.95), col = col.bp[i])
      add.outlier(x = 1 - frac.internal2.AGMT_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 5 + c(0.2, 0, -0.2)[i], pch = "|", col = col.outlier[i], cex = 0.6)
    })
    
    text(x = 0.02, y = 4.4, "Mean warming pattern", pos = 4, col = "grey25", cex = 0.8)
    sapply(1:3, FUN=function(i) {
      add.boxplot(cur.data = 1 - frac.internal2.EOF_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 4 + c(0.2, 0, -0.2)[i], outline = F, horizontal = T, boxwex = 0.1, whisker.quantiles = c(0.05, 0.95), col = col.bp[i])
      add.outlier(x = 1 - frac.internal2.EOF_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 4 + c(0.2, 0, -0.2)[i], pch = "|", col = col.outlier[i], cex = 0.6)
    })
    
    text(x = 0.02, y = 3.4, "Ridge", pos = 4, col = "grey25", cex = 0.8)
    sapply(1:3, FUN=function(i) {
      add.boxplot(cur.data = 1 - frac.internal2.ridge_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 3 + c(0.2, 0, -0.2)[i], outline = F, horizontal = T, boxwex = 0.1, whisker.quantiles = c(0.05, 0.95), col = col.bp[i])
      add.outlier(x = 1 - frac.internal2.ridge_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 3 + c(0.2, 0, -0.2)[i], pch = "|", col = col.outlier[i], cex = 0.6)
    })
    
    text(x = 0.02, y = 2.4, "Anchor (0.5/0.5)", pos = 4, col = "grey25", cex = 0.8)
    sapply(1:3, FUN=function(i) {
      add.boxplot(cur.data = 1 - frac.internal2.opti_0.5_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 2 + c(0.2, 0, -0.2)[i], outline = F, horizontal = T, boxwex = 0.1, whisker.quantiles = c(0.05, 0.95), col = col.bp[i])
      add.outlier(x = 1 - frac.internal2.opti_0.5_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 2 + c(0.2, 0, -0.2)[i], pch = "|", col = col.outlier[i], cex = 0.6)
    })
    
    text(x = 0.02, y = 1.4, "Anchor (0.1/0.9)", pos = 4, col = "grey25", cex = 0.8)
    sapply(1:3, FUN=function(i) {
      add.boxplot(cur.data = 1 - frac.internal2.opti_0.1_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 1 + c(0.2, 0, -0.2)[i], outline = F, horizontal = T, boxwex = 0.1, whisker.quantiles = c(0.05, 0.95), col = col.bp[i])
      add.outlier(x = 1 - frac.internal2.opti_0.1_EOF10[[c(4, 2, 1)[i]]][[t.ix]], at = 1 + c(0.2, 0, -0.2)[i], pch = "|", col = col.outlier[i], cex = 0.6)
    })
  }
  dev.off()
}