R/results.R

In rbiogeme: Interface to 'Biogeme'

validate_estimation_controls <- function(controls) {
  if (is.null(controls)) {
    return(list())
  }
  if (is.list(controls) && length(controls) == 0L) {
    return(list())
  }
  if (!is.list(controls) || is.null(names(controls)) || anyNA(names(controls)) ||
      any(!nzchar(names(controls))) || anyDuplicated(names(controls))) {
    stop("controls must be a named list.", call. = FALSE)
  }
  controls
}

# Persistent native files must have an explicit destination. In particular,
# never turn a missing destination into getwd() or a package-directory path.
biogeme_output_control_names <- c(
  "generate_html", "generate_yaml", "generate_netcdf", "generate_pickle",
  "save_iterations"
)

biogeme_output_requested <- function(controls, default_controls = character()) {
  explicit_output <- any(vapply(
    biogeme_output_control_names,
    function(name) isTRUE(controls[[name]]),
    logical(1)
  ))
  default_output <- any(vapply(
    default_controls,
    function(name) !name %in% names(controls) || isTRUE(controls[[name]]),
    logical(1)
  ))
  explicit_output || default_output
}

biogeme_validate_output_directory <- function(output_directory) {
  if (!is.character(output_directory) || length(output_directory) != 1L ||
      is.na(output_directory) || !nzchar(output_directory)) {
    stop(
      "output_directory must be one non-empty path supplied explicitly by the user.",
      call. = FALSE
    )
  }
  output_directory <- path.expand(output_directory)
  if (file.exists(output_directory) && !dir.exists(output_directory)) {
    stop("output_directory must identify a directory: ", output_directory, call. = FALSE)
  }
  dir.create(output_directory, recursive = TRUE, showWarnings = FALSE)
  if (!dir.exists(output_directory)) {
    stop("Cannot create output_directory: ", output_directory, call. = FALSE)
  }
  normalizePath(output_directory, mustWork = TRUE)
}

biogeme_prepare_output_controls <- function(
    controls,
    yaml_file_name = NULL,
    default_controls = character()
) {
  controls <- validate_estimation_controls(controls)
  if (!is.null(yaml_file_name)) controls$generate_yaml <- TRUE
  output_directory <- controls$output_directory
  output_requested <- biogeme_output_requested(controls, default_controls)
  yaml_only <- isTRUE(controls$generate_yaml) &&
    !any(vapply(
      setdiff(biogeme_output_control_names, "generate_yaml"),
      function(name) isTRUE(controls[[name]]),
      logical(1)
    ))
  if (output_requested && is.null(output_directory) &&
      !(yaml_only && !is.null(yaml_file_name))) {
    stop(
      paste0(
        "An explicit output_directory is required whenever native Biogeme ",
        "outputs are requested. Use a directory you selected, such as ",
        "tempdir() in tests and vignettes."
      ),
      call. = FALSE
    )
  }
  if (!is.null(output_directory)) {
    controls$output_directory <- biogeme_validate_output_directory(output_directory)
  }
  controls
}

validate_starting_values <- function(starting_values) {
  if (is.null(starting_values)) {
    return(NULL)
  }
  if (!is.numeric(starting_values) || is.null(names(starting_values)) ||
      anyNA(names(starting_values)) || any(!nzchar(names(starting_values))) ||
      anyNA(starting_values) || any(!is.finite(starting_values))) {
    stop(
      "starting_values must be a named finite numeric vector.",
      call. = FALSE
    )
  }
  as.numeric(starting_values) |> setNames(names(starting_values))
}

as_biogeme_fit <- function(raw_results, model) {
  result <- reticulate::py_to_r(raw_results)
  if (!is.list(result) || is.null(result$beta_names) || is.null(result$beta_values)) {
    biogeme_abort(
      "Biogeme returned an invalid estimation result.",
      class = "biogeme_result_error",
      operation = "estimation result conversion",
      suggestion = "enable debug mode and report the complete result structure"
    )
  }
  result$beta_names <- as.character(result$beta_names)
  result$beta_values <- as.numeric(result$beta_values)
  names(result$beta_values) <- result$beta_names
  for (field in c("lower_bounds", "upper_bounds")) {
    if (!is.null(result[[field]])) {
      result[[field]] <- vapply(
        result[[field]],
        function(value) if (is.null(value)) NA_real_ else as.numeric(value),
        numeric(1)
      )
    }
  }
  for (field in c("standard_errors", "t_statistics", "p_values")) {
    if (!is.null(result[[field]])) {
      result[[field]] <- as.numeric(unlist(result[[field]], use.names = FALSE))
      names(result[[field]]) <- result$beta_names
    }
  }
  if (!is.null(result$variance_covariance)) {
    n_parameters <- length(result$beta_names)
    result$variance_covariance <- matrix(
      as.numeric(unlist(result$variance_covariance, use.names = FALSE)),
      nrow = n_parameters,
      ncol = n_parameters
    )
    dimnames(result$variance_covariance) <- list(result$beta_names, result$beta_names)
  }
  result$model <- model
  structure(result, class = c("biogeme_fit", "biogeme_result"))
}

# Convert the native Bayesian payload into an R-only result object. Posterior
# draws remain in the native NetCDF file; the R object contains the native
# summary and paths, not a live PyMC or ArviZ object.
as_biogeme_bayesian_fit <- function(raw_results, model) {
  result <- raw_results
  if (!is.list(result) || is.null(result$summary) || is.null(result$beta_names)) {
    biogeme_abort(
      "Biogeme returned an invalid Bayesian estimation result.",
      class = "biogeme_result_error",
      operation = "Bayesian estimation result conversion",
      suggestion = "enable debug mode and report the complete result structure"
    )
  }
  result$beta_names <- as.character(result$beta_names)
  result$beta_values <- as.numeric(unlist(result$beta_values, use.names = FALSE))
  names(result$beta_values) <- result$beta_names
  result$posterior_draws <- as.integer(result$posterior_draws)
  result$chains <- as.integer(result$chains)
  result$draws <- as.integer(result$draws)
  result$model <- model
  structure(result, class = c("biogeme_bayesian_fit", "biogeme_result"))
}

