R/mdcev.R

In rbiogeme: Interface to 'Biogeme'

#' Construct an MDCEV model specification
#'
#' This is a declarative model object.  The alternative utilities, shape
#' parameters, and observed consumptions remain neutral R expression trees
#' until the Python bridge sends the complete specification to the native
#' Biogeme engine. The likelihood and forecasting algorithms
#' are therefore provided by Biogeme itself.
#'
#' @param database A [biogeme_database()] object.
#' @param model_type One of `"gamma_profile"`, `"generalized"`,
#'   `"translated"`, or `"non_monotonic"`.
#' @param baseline_utilities Non-empty named list of baseline utility
#'   expressions, keyed by integer alternative code.
#' @param gamma_parameters Named list of gamma expressions, keyed by the same
#'   alternatives.  One value may be `NULL` to declare an outside good.
#' @param alpha_parameters Named list of alpha expressions. Required for the
#'   generalized, translated, and non-monotonic variants.
#' @param mu_utilities Named list of non-monotonic utility expressions.
#' @param scale_parameter Optional scale expression.
#' @param prices Optional named list of price expressions. Supported by the
#'   gamma-profile and generalized native classes.
#' @param weights Optional observation-weight expression.
#' @param number_of_chosen_alternatives Expression containing the number of
#'   goods chosen in each observation.
#' @param consumed_quantities Named list of observed consumption expressions.
#' @param subset Optional logical expression selecting observations to remove
#'   when it is false, using the same semantics as [biogeme_model()].
#' @param control Optional [biogeme_control()] object.
#' @return An object of class `biogeme_mdcev_model`.
#' @export
biogeme_mdcev_model <- function(
    database,
    model_type = c("gamma_profile", "generalized", "translated", "non_monotonic"),
    baseline_utilities,
    gamma_parameters,
    alpha_parameters = NULL,
    mu_utilities = NULL,
    scale_parameter = NULL,
    prices = NULL,
    weights = NULL,
    number_of_chosen_alternatives,
    consumed_quantities,
    subset = NULL,
    control = NULL
) {
  validate_biogeme_database(database)
  model_type <- match.arg(model_type)

  validate_map <- function(values, argument, allow_null = FALSE) {
    if (!is.list(values) || length(values) == 0L || is.null(names(values)) ||
        anyNA(names(values)) || any(!nzchar(names(values))) ||
        anyDuplicated(names(values))) {
      stop(argument, " must be a non-empty named list.", call. = FALSE)
    }
    numeric_keys <- suppressWarnings(as.numeric(names(values)))
    if (anyNA(numeric_keys) || any(!is.finite(numeric_keys)) ||
        any(numeric_keys != floor(numeric_keys)) ||
        anyDuplicated(as.integer(numeric_keys))) {
      stop(argument, " names must be unique finite integer codes.", call. = FALSE)
    }
    result <- lapply(values, function(value) {
      if (is.null(value)) {
        if (!isTRUE(allow_null)) {
          stop(argument, " cannot contain NULL values.", call. = FALSE)
        }
        return(NULL)
      }
      as_biogeme_expression(value)
    })
    names(result) <- as.character(as.integer(numeric_keys))
    result
  }

  baseline_utilities <- validate_map(baseline_utilities, "baseline_utilities")
  gamma_parameters <- validate_map(
    gamma_parameters,
    "gamma_parameters",
    allow_null = TRUE
  )
  alternatives <- names(baseline_utilities)
  if (!setequal(names(gamma_parameters), alternatives)) {
    stop("gamma_parameters must contain exactly the baseline alternatives.", call. = FALSE)
  }
  if (sum(vapply(gamma_parameters, is.null, logical(1))) > 1L) {
    stop("gamma_parameters may contain at most one NULL outside-good entry.", call. = FALSE)
  }

  requires_alpha <- model_type %in% c("generalized", "translated", "non_monotonic")
  if (requires_alpha && is.null(alpha_parameters)) {
    stop(model_type, " requires alpha_parameters.", call. = FALSE)
  }
  if (!requires_alpha && !is.null(alpha_parameters)) {
    stop("alpha_parameters is only used by the generalized, translated, and non_monotonic variants.", call. = FALSE)
  }
  if (!is.null(alpha_parameters)) {
    alpha_parameters <- validate_map(alpha_parameters, "alpha_parameters")
    if (!setequal(names(alpha_parameters), alternatives)) {
      stop("alpha_parameters must contain exactly the baseline alternatives.", call. = FALSE)
    }
  }

  if (identical(model_type, "non_monotonic") && is.null(mu_utilities)) {
    stop("non_monotonic requires mu_utilities.", call. = FALSE)
  }
  if (!identical(model_type, "non_monotonic") && !is.null(mu_utilities)) {
    stop("mu_utilities is only used by the non_monotonic variant.", call. = FALSE)
  }
  if (!is.null(mu_utilities)) {
    mu_utilities <- validate_map(mu_utilities, "mu_utilities")
    if (!setequal(names(mu_utilities), alternatives)) {
      stop("mu_utilities must contain exactly the baseline alternatives.", call. = FALSE)
    }
  }

  prices_allowed <- model_type %in% c("gamma_profile", "generalized")
  if (!prices_allowed && !is.null(prices)) {
    stop("prices are supported only by the gamma_profile and generalized variants.", call. = FALSE)
  }
  if (!is.null(prices)) {
    prices <- validate_map(prices, "prices")
    if (!setequal(names(prices), alternatives)) {
      stop("prices must contain exactly the baseline alternatives.", call. = FALSE)
    }
  }

  if (!is.null(scale_parameter)) scale_parameter <- as_biogeme_expression(scale_parameter)
  if (!is.null(weights)) weights <- as_biogeme_expression(weights)
  number_of_chosen_alternatives <- as_biogeme_expression(number_of_chosen_alternatives)
  consumed_quantities <- validate_map(consumed_quantities, "consumed_quantities")
  if (!setequal(names(consumed_quantities), alternatives)) {
    stop("consumed_quantities must contain exactly the baseline alternatives.", call. = FALSE)
  }
  if (!is.null(subset)) subset <- as_biogeme_expression(subset)
  if (!is.null(control) && !is.list(control)) {
    stop("control must be a list returned by biogeme_control().", call. = FALSE)
  }

  expressions <- c(
    unname(baseline_utilities),
    Filter(Negate(is.null), unname(gamma_parameters)),
    if (is.null(alpha_parameters)) list() else unname(alpha_parameters),
    if (is.null(mu_utilities)) list() else unname(mu_utilities),
    if (is.null(prices)) list() else unname(prices),
    list(scale_parameter, weights, number_of_chosen_alternatives, subset),
    unname(consumed_quantities)
  )
  expressions <- Filter(Negate(is.null), expressions)
  validate_expression_variables(expressions, database, argument = "model expression")

  parameters <- list()
  for (expression in expressions) {
    parameters <- collect_biogeme_betas(expression, parameters)
  }

  structure(
    list(
      kind = "mdcev",
      database = database,
      model_type = model_type,
      baseline_utilities = baseline_utilities,
      gamma_parameters = gamma_parameters,
      alpha_parameters = alpha_parameters,
      mu_utilities = mu_utilities,
      scale_parameter = scale_parameter,
      prices = prices,
      weights = weights,
      number_of_chosen_alternatives = number_of_chosen_alternatives,
      consumed_quantities = consumed_quantities,
      subset = subset,
      control = control,
      parameters = parameters
    ),
    class = c("biogeme_mdcev_model", "biogeme_model")
  )
}

