inst/examples/swissmetro/plot_b05d_normal_mixture_all_algos.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b05d. Normal mixture estimated with several native algorithms/settings
#
# This example runs the same 18 combinations as the native Python example.
# Every run uses the same symbolic Monte Carlo likelihood and changes only
# native Biogeme optimizer controls. Results are summarized in a CSV file.

library(rbiogeme)

# The shared helper contains command-line parsing and data preparation. The
# complete random-coefficient 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_b05d_normal_mixture_model <- function(
    database,
    number_of_draws = 10000L,
    seed = 1223L
) {
  # These names and starting values are identical to native b05a/b05d. The
  # Swissmetro ASC is fixed at zero to identify the utility scale.
  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_cost <- biogeme_beta("b_cost", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1)

  # draw() is a symbolic native Draws node, not an R random number. The
  # complete expression tree is compiled once for each native estimation.
  b_time_rnd <- b_time + b_time_s * draw("b_time_rnd", "NORMAL")
  utilities <- list(
    `1` = asc_train + b_time_rnd * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time_rnd * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time_rnd * 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")
  )
  conditional_probability <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )
  draws <- biogeme_draws(
    name = "b_time_rnd",
    draw_type = "NORMAL",
    number_of_draws = number_of_draws,
    seed = seed
  )

  biogeme_model(
    database = database,
    formula = log(monte_carlo(conditional_probability)),
    draws = draws
  )
}

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. The --data, --python,
# --output, --draws, and --seed options work from any current working
# directory.
prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b05d_normal_mixture_all_algos"
)

number_of_draws <- if (!is.null(prepared$options$draws) && nzchar(prepared$options$draws)) {
  example_integer(prepared$options$draws, "draws")
} else {
  10000L
}
seed <- if (!is.null(prepared$options$seed) && nzchar(prepared$options$seed)) {
  example_integer(prepared$options$seed, "seed")
} else {
  1223L
}

# The native example uses itertools.product(
# [True, False], [0.1, 1.0, 10.0], [0.0, 0.5, 1.0]). This expand.grid order
# produces the same 18 rows and therefore the same model-name sequence.
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
)

# Always estimate afresh. The native script can recycle saved YAML files;
# this R example removes only the exact artifacts it can create and never
# uses estimate_or_load().
summary_file <- file.path(prepared$output, "05d_normal_mixture_all_algos.csv")
stale_files <- c(
  summary_file,
  list.files(
    prepared$output,
    pattern = "^(b05normal_mixture_algo_.*\\.(yaml|html)|__b05normal_mixture_algo_.*\\.iter)$",
    all.files = FALSE,
    full.names = TRUE
  )
)
stale_files <- unique(stale_files[file.exists(stale_files)])
if (length(stale_files) > 0L) unlink(stale_files, force = TRUE)

database <- swissmetro_data(prepared$data)
model <- build_b05d_normal_mixture_model(
  database,
  number_of_draws = number_of_draws,
  seed = seed
)
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)
  )
  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,
        number_of_draws = number_of_draws,
        seed = seed,
        infeasible_cg = infeasible_cg,
        initial_radius = initial_radius,
        second_derivatives_percentage = second_derivatives,
        analytical_hessian_mode = "automatic",
        generate_html = FALSE,
        generate_yaml = FALSE,
        save_iterations = FALSE
      )
      current_fit <- estimate(
        model,
        model_name = native_model_name,
        control = controls
      )
      # Native b05d repeats the first estimation to warm up Python before
      # comparing optimization times. Preserve that workflow exactly.
      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)