#' Estimate a model with native Bayesian inference
#'
#' The complete R expression tree is compiled once, then native Biogeme
#' builds and samples the PyMC model. The returned R object contains the native
#' posterior summary and the paths of any generated NetCDF/YAML/HTML files;
#' posterior draws are not represented as manually managed Python objects.
#'
#' NetCDF and YAML output are enabled by default for this operation. Additional
#' Bayesian controls such as `bayesian_draws`, `warmup`, `chains`, `target_accept`,
#' `calculate_likelihood`, `calculate_waic`, `calculate_loo`, and
#' `mcmc_sampling_strategy` are passed to native Biogeme through `controls`.
#' Because those default Bayesian files are persistent, supply an explicit
#' `output_directory` through [biogeme_control()] (or disable both output types
#' explicitly).
#'
#' @param model A `biogeme_model`.
#' @param model_name Native Biogeme model name.
#' @param controls Named native Biogeme controls.
#' @param starting_values Optional named numeric vector of starting values.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_bayesian_fit` containing the native
#'   posterior summary and output paths.
#' @export
bayesian_estimate <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    starting_values = NULL,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "Bayesian model validation",
      suggestion = "construct the model with biogeme_model()"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_model") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(
    controls,
    default_controls = c("generate_netcdf", "generate_yaml")
  )
  starting_values <- validate_starting_values(starting_values)
  raw_results <- biogeme_bayesian_estimate_model(
    model = model,
    model_name = model_name,
    controls = controls,
    starting_values = starting_values
  )
  as_biogeme_bayesian_fit(raw_results, model)
}

#' Estimate a Biogeme model
#'
#' @param model A `biogeme_model`, including a model created by
#'   [logit_model()] or another specialized constructor.
#' @param model_name Output/model name used by Biogeme.
#' @param controls Named list of Biogeme controls.
#' @param starting_values Optional named numeric vector of starting values.
#' @param run_bootstrap Whether to run bootstrap re-estimation.
#' @param yaml_file_name Optional path for Biogeme's standard YAML output.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_fit`.
#' @details
#' `estimate()` always requests a fresh native estimation. It does not load a
#' previous YAML or iteration file implicitly. Use [estimate_or_load()] when
#' explicit result loading or recycling is required.
#' @examples
#' \dontrun{
#' fit <- estimate(
#'   model,
#'   model_name = "demo",
#'   control = biogeme_control(
#'     output_directory = tempfile("rbiogeme-estimate-"),
#'     generate_html = FALSE,
#'     generate_yaml = FALSE,
#'     save_iterations = FALSE
#'   )
#' )
#' summary(fit)
#' }
#' @export
estimate <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE,
    yaml_file_name = NULL,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "model validation",
      suggestion = "construct the model with biogeme_model(), logit_model(), or another model constructor"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_model") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  if (!is.null(yaml_file_name)) {
    yaml_file_name <- normalizePath(
      path.expand(validate_name(yaml_file_name, "yaml_file_name")),
      mustWork = FALSE
    )
  }
  controls <- biogeme_prepare_output_controls(
    controls,
    yaml_file_name = yaml_file_name
  )
  starting_values <- validate_starting_values(starting_values)
  if (!is.logical(run_bootstrap) || length(run_bootstrap) != 1L || is.na(run_bootstrap)) {
    stop("run_bootstrap must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.null(yaml_file_name)) controls$generate_yaml <- TRUE
  raw_results <- tryCatch(
    biogeme_estimate_model(
      model = model,
      model_name = model_name,
      controls = controls,
      starting_values = starting_values,
      run_bootstrap = isTRUE(run_bootstrap),
      yaml_file_name = yaml_file_name
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme estimation",
        suggestion = "check starting values, controls, data scaling, and the model's numerical conditioning"
      )
    }
  )
  tryCatch(
    as_biogeme_fit(raw_results, model),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "estimation result conversion",
        suggestion = "enable debug mode and report the complete result structure"
      )
    }
  )
}

#' Estimate one named native catalog configuration
#'
#' The configuration identifier is resolved by native
#' `BIOGEME.from_configuration()`.  R only supplies the compiled expression
#' tree and receives ordinary serialized estimation results; no R expression
#' or callback is evaluated during estimation.
#'
#' @param model A `biogeme_model` containing catalog expressions.
#' @param configuration_id Exact native configuration identifier, for example
#'   `"model_catalog:logit;train_tt_catalog:linear"`.
#' @param model_name Native Biogeme model name.
#' @param controls Named list of 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; an alias for `controls`.
#' @return An object of class `biogeme_fit` containing native results.
#' @export
estimate_configuration <- function(
    model,
    configuration_id,
    model_name = "rbiogeme_configuration",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE,
    yaml_file_name = NULL,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "selected configuration validation",
      suggestion = "construct the catalog model with biogeme_model()"
    )
  }
  if (!is.character(configuration_id) || length(configuration_id) != 1L ||
      is.na(configuration_id) || !nzchar(configuration_id)) {
    stop("configuration_id must be one non-empty character string.", call. = FALSE)
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_configuration") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  if (!is.null(yaml_file_name)) {
    yaml_file_name <- normalizePath(
      path.expand(validate_name(yaml_file_name, "yaml_file_name")),
      mustWork = FALSE
    )
  }
  controls <- biogeme_prepare_output_controls(
    controls,
    yaml_file_name = yaml_file_name
  )
  starting_values <- validate_starting_values(starting_values)
  if (!is.logical(run_bootstrap) || length(run_bootstrap) != 1L || is.na(run_bootstrap)) {
    stop("run_bootstrap must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.null(yaml_file_name)) controls$generate_yaml <- TRUE
  raw_results <- biogeme_estimate_configuration_model(
    model = model,
    configuration_id = configuration_id,
    model_name = model_name,
    controls = controls,
    starting_values = starting_values,
    run_bootstrap = isTRUE(run_bootstrap),
    yaml_file_name = yaml_file_name
  )
  tryCatch(
    as_biogeme_fit(raw_results, model),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "selected configuration result conversion",
        suggestion = "enable debug mode and report the complete native result structure"
      )
    }
  )
}

