R/tune.R

In flr/mse: Tools for Running Management Strategy Evaluations using FLR

# tune.R - DESC
# /tune.R

# Copyright European Union, 2018
# Author: Iago Mosqueira (EC JRC) <iago.mosqueira@ec.europa.eu>
#
# Distributed under the terms of the European Union Public Licence (EUPL) V.1.1.

# tunebisect: BISECTION tuning {{{

#' 
#' @examples
#' # dataset contains both OM (FLom) and OEM (FLoem)
#' data(sol274)
#' # choose sa and hcr
#' control <- mpCtrl(list(
#'   est = mseCtrl(method=perfect.sa),
#'   hcr = mseCtrl(method=hockeystick.hcr, args=list(lim=0,
#'   trigger=41500, target=0.27))))
#' # load statistics
#' data(statistics)
#' tun <- tunebisect(om, oem=oem, control=control, args=list(iy=2021, fy=2035),
#' tune=list(target=c(0.15, 0.35)),
#' metrics=list(SB=ssb), statistic=statistics['PSBMSY'], years=2025:2034)
#' # Plot tuned MP
#' plot(om, tun)

tunebisect <- function(om, oem="missing", control, metrics, statistic, tune,
  args, prob=0.5, tol=0.01, maxit=12, verbose=TRUE,
  years=ac(seq(args$iy + 1, args$fy - 1)), ...) {
  
  # args
  args$fy <- if(is.null(args$fy)) dims(om)$maxyear else args$fy

  # CHECK years
  if(!all(unique(unlist(years)) %in% seq(args$iy, args$fy)))
    stop("Years for statistic computation 'years' outside of 'args' year range (iy:fy).")

  # CHECK that tune names match args(control$hcr)
  if(!names(tune) %in% formalArgs(method(control$hcr)))
    stop("Element in tune must be a argument of method(control$hcr).")

  # CHECK length(tune)  == 1
  if(length(tune) > 1)
    stop(paste0("'tunebisect' can only tune for a single HCR argument, got ",
      names(tune), "."))

  # CHECK length(tune[[1]]) == 2
  if(length(tune[[1]]) != 2)
    stop("Range of argument for tuning has more than 2 values, need min and max.")

  # CHECK statistic: single statistic as length 1 list
  if(length(statistic) != 1) {
    stop("'statistic' must be a named list of length 1 with formula, 'name', and 'desc'")
  }

  # CHECK 0 < prob < 1
  if(!(prob >=0 & prob <=1))
    stop("prob must be a value between 0 and 1.")

  # --- RUN at min

  cmin <- control
  cmin$hcr@args[names(tune)] <- lapply(tune, '[', 1)

  # PRINT at top
  if(verbose)
    message(paste0("[1] ", names(tune), ": ",
      unlist(cmin$hcr@args[names(tune)])))

  rmin <- mp(om, oem=oem, ctrl=cmin, args=args, scenario=paste0("min"),
    verbose=FALSE, ...)

  pmin <- performance(rmin, metrics=metrics, 
    statistics=statistic, refpts=refpts(om), probs=NULL, years=years)
  obmin <- mean(pmin$data, na.rm=TRUE) - prob

  # PRINT result
  if(verbose)
    message(paste0("[1] diff: ", format(obmin, digits=2)))

  # CHECK cmin result
  if(isTRUE(all.equal(obmin, 0, tolerance=tol)))
    return(rmin)
  
  # --- RUN at max

  cmax <- control
  cmax$hcr@args[names(tune)] <- lapply(tune, '[', 2)

  # PRINT at top
  if(verbose)
    message(paste0("[2] ", names(tune), ": ",
      unlist(cmax$hcr@args[names(tune)])))

  rmax <- mp(om, oem=oem, ctrl=cmax, args=args, scenario=paste0("max"),
    verbose=FALSE, ...)
  
  pmax <- performance(rmax, metrics=metrics,
    statistic=statistic, refpts=refpts(om), probs=NULL, years=years)
  obmax <- mean(pmax$data, na.rm=TRUE) - prob

  # PRINT result
  if(verbose)
    message(paste0("[2] diff: ", format(obmax, digits=2)))
  
  # CHECK cmax result
  if(isTRUE(all.equal(obmax, 0, tolerance=tol)))
    return(rmax)

  # CHECK range includes 0
  # TODO: PROPOSE wider range, min * 0.5, max * 1.5, and restart (recursive?).
  if((obmin * obmax) > 0) {
    warning("Range of hcr param(s) cannot achieve requested tuning objective probability")
    return(list(min=rmin, max=rmax))
  }

  # --- LOOP bisecting

  count <- 0
  while(count <= maxit) {

    # RUN at mid
    cmid <- control
    cmid$hcr@args[[names(tune)]] <-
      (cmin$hcr@args[[names(tune)]] + cmax$hcr@args[[names(tune)]]) / 2

    # PRINT at top
    if(verbose)
      message(paste0("[", count + 2, "] ", names(tune), ": ",
        format(unlist(cmid$hcr@args[names(tune)]), digits=3)))

    rmid <- mp(om, oem=oem, ctrl=cmid, args=args, scenario=paste0("mid"),
      verbose=FALSE, ...)

    pmid <- performance(rmid, metrics=metrics, 
      statistics=statistic, refpts=refpts(om), probs=NULL, years=years)
    obmid <- mean(pmid$data, na.rm=TRUE) - prob

    # PRINT result
    if(verbose)
      message(paste0("[", count + 2, "] diff: ", format(obmid, digits=2)))
  
    # CHECK and RETURN cmid result
    if(isTRUE(all.equal(obmid, 0, tolerance=tol))) {
      return(rmid)
    }

    # TEST LEFT
    if((obmin * obmid) < 0) {

      # SET max as new mid
      cmax <- cmid
      obmax <- obmid

    } else {
      
      # SET min as new mid
      cmin <- cmid
      obmin <- obmid
    }

    count <- count + 1
  }

  stop("Solution not found within 'maxit', check 'range', 'maxit' or 'tol'.")

} # }}}

