############################### ####################################
## Compares Qbf estimated values with hydrologic and process-based 
## as well as hydrologic-based Qbf predictors.
## Creates log-linear models and calculates goodness of fit metrics.
## Also creates most plots used in Chapter 3 
####################################################################

# Clear data
rm(list = ls())

# set working dir
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")
#library(wq) # Sens Slope
library(MASS) # RLM

source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")

# Import tables from various analyses. Should all have same # of rows
yield.metrics <- read.table("results_yield_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
flow.metrics <- read.table("results_flow_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
site.info <- read.table("site_summary.csv", stringsAsFactors = F, header = TRUE, sep = ",")
reg.results <- read.table("results_bed_regression.csv", stringsAsFactors = F, header = TRUE, sep = ",")
peaks <- read.table("results_peaks.csv", stringsAsFactors = F, header = TRUE, sep = ",")

# 
# yield.metrics$Qeff.spread <- round ((yield.metrics$Qeff.75 - yield.metrics$Qeff.25) 
#                                        / yield.metrics$Qeff.e,3)

# # Changes name so that tables can be merged based on "site.no"
# colnames(reg.results)[1:2] <- c("site.no", "site.nm")

# Check to see if tables match up: Should all be TRUE
sum( c(
  sum(yield.metrics$site.no == flow.metrics$site.no) == nrow(yield.metrics),
  sum(yield.metrics$site.no == site.info$site.no) == nrow(yield.metrics),
  sum(yield.metrics$site.no == reg.results$site.no) == nrow(yield.metrics), 
  sum(yield.metrics$site.no == peaks$site.no) == nrow(yield.metrics) ) ) == 4

# Calc some new, normalized flow metrics
# NOTE: site.info$Q1.5 contains Q1.5 calc'd from inst. peak annual max series
# with some Q1.5's calc'd from daily avg flow data where peak flow data not available
# Qbf <- ifelse(!is.na(site.info$Qbf), round(site.info$Qbf, 3), round(site.info$Q1.5, 3) )
Q1.5 <- site.info$Q1.5
Q2 <- peaks$Q2 * 0.028316847
Q2[which(is.na(Q2))] <- flow.metrics$Q2[which(is.na(Q2))]
Qavg <- peaks$Qavg * 0.028316847
Qbf <- site.info$Qbf
Qbf.mean <- round( Qbf / flow.metrics$mean, 3)
Qeff.mean <- round(yield.metrics$Qeff.e / flow.metrics$mean, 3)
Qeff.Q1.5 <- round(yield.metrics$Qeff.e / Q1.5, 3)
Q1.5.mean <- round(Q1.5 / flow.metrics$mean, 3)
Q1.5.Qavg <- round(Q1.5 / Qavg, 3)
Qeff.Qbf <-  round(yield.metrics$Qeff.e / Qbf, 3)
da.km2 <- round(site.info$da.km2, 0)

## Compare estimates of yield with Qbf
Q.25 <- yield.metrics$Q.25
Q.50 <- yield.metrics$Q.50
Q.75 <- yield.metrics$Q.75
Qeff <- yield.metrics$Qeff.e


# Remove pairs with really large differences (Percent error > 600%, squared error > 1E06)
Qbf.metrics <- data.frame(Qbf = Qbf, Qeff = Qeff, Q1.5 = Q1.5, Q2 = Q2, 
                          Q.25 = Q.25, Q.50 = Q.50, Q.75 = Q.75, 
                          sqr.err = (Q1.5 - Qbf)^2, pct.err = (Q1.5 - Qbf)/Qbf,
                          beta = yield.metrics$beta, cv = flow.metrics$cv, 
                          skewness = flow.metrics$skewness, RL = flow.metrics$yrs, 
                          DA = site.info$da.km2)

Qbf.metrics <- subset(Qbf.metrics, pct.err < 4 | is.na(pct.err))
Qbf.metrics <- subset(Qbf.metrics, sqr.err < 1E5 | is.na(pct.err))

Q.25 <- Qbf.metrics$Q.25
Q.50 <- Qbf.metrics$Q.50
Q.75 <- Qbf.metrics$Q.75
Qeff <- Qbf.metrics$Qeff
Qbf <- Qbf.metrics$Qbf
Q1.5 <- Qbf.metrics$Q1.5
Q2 <- Qbf.metrics$Q2

lm_Q25 <- lm(Qbf ~ Q.25)
a.25 <- exp(coefficients(lm_Q25)[1])
b.25 <- coefficients(lm_Q25)[2]
r2.25 <- round(summary(lm_Q25)$r.squared, 2)

lm_Q50 <- lm(Qbf ~ Q.50)
log.slope.Q50 <- round(summary(lm(log(Qbf) ~ log(Q.50)))$coef[[2]],2)
a.50 <- exp(coefficients(lm_Q50)[1])
b.50 <- coefficients(lm_Q50)[2]
r2.50 <- round(summary(lm_Q50)$r.squared, 2)
SSE.Q50 <- signif(sum((Q.50 - Qbf)^2, na.rm = T),2)
RMSE.Q50 <- signif(sqrt( sum((Q.50 - Qbf)^2, na.rm = T) / sum(!is.na(Q.50 - Qbf)) ), 2)
MAPE.Q50 <- signif(sum(abs(Q.50 - Qbf)/Qbf, na.rm = T) / sum(!is.na(Q.50 - Qbf)) * 100, 3)
PE.Q50 <- (Q.50 - Qbf)/Qbf * 100
MPE.Q50 <- mean((Q.50 - Qbf)/Qbf * 100, na.rm = T)
sum_Q50 <- summary(lm_Q50)
slope.Q50 <- round(sum_Q50$coef[[2]],2)
se_slp.Q50 <- round(sum_Q50$coef[[4]],3)

lm_Q75 <- lm(Qbf ~ Q.75)
a.75 <- exp(coefficients(lm_Q75)[1])
b.75 <- coefficients(lm_Q75)[2]
r2.75 <- round(summary(lm_Q75)$r.squared, 2)

lm_Qeff <- lm(Qbf ~ Qeff)
log.slope.Qeff <- round(summary(lm(log(Qbf) ~ log(Qeff)))$coef[[2]],2)
a.eff <- exp(coefficients(lm_Qeff)[1])
b.eff <- coefficients(lm_Qeff)[2]
r2.eff <- round(summary(lm_Qeff)$r.squared, 2)
SSE.Qeff <- signif(sum((Qeff - Qbf)^2, na.rm = T),2)
RMSE.Qeff <- signif(sqrt( sum((Qeff - Qbf)^2, na.rm = T) / sum(!is.na(Qeff - Qbf)) ), 2)
MAPE.Qeff <- signif(sum(abs(Qeff - Qbf)/Qbf, na.rm = T) / sum(!is.na(Qeff - Qbf)) * 100, 3)
PE.Qeff <- (Qeff - Qbf)/Qbf * 100
MPE.Qeff <- mean((Qeff - Qbf)/Qbf * 100, na.rm = T)
sum_Qeff <- summary(lm_Qeff)
slope.Qeff <- round(sum_Qeff$coef[[2]],2)
se_slp.Qeff <- round(sum_Qeff$coef[[4]],3)

lm_Q1.5 <- lm(Qbf ~ Q1.5)
r2.Q1.5 <- round(summary(lm_Q1.5)$r.squared, 2)
SSE.Q1.5 <- signif(sum((Q1.5 - Qbf)^2, na.rm = T), 2)
RMSE.Q1.5 <- signif(sqrt( sum((Q1.5 - Qbf)^2, na.rm = T) / sum(!is.na(Q1.5 - Qbf)) ), 2)
MAPE.Q1.5 <- signif(sum(abs(Q1.5 - Qbf)/Qbf, na.rm = T) / sum(!is.na(Q1.5 - Qbf)) * 100, 3)
PE.Q1.5 <- (Q1.5 - Qbf)/Qbf * 100
MPE.Q1.5 <- mean((Q1.5 - Qbf)/Qbf * 100, na.rm = T)
sum_Q1.5 <- summary(lm_Q1.5)
slope.Q1.5 <- round(sum_Q1.5$coef[[2]],2)
log.slope.Q1.5 <- round(summary(lm(log(Qbf) ~ log(Q1.5)))$coef[[2]],2)
se_slp.Q1.5 <- round(sum_Q1.5$coef[[4]],3)

lm_Q2 <- lm(Qbf ~ Q2)
r2.Q2 <- round(summary(lm_Q2)$r.squared, 2)
SSE.Q2 <- signif(sum((Q2 - Qbf)^2, na.rm = T),2)
RMSE.Q2 <- signif(sqrt( sum((Q2 - Qbf)^2, na.rm = T) / sum(!is.na(Q2 - Qbf)) ), 2)
MAPE.Q2 <- signif(sum(abs(Q2 - Qbf)/Qbf, na.rm = T) / sum(!is.na(Q2 - Qbf)) * 100, 3)
PE.Q2 <- (Q2 - Qbf)/Qbf * 100
MPE.Q2 <- mean((Q2 - Qbf)/Qbf * 100, na.rm = T)
sum_Q2 <- summary(lm_Q2)
slope.Q2 <- round(sum_Q2$coef[[2]],2)
log.slope.Q2 <- round(summary(lm(log(Qbf) ~ log(Q2)))$coef[[2]],2)
se_slp.Q2 <- round(sum_Q2$coef[[4]],3)

#### Characterize Pct Err as function of CV, and magnitude of Qbf, as well as beta

PE.table <- cbind(Qbf.metrics, PE.Qeff = PE.Qeff, PE.Qs50 = PE.Q50, 
                  PE.Q1.5 = PE.Q1.5, PE.Q2 = PE.Q2)

#setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/Figures")
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Writing/Qbf Paper/JHE Submission/Figs")
#pdf("PE_beta_skew_DA.pdf", width = 8, height = 3)
postscript("PE_beta_skew_DA.eps", width = 6.5, height = 3, paper = "special")
  #par(pty = "s", mar = c(2, 0, 0, 0), mfrow = c(1,3), oma = c(0, 4.5, .5, .5))
  par(pty = "s", mar = c(4.5, 0, 0, 0), mfrow = c(1,3), oma = c(2, 4.5, 0.5, 0.5), cex.lab = 1.2)
  
  ### Beta ###
  plot(PE.table$beta, PE.table$PE.Q1.5, pch = 0, lwd = 1.5, col = "royalblue",
       xlab = "",#bquote(beta),
       ylab = "Percent Error",
       cex.axis = 1.25,
       xlim = c(0,5))
  # axis(2, at = c(-60, -40, -20, 0, 20, 40, 60, 80, 100), 
  #      labels = c(-60, -40, -20, 0, 20, 40, 60, 80, 100),
  #      cex.axis  = 0.8,
  #      las = 2)
  abline(0,0, col = "gray40", lwd = 1.5)
  abline(MPE.Qeff, 0, col = "firebrick", lty = 6, lwd = 1.5)
  abline(MPE.Q50, 0, lty = 2, lwd = 1.5)
  abline(MPE.Q1.5, 0, col = "royalblue", lty = 4, lwd = 1.5)
  points(PE.table$beta, PE.table$PE.Qeff, pch = 2, lwd = 1.5, col = "firebrick")
  points(PE.table$beta, PE.table$PE.Qs50, pch = 20, lwd = 1.5)
  mtext("Percent Error", side = 2, line = 3, cex = 1.25)
    
  ### Skewness ###
#   skewness <- PE.table$skewness[which(PE.table$skewness < 7)]
#   PE.Qs50 <- PE.table$PE.Qs50[which(PE.table$skewness < 7)]
#   PE.Q1.5 <- PE.table$PE.Q1.5[which(PE.table$skewness < 7)]
#   PE.Qeff <- PE.table$PE.Qeff[which(PE.table$skewness > 7)]
  
  plot(PE.table$skewness, PE.table$PE.Q1.5, pch = 0, 
       xlim = c(0.5, 6),
       xlab = "",#"Skewness",
       yaxt = "n",
       ylab = "",
       cex.axis = 1.25,
       col = "royalblue")
  abline(0,0, col = "gray40", lwd = 1.5)
  points(PE.table$skewness, PE.table$PE.Qeff, pch = 2, lwd = 1.5, col = "firebrick")
  points(PE.table$skewness, PE.table$PE.Qs50, pch = 20, lwd = 1.5)
  abline(MPE.Qeff, 0, col = "firebrick", lty = 6, lwd = 1.5)
  abline(MPE.Q50, 0, lty = 2, lwd = 1.5)
  abline(MPE.Q1.5, 0, col = "royalblue", lty = 4, lwd = 1.5)

  plot(log10(PE.table$DA), PE.table$PE.Q1.5, pch = 0, 
       xlim = log10(c(10, 10000)),
       xaxt = "n",
       xlab = "",#"Skewness",
       yaxt = "n",
       ylab = "",
       cex.axis = 1.25,
       col = "royalblue")
  abline(0,0, col = "gray40", lwd = 1.5)
  points(log10(PE.table$DA), PE.table$PE.Qeff, pch = 2, lwd = 1.5, col = "firebrick")
  points(log10(PE.table$DA), PE.table$PE.Qs50, pch = 20, lwd = 1.5)
  abline(MPE.Qeff, 0, col = "firebrick", lty = 6, lwd = 1.5)
  abline(MPE.Q50, 0, lty = 2, lwd = 1.5)
  abline(MPE.Q1.5, 0, col = "royalblue", lty = 4, lwd = 1.5)
  minor.ticks.axis(1,9, mn = 1, mx = 4, cex = 1.25)
  text(2.5, 140, labels = "Coarse Bed Sites", font = 1, adj = c(0,0), cex = 1.25)

dev.off()