# Convert the serialized native pandas table returned by the catalog bridge.
as_biogeme_catalog_table <- 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 catalog summary table.", call. = FALSE)
  }
  index <- as.character(value$index)
  columns <- as.character(value$columns)
  rows <- if (is.null(value$data)) list() else value$data
  if (length(rows) == 0L) {
    table <- as.data.frame(
      matrix(nrow = 0L, ncol = length(columns)),
      stringsAsFactors = FALSE
    )
  } else {
    table <- as.data.frame(
      do.call(rbind, lapply(rows, function(row) as.character(unlist(row)))),
      stringsAsFactors = FALSE,
      check.names = FALSE
    )
  }
  if (length(columns) > 0L) names(table) <- columns
  rownames(table) <- index
  table
}

#' Estimate every specification in a native Biogeme catalog
#'
#' Catalog choices, synchronized controllers, estimation, result summaries,
#' Pareto selection, and parameter comparison are all handled by native
#' Biogeme. The returned object only contains ordinary R representations of
#' those native results.
#'
#' @param model A `biogeme_model` containing one or more catalog expressions.
#' @param model_name Native Biogeme model name prefix.
#' @param controls Named list of native Biogeme controls.
#' @param quick_estimate Whether to use native quick estimation.
#' @param run_bootstrap Whether to run native bootstrap re-estimation.
#' @param force Whether to estimate afresh rather than recycle native files.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_catalog_fit` containing native model
#'   results, summaries, Pareto-optimal specifications, and LaTeX comparison.
#' @export
estimate_catalog <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list(),
    quick_estimate = FALSE,
    run_bootstrap = FALSE,
    force = TRUE,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "catalog model validation",
      suggestion = "construct the model with biogeme_model()"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.logical(quick_estimate) || length(quick_estimate) != 1L || is.na(quick_estimate)) {
    stop("quick_estimate must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.logical(run_bootstrap) || length(run_bootstrap) != 1L || is.na(run_bootstrap)) {
    stop("run_bootstrap must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.logical(force) || length(force) != 1L || is.na(force)) {
    stop("force must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_catalog") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(controls)
  raw <- tryCatch(
    biogeme_estimate_catalog_model(
      model = model,
      model_name = model_name,
      controls = controls,
      quick_estimate = quick_estimate,
      run_bootstrap = run_bootstrap,
      force = force
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme catalog estimation",
        suggestion = "check the catalog controllers, model expressions, controls, and data scaling"
      )
    }
  )
  raw <- reticulate::py_to_r(raw)
  converted_results <- lapply(raw$results, as_biogeme_fit, model = model)
  structure(
    list(
      results = converted_results,
      summary = as_biogeme_catalog_table(raw$summary),
      description = unlist(raw$description, use.names = TRUE),
      non_dominated = as.character(raw$non_dominated),
      non_dominated_summary = as_biogeme_catalog_table(raw$non_dominated_summary),
      non_dominated_description = unlist(raw$non_dominated_description, use.names = TRUE),
      latex = as.character(raw$latex),
      model = model
    ),
    class = c("biogeme_catalog_fit", "biogeme_result")
  )
}

#' Run native Biogeme assisted specification
#'
#' The complete catalog search, quick-estimation objective evaluation, Pareto
#' persistence, and final re-estimation are delegated to native Biogeme. The
#' R interface accepts named objective and validity descriptors rather than
#' exposing Python callback objects.  The selected validity rule is created and
#' executed inside the Python bridge between native estimates.
#'
#' @param model A biogeme_model containing catalog expressions.
#' @param objectives Native objective preset. Supported values are
#'   "loglikelihood_dimension" and "aic_bic_dimension".
#' @param validity Optional native validity-rule name. Currently supported is
#'   `"negative_time_cost"`, matching the predicate in the native assisted
#'   Swissmetro example.
#' @param pareto_file_name Explicit path of the native Pareto checkpoint file.
#'   It has no default because native checkpoints are persistent files.
#' @param model_name Native Biogeme model name prefix.
#' @param controls Named list of native Biogeme controls.
#' @param force Whether to remove the named Pareto checkpoint and start fresh.
#' @param control Optional biogeme_control() object; an alias for controls.
#' @return An object of class biogeme_assisted_fit containing native final
#'   results, the summary table, descriptions, and Pareto metadata.
#' @export
assisted_specification <- function(
    model,
    objectives = "loglikelihood_dimension",
    pareto_file_name = NULL,
    model_name = "rbiogeme_assisted",
    controls = list(),
    force = TRUE,
    control = NULL,
    validity = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "assisted specification validation",
      suggestion = "construct the model with biogeme_model() and catalog expressions"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.character(objectives) || length(objectives) != 1L ||
      is.na(objectives) || !nzchar(objectives)) {
    stop("objectives must be one non-empty character string.", call. = FALSE)
  }
  if (!is.null(validity) && (!is.character(validity) || length(validity) != 1L ||
      is.na(validity) || !nzchar(validity))) {
    stop("validity must be NULL or one non-empty character string.", call. = FALSE)
  }
  if (is.null(pareto_file_name)) {
    stop(
      paste0(
        "pareto_file_name must be supplied explicitly; persistent native ",
        "outputs are never written to the working directory by default."
      ),
      call. = FALSE
    )
  }
  if (!is.character(pareto_file_name) ||
      length(pareto_file_name) != 1L ||
      is.na(pareto_file_name) || !nzchar(pareto_file_name)) {
    stop("pareto_file_name must be one non-empty character string.", call. = FALSE)
  }
  if (!is.logical(force) || length(force) != 1L || is.na(force)) {
    stop("force must be one non-missing logical value.", call. = FALSE)
  }
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_assisted") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(controls)
  pareto_file_name <- normalizePath(path.expand(pareto_file_name), mustWork = FALSE)
  raw <- tryCatch(
    biogeme_assisted_specification_model(
      model = model,
      model_name = model_name,
      controls = controls,
      objectives = objectives,
      pareto_file_name = pareto_file_name,
      force = force,
      validity = validity
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme assisted specification",
        suggestion = "check the catalog expressions, objective preset, controls, and Pareto path"
      )
    }
  )
  raw <- reticulate::py_to_r(raw)
  converted_results <- lapply(raw$results, as_biogeme_fit, model = model)
  structure(
    list(
      results = converted_results,
      summary = as_biogeme_catalog_table(raw$summary),
      description = unlist(raw$description, use.names = TRUE),
      pareto_file_name = raw$pareto_file_name,
      pareto_statistics = as.character(raw$pareto_statistics),
      model = model
    ),
    class = c("biogeme_assisted_fit", "biogeme_result")
  )
}