# tunegrid: GRID tuning {{{

#' @examples
#' data(sol274)
#' # tes <- tunegrid(om, oem, control, metric=list(SB=ssb), statistic=stats['S8'],
#' #  grid=grid, args=mseargs, years=list(2030:2039), iters=1:100)

tunegrid <- function(om, oem="missing", control, metric, statistic, grid, args,
  years=ac(seq(args$iy + 1, args$fy)), parallel=FALSE, verbose=FALSE, 
  iters=dims(om)$iter, ...) {

  if(!missing(iters)) {
    om <- iter(om, iters)
    oem <- iter(oem, iters)
  }

  fixed <- as.data.frame(args(control$hcr)[!names(args(control$hcr)) %in% names(grid)])

  # dopar LOOP over grid

  res <- foreach(i=seq(dim(grid)[1]), .errorhandling = "remove") %dofuture% {

    # TODO ADD progressr

    cat("[", i, "]\n")

    # ASSIGN new control$hcr@args values
    args(control$hcr)[names(grid)] <- grid[i,]

    # RUN mp
    run <- mp(om=om, oem=oem, ctrl=control,
      args=args, parallel=parallel, verbose=verbose, ...)

    # COMPUTE performance & BIND to grid
    cbind(performance(run, metric=metric, statistic=statistic,
      years=years), grid[i,, drop=FALSE], fixed)
  }

  # TODO ADD empty row for missing or use try()and .combine

  res <- rbindlist(res, idcol="run")

  return(res)
}

# }}}

# bisect {{{

