R/bridge.R

In rbiogeme: Interface to 'Biogeme'

#' Locate and import the package's Python bridge.
#'
#' This is intentionally an internal helper.  Users interact with the
#' high-level estimation functions; keeping the import in one place makes the
#' Python boundary easy to test and keeps all numerical work on the Python side.
#' @keywords internal
biogeme_bridge <- function() {
  biogeme_python()
  bridge_path <- system.file("python", package = "rbiogeme")
  if (!nzchar(bridge_path)) {
    # This fallback is useful while developing the package from a source tree
    # (an installed package normally takes the first branch).
    bridge_path <- file.path(getwd(), "inst", "python")
  }
  if (!dir.exists(bridge_path)) {
    biogeme_abort(
      "The rbiogeme Python bridge is not installed.",
      class = "biogeme_environment_error",
      operation = "rbiogeme Python bridge loading",
      suggestion = "reinstall the package so that inst/python/biogeme_bridge.py is included"
    )
  }
  tryCatch(
    reticulate::import_from_path(
      module = "biogeme_bridge",
      path = bridge_path,
      convert = FALSE
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_environment_error",
        operation = "rbiogeme Python bridge loading",
        suggestion = "reinstall rbiogeme so that inst/python/biogeme_bridge.py is included"
      )
    }
  )
}

