test_that("Simple mnl example 1", {
  #' Simple MNL log-likelihood function assuming linear utility
  #' @param b Numeric K-long vector of parameters to be estimated.
  #' @param Y Numeric NxJ matrix. \code{Y[n,j]=1} if alt j selected in obs n.
  #' @param X List of J NxK matrices with explanatory variables.
  #' @return an N-vector of choice probabilities.

  mnl <- function(b, Y, X){
    ### Check input
    test <- is.vector(b) && is.matrix(Y) && is.list(X) && all(sapply(X, is.matrix))
    if(!test) stop("Arguments 'b', 'Y', 'X' must be a vector, matrix, ",
                   "and list of matrices, respectively.")
    N <- nrow(Y)
    K <- length(b)
    J <- length(X)
    test <- all(dim(Y)==c(N,J)) && all(sapply(X, nrow)==N) && all(sapply(X, ncol)==K)
    if(!test) stop("Dimensions of arguments 'b', 'Y', 'X' do not match.")

    ### Calculate and return MNL loglikelihood
    eV <- sapply(X, function(x) exp(x%*%b))
    p  <- rowSums(Y*eV)/rowSums(eV)
    return( p )
  }

  ### Generate synthetic data
  N <- 1000
  K <- 3
  J <- 4
  set.seed(27)
  X <- list()
  for(j in 1:J) X[[j]] <- matrix(runif(N*K), nrow=N, ncol=K)
  b<- round(runif(K, min=-1, max=1), 2)
  U <- sapply(X, function(x) x%*%b + -log(-log(runif(N))))
  Y <- U==apply(U, MARGIN=1, FUN=max)

  ### Create starting values for estimation
  # b0<- setNames(rep(0, K), paste0("b", 1:K))
  b0 <- rep(0, K)

  ### Estimate using bgw
  mnl2 <- function(b) mnl(b, Y, X) # necessary so Y and X do not need to be given

  # First test example.  CRAN reviewers prefer non-verbose output, so
  # we impose this on the initial example.
  # We do not prefer this as the default...
  test_settings = list()
  test_settings[["modelName"]] <- "Simple mnl example 1 - non-verbose"
  test_settings[["silent"]] <- TRUE
  model <- bgw_mle(calcR=mnl2, betaStart=b0, bgw_settings=test_settings)
  #
  # Print information indicating what happened.
  cat(model$bgw_settings[["modelName"]])
  cat("\n")
  cat("\nOutcome of estimation search:",model$message,"\n")
  cat("Number of iterations: ",model$iterations,"\n")
  cat("\nMLE parameter estimate:\n")
  cat(model$estimate)
  cat("\n")
  cat("\nEstimated t-ratios:\n")
  cat(model$tstatBGW)
  cat("\n")

  expect_equal(model$maximum, -1361.97, tolerance = 0.01)
})

test_that("Simple mnl model 1 - non-verbose", {
  #' Simple MNL loglikelihood function assuming linear utility
  #' @param b Numeric K-long vector of parameters to be estimated.
  #' @param Y Numeric NxJ matrix. \code{Y[n,j]=1} if alt j selected in obs n.
  #' @param X List of J NxK matrices with explanatory variables.
  #' @return an N-vector of choice probabilities.

  mnl <- function(b, Y, X){
    ### Check input
    test <- is.vector(b) && is.matrix(Y) && is.list(X) && all(sapply(X, is.matrix))
    if(!test) stop("Arguments 'b', 'Y', 'X' must be a vector, matrix, ",
                   "and list of matrices, respectively.")
    N <- nrow(Y)
    K <- length(b)
    J <- length(X)
    test <- all(dim(Y)==c(N,J)) && all(sapply(X, nrow)==N) && all(sapply(X, ncol)==K)
    if(!test) stop("Dimensions of arguments 'b', 'Y', 'X' do not match.")

    ### Calculate and return MNL loglikelihood
    eV <- sapply(X, function(x) exp(x%*%b))
    p  <- rowSums(Y*eV)/rowSums(eV)
    return( p )
  }

  ### Generate synthetic data
  N <- 1000
  K <- 3
  J <- 4
  set.seed(27)
  X <- list()
  for(j in 1:J) X[[j]] <- matrix(runif(N*K), nrow=N, ncol=K)
  b<- round(runif(K, min=-1, max=1), 2)
  U <- sapply(X, function(x) x%*%b + -log(-log(runif(N))))
  Y <- U==apply(U, MARGIN=1, FUN=max)

  ### Create starting values for estimation
  b0<- setNames(rep(0, K), paste0("b", 1:K))

  ### Estimate using bgw
  mnl2 <- function(b) mnl(b, Y, X) # necessary so Y and X do not need to be given

  test_settings = list()
  test_settings[["printLevel"]] <- 2L
  model <- bgw_mle(calcR=mnl2, betaStart=b0, bgw_settings = test_settings)

  expect_equal(model$maximum, -1361.97, tolerance = 0.01)
})