#' @rdname assisted_specification
#' @export
run_assisted_specification <- assisted_specification

#' Count catalog specifications using native Biogeme
#'
#' @param model A `biogeme_model` containing catalog expressions.
#' @param model_name Native Biogeme model name.
#' @param controls Named list of native Biogeme controls.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return The native number of catalog configurations, or `NULL` when native
#'   Biogeme cannot determine a count.
#' @export
count_number_of_specifications <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list(),
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "catalog specification count validation",
      suggestion = "construct the model with biogeme_model() and catalog expressions"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_catalog") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- validate_estimation_controls(controls)
  raw <- tryCatch(
    biogeme_count_number_of_specifications_model(
      model = model,
      model_name = model_name,
      controls = controls
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "native catalog specification count",
        suggestion = "check the catalog controllers and model expression"
      )
    }
  )
  value <- reticulate::py_to_r(raw)
  if (is.null(value)) NULL else as.integer(value)
}

#' Return catalog configuration identifiers using native Biogeme
#'
#' @param model A `biogeme_model` containing catalog expressions.
#' @param model_name Native Biogeme model name.
#' @param controls Named list of native Biogeme controls.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return Character vector of exact native configuration identifiers.
#' @export
catalog_configuration_ids <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list(),
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "catalog configuration discovery validation",
      suggestion = "construct a model with biogeme_model() and catalog expressions"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_catalog") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- validate_estimation_controls(controls)
  raw <- tryCatch(
    biogeme_catalog_configuration_ids_model(
      model = model,
      model_name = model_name,
      controls = controls
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "native catalog configuration discovery",
        suggestion = "check the catalog controllers and model expression"
      )
    }
  )
  value <- reticulate::py_to_r(raw)
  if (is.null(value)) character() else as.character(value)
}

#' Re-estimate the Pareto-optimal models saved by native Biogeme
#'
#' The Pareto file is read by native `ParetoPostProcessing`. Re-estimation,
#' summary compilation, model descriptions, and optional Pareto plotting all
#' remain native operations; the returned object contains ordinary R values.
#'
#' @param model A `biogeme_model` containing the catalog expressions used to
#'   create the Pareto file.
#' @param pareto_file_name Path to an existing native Pareto file.
#' @param model_name Native Biogeme model-name prefix for re-estimation.
#' @param controls Named list of native Biogeme controls.
#' @param recycle Whether native Biogeme may recycle complete estimation files.
#'   The default is `FALSE`, matching the Python example.
#' @param plot_file_name Optional path where native Pareto plot output is saved.
#' @param objective_x Zero-based native Pareto objective index for the plot.
#' @param objective_y Zero-based native Pareto objective index for the plot.
#' @param label_x Optional native plot x-axis label.
#' @param label_y Optional native plot y-axis label.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_pareto_fit` containing native results,
#'   the compiled summary, descriptions, Pareto statistics, and plot path.
#' @export
pareto_post_processing <- function(
    model,
    pareto_file_name,
    model_name = "rbiogeme_pareto",
    controls = list(),
    recycle = FALSE,
    plot_file_name = NULL,
    objective_x = 0L,
    objective_y = 1L,
    label_x = NULL,
    label_y = NULL,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "Pareto post-processing validation",
      suggestion = "construct the catalog model used to create the Pareto file"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.character(pareto_file_name) || length(pareto_file_name) != 1L ||
      is.na(pareto_file_name) || !nzchar(pareto_file_name)) {
    stop("pareto_file_name must be one non-empty character string.", call. = FALSE)
  }
  pareto_file_name <- normalizePath(path.expand(pareto_file_name), mustWork = TRUE)
  if (!file.exists(pareto_file_name) || dir.exists(pareto_file_name)) {
    stop("pareto_file_name must identify an existing Pareto file.", call. = FALSE)
  }
  if (!is.logical(recycle) || length(recycle) != 1L || is.na(recycle)) {
    stop("recycle must be one non-missing logical value.", call. = FALSE)
  }
  for (argument in list(objective_x = objective_x, objective_y = objective_y)) {
    value <- argument[[1L]]
    if (!is.numeric(value) || length(value) != 1L || is.na(value) ||
        !is.finite(value) || value < 0 || value != floor(value)) {
      stop("Pareto objective indices must be non-negative integers.", call. = FALSE)
    }
  }
  validate_plot_label <- function(value, name) {
    if (!is.null(value) && (!is.character(value) || length(value) != 1L ||
        is.na(value) || !nzchar(value))) {
      stop(name, " must be NULL or one non-empty character string.", call. = FALSE)
    }
    value
  }
  label_x <- validate_plot_label(label_x, "label_x")
  label_y <- validate_plot_label(label_y, "label_y")
  if (!is.null(plot_file_name)) {
    if (!is.character(plot_file_name) || length(plot_file_name) != 1L ||
        is.na(plot_file_name) || !nzchar(plot_file_name)) {
      stop("plot_file_name must be NULL or one non-empty character string.", call. = FALSE)
    }
    plot_file_name <- normalizePath(path.expand(plot_file_name), mustWork = FALSE)
  }
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_pareto") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(controls)
  raw <- tryCatch(
    biogeme_pareto_post_processing_model(
      model = model,
      model_name = model_name,
      controls = controls,
      pareto_file_name = pareto_file_name,
      recycle = recycle,
      plot_file_name = plot_file_name,
      objective_x = as.integer(objective_x),
      objective_y = as.integer(objective_y),
      label_x = label_x,
      label_y = label_y
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme Pareto post-processing",
        suggestion = "check the Pareto file, catalog expressions, controls, and plot options"
      )
    }
  )
  raw <- reticulate::py_to_r(raw)
  converted_results <- lapply(raw$results, as_biogeme_fit, model = model)
  structure(
    list(
      results = converted_results,
      summary = as_biogeme_catalog_table(raw$summary),
      description = unlist(raw$description, use.names = TRUE),
      pareto_file_name = raw$pareto_file_name,
      pareto_statistics = as.character(raw$pareto_statistics),
      plot_file_name = raw$plot_file_name,
      model = model
    ),
    class = c("biogeme_pareto_fit", "biogeme_result")
  )
}

