R/Build_TMB_Fn.R

In nwfsc-assess/geostatistical_delta-GLMM: Package for conducting spatial estimation of CPUE index standardization (fisheries, ecology, etc)

#' Build TMB object for geostatistical delta-GLMM
#'
#' \code{Build_TMB_Fn} builds a tagged list with everything necessary to run or interpret inputs to a geostatistical delta-GLMM
#'
#' @param TmbData, a tagged list of data inputs generated by \code{Data_Fn}
#' @param VesselConfig, a vector of form c("Vessel"=0,"VesselYear"=0) turning off vessel or vessel-year effects
#' @param Q_Config, a boolean whether to estimate catchability covariates (Q_ik in the TmbData input) or not
#' @param CovConfig, a boolean whether to estimate density covariates (X_xj and X_xtp in the TmbData input) or not
#' @param Method Spatial method used for estimation (determines bounds for logkappa)
#' @param ConvergeTol, OPTIONAL override for TMB convergence criteria
#' @param Use_REML, OPTIONAL boolean whether to use maximum marginal likelihood or restricted maximum likelihood (termed "REML")
#' @param loc_x OPTIONAL, location for each sample used to generate plausible bounds for scale parameter
#' @param Parameters OPTIONAL, a tagged list of starting parameters
#' @param Random OPTIONAL, a character vector of random effects
#' @param Map OPTIONAL, a tagged list of parameters to either mirror or turn off
#' @param DiagnosticDir OPTIONAL, a directory where diagonstic runtime information should be stored
#' @param TmbDir OPTIONAL, a directory where the CPP file for the VAST model can be found locally
#' @param RunDir OPTIONAL, a directory where the CPP file is copied, copiled, and run (must have write privileges or else the function will crash)
#' @param silent Boolean, whether TMB should run silently
#' @inheritParams Data_Fn

#' @return Tagged list containing objects for running a VAST model
#' \describe{
#'   \item{Obj}{The built TMB object}
#'   \item{Upper}{A vector of upper bounds for parameters, optionally for use during optimization}
#'   \item{Lower}{A vector of lower bounds for parameters, optionally for use during optimization}
#'   \item{Parameters}{A tagged list of parameter starting values used when building Obj, which can be extracted, modified, and then put back into \code{Build_TMB_Fn} to define different starting values}
#'   \item{Map}{A taggged list of parameters to be turned off or mirrored, for similar use as Parameters}
#'   \item{Random}{A character vector of random effects, for similar use as Parameters}
#' }

#' @importFrom grDevices colorRampPalette dev.off jpeg png rainbow rgb tiff
#' @importFrom graphics abline arrows axis box hist legend lines matplot mtext par plot plot.new points polygon rect title
#' @importFrom stats cor dist na.omit plogis pnorm qlogis quantile rbinom rgamma rlnorm rnorm runif sd var
#' @importFrom utils capture.output write.csv write.table