biogeme_expression_ir <- function(expression) {
  expression <- as_biogeme_expression(expression)
  kind <- expression$kind
  if (identical(kind, "numeric")) {
    return(list(kind = "numeric", value = as.numeric(expression$value)))
  }
  if (identical(kind, "variable")) {
    return(list(kind = "variable", name = expression$name))
  }
  if (identical(kind, "beta")) {
    return(list(
      kind = "beta",
      name = expression$name,
      start = as.numeric(expression$start),
      lower = expression$lower,
      upper = expression$upper,
      fixed = isTRUE(expression$fixed),
      sigma_prior = as.numeric(expression$sigma_prior %||% 5),
      prior = if (is.null(expression$prior)) NULL else {
        if (!inherits(expression$prior, "biogeme_prior")) {
          stop("beta prior must be a biogeme_prior object.", call. = FALSE)
        }
        list(
          distribution = expression$prior$distribution,
          sigma = as.numeric(expression$prior$sigma),
          nu = as.numeric(expression$prior$nu)
        )
      }
    ))
  }
  if (identical(kind, "binary")) {
    return(list(
      kind = "binary",
      operator = expression$operator,
      left = biogeme_expression_ir(expression$left),
      right = biogeme_expression_ir(expression$right)
    ))
  }
  if (identical(kind, "unary")) {
    return(list(
      kind = "unary",
      operator = expression$operator,
      operand = biogeme_expression_ir(expression$operand)
    ))
  }
  if (identical(kind, "function")) {
    return(list(
      kind = "function",
      operator = expression$operator,
      args = lapply(expression$args, biogeme_expression_ir),
      attributes = expression$attributes
    ))
  }
  if (identical(kind, "elem")) {
    mapping <- lapply(seq_along(expression$mapping), function(index) {
      list(
        key = as.integer(expression$keys[[index]]),
        expression = biogeme_expression_ir(expression$mapping[[index]])
      )
    })
    return(list(
      kind = "elem",
      mapping = mapping,
      index = biogeme_expression_ir(expression$index)
    ))
  }
  if (identical(kind, "linear_utility")) {
    return(list(
      kind = kind,
      terms = lapply(expression$terms, function(term) {
        list(
          beta = biogeme_expression_ir(term$beta),
          x = biogeme_expression_ir(term$x)
        )
      })
    ))
  }
  if (identical(kind, "segmented_beta")) {
    return(list(
      kind = kind,
      beta = biogeme_expression_ir(expression$beta),
      segmentations = lapply(expression$segmentations, biogeme_segmentation_ir),
      prefix = expression$prefix
    ))
  }
  if (identical(kind, "catalog")) {
    return(list(
      kind = kind,
      name = expression$name,
      controller = list(
        name = expression$controller$name,
        specification_names = expression$controller$specification_names,
        selected_name = expression$controller$selected_name
      ),
      expressions = lapply(seq_along(expression$expressions), function(index) {
        list(
          name = names(expression$expressions)[[index]],
          expression = biogeme_expression_ir(expression$expressions[[index]])
        )
      })
    ))
  }
  if (kind %in% c("logit_probability", "logit_log_probability")) {
    mapping_ir <- function(expressions) {
      if (is.null(expressions)) return(NULL)
      lapply(seq_along(expressions), function(index) {
        list(
          key = as.integer(names(expressions)[[index]]),
          expression = biogeme_expression_ir(expressions[[index]])
        )
      })
    }
    return(list(
      kind = kind,
      utilities = mapping_ir(expression$utilities),
      availability = mapping_ir(expression$availability),
      alternative = if (is_biogeme_expression(expression$alternative)) {
        biogeme_expression_ir(expression$alternative)
      } else {
        as.integer(expression$alternative)
      }
    ))
  }
  if (kind %in% c("ordered_logit", "ordered_probit")) {
    return(list(
      kind = kind,
      eta = biogeme_expression_ir(expression$eta),
      cutpoints = lapply(expression$cutpoints, biogeme_expression_ir),
      alternative = biogeme_expression_ir(expression$alternative),
      categories = expression$categories,
      neutral_labels = expression$neutral_labels,
      enforce_order = isTRUE(expression$enforce_order),
      eps = as.numeric(expression$eps)
    ))
  }
  if (kind %in% c("nested_probability", "nested_logit", "nested_logit_mev_mu")) {
    mapping_ir <- function(expressions) {
      if (is.null(expressions)) return(NULL)
      lapply(seq_along(expressions), function(index) {
        list(
          key = as.integer(names(expressions)[[index]]),
          expression = biogeme_expression_ir(expressions[[index]])
        )
      })
    }
    return(list(
      kind = kind,
      utilities = mapping_ir(expression$utilities),
      availability = mapping_ir(expression$availability),
      nests = biogeme_nested_nests_ir(expression$nests),
      scale_parameter = if (kind == "nested_logit_mev_mu") {
        biogeme_expression_ir(expression$scale_parameter)
      } else {
        NULL
      },
      alternative = if (is_biogeme_expression(expression$alternative)) {
        biogeme_expression_ir(expression$alternative)
      } else {
        as.integer(expression$alternative)
      }
    ))
  }
  if (identical(kind, "nested_endogenous_sampling")) {
    mapping_ir <- function(expressions) {
      if (is.null(expressions)) return(NULL)
      lapply(seq_along(expressions), function(index) {
        list(
          key = as.integer(names(expressions)[[index]]),
          expression = biogeme_expression_ir(expressions[[index]])
        )
      })
    }
    return(list(
      kind = kind,
      utilities = mapping_ir(expression$utilities),
      availability = mapping_ir(expression$availability),
      nests = biogeme_nested_nests_ir(expression$nests),
      correction = mapping_ir(expression$correction),
      alternative = if (is_biogeme_expression(expression$alternative)) {
        biogeme_expression_ir(expression$alternative)
      } else {
        as.integer(expression$alternative)
      }
    ))
  }
  if (kind %in% c(
      "cross_nested_logit",
      "cross_nested_logit_mu",
      "cross_nested_probability"
  )) {
    mapping_ir <- function(expressions) {
      if (is.null(expressions)) return(NULL)
      lapply(seq_along(expressions), function(index) {
        list(
          key = as.integer(names(expressions)[[index]]),
          expression = biogeme_expression_ir(expressions[[index]])
        )
      })
    }
    return(list(
      kind = kind,
      utilities = mapping_ir(expression$utilities),
      availability = mapping_ir(expression$availability),
      nests = biogeme_cross_nested_nests_ir(expression$nests),
      scale_parameter = if (kind == "cross_nested_logit_mu") {
        biogeme_expression_ir(expression$scale_parameter)
      } else {
        NULL
      },
      alternative = if (is_biogeme_expression(expression$alternative)) {
        biogeme_expression_ir(expression$alternative)
      } else {
        as.integer(expression$alternative)
      }
    ))
  }
  if (kind %in% c("draw", "random_variable")) {
    result <- list(kind = kind, name = expression$name)
    if (identical(kind, "draw")) result$draw_type <- expression$draw_type
    return(result)
  }
  if (identical(kind, "distributed_parameter")) {
    return(list(
      kind = kind,
      name = expression$name,
      operand = biogeme_expression_ir(expression$operand)
    ))
  }
  if (identical(kind, "integrate_normal")) {
    return(list(
      kind = kind,
      name = expression$name,
      number_of_quadrature_points = as.integer(expression$number_of_quadrature_points),
      operand = biogeme_expression_ir(expression$operand)
    ))
  }
  if (kind %in% c("monte_carlo", "panel_likelihood_trajectory")) {
    return(list(kind = kind, operand = biogeme_expression_ir(expression$operand)))
  }
  if (identical(kind, "derive")) {
    return(list(
      kind = kind,
      name = expression$name,
      operand = biogeme_expression_ir(expression$operand)
    ))
  }
  if (identical(kind, "boxcox")) {
    return(list(
      kind = kind,
      left = biogeme_expression_ir(expression$left),
      right = biogeme_expression_ir(expression$right)
    ))
  }
  if (identical(kind, "piecewise")) {
    return(list(
      kind = kind,
      operand = biogeme_expression_ir(expression$operand),
      thresholds = expression$thresholds,
      betas = if (is.null(expression$betas)) NULL else
        lapply(expression$betas, biogeme_expression_ir),
      transform = expression$transform
    ))
  }
  stop("Unsupported Biogeme expression kind: ", kind, call. = FALSE)
}