#' Check analytical derivatives against native finite differences
#'
#' This delegates to the public `BIOGEME.check_derivatives()` operation. The
#' complete expression tree is compiled first; no R callback is evaluated by
#' the derivative checker.
#'
#' @param model A `biogeme_model`.
#' @param model_name Native Biogeme model name.
#' @param controls Named native Biogeme controls.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @param verbose Whether native Biogeme should print the comparison.
#' @return A list containing the native function value, analytical and finite
#'   difference gradients and Hessians, and their error vectors.
#' @export
check_derivatives <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    control = NULL,
    verbose = FALSE
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "derivative-check model validation",
      suggestion = "construct the model with logit_model() or biogeme_model()"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.logical(verbose) || length(verbose) != 1L || is.na(verbose)) {
    stop("verbose must be one non-missing logical value.", call. = FALSE)
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- validate_estimation_controls(controls)
  biogeme_check_derivatives_model(
    model = model,
    model_name = model_name,
    controls = controls,
    verbose = verbose
  )
}

#' Profile native JAX formula evaluation
#'
#' The complete expression tree is compiled once, then native Biogeme's
#' `CompiledFormulaEvaluator` is used for each requested evaluation mode. Each
#' mode is evaluated twice, matching the profiling example's first-call and
#' steady-state measurements. R receives only ordinary environment and timing
#' summaries; no Python evaluator or profiler object is exposed.
#'
#' @param model A `biogeme_model`.
#' @param beta_values Named finite numeric vector of parameter values. If
#'   `NULL`, the model's declared starting values are used.
#' @param cases A list of named case lists. Each case must contain a non-empty
#'   `label` and optional logical `gradient`, `hessian`, and `bhhh` flags.
#' @param model_name Native Biogeme model name.
#' @param controls Named native Biogeme controls.
#' @param numerically_safe Whether native formula evaluation should use
#'   Biogeme's numerical-safety transformations.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_jax_profile` containing the native JAX
#'   environment and one serialized profile for each case.
#' @export
profile_jax <- function(
    model,
    beta_values = NULL,
    cases = list(
      list(label = "Function only", gradient = FALSE, hessian = FALSE, bhhh = FALSE),
      list(label = "Function + gradient", gradient = TRUE, hessian = FALSE, bhhh = FALSE),
      list(label = "Function + gradient + Hessian", gradient = TRUE, hessian = TRUE, bhhh = FALSE),
      list(label = "Function + gradient + BHHH", gradient = TRUE, hessian = FALSE, bhhh = TRUE)
    ),
    model_name = "rbiogeme_model",
    controls = list(),
    numerically_safe = FALSE,
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "JAX profiling model validation",
      suggestion = "construct the model with biogeme_model()"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  controls <- validate_estimation_controls(controls)
  if (is.null(beta_values)) {
    beta_values <- vapply(model$parameters, function(parameter) parameter$start, numeric(1))
    names(beta_values) <- names(model$parameters)
  }
  if (!is.numeric(beta_values) || is.null(names(beta_values)) ||
      length(beta_values) == 0L || anyNA(names(beta_values)) ||
      any(!nzchar(names(beta_values))) || anyDuplicated(names(beta_values)) ||
      anyNA(beta_values) || any(!is.finite(beta_values))) {
    stop("beta_values must be a named non-empty finite numeric vector.", call. = FALSE)
  }
  if (!is.list(cases) || length(cases) == 0L) {
    stop("cases must be a non-empty list.", call. = FALSE)
  }
  cases <- lapply(unname(cases), function(case) {
    if (!is.list(case) || is.null(case$label) ||
        !is.character(case$label) || length(case$label) != 1L ||
        is.na(case$label) || !nzchar(case$label)) {
      stop("Each profiling case must have one non-empty character label.", call. = FALSE)
    }
    flags <- lapply(c("gradient", "hessian", "bhhh"), function(flag) {
      value <- case[[flag]]
      if (is.null(value)) value <- FALSE
      if (!is.logical(value) || length(value) != 1L || is.na(value)) {
        stop("Profiling case flags must be one non-missing logical value.", call. = FALSE)
      }
      isTRUE(value)
    })
    names(flags) <- c("gradient", "hessian", "bhhh")
    c(list(label = case$label), flags)
  })
  if (!is.logical(numerically_safe) || length(numerically_safe) != 1L ||
      is.na(numerically_safe)) {
    stop("numerically_safe must be one non-missing logical value.", call. = FALSE)
  }
  if (identical(model_name, "rbiogeme_model") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  raw <- biogeme_profile_jax_model(
    model = model,
    model_name = model_name,
    controls = controls,
    beta_values = beta_values,
    cases = cases,
    numerically_safe = numerically_safe
  )
  if (!is.list(raw) || is.null(raw$environment) || is.null(raw$cases)) {
    biogeme_abort(
      "Biogeme returned an invalid JAX profiling result.",
      class = "biogeme_result_error",
      operation = "JAX profiling result conversion",
      suggestion = "enable debug mode and report the complete result structure"
    )
  }
  structure(
    list(
      environment = raw$environment,
      cases = raw$cases,
      beta_values = beta_values,
      model = model
    ),
    class = c("biogeme_jax_profile", "biogeme_result")
  )
}

#' Run native post-estimation Monte Carlo draw-stability diagnostics
#'
#' The model is compiled once and the supplied fixed estimate is passed to
#' native `BIOGEME.check_monte_carlo_stability()`. Native Biogeme evaluates the
#' objective and gradient at fresh draw designs and writes its YAML checkpoint
#' and Markdown report.
#'
#' @param model A `biogeme_model`.
#' @param fit A `biogeme_fit` obtained from the same model.
#' @param model_name Native Biogeme model name.
#' @param controls Named native Biogeme controls.
#' @param output_directory Explicit directory for the native diagnostic files.
#'   It has no default: the operation refuses to write into the working
#'   directory when this argument is omitted.
#' @param basename Optional native diagnostic filename prefix.
#' @param resume Whether native Biogeme may resume a compatible checkpoint.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_monte_carlo_diagnostic` containing the
#'   native status, conclusion, recommendation, data, and output paths.
#' @export
check_monte_carlo_stability <- function(
    model,
    fit,
    model_name = "rbiogeme_model",
    controls = list(),
    output_directory = NULL,
    basename = NULL,
    resume = TRUE,
    control = NULL
) {
  if (!inherits(model, "biogeme_model") || !inherits(fit, "biogeme_fit")) {
    biogeme_abort(
      "model must be a biogeme_model and fit must be a biogeme_fit.",
      class = "biogeme_specification_error",
      operation = "Monte Carlo diagnostic validation",
      suggestion = "estimate the model with estimate() and pass that fit"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  output_directory <- biogeme_validate_output_directory(output_directory)
  if (!is.null(basename) && (!is.character(basename) || length(basename) != 1L ||
      is.na(basename) || !nzchar(basename))) {
    stop("basename must be NULL or one non-empty character string.", call. = FALSE)
  }
  if (!is.logical(resume) || length(resume) != 1L || is.na(resume)) {
    stop("resume must be one non-missing logical value.", call. = FALSE)
  }
  if (identical(model_name, "rbiogeme_model") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(controls)
  raw <- tryCatch(
    biogeme_check_monte_carlo_stability_model(
      model = model,
      fit = fit,
      model_name = model_name,
      controls = controls,
      output_directory = output_directory,
      basename = basename,
      resume = resume
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme Monte Carlo stability diagnostic",
        suggestion = "check the fitted model, diagnostic controls, and output directory"
      )
    }
  )
  raw <- reticulate::py_to_r(raw)
  structure(
    list(
      execution_status = as.character(raw$execution_status),
      diagnostic_conclusion = as.character(raw$diagnostic_conclusion),
      recommendation = as.character(raw$recommendation),
      data = raw$data,
      yaml_file = as.character(raw$yaml_file),
      markdown_file = as.character(raw$markdown_file),
      model = model,
      fit = fit
    ),
    class = c("biogeme_monte_carlo_diagnostic", "biogeme_result")
  )
}

#' Estimate a Biogeme model using the native quick-estimation operation
#'
#' `quick_estimate()` delegates to native `BIOGEME.quick_estimate()`. It
#' returns parameter values and the statistics that native Biogeme makes
#' available without the ordinary post-estimation derivative calculations.
#' Use [save_results()] when a standard YAML result is required.
#'
#' @param model A `biogeme_model`.
#' @param model_name Output/model name used by Biogeme.
#' @param controls Named list of Biogeme controls.
#' @param control Optional [biogeme_control()] object; an alias for `controls`.
#' @return An object of class `biogeme_fit`.
#' @export
quick_estimate <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    control = NULL
) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "model validation",
      suggestion = "construct the model with biogeme_model(), logit_model(), or another model constructor"
    )
  }
  if (!is.null(control)) {
    if (!is.null(controls) && length(controls) > 0L) {
      stop("Provide either controls or control, not both.", call. = FALSE)
    }
    controls <- control
  } else if (length(controls) == 0L && !is.null(model$control)) {
    controls <- model$control
  }
  controls <- if (is.null(controls)) list() else controls
  if (!is.list(controls)) stop("controls must be a named list.", call. = FALSE)
  if (identical(model_name, "rbiogeme_model") && !is.null(controls$model_name)) {
    model_name <- controls$model_name
  }
  model_name <- validate_name(model_name, "model_name")
  controls <- biogeme_prepare_output_controls(controls)
  raw_results <- tryCatch(
    biogeme_quick_estimate_model(
      model = model,
      model_name = model_name,
      controls = controls
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme quick estimation",
        suggestion = "check starting values, controls, data scaling, and the model's numerical conditioning"
      )
    }
  )
  tryCatch(
    as_biogeme_fit(raw_results, model),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "quick-estimation result conversion",
        suggestion = "enable debug mode and report the complete result structure"
      )
    }
  )
}