#' Bisection search for a forecast target providing a given performance statistic value.
#'
#' The plain bisection algorithm (Burden & Douglas, 1985) is employed here to
#' find the value of a given forecast target quantity (e.g. `fbar`) for which
#' a selected value of a performance statistic is obtained over a chosen period.
#' @references {Burden, Richard L.; Faires, J. Douglas (1985), "2.1 The Bisection Algorithm", Numerical Analysis (3rd ed.), PWS Publishers, ISBN 0-87150-857-5}
#' @param stock
#' @param sr
#' @param metrics
#' @param statistic
#' @param years
#' @param tune
#' @param prob
#' @param tol
#' @param maxit
#' @examples
#' data(ple4)
#' stock <- propagate(stf(ple4, end=2118), 100)
#' srr <- predictModel(model=rec~a*ssb*exp(-b*ssb), params=FLPar(a=5.20, b=1.65e-6))
#' # GENERATE SRR deviances
#' devs <- ar1rlnorm(rho=0.4, 2018:2118, iters=100, meanlog=0, sdlog=0.5)
#' # DEFINE Fp05 statistic
#' statistic <- list(FP05=list(~yearMeans((SB/SBlim) < 1), name="P.05",
#'   desc="ICES P.05"))
#' # CALL bisect over 100 years, Fp.05 calculated over last 50.
#' fp05fwd <- bisect(stock, sr=srr, deviances=devs, metrics=list(SB=ssb), 
#'   refpts=FLPar(SBlim=150000), statistic=statistic, years=2018:2118,
#'   pyears=2069:2118, tune=list(fbar=c(0.1, 1)), prob=0.05)

bisect <- function(stock, sr, deviances=rec(stock) %=% 1, metrics, refpts,
  statistic, years, pyears=years, tune, prob, tol=0.01, maxit=15, verbose=TRUE) {

  # CHOOSE method

  if(names(tune)[1] %in% c("f", "fbar"))
    foo <- ffwd
  else
    foo <- fwd

  # --- RUN at min

  cmin <- fwdControl(year=years, quant=names(tune)[1], value=unlist(tune)[1])

  # PRINT at top
  if(verbose)
    cat(paste0("[1] ", names(tune), ": ", unlist(tune)[1]))
  rmin <- foo(stock, sr=sr, control=cmin, deviances=deviances)
  
  pmin <- performance(rmin, metrics=metrics, 
    statistics=statistic, refpts=refpts, probs=NULL, years=pyears)

  obmin <- mean(pmin$data, na.rm=TRUE) - prob

  if(verbose)
    cat(" - prob:", mean(pmin$data, na.rm=TRUE), " - diff: ", obmin, "\n")
  
  # CHECK cmin result
  if(isTRUE(all.equal(obmin, 0, tolerance=tol)))
    return(rmin)
  
  # --- RUN at max

  cmax <- fwdControl(year=years, quant=names(tune)[1], value=unlist(tune)[2])

  # PRINT at top
  if(verbose)
    cat(paste0("[2] ", names(tune), ": ", unlist(tune)[2]))

  rmax <- foo(stock, sr=sr, control=cmax, deviances=deviances)
  
  pmax <- performance(rmax, metrics=metrics, 
    statistics=statistic, refpts=refpts, probs=NULL, years=pyears)
  obmax <- mean(pmax$data, na.rm=TRUE) - prob

  if(verbose)
    cat(" - prob:", mean(pmax$data, na.rm=TRUE), " - diff: ", obmax, "\n")
  
  # CHECK cmax result
  if(isTRUE(all.equal(obmax, 0, tolerance=tol)))
    return(rmax)
 
  # --- CHECK range includes 0
  if((obmin * obmax) > 0) {
    warning("Range of hcr param(s) cannot achieve requested tuning objective probability")
    return(list(min=rmin, max=rmax))
  } 

  # --- LOOP bisecting

  count <- 0
  while(count <= maxit) {

    # RUN at mid
    cmid <- control
    cmid <- fwdControl(year=years, quant=names(tune)[1],
      value=(cmin$value + cmax$value) / 2)

    # PRINT at mid
    if(verbose)
      cat(paste0("[", count + 3, "] ", names(tune), ": ", cmid$value[1]))

    rmid <- foo(stock, sr=sr, control=cmid, deviances=deviances)

    pmid <- performance(rmid, metrics=metrics, 
      statistics=statistic, refpts=refpts, probs=NULL, years=pyears)
    obmid <- mean(pmid$data, na.rm=TRUE) - prob

    if(verbose)
      cat(" - prob:", mean(pmid$data, na.rm=TRUE), " - diff: ", obmid, "\n")

    # CHECK and RETURN cmid result
    if(isTRUE(all.equal(obmid, 0, tolerance=tol))) {
      return(rmid)
    }

    # TEST LEFT
    if((obmin * obmid) < 0) {

      # SET max as new mid
      cmax <- cmid
      obmax <- obmid
      if(isTRUE(all.equal(cmin$value[1], cmid$value[1], tolerance=tol))) {
        return(rmid)
      }
    } else {

      # SET min as new mid
      cmin <- cmid
      obmin <- obmid
      if(isTRUE(all.equal(cmid$value[1], cmax$value[1], tolerance=tol))) {
        return(rmid)
      }
    }
    count <- count + 1
  }

  warning("Solution not found within 'maxit', check 'range', 'maxit' or 'tol'.")

  return(rmid)

} # }}}