biogeme_expression_mapping <- function(expressions) {
  if (is.null(expressions)) {
    return(NULL)
  }
  lapply(expressions, biogeme_expression_ir)
}

biogeme_mdcev_mapping_ir <- function(expressions, argument, allow_null = FALSE) {
  if (is.null(expressions)) return(NULL)
  if (!is.list(expressions) || is.null(names(expressions))) {
    stop(argument, " must be a named list.", call. = FALSE)
  }
  lapply(seq_along(expressions), function(index) {
    value <- expressions[[index]]
    list(
      key = as.integer(names(expressions)[[index]]),
      expression = if (is.null(value)) {
        if (!isTRUE(allow_null)) {
          stop(argument, " cannot contain NULL values.", call. = FALSE)
        }
        NULL
      } else {
        biogeme_expression_ir(value)
      }
    )
  })
}

biogeme_segmentation_ir <- function(segmentation) {
  if (!inherits(segmentation, "biogeme_segmentation")) {
    stop("segmentation must be a biogeme_segmentation object.", call. = FALSE)
  }
  list(
    variable = segmentation$variable,
    mapping = lapply(seq_along(segmentation$mapping), function(index) {
      list(
        key = as.integer(names(segmentation$mapping)[[index]]),
        label = unname(segmentation$mapping[[index]])
      )
    }),
    reference = segmentation$reference
  )
}

biogeme_nested_nests_ir <- function(nests) {
  if (!inherits(nests, "biogeme_nested_nests")) {
    stop("nests must be a biogeme_nested_nests object.", call. = FALSE)
  }
  list(
    choice_set = as.integer(nests$choice_set),
    nests = lapply(nests$nests, function(nest) {
      list(
        nest_parameter = biogeme_expression_ir(nest$nest_parameter),
        alternatives = as.list(as.integer(nest$alternatives)),
        name = nest$name
      )
    })
  )
}

biogeme_cross_nested_nests_ir <- function(nests) {
  if (!inherits(nests, "biogeme_cross_nested_nests")) {
    stop("nests must be a biogeme_cross_nested_nests object.", call. = FALSE)
  }
  list(
    choice_set = as.integer(nests$choice_set),
    sparse = isTRUE(nests$sparse),
    nests = lapply(nests$nests, function(nest) {
      list(
        nest_parameter = biogeme_expression_ir(nest$nest_parameter),
        allocation = lapply(seq_along(nest$allocation), function(index) {
          list(
            key = as.integer(names(nest$allocation)[[index]]),
            expression = biogeme_expression_ir(nest$allocation[[index]])
          )
        }),
        name = nest$name
      )
    })
  )
}

biogeme_sampling_partition_ir <- function(partition, argument) {
  if (!inherits(partition, "biogeme_sampling_partition")) {
    stop(argument, " must be a biogeme_sampling_partition object.", call. = FALSE)
  }
  list(
    segments = lapply(partition$segments, as.list),
    sample_sizes = as.list(as.integer(partition$sample_sizes)),
    full_set = as.list(as.integer(partition$full_set))
  )
}