#' Estimate or explicitly load a standard Biogeme result
#'
#' `force = FALSE` permits loading the named YAML file; `force = TRUE` always
#' performs fresh estimation. The default `estimate()` path never recycles an
#' old result file.
#' @param model A `biogeme_model`.
#' @param yaml_file_name YAML path used for loading or saving.
#' @param force Whether to force fresh estimation.
#' @param model_name Native Biogeme model name.
#' @param controls Named native controls.
#' @param starting_values Optional named starting values.
#' @param run_bootstrap Whether to run bootstrap estimation.
#' @return A `biogeme_fit` object.
#' @export
estimate_or_load <- function(
    model,
    yaml_file_name,
    force = FALSE,
    model_name = "rbiogeme_model",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE
) {
  if (!inherits(model, "biogeme_model")) stop("model must be a Biogeme model.", call. = FALSE)
  yaml_file_name <- normalizePath(
    path.expand(validate_name(yaml_file_name, "yaml_file_name")),
    mustWork = FALSE
  )
  if (!is.logical(force) || length(force) != 1L || is.na(force)) {
    stop("force must be one non-missing logical value.", call. = FALSE)
  }
  controls <- biogeme_prepare_output_controls(
    controls,
    yaml_file_name = yaml_file_name
  )
  model_name <- validate_name(model_name, "model_name")
  compiled <- biogeme_compile_model(model)
  raw <- tryCatch(
    biogeme_bridge()$estimate_or_load_biogeme(
      compiled_model = compiled,
      model_name = model_name,
      controls = if (length(controls) == 0L) NULL else reticulate::r_to_py(controls),
      starting_values = if (is.null(starting_values)) NULL else reticulate::r_to_py(starting_values),
      run_bootstrap = isTRUE(run_bootstrap),
      yaml_file_name = yaml_file_name,
      force = isTRUE(force)
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme estimate-or-load",
        suggestion = "set force = TRUE for fresh estimation or provide a valid YAML result"
      )
    }
  )
  as_biogeme_fit(raw, model)
}

