Investigating short duration Areal Reduction Factors https://tonyladson.wordpress.com/2020/04/19/arr2019-areal-reduction-factors-wobbles-in-short-duration-arfs/

ARF_sd.R

library(tidyverse)
library(optimx)

devtools::source_url("https://gist.githubusercontent.com/TonyLadson/fc870cf7ebfe39ea3d1a812bcc53c8fb/raw/d8112631a92a32be749cabe334a22931c035711e/ARF2019.R?raw=TRUE")


#source(file.path('ARR2019_ARF', "ARF_2019.R"))
#source(file.path('ARR2019_ARF', "ARF_tests.R")) # Check that we pass tests


# Functions and constants ------------------------------------------------------



log_breaks = function(maj, radix=10) {
  function(x) {
    minx         = floor(min(logb(x,radix), na.rm=T)) - 1
    maxx         = ceiling(max(logb(x,radix), na.rm=T)) + 1
    n_major      = maxx - minx + 1
    major_breaks = seq(minx, maxx, by=1)
    if (maj) {
      breaks = major_breaks
    } else {
      steps = logb(1:(radix-1),radix)
      breaks = rep(steps, times=n_major) +
        rep(major_breaks, each=radix-1)
    }
    radix^breaks
  }
  
  
  
  
}
scale_x_log_eng = function(..., radix=10) {
  scale_x_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}
scale_y_log_eng = function(..., radix=10) {
  scale_y_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}


params <-
  structure(list(`East Coast North` =     c(0.327,  0.241, 0.448, 0.36,  0.00096,   0.48,  -0.21,   0.012,   -0.0013), 
                 `Semi-arid Inland QLD` = c(0.159,  0.283, 0.25,  0.308, 7.3e-07,   1,      0.039,  0,        0), 
                 Tasmania =               c(0.0605, 0.347, 0.2,   0.283, 0.00076,   0.347,  0.0877, 0.012,   -0.00033), 
                 `SW WA` =                c(0.183,  0.259, 0.271, 0.33,  3.845e-06, 0.41,   0.55,   0.00817, -0.00045), 
                 `Central NSW` =          c(0.265,  0.241, 0.505, 0.321, 0.00056,   0.414, -0.021,  0.015,   -0.00033), 
                 `SE Coast` =             c(0.06,   0.361, 0,     0.317, 8.11e-05,  0.651,  0,      0,        0), 
                 `Southern Semi-arid` =   c(0.254,  0.247, 0.403, 0.351, 0.0013,    0.302,  0.058,  0,        0), 
                 `Southern Temperate` =   c(0.158,  0.276, 0.372, 0.315, 0.000141,  0.41,   0.15,   0.01,    -0.0027), 
                 `Northern Coastal` =     c(0.326,  0.223, 0.442, 0.323, 0.0013,    0.58,  -0.374,  0.013,   -0.0015), 
                 `Inland Arid` =          c(0.297,  0.234, 0.449, 0.344, 0.00142,    0.216, 0.129,  0,        0)), 
            class = "data.frame", row.names = c("a", "b", "c", "d", "e", "f", "g", "h", "i"))




# Plots -------------------------------------------------------------------


## Plots against duration


AEP <- c(0.005)
area <- c(5, 10, 100, 1000, 10000, 20000)
duration <- 10^seq(0,4,length.out = 500)





xdf <- expand_grid(area = area, AEP = AEP, duration = duration)