#' @export
Build_TMB_Fn <-
function( TmbData, Version, VesselConfig=c("Vessel"=0,"VesselYear"=0), Q_Config=TRUE, CovConfig=TRUE,
  RhoConfig=c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0), Method="Mesh",
  ConvergeTol=1, Use_REML=FALSE, loc_x=NULL, Parameters="generate", Random="generate", Map="generate",
  DiagnosticDir=NULL, TmbDir=system.file("executables",package="SpatialDeltaGLMM"), RunDir=getwd(), silent=FALSE ){

  # Compile TMB software
  if( !file.exists(paste0(TmbDir,"/",Version,".cpp")) ){
    stop("Make sure that ",TmbDir," contains ",Version,".cpp")
  }
  file.copy( from=paste0(TmbDir,"/",Version,".cpp"), to=paste0(RunDir,"/",Version,".cpp"), overwrite=FALSE)
  setwd( RunDir )
  compile( paste0(Version,".cpp") )

  # Local functions
  boundsifpresent_fn = function( par, map, name, lower, upper, bounds ){
    if( name %in% names(par) ){
      bounds[grep(name,names(par)),c('Lower','Upper')] = c(lower,upper)
    }
    return( bounds )
  }

  # Parameters
  if( length(Parameters)==1 && Parameters=="generate" ) Parameters = Param_Fn( Version=Version, DataList=TmbData, RhoConfig=RhoConfig )

  # Which parameters are turned off
  if( length(Map)==1 && Map=="generate" ) Map = Make_Map( Version=Version, TmbData=TmbData, VesselConfig=VesselConfig, CovConfig=CovConfig, Q_Config=Q_Config, RhoConfig=RhoConfig, Aniso=TmbData[['Options_vec']]['Aniso'])
  if( "hyperparameters_z"%in%names(Parameters) && TmbData$Options_vec['AreaAbundanceCurveTF']==0 ) Map[["hyperparameters_z"]] = factor( rep(NA,length(Parameters[["hyperparameters_z"]])) )

  # Which are random
  if( length(Random)==1 && Random=="generate" ){
    Random = c("Epsiloninput1_st", "Omegainput1_s", "Epsiloninput2_st", "Omegainput2_s", "nu1_v", "nu2_v", "nu1_vt", "nu2_vt")
    if( RhoConfig[["Beta1"]]!=0 ) Random = c(Random, "beta1_t")
    if( RhoConfig[["Beta2"]]!=0 ) Random = c(Random, "beta2_t")
    if( Use_REML==TRUE ){
      Random = union(Random, c("beta1_t","gamma1_j","gamma1_tp","lambda1_k","beta2_t","gamma2_j","gamma2_tp","lambda2_k"))
      Random = Random[which(Random %in% names(Parameters))]
    }
    # Avoid problems with mapping
    Random = setdiff(Random, names(Map))
    if( length(Random)==0) Random = NULL
  }

  # Build object
  dyn.load( paste0(RunDir,"/",TMB::dynlib(Version)) ) # random=Random,
  Obj <- MakeADFun(data=TmbData, parameters=Parameters, hessian=FALSE, map=Map, random=Random, DLL=Version, silent=silent)  #
  Obj$control <- list(trace=1, parscale=1, REPORT=1, reltol=1e-12, maxit=100)

  # Diagnostic functions (optional)
  if( !is.null(DiagnosticDir) ){
    Obj$gr_orig = Obj$gr
    Obj$fn_orig = Obj$fn
    Obj$fn = function( vec ){
      capture.output( matrix(vec,ncol=1,dimnames=list(names(Obj$par),NULL)), file=paste0(DiagnosticDir,"fn.txt") )
      write.table( matrix(vec,nrow=1), row.names=FALSE, sep=",", col.names=FALSE, append=TRUE, file=paste0(DiagnosticDir,"trace.csv"))
      return( Obj$fn_orig(vec) )
    }
    Obj$gr = function( vec ){
      capture.output( matrix(vec,ncol=1,dimnames=list(names(Obj$par),NULL)), file=paste0(DiagnosticDir,"gr.txt") )
      return( Obj$gr_orig(vec) )
    }
    write.table( matrix(Obj$par,nrow=1), row.names=FALSE, sep=",", col.names=FALSE, file=paste0(DiagnosticDir,"trace.csv"))
  }

  # Declare upper and lower bounds for parameter search
  Bounds = matrix( NA, ncol=2, nrow=length(Obj$par), dimnames=list(names(Obj$par),c("Lower","Upper")) )
  Bounds[,'Lower'] = rep(-50, length(Obj$par))
  Bounds[,'Upper'] = rep( 50, length(Obj$par))
  Bounds[grep("logsigmaV",names(Obj$par)),'Lower'] = log(0.01)
  Bounds[grep("logsigmaVT",names(Obj$par)),'Lower'] = log(0.01)
  Bounds[grep("logtau",names(Obj$par)),'Upper'] = 10   # Version < v2i
  Bounds[grep("logeta",names(Obj$par)),'Upper'] = log(1/(1e-2*sqrt(4*pi))) # Version >= v2i: Lower bound on margSD = 1e-4
  Bounds[grep("SigmaM",names(Obj$par)),'Upper'] = 10 # ZINB can crash if it gets > 20
  if( !is.null(loc_x) && !is.na(TmbData$Options_vec['Method']) && TmbData$Options_vec['Method']==0 && Method!="Spherical_mesh" ){
    Dist = dist(loc_x)
    Bounds[grep("logkappa",names(Obj$par)),'Lower'] = log( sqrt(8)/max(Dist) ) # Range = nu*sqrt(8)/kappa
    Bounds[grep("logkappa",names(Obj$par)),'Upper'] = log( sqrt(8)/min(Dist) ) # Range = nu*sqrt(8)/kappa
  }
  if( !is.na(TmbData$Options_vec['Method']) && TmbData$Options_vec['Method']==1 && Method!="Spherical_mesh" ){
    Bounds[grep("logkappa",names(Obj$par)),'Upper'] = log(0.9999) # Has to be negative because exp(logkappa)=Rho < 1
  }
  Bounds = boundsifpresent_fn( par=Obj$par, name="gamma1", lower=-20, upper=20, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="gamma2", lower=-20, upper=20, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="lambda1", lower=-20, upper=20, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="lambda2", lower=-20, upper=20, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Beta_rho1", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Beta_rho2", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho1", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho2", lower=-0.99, upper=0.99, bounds=Bounds)

  # Change convergence tolerance
  Obj$env$inner.control$step.tol <- c(1e-8,1e-12,1e-15)[ConvergeTol] # Default : 1e-8  # Change in parameters limit inner optimization
  Obj$env$inner.control$tol10 <- c(1e-6,1e-8,1e-12)[ConvergeTol]  # Default : 1e-3     # Change in pen.like limit inner optimization
  Obj$env$inner.control$grad.tol <- c(1e-8,1e-12,1e-15)[ConvergeTol] # # Default : 1e-8  # Maximum gradient limit inner optimization

  # Print number of parameters
  message("Number of fixed and random effects:")
  print( table(names(Obj$env$last.par)) )

  # Return stuff
  Return = list("Obj"=Obj, "Upper"=Bounds[,'Upper'], "Lower"=Bounds[,'Lower'], "Parameters"=Parameters, "Map"=Map, "Random"=Random)
  return( Return )
}