biogeme_cross_variables_ir <- function(combined_variables) {
  lapply(combined_variables, function(item) {
    if (!inherits(item, "biogeme_cross_variable")) {
      stop("combined_variables must contain cross_variable() objects.", call. = FALSE)
    }
    list(name = item$name, formula = biogeme_expression_ir(item$formula))
  })
}

biogeme_draw_metadata_ir <- function(draws) {
  if (is.null(draws)) return(NULL)
  if (inherits(draws, "biogeme_draws")) draws <- list(draws)
  if (!is.list(draws)) stop("draws must be a biogeme_draws object or list of them.", call. = FALSE)
  lapply(draws, function(item) {
    if (!inherits(item, "biogeme_draws")) stop("draws must contain biogeme_draws objects.", call. = FALSE)
    list(
      name = item$name,
      draw_type = item$draw_type,
      number_of_draws = item$number_of_draws,
      seed = item$seed,
      matrix = item$matrix,
      generator = item$generator
    )
  })
}

biogeme_parameter_overrides_ir <- function(overrides) {
  if (is.null(overrides)) return(NULL)
  if (!is.list(overrides) || is.null(names(overrides)) ||
      anyNA(names(overrides)) || any(!nzchar(names(overrides))) ||
      anyDuplicated(names(overrides))) {
    stop("parameter_overrides must be a named list of expressions.", call. = FALSE)
  }
  setNames(lapply(overrides, biogeme_expression_ir), names(overrides))
}

biogeme_database_operations_ir <- function(database) {
  validate_biogeme_database(database)
  operations <- list()
  if (length(database$derived_variables) > 0L) {
    for (name in names(database$derived_variables)) {
      operations[[length(operations) + 1L]] <- list(
        kind = "define_variable",
        name = name,
        expression = biogeme_expression_ir(database$derived_variables[[name]])
      )
    }
  }
  if (length(database$filters) > 0L) {
    for (condition in database$filters) {
      operations[[length(operations) + 1L]] <- list(
        kind = "remove",
        expression = biogeme_expression_ir(condition)
      )
    }
  }
  operations
}

biogeme_materialize_database <- function(database) {
  validate_biogeme_database(database)
  if (isTRUE(database$materialized) && length(database$derived_variables) == 0L &&
      length(database$filters) == 0L) {
    return(database)
  }
  result <- tryCatch(
    biogeme_bridge()$materialize_database(
      name = database$name,
      data = reticulate::r_to_py(database$data),
      operations = reticulate::r_to_py(biogeme_database_operations_ir(database)),
      panel_id = database$panel_id,
      row_ids = reticulate::r_to_py(as.character(database$row_ids))
    ),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "native database materialization",
        suggestion = "check derived expressions, filter conditions, and the panel identifier"
      )
    }
  )
  converted <- reticulate::py_to_r(result)
  materialized <- database
  materialized$data <- validate_biogeme_data(as.data.frame(converted$data, check.names = FALSE))
  materialized$row_ids <- as.character(converted$row_ids)
  materialized$derived_variables <- list()
  materialized$filters <- list()
  materialized$filtered_row_count <- as.integer(converted$filtered_row_count)
  materialized$number_of_excluded_data <- as.integer(converted$number_of_excluded_data)
  materialized$materialized <- TRUE
  materialized
}

#' Materialize lazy database operations through native Biogeme
#' @param database A `biogeme_database` object.
#' @return A database with native derived columns and filters applied.
#' @export
biogeme_database_materialize <- biogeme_materialize_database

