R/Movement_calcs.R

In tcarruth/MSEtool: Management Strategy Evaluation Toolkit

# Movement matrix calculation functions
#' Calculates movement matrices from user inputs
#'
#' A wrapper function for \link{makemov} used to generate movement matrices for a DLMtool 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.)
#' @param mov A four-dimensional array of dimension \code{c(nsim, maxage, nareas, nareas)}
#' specifying movement in the operating model.
#' @param age An age from 1 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"})
#' @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
#' @import TMB
#' @useDynLib MSEtool
#' @examples
#' \donttest{
#' 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

  mov<-array(NA,c(nsim,maxage,nareas,nareas))

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

  if(length(prob)==1){

    prob_s<-ilogit(matrix(rnorm(nsim,logit(prob),probE),nrow=nsim))

  }else if(length(prob)==length(dist)&is.na(prob2)){

    prob_s<-ilogit(t(matrix(rnorm(nareas*nsim,logit(prob),probE),nrow=nareas)))

  }else if(length(prob)==length(dist)&length(prob)==length(prob2)){

    prob_s<-t(matrix(runif(nareas*nsim,prob,prob2),nrow=nareas))

  }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")
  }

  for(i in 1:nsim){

    movt<-makemov(fracs=dist_s[i,],prob=prob_s[i,])
    mov[i,,,]<-rep(movt,each=maxage)
  } # nsim

  OM@cpars$mov<-mov
  if(figure)plot_mov(mov)
  OM

}
# myOM<-simmov(testOM, dist=c(0.1,0.8,0.05,0.05),prob=c(0.6,0.5,0.7,0.85))


#' 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 a DLMtool 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 makemov
#' @import TMB
#' @importFrom stats nlminb
#' @useDynLib MSEtool
#' @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)

  if(length(prob)==1) data <- list(model = "grav",fracs = fracs, prob = prob, nareas = nareas)
  if(length(prob)==nareas) data <- list(model = "grav_Pbyarea",fracs = fracs, prob = prob, nareas = nareas)

  if(length(prob)==1)params <- list(log_visc = 0,log_grav = rep(0,nareas-1))
  if(length(prob)==nareas)params <- list(log_visc = rep(0,nareas),log_grav = rep(0,nareas-1))

  info <- list(data = data, params = params)

  obj <- MakeADFun(data = info$data, parameters = info$params, DLL = "MSEtool", silent = TRUE)
  opt <- nlminb(start = obj$par, objective = obj$fn, gradient = obj$gr)

  if(length(prob)==1)params.new<-list(log_visc=opt$par[1],log_grav=opt$par[2:nareas])
  if(length(prob)==nareas)params.new<-list(log_visc=opt$par[1:nareas],log_grav=opt$par[nareas+(1:(nareas-1))])

  info <- list(data = data, params = params.new)
  obj.new<-MakeADFun(data = info$data, parameters = info$params, DLL = "MSEtool", silent = TRUE)

  obj.new$report(obj.new$env$last.par.best)$mov
  #validateTMB(obj.new)

}

# #' 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
# #' @import TMB
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)

}

#' @describeIn simmov Plotting function.
#' @export
plot_mov <- function(mov, age = 1, type = c("matrix", "all")) {
  type <- match.arg(type)
  nsim <- dim(mov)[1]
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par))

  if(type == "matrix") {
    mov_at_age <- mov[ , age, , ] # dimension nsim, narea, narea
    mean_mov <- apply(mov_at_age, c(2, 3), mean)
    nareas <- nrow(mean_mov)

    plot(NULL, NULL, xlab = "Area (to)", ylab = "Area (from)", 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))
    for(i in 1:nareas) {
      for(j in 1:nareas) {
        text(j, i, labels = paste0(signif(mean_mov[i, j], 2), "\n(", signif(min(mov_at_age[, i, j]), 2), "-", signif(max(mov_at_age[, i, j]), 2), ")"))
      }
    }
    title(paste("Movement matrix at age", age, "\n(mean and range across", nsim, "simulations)"))
  }

  if(type == "all") {
    mean_mov <- apply(mov, 2:4, mean)
    maxage <- dim(mean_mov)[1]
    nareas <- dim(mean_mov)[2]

    n_plots <- nareas + 1
    ylim <- c(0, max(1, 1.1 * max(mean_mov)))
    color.vec <- rich.colors(nareas)

    par(mfrow = c(ceiling(n_plots/3), 3), mar = c(2, 5, 1, 1), oma = c(5, 0, 2, 0))
    # Prob staying
    plot(NULL, NULL, xlim = c(1, maxage), ylim = ylim, ylab = "Mean probability of staying")
    abline(h = 0, col = "grey")
    for(i in 1:nareas) lines(1:maxage, mean_mov[, i, i], col = color.vec[i], typ = 'o', pch = 16)
    legend("topleft", paste("Area", 1:nareas), col = color.vec, pch = 16, lty = 1)


    # Each figure shows movement from one area to all areas
    for(j in 1:nareas) {
      plot(NULL, NULL, xlim = c(1, maxage), ylim = ylim, ylab = paste("Mean movement probability\nfrom Area", j))
      abline(h = 0, col = "grey")
      for(i in 1:nareas) lines(1:maxage, mean_mov[, j, i], col = color.vec[i], typ = 'o', pch = 16)
    }
    mtext("Age", side = 1, line = 2, outer = TRUE)
  }

  invisible()
}