p <- xdf %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Tasmania', params = params))  %>% 
  ggplot(aes(x = duration, y = ARF_calc, colour = factor(area))) +
  geom_line(size = 0.2) +
  
  geom_vline(xintercept = 720, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  geom_vline(xintercept = 1440, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  scale_x_log_eng(name = 'Duration (min)', labels = scales::label_comma(accuracy = 1)) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_hue(name = bquote('Area ('~km^2*')'))+
  labs(title = 'Tasmania',
       subtitle = "AEP = 0.5%") +
  theme_grey(base_size = 5) +
  theme(legend.position = c(0.6, 0.3)) + 
  theme(legend.key.size = unit(2, 'mm')) +
  annotate(geom = 'text',x = 400, y = 0.02, label = 'Short duration', hjust = 'right', size = 1.5) +
  annotate(geom = 'text',x = 2640, y = 0.02, label = 'Long duration', hjust = 'left', size = 1.5) +
  annotate(geom = "segment", x = 720, y = 0.02, xend = 420, yend = 0.02,
               arrow = arrow(angle = 20, length = unit(1, 'mm'), ends = "last", type = "closed" ), 
           size = 0.2) +
  annotate(geom = "segment", x = 1440, y = 0.02, xend = 2440 , yend = 0.02,
           arrow = arrow(angle = 20, length = unit(1, 'mm'), ends = "last", type = "closed" ),
           size = 0.2)
  
# These plots look odd on-screen but are fine when saved and loaded into a blog
# The warning messages relate to attempting to calculate short duration ARFs for areas
# larger than 1000 km-squared.

p

ggsave(file.path('figure','ARF-duration-smallAEP.png'), plot = p, width = 4, height = 3)

# Derivative of the Short Duration ARF equation
# First Derivative with respect to duration

DARF_short <- function(A, D, P){
  
  (0.10332 * (A^0.265 - 0.190655 * log(D)))/(D^1.36) + 
    (0.0002825 *A^0.226 * (log(P)/log(10) + 0.3))/D^0.875 - 
    9.46938e-7*10^(-0.0000145833 * (D - 180)^2) * A^0.213 * (D - 180) * (log(P)/log(10) + 0.3) + 
    0.0547181/D^1.36
}  


# Second derivative with respect to duration

D2ARF_short <- function(A, D, P){
  #Second derivative
  -(0.140515 * (A^0.265 - 0.190655 * log(D)))/D^2.36 - 
    (0.000247188 * A^0.226 * (log(P)/log(10) + 0.3))/D^1.875 + 
    0.0141 * A^0.213 * (4.5103e-9 * 10^(-0.0000145833 * (D - 180)^2) * (D - 180)^2 - 0.0000671587 * 10^(-0.0000145833 * (D - 180)^2)) * (log(P)/log(10) + 0.3) - 
    0.0941151/D^2.36
}



# Create plot that shows the short duration ARF and its derivative

duration <- 10^seq(0,4,length.out = 500)
area <- c(100)
aep <- 0.005




p1 <- expand_grid(duration, area, aep) %>% 
  rowwise() %>% 
  mutate(ARF_s = ARF_short(area = area, duration = duration, aep = aep) ) %>% 
  ggplot(aes(x = duration, y = ARF_s)) +
  geom_line(size = 0.3) +
  geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
  scale_x_continuous(name = 'Duration (min)', trans = 'log10', limits = c(NA, 720), breaks = c(1, 10, 100, 180, 720)) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  labs(title = 'Tasmania',
       subtitle=expression(paste("0.5% AEP, Area = 100 ", km^2))) +
  theme_grey(base_size = 5)



p1



xdf_ARFs <- expand_grid(duration, area, aep) %>% 
  rowwise() %>% 
  mutate(ARF_s = ARF_short(area = area, duration = duration, aep = aep))  


p2 <- expand_grid(duration, area, aep) %>% 
  rowwise() %>% 
  mutate(DARF_s = DARF_short(A = area, D = duration, P = aep)) %>% 
  ggplot(aes(duration, DARF_s)) +
  geom_line(colour = 'red', size = 0.2) +
  geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
  #geom_line(data = xdf_ARFs, aes(duration, ARF_s)) +
  scale_x_continuous(name = 'Duration (min)', trans = 'log10', limits = c(NA, 720), breaks = c(1, 10, 100, 180, 720)) +
  scale_y_continuous(name = 'Derivative(Areal Reduction Factor)', trans = 'log10') +
  theme_grey(base_size = 5)

p2

gridExtra::grid.arrange(p1, p2)

p <- gridExtra::arrangeGrob(p1, p2)
ggsave(file.path('figure','ARF-derivative.png'), plot = p, width = 4, height = 3)



# Breaking the short duration equation up into three parts

# Constants for short duration ARF

a <- 0.287	
b <- 0.265
c <- 0.439
d <- 0.36
e <- 0.00226	
f <- 0.226
g <- 0.125
h <- 0.0141
i <- -0.021
j <- 0.213


ARF_short <- function(area, duration, aep) {
  
  
  
  min(1, (1 - a*(area^b - c*log10(duration)) * duration^(-d) +
            e*area^f*duration^g * (0.3 + log10(aep)) + 
            h * area^j * 10^(i* (1/1440) * (duration - 180)^2)  * (0.3 + log10(aep))))
}

ARF_short1 <- function(area, duration, aep){
  
  -a*(area^b - c*log10(duration)) * duration^(-d)
}


ARF_short2 <- function(area, duration, aep){
  
  e*area^f*duration^g * (0.3 + log10(aep)) 
  
  
}  



ARF_short3 <- function(area, duration, aep){
  
  h * area^j * 10^(i* (1/1440) * (duration - 180)^2)  * (0.3 + log10(aep))
  
}


# Test the equations
ARF_short1(area = 1000, duration = 200, aep = 0.01)
ARF_short2(area = 1000, duration = 200, aep = 0.01)
ARF_short3(area = 1000, duration = 200, aep = 0.01)

# Plot the three parts

duration <- 10^seq(0,4,length.out = 500)

p <- expand_grid(area = 1000, duration = duration, aep = 0.005) %>% 
  rowwise() %>% 
  mutate(s1 = ARF_short1(area, duration, aep)) %>% 
  mutate(s2 = ARF_short2(area, duration, aep)) %>% 
  mutate(s3 = ARF_short3(area, duration, aep)) %>% 
  mutate(s4 = 1+s1+s2+s3) %>% 
  pivot_longer(cols = c(s1, s2, s3, s4), names_to = 'component') %>% 
  ggplot(aes(x = duration, y = value, colour = component)) +
  geom_line(size = 0.2) +
  geom_vline(xintercept = 180, linetype = 'dotted', size = 0.2) +
  scale_colour_hue(name = 'Eqn parts', labels = c('Part 1', 'Part 2', 'Part 3', 'Sum')) +
  scale_x_continuous(name = 'Duration (min)', 
                     trans = "log10", 
                     breaks = c(1, 10, 100, 180, 720), 
                     limits = c(NA, 720)) +
  #theme(legend.position = c(0.65, 0.2)) +
  theme_gray(base_size = 5) +
  theme(legend.key.size = unit(2, 'mm'), legend.position = c(0.65, 0.2)) 

p  

ggsave(file.path('figure','ARF-short_parts.png'), plot = p, width = 4, height = 3)


# test of the (D-180)^2 part
T1 <- 10^(-0.021* (1/1440)*(D-180)^2)
D <- 100:300
plot(T1)


## Plots against aep


# Make a sequence of AEP values evenly spaced in as quantiles of the normal distribution
z_seq <- seq(qnorm(0.5), qnorm(0.005), length.out = 200)
AEP <- pnorm(z_seq)
area <- c(5, 10, 100, 1000, 10000, 20000)
duration <- c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080)

aep.labs <-  c(0.5, 1, 2, 5, 10, 20, 50)
names(aep.labs) <- c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)
aep_breaks = c(0.005, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5)

