R/expressions.R

In rbiogeme: Interface to 'Biogeme'

#' Create a Biogeme expression node
#'
#' This is an internal constructor. Expressions are represented in R until a
#' complete model is compiled by the Python bridge.
#'
#' @param kind Expression-node kind.
#' @param ... Node attributes.
#' @return An object of class `biogeme_expression`.
#' @keywords internal
new_biogeme_expression <- function(kind, ...) {
  structure(
    c(list(kind = kind), list(...)),
    class = "biogeme_expression"
  )
}

is_biogeme_expression <- function(x) {
  inherits(x, "biogeme_expression")
}

as_biogeme_expression <- function(x) {
  if (is_biogeme_expression(x)) {
    return(x)
  }
  if (inherits(x, "biogeme_draws")) {
    return(draw(x$name, x$draw_type))
  }
  if (is.numeric(x) && length(x) == 1L && !is.na(x)) {
    return(new_biogeme_expression("numeric", value = as.numeric(x)))
  }
  stop(
    "Expressions can only contain scalar numeric values or Biogeme expressions.",
    call. = FALSE
  )
}

validate_name <- function(name, argument) {
  if (!is.character(name) || length(name) != 1L || is.na(name) || !nzchar(name)) {
    stop(argument, " must be a single non-empty character string.", call. = FALSE)
  }
  name
}

validate_optional_bound <- function(value, argument) {
  if (is.null(value)) {
    return(NULL)
  }
  if (!is.numeric(value) || length(value) != 1L || is.na(value)) {
    stop(argument, " must be NULL or one numeric value.", call. = FALSE)
  }
  as.numeric(value)
}

