R/model.R

In rbiogeme: Interface to 'Biogeme'

#' Create a generic Biogeme model
#'
#' The formula is a neutral R expression. It is compiled once into a native
#' Biogeme expression graph before any estimation or simulation starts.
#'
#' @param database A `biogeme_database` object.
#' @param formula A log-likelihood expression, or a named list of formulas.
#'   The names `log_like` and `loglike` identify the likelihood formula.
#' @param weight Optional observation-weight expression.
#' @param probability Optional probability expression retained for simulation
#'   workflows.
#' @param simulations Optional named list of expressions to simulate.
#' @param panel_trajectory If `TRUE`, aggregate the likelihood using native
#'   `PanelLikelihoodTrajectory`.
#' @param draws Optional draw metadata object or list of draw metadata.
#' @param subset Optional logical expression selecting observations to retain.
#' @param availability Optional named list of availability expressions. This
#'   metadata is used by native post-estimation operations such as the null
#'   log-likelihood calculation for generic catalog models.
#' @param parameter_overrides Optional named list of native parameter
#'   replacements, keyed by the original Beta names.
#' @param control Optional [biogeme_control()] object stored with the model.
#' @return An object of class `biogeme_model`.
#' @details
#' `formula` is the native log-likelihood expression. A named list may contain
#' `log_like` (or `loglike`) and additional expressions, but a generic model
#' must still provide a likelihood, probability, or simulation expressions.
#' `simulations` is a named list evaluated only when [simulate()] is called.
#' @examples
#' database <- biogeme_database("demo", data.frame(choice = c(1, 2), x = c(1, 2)))
#' probability <- logit_probability(
#'   utilities = list(`1` = 0, `2` = biogeme_beta("b") * variable("x")),
#'   alternative = variable("choice")
#' )
#' model <- biogeme_model(
#'   database,
#'   formula = logzero(probability),
#'   simulations = list(probability = probability)
#' )
#' model
#' @export
biogeme_model <- function(
    database,
    formula = NULL,
    weight = NULL,
    probability = NULL,
    simulations = NULL,
    panel_trajectory = FALSE,
    draws = NULL,
    subset = NULL,
    parameter_overrides = NULL,
    control = NULL,
    availability = NULL
) {
  validate_biogeme_database(database)
  if (!is.null(formula) && !is_biogeme_expression(formula) && !is.list(formula)) {
    stop("formula must be a Biogeme expression or a named list of expressions.", call. = FALSE)
  }
  if (is.list(formula) && !is_biogeme_expression(formula)) {
    if (is.null(names(formula)) || anyNA(names(formula)) || any(!nzchar(names(formula)))) {
      stop("A formula list must be named.", call. = FALSE)
    }
    formula <- lapply(formula, as_biogeme_expression)
    if (!any(names(formula) %in% c("log_like", "loglike")) && is.null(simulations)) {
      stop("A formula list must contain a 'log_like' or 'loglike' expression.", call. = FALSE)
    }
  } else if (!is.null(formula)) {
    formula <- as_biogeme_expression(formula)
  }
  if (is.null(formula) && is.null(simulations) && is.null(probability)) {
    stop("Provide a formula, probability, or simulations.", call. = FALSE)
  }
  if (!is.null(weight)) weight <- as_biogeme_expression(weight)
  if (!is.null(probability)) probability <- as_biogeme_expression(probability)
  if (!is.null(availability)) {
    if (!is.list(availability) || is.null(names(availability)) ||
        anyNA(names(availability)) || any(!nzchar(names(availability))) ||
        anyDuplicated(names(availability)) || length(availability) == 0L) {
      stop("availability must be a non-empty named list of expressions.", call. = FALSE)
    }
    availability <- lapply(availability, as_biogeme_expression)
  }
  if (!is.null(simulations)) {
    if (!is.list(simulations) || is.null(names(simulations)) ||
        anyNA(names(simulations)) || any(!nzchar(names(simulations))) ||
        anyDuplicated(names(simulations)) || length(simulations) == 0L) {
      stop("simulations must be a non-empty named list of expressions.", call. = FALSE)
    }
    simulations <- lapply(simulations, as_biogeme_expression)
  }
  if (!is.logical(panel_trajectory) || length(panel_trajectory) != 1L || is.na(panel_trajectory)) {
    stop("panel_trajectory must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.null(subset)) subset <- as_biogeme_expression(subset)
  if (!is.null(parameter_overrides)) {
    if (!is.list(parameter_overrides) || is.null(names(parameter_overrides)) ||
        anyNA(names(parameter_overrides)) || any(!nzchar(names(parameter_overrides))) ||
        anyDuplicated(names(parameter_overrides))) {
      stop("parameter_overrides must be a named list of expressions.", call. = FALSE)
    }
    parameter_overrides <- lapply(parameter_overrides, as_biogeme_expression)
  }
  if (!is.null(control) && !is.list(control)) {
    stop("control must be a list returned by biogeme_control().", call. = FALSE)
  }

  all_expressions <- list()
  if (is_biogeme_expression(formula)) all_expressions <- c(all_expressions, list(formula))
  if (is.list(formula) && !is_biogeme_expression(formula)) all_expressions <- c(all_expressions, formula)
  all_expressions <- c(
    all_expressions,
    list(weight, probability, subset),
    availability,
    simulations,
    parameter_overrides
  )
  all_expressions <- Filter(Negate(is.null), all_expressions)
  validate_expression_variables(all_expressions, database, argument = "model expression")

  parameters <- list()
  base_expressions <- c(
    if (is_biogeme_expression(formula)) list(formula) else formula,
    list(weight, probability, subset),
    availability,
    simulations
  )
  base_expressions <- Filter(Negate(is.null), base_expressions)
  for (expression in base_expressions) {
    parameters <- collect_biogeme_betas(expression, parameters)
  }
  if (!is.null(parameter_overrides)) {
    for (name in names(parameter_overrides)) parameters[[name]] <- NULL
    for (expression in parameter_overrides) {
      parameters <- collect_biogeme_betas(expression, parameters)
    }
  }
  structure(
    list(
      kind = "generic",
      database = database,
      formula = formula,
      log_likelihood = if (is_biogeme_expression(formula)) formula else NULL,
      weight = weight,
      probability = probability,
      simulations = simulations,
      panel_trajectory = isTRUE(panel_trajectory),
      draws = draws,
      subset = subset,
      availability = availability,
      parameter_overrides = parameter_overrides,
      control = control,
      parameters = parameters
    ),
    class = "biogeme_model"
  )
}

#' Define estimation and simulation controls
#'
#' Unspecified fields are omitted so native Biogeme defaults remain in force.
#' Additional named arguments are passed through to the native bridge.
#' @param model_name Optional native Biogeme model name.
#' @param output_directory Optional explicit directory for native files. It is
#'   required whenever a control requests native HTML, YAML, NetCDF, pickle, or
#'   iteration output; no operation writes to the working directory by default.
#' @param seed Optional random seed.
#' @param numerically_safe Optional numerical-safety flag.
#' @param second_derivatives Optional second-derivative policy.
#' @param second_derivatives_percentage Optional native percentage of
#'   iterations using analytical second derivatives.
#' @param optimization_algorithm Optional native algorithm name.
#' @param number_of_draws Optional number of simulation draws.
#' @param draw_type Optional default draw type metadata.
#' @param draw_seed Optional draw seed (mapped to native `seed`).
#' @param bootstrap_samples Optional bootstrap sample count.
#' @param user_notes Optional notes included in native reports/results.
#' @param variance_covariance_type Optional result covariance selection, such
#'   as `"BHHH"` or `"Bootstrap"`.
#' @param save_iterations Optional iteration-file persistence flag.
#' @param generate_html Optional HTML-report flag.
#' @param generate_yaml Optional YAML-report flag.
#' @param validation_folds Optional validation-fold count.
#' @param ... Additional native Biogeme controls.
#' @return A named list of controls.
#' @details
#' Unspecified values are omitted. This lets native Biogeme apply its normal
#' defaults. The control object can be stored in a model or passed to an
#' estimation or simulation operation; it is not a Python object.
#' @examples
#' biogeme_control(
#'   seed = 1234,
#'   generate_html = FALSE,
#'   generate_yaml = FALSE,
#'   save_iterations = FALSE
#' )
#' @export
biogeme_control <- function(
    model_name = NULL,
    output_directory = NULL,
    seed = NULL,
    numerically_safe = NULL,
    second_derivatives = NULL,
    second_derivatives_percentage = NULL,
    optimization_algorithm = NULL,
    number_of_draws = NULL,
    draw_type = NULL,
    draw_seed = NULL,
    bootstrap_samples = NULL,
    user_notes = NULL,
    variance_covariance_type = NULL,
    save_iterations = NULL,
    generate_html = NULL,
    generate_yaml = NULL,
    validation_folds = NULL,
    ...
) {
  result <- list(...)
  add <- function(name, value) {
    if (!is.null(value)) result[[name]] <<- value
  }
  add("model_name", model_name)
  add("output_directory", output_directory)
  add("seed", if (is.null(draw_seed)) seed else draw_seed)
  add("numerically_safe", numerically_safe)
  add("calculating_second_derivatives", second_derivatives)
  add("second_derivatives", second_derivatives_percentage)
  add("optimization_algorithm", optimization_algorithm)
  add("number_of_draws", number_of_draws)
  add("draw_type", draw_type)
  add("bootstrap_samples", bootstrap_samples)
  add("user_notes", user_notes)
  add("variance_covariance_type", variance_covariance_type)
  add("save_iterations", save_iterations)
  add("generate_html", generate_html)
  add("generate_yaml", generate_yaml)
  add("validation_folds", validation_folds)
  if (length(result) > 0L && (is.null(names(result)) || any(!nzchar(names(result))))) {
    stop("Additional controls must be named.", call. = FALSE)
  }
  class(result) <- c("biogeme_control", "list")
  result
}

#' @export
print.biogeme_model <- function(x, ...) {
  cat("Biogeme generic model with ", length(x$parameters), " parameter definition(s)\n", sep = "")
  if (!is.null(x$database$panel_id)) cat("Panel identifier: ", x$database$panel_id, "\n", sep = "")
  invisible(x)
}

#' Validate a model before estimation
#'
#' Model construction, database operations, and native Biogeme model creation
#' are checked without running an optimizer or writing estimation results.
#' The returned diagnostics are ordinary R values; the native model object is
#' created and discarded inside the bridge.
#'
#' @param model A `biogeme_model`.
#' @param control Optional [biogeme_control()] object. It is used only when
#'   native model construction needs a control value.
#' @return An object of class `biogeme_model_validation` containing the native
#'   validation status, database information, formula names, and parameter
#'   count.
#' @details
#' This is specification validation, not cross-validation. Use [validate()]
#' after estimation when the goal is out-of-sample fold evaluation.
#' @examples
#' \dontrun{
#' database <- biogeme_database(
#'   "demo",
#'   data.frame(choice = c(1, 2), x = c(1, 2))
#' )
#' model <- logit_model(
#'   database,
#'   choice = "choice",
#'   utilities = list(`1` = 0, `2` = biogeme_beta("b") * variable("x"))
#' )
#' validate_model(model)
#' }
#' @export
validate_model <- function(model, control = NULL) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "Biogeme model validation",
      suggestion = "construct the model with logit_model() or biogeme_model()"
    )
  }
  if (!is.null(control)) {
    if (!is.list(control)) stop("control must be a named list.", call. = FALSE)
    controls <- control
  } else {
    controls <- model$control %||% list()
  }
  controls <- validate_estimation_controls(controls)
  model_name <- if (!is.null(controls$model_name)) controls$model_name else {
    "rbiogeme_validation"
  }
  model_name <- validate_name(model_name, "model_name")
  result <- biogeme_validate_model(
    model = model,
    model_name = model_name,
    controls = controls
  )
  structure(result, class = c("biogeme_model_validation", "list"))
}

#' @export
print.biogeme_model_validation <- function(x, ...) {
  cat("Biogeme model validation: ", if (isTRUE(x$valid)) "valid" else "invalid", "\n", sep = "")
  cat("Database: ", x$database_name, " (", x$number_of_rows, " rows", sep = "")
  if (isTRUE(x$is_panel)) cat(", panel")
  cat(")\n", sep = "")
  cat("Parameters: ", x$number_of_parameters, "\n", sep = "")
  cat("Formulas: ", paste(x$formula_names, collapse = ", "), "\n", sep = "")
  invisible(x)
}