R/Movement_calcs.R

In MSEtool: Management Strategy Evaluation Toolkit

#' Calculates movement matrices from user inputs
#'
#' A wrapper function for \link{makemov} used to generate movement matrices for the operating model.
#' Calculates a movement matrix from user-specified unfished stock biomass fraction in each area and probability of
#' staying in the area in each time step.
#' @param OM Operating model, an object of class \linkS4class{OM}.
#' @param dist A vector of fractions of unfished stock in each area. The length of this vector will determine the
#' number of areas (\code{nareas}) in the OM.
#' @param prob Mean probability of staying across all areas (single value) or a vector of the probability of
#' individuals staying in each area (same length as dist).
#' @param distE Logit (normal) St.Dev error for sampling stock fractions from the fracs vector
#' @param probE Logit (normal) St.Dev error for sampling desired probability of staying either by area
#' (prob is same length as dist) or the mean probability of staying (prob is a single number).
#' @param prob2 Optional vector as long as prob and dist. Upper bounds on uniform sampling of
#' probability of staying, lower bound is prob.
#' @param figure Logical to indicate if the movement matrix will be plotted (mean values and range
#' across \code{OM@@nsim} simulations.)
#' @return The operating model \code{OM} with movement parameters in slot \code{cpars}.
#' The \code{mov} array is of dimension \code{nsim}, \code{maxage}, \code{nareas}, \code{nareas}.
#' @note Array \code{mov} is age-specific, but currently the movement generated by \code{simmov} is
#' independent of age.
#' @author T. Carruthers and Q. Huynh
#' @export
#' @examples
#' \dontrun{
#' movOM_5areas <- simmov(testOM, dist = c(0.01,0.1,0.2,0.3,0.39), prob = c(0.1,0.6,0.6,0.7,0.9))
#' movOM_5areas@@cpars$mov[1, 1, , ] # sim 1, age 1, movement from areas in column i to areas in row j
#' plot_mov(movOM_5areas@@cpars$mov)
#' plot_mov(movOM_5areas@@cpars$mov, type = "all")
#'
#' }
#' @describeIn simmov Estimation function for creating movement matrix.
simmov <- function(OM, dist = c(0.1, 0.2, 0.3, 0.4), prob = 0.5, distE = 0.1, probE = 0.1, prob2 = NA, figure = TRUE) {

  nareas <- length(dist)
  if(nareas < 2) stop("Error: nareas, i.e., length(dist), is less than 2.")
  nsim <- OM@nsim
  maxage <- OM@maxage

  dist_s <- rnorm(nareas*nsim, log(dist), distE) %>% matrix(nrow = nareas) %>% t() %>% ilogitm()

  if(length(prob) == 1) {
    prob_s <- rnorm(nsim, logit(prob), probE) %>% matrix(nrow = nsim) %>% ilogit()
  } else if(length(prob) == length(dist) && is.na(prob2)) {
    prob_s <- rnorm(nareas*nsim, logit(prob), probE) %>% matrix(nrow = nareas) %>% t() %>% ilogit()
  } else if(length(prob) == length(dist) && length(prob) == length(prob2)) {
    prob_s <- runif(nareas*nsim, prob, prob2) %>% matrix(nrow = nareas) %>% t()
  } else{
    stop("Error: either prob wasn't of length 1, or prob wasn't of length dist or prob 2 wasn't the same length as prob and dist.
            You have three options:
            (1) provide one value for prob which represents mean probability of staying across all areas sampled for each simulation with probE logit error
            (2) provide nareas values of prob which represent probability of staying across all areas sampled for each simulation with probE logit error
            (3) provide nareas values of prob and prob2 which are the upper and lower bounds for sampling uniform probability of staying for each area")
  }

  movt <- lapply(1:nsim, function(i) makemov(fracs = dist_s[i, ], prob = prob_s[i, ]) %>% array(c(nareas, nareas, maxage + 1)))
  mov <- simplify2array(movt) %>% aperm(c(4, 3, 1, 2))

  OM@cpars$mov <- mov
  if(figure) plot_mov(mov, type = "matrix")

  return(OM)
}

opt_mov <- function(x, fracs, prob, nareas) {
  grav_out <- grav(log_visc = x[1:length(prob)], log_grav = x[(length(prob)+1):length(x)],
                   fracs = fracs, nareas = nareas)

  nll_dist <- dnorm(log(grav_out$idist), log(fracs), 0.1, TRUE)
  if(length(prob) == 1) {
    nll_stay <- dnorm(log(grav_out$psum), log(prob), 0.1, TRUE)
  } else {
    nll_stay <- dnorm(log(diag(grav_out$mov)), log(prob), 0.1, TRUE)
  }
  nll <- c(nll_dist, nll_stay) %>% sum()
  return(-nll)
}

#' Calculates movement matrices from user inputs for fraction in each area (fracs) and probability of staying in areas (prob)
#'
#' @description A function for calculating a movement matrix from user specified unfished stock biomass fraction in each area.
#' Used by \link{simmov} to generate movement matrices for an operating model.
#' @param fracs A vector nareas long of fractions of unfished stock biomass in each area
#' @param prob A vector of the probability of individuals staying in each area or a single value for the mean probability of staying among all areas
#' @author T. Carruthers
#' @export
#' @seealso \link{simmov}
makemov <- function(fracs = c(0.1, 0.2, 0.3, 0.4), prob = c(0.5, 0.8, 0.9, 0.95)) {

  nareas <- length(fracs)
  nprob <- length(prob)

  opt <- stats::nlminb(rep(0, nprob + nareas - 1), opt_mov, fracs = fracs, prob = prob, nareas = nareas,
                control = list(iter.max = 5e3, eval.max = 1e4))
  mov <- grav(log_visc = opt$par[1:length(prob)], log_grav = opt$par[(length(prob)+1):length(opt$par)],
              fracs = fracs, nareas = nareas)$mov
  return(mov)
}


# #' Checks the TMB equations again estimated movement matrices and also checks fit
# #'
# #' @description A function for calculating a movement matrix from user specified unfished stock biomass fraction in each area
# #' @param obj A list object arising from MakeADFun from either grav.h or grav_Pbyarea.h
# #' @author T. Carruthers
# #' @export validateTMB
validateTMB <- function(obj) {

  log_grav<-obj$report()$log_grav
  log_visc<-obj$report()$log_visc
  nareas<-length(log_grav)+1

  grav<-array(0,c(nareas,nareas))
  grav[,2:nareas]<-rep(log_grav,each=nareas)
  grav[cbind(1:nareas,1:nareas)]<-grav[cbind(1:nareas,1:nareas)]+log_visc

  mov<-exp(grav)/apply(exp(grav),1,sum)

  idist<-rep(1/nareas,nareas)
  for(i in 1:50)idist<-apply(idist*mov,2,sum)

  print(obj$report()$idist)
  print(idist)
  print(obj$report()$fracs)

  print(obj$report()$mov)
  print(mov)

}

#' @param mov A four-dimensional array of dimension \code{c(nsim, maxage, nareas, nareas)} or a five-dimensional
#' array of dimension \code{c(nsim, maxage, nareas, nareas, nyears + proyears)} specifying movement 
#' in the operating model.
#' @param age An age from 0 to maxage for the movement-at-age matrix figure when \code{type = "matrix"}.
#' @param type Whether to plot a movement matrix for a single age (\code{"matrix"}) or the
#' full movement versus age figure (\code{"all"})
#' @param year If \code{mov} is a 5-dimensional array, the year (from 1 to nyears + proyears) for which 
#' to plot movement.
#' @param qval The quantile to plot or report the range of values among simulations.
#' @describeIn simmov Plotting function.
#' @export
plot_mov <- function(mov, age = 1, type = c("matrix", "all"), year = 1, qval = 0.9) {
  type <- match.arg(type)
  nsim <- dim(mov)[1]
  maxage <- dim(mov)[2] - 1
  nareas <- dim(mov)[3]
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par))
  
  if(length(dim(mov)) == 5) mov <- mov[, , , , year]

  if(type == "matrix") {
    mov_at_age <- mov[ , age + 1, , ] # dimension nsim, narea, narea
    q_mov <- apply(mov_at_age, c(2, 3), quantile, probs = c(0.5 * (1 - qval), 0.5, qval + 0.5 * (1 - qval)))
    
    plot(NULL, NULL, xlab = "Area (from)", ylab = "Area (to)", 
         xaxp = c(1, nareas, nareas - 1), yaxp = c(1, nareas, nareas - 1),
         xaxs = "i", yaxs = "i",
         xlim = c(0.5, nareas + 0.5), ylim = c(0.5, nareas + 0.5))
    abline(h = 1.5:(nareas - 0.5), v = 1.5:(nareas - 0.5))
    
    if(identical(q_mov[1, , ], q_mov[3, , ])) {
      for(i in 1:nareas) {
        for(j in 1:nareas) {
          text(i, j, labels = signif(q_mov[2, i, j], 2))
        }
      }
      title(paste0("Movement matrix at age ", age))
    } else {
      for(i in 1:nareas) {
        for(j in 1:nareas) {
          text(i, j, labels = paste0(signif(q_mov[2, i, j], 2), "\n(", signif(q_mov[1, i, j], 2), "-", signif(q_mov[3, i, j], 2), ")"))
        }
      }
      title(paste0("Movement matrix at age ", age, "\n(Median and ", 100 * qval, "% quantile over ", nsim, " simulations)"))
    }
    
    
  }

  if(type == "all") {
    q_mov <- apply(mov, 2:4, quantile, probs = c(0.5 * (1 - qval), 0.5, qval + 0.5 * (1 - qval)))

    n_plots <- nareas + 1
    ylim <- c(0, max(1, 1.1 * max(q_mov)))
    color.vec <- grDevices::rainbow(nareas)
    
    if(nareas <= 3) {
      layout_matrix <- matrix(1:4, 2, 2, byrow = TRUE) %>%
        cbind(rep(0, 2))
      layout(layout_matrix + 1, widths = c(1, 1, 0.5))
    } else {
      layout_matrix <- matrix(seq_len(ceiling(n_plots/3) * 3), ceiling(n_plots/3), 3, byrow = TRUE) %>%
        cbind(rep(0, ceiling(n_plots/3)))
      layout(layout_matrix + 1, widths = c(1, 1, 1, 0.8))
    }
    par(mar = c(2, 4, 1, 1), oma = c(4, 0, 0, 0))
    
    # Legend
    plot(1, 1, typ = "n", axes = FALSE, xlab = "", ylab = "")
    legend("left", paste("Area", 1:nareas), title = "Movement to:",
           col = color.vec, pch = 16, lty = 1, bty = "n")
    
    # Prob staying
    plot(NULL, NULL, xlim = c(0, maxage), ylim = ylim, ylab = "Probability of staying")
    abline(h = 0, col = "grey")
    lapply(1:nareas, function(i) {
      polygon(c(0:maxage, maxage:0), c(q_mov[1, , i, i], rev(q_mov[3, , i, i])), 
              col = makeTransparent(color.vec[i], 80), border = NA)
      lines(0:maxage, q_mov[2, , i, i], col = color.vec[i], lwd = 3, pch = 16)
    })
    
    # Each figure shows movement from one area to all areas
    for(j in 1:nareas) {
      plot(NULL, NULL, xlim = c(0, maxage), ylim = ylim, ylab = paste("Movement from Area", j))
      abline(h = 0, col = "grey")
      lapply(1:nareas, function(i) {
        polygon(c(0:maxage, maxage:0), c(q_mov[1, , j, i], rev(q_mov[3, , j, i])), 
                col = makeTransparent(color.vec[i], 80), border = NA)
        lines(0:maxage, q_mov[2, , j, i], col = color.vec[i], lwd = 3, pch = 16)
      })
    }
    mtext("Age", side = 1, line = 2, outer = TRUE)
    
  }

  invisible()
}