R/apollo_lc.R

In apollo: Tools for Choice Model Estimation and Application

#' Calculates the likelihood of a latent class model
#' 
#' Given within class probabilities, and class allocation probabilities, calculates the probabilities of an Exploded Logit model and can also perform other operations based on the value of the \code{functionality} argument.
#' 
#' @param lc_settings List. Contains settings for this function. User input is required for all settings except those with a default or marked as optional. 
#'                  \itemize{
#'                    \item \strong{inClassProb}: List of probabilities. Conditional likelihood for each class. One element per class, in the same order as \code{classProb}.
#'                    \item \strong{classProb}: List of probabilities. Allocation probability for each class. One element per class, in the same order as \code{inClassProb}.
#'                    \item \strong{componentName}: Character. Name given to model component (optional).
#'                  }
#' @param apollo_inputs List grouping most common inputs. Created by function \link{apollo_validateInputs}.
#' @param functionality Character. Setting instructing Apollo what processing to apply to the likelihood function. This is in general controlled by the functions that call \code{apollo_probabilities}, though the user can also call \code{apollo_probabilities} manually with a given functionality for testing/debugging. Possible values are:
#'                      \itemize{
#'                        \item \strong{\code{"components"}}: For further processing/debugging, produces likelihood for each model component (if multiple components are present), at the level of individual draws and observations.
#'                        \item \strong{\code{"conditionals"}}: For conditionals, produces likelihood of the full model, at the level of individual inter-individual draws.
#'                        \item \strong{\code{"estimate"}}: For model estimation, produces likelihood of the full model, at the level of individual decision-makers, after averaging across draws.
#'                        \item \strong{\code{"gradient"}}: For model estimation, produces analytical gradients of the likelihood, where possible.
#'                        \item \strong{\code{"output"}}: Prepares output for post-estimation reporting.
#'                        \item \strong{\code{"prediction"}}: For model prediction, produces probabilities for individual alternatives and individual model components (if multiple components are present) at the level of an observation, after averaging across draws.
#'                        \item \strong{\code{"preprocess"}}: Prepares likelihood functions for use in estimation.
#'                        \item \strong{\code{"raw"}}: For debugging, produces probabilities of all alternatives and individual model components at the level of an observation, at the level of individual draws.
#'                        \item \strong{\code{"report"}}: Prepares output summarising model and choiceset structure.
#'                        \item \strong{\code{"shares_LL"}}: Produces overall model likelihood with constants only.
#'                        \item \strong{\code{"validate"}}: Validates model specification, produces likelihood of the full model, at the level of individual decision-makers, after averaging across draws.
#'                        \item \strong{\code{"zero_LL"}}: Produces overall model likelihood with all parameters at zero.
#'                      }
#' @return The returned object depends on the value of argument \code{functionality} as follows.
#'         \itemize{
#'           \item \strong{\code{"components"}}: Returns nothing.
#'           \item \strong{\code{"conditionals"}}: Same as \code{"estimate"}
#'           \item \strong{\code{"estimate"}}: vector/matrix/array. Returns the probabilities for the chosen alternative for each observation.
#'           \item \strong{\code{"gradient"}}: List containing the likelihood and gradient of the model component.
#'           \item \strong{\code{"output"}}: Same as \code{"estimate"} but also writes summary of input data to internal Apollo log.
#'           \item \strong{\code{"prediction"}}: List of vectors/matrices/arrays. Returns a list with the probabilities for all models components, for each class.
#'           \item \strong{\code{"preprocess"}}: Returns a list with pre-processed inputs, based on \code{lc_settings}.
#'           \item \strong{\code{"raw"}}: Same as \code{"prediction"}
#'           \item \strong{\code{"report"}}: Class allocation overview.
#'           \item \strong{\code{"shares_LL"}}: Same as \code{"estimate"}
#'           \item \strong{\code{"validate"}}: Same as \code{"estimate"}, but also runs a set of tests on the given arguments.
#'           \item \strong{\code{"zero_LL"}}: Same as \code{"estimate"}
#'         }
#' @importFrom utils capture.output
#' @export
apollo_lc <- function(lc_settings, apollo_inputs, functionality){
  modelType = 'LC'
  if(is.null(lc_settings[["componentName"]])){
    lc_settings[["componentName"]] = ifelse(!is.null(lc_settings[['componentName2']]),
                                            lc_settings[['componentName2']], modelType)
    test <- functionality=="validate" && lc_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 "', lc_settings[["componentName"]], '" by default.'))
  }
  
  # ############################### #
  #### Load or do pre-processing ####
  # ############################### #
  if(!is.null(apollo_inputs[[paste0(lc_settings$componentName, "_settings")]]) && (functionality!="preprocess")){
    # Load lc_settings from apollo_inputs
    tmp <- apollo_inputs[[paste0(lc_settings$componentName, "_settings")]]
    # If there is no V inside the loaded mnl_settings, restore the one received as argument
    if(is.null(tmp$inClassProb)) tmp$inClassProb <- lc_settings$inClassProb
    if(is.null(tmp$classProb  )) tmp$classProb   <- lc_settings$classProb
    lc_settings <- tmp
    rm(tmp)
  } else {
    if(functionality=="preprocess"){
      # Do preprocessing
      countOutcomes <- function(L){ # Count number of outcomes inside EACH class
        if(!is.list(L)) return(0)
        if("rows" %in% names(L)) return(sum(L$rows))
        for(i in names(L)) return( sum(sapply(L, countOutcomes)) )
      }
      lc_settings <- list(componentName = lc_settings$componentName,
                          LCNObs        = max(sapply(lc_settings$inClassProb, countOutcomes)), 
                          LCCompNames   = names(lc_settings$inClassProb),
                          modelType     = 'LC',
                          gradient      = TRUE, # ADD CHECK
                          hessian       = TRUE, # ADD CHECK
                          classAlloc_settings = lc_settings$classProb)
    }
    # Diagnostics function
    lc_settings$lc_diagnostics <- function(inputs, apollo_inputs, data=TRUE, param=TRUE){
      class_summary = matrix(0, nrow=length(inputs$classProb), ncol=1, dimnames=list(names(inputs$inClassProb), "Mean prob."))
      if(is.null(rownames(class_summary))) rownames(class_summary) <- paste0("Class_", 1:length(inputs$classProb))
      for(cc in 1:length(inputs$classProb)) class_summary[cc,1] <- mean(inputs$classProb[[cc]])
      if(!apollo_inputs$silent & param){
        apollo_print('\n')
        apollo_print(paste0('Summary of class allocation for ', toupper(inputs$modelType), ' model component ', 
                            ifelse(inputs$componentName=='model', '', inputs$componentName), ':'))
        apollo_print(class_summary)
      }
      return(invisible(TRUE))
    }
    if(functionality=="preprocess") return(lc_settings)
  }
  
  # ############################################### #
  #### functionalities with untransformed return ####
  # ############################################### #
  
  if(functionality=="components") return(NULL)
  
  ### Special case for EM estimation
  test <- !is.null(apollo_inputs$EM) && apollo_inputs$EM 
  test <- test && !is.null(apollo_inputs$class_specific) && apollo_inputs$class_specific>0
  if(test) return(lc_settings$inClassProb[[1]])
  rm(test)

  # #################### #
  #### General checks ####
  # #################### #
  if(is.null(lc_settings[["componentName"]])){
    modelType <- 'LC'
    lc_settings[["componentName"]] = ifelse(!is.null(lc_settings[['componentName2']]),
                                            lc_settings[['componentName2']], modelType)
    test <- functionality=="validate" && lc_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 "',
                                 lc_settings[["componentName"]],'" by default.'))
  }
  if(is.null(lc_settings[["inClassProb"]])) stop('SYNTAX ISSUE - The lc_settings list for model component "',
                                                 lc_settings$componentName,'" needs to include an object called "inClassProb"!')
  if(is.null(lc_settings[["classProb"]])) stop('SYNTAX ISSUE - The lc_settings list for model component "',
                                               lc_settings$componentName,'" needs to include an object called "classProb"!')
  inClassProb    = lc_settings[["inClassProb"]]
  classProb      = lc_settings[["classProb"]]
  apollo_control = apollo_inputs[["apollo_control"]]
  if(apollo_control$workInLogs && apollo_control$mixing) stop('INTERNAL ISSUE - The settings "workInLogs" and "mixing" in "apollo_control" ',
                                                              'cannot be used together for latent class models.')
  
  # Check that inClassProb is not a list of lists (if its is, use 'model' component or fail)
  if(!is.list(inClassProb)) stop('INPUT ISSUE - Setting "inClassProb" inside "lc_settings" must be a list.')
  for(i in 1:length(inClassProb)) if(is.list(inClassProb[[i]]) && functionality %in% c("conditionals","estimate","validate","zero_LL", "shares_LL", "output")){
    test <- is.null(inClassProb[[i]]$model)
    if(test) stop(paste0('INPUT ISSUE - At least one element inside "inClassProb" setting is a list. It should be a numeric vector, matrix',
                         ' or array. If you are using apollo_combineModels inside each class, try using it with the argument', 
                         ' apollo_combineModel(..., asList=FALSE)')) else inClassProb[[i]] <- inClassProb[[i]]$model
  }
  
  # Count number of rows in classProb
  Dim1 <- function(x){if(is.array(x)) return(dim(x)[[1]]) else return(length(x))}
  nRowsClassProb <- Dim1(classProb[[1]])
  if(functionality %in% c("gradient", "hessian")) nRowsClassProb <- Dim1(classProb$like[[1]])
  
  # Count number of rows in inClassProb
  nRowsInClassProb <- Dim1(inClassProb[[1]])
  if(functionality %in% c("gradient", "hessian")) nRowsInClassProb <- Dim1(inClassProb[[1]]$like)
  
  # Count number of observations per individual
  indivID <- get(apollo_control$indivID)
  nObsPerIndiv <- unique(indivID)
  nObsPerIndiv <- setNames(sapply(as.list(nObsPerIndiv),
                                  function(x) sum(indivID==x)), nObsPerIndiv)
  
  # Function to resize classProbs
  resize <- function(x, newL){
    # If it is a list, apply function recursively
    if(is.list(x)) return( lapply(x, resize, newL=newL) )
    # If no change is needed, return as is
    if(is.array(x) && length(x)==1) x <- as.vector(x)
    if(is.array(x)) oldL <- dim(x)[1] else oldL <- length(x)
    isCube <- is.array(x) && length(dim(x))==3
    if(is.matrix(x) && ncol(x)==1) x <- as.vector(x)
    if(oldL==1 | oldL==newL) return(x)
    # If longer
    if(oldL>newL){
      if(isCube){
        x <- apply(x, MARGIN=c(2,3), FUN=rowsum, group=c(1,1,2), reorder=FALSE)
      } else x <- rowsum(x, group=indivID, reorder=FALSE)
      x <- x/nObsPerIndiv
    } 
    # If shorter
    if(oldL<newL){
      if(isCube){
        tmp <- array(0, dim=c(newL, dim(x)[2:3]))
        tmp[1:nObsPerIndiv[1],,] <- x[1,,]
        if(length(nObsPerIndiv)>1) for(n in 2:length(nObsPerIndiv)){
          tmp[sum(nObsPerIndiv[1:(n-1)]):sum(nObsPerIndiv[1:n]),,] <- x[n,,]
        }; x <- tmp; rm(tmp)
      } else if(is.matrix(x)){
        tmp <- matrix(0, newL, ncol(x))
        tmp[1:nObsPerIndiv[1],] <- x[1,]
        if(length(nObsPerIndiv)>1) for(n in 2:length(nObsPerIndiv)){
          tmp[sum(nObsPerIndiv[1:(n-1)]):sum(nObsPerIndiv[1:n]),] <- x[n,]
        }; x <- tmp; rm(tmp)
      } else x <- rep(x, times=nObsPerIndiv)
    }
    # Simplify x to a vector if appropriate
    if(is.matrix(x) && ncol(x)==1) x <- as.vector(x)
    return(x)
  }
  
  # ############## #
  #### Validate ####
  # ############## #
  if(functionality=="validate"){
    
    # Validate inputs
    if(!apollo_inputs$apollo_control$noValidation){
      if(length(inClassProb)!=length(classProb)) stop("INPUT ISSUE - Arguments 'inClassProb' and 'classProb' for model component \"",
                                                      lc_settings$componentName,"\" must have the same length.")
      for(cc in 1:length(classProb)) if(sum(inClassProb[[cc]])==0) stop('CALCULATION ISSUE - Within class probabilities for model component "',
                                                                        lc_settings$componentName, 
                                                                        '" are zero for every individual for class ',
                                                                        ifelse(is.numeric(names(inClassProb)[cc]),cc,names(inClassProb)[cc]))
      classprobsum = Reduce("+",classProb)
      if(any(round(classprobsum,10)!=1)) stop("SPECIFICATION ISSUE - Class allocation probabilities in 'classProb' for model component \"",
                                              lc_settings$componentName,"\" must sum to 1.")
    }
    
    # Print diagnostics
    if(!apollo_inputs$apollo_control$noDiagnostics) lc_settings$lc_diagnostics(lc_settings, apollo_inputs)
    
    # Check if classProb needs to be averaged to match the dimensionality of inClassProb, warn if so
    test <- nRowsInClassProb!=nRowsClassProb && nRowsClassProb!=1 && nRowsInClassProb < nRowsClassProb && !apollo_inputs$silent
    #if(test) apollo_print(paste0('Class probabilities for model component "', lc_settings$componentName, 
    #                             '" averaged across observations of each individual.'))
    if(test) apollo_print(paste0('The class allocation probabilities for model component "', lc_settings$componentName, 
                                 '" are calculated at the observation level in \'apollo_lcPars\', but are used in',
                                 ' \'apollo_probabilities\' to multiply within class probabilities that are at the',
                                 ' individual level. Apollo will average the class allocation probabilities across',
                                 ' observations for the same individual level before using them to multiply the',
                                 ' within-class probabilities. If your class allocation probabilities are',
                                 ' constant across choice situations for the same individual, then this is of no concern.',
                                 ' If your class allocation probabilities however vary across choice tasks, then you should',
                                 ' change your model specification in \'apollo_probabilities\' to only call \'apollo_panelProd\'',
                                 ' after calling \'apollo_lc\'.'))
    
    testL = apollo_lc(lc_settings, apollo_inputs, functionality="estimate")
    return(testL)
  } 
  
  
  # ######################################################## #
  #### Estimate, conditionals, output, zero_LL, shares_LL ####
  # ######################################################## #
  if(functionality %in% c("estimate", "conditionals", "output", "zero_LL", "shares_LL")){
    # Match the number of rows in each list
    classProb <- resize(classProb, nRowsInClassProb)
    # # The dimension of classProb is changed if necessary, dim(inClassProb) remains the same
    # if( nRowsInClassProb!=nRowsClassProb & nRowsClassProb!=1 ){
    #     
    #   if( nRowsInClassProb < nRowsClassProb ){
    #     # If classProb is calculated at the obs level while the inClassProb is at the indiv level
    #     classProb <- lapply(classProb, rowsum, group=indivID, reorder=FALSE)
    #     classProb <- lapply(classProb, function(x) if(is.matrix(x) && ncol(x)==1) return(as.vector(x)) else return(x))
    #     classProb <- lapply(classProb, '/', nObsPerIndiv)
    #     #if(!apollo_inputs$silent) apollo_print(paste0('Class probabilities for model component "', lc_settings$componentName, 
    #     #                                              '" averaged across observations of each individual.'))
    #   } else {
    #     # If inClassProb is calculated at the obs level while the classProb is at the indiv level
    #     classProb <- lapply(classProb, rep, nObsPerIndiv)
    #   }
    # }
    
    # Calculate inClassProb*classProb
    nIndiv <- length(unique(get(apollo_control$indivID)))                                      ## FIX DP 18/05/2020
    checkForPanel <- function(x) (is.vector(x) && length(x)==nIndiv) || (is.array(x) && dim(x)[1]==nIndiv) ## FIX DP 13/08/2021
    panelProdApplied <- sapply(inClassProb, checkForPanel)                                                 ## FIX DP 13/08/2021
    if(apollo_control$workInLogs && all(panelProdApplied)){                                    ## FIX DP 18/05/2020
      # Assuming panelProd was applied inside each class, inClassProb is log(P[ns])            ## FIX DP 18/05/2020
      Bbar <- Reduce("+", inClassProb)/length(inClassProb)
      Pout <- mapply(function(pi, lnP) pi*exp(lnP - Bbar), classProb, inClassProb, SIMPLIFY=FALSE)
      Pout <- Bbar + log( Reduce("+", Pout) )
    } else {
      Pout <- mapply(function(p,pi) pi*p, inClassProb, classProb, SIMPLIFY=FALSE)
      Pout <- Reduce('+', Pout)
    }
    return( Pout )
  }
  
  # ##################### #
  #### Prediction, raw ####
  # ##################### #
  if(functionality %in% c("prediction", "raw")){
    # Match the number of rows in each list
    classProb <- resize(classProb, nRowsInClassProb)
    # # The dimension of classProb is changed if necessary, dim(inClassProb) remains the same
    # if( nRowsInClassProb!=nRowsClassProb  & nRowsClassProb!=1 ){
    #   
    #   if( nRowsInClassProb < nRowsClassProb ){
    #     # If classProb is calculated at the obs level while the inClassProb is at the indiv level
    #     stop("SPECIFICATION ISSUE - Class-probability variable for model component \"",lc_settings$componentName,"\" has more elements than in-class-probability.")
    #   } else {
    #     # If inClassProb is calculated at the obs level while the classProb is at the indiv level
    #     S=length(classProb)
    #     for(s in 1:S){
    #       isMat <- is.matrix(classProb[[s]])
    #       isCub <- is.array(classProb[[s]]) && !isMat && length(dim(classProb[[s]]))==3
    #       if(isCub){
    #         # first average out over intra
    #         classProb[[s]]=colSums(aperm(classProb[[s]], perm=c(3,1,2)))/dim(classProb[[s]])[3]
    #       } 
    #       # now check what's left
    #       isVec <- is.vector(classProb[[s]])
    #       isMat <- is.matrix(classProb[[s]])
    #       if(isVec) classProb[[s]]=rep(classProb[[s]],nObsPerIndiv)
    #       if(isMat){
    #         tmp <- matrix(0, nrow=sum(nObsPerIndiv), ncol=ncol(classProb[[s]]))
    #         for(n in 1:length(nObsPerIndiv)){
    #           a <- ifelse(n==1, 1, sum(nObsPerIndiv[1:(n-1)]) + 1)
    #           b <- a + nObsPerIndiv[n] - 1 
    #           tmp[a:b,] <- rep(classProb[[s]][n,], each=nObsPerIndiv[n])
    #         }
    #         classProb[[s]] <- tmp
    #       } 
    #     }
    #     
    #   }
    # }
    
    nClass<- length(classProb)
    for(s in 1:length(inClassProb)){
      if(!is.list(inClassProb[[s]])) inClassProb[[s]]=list(inClassProb[[s]])  
    }
    
    nAlts <- length(inClassProb[[1]])
    
    Pout <- vector(mode="list", length=nAlts)
    for(i in 1:nAlts){
      Pout[[i]] <- 0*inClassProb[[1]][[1]]
      for(k in 1:nClass) Pout[[i]] <- Pout[[i]] + inClassProb[[k]][[i]]*classProb[[k]]
    }
    names(Pout)=names(inClassProb[[1]])
    return(Pout)
    
  }
  
  
  # ############## #
  #### Gradient ####
  # ############## #
  
  if(functionality==("gradient")){
    # Initial checks
    if(any(sapply(inClassProb, function(p) is.null(p$like) || is.null(p$grad)))) stop("INTERNAL ISSUE - Some components in inClassProb are missing the 'like' and/or 'grad' elements")
    if(is.null(classProb$like) || is.null(classProb$grad)) stop("INTERNAL ISSUE - Some components in classProb are missing the 'like' and/or 'grad' elements")
    if(apollo_inputs$apollo_control$workInLogs && apollo_inputs$apollo_control$analyticGrad) stop("INCORRECT FUNCTION/SETTING USE - workInLogs cannot be used in conjunction with analyticGrad")
    K <- length(inClassProb[[1]]$grad) # number of parameters
    if(length(classProb$grad)!=K) stop("INTERNAL ISSUE - Dimensions of gradients not consistent between classProb and inClassProb")
    if(any(sapply(inClassProb, function(p) length(p$grad))!=K)) stop("INTERNAL ISSUE - Dimensions of gradients from different components in inClassProb imply different number of parameters")
    if(length(classProb$grad)!=K) stop("INTERNAL ISSUE - Dimensions of gradients from different components in classProb imply different number of parameters")
    S <- length(classProb$like) # number of classes
    apollo_beta <- tryCatch(get("apollo_beta", envir=parent.frame(), inherits=TRUE),
                            error=function(e) stop("INTERNAL ISSUE - Gradient ",
                                                   "could not fetch apollo_beta"))
    parName <- names(apollo_beta)[!(names(apollo_beta) %in% apollo_inputs$apollo_fixed)]
    
    # Match the number of rows in each list
    classProb <- resize(classProb, nRowsInClassProb)
    inClassProb <- resize(inClassProb, nRowsInClassProb)
    # # The dimension of classProb is changed if necessary, dim(inClassProb) remains the same
    # if( nRowsClassProb!=1 && nRowsInClassProb < nRowsClassProb){
    #   classProb$like <- lapply(classProb$like, \(l) rowsum(l, group=indivID, reorder=FALSE)/nObsPerIndiv)
    #   for(k in 1:K) classProb$grad[[k]] <- lapply(classProb$grad[[k]], \(g) rowsum(g, group=indivID, reorder=FALSE)/nObsPerIndiv)
    # }
    # if( nRowsClassProb!=1 && nRowsInClassProb > nRowsClassProb){
    #   classProb$like <- lapply(classProb$like, rep, times=nObsPerIndiv)
    #   for(k in 1:K) classProb$grad[[k]] <- lapply(classProb$grad[[k]], rep, nObsPerIndiv)
    # }
    
    ### Calculate actual likelihood and gradient
    nClass <- length(inClassProb)
    Pout = list(like=0, grad=setNames(vector(mode="list", length=K), parName))
    for(k in 1:K) Pout$grad[[k]] <- 0
    for(i in 1:nClass){
      Pout$like <- Pout$like + inClassProb[[i]]$like*classProb$like[[i]]
      for(k in 1:K) Pout$grad[[k]] <- Pout$grad[[k]] + 
          inClassProb[[i]]$grad[[k]]*classProb$like[[i]] + inClassProb[[i]]$like*classProb$grad[[k]][[i]]
      if(is.array(Pout$grad[[k]])) rownames(Pout$grad[[k]]) <- NULL
      if(is.vector(Pout$grad[[k]])) names(Pout$grad[[k]]) <- NULL
    }
    
    return(Pout)
  }
  
  # ############# #
  #### Hessian ####
  # ############# #
  
  if(functionality==("hessian")){
    ### TO DO
     # Checks
     # No need to check for rows
    K <- length(inClassProb[[1]]$grad) # number of parameters
    S <- length(classProb$like) # number of classes
    
    # Match the number of rows in each list
    classProb <- resize(classProb, nRowsInClassProb)
    # # The dimension of classProb is changed if necessary, dim(inClassProb) remains the same
    # if( nRowsClassProb!=1 && nRowsInClassProb < nRowsClassProb){
    #   classProb$like <- lapply(classProb$like, \(l) rowsum(l, group=indivID, reorder=FALSE)/nObsPerIndiv)
    #   for(k in 1:K){
    #     classProb$grad[[k]] <- lapply(classProb$grad[[k]], \(g) rowsum(g, group=indivID, reorder=FALSE)/nObsPerIndiv)
    #     for(k2 in 1:K) classProb$hess[[k]][[k2]] <- lapply(classProb$hess[[k]][[k2]], \(h) rowsum(h, group=indivID, reorder=FALSE)/nObsPerIndiv)
    #   }; rm(k, k2)
    # }
    # if( nRowsClassProb!=1 && nRowsInClassProb > nRowsClassProb){
    #   classProb$like <- lapply(classProb$like, rep, times=nObsPerIndiv)
    #   for(k in 1:K){
    #     classProb$grad[[k]] <- lapply(classProb$grad[[k]], rep, times=nObsPerIndiv)
    #     for(k2 in 1:K) classProb$hess[[k]][[k2]] <- lapply(classProb$hess[[k]][[k2]], rep, times=nObsPerIndiv)
    #   }; rm(k, k2)
    # }
    
    # Calculate likelihood and gradient
    nClass <- length(classProb$like)
    apollo_beta <- tryCatch(get("apollo_beta", envir=parent.frame(), inherits=TRUE),
                            error=function(e) stop("INTERNAL ISSUE - Gradient ",
                                                   "could not fetch apollo_beta"))
    parName <- names(apollo_beta)[!(names(apollo_beta) %in% apollo_inputs$apollo_fixed)]
    L <- 0
    G <- setNames(vector(mode="list", length=K), parName)
    for(k in 1:K) G[[k]] <- 0
    for(s in 1:S){
      L <- L + inClassProb[[s]]$like*classProb$like[[s]]
      if(is.array(L)) rownames(L) <- NULL else names(L) <- NULL
      for(k in 1:K) G[[k]] <- G[[k]] + 
          inClassProb[[s]]$grad[[k]]*classProb$like[[s]] + inClassProb[[s]]$like*classProb$grad[[k]][[s]]
      if(is.array(G[[k]])) rownames(G[[k]]) <- NULL else names(G[[k]]) <- NULL
    }
    
    # Calculate the hessian
    H <- setNames(vector(mode="list", length=K), parName)
    for(k1 in 1:K){
      H[[k1]] <- setNames(vector(mode="list", length=K), parName)
      for(k2 in 1:K){
        H[[k1]][[k2]] <- 0
        for(s in 1:S) H[[k1]][[k2]] <- H[[k1]][[k2]] + 
            classProb$hess[[k1]][[k2]][[s]]*inClassProb[[s]]$like + 
            classProb$grad[[k1]][[s]]*inClassProb[[s]]$grad[[k2]] + 
            classProb$grad[[k2]][[s]]*inClassProb[[s]]$grad[[k1]] + 
            classProb$like[[s]]*inClassProb[[s]]$hess[[k1]][[k2]]
        if(is.array(H[[k1]][[k2]])) rownames(H[[k1]][[k2]]) <- NULL else 
          names(H[[k1]][[k2]]) <- NULL
      }
    }
    
    return(list(like=L, grad=G, hess=H))
  }
  
  # ############ #
  #### Report ####
  # ############ #
  if(functionality=='report'){
    P <- list()
    apollo_inputs$silent <- FALSE
    P$data  <- capture.output(lc_settings$lc_diagnostics(lc_settings, apollo_inputs, param=FALSE))
    P$param <- capture.output(lc_settings$lc_diagnostics(lc_settings, apollo_inputs, data =FALSE))
    return(P)
  }
  
}