# computeFp05 {{{

#' Calculates the Fbar value giving a maximum probability of ssb being below Blim of 5%
#'
#' @param stock An FLStock over which the calculation is carried out.
#' @param sr The stock-recruits relationship to use in fwd.
#' @param SBlim
#' @param range
#' @param nyears
#' @param sigmaR
#' @param rho
#' @param its Number of iterations
#' @param verbose Should progress be shown, TRUE.
#' @examples
#' data(ple4)
#' sr <- predictModel(model=bevholt, params=FLPar(a=1.4e6, b=1.5e5))
#' fp05 <- computeFp05(ple4, sr, SBlim=150000, its=300, range=c(0.40, 0.50))
#' # RUN projection for obtained Fp.05 value
#' proj <- fwd(propagate(stf(ple4, nyears=100), 300), sr=sr,
#'   fbar=FLQuant(fp05, dimnames=list(year=2018:2117)),
#'   deviances=ar1rlnorm(rho=0.43, years=2018:2117, iters=300, meanlog=0,
#'   sdlog=0.5))
#' plot(ssb(proj), prob=c(0.01, 0.25, 0.50, 0.75, 0.99)) +
#'   geom_hline(yintercept=150000)

computeFp05 <- function(stock, sr, SBlim, range=c(0.01, 0.75), nyears=3,
  sigmaR=0.5, rho=0.43, its=500, verbose=TRUE) {

  years <- seq(dims(stock)$maxyear + 1, length=100)
  pyears <- years[51:100]

  # GENERATE SRR deviances
  devs <- ar1rlnorm(rho=rho, years=years, iters=its, meanlog=0, sdlog=sigmaR)

  # EXTEND & PROPAGATE
  
  stock <- stf(stock, end=years[100], wts.nyears=nyears)
  stock <- propagate(stock, its)

  statistic <- list(FP05=list(~apply(iterMeans((SB/SBlim) < 1), c(1, 3:6), max),
    name="P.05", desc="ICES P.05"))

  res <- bisect(stock, sr=sr, refpts=FLPar(SBlim=SBlim), deviances=devs,
    metrics=list(SB=ssb), statistic=statistic, years=years, pyears=pyears, 
    tune=list(fbar=range), prob=0.05, tol=0.01, verbose=verbose)

  if(is(res, "list"))
    stop("FP.05 cannot be found within the range provided.")

  fp05 <- mean(fbar(res)[,100])

  return(c(Fp05=fp05))
}
# }}}