# Compile the R model into Python objects.  The returned value is a Python
# dictionary and is intentionally not converted back to R: it contains native
# Biogeme expressions that must stay on the Python side.
biogeme_compile_model <- function(model, database_override = NULL) {
  if (!inherits(model, "biogeme_model")) {
    biogeme_abort(
      "model must be a Biogeme model.",
      class = "biogeme_specification_error",
      operation = "model compilation",
      suggestion = "construct the model with logit_model() or biogeme_model()"
    )
  }
  tryCatch(
    {
      bridge <- biogeme_bridge()
      database_spec <- if (is.null(database_override)) model$database else {
        validate_biogeme_database(database_override)
        database_override
      }
      if (!is.null(model$subset)) {
        database_spec <- biogeme_database_remove(database_spec, !model$subset)
      }
      data <- reticulate::r_to_py(database_spec$data)
      database <- bridge$build_database(
        name = database_spec$name,
        data = data,
        operations = reticulate::r_to_py(biogeme_database_operations_ir(database_spec)),
        panel_id = database_spec$panel_id,
        row_ids = reticulate::r_to_py(as.character(database_spec$row_ids))
      )
      if (inherits(model, "biogeme_mdcev_model")) {
        return(bridge$compile_mdcev_model(
          database = database,
          model_name = "rbiogeme_mdcev",
          model_type = model$model_type,
          baseline_utilities = reticulate::r_to_py(
            biogeme_mdcev_mapping_ir(model$baseline_utilities, "baseline_utilities")
          ),
          gamma_parameters = reticulate::r_to_py(
            biogeme_mdcev_mapping_ir(
              model$gamma_parameters,
              "gamma_parameters",
              allow_null = TRUE
            )
          ),
          alpha_parameters = if (is.null(model$alpha_parameters)) NULL else
            reticulate::r_to_py(
              biogeme_mdcev_mapping_ir(model$alpha_parameters, "alpha_parameters")
            ),
          mu_utilities = if (is.null(model$mu_utilities)) NULL else
            reticulate::r_to_py(
              biogeme_mdcev_mapping_ir(model$mu_utilities, "mu_utilities")
            ),
          scale_parameter = if (is.null(model$scale_parameter)) NULL else
            reticulate::r_to_py(biogeme_expression_ir(model$scale_parameter)),
          prices = if (is.null(model$prices)) NULL else
            reticulate::r_to_py(
              biogeme_mdcev_mapping_ir(model$prices, "prices")
            ),
          weights = if (is.null(model$weights)) NULL else
            reticulate::r_to_py(biogeme_expression_ir(model$weights)),
          number_of_chosen_alternatives = reticulate::r_to_py(
            biogeme_expression_ir(model$number_of_chosen_alternatives)
          ),
          consumed_quantities = reticulate::r_to_py(
            biogeme_mdcev_mapping_ir(
              model$consumed_quantities,
              "consumed_quantities"
            )
          )
        ))
      }
      if (inherits(model, "biogeme_logit_model")) {
      utilities <- reticulate::r_to_py(biogeme_expression_mapping(model$utilities))
      availability <- if (is.null(model$availability)) NULL else
        reticulate::r_to_py(biogeme_expression_mapping(model$availability))
      alternative_codes <- reticulate::r_to_py(as.list(model$alternative_codes))
        weight <- if (is.null(model$weight)) NULL else
          reticulate::r_to_py(biogeme_expression_ir(model$weight))
        return(bridge$compile_logit_model(
          database = database,
          choice = model$choice,
          utilities = utilities,
          availability = availability,
          alternative_codes = alternative_codes,
          weight = weight,
          formula = NULL,
          simulations = if (is.null(model$simulations)) NULL else
            reticulate::r_to_py(biogeme_expression_mapping(model$simulations)),
          draw_metadata = reticulate::r_to_py(biogeme_draw_metadata_ir(model$draws))
        ))
      }
      if (inherits(model, "biogeme_nested_logit_model")) {
        utilities <- reticulate::r_to_py(biogeme_expression_mapping(model$utilities))
        availability <- if (is.null(model$availability)) NULL else
          reticulate::r_to_py(biogeme_expression_mapping(model$availability))
        alternative_codes <- reticulate::r_to_py(as.list(model$alternative_codes))
        return(bridge$compile_nested_logit_model(
          database = database,
          choice = model$choice,
          utilities = utilities,
          availability = availability,
          nests = reticulate::r_to_py(biogeme_nested_nests_ir(model$nests)),
          alternative_codes = alternative_codes,
          scale_parameter = if (is.null(model$scale_parameter)) NULL else
            reticulate::r_to_py(biogeme_expression_ir(model$scale_parameter)),
          weight = NULL,
          simulations = if (is.null(model$simulations)) NULL else
            reticulate::r_to_py(biogeme_expression_mapping(model$simulations)),
          draw_metadata = reticulate::r_to_py(biogeme_draw_metadata_ir(model$draws))
        ))
      }
      if (inherits(model, "biogeme_cross_nested_logit_model")) {
        utilities <- reticulate::r_to_py(biogeme_expression_mapping(model$utilities))
        availability <- if (is.null(model$availability)) NULL else
          reticulate::r_to_py(biogeme_expression_mapping(model$availability))
        alternative_codes <- reticulate::r_to_py(as.list(model$alternative_codes))
        return(bridge$compile_cross_nested_logit_model(
          database = database,
          choice = model$choice,
          utilities = utilities,
          availability = availability,
          nests = reticulate::r_to_py(biogeme_cross_nested_nests_ir(model$nests)),
          alternative_codes = alternative_codes,
          scale_parameter = if (is.null(model$scale_parameter)) NULL else
            reticulate::r_to_py(biogeme_expression_ir(model$scale_parameter)),
          weight = NULL,
          simulations = if (is.null(model$simulations)) NULL else
            reticulate::r_to_py(biogeme_expression_mapping(model$simulations)),
          draw_metadata = reticulate::r_to_py(biogeme_draw_metadata_ir(model$draws))
        ))
      }
      formula <- model$formula
      if (is.null(formula)) formula <- model$log_likelihood
      simulations <- model$simulations
      bridge$compile_generic_model(
        database = database,
        formula = if (is.null(formula)) NULL else reticulate::r_to_py(
          if (is.list(formula) && !is_biogeme_expression(formula))
            lapply(formula, biogeme_expression_ir)
          else biogeme_expression_ir(formula)
        ),
        probability = if (is.null(model$probability)) NULL else
          reticulate::r_to_py(biogeme_expression_ir(model$probability)),
        simulations = if (is.null(simulations)) NULL else
          reticulate::r_to_py(biogeme_expression_mapping(simulations)),
        weight = if (is.null(model$weight)) NULL else
          reticulate::r_to_py(biogeme_expression_ir(model$weight)),
        availability = if (is.null(model$availability)) NULL else
          reticulate::r_to_py(
            biogeme_mdcev_mapping_ir(model$availability, "availability")
          ),
        panel_trajectory = isTRUE(model$panel_trajectory),
        draw_metadata = reticulate::r_to_py(biogeme_draw_metadata_ir(model$draws)),
        parameter_overrides = reticulate::r_to_py(
          biogeme_parameter_overrides_ir(model$parameter_overrides)
        )
      )
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "model compilation",
        suggestion = "check the database columns, expressions, alternatives, and availability definitions"
      )
    }
  )
}

