R/apollo_classAlloc.R

In byu-transpolab/apollo-byu: Tools for Choice Model Estimation and Application

#' Calculates Multinomial Logit probabilities
#'
#' Calculates probabilities of a Multinomial Logit model.
#'
#' @param classAlloc_settings List of inputs of the MNL model. It should contain the following.
#'                     \itemize{
#'                       \item \strong{\code{V}}: Named list of deterministic utilities . Utilities of the alternatives. Names of elements must match those in \code{avail}, if provided.
#'                       \item \strong{\code{avail}}: Named list of numeric vectors or scalars. Availabilities of alternatives, one element per alternative. Names of elements must match those in \code{V}. Values can be 0 or 1.
#'                       \item \strong{\code{rows}}: Boolean vector. Consideration of rows in the likelihood calculation, FALSE to exclude. Length equal to the number of observations (nObs). Default is \code{"all"}, equivalent to \code{rep(TRUE, nObs)}.
#'                       \item \strong{\code{componentName}}: Character. Name given to model component.
#'                     }
#' @param functionality Character. Can take different values depending on desired output.
#'                      \itemize{
#'                        \item \code{"estimate"}: Used for model estimation.
#'                        \item \code{"prediction"}: Used for model predictions.
#'                        \item \code{"validate"}: Used for validating input.
#'                        \item \code{"zero_LL"}: Used for calculating null likelihood.
#'                        \item \code{"conditionals"}: Used for calculating conditionals.
#'                        \item \code{"output"}: Used for preparing output after model estimation.
#'                        \item \code{"raw"}: Used for debugging.
#'                      }
#' @return The returned object depends on the value of argument \code{functionality} as follows.
#'         \itemize{
#'           \item \strong{\code{"estimate"}}: vector/matrix/array. Returns the probabilities for the chosen alternative for each observation.
#'           \item \strong{\code{"prediction"}}: List of vectors/matrices/arrays. Returns a list with the probabilities for all alternatives, with an extra element for the probability of the chosen alternative.
#'           \item \strong{\code{"validate"}}: Same as \code{"estimate"}, but it also runs a set of tests to validate the function inputs.
#'           \item \strong{\code{"zero_LL"}}: vector/matrix/array. Returns the probability of the chosen alternative when all parameters are zero.
#'           \item \strong{\code{"conditionals"}}: Same as \code{"estimate"}
#'           \item \strong{\code{"output"}}: Same as \code{"estimate"} but also writes summary of input data to internal Apollo log.
#'           \item \strong{\code{"raw"}}: Same as \code{"prediction"}
#'         }
#' @export
#' @importFrom utils capture.output
apollo_classAlloc <- function(classAlloc_settings){
  # Fetch functionality
  functionality = tryCatch(get('functionality'), parent.frame(), inherits=TRUE, error=function(e) return('estimate'))
  
  # Fetch apollo_inputs
  apollo_inputs = tryCatch(get("apollo_inputs", parent.frame(), inherits=FALSE),
                           error=function(e) return( list(apollo_control=list(cpp=FALSE, silent=FALSE, analyticGrad=TRUE),
                                                          silent=FALSE) ))
  
  ### Set or extract componentName
  modelType = "classAlloc"
  if(is.null(classAlloc_settings[["componentName"]])){
    classAlloc_settings[["componentName"]] = ifelse(!is.null(classAlloc_settings[['componentName2']]),
                                             classAlloc_settings[['componentName2']], modelType)
    test <- functionality=="validate" && classAlloc_settings[["componentName"]]!='model' && !apollo_inputs$silent
    if(test) apollo_print(paste0('Apollo found a model component of type ', modelType, ' without a componentName.', 
                                 ' The name was set to "', classAlloc_settings[["componentName"]], '" by default.'))
  }
  ### Check for duplicated modelComponent name
  if(functionality=="validate"){
    apollo_modelList <- tryCatch(get("apollo_modelList", envir=parent.frame(), inherits=FALSE), error=function(e) c())
    apollo_modelList <- c(apollo_modelList, classAlloc_settings$componentName)
    if(anyDuplicated(apollo_modelList)) stop("Duplicated componentName found (", classAlloc_settings$componentName,
                                             "). Names must be different for each component.")
    assign("apollo_modelList", apollo_modelList, envir=parent.frame())
  }
  
  # ############################### #
  #### Load or do pre-processing ####
  # ############################### #
  if( !is.null(apollo_inputs[[paste0(classAlloc_settings$componentName, "_settings")]]) && (functionality!="preprocess") ){
    # Load classAlloc_settings from apollo_inputs
    tmp <- apollo_inputs[[paste0(classAlloc_settings$componentName, "_settings")]]
    # If there is no V inside the loaded classAlloc_settings, restore the one received as argument
    if(is.null(tmp$V)) tmp$V <- classAlloc_settings$V
    classAlloc_settings <- tmp
    rm(tmp)
    
  } else { 
    ### Do pre-processing
    # Do pre-processing common to most models
    classAlloc_settings <- apollo_preprocess(inputs = classAlloc_settings, modelType, functionality, apollo_inputs)
    
    # Determine which mnl likelihood to use (R or C++)
    if(apollo_inputs$apollo_control$cpp & !apollo_inputs$silent) apollo_print("No C++ optimisation for classAlloc available")
    # Using R likelihood
    classAlloc_settings$probs_MNL <- function(classAlloc_settings){
      # Set utility of unavailable alternatives to 0 to avoid numerical issues (eg attributes = -999)
      classAlloc_settings$V <- mapply(function(v,a) apollo_setRows(v, !a, 0), 
                               classAlloc_settings$V, classAlloc_settings$avail, SIMPLIFY=FALSE)
      # subtract the maxV
      maxV <- do.call(pmax, classAlloc_settings$V)
      classAlloc_settings$V <- lapply(classAlloc_settings$V, "-", maxV)
      # calculate probabilities of all alternatives
      classAlloc_settings$V <- lapply(X=classAlloc_settings$V, FUN=exp)
      classAlloc_settings$V <- mapply('*', classAlloc_settings$V, classAlloc_settings$avail, SIMPLIFY = FALSE)
      denom = Reduce('+',classAlloc_settings$V)
      P <- lapply(classAlloc_settings$V, "/", denom)
      return(P)
    }
    
    # Construct necessary input for gradient (including gradient of utilities)
    apollo_beta <- tryCatch(get("apollo_beta", envir=parent.frame(), inherits=TRUE),
                            error=function(e) return(NULL))
    test <- !is.null(apollo_beta) && (functionality %in% c("preprocess", "gradient"))
    test <- test && all(sapply(classAlloc_settings$V, is.function))
    test <- test && apollo_inputs$apollo_control$analyticGrad
    classAlloc_settings$gradient <- FALSE
    if(test){
      classAlloc_settings$dV       <- apollo_dVdB(apollo_beta, apollo_inputs, classAlloc_settings$V)
      classAlloc_settings$gradient <- !is.null(classAlloc_settings$dV)
    }; rm(test)
    
    # Return classAlloc_settings if pre-processing
    if(functionality=="preprocess"){
      # Remove things that change from one iteration to the next
      classAlloc_settings$V      <- NULL
      return(classAlloc_settings)
    }
  }
  
  # ############################################ #
  #### Transform V into numeric and drop rows ####
  # ############################################ #
  
  ### Execute V (makes sure we are now working with vectors/matrices/arrays and not functions)
  if(any(sapply(classAlloc_settings$V, is.function))){
    classAlloc_settings$V = lapply(classAlloc_settings$V, function(f) if(is.function(f)) f() else f )
  } 
  classAlloc_settings$V <- lapply(classAlloc_settings$V, function(v) if(is.matrix(v) && ncol(v)==1) as.vector(v) else v)
  
  ### Drop rows from V if necessary
  if(!all(classAlloc_settings$rows)) classAlloc_settings$V <- lapply(classAlloc_settings$V, apollo_keepRows, r=classAlloc_settings$rows)
  # No need to drop rows in avail, as it was already filtered durin pre-processing
  
  # ############################## #
  #### functionality="validate" ####
  # ############################## #
  
  if(functionality=="validate"){
    if(!apollo_inputs$apollo_control$noValidation) apollo_validate(classAlloc_settings, modelType, 
                                                                   functionality, apollo_inputs)
    
    # No diagnose for the class allocation component
    #if(!apollo_inputs$apollo_control$noDiagnostics) apollo_diagnostics(classAlloc_settings, modelType, apollo_inputs)
    
    testL = classAlloc_settings$probs_MNL(classAlloc_settings)
    if(all(testL==0)) stop("All observations have zero probability at starting value for model component \"",classAlloc_settings$componentName,"\"")
    if(any(testL==0) && !apollo_inputs$silent && apollo_inputs$apollo_control$debug) apollo_print(paste0("Some observations have zero probability at starting value for model component \"",classAlloc_settings$componentName,"\"", sep=""))
    return(invisible(testL))
  }
  
  # ############################## #
  #### functionality="zero_LL" ####
  # ############################## #
  
  if(functionality=="zero_LL") return(NA)
  
  # ################################################################# #
  #### functionality="estimate/prediction/conditionals/raw/output" ####
  # ################################################################# #
  
  if(functionality %in% c("estimate","conditionals", "components", "output", "prediction", "raw")){
    P <- classAlloc_settings$probs_MNL(classAlloc_settings)
    if(!all(classAlloc_settings$rows)) P <- lapply(P, apollo_insertRows, r=classAlloc_settings$rows, val=1)
    P <- lapply(P, apollo_firstRow, apollo_inputs=apollo_inputs)
    return(P)
  }
  
  # ############################## #
  #### functionality="gradient" ####
  # ############################## #
  
  if(functionality=="gradient"){
    # Verify everything necesary is available
    if(is.null(classAlloc_settings$dV) || !all(sapply(classAlloc_settings$dV, is.function))) stop("Analytical gradient could not be calculated. Please set apollo_control$analyticGrad=FALSE.")
    apollo_beta <- tryCatch(get("apollo_beta", envir=parent.frame(), inherits=TRUE),
                            error=function(e) stop("apollo_mnl could not fetch apollo_beta for gradient estimation."))
    if(is.null(apollo_inputs$database)) stop("apollo_mnl could not fetch apollo_inputs$database for gradient estimation.")
    
    # Calculate probabilities and derivatives of utilities for all alternatives
    P <- classAlloc_settings$probs_MNL(classAlloc_settings)
    e <- list2env(c(as.list(apollo_beta), apollo_inputs$database, list(apollo_inputs=apollo_inputs)), hash=TRUE)
    for(i in 1:length(classAlloc_settings$dV)) environment(classAlloc_settings$dV[[i]]) <- e
    dV<- lapply(classAlloc_settings$dV, function(dv) dv())
    if(!all(classAlloc_settings$rows)) for(i in 1:length(dV)) dV[[i]] <- lapply(dV[[i]], apollo_keepRows, classAlloc_settings$rows)
    for(i in 1:classAlloc_settings$nAlt) dV[[i]] <- lapply(dV[[i]], 
                                                    function(dvik){ # Make dV=0 for unavailable alternatives
                                                      test <- length(dvik)==1 && length(classAlloc_settings$avail[[i]])>1
                                                      if(test) dvik <- rep(dvik, classAlloc_settings$nObs)
                                                      dvik[!classAlloc_settings$avail[[i]]] <- 0
                                                      return(dvik)
                                                    })
    
    # Calculate gradient for each alternative
    K   <- length(dV[[1]])
    Pd <- list()
    for(k in 1:K) Pd[[k]] <- Reduce('+', mapply('*', P, lapply(dV, function(dv) dv[[k]]), SIMPLIFY=FALSE))
    G <- list()
    for(i in 1:classAlloc_settings$nAlt){
      G[[i]] <- mapply(function(dvik, pdk) dvik - pdk, dV[[i]], Pd, SIMPLIFY=FALSE)
      G[[i]] <- mapply('*', P, G[[i]], SIMPLIFY=FALSE)
    }
    
    # Restore rows and return
    if(!all(classAlloc_settings$rows)){
      P <- lapply(P, apollo_insertRows, r=classAlloc_settings$rows, val=1)
      for(i in 1:classAlloc_settings$nAlt) G[[i]] <- lapply(G[[i]], apollo_insertRows, r=classAlloc_settings$rows, val=0)
    }
    
    # Change order and return
    ans <- list()
    for(i in 1:classAlloc_settings$nAlt) ans[[i]] <- list(like=P[[i]], grad=G[[i]])
    return(ans)
  }
  
  # ############ #
  #### Report ####
  # ############ #
  if(functionality=='report') return(list(data='', param=''))
  
}