R/helper_functions.R

In MDgof: Various Methods for the Goodness-of-Fit Problem in D>1 Dimensions

# Internal data/argument validators -----------------------------------------

is_discrete_data <- function(x) {
  out <- tryCatch(prepare_discrete_data(x), error=function(e) NULL)
  !is.null(out)
}

prepare_discrete_data <- function(x) {
  if(is.data.frame(x)) x <- as.matrix(x)
  if(!is.matrix(x) || !is.numeric(x) || ncol(x) != 3L || nrow(x) < 1L)
    return(NULL)
  if(any(!is.finite(x))) return(NULL)

  # With standard names, honor them regardless of the physical column order.
  standard_names <- c("x", "y", "counts")
  if(!is.null(colnames(x)) && all(standard_names %in% colnames(x))) {
    out <- x[, standard_names, drop=FALSE]
    counts <- out[, 3L]
    if(any(counts < 0) ||
       any(abs(counts - round(counts)) > sqrt(.Machine$double.eps))) return(NULL)

    # MDgof's discrete representation contains one row for every x-y bin pair.
    if(nrow(out) != length(unique(out[, 1L])) * length(unique(out[, 2L])))
      return(NULL)
    if(anyDuplicated(out[, 1:2, drop=FALSE])) return(NULL)
    return(out)
  }

  # Otherwise guess the count column. Counts typically have the largest number
  # of distinct values, so examine columns in that order. A candidate is only
  # accepted if the remaining two columns form the complete rectangular grid
  # required by MDgof's discrete representation.
  n_unique <- vapply(seq_len(ncol(x)), function(j) length(unique(x[, j])), integer(1))
  candidates <- order(n_unique, decreasing=TRUE)

  for(count_col in candidates) {
    counts <- x[, count_col]
    if(any(counts < 0) ||
       any(abs(counts - round(counts)) > sqrt(.Machine$double.eps))) next

    support_cols <- setdiff(seq_len(3L), count_col)
    support <- x[, support_cols, drop=FALSE]
    if(nrow(support) != length(unique(support[, 1L])) * length(unique(support[, 2L])))
      next
    if(anyDuplicated(support)) next

    out <- cbind(support, counts)
    colnames(out) <- c("x", "y", "counts")
    return(out)
  }

  NULL
}

prepare_gof_input <- function(x) {
  if(is.data.frame(x)) x <- as.matrix(x)
  discrete <- prepare_discrete_data(x)
  if(!is.null(discrete)) return(list(x=discrete, Continuous=FALSE))

  if(!is.matrix(x) || !is.numeric(x) || nrow(x) < 1L || ncol(x) < 1L)
    stop("x must be a non-empty numeric matrix (or data frame coercible to one).", call.=FALSE)
  if(any(!is.finite(x)))
    stop("x must contain only finite numeric values.", call.=FALSE)

  list(x=x, Continuous=TRUE)
}

standardize_discrete_generator <- function(fun) {
  n <- length(formals(fun))
  if(n == 0L) {
    return(function() {
      out <- prepare_discrete_data(fun())
      if(is.null(out)) stop("A discrete-data generator returned data that are not in a valid MDgof discrete format.", call.=FALSE)
      out
    })
  }
  if(n == 1L) {
    return(function(p) {
      out <- prepare_discrete_data(fun(p))
      if(is.null(out)) stop("A discrete-data generator returned data that are not in a valid MDgof discrete format.", call.=FALSE)
      out
    })
  }
  fun
}

validate_integer_scalar <- function(x, name, min=0L) {
  if(length(x) != 1L || is.na(x) || !is.finite(x) || x < min ||
     abs(x - round(x)) > sqrt(.Machine$double.eps)) {
    stop(name, " must be a single integer >= ", min, ".", call.=FALSE)
  }
  invisible(as.integer(round(x)))
}

validate_probability <- function(x, name="alpha") {
  if(length(x) != 1L || is.na(x) || !is.finite(x) || x <= 0 || x >= 1) {
    stop(name, " must be a single number strictly between 0 and 1.", call.=FALSE)
  }
  invisible(x)
}

validate_ts_output <- function(x) {
  if(!is.numeric(x) || length(x) < 1L) stop("TS must return a non-empty numeric vector.", call.=FALSE)
  nm <- names(x)
  if(is.null(nm) || anyNA(nm) || any(!nzchar(nm))) stop("TS must return a named numeric vector.", call.=FALSE)
  if(anyDuplicated(nm)) stop("TS must return a vector with unique names.", call.=FALSE)
  invisible(x)
}

#' determine number of processors
#' @param maxProcessor use this if given
#' @param dta data set used to estimate computation time
#' @param TS test statistic
#' @param typeTS format of TS
#' @param TSextra optional
#' @param B number of simulation runs
#' @param useSingleProcessor =FALSE, no parallel processing
#' @return a number
#' @keywords internal
makemaxProcessor=function(maxProcessor, dta, TS,
                          typeTS, TSextra, B,
                          useSingleProcessor=FALSE) {
  if(useSingleProcessor) return(1L)
  if(!missing(maxProcessor)) {
    validate_integer_scalar(maxProcessor, "maxProcessor", 1L)
    return(as.integer(maxProcessor))
  }
  m <- parallel::detectCores(logical=FALSE)
  if(is.na(m) || m < 1L) m <- 1L
  maxProcessor <- max(1L, as.integer(m) - 1L)
  if(maxProcessor > 1L) {
    tm <- timecheck(dta, TS, typeTS, TSextra)
    if(is.finite(tm) && tm * B < 20) {
      maxProcessor <- 1L
      message("maxProcessor set to 1 because parallel overhead would dominate this short computation")
    } else {
      message("Using ", maxProcessor, " cores.")
    }
  }
  maxProcessor
}