# Run native model construction and validation without estimation. The bridge
# returns only serialisable diagnostics; the native BIOGEME object stays in
# Python and is discarded after the check.
biogeme_validate_model <- function(
    model,
    model_name = "rbiogeme_validation",
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$validate_model_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls
      ))
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "Biogeme model validation",
        suggestion = "check the database columns, expressions, alternatives, availability, and panel declaration"
      )
    }
  )
}

# Run native Biogeme's analytical-versus-finite-difference derivative check.
# The returned object is converted to ordinary R vectors only after the
# native check has completed.
biogeme_check_derivatives_model <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    verbose = FALSE
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      result <- reticulate::py_to_r(bridge$check_derivatives_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        verbose = isTRUE(verbose)
      ))
      vector_fields <- c(
        "analytical_gradient",
        "finite_differences_gradient",
        "errors_gradient"
      )
      matrix_fields <- c(
        "analytical_hessian",
        "finite_differences_hessian",
        "errors_hessian"
      )
      for (field in vector_fields) {
        result[[field]] <- as.numeric(unlist(result[[field]], use.names = FALSE))
      }
      for (field in matrix_fields) {
        value <- result[[field]]
        if (!is.matrix(value)) {
          rows <- lapply(value, function(row) {
            as.numeric(unlist(row, use.names = FALSE))
          })
          value <- do.call(rbind, rows)
        }
        result[[field]] <- unname(value)
      }
      result[["function"]] <- as.numeric(result[["function"]])
      result
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme derivative check",
        suggestion = "check the model expression and its starting values"
      )
    }
  )
}

# Run native JAX execution profiling for a compiled model.  The evaluator and
# profiler remain Python-side objects; only ordinary environment and timing
# summaries are converted back to R.
biogeme_profile_jax_model <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    beta_values,
    cases,
    numerically_safe = FALSE
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$profile_jax_biogeme(
        compiled_model = compiled,
        beta_values = reticulate::r_to_py(as.list(beta_values)),
        cases = reticulate::r_to_py(cases),
        model_name = model_name,
        controls = controls,
        numerically_safe = isTRUE(numerically_safe)
      ))
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "native JAX profiling",
        suggestion = "check the model, beta values, profiling cases, and JAX support"
      )
    }
  )
}