#' @export
print.biogeme_fit <- function(x, ...) {
  cat("Biogeme estimation results for model ", x$model_name, "\n", sep = "")
  cat("Final log likelihood: ", format(x$final_log_likelihood), "\n", sep = "")
  cat("Converged: ", isTRUE(x$convergence), "\n", sep = "")
  if (!is.null(x$termination_reason)) {
    cat("Termination: ", x$termination_reason, "\n", sep = "")
  }
  if (!is.null(x$beta_values)) {
    print(coef(x))
  }
  invisible(x)
}

#' @export
print.biogeme_catalog_fit <- function(x, ...) {
  cat("Biogeme catalog estimation results: ", length(x$results), " models\n", sep = "")
  print(x$summary)
  invisible(x)
}

#' @export
print.biogeme_assisted_fit <- function(x, ...) {
  cat("Biogeme assisted specification results: ", length(x$results), " models\n", sep = "")
  print(x$summary)
  invisible(x)
}

#' @export
print.biogeme_pareto_fit <- function(x, ...) {
  cat("Biogeme Pareto post-processing results: ", length(x$results), " models\n", sep = "")
  if (length(x$pareto_statistics) > 0L) {
    cat(paste(x$pareto_statistics, collapse = "\n"), "\n", sep = "")
  }
  print(x$summary)
  invisible(x)
}

#' @export
print.biogeme_monte_carlo_diagnostic <- function(x, ...) {
  cat("Biogeme Monte Carlo diagnostic: ", x$execution_status, "\n", sep = "")
  cat("Conclusion: ", x$diagnostic_conclusion, "\n", sep = "")
  cat("Recommendation: ", x$recommendation, "\n", sep = "")
  cat("YAML checkpoint: ", x$yaml_file, "\n", sep = "")
  cat("Markdown report: ", x$markdown_file, "\n", sep = "")
  invisible(x)
}

#' @export
print.biogeme_jax_profile <- function(x, ...) {
  cat("Native JAX profile results\n")
  if (length(x$environment) > 0L) {
    cat("JAX environment:\n")
    for (name in names(x$environment)) {
      value <- x$environment[[name]]
      if (is.null(value)) value <- "<unset>"
      cat("  ", name, ": ", paste(value, collapse = ", "), "\n", sep = "")
    }
  }
  for (case in x$cases) {
    cat("\n=== ", case$label, " ===\n", sep = "")
    cat(case$summary, "\n", sep = "")
  }
  invisible(x)
}

biogeme_coefficient_table <- function(object) {
  values <- coef(object)
  standard_errors <- object$standard_errors
  t_statistics <- object$t_statistics
  p_values <- object$p_values
  if (is.null(standard_errors)) {
    standard_errors <- rep(NA_real_, length(values))
    t_statistics <- standard_errors
    p_values <- standard_errors
  }
  data.frame(
    Estimate = unname(values),
    `Std. Error` = unname(standard_errors),
    `t value` = unname(t_statistics),
    `Pr(>|t|)` = unname(p_values),
    row.names = names(values),
    check.names = FALSE
  )
}

#' @export
summary.biogeme_fit <- function(object, ...) {
  structure(
    list(
      coefficients = biogeme_coefficient_table(object),
      model_name = object$model_name,
      final_log_likelihood = object$final_log_likelihood,
      convergence = object$convergence,
      termination_reason = object$termination_reason,
      number_of_iterations = object$number_of_iterations,
      optimization_time = object$optimization_time,
      optimization_messages = object$optimization_messages,
      biogeme_version = object$biogeme_version,
      python_version = object$python_version,
      derivatives_available = isTRUE(object$derivatives_available) &&
        !is.null(object$variance_covariance)
    ),
    class = "summary.biogeme_fit"
  )
}

#' @export
print.summary.biogeme_fit <- function(x, ...) {
  cat("Biogeme estimation results for model ", x$model_name, "\n", sep = "")
  cat("Final log likelihood: ", format(x$final_log_likelihood), "\n", sep = "")
  cat("Converged: ", isTRUE(x$convergence), "\n\n", sep = "")
  if (!is.null(x$termination_reason)) {
    cat("Termination: ", x$termination_reason, "\n\n", sep = "")
  }
  print(x$coefficients)
  if (!isTRUE(x$derivatives_available)) {
    cat("Covariance and coefficient tests are unavailable because derivatives were not calculated.\n")
  }
  invisible(x)
}

#' @export
print.biogeme_bayesian_fit <- function(x, ...) {
  print(summary(x))
  invisible(x)
}