#' Create list with needed info
#' 
#' This function creates a list with info  needed in various parts of the package
#' 
#' @param  x data set
#' @param  Continuous =TRUE, is data continuous?
#' @param  pnull  cdf under the null hypothesis
#' @param  rnull  routine to generate data under the null hypothesis
#' @param  phat  =function(x) -99, function to estimate parameters from the data, or -99 if no parameters are estimated
#' @param  dnull =function(x) -99, density function under the null hypothesis, if available, or -99 if missing
#' @param  Ranges Range of variables
#' @param  TSextra (optional) list passed to TS, if needed.
#' @return A list containing the null-model functions and auxiliary information used by test statistics.
#' @keywords internal
makeTSextra <- function(x, Continuous, pnull, rnull, phat=function(x) -99, 
                        dnull=function(x) -99, Ranges, TSextra) {
  if(missing(TSextra)) {
    TSextra=list(pnull=pnull, rnull=rnull, phat=phat, Range=Ranges,
                 dnull=dnull, ripleyK=RipleyK, Continuous=Continuous)
    if(Continuous) {
      Eval=gen_eval(rnull, phat(x), 100) #Needed for Bakshaev_Rudzkis
      TSextra=c(TSextra, list(Eval=Eval))
    }   
  }   
  else {
    if(!is.list(TSextra)) {
      stop("TSextra must be a list when supplied.", call.=FALSE)
    } 
    TSextra = c(TSextra, pnull=pnull, rnull=rnull, phat=phat, 
                dnull=dnull, ripleyK=RipleyK, 
                Range=Ranges, Continuous=Continuous)
    if(Continuous) {
      Eval=gen_eval(rnull, phat(x), 100) #Needed for Bakshaev_Rudzkis
      TSextra=c(TSextra, list(Eval=Eval))
    }
  }  
  NoDensity=FALSE
  if(length(formals(dnull))==1 && dnull(x)[1]<0) NoDensity=TRUE
  if(ncol(x)>2) NoDensity=TRUE
  TSextra=c(TSextra, NoDensity=NoDensity)
  TSextra$knn=function(x) -99
  if(!NoDensity) {
    if(Continuous) {
      TSextra$knn=function(x) FNN::get.knn(x, 1)
      TSextra$nn=FNN::get.knn(x, 1)
      TSextra$scf=function(x, p=0) {
        if(length(formals(dnull))==1) f=function(x) dnull(x)
        if(length(formals(dnull))==2) f=function(x) dnull(x, p)
        if(!is.matrix(x)) return(grad_vec(x, f))
        grad_mat(x, f)
      }
    }
  }
  return(TSextra)
  
}

#' maketypeTS
#' 
#' find typeTS and TS
#' 
#' @param TS user-supplied test statistic, if not missing
#' @param Continuous is data continuous?
#' @return a list
#' @keywords internal
maketypeTS <- function(TS, Continuous) { 
  useSingleProcessor <- FALSE
  
  # If no custom TS was supplied, use package default test statistic
  if(missing(TS)) {
    NewTS <- FALSE
    typeTS <- 2L
    
    if(Continuous)
      TS <- TS_cont
    else
      TS <- TS_disc
    
  } else {
    if(!is.function(TS))
      stop("TS must be a function.", call. = FALSE)
    
    # Custom TS written through .Call may not be safe for parallel execution.
    if(any(grepl("^\\s*\\.Call", deparse(body(TS))))) {
      message(
        "Parallel processing is not possible if custom TS is written in C++. ",
        "Switching to single processor"
      )
      useSingleProcessor <- TRUE
    }
    
    NewTS <- TRUE
    nf <- length(formals(TS))
    
    if(Continuous) {
      
      # Continuous custom statistics retain the established interface:
      # TS(dta, pnull, param) or
      # TS(dta, pnull, param, TSextra).
      if(!nf %in% c(3L, 4L))
        stop(
          paste0(
            "For continuous data, TS must have 3 or 4 arguments: ",
            "dta, pnull, param, and optional TSextra."
          ),
          call. = FALSE
        )
      
      typeTS <- if(nf == 3L) 1L else 2L
      
    } else if(nf == 1L) {
      
      # Simple discrete interface: TS(dta)
      TS0 <- TS
      
      TS1 <- function(dta, pnull, param) {
        TS0(dta)
      }
      
      TS <- TS1
      typeTS <- 1L
      
    } else if(nf == 2L) {
      
      # Simple discrete interface: TS(dta, TSextra)
      TS0 <- TS
      
      TS2 <- function(dta, pnull, param, TSextra) {
        TS0(dta, TSextra)
      }
      
      TS <- TS2
      typeTS <- 2L
      
    } else {
      
      # Preserve the previous discrete custom-statistic signatures:
      # TS(dta, pnull, param)
      # TS(dta, pnull, param, TSextra)
      if(!nf %in% c(3L, 4L))
        stop(
          paste0(
            "For discrete data, TS must have 1 or 2 arguments ",
            "(dta and optional TSextra), or the legacy 3 or 4 arguments."
          ),
          call. = FALSE
        )
      
      typeTS <- if(nf == 3L) 1L else 2L
    }
  }
  
  list(
    TS = TS,
    typeTS = typeTS,
    NewTS = NewTS,
    useSingleProcessor = useSingleProcessor
  )
}