inst/examples/swissmetro/plot_b01e_logit_all_algos.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b01e. Logit estimation with several native algorithms/settings
#
# This example estimates the Swissmetro logit model for every combination of
# three native Biogeme controls. It records the result of each run in a CSV
# summary, including failures, as in the corresponding Python example.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# complete model specification remains in this script.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "example_utils.R"))

build_b01e_model <- function(database) {
  # biogeme_beta() creates symbolic parameters. The Swissmetro ASC is fixed
  # at zero so the model is identified.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)

  # variable() and arithmetic construct a symbolic utility tree. The derived
  # scaled columns were prepared by swissmetro_data(), matching the native
  # Swissmetro data-processing module.
  logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm + b_time * variable("SM_TT_SCALED") +
        b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * variable("CAR_TT_SCALED") +
        b_cost * variable("CAR_CO_SCALED")
    ),
    availability = list(
      `1` = variable("TRAIN_AV_SP"),
      `2` = variable("SM_AV"),
      `3` = variable("CAR_AV_SP")
    )
  )
}

format_number <- function(value) {
  formatC(value, format = "f", digits = 1)
}

format_optimization_time <- function(seconds) {
  if (is.null(seconds) || length(seconds) == 0L || is.na(seconds)) {
    return(NA_character_)
  }
  sprintf("%.6f seconds", as.numeric(seconds))
}

# prepare_swissmetro_example() is defined in example_utils.R. It parses the
# command line, validates the data/Python paths, configures the bridge, reads
# the data, and creates a fresh output directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b01e_logit_all_algos"
)

# The product order matches itertools.product(
#   [TRUE, FALSE], [0.1, 1.0, 10.0], [0.0, 0.5, 1.0]
# ) in the native Python example.
settings_grid <- expand.grid(
  second_derivatives = c(0.0, 0.5, 1.0),
  initial_radius = c(0.1, 1.0, 10.0),
  infeasible_cg = c(TRUE, FALSE),
  KEEP.OUT.ATTRS = FALSE,
  stringsAsFactors = FALSE
)

# This script always estimates afresh. Remove only the exact summary file it
# writes so an existing output directory cannot make an old summary appear to
# be the result of the current run.
summary_file <- file.path(prepared$output, "01e_logit_all_algos.csv")
if (file.exists(summary_file)) unlink(summary_file, force = TRUE)

database <- swissmetro_data(prepared$data)
model <- build_b01e_model(database)
summary_rows <- vector("list", nrow(settings_grid))
first <- TRUE

for (index in seq_len(nrow(settings_grid))) {
  settings <- settings_grid[index, , drop = FALSE]
  infeasible_cg <- isTRUE(settings$infeasible_cg)
  initial_radius <- as.numeric(settings$initial_radius)
  second_derivatives <- as.numeric(settings$second_derivatives)
  suffix <- paste0(
    "cg_", infeasible_cg,
    "_radius_", format_number(initial_radius),
    "_second_deriv_", format_number(second_derivatives)
  )
  # Preserve the native Python model-name convention, including its original
  # prefix, because model names are part of the equivalence contract.
  native_model_name <- paste0("b05normal_mixture_algo_", suffix)
  result_data <- data.frame(
    InfeasibleCG = infeasible_cg,
    InitialRadius = initial_radius,
    SecondDerivatives = second_derivatives,
    Status = "Success",
    LogLikelihood = NA_real_,
    GradientNorm = NA_real_,
    `Number of draws` = NA_real_,
    `Optimization time` = NA_character_,
    TerminationCause = NA_character_,
    check.names = FALSE,
    stringsAsFactors = FALSE
  )

  message(sprintf("Running %d/%d: %s", index, nrow(settings_grid), suffix))
  fit <- tryCatch(
    {
      controls <- biogeme_control(
    output_directory = prepared$output,
        infeasible_cg = infeasible_cg,
        initial_radius = initial_radius,
        second_derivatives_percentage = second_derivatives,
        generate_html = FALSE,
        generate_yaml = FALSE,
        save_iterations = FALSE
      )
      current_fit <- estimate(
        model,
        model_name = native_model_name,
        control = controls
      )
      # The native example estimates the first configuration twice to warm up
      # Python code before timing the remaining configurations.
      if (first) {
        current_fit <- estimate(
          model,
          model_name = native_model_name,
          control = controls
        )
        first <- FALSE
      }
      current_fit
    },
    error = function(error) error
  )

  if (inherits(fit, "error")) {
    result_data$Status <- "Failed"
    result_data$TerminationCause <- conditionMessage(fit)
  } else {
    result_data$LogLikelihood <- as.numeric(fit$final_log_likelihood)
    result_data$GradientNorm <- if (is.null(fit$gradient_norm)) {
      NA_real_
    } else {
      as.numeric(fit$gradient_norm)
    }
    result_data$`Number of draws` <- if (is.null(fit$number_of_draws)) {
      NA_real_
    } else {
      as.numeric(fit$number_of_draws)
    }
    result_data$`Optimization time` <- format_optimization_time(fit$optimization_time)
    result_data$TerminationCause <- fit$termination_reason
  }
  summary_rows[[index]] <- result_data
}

summary <- do.call(rbind, summary_rows)
print(summary)
write.csv(summary, summary_file, row.names = FALSE, quote = TRUE)
message("Summary reported in file ", summary_file)

invisible(summary)