#' Create a Biogeme parameter
#'
#' @param name Parameter name.
#' @param start Initial value.
#' @param lower Optional lower bound.
#' @param upper Optional upper bound.
#' @param fixed Logical value indicating whether the parameter is fixed.
#' @param sigma_prior Standard deviation of the native default normal prior.
#' @param prior Optional declarative prior from [biogeme_prior()]. It is
#'   compiled into a native PyMC prior factory; R callbacks are not accepted.
#' @return A parameter expression.
#' @details
#' Parameter names are preserved exactly by the bridge and are therefore part
#' of the equivalence contract with native Biogeme. Bounds and `fixed` are
#' passed to the native parameter definition.
#' @examples
#' b_cost <- biogeme_beta("b_cost", start = -1, upper = 0)
#' b_cost
#' @export
biogeme_beta <- function(
    name,
    start = 0,
    lower = NULL,
    upper = NULL,
    fixed = FALSE,
    sigma_prior = 5,
    prior = NULL
) {
  name <- validate_name(name, "name")
  if (!is.numeric(start) || length(start) != 1L || is.na(start) || !is.finite(start)) {
    stop("start must be one finite numeric value.", call. = FALSE)
  }
  lower <- validate_optional_bound(lower, "lower")
  upper <- validate_optional_bound(upper, "upper")
  if (!is.null(lower) && !is.null(upper) && lower > upper) {
    stop("lower must not be greater than upper.", call. = FALSE)
  }
  if (!is.null(lower) && start < lower) {
    stop("start must not be smaller than lower.", call. = FALSE)
  }
  if (!is.null(upper) && start > upper) {
    stop("start must not be greater than upper.", call. = FALSE)
  }
  if (!is.logical(fixed) || length(fixed) != 1L || is.na(fixed)) {
    stop("fixed must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.numeric(sigma_prior) || length(sigma_prior) != 1L ||
      is.na(sigma_prior) || !is.finite(sigma_prior) || sigma_prior <= 0) {
    stop("sigma_prior must be one positive finite numeric value.", call. = FALSE)
  }
  if (!is.null(prior) && !inherits(prior, "biogeme_prior")) {
    stop("prior must be NULL or a biogeme_prior() object.", call. = FALSE)
  }

  structure(
    new_biogeme_expression(
      "beta",
      name = name,
      start = as.numeric(start),
      lower = lower,
      upper = upper,
      fixed = isTRUE(fixed),
      sigma_prior = as.numeric(sigma_prior),
      prior = prior
    ),
    class = c("biogeme_beta", "biogeme_expression")
  )
}

#' Declare a native Bayesian prior
#'
#' This descriptor is intentionally data-only. The Python bridge maps it to a
#' PyMC distribution before native Biogeme starts estimation, so no R
#' function or callback is evaluated by the sampler.
#'
#' @param distribution Prior distribution. Currently `"normal"` and
#'   `"student_t"` are supported.
#' @param sigma Positive scale parameter.
#' @param nu Positive degrees of freedom for a Student-t prior.
#' @return A declarative `biogeme_prior` object.
#' @export
biogeme_prior <- function(
    distribution = c("normal", "student_t"),
    sigma = 5,
    nu = 5
) {
  distribution <- match.arg(distribution)
  if (!is.numeric(sigma) || length(sigma) != 1L || is.na(sigma) ||
      !is.finite(sigma) || sigma <= 0) {
    stop("sigma must be one positive finite numeric value.", call. = FALSE)
  }
  if (!is.numeric(nu) || length(nu) != 1L || is.na(nu) ||
      !is.finite(nu) || nu <= 0) {
    stop("nu must be one positive finite numeric value.", call. = FALSE)
  }
  structure(
    list(
      distribution = distribution,
      sigma = as.numeric(sigma),
      nu = as.numeric(nu)
    ),
    class = "biogeme_prior"
  )
}

#' Create a Biogeme data variable
#'
#' @param name Column name in the Biogeme database.
#' @return A variable expression.
#' @details
#' `variable()` creates a symbolic reference. It does not extract an R column
#' or calculate a vector immediately. Use it in arithmetic and model
#' expressions after the corresponding database column has been defined.
#' @examples
#' x <- variable("income")
#' b <- biogeme_beta("b_income")
#' b * x + 1
#' @export
variable <- function(name) {
  structure(
    new_biogeme_expression("variable", name = validate_name(name, "name")),
    class = c("biogeme_variable", "biogeme_expression")
  )
}

#' Define one coefficient-variable term for a linear utility
#'
#' @param beta A `biogeme_beta` expression.
#' @param x A `biogeme_variable` expression.
#' @return A linear-utility term.
#' @export
linear_term <- function(beta, x) {
  beta <- as_biogeme_expression(beta)
  x <- as_biogeme_expression(x)
  if (!identical(beta$kind, "beta") || !identical(x$kind, "variable")) {
    stop("A linear term must contain one Beta and one data variable.", call. = FALSE)
  }
  structure(list(beta = beta, x = x), class = "biogeme_linear_term")
}

#' @rdname linear_term
#' @export LinearTermTuple
LinearTermTuple <- linear_term

#' Construct a native linear utility expression
#'
#' @param terms A non-empty list of `linear_term()` objects.
#' @return A Biogeme expression compiled to native `LinearUtility`.
#' @export
linear_utility <- function(terms) {
  if (!is.list(terms) || length(terms) == 0L ||
      any(!vapply(terms, inherits, logical(1), what = "biogeme_linear_term"))) {
    stop("terms must be a non-empty list of linear_term() objects.", call. = FALSE)
  }
  new_biogeme_expression("linear_utility", terms = terms)
}

#' @rdname linear_utility
#' @export LinearUtility
LinearUtility <- linear_utility

#' Create a neutral controller for one or more expression catalogs
#'
#' Catalog controllers select the same named specification in every catalog
#' that shares the controller. The object remains an R specification until the
#' complete expression is compiled to native Biogeme.
#'
#' @param name Controller name.
#' @param specification_names Ordered, unique catalog specification names.
#' @param selected_name Optional specification selected before native compilation.
#' @return A neutral catalog controller.
#' @export
biogeme_catalog_controller <- function(
    name,
    specification_names,
    selected_name = NULL
) {
  name <- validate_name(name, "name")
  if (grepl("[;:]", name)) {
    stop("name must not contain ';' or ':'.", call. = FALSE)
  }
  if (!is.character(specification_names) || length(specification_names) == 0L ||
      anyNA(specification_names) || any(!nzchar(specification_names)) ||
      any(grepl("[;:]", specification_names)) ||
      anyDuplicated(specification_names)) {
    stop(
      "specification_names must be a non-empty character vector with unique names.",
      call. = FALSE
    )
  }
  if (!is.null(selected_name)) {
    if (!is.character(selected_name) || length(selected_name) != 1L ||
        is.na(selected_name) || !nzchar(selected_name) ||
        !selected_name %in% specification_names) {
      stop(
        "selected_name must be NULL or one of specification_names.",
        call. = FALSE
      )
    }
  }
  structure(
    list(
      name = name,
      specification_names = as.character(specification_names),
      selected_name = selected_name
    ),
    class = "biogeme_catalog_controller"
  )
}

#' @rdname biogeme_catalog_controller
#' @export
catalog_controller <- biogeme_catalog_controller

#' Create a native-compatible expression catalog
#'
#' @param name Catalog name.
#' @param expressions Non-empty named list of alternative expressions.
#' @param controller Optional [biogeme_catalog_controller()] shared by catalogs.
#' @return A catalog expression compiled to native Biogeme.
#' @export
biogeme_catalog <- function(name, expressions, controller = NULL) {
  name <- validate_name(name, "name")
  if (grepl("[;:]", name)) {
    stop("name must not contain ';' or ':'.", call. = FALSE)
  }
  if (!is.list(expressions) || length(expressions) == 0L ||
      is.null(names(expressions)) || anyNA(names(expressions)) ||
      any(!nzchar(names(expressions))) || any(grepl("[;:]", names(expressions))) ||
      anyDuplicated(names(expressions))) {
    stop(
      "expressions must be a non-empty named list with unique names.",
      call. = FALSE
    )
  }
  expressions <- lapply(expressions, as_biogeme_expression)
  specification_names <- names(expressions)
  if (is.null(controller)) {
    controller <- biogeme_catalog_controller(name, specification_names)
  }
  if (!inherits(controller, "biogeme_catalog_controller")) {
    stop("controller must be a biogeme_catalog_controller object.", call. = FALSE)
  }
  if (!identical(controller$specification_names, specification_names)) {
    stop(
      "controller specification_names must match the catalog expression names.",
      call. = FALSE
    )
  }
  node <- new_biogeme_expression(
    "catalog",
    name = name,
    expressions = expressions,
    controller = controller
  )
  names(node$expressions) <- specification_names
  node
}

#' @rdname biogeme_catalog
#' @export
catalog <- biogeme_catalog

#' Build synchronized generic/alternative-specific catalogs
#'
#' This is the neutral R counterpart of native Biogeme's
#' `generic_alt_specific_catalogs()`.  Each returned element is a named list
#' of catalogs, one for every alternative.  The catalogs share one native
#' controller with the exact specifications `generic` and `altspec`; when
#' segmentations are supplied, the segmentation catalogs are nested below
#' that controller in the same order as native Biogeme.
#'
#' @param generic_name Shared name used for the segmentation and
#'   generic/alternative-specific controllers.
#' @param beta_parameters Non-empty list of `biogeme_beta()` expressions.
#' @param alternatives At least two unique alternative names.
#' @param potential_segmentations Optional non-empty list of segmentation
#'   objects.
#' @param maximum_number Maximum number of potential segmentations in one
#'   option.
#' @return A list of named lists of catalog expressions, one list per Beta.
#' @export
biogeme_generic_alt_specific_catalogs <- function(
    generic_name,
    beta_parameters,
    alternatives,
    potential_segmentations = NULL,
    maximum_number = 5L
) {
  generic_name <- validate_name(generic_name, "generic_name")
  if (grepl("[;:]", generic_name)) {
    stop("generic_name must not contain ';' or ':'.", call. = FALSE)
  }
  if (!is.list(beta_parameters) || length(beta_parameters) == 0L ||
      any(!vapply(beta_parameters, function(x) identical(x$kind, "beta"), logical(1)))) {
    stop("beta_parameters must be a non-empty list of Beta expressions.", call. = FALSE)
  }
  if (!is.character(alternatives) || length(alternatives) < 2L ||
      anyNA(alternatives) || any(!nzchar(alternatives)) ||
      any(grepl("[;:]", alternatives)) || anyDuplicated(alternatives)) {
    stop(
      "alternatives must contain at least two unique non-empty names.",
      call. = FALSE
    )
  }
  if (!is.null(potential_segmentations)) {
    if (!is.list(potential_segmentations) || length(potential_segmentations) == 0L ||
        any(!vapply(
          potential_segmentations,
          inherits,
          logical(1),
          what = "biogeme_segmentation"
        ))) {
      stop(
        "potential_segmentations must be NULL or a non-empty list of segmentation objects.",
        call. = FALSE
      )
    }
  }
  if (!is.numeric(maximum_number) || length(maximum_number) != 1L ||
      is.na(maximum_number) || maximum_number < 0 ||
      maximum_number != floor(maximum_number)) {
    stop("maximum_number must be one non-negative integer.", call. = FALSE)
  }
  maximum_number <- as.integer(maximum_number)

  # Native SegmentedParameters creates the generic Beta first, followed by
  # one alternative-specific Beta block per alternative.  Keeping that order
  # is important because the generated parameter names are part of the
  # equivalence contract.
  all_beta_parameters <- beta_parameters
  for (alternative in alternatives) {
    all_beta_parameters <- c(
      all_beta_parameters,
      lapply(beta_parameters, function(beta) {
        # Native SegmentedParameters copies the value, bounds, and fixed
        # status, while newly created alternative-specific Betas use native
        # default prior metadata.
        biogeme_beta(
          name = paste0(beta$name, "_", alternative),
          start = beta$start,
          lower = beta$lower,
          upper = beta$upper,
          fixed = beta$fixed
        )
      })
    )
  }

  if (is.null(potential_segmentations)) {
    segmented_catalogs <- all_beta_parameters
  } else {
    segmented_catalogs <- biogeme_segmentation_catalogs(
      generic_name = generic_name,
      beta_parameters = all_beta_parameters,
      potential_segmentations = potential_segmentations,
      maximum_number = maximum_number
    )
  }

  generic_specific_controller <- biogeme_catalog_controller(
    paste0(generic_name, "_gen_altspec"),
    c("generic", "altspec")
  )
  number_of_parameters <- length(beta_parameters)
  result <- vector("list", number_of_parameters)
  names(result) <- vapply(beta_parameters, function(beta) beta$name, character(1))

  for (parameter_index in seq_len(number_of_parameters)) {
    generic_expression <- segmented_catalogs[[parameter_index]]
    alternative_expressions <- lapply(seq_along(alternatives), function(index) {
      segmented_catalogs[[number_of_parameters * index + parameter_index]]
    })
    catalogs_for_parameter <- lapply(seq_along(alternatives), function(index) {
      alternative <- alternatives[[index]]
      biogeme_catalog(
        name = paste0(
          beta_parameters[[parameter_index]]$name,
          "_",
          alternative,
          "_gen_altspec"
        ),
        expressions = list(
          generic = generic_expression,
          altspec = alternative_expressions[[index]]
        ),
        controller = generic_specific_controller
      )
    })
    names(catalogs_for_parameter) <- alternatives
    result[[parameter_index]] <- catalogs_for_parameter
  }
  result
}

#' @rdname biogeme_generic_alt_specific_catalogs
#' @export
generic_alt_specific_catalogs <- biogeme_generic_alt_specific_catalogs

#' Build synchronized catalogs for possible parameter segmentations
#'
#' This is the structural counterpart of native Biogeme's
#' `segmentation_catalogs()`. It only enumerates catalog names and builds
#' neutral expression nodes; native Biogeme still creates the segmented
#' parameters and evaluates the selected specification.
#'
#' @param generic_name Shared controller name.
#' @param beta_parameters Non-empty list of `biogeme_beta()` expressions.
#' @param potential_segmentations Non-empty list of segmentation objects.
#' @param maximum_number Maximum number of segmentations in one option.
#' @param selected_name Optional segmentation specification selected before
#'   native compilation.
#' @return A list of synchronized catalog expressions, one per beta.
#' @export
biogeme_segmentation_catalogs <- function(
    generic_name,
    beta_parameters,
    potential_segmentations,
    maximum_number,
    selected_name = NULL
) {
  generic_name <- validate_name(generic_name, "generic_name")
  if (!is.list(beta_parameters) || length(beta_parameters) == 0L ||
      any(!vapply(beta_parameters, function(x) identical(x$kind, "beta"), logical(1)))) {
    stop("beta_parameters must be a non-empty list of Beta expressions.", call. = FALSE)
  }
  if (!is.list(potential_segmentations) || length(potential_segmentations) == 0L ||
      any(!vapply(potential_segmentations, inherits, logical(1), what = "biogeme_segmentation"))) {
    stop(
      "potential_segmentations must be a non-empty list of segmentation objects.",
      call. = FALSE
    )
  }
  if (!is.numeric(maximum_number) || length(maximum_number) != 1L ||
      is.na(maximum_number) || maximum_number < 0 ||
      maximum_number != floor(maximum_number)) {
    stop("maximum_number must be one non-negative integer.", call. = FALSE)
  }
  maximum_number <- as.integer(maximum_number)
  if (!is.null(selected_name)) {
    if (!is.character(selected_name) || length(selected_name) != 1L ||
        is.na(selected_name) || !nzchar(selected_name)) {
      stop("selected_name must be NULL or one non-empty string.", call. = FALSE)
    }
  }

  # expand.grid varies its first column fastest. Reversing the columns gives
  # the same False/True product order as itertools.product in native Python.
  n_segmentations <- length(potential_segmentations)
  combinations <- expand.grid(
    rep(list(c(FALSE, TRUE)), n_segmentations),
    KEEP.OUT.ATTRS = FALSE,
    stringsAsFactors = FALSE
  )
  combinations <- combinations[n_segmentations:1L]
  keep <- rowSums(combinations) <= maximum_number
  combinations <- combinations[keep, , drop = FALSE]
  specification_names <- vapply(seq_len(nrow(combinations)), function(row) {
    selected <- as.logical(combinations[row, ])
    if (!any(selected)) return("no_seg")
    paste(
      vapply(potential_segmentations[selected], function(segmentation) segmentation$variable, character(1)),
      collapse = "-"
    )
  }, character(1))
  controller <- biogeme_catalog_controller(
    generic_name,
    specification_names,
    selected_name = selected_name
  )

  lapply(beta_parameters, function(beta) {
    options <- lapply(seq_len(nrow(combinations)), function(row) {
      selected <- as.logical(combinations[row, ])
      if (!any(selected)) {
        beta
      } else {
        segment_beta(beta, potential_segmentations[selected])
      }
    })
    names(options) <- specification_names
    biogeme_catalog(
      name = paste0("segmented_", beta$name),
      expressions = options,
      controller = controller
    )
  })
}

#' @rdname biogeme_segmentation_catalogs
#' @export
segmentation_catalogs <- biogeme_segmentation_catalogs

#' Define a discrete parameter segmentation
#'
#' @param variable A data-variable name or expression.
#' @param mapping Named mapping from integer values to segment names.
#' @param reference Optional reference segment name.
#' @return A `biogeme_segmentation` specification.
#' @export
biogeme_segmentation <- function(variable, mapping, reference = NULL) {
  variable_expression <- if (is.character(variable)) {
    structure(
      new_biogeme_expression("variable", name = validate_name(variable[[1L]], "variable")),
      class = c("biogeme_variable", "biogeme_expression")
    )
  } else {
    as_biogeme_expression(variable)
  }
  if (!identical(variable_expression$kind, "variable")) {
    stop("variable must be a database variable name or expression.", call. = FALSE)
  }
  if (!is.atomic(mapping) && !is.list(mapping)) {
    stop("mapping must be a named vector or list.", call. = FALSE)
  }
  mapping <- unlist(mapping, use.names = TRUE)
  if (length(mapping) == 0L || is.null(names(mapping)) ||
      anyNA(names(mapping)) || any(!nzchar(names(mapping)))) {
    stop("mapping must be non-empty and named by integer values.", call. = FALSE)
  }
  keys <- suppressWarnings(as.integer(names(mapping)))
  labels <- as.character(unname(mapping))
  if (anyNA(keys) || anyDuplicated(keys) || anyNA(labels) || any(!nzchar(labels)) ||
      anyDuplicated(labels)) {
    stop("mapping keys must be unique integers and labels must be unique non-empty strings.", call. = FALSE)
  }
  if (is.null(reference)) reference <- labels[[1L]]
  reference <- validate_name(reference, "reference")
  if (!reference %in% labels) {
    stop("reference must be one of the mapping labels.", call. = FALSE)
  }
  structure(
    list(
      variable = variable_expression$name,
      mapping = stats::setNames(labels, as.character(keys)),
      reference = reference
    ),
    class = "biogeme_segmentation"
  )
}

#' Generate a native-compatible segmentation specification from a database
#'
#' @param database A `biogeme_database` object.
#' @param variable A database variable name or expression.
#' @param mapping Named mapping from integer values to segment names.
#' @param reference Optional reference segment name.
#' @return A `biogeme_segmentation` specification.
#' @export
biogeme_database_segmentation <- function(database, variable, mapping, reference = NULL) {
  validate_biogeme_database(database)
  variable_expression <- if (is.character(variable)) {
    structure(
      new_biogeme_expression("variable", name = validate_name(variable[[1L]], "variable")),
      class = c("biogeme_variable", "biogeme_expression")
    )
  } else {
    as_biogeme_expression(variable)
  }
  if (!biogeme_database_has_column(database, variable_expression$name)) {
    stop("Segmentation variable is not present in the database.", call. = FALSE)
  }
  biogeme_segmentation(variable_expression, mapping, reference)
}

#' @rdname biogeme_database_segmentation
#' @export
database_generate_segmentation <- biogeme_database_segmentation

#' Create a segmented parameter expression
#'
#' The generated native parameters follow Biogeme's public naming contract:
#' `<base>_ref` for the reference and `<base>_diff_<segment>` for differences.
#'
#' @param beta A `biogeme_beta` expression.
#' @param segmentations A non-empty list of `biogeme_segmentation` objects.
#' @param prefix Native segmentation expression prefix.
#' @return A Biogeme expression compiled through native `Segmentation`.
#' @export
segment_beta <- function(beta, segmentations, prefix = "segmented") {
  beta <- as_biogeme_expression(beta)
  if (!identical(beta$kind, "beta")) {
    stop("beta must be a biogeme_beta expression.", call. = FALSE)
  }
  if (!is.list(segmentations) || length(segmentations) == 0L ||
      any(!vapply(segmentations, inherits, logical(1), what = "biogeme_segmentation"))) {
    stop("segmentations must be a non-empty list of biogeme_segmentation objects.", call. = FALSE)
  }
  structure(
    new_biogeme_expression(
      "segmented_beta",
      beta = beta,
      segmentations = segmentations,
      prefix = validate_name(prefix, "prefix")
    ),
    class = c("biogeme_segmented_beta", "biogeme_expression")
  )
}

#' @rdname segment_beta
#' @export
segmented_beta <- segment_beta

biogeme_segmented_beta_definitions <- function(expression) {
  if (!is_biogeme_expression(expression) || !identical(expression$kind, "segmented_beta")) {
    stop("expression must be a segmented_beta expression.", call. = FALSE)
  }
  base <- expression$beta
  definitions <- list(
    biogeme_beta(
      paste0(base$name, "_ref"),
      start = base$start,
      lower = base$lower,
      upper = base$upper,
      fixed = base$fixed
    )
  )
  for (segmentation in expression$segmentations) {
    labels <- unname(segmentation$mapping)
    labels <- labels[labels != segmentation$reference]
    for (label in labels) {
      definitions[[length(definitions) + 1L]] <- biogeme_beta(
        paste0(base$name, "_diff_", label),
        start = base$start,
        fixed = base$fixed
      )
    }
  }
  definitions
}

binary_expression <- function(operator, left, right) {
  if (!operator %in% c("+", "-", "*", "/", "^", "==", "!=", "<", "<=", ">", ">=", "&", "|")) {
    stop("Unsupported Biogeme binary operator: ", operator, call. = FALSE)
  }
  new_biogeme_expression(
    "binary",
    operator = operator,
    left = as_biogeme_expression(left),
    right = as_biogeme_expression(right)
  )
}

#' @export
#' @method Ops biogeme_expression
Ops.biogeme_expression <- function(e1, e2 = NULL) {
  operator <- .Generic
  if (operator == "!") {
    return(new_biogeme_expression("unary", operator = "!", operand = e1))
  }
  if (is.null(e2)) {
    stop("A second operand is required for Biogeme operator '", operator, "'.", call. = FALSE)
  }
  binary_expression(operator, e1, e2)
}

#' @export
#' @method + biogeme_expression
`+.biogeme_expression` <- function(e1, e2) {
  binary_expression("+", e1, e2)
}

#' @export
#' @method - biogeme_expression
`-.biogeme_expression` <- function(e1, e2 = NULL) {
  if (is.null(e2)) {
    return(new_biogeme_expression("unary", operator = "-", operand = e1))
  }
  binary_expression("-", e1, e2)
}

#' @export
#' @method * biogeme_expression
`*.biogeme_expression` <- function(e1, e2) {
  binary_expression("*", e1, e2)
}

#' @export
#' @method / biogeme_expression
`/.biogeme_expression` <- function(e1, e2) {
  binary_expression("/", e1, e2)
}

#' @export
#' @method ^ biogeme_expression
`^.biogeme_expression` <- function(e1, e2) {
  binary_expression("^", e1, e2)
}

#' @export
#' @method log biogeme_expression
log.biogeme_expression <- function(x, base = exp(1)) {
  if (!is.numeric(base) || length(base) != 1L || is.na(base) || base != exp(1)) {
    stop("Only the natural logarithm is supported.", call. = FALSE)
  }
  new_biogeme_expression("unary", operator = "log", operand = as_biogeme_expression(x))
}

#' @export
#' @method exp biogeme_expression
exp.biogeme_expression <- function(x) {
  new_biogeme_expression("unary", operator = "exp", operand = as_biogeme_expression(x))
}

#' Create a Biogeme function node
#'
#' This constructor is intentionally generic so that new public Biogeme
#' functions can be added without evaluating the expression in R.
#' @param operator Native Biogeme function/operator name.
#' @param args A list of scalar values or Biogeme expressions.
#' @return A Biogeme expression.
#' @keywords internal
biogeme_function <- function(operator, args = list(), ...) {
  operator <- validate_name(operator, "operator")
  if (!is.list(args)) {
    stop("args must be a list.", call. = FALSE)
  }
  new_biogeme_expression(
    "function",
    operator = operator,
    args = lapply(args, as_biogeme_expression),
    attributes = list(...)
  )
}

#' A numerically safe logarithm that returns zero at zero
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
logzero <- function(x) {
  biogeme_function("logzero", list(x))
}

#' Construct a native ordered-response log-likelihood expression
#'
#' The observed response is matched to `categories`; the intervals are
#' defined by the ordered `cutpoints`.  This internal helper selects the
#' native logistic or normal cumulative distribution.  No probability is
#' evaluated in R.
#'
#' @param eta Latent-index expression.
#' @param cutpoints Non-empty list of cutpoint expressions.
#' @param alternative Observed response expression.
#' @param categories Optional ordered numeric labels. If omitted, native
#'   Biogeme uses `0:(length(cutpoints))`.
#' @param neutral_labels Optional numeric labels whose contribution is one.
#' @param enforce_order Whether native Biogeme should enforce ordered
#'   cutpoints during evaluation.
#' @param eps Positive numerical lower bound used by native Biogeme.
#' @param distribution Native ordered-response distribution.
#' @return A Biogeme expression compiled to native Biogeme.
#' @keywords internal
ordered_response_log_probability <- function(
    distribution = c("logit", "probit"),
    eta,
    cutpoints,
    alternative,
    categories = NULL,
    neutral_labels = numeric(),
    enforce_order = TRUE,
    eps = 1e-12
) {
  distribution <- match.arg(distribution)
  eta <- as_biogeme_expression(eta)
  alternative <- as_biogeme_expression(alternative)
  if (!is.list(cutpoints) || length(cutpoints) == 0L) {
    stop("cutpoints must be a non-empty list of expressions.", call. = FALSE)
  }
  cutpoints <- lapply(cutpoints, as_biogeme_expression)
  n_categories <- length(cutpoints) + 1L
  if (!is.null(categories)) {
    if (!is.numeric(categories) || length(categories) != n_categories ||
        anyNA(categories) || any(!is.finite(categories))) {
      stop(
        "categories must be NULL or a finite numeric vector of length ",
        n_categories,
        ".",
        call. = FALSE
      )
    }
    if (anyDuplicated(categories)) {
      stop("categories must contain unique labels.", call. = FALSE)
    }
    categories <- as.numeric(categories)
  }
  if (!is.numeric(neutral_labels) || anyNA(neutral_labels) ||
      any(!is.finite(neutral_labels))) {
    stop("neutral_labels must be a finite numeric vector.", call. = FALSE)
  }
  neutral_labels <- as.numeric(neutral_labels)
  if (!is.logical(enforce_order) || length(enforce_order) != 1L || is.na(enforce_order)) {
    stop("enforce_order must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.numeric(eps) || length(eps) != 1L || is.na(eps) ||
      !is.finite(eps) || eps <= 0) {
    stop("eps must be one positive finite numeric value.", call. = FALSE)
  }
  new_biogeme_expression(
    paste0("ordered_", distribution),
    eta = eta,
    cutpoints = cutpoints,
    alternative = alternative,
    categories = categories,
    neutral_labels = neutral_labels,
    enforce_order = isTRUE(enforce_order),
    eps = as.numeric(eps)
  )
}

#' Construct a native ordered-logit log-likelihood expression
#'
#' @param eta Latent-index expression.
#' @param cutpoints Non-empty list of cutpoint expressions.
#' @param alternative Observed response expression.
#' @param categories Optional ordered numeric labels.
#' @param neutral_labels Optional numeric labels whose contribution is one.
#' @param enforce_order Whether native Biogeme should enforce ordered cutpoints.
#' @param eps Positive numerical lower bound used by native Biogeme.
#' @return A Biogeme expression compiled to native `OrderedLogLogit`.
#' @export
ordered_logit_log_probability <- function(
    eta,
    cutpoints,
    alternative,
    categories = NULL,
    neutral_labels = numeric(),
    enforce_order = TRUE,
    eps = 1e-12
) {
  ordered_response_log_probability(
    distribution = "logit",
    eta = eta,
    cutpoints = cutpoints,
    alternative = alternative,
    categories = categories,
    neutral_labels = neutral_labels,
    enforce_order = enforce_order,
    eps = eps
  )
}

#' Construct a native ordered-probit log-likelihood expression
#'
#' Ordered probit uses the native normal CDF with the same cutpoint and
#' category semantics as [ordered_logit_log_probability()].
#'
#' @param eta Latent-index expression.
#' @param cutpoints Non-empty list of cutpoint expressions.
#' @param alternative Observed response expression.
#' @param categories Optional ordered numeric labels.
#' @param neutral_labels Optional numeric labels whose contribution is one.
#' @param enforce_order Whether native Biogeme should enforce ordered cutpoints.
#' @param eps Positive numerical lower bound used by native Biogeme.
#' @return A Biogeme expression compiled to native `OrderedLogProbit`.
#' @export
ordered_probit_log_probability <- function(
    eta,
    cutpoints,
    alternative,
    categories = NULL,
    neutral_labels = numeric(),
    enforce_order = TRUE,
    eps = 1e-12
) {
  ordered_response_log_probability(
    distribution = "probit",
    eta = eta,
    cutpoints = cutpoints,
    alternative = alternative,
    categories = categories,
    neutral_labels = neutral_labels,
    enforce_order = enforce_order,
    eps = eps
  )
}

#' Normal cumulative distribution function
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
normal_cdf <- function(x) {
  biogeme_function("normal_cdf", list(x))
}

#' Normal probability density function
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
normal_pdf <- function(x) {
  biogeme_function("normal_pdf", list(x))
}

#' Numerically safe exponential
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
safe_exp <- function(x) {
  biogeme_function("safe_exp", list(x))
}

#' @rdname logzero
#' @export
safe_log <- logzero

#' Square root of a Biogeme expression
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
sqrt.biogeme_expression <- function(x) {
  biogeme_function("sqrt", list(x))
}

#' Absolute value of a Biogeme expression
#' @param x A scalar or Biogeme expression.
#' @return A Biogeme expression.
#' @export
abs.biogeme_expression <- function(x) {
  biogeme_function("abs", list(x))
}

#' Minimum of Biogeme expressions
#' @param ... Scalar values or Biogeme expressions.
#' @return A Biogeme expression.
#' @export
biogeme_min <- function(...) {
  args <- list(...)
  if (length(args) == 0L) {
    stop("biogeme_min requires at least one argument.", call. = FALSE)
  }
  biogeme_function("min", args)
}

#' Maximum of Biogeme expressions
#' @param ... Scalar values or Biogeme expressions.
#' @return A Biogeme expression.
#' @export
biogeme_max <- function(...) {
  args <- list(...)
  if (length(args) == 0L) {
    stop("biogeme_max requires at least one argument.", call. = FALSE)
  }
  biogeme_function("max", args)
}

#' @export
#' @method Summary biogeme_expression
Summary.biogeme_expression <- function(..., na.rm = FALSE) {
  if (!.Generic %in% c("min", "max")) {
    stop("Unsupported Summary operation for a Biogeme expression: ", .Generic, call. = FALSE)
  }
  if (isTRUE(na.rm)) stop("na.rm is not meaningful for symbolic expressions.", call. = FALSE)
  if (.Generic == "min") biogeme_min(...) else biogeme_max(...)
}

#' Select an expression from a numeric mapping
#' @param mapping Named list whose names are integer keys.
#' @param index Numeric or Biogeme expression used as the key.
#' @return A Biogeme expression.
#' @export
Elem <- function(mapping, index) {
  if (!is.list(mapping) || is.null(names(mapping)) || length(mapping) == 0L ||
      anyNA(names(mapping)) || any(!nzchar(names(mapping)))) {
    stop("mapping must be a non-empty named list.", call. = FALSE)
  }
  keys <- suppressWarnings(as.integer(names(mapping)))
  if (anyNA(keys) || anyDuplicated(keys)) {
    stop("mapping names must be unique integer keys.", call. = FALSE)
  }
  node <- new_biogeme_expression(
    "elem",
    mapping = lapply(mapping, as_biogeme_expression),
    keys = keys,
    index = as_biogeme_expression(index)
  )
  names(node$mapping) <- as.character(keys)
  node
}

#' @rdname Elem
#' @export
elem <- Elem

#' Create a named Biogeme draw node
#' @param name Draw variable name.
#' @param draw_type Biogeme draw type, such as `NORMAL` or `UNIFORM_HALTON2`.
#' @return A Biogeme expression.
#' @export
draw <- function(name, draw_type = "NORMAL") {
  new_biogeme_expression(
    "draw",
    name = validate_name(name, "name"),
    draw_type = validate_name(draw_type, "draw_type")
  )
}

#' Declare first-class draw metadata for a model
#' @param name Draw variable name.
#' @param draw_type Native draw type.
#' @param number_of_draws Optional draw count.
#' @param seed Optional draw seed.
#' @param matrix Optional user-supplied numeric draw matrix.
#' @param generator Optional bridge generator name. Supported bridge-owned
#'   generators include `"TRIANGULAR"`, `"HALTON13"`, and
#'   `"HALTON13_ANTI"`. The latter two match the custom base-13, skip-10
#'   generators used by the Monte Carlo examples.
#' @return A `biogeme_draws` object.
#' @export
biogeme_draws <- function(
    name,
    draw_type = "NORMAL",
    number_of_draws = NULL,
    seed = NULL,
    matrix = NULL,
    generator = NULL
) {
  name <- validate_name(name, "name")
  draw_type <- validate_name(draw_type, "draw_type")
  if (!is.null(number_of_draws) &&
      (!is.numeric(number_of_draws) || length(number_of_draws) != 1L ||
       is.na(number_of_draws) || number_of_draws < 1 || number_of_draws != floor(number_of_draws))) {
    stop("number_of_draws must be NULL or one positive integer.", call. = FALSE)
  }
  if (!is.null(seed) && (!is.numeric(seed) || length(seed) != 1L || is.na(seed) || !is.finite(seed))) {
    stop("seed must be NULL or one finite numeric value.", call. = FALSE)
  }
  if (!is.null(matrix) && (!is.numeric(matrix) || is.null(dim(matrix)) || length(dim(matrix)) != 2L)) {
    stop("matrix must be NULL or a numeric matrix.", call. = FALSE)
  }
  if (!is.null(generator) &&
      (!is.character(generator) || length(generator) != 1L ||
       is.na(generator) || !nzchar(generator))) {
    stop("generator must be NULL or one non-empty character string.", call. = FALSE)
  }
  structure(
    list(
      name = name,
      draw_type = draw_type,
      number_of_draws = if (is.null(number_of_draws)) NULL else as.integer(number_of_draws),
      seed = if (is.null(seed)) NULL else as.numeric(seed),
      matrix = matrix,
      generator = generator
    ),
    class = "biogeme_draws"
  )
}

#' Create a random variable for native numerical integration
#' @param name Random variable name.
#' @return A Biogeme expression.
#' @export
random_variable <- function(name) {
  new_biogeme_expression("random_variable", name = validate_name(name, "name"))
}

#' Store a simulated individual-level parameter in Bayesian results
#'
#' `distributed_parameter()` is a symbolic wrapper around native Biogeme's
#' `DistributedParameter` expression. It preserves the named variable in the
#' Bayesian output while its child expression supplies the value used by the
#' likelihood. The wrapper is compiled once with the rest of the expression
#' tree; it is not evaluated by R during sampling.
#' @param name Name used for the stored native variable.
#' @param expression Child expression, usually a location parameter plus a
#'   scale parameter multiplied by [draw()].
#' @return A Biogeme expression.
#' @export
distributed_parameter <- function(name, expression) {
  new_biogeme_expression(
    "distributed_parameter",
    name = validate_name(name, "name"),
    operand = as_biogeme_expression(expression)
  )
}

#' Integrate an expression over a standard normal random variable
#' @param expression Integrand expression.
#' @param name Random variable name.
#' @param number_of_quadrature_points Number of quadrature points.
#' @return A Biogeme expression.
#' @export
integrate_normal <- function(expression, name, number_of_quadrature_points = 30L) {
  if (!is.numeric(number_of_quadrature_points) || length(number_of_quadrature_points) != 1L ||
      is.na(number_of_quadrature_points) || number_of_quadrature_points < 1 ||
      number_of_quadrature_points != floor(number_of_quadrature_points)) {
    stop("number_of_quadrature_points must be one positive integer.", call. = FALSE)
  }
  new_biogeme_expression(
    "integrate_normal",
    operand = as_biogeme_expression(expression),
    name = validate_name(name, "name"),
    number_of_quadrature_points = as.integer(number_of_quadrature_points)
  )
}

#' Monte Carlo average of an expression
#' @param expression Expression to integrate.
#' @return A Biogeme expression.
#' @export
monte_carlo <- function(expression) {
  new_biogeme_expression("monte_carlo", operand = as_biogeme_expression(expression))
}

#' Aggregate an observation-level likelihood over a panel trajectory
#' @param expression Observation-level likelihood or probability expression.
#' @return A Biogeme expression.
#' @export
panel_likelihood_trajectory <- function(expression) {
  new_biogeme_expression(
    "panel_likelihood_trajectory",
    operand = as_biogeme_expression(expression)
  )
}

#' Symbolically differentiate an expression with respect to a named variable
#' @param expression Expression to differentiate.
#' @param name Variable or parameter name.
#' @return A Biogeme expression.
#' @export
derive <- function(expression, name) {
  new_biogeme_expression(
    "derive",
    operand = as_biogeme_expression(expression),
    name = validate_name(name, "name")
  )
}

#' @rdname derive
#' @export
Derive <- derive

#' Box--Cox transformation
#' @param x Expression to transform.
#' @param lambda Box--Cox exponent expression.
#' @return A Biogeme expression.
#' @export
boxcox <- function(x, lambda) {
  new_biogeme_expression(
    "boxcox",
    left = as_biogeme_expression(x),
    right = as_biogeme_expression(lambda)
  )
}

validate_piecewise_thresholds <- function(thresholds) {
  if (!is.numeric(thresholds) && !is.list(thresholds)) {
    stop("thresholds must be numeric, with optional NULL endpoints.", call. = FALSE)
  }
  thresholds <- as.list(thresholds)
  if (length(thresholds) < 2L || all(vapply(thresholds, is.null, logical(1)))) {
    stop("thresholds must contain at least two values and not be all NULL.", call. = FALSE)
  }
  if (any(vapply(thresholds[-c(1L, length(thresholds))], is.null, logical(1)))) {
    stop("Only the first and last piecewise thresholds may be NULL.", call. = FALSE)
  }
  values <- vapply(thresholds, function(value) {
    if (is.null(value)) return(NA_real_)
    if (!is.numeric(value) || length(value) != 1L || is.na(value) || !is.finite(value)) {
      stop("Piecewise thresholds must be finite numeric values or NULL.", call. = FALSE)
    }
    as.numeric(value)
  }, numeric(1))
  finite <- values[is.finite(values)]
  if (length(finite) > 1L && any(diff(finite) <= 0)) {
    stop("Finite piecewise thresholds must be strictly increasing.", call. = FALSE)
  }
  thresholds
}

#' Piecewise-linear expression using native Biogeme naming rules
#' @param expression Variable expression.
#' @param thresholds Piecewise thresholds; only endpoints may be NULL.
#' @param betas Optional list of parameter expressions.
#' @param transform Whether to return the full formula or the transformed
#'   variable form.
#' @return A Biogeme expression.
#' @export
piecewise <- function(expression, thresholds, betas = NULL, transform = c("formula", "variable")) {
  transform <- match.arg(transform)
  thresholds <- validate_piecewise_thresholds(thresholds)
  expected_betas <- if (transform == "variable") length(thresholds) - 2L else length(thresholds) - 1L
  if (!is.null(betas)) {
    if (!is.list(betas) || length(betas) != expected_betas) {
      stop("betas must contain one expression per piecewise interval.", call. = FALSE)
    }
    betas <- lapply(betas, as_biogeme_expression)
  }
  new_biogeme_expression(
    "piecewise",
    operand = as_biogeme_expression(expression),
    thresholds = thresholds,
    betas = betas,
    transform = transform
  )
}

format_biogeme_expression <- function(x) {
  if (!is_biogeme_expression(x)) {
    return(format(x))
  }
  if (identical(x$kind, "numeric")) {
    return(format(x$value))
  }
  if (identical(x$kind, "variable")) {
    return(x$name)
  }
  if (identical(x$kind, "beta")) {
    return(x$name)
  }
  if (identical(x$kind, "unary")) {
    operand <- format_biogeme_expression(x$operand)
    if (identical(x$operator, "-")) {
      return(paste0("(-", operand, ")"))
    }
    return(paste0(x$operator, "(", operand, ")"))
  }
  if (identical(x$kind, "function")) {
    args <- vapply(x$args, format_biogeme_expression, character(1))
    return(paste0(x$operator, "(", paste(args, collapse = ", "), ")"))
  }
  if (identical(x$kind, "elem")) {
    entries <- vapply(seq_along(x$mapping), function(i) {
      paste0(x$keys[[i]], ": ", format_biogeme_expression(x$mapping[[i]]))
    }, character(1))
    return(paste0("Elem({", paste(entries, collapse = ", "), "}, ",
      format_biogeme_expression(x$index), ")"))
  }
  if (identical(x$kind, "catalog")) {
    return(paste0(
      "Catalog(", x$name, ": ",
      paste(names(x$expressions), collapse = " | "), ")"
    ))
  }
  if (identical(x$kind, "draw")) {
    return(paste0("Draws(\"", x$name, "\", \"", x$draw_type, "\")"))
  }
  if (identical(x$kind, "random_variable")) {
    return(paste0("RandomVariable(\"", x$name, "\")"))
  }
  if (identical(x$kind, "distributed_parameter")) {
    return(paste0(
      "DistributedParameter(\"", x$name, "\", ",
      format_biogeme_expression(x$operand), ")"
    ))
  }
  if (identical(x$kind, "integrate_normal")) {
    return(paste0("IntegrateNormal(", format_biogeme_expression(x$operand), ", \"",
      x$name, "\")"))
  }
  if (identical(x$kind, "monte_carlo")) {
    return(paste0("MonteCarlo(", format_biogeme_expression(x$operand), ")"))
  }
  if (identical(x$kind, "panel_likelihood_trajectory")) {
    return(paste0("PanelLikelihoodTrajectory(", format_biogeme_expression(x$operand), ")"))
  }
  if (identical(x$kind, "derive")) {
    return(paste0("Derive(", format_biogeme_expression(x$operand), ", \"", x$name, "\")"))
  }
  if (identical(x$kind, "boxcox")) {
    return(paste0("BoxCox(", format_biogeme_expression(x$left), ", ",
      format_biogeme_expression(x$right), ")"))
  }
  if (identical(x$kind, "piecewise")) {
    thresholds <- vapply(x$thresholds, function(value) {
      if (is.null(value)) return("NULL")
      format(value)
    }, character(1))
    return(paste0("piecewise(", format_biogeme_expression(x$operand), ", c(",
      paste(thresholds, collapse = ", "), "))"))
  }
  if (identical(x$kind, "binary")) {
    left <- format_biogeme_expression(x$left)
    right <- format_biogeme_expression(x$right)
    return(paste0("(", left, " ", x$operator, " ", right, ")"))
  }
  if (identical(x$kind, "linear_utility")) {
    terms <- vapply(x$terms, function(term) {
      paste0(
        format_biogeme_expression(term$beta), " * ",
        format_biogeme_expression(term$x)
      )
    }, character(1))
    return(paste(terms, collapse = " + "))
  }
  if (identical(x$kind, "segmented_beta")) {
    return(paste0(x$beta$name, " segmented by ", length(x$segmentations), " variable(s)"))
  }
  if (x$kind %in% c("logit_probability", "logit_log_probability")) {
    alternative <- if (is_biogeme_expression(x$alternative)) {
      format_biogeme_expression(x$alternative)
    } else {
      as.character(x$alternative)
    }
    return(paste0(
      if (x$kind == "logit_log_probability") "logit log probability" else "logit probability",
      " for alternative ", alternative
    ))
  }
  if (x$kind %in% c("ordered_logit", "ordered_probit")) {
    return(paste0(
      "ordered ",
      if (x$kind == "ordered_logit") "logit" else "probit",
      " log probability for response ",
      format_biogeme_expression(x$alternative)
    ))
  }
  if (x$kind %in% c(
      "nested_logit",
      "nested_logit_mev_mu",
      "nested_probability",
      "nested_endogenous_sampling",
      "cross_nested_logit",
      "cross_nested_logit_mu",
      "cross_nested_probability"
  )) {
    alternative <- if (is_biogeme_expression(x$alternative)) {
      format_biogeme_expression(x$alternative)
    } else {
      as.character(x$alternative)
    }
    return(paste0("nested logit probability for alternative ", alternative))
  }
  "<unknown Biogeme expression>"
}

#' @export
format.biogeme_expression <- function(x, ...) {
  format_biogeme_expression(x)
}

#' @export
print.biogeme_expression <- function(x, ...) {
  cat(format_biogeme_expression(x), "\n")
  invisible(x)
}

#' Collect parameter names from an expression
#'
#' @param expression A Biogeme expression.
#' @return A character vector of parameter names.
#' @keywords internal
collect_biogeme_parameters <- function(expression) {
  if (!is_biogeme_expression(expression)) {
    stop("expression must be a Biogeme expression.", call. = FALSE)
  }
  if (identical(expression$kind, "beta")) {
    return(expression$name)
  }
  if (identical(expression$kind, "segmented_beta")) {
    names <- paste0(expression$beta$name, "_ref")
    for (segmentation in expression$segmentations) {
      labels <- unname(segmentation$mapping)
      labels <- labels[labels != segmentation$reference]
      names <- c(names, paste0(expression$beta$name, "_diff_", labels))
    }
    return(names)
  }
  children <- biogeme_expression_children(expression)
  if (length(children) == 0L) return(character())
  unique(unlist(lapply(children, collect_biogeme_parameters), use.names = FALSE))
}

biogeme_expression_children <- function(expression) {
  if (!is_biogeme_expression(expression)) {
    stop("expression must be a Biogeme expression.", call. = FALSE)
  }
  kind <- expression$kind
  if (kind == "binary") return(list(expression$left, expression$right))
  if (kind == "linear_utility") {
    return(unlist(lapply(expression$terms, function(term) list(term$beta, term$x)), recursive = FALSE))
  }
  if (kind == "catalog") return(unname(expression$expressions))
  if (kind == "segmented_beta") return(list())
  if (kind %in% c("logit_probability", "logit_log_probability")) {
    return(c(
      unname(expression$utilities),
      unname(expression$availability),
      if (is_biogeme_expression(expression$alternative)) list(expression$alternative) else list()
    ))
  }
  if (kind %in% c("ordered_logit", "ordered_probit")) {
    return(c(
      list(expression$eta),
      expression$cutpoints,
      list(expression$alternative)
    ))
  }
  if (kind %in% c("nested_probability", "nested_logit", "nested_logit_mev_mu")) {
    return(c(
      unname(expression$utilities),
      unname(expression$availability),
      if (is_biogeme_expression(expression$alternative)) list(expression$alternative) else list(),
      lapply(expression$nests$nests, function(nest) nest$nest_parameter),
      if (kind == "nested_logit_mev_mu") list(expression$scale_parameter) else list()
    ))
  }
  if (kind == "nested_endogenous_sampling") {
    return(c(
      unname(expression$utilities),
      unname(expression$availability),
      unname(expression$correction),
      if (is_biogeme_expression(expression$alternative)) list(expression$alternative) else list(),
      lapply(expression$nests$nests, function(nest) nest$nest_parameter)
    ))
  }
  if (kind %in% c(
      "cross_nested_logit",
      "cross_nested_logit_mu",
      "cross_nested_probability"
  )) {
    return(c(
      unname(expression$utilities),
      unname(expression$availability),
      if (is_biogeme_expression(expression$alternative)) list(expression$alternative) else list(),
      unlist(lapply(expression$nests$nests, function(nest) {
        c(list(nest$nest_parameter), unname(nest$allocation))
      }), recursive = FALSE),
      if (kind == "cross_nested_logit_mu") list(expression$scale_parameter) else list()
    ))
  }
  if (kind %in% c("unary", "derive", "integrate_normal", "monte_carlo",
                  "panel_likelihood_trajectory", "distributed_parameter")) {
    return(list(expression$operand))
  }
  if (kind == "function") return(expression$args)
  if (kind == "elem") return(c(list(expression$index), unname(expression$mapping)))
  if (kind == "boxcox") return(list(expression$left, expression$right))
  if (kind == "piecewise") {
    return(c(list(expression$operand), if (is.null(expression$betas)) list() else expression$betas))
  }
  list()
}

#' Collect data-variable names from an expression
#' @param expression A Biogeme expression.
#' @return A character vector.
#' @keywords internal
collect_biogeme_variables <- function(expression) {
  if (!is_biogeme_expression(expression)) {
    stop("expression must be a Biogeme expression.", call. = FALSE)
  }
  result <- if (identical(expression$kind, "variable")) expression$name else character()
  if (identical(expression$kind, "segmented_beta")) {
    return(vapply(expression$segmentations, function(segmentation) segmentation$variable, character(1)))
  }
  children <- biogeme_expression_children(expression)
  if (length(children) > 0L) {
    result <- c(result, unlist(lapply(children, collect_biogeme_variables), use.names = FALSE))
  }
  unique(result)
}