#' @export
summary.biogeme_bayesian_fit <- function(object, ...) {
  result <- object$summary
  class(result) <- c("summary.biogeme_bayesian_fit", "list")
  result
}

#' @export
print.summary.biogeme_bayesian_fit <- function(x, ...) {
  cat("Bayesian estimation results for model ", x$model_name, "\n", sep = "")
  cat("Chains: ", x$chains, "; draws per chain: ", x$draws, "\n", sep = "")
  if (!is.null(x$parameters) && length(x$parameters) > 0L) {
    parameters <- do.call(
      rbind,
      lapply(x$parameters, function(parameter) {
        data.frame(
          Mean = parameter$mean,
          Median = parameter$median,
          `Std. Error` = parameter$std_err,
          `HDI low` = parameter$hdi_low,
          `HDI high` = parameter$hdi_high,
          check.names = FALSE
        )
      })
    )
    rownames(parameters) <- names(x$parameters)
    print(parameters)
  }
  invisible(x)
}

#' Report variables stored in native Bayesian results
#'
#' The report is serialized by the Python bridge from native Biogeme's
#' `BayesianResultsSummary.report_stored_variables()` method. R receives only
#' an ordinary data frame; posterior and ArviZ objects are not exposed.
#' @param x A `biogeme_bayesian_fit` object.
#' @return A data frame with the native result group, variable, dimensions, and
#'   shape for each stored variable.
#' @export
bayesian_stored_variables <- function(x) {
  if (!inherits(x, "biogeme_bayesian_fit")) {
    stop("x must be a biogeme_bayesian_fit object.", call. = FALSE)
  }
  report <- x$summary$stored_variables_report
  empty <- data.frame(
    group = character(),
    variable = character(),
    dims = character(),
    shape = character(),
    stringsAsFactors = FALSE
  )
  if (is.null(report) || length(report) == 0L) return(empty)
  as_text <- function(value) {
    if (is.null(value)) return(NA_character_)
    if (is.list(value)) value <- unlist(value, use.names = FALSE)
    paste(as.character(value), collapse = ", ")
  }
  rows <- lapply(report, function(item) {
    if (!is.list(item)) stop("stored variable report entries must be lists.", call. = FALSE)
    data.frame(
      group = as_text(item$group),
      variable = as_text(item$variable),
      dims = as_text(item$dims),
      shape = as_text(item$shape),
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

#' @export
coef.biogeme_fit <- function(object, ...) {
  values <- object$beta_values
  names(values) <- object$beta_names
  values
}

#' @export
coef.biogeme_bayesian_fit <- function(object, ...) {
  values <- object$beta_values
  names(values) <- object$beta_names
  values
}

#' @export
vcov.biogeme_fit <- function(object, ...) {
  if (is.null(object$variance_covariance)) {
    biogeme_abort(
      "The variance-covariance matrix is unavailable because derivatives were not calculated.",
      class = "biogeme_result_error",
      operation = "variance-covariance extraction",
      suggestion = "estimate with derivative calculation enabled"
    )
  }
  object$variance_covariance
}

#' @export
logLik.biogeme_fit <- function(object, ...) {
  structure(
    as.numeric(object$final_log_likelihood),
    df = length(object$beta_names),
    nobs = object$sample_size,
    class = "logLik"
  )
}

#' @export
nobs.biogeme_fit <- function(object, use.fallback = TRUE, ...) {
  as.integer(object$sample_size)
}

#' @export
nobs.biogeme_bayesian_fit <- function(object, use.fallback = TRUE, ...) {
  as.integer(object$summary$number_of_observations)
}

#' Save standard Biogeme YAML results
#'
#' @param fit A `biogeme_fit` object.
#' @param filename Destination YAML file.
#' @return `fit`, invisibly.
#' @export
save_results <- function(fit, filename) {
  if (!inherits(fit, "biogeme_fit")) {
    biogeme_abort(
      "fit must be a biogeme_fit object.",
      class = "biogeme_result_error",
      operation = "saving Biogeme results",
      suggestion = "pass the object returned by estimate()"
    )
  }
  filename <- validate_name(filename, "filename")
  raw <- fit
  raw$model <- NULL
  tryCatch(
    biogeme_bridge()$save_results_yaml(
      result = reticulate::r_to_py(raw),
      filename = filename
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "saving Biogeme results",
        suggestion = "check that the result object is complete and that the destination is writable"
      )
    }
  )
  invisible(fit)
}

#' Return native Biogeme general estimation statistics
#'
#' @param fit A `biogeme_fit` object.
#' @return A named list of native Biogeme statistics.
#' @export
biogeme_general_statistics <- function(fit) {
  if (!inherits(fit, "biogeme_fit")) {
    biogeme_abort(
      "fit must be a biogeme_fit object.",
      class = "biogeme_result_error",
      operation = "reading Biogeme general statistics",
      suggestion = "pass the object returned by estimate() or quick_estimate()"
    )
  }
  if (is.null(fit$general_statistics)) {
    biogeme_abort(
      "General statistics are unavailable in this result.",
      class = "biogeme_result_error",
      operation = "reading Biogeme general statistics",
      suggestion = "run the model through the native Biogeme bridge"
    )
  }
  fit$general_statistics
}

#' Read standard Biogeme YAML results
#'
#' @param filename YAML file generated by Biogeme.
#' @return A `biogeme_fit` object.
#' @export
read_results <- function(filename) {
  filename <- validate_name(filename, "filename")
  if (!file.exists(filename)) {
    biogeme_abort(
      paste0("Results file does not exist: ", filename),
      class = "biogeme_result_error",
      operation = "reading Biogeme results",
      suggestion = "provide the path to an existing Biogeme YAML file"
    )
  }
  tryCatch(
    {
      raw_results <- biogeme_bridge()$read_results_yaml(filename)
      as_biogeme_fit(raw_results, model = NULL)
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_result_error",
        operation = "reading Biogeme results",
        suggestion = "check that the YAML file was generated by Biogeme and is not truncated"
      )
    }
  )
}