#' @rdname biogeme_mdcev_model
#' @export
mdcev_model <- biogeme_mdcev_model

#' @export
print.biogeme_mdcev_model <- function(x, ...) {
  cat(
    "Biogeme MDCEV model (", x$model_type, ") with ",
    length(x$baseline_utilities), " alternatives and ",
    length(x$parameters), " parameter definition(s)\n",
    sep = ""
  )
  invisible(x)
}

validate_mdcev_operation_model <- function(model, operation) {
  if (!inherits(model, "biogeme_mdcev_model")) {
    biogeme_abort(
      "model must be a biogeme_mdcev_model.",
      class = "biogeme_specification_error",
      operation = operation,
      suggestion = "construct the model with biogeme_mdcev_model()"
    )
  }
  invisible(model)
}

#' Estimate an MDCEV model with native Biogeme
#'
#' The MDCEV likelihood is generated by the native Biogeme MDCEV class after
#' compilation.  This convenience wrapper has the same fresh-estimation
#' semantics as [estimate()].
#'
#' @param model A [biogeme_mdcev_model()] object.
#' @param model_name Native Biogeme model name.
#' @param controls Named native Biogeme controls.
#' @param starting_values Optional named numeric vector of starting values.
#' @param run_bootstrap Whether to run native bootstrap re-estimation.
#' @param yaml_file_name Optional path for native YAML output.
#' @param control Optional [biogeme_control()] object.
#' @return A `biogeme_fit` object.
#' @export
mdcev_estimate <- function(
    model,
    model_name = "rbiogeme_mdcev",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE,
    yaml_file_name = NULL,
    control = NULL
) {
  validate_mdcev_operation_model(model, "MDCEV estimation")
  estimate(
    model = model,
    model_name = model_name,
    controls = controls,
    starting_values = starting_values,
    run_bootstrap = run_bootstrap,
    yaml_file_name = yaml_file_name,
    control = control
  )
}

