R/make_model.R

In James-Thorson/VAST: Vector-Autoregressive Spatio-Temporal (VAST) Model

#' Build TMB object for VAST model
#'
#' \code{make_model} builds a tagged list with everything necessary to run or interpret inputs to a VAST model
#'
#' @inheritParams make_data
#' @inheritParams make_map
#' @inheritParams TMB::MakeADFun
#' @inheritParams TMB::compile
#' @param TmbData a tagged list of data inputs generated by \code{Data_Fn}
#' @param Method Spatial method used for estimation (determines bounds for logkappa)
#' @param ConvergeTol Integer specifying override for TMB convergence criteria (OPTIONAL)
#' @param Use_REML boolean whether to use maximum marginal likelihood or restricted maximum likelihood (termed "REML")  (OPTIONAL)
#' @param loc_x location for each sample used to generate plausible bounds for scale parameter  (OPTIONAL)
#' @param Parameters a tagged list of starting parameters  (OPTIONAL)
#' @param Random a character vector of random effects  (OPTIONAL)
#' @param Map a tagged list of parameters to either mirror or turn off, using standard TMB interface.
#'        This input is useful, e.g., to build a model without estimating parameters, extracting Map
#'        from the list of outputs, modifying it manually, and then passing it explicitly, \code{make_model(Map=NewMap)} (OPTIONAL)
#' @param DiagnosticDir a directory where diagonstic runtime information should be stored (OPTIONAL)
#' @param TmbDir a directory where the CPP file for the VAST model can be found locally (OPTIONAL)
#' @param RunDir a directory where model results are written; by default uses the working directory (OPTIONAL)
#' @param CompileDir a directory where the CPP file is copied, copiled, and run
#'        (must have write privileges or else the function will crash); by default uses \code{TmbDir} (OPTIONAL)
#' @param build_model Boolean indicating whether to build the model, \code{build_model=TRUE}, or simply build the inputs, \code{build_model=FALSE}

#' @return Object of class \code{make_model}, 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{make_model} to define different starting values}
#'   \item{Map}{A tagged 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 graphics axis box hist image mtext par plot text title
#' @importFrom grDevices colorRampPalette dev.off png
#' @importFrom stats cov cov2cor na.omit nlminb qlogis rnorm var weighted.mean
#' @importFrom utils capture.output