duration.labs <-  str_c(c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080), ' minute')
names(duration.labs) <- c(1, 5, 10, 100, 180, 720, 1000, 1440, 10080)




xdf <- expand_grid(area = area, AEP = AEP, duration = duration)



p <- xdf %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Northern Coastal', params = params)) %>% 
  ggplot(aes(x = 1-AEP, y = ARF_calc, colour = factor(area))) +
  geom_line(size = 0.2) +
  scale_x_continuous(name = 'AEP(%)', trans = 'probit', breaks = 1-aep_breaks, labels = aep.labs) +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_hue(name = expression(paste("Area ", km^2))) +
  facet_wrap(~duration, labeller = labeller(duration = duration.labs)) +
  theme_grey(base_size = 5) +
  theme(legend.key.size = unit(1.5, 'mm'), legend.position = c(0.8,0.12)) 
  
p
ggsave(file.path('figure','ARF-AEP_facet.png'), plot = p, width = 4, height = 3)
  
# What combination of area and duration provides the greatest change in ARF as AEP goes from 50% to 0.05%
# Note that Short duration ARFs are not a function of region.


#https://stackoverflow.com/questions/44224761/getting-error-using-optim-function-with-bounds-in-r
# Uses penalties to enforce bounds to keep area between 0 and 1000
# and duration between 1 and 720

# The objective is the difference between the ARF at a low AEP and the ARF and a high AEP
# Find where this is a maximum

Obj <- function(par){
  area <- par[1]
  duration <- par[2]
  if(area > 1000) return(100 + sum(par^2))
  if(area < 0) return(100 + sum(par^2))
  if(duration > 720) return(100 + sum(par^2))
  if(duration < 1)return(100 + sum(par^2))
  -ARF_short(area, duration, aep = 0.5) +
    ARF_short(area, duration, aep = 0.005)
}

par.init <- (c(10, 10))
Obj(par.init)
optimx(par = par.init, fn = Obj, method = 'Nelder-Mead')

# I tried this using several different starting points and got similar results.

# > optimx(par = par.init, fn = Obj, method = 'Nelder-Mead')
#               p1       p2      value fevals gevals niter convcode  kkt1  kkt2 xtime
# Nelder-Mead 1000 183.4142 -0.1640775    161     NA    NA        0 FALSE FALSE 0.001


# i.e. 1000 km-squared catchment and a duration of 183.4 min