biogeme_check_monte_carlo_stability_model <- function(
    model,
    fit,
    model_name = "rbiogeme_model",
    controls = list(),
    output_directory = NULL,
    basename = NULL,
    resume = TRUE
) {
  compiled <- biogeme_compile_model(model)
  fit_for_bridge <- fit
  fit_for_bridge$model <- NULL
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$check_monte_carlo_stability_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        estimation_results = reticulate::r_to_py(fit_for_bridge),
        output_directory = output_directory,
        basename = basename,
        resume = isTRUE(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"
      )
    }
  )
}

# Generate sampled alternatives and estimate the resulting native model.  The
# sampling context, generated data, and sampled likelihood all remain on the
# Python side; R receives only serialized results and metadata.
biogeme_estimate_sampled_alternatives_model <- function(
    model,
    model_name = "rbiogeme_sampled",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE
) {
  if (!inherits(model, "biogeme_sampled_alternatives_model")) {
    stop("model must be a sampled_alternatives_model object.", call. = FALSE)
  }
  tryCatch(
    {
      bridge <- biogeme_bridge()
      raw <- bridge$estimate_sampled_alternatives_biogeme(
        alternatives = reticulate::r_to_py(model$alternatives),
        individuals = reticulate::r_to_py(model$individuals),
        choice_column = model$choice_column,
        id_column = model$id_column,
        utility = reticulate::r_to_py(biogeme_expression_ir(model$utility)),
        combined_variables = reticulate::r_to_py(
          biogeme_cross_variables_ir(model$combined_variables)
        ),
        partition = reticulate::r_to_py(
          biogeme_sampling_partition_ir(model$partition, "partition")
        ),
        sample_file_name = model$biogeme_file_name,
        model_type = model$model_type,
        mev_partition = if (is.null(model$mev_partition)) NULL else
          reticulate::r_to_py(
            biogeme_sampling_partition_ir(model$mev_partition, "mev_partition")
          ),
        mev_sample_sizes = if (is.null(model$mev_sample_sizes)) NULL else
          reticulate::r_to_py(as.list(as.integer(model$mev_sample_sizes))),
        nests = if (is.null(model$nests)) NULL else if (
          inherits(model$nests, "biogeme_nested_nests")
        ) {
          reticulate::r_to_py(biogeme_nested_nests_ir(model$nests))
        } else {
          reticulate::r_to_py(biogeme_cross_nested_nests_ir(model$nests))
        },
        model_name = model_name,
        controls = if (is.null(controls) || 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)
      )
      reticulate::py_to_r(raw)
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "sampled-alternative Biogeme estimation",
        suggestion = "check the alternative IDs, sampling partitions, and model specification"
      )
    }
  )
}

# Internal helper used by the public estimation function.
biogeme_estimate_model <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE,
    yaml_file_name = NULL
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$estimate_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        starting_values = if (is.null(starting_values)) NULL else reticulate::r_to_py(starting_values),
        run_bootstrap = 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"
      )
    }
  )
}

biogeme_estimate_configuration_model <- function(
    model,
    configuration_id,
    model_name = "rbiogeme_configuration",
    controls = list(),
    starting_values = NULL,
    run_bootstrap = FALSE,
    yaml_file_name = NULL
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$estimate_configuration_biogeme(
        compiled_model = compiled,
        configuration_id = configuration_id,
        model_name = model_name,
        controls = 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
      )
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme selected configuration estimation",
        suggestion = "check the configuration ID, catalog controllers, and model expression"
      )
    }
  )
}

biogeme_native_parameter_names_model <- function(
    model,
    configuration_id = NULL,
    model_name = "rbiogeme_parameters",
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$native_parameter_names_biogeme(
        compiled_model = compiled,
        configuration_id = configuration_id,
        model_name = model_name,
        controls = controls
      ))
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_specification_error",
        operation = "native parameter-name discovery",
        suggestion = "check the configuration ID and model expression"
      )
    }
  )
}

biogeme_catalog_configuration_ids_model <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$catalog_configuration_ids_biogeme(
        compiled_model = compiled,
        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"
      )
    }
  )
}