#' @export
make_model <-
function( TmbData,
          Version,
          RhoConfig = c("Beta1" = 0,"Beta2" = 0,"Epsilon1" = 0,"Epsilon2" = 0),
          Method = "Mesh",
          Npool = 0,
          ConvergeTol = 1,
          Use_REML = FALSE,
          loc_x = NULL,
          Parameters = "generate",
          Random = "generate",
          Map = "generate",
          DiagnosticDir = NULL,
          TmbDir = system.file("executables",package = "VAST"),
          RunDir = getwd(),
          CompileDir = TmbDir,
          build_model = TRUE,
          framework = "TMBad",
          intern = FALSE,
          inner.control = list(sparse=TRUE, lowrank=TRUE, trace=FALSE),
          supernodal = FALSE,
          flags = "" ){

  # Extract Options and Options_vec (depends upon version)
  if( all(c("Options","Options_vec") %in% names(TmbData)) ){
    Options_vec = TmbData$Options_vec
    Options = TmbData$Options
  }
  if( "Options_list" %in% names(TmbData) ){
    Options_vec = TmbData$Options_list$Options_vec
    Options = TmbData$Options_list$Options
  }

  # Augment objects in TmbData (to deal with backwards compatibility)
  if( !("n_e" %in% names(TmbData)) ){
    TmbData[["n_e"]] = TmbData$n_c
  }
  if( !("ObsModel_ez" %in% names(TmbData)) ){
    TmbData[["ObsModel_ez"]] = rep(1,TmbData[["n_e"]]) %o% TmbData$ObsModel
  }
  if( !("c_iz" %in% names(TmbData)) ){
    TmbData[["c_iz"]] = matrix( TmbData$c_i, ncol=1 )
  }
  if( !("e_i" %in% names(TmbData)) ){
    TmbData[["e_i"]] = TmbData$c_iz[,1]
  }
  if( !("t_iz" %in% names(TmbData)) ){
    TmbData[["t_iz"]] = matrix( TmbData$t_i, ncol=1 )
  }

  # Save package version info
  capture.output( packageDescription("VAST"), file=paste0(RunDir,"/packageDescription.txt") )
  capture.output( packageDescription("FishStatsUtils"), file=paste0(RunDir,"/packageDescription.txt"), append=TRUE )

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

  # Which parameters are turned off
  if( length(Map)==1 && Map=="generate" ){
    Map = make_map( DataList = TmbData,
                    TmbParams = Parameters,
                    RhoConfig = RhoConfig,
                    Npool = Npool )
  }else{
    warning( "Please carefully check starting values for all parameters to ensure that mapping off parameters will work as expected.")
  }

  # Which are random
    # Using redundant parameter names from all past versions, and then eliminating redundancies by checking against contents of Parameters
  if( length(Random)==1 && Random=="generate" ){
    Random = c("Epsiloninput1_sct", "Omegainput1_sc", "Epsiloninput1_sft", "Omegainput1_sf", "eta1_vf", "Xiinput1_scp", "Phiinput1_sk", "Epsiloninput1_sff",
      "Epsiloninput2_sct", "Omegainput2_sc", "Epsiloninput2_sft", "Omegainput2_sf", "eta2_vf", "Xiinput2_scp", "Phiinput2_sk", "Epsiloninput2_sff",
      "delta_i")
    if( RhoConfig[["Beta1"]]%in%c(1,2,4) ) Random = c(Random, "beta1_ct", "beta1_ft")
    if( RhoConfig[["Beta2"]]%in%c(1,2,4) ) Random = c(Random, "beta2_ct", "beta2_ft")
    if( Use_REML==TRUE ){
      Random = union(Random, c("beta1_ct","beta1_ft","gamma1_j","gamma1_tp","gamma1_ctp","lambda1_k","gamma1_cp",
        "beta2_ct","beta2_ft","gamma2_j","gamma2_tp","gamma2_ctp","lambda2_k","gamma2_cp"))
    }
    # Avoid problems with mapping
    Random = Random[which(Random %in% names(Parameters))]
    if( length(Random)==0) Random = NULL
  }

  # Bundle for debugging
  #Save = list("Map"=Map, "Data"=TmbData, "Parameters"=Parameters, "Random"=Random)
  #return(Save)
  #on.exit( return(Save) )
  #save(Save, file=paste0(RunDir,"/Save.RData"))

  #
  if( build_model==FALSE ){
    Return = list("Map"=Map, "Data"=TmbData, "Parameters"=Parameters, "Random"=Random)
    return( Return )
  }

   # Compile TMB software
  #dyn.unload( paste0(RunDir,"/",dynlib(TMB:::getUserDLL())) ) # random=Random,
  Version_framework = paste0( Version, "_", framework )
  file.copy( from=paste0(TmbDir,"/",Version,".cpp"),
    to=paste0(CompileDir,"/",Version_framework,".cpp"),
    overwrite=FALSE)
  origwd = getwd()
  on.exit(setwd(origwd),add=TRUE)
  setwd( CompileDir )
  # SEE https://github.com/kaskr/adcomp/issues/321 for flags argument
  if( "framework" %in% formalArgs(TMB::compile)){
    TMB::compile( file = paste0(Version_framework,".cpp"),
                  framework = framework,
                  flags = flags,
                  supernodal = supernodal )
  }else{
    TMB::compile( file = paste0(Version_framework,".cpp"),
                  flags = flags )
  }

  # Build object
  dyn.load( paste0(CompileDir,"/",TMB::dynlib(Version_framework)) ) # random=Random,
  if( ("framework" %in% formalArgs(TMB::compile)) && !is.null(framework) && (framework=="TMBad") ){
    Obj <- TMB::MakeADFun(
      data = lapply(TmbData,strip_units),
      parameters = Parameters,
      hessian = FALSE,
      map = Map,
      random = Random,
      inner.method = "newton",
      DLL = Version_framework,
      intern = intern,
      inner.control = inner.control )  #
  }else{
    Obj <- TMB::MakeADFun(
      data = lapply(TmbData,strip_units),
      parameters = Parameters,
      hessian = FALSE,
      map = Map,
      random = Random,
      inner.method = "newton",
      DLL = Version_framework )  #
  }
  Obj$control <- list(parscale=1, REPORT=1, reltol=1e-12, maxit=100)

  # Add normalization in
  if( Options['normalize_GMRF_in_CPP']==FALSE ){
    message("Normalizing GMRF in R using `TMB::normalize` feature")
    Obj = TMB::normalize(Obj, flag="include_data", value=FALSE)
  }

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

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

  # 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(-Inf, length(Obj$par))
  Bounds[,'Upper'] = rep( Inf, length(Obj$par))
  Bounds[grep("SigmaM",names(Obj$par)),'Upper'] = 10 # ZINB can crash if it gets > 20
  if( any(TmbData$ObsModel_ez[1,]==8) ) Bounds[grep("SigmaM",names(Obj$par)),'Upper'] = 3 # Tweedie can crash if logSigmaM gets too high
  if( !is.null(loc_x) && !is.na(Options_vec['Method']) && Options_vec['Method']==0 && Method!="Spherical_mesh" ){
    Dist = stats::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(Options_vec['Method']) && Options_vec['Method']==1 && Method!="Spherical_mesh" ){
    Bounds[grep("logkappa",names(Obj$par)),'Upper'] = log(0.9999) # Must be negative, so that Rho<1
  }
  Bounds = boundsifpresent_fn( par=Obj$par, name="ln_H_input", lower=-5, upper=5, bounds=Bounds)
  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="Beta_rho1_f", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Beta_rho2_f", 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)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho1_f", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho2_f", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="rho_c1", lower=-0.99, upper=0.99, bounds=Bounds)
  Bounds = boundsifpresent_fn( par=Obj$par, name="rho_c2", lower=-0.99, upper=0.99, bounds=Bounds)
  if( ("n_f_input"%in%names(TmbData)) && TmbData[["n_f_input"]]==0 ){
    Bounds = boundsifpresent_fn( par=Obj$par, name="L1_z", lower=c(-Inf,-0.99), upper=c(Inf,0.99), bounds=Bounds)
    Bounds = boundsifpresent_fn( par=Obj$par, name="L2_z", lower=c(-Inf,-0.99), upper=c(Inf,0.99), bounds=Bounds)
  }
  if( ("OverdispersionConfig"%in%names(TmbData)) ){
    if( TmbData[["OverdispersionConfig"]][1]==0 ) Bounds = boundsifpresent_fn( par=Obj$par, name="L1_z", lower=c(-Inf,-0.99), upper=c(Inf,0.99), bounds=Bounds)
    if( TmbData[["OverdispersionConfig"]][2]==0 ) Bounds = boundsifpresent_fn( par=Obj$par, name="L2_z", lower=c(-Inf,-0.99), upper=c(Inf,0.99), bounds=Bounds)
  }
  #for(i in 1:4){
  #  if( TmbData[["FieldConfig"]][i]==0 ){
  #    Bounds = boundsifpresent_fn( par=Obj$par, name=c("L_omega1_z","L_epsilon1_z","L_omega2_z","L_epsilon2_z")[i], lower=c(-Inf,-0.99), upper=c(Inf,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
  ThorsonUtilities::list_parameters( Obj )

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

#' Print model object generated by \code{\link{VAST}}
#'
#' @title Print model object
#' @param x Output from \code{\link{make_model}}
#' @param ... Not used
#' @return NULL
#' @method print make_model
#' @export
print.make_model <- function(x, ...)
{
  cat("make_model(.) result\n")
  cat("Starting values...\n")
  print( x$Obj$par )

  return(invisible(NULL))
}