as_mdcev_forecast_data_frame <- function(value) {
  value <- reticulate::py_to_r(value)
  if (!is.list(value) || is.null(value$index) || is.null(value$columns)) {
    stop("Biogeme returned an invalid MDCEV forecast table.", call. = FALSE)
  }
  index <- as.character(value$index)
  columns <- as.character(value$columns)
  rows <- value$data
  if (length(rows) == 0L) {
    result <- as.data.frame(
      setNames(replicate(length(columns), numeric(), simplify = FALSE), columns),
      check.names = FALSE
    )
  } else {
    result <- as.data.frame(
      setNames(lapply(seq_along(columns), function(column) {
        unlist(lapply(rows, function(row) row[[column]]), use.names = FALSE)
      }), columns),
      check.names = FALSE
    )
  }
  rownames(result) <- index
  result
}

validate_mdcev_positive_integer <- function(value, argument) {
  if (!is.numeric(value) || length(value) != 1L || is.na(value) ||
      !is.finite(value) || value <= 0 || value != floor(value)) {
    stop(argument, " must be one positive integer.", call. = FALSE)
  }
  as.integer(value)
}

validate_mdcev_seed <- function(value) {
  if (!is.numeric(value) || length(value) != 1L || is.na(value) ||
      !is.finite(value) || value < 0 || value != floor(value)) {
    stop("seed must be one non-negative integer.", call. = FALSE)
  }
  as.integer(value)
}

validate_mdcev_budget <- function(value, argument = "total_budget") {
  if (!is.numeric(value) || length(value) != 1L || is.na(value) ||
      !is.finite(value) || value <= 0) {
    stop(argument, " must be one positive finite number.", call. = FALSE)
  }
  as.numeric(value)
}

validate_mdcev_tolerance <- function(value, argument) {
  if (!is.numeric(value) || length(value) != 1L || is.na(value) ||
      !is.finite(value) || value <= 0) {
    stop(argument, " must be one positive finite number.", call. = FALSE)
  }
  as.numeric(value)
}

validate_mdcev_epsilons <- function(epsilons, number_of_alternatives) {
  if (!is.list(epsilons) || length(epsilons) == 0L) {
    stop("epsilons must be a non-empty list of numeric matrices.", call. = FALSE)
  }
  result <- lapply(epsilons, function(value) {
    if (!is.matrix(value) || !is.numeric(value) || length(dim(value)) != 2L ||
        ncol(value) != number_of_alternatives || anyNA(value) ||
        any(!is.finite(value))) {
      stop(
        "Each epsilon entry must be a finite numeric matrix with ",
        number_of_alternatives, " columns.",
        call. = FALSE
      )
    }
    value
  })
  result
}