biogeme_bayesian_estimate_model <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list(),
    starting_values = NULL
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$bayesian_estimate_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        starting_values = if (is.null(starting_values)) NULL else
          reticulate::r_to_py(starting_values)
      ))
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Bayesian Biogeme estimation",
        suggestion = "check the model expression, Bayesian controls, starting values, and Python Bayesian dependencies"
      )
    }
  )
}

biogeme_simulate_bayesian_model <- function(
    model,
    bayesian_results_file,
    model_name = "rbiogeme_bayesian_simulation",
    percentage_of_draws_to_use = 10,
    lower_quantile = 0.025,
    upper_quantile = 0.975,
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      reticulate::py_to_r(bridge$simulate_bayesian_biogeme(
        compiled_model = compiled,
        bayesian_results_file = bayesian_results_file,
        model_name = model_name,
        percentage_of_draws_to_use = as.numeric(percentage_of_draws_to_use),
        lower_quantile = as.numeric(lower_quantile),
        upper_quantile = as.numeric(upper_quantile),
        controls = controls
      ))
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Bayesian Biogeme simulation",
        suggestion = "check the NetCDF result path, simulation expressions, and posterior-draw controls"
      )
    }
  )
}

biogeme_posterior_mean_by_observation <- function(
    bayesian_results_file,
    variable_name,
    controls = list()
) {
  if (!is.character(bayesian_results_file) || length(bayesian_results_file) != 1L ||
      is.na(bayesian_results_file) || !nzchar(bayesian_results_file)) {
    stop("bayesian_results_file must be one non-empty path.", call. = FALSE)
  }
  if (!is.character(variable_name) || length(variable_name) != 1L ||
      is.na(variable_name) || !nzchar(variable_name)) {
    stop("variable_name must be one non-empty string.", call. = FALSE)
  }
  controls <- if (is.null(controls) || length(controls) == 0L) {
    NULL
  } else {
    reticulate::r_to_py(controls)
  }
  tryCatch(
    reticulate::py_to_r(biogeme_bridge()$posterior_mean_by_observation_biogeme(
      bayesian_results_file = bayesian_results_file,
      variable_name = variable_name,
      controls = controls
    )),
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Bayesian posterior observation summary",
        suggestion = "check the NetCDF result path and the stored variable name"
      )
    }
  )
}

biogeme_estimate_catalog_model <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list(),
    quick_estimate = FALSE,
    run_bootstrap = FALSE,
    force = TRUE
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$estimate_catalog_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        quick_estimate = isTRUE(quick_estimate),
        run_bootstrap = isTRUE(run_bootstrap),
        force = isTRUE(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"
      )
    }
  )
}

biogeme_assisted_specification_model <- function(
    model,
    model_name = "rbiogeme_assisted",
    controls = list(),
    objectives = "loglikelihood_dimension",
    pareto_file_name = NULL,
    force = TRUE,
    validity = NULL
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$assisted_specification_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        objectives = objectives,
        pareto_file_name = pareto_file_name,
        force = isTRUE(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"
      )
    }
  )
}

biogeme_count_number_of_specifications_model <- function(
    model,
    model_name = "rbiogeme_catalog",
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$count_number_of_specifications_biogeme(
        compiled_model = compiled,
        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"
      )
    }
  )
}

biogeme_pareto_post_processing_model <- function(
    model,
    model_name = "rbiogeme_pareto",
    controls = list(),
    pareto_file_name = NULL,
    recycle = FALSE,
    plot_file_name = NULL,
    objective_x = 0L,
    objective_y = 1L,
    label_x = NULL,
    label_y = NULL
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$pareto_post_processing_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls,
        pareto_file_name = pareto_file_name,
        recycle = isTRUE(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 model data"
      )
    }
  )
}

biogeme_quick_estimate_model <- function(
    model,
    model_name = "rbiogeme_model",
    controls = list()
) {
  compiled <- biogeme_compile_model(model)
  tryCatch(
    {
      bridge <- biogeme_bridge()
      controls <- if (is.null(controls) || length(controls) == 0L) {
        NULL
      } else {
        reticulate::r_to_py(controls)
      }
      bridge$quick_estimate_biogeme(
        compiled_model = compiled,
        model_name = model_name,
        controls = controls
      )
    },
    error = function(error) {
      biogeme_rethrow(
        error,
        class = "biogeme_estimation_error",
        operation = "Biogeme quick estimation",
        suggestion = "check the model specification, starting values, and data scaling"
      )
    }
  )
}