test_that("Using settings to turn off VC calculation for mnl", {
  #' Simple MNL loglikelihood function assuming linear utility
  #' @param b Numeric K-long vector of parameters to be estimated.
  #' @param Y Numeric NxJ matrix. \code{Y[n,j]=1} if alt j selected in obs n.
  #' @param X List of J NxK matrices with explanatory variables.

  mnl <- function(b, Y, X){
    ### Check input
    test <- is.vector(b) && is.matrix(Y) && is.list(X) && all(sapply(X, is.matrix))
    if(!test) stop("Arguments 'b', 'Y', 'X' must be a vector, matrix, ",
                   "and list of matrices, respectively.")
    N <- nrow(Y)
    K <- length(b)
    J <- length(X)
    test <- all(dim(Y)==c(N,J)) && all(sapply(X, nrow)==N) && all(sapply(X, ncol)==K)
    if(!test) stop("Dimensions of arguments 'b', 'Y', 'X' do not match.")

    ### Calculate and return MNL loglikelihood
    eV <- sapply(X, function(x) exp(x%*%b))
    p  <- rowSums(Y*eV)/rowSums(eV)
    return( p )
  }

  ### Generate synthetic data
  N <- 1000
  K <- 3
  J <- 4
  set.seed(27)
  X <- list()
  for(j in 1:J) X[[j]] <- matrix(runif(N*K), nrow=N, ncol=K)
  b<- round(runif(K, min=-1, max=1), 2)
  U <- sapply(X, function(x) x%*%b + -log(-log(runif(N))))
  Y <- U==apply(U, MARGIN=1, FUN=max)

  ### Create starting values for estimation
  b0<- setNames(rep(0, K), paste0("b", 1:K))

  ### Estimate using bgw
  mnl2 <- function(b) mnl(b, Y, X) # necessary so Y and X do not need to be given

  test_settings = list()
  test_settings[["printLevel"]] <- 2L
  test_settings[["vcHessianMethod"]] <- "none"
  model <- bgw_mle(calcR=mnl2, betaStart=b0, bgw_settings = test_settings)

  expect_equal(model$maximum, -1361.97, tolerance = 0.01)
})

# ---- The R / Fortran interface ------------------------------------------------

set.seed(1)
int_n  <- 300
int_z1 <- rnorm(int_n)
int_z2 <- rnorm(int_n)
int_y  <- rbinom(int_n, 1, plogis(0.8*int_z1 - 0.5*int_z2))
int_R  <- function(b) {
  b <- unname(b)
  p <- plogis(b[1]*int_z1 + b[2]*int_z2)
  ifelse(int_y == 1, p, 1 - p)
}
int_J  <- function(b) {
  b <- unname(b)
  p <- plogis(b[1]*int_z1 + b[2]*int_z2)
  cbind(int_z1*(int_y - p), int_z2*(int_y - p))
}

test_that("the Fortran does not overwrite an unnamed betaStart when calcJ is supplied", {
  b <- c(0.1, 0.2)
  m <- bgw_mle(int_R, b, calcJ = int_J, bgw_settings = list(silent = TRUE),
               calcJ_type = "log_probability")
  expect_identical(b, c(0.1, 0.2))
  expect_identical(m$betaStart, c(0.1, 0.2))
  expect_false(isTRUE(all.equal(unname(m$estimate), c(0.1, 0.2))))
})

test_that("each beta given to calcR is left unchanged by later iterations", {
  seen <- list()
  remember <- function(b) {
    seen[[length(seen) + 1]] <<- b
    int_R(b)
  }
  m <- bgw_mle(remember, c(0.1, 0.2), calcJ = int_J, bgw_settings = list(silent = TRUE),
               calcJ_type = "log_probability")
  expect_identical(seen[[1]], c(0.1, 0.2))
  expect_gt(length(unique(seen)), 1)
})

test_that("limits too large for a Fortran integer are rejected", {
  for (limit in c(1e10, Inf)) {
    expect_error(bgw_mle(int_R, c(a = 0.1, b = 0.2),
                         bgw_settings = list(silent = TRUE, maxIterations = limit)),
                 "not valid")
    expect_error(bgw_mle(int_R, c(a = 0.1, b = 0.2),
                         bgw_settings = list(silent = TRUE, maxFunctionEvals = limit)),
                 "not valid")
  }
})

test_that("a model too large for a Fortran integer is rejected with a clear message", {
  expect_error(bgw:::bgw_mle_setup(2, 1e10, TRUE), "too large for BGW")
})