#' Generate native MDCEV error-term draws
#'
#' @param model A [biogeme_mdcev_model()] object.
#' @param number_of_observations Number of database observations.
#' @param number_of_draws Number of Gumbel draws per observation.
#' @param seed Optional native NumPy seed.
#' @return A list of numeric matrices, one matrix per observation.
#' @export
mdcev_generate_epsilons <- function(
    model,
    number_of_observations,
    number_of_draws,
    seed = NULL
) {
  validate_mdcev_operation_model(model, "MDCEV epsilon generation")
  number_of_observations <- validate_mdcev_positive_integer(
    number_of_observations, "number_of_observations"
  )
  number_of_draws <- validate_mdcev_positive_integer(number_of_draws, "number_of_draws")
  if (!is.null(seed)) seed <- validate_mdcev_seed(seed)
  compiled <- biogeme_compile_model(model)
  result <- tryCatch({
    bridge <- biogeme_bridge()
    raw <- if (is.null(seed)) {
      bridge$generate_epsilons_mdcev_biogeme(
        compiled_model = compiled,
        number_of_observations = number_of_observations,
        number_of_draws = number_of_draws,
        seed = NULL
      )
    } else {
      bridge$generate_epsilons_mdcev_biogeme(
        compiled_model = compiled,
        number_of_observations = number_of_observations,
        number_of_draws = number_of_draws,
        seed = seed
      )
    }
    lapply(reticulate::py_to_r(raw), function(value) {
      matrix(as.numeric(unlist(value, use.names = FALSE)), ncol = length(model$baseline_utilities), byrow = TRUE)
    })
  }, error = function(error) {
    biogeme_rethrow(
      error,
      class = "biogeme_simulation_error",
      operation = "native MDCEV epsilon generation",
      suggestion = "check the observation count, draw count, and Python Biogeme environment"
    )
  })
  result
}

#' Forecast MDCEV consumption using native Biogeme algorithms
#'
#' @param model A [biogeme_mdcev_model()] object.
#' @param fit A fresh [mdcev_estimate()] result.
#' @param database Optional database used only for forecasting. This permits
#'   forecasting a native row subset without changing the estimation fit.
#' @param total_budget Positive total budget.
#' @param epsilons Optional list of Gumbel-draw matrices. If omitted, native
#'   draws are generated using `number_of_draws` and `seed`.
#' @param number_of_draws Number of draws when `epsilons` is omitted.
#' @param seed Optional native NumPy seed when draws are generated.
#' @param brute_force Whether to use native brute-force optimization.
#' @param tolerance_dual Native dual-variable tolerance.
#' @param tolerance_budget Native budget-constraint tolerance.
#' @return A `biogeme_mdcev_forecast` object containing R data frames.
#' @export
mdcev_forecast <- function(
    model,
    fit,
    database = NULL,
    total_budget,
    epsilons = NULL,
    number_of_draws = NULL,
    seed = NULL,
    brute_force = FALSE,
    tolerance_dual = 1e-10,
    tolerance_budget = 1e-10
) {
  validate_mdcev_operation_model(model, "MDCEV forecasting")
  if (!inherits(fit, "biogeme_fit")) {
    stop("fit must be a biogeme_fit returned by mdcev_estimate().", call. = FALSE)
  }
  if (!is.null(database)) validate_biogeme_database(database)
  total_budget <- validate_mdcev_budget(total_budget)
  tolerance_dual <- validate_mdcev_tolerance(tolerance_dual, "tolerance_dual")
  tolerance_budget <- validate_mdcev_tolerance(tolerance_budget, "tolerance_budget")
  if (!is.logical(brute_force) || length(brute_force) != 1L || is.na(brute_force)) {
    stop("brute_force must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.null(number_of_draws)) {
    number_of_draws <- validate_mdcev_positive_integer(number_of_draws, "number_of_draws")
  }
  if (!is.null(seed)) seed <- validate_mdcev_seed(seed)
  if (is.null(epsilons) && is.null(number_of_draws)) {
    stop("number_of_draws is required when epsilons is NULL.", call. = FALSE)
  }
  if (!is.null(epsilons)) {
    epsilons <- validate_mdcev_epsilons(epsilons, length(model$baseline_utilities))
  }

  compiled <- biogeme_compile_model(model, database_override = database)
  fit_for_bridge <- fit
  fit_for_bridge$model <- NULL
  raw <- tryCatch({
    bridge <- biogeme_bridge()
    epsilon_payload <- if (is.null(epsilons)) NULL else reticulate::r_to_py(epsilons)
    result <- bridge$forecast_mdcev_biogeme(
      compiled_model = compiled,
      estimation_results = reticulate::r_to_py(fit_for_bridge),
      total_budget = total_budget,
      epsilons = epsilon_payload,
      number_of_draws = if (is.null(number_of_draws)) NULL else number_of_draws,
      seed = if (is.null(seed)) NULL else seed,
      brute_force = isTRUE(brute_force),
      tolerance_dual = tolerance_dual,
      tolerance_budget = tolerance_budget
    )
    reticulate::py_to_r(result)
  }, error = function(error) {
    biogeme_rethrow(
      error,
      class = "biogeme_simulation_error",
      operation = "native MDCEV forecasting",
      suggestion = "check the fit, epsilon dimensions, budget, and forecast tolerances"
    )
  })
  forecasts <- lapply(raw$forecasts, as_mdcev_forecast_data_frame)
  descriptions <- lapply(raw$descriptions, as_mdcev_forecast_data_frame)
  structure(
    list(
      values = forecasts,
      descriptions = descriptions,
      total_budget = total_budget,
      brute_force = isTRUE(brute_force),
      number_of_draws = as.integer(raw$number_of_draws),
      number_of_observations = length(forecasts)
    ),
    class = "biogeme_mdcev_forecast"
  )
}

#' Return native pandas description tables for an MDCEV forecast
#'
#' @param forecast A `biogeme_mdcev_forecast` object.
#' @return A list of ordinary R data frames corresponding to native
#'   pandas `DataFrame.describe()` output.
#' @export
mdcev_forecast_describe <- function(forecast) {
  if (!inherits(forecast, "biogeme_mdcev_forecast")) {
    stop("forecast must be a biogeme_mdcev_forecast object.", call. = FALSE)
  }
  forecast$descriptions
}

#' @export
print.biogeme_mdcev_forecast <- function(x, ...) {
  cat(
    "Native MDCEV forecast for ", x$number_of_observations,
    " observations and ", x$number_of_draws, " draws",
    if (isTRUE(x$brute_force)) " (brute force)\n" else " (analytical)\n",
    sep = ""
  )
  invisible(x)
}

#' @export
as.data.frame.biogeme_mdcev_forecast <- function(x, ...) {
  if (length(x$values) != 1L) {
    stop("A multi-observation MDCEV forecast cannot be represented as one data frame.", call. = FALSE)
  }
  x$values[[1L]]
}

#' Validate the two native MDCEV forecasting algorithms
#'
#' @param model A [biogeme_mdcev_model()] object.
#' @param fit A fresh [mdcev_estimate()] result.
#' @param database Optional database used only for forecast validation.
#' @param total_budget Positive total budget.
#' @param epsilons Optional list of Gumbel-draw matrices.
#' @param number_of_draws Number of native draws when `epsilons` is omitted.
#' @param seed Optional native NumPy seed when draws are generated.
#' @param tolerance_dual Native dual-variable tolerance.
#' @param tolerance_budget Native budget-constraint tolerance.
#' @return `TRUE` when native validation completes without an exception.
#' @export
mdcev_validate_forecast <- function(
    model,
    fit,
    database = NULL,
    total_budget,
    epsilons = NULL,
    number_of_draws = NULL,
    seed = NULL,
    tolerance_dual = 1e-13,
    tolerance_budget = 1e-13
) {
  validate_mdcev_operation_model(model, "MDCEV forecast validation")
  if (!inherits(fit, "biogeme_fit")) {
    stop("fit must be a biogeme_fit returned by mdcev_estimate().", call. = FALSE)
  }
  if (!is.null(database)) validate_biogeme_database(database)
  total_budget <- validate_mdcev_budget(total_budget)
  tolerance_dual <- validate_mdcev_tolerance(tolerance_dual, "tolerance_dual")
  tolerance_budget <- validate_mdcev_tolerance(tolerance_budget, "tolerance_budget")
  if (!is.null(number_of_draws)) {
    number_of_draws <- validate_mdcev_positive_integer(number_of_draws, "number_of_draws")
  }
  if (!is.null(seed)) seed <- validate_mdcev_seed(seed)
  if (is.null(epsilons) && is.null(number_of_draws)) {
    stop("number_of_draws is required when epsilons is NULL.", call. = FALSE)
  }
  if (!is.null(epsilons)) {
    epsilons <- validate_mdcev_epsilons(epsilons, length(model$baseline_utilities))
  }

  compiled <- biogeme_compile_model(model, database_override = database)
  fit_for_bridge <- fit
  fit_for_bridge$model <- NULL
  tryCatch({
    bridge <- biogeme_bridge()
    epsilon_payload <- if (is.null(epsilons)) NULL else reticulate::r_to_py(epsilons)
    result <- bridge$validate_forecast_mdcev_biogeme(
      compiled_model = compiled,
      estimation_results = reticulate::r_to_py(fit_for_bridge),
      total_budget = total_budget,
      epsilons = epsilon_payload,
      number_of_draws = if (is.null(number_of_draws)) NULL else number_of_draws,
      seed = if (is.null(seed)) NULL else seed,
      tolerance_dual = tolerance_dual,
      tolerance_budget = tolerance_budget
    )
    isTRUE(reticulate::py_to_r(result)$validated)
  }, error = function(error) {
    biogeme_rethrow(
      error,
      class = "biogeme_simulation_error",
      operation = "native MDCEV forecast validation",
      suggestion = "check the fit, epsilon dimensions, budget, and forecast tolerances"
    )
  })
}

#' Return native Biogeme's compact MDCEV estimation summary
#'
#' @param fit A `biogeme_fit` returned by [mdcev_estimate()].
#' @return The native Biogeme summary string.
#' @export
mdcev_short_summary <- function(fit) {
  if (!inherits(fit, "biogeme_fit")) {
    stop("fit must be a biogeme_fit returned by mdcev_estimate().", call. = FALSE)
  }
  fit_for_bridge <- fit
  fit_for_bridge$model <- NULL
  tryCatch(
    as.character(reticulate::py_to_r(
      biogeme_bridge()$mdcev_short_summary_biogeme(
        estimation_results = reticulate::r_to_py(fit_for_bridge)
      )
    )),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "native MDCEV short summary",
        suggestion = "check that fit is a complete native Biogeme estimation result"
      )
    }
  )
}

#' Return native Biogeme's estimated-parameter table for an MDCEV fit
#'
#' @param fit A `biogeme_fit` returned by [mdcev_estimate()].
#' @param variance_covariance_type Optional native covariance type, such as
#'   `"BHHH"` or `"Bootstrap"`.
#' @return A named list of ordinary R data frames returned by native Biogeme.
#' @export
mdcev_parameter_table <- function(fit, variance_covariance_type = NULL) {
  if (!inherits(fit, "biogeme_fit")) {
    stop("fit must be a biogeme_fit returned by mdcev_estimate().", call. = FALSE)
  }
  fit_for_bridge <- fit
  fit_for_bridge$model <- NULL
  raw <- tryCatch(
    biogeme_bridge()$mdcev_parameter_table_biogeme(
      estimation_results = reticulate::r_to_py(fit_for_bridge),
      variance_covariance_type = variance_covariance_type
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "native MDCEV parameter table",
        suggestion = "check the covariance type and native estimation result"
      )
    }
  )
  raw <- reticulate::py_to_r(raw)
  lapply(raw, as_mdcev_forecast_data_frame)
}