inst/examples/swissmetro/plot_b21c_process_pareto.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# b21c. Re-estimate the Pareto-optimal models
#
# This example mirrors plot_b21c_process_pareto.py. The complete b21b model
# specification is included below so this R script is self-contained. Native
# Biogeme reads the Pareto file, re-estimates the selected configurations,
# compiles the summary, and creates the optional Pareto plot.

library(rbiogeme)

# 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 the output directory. It does not define this model.
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_b21c_multiple_models_model <- function(database) {
  # These parameters and starting values match native b21b.
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)

  # The explicit noGA reference preserves native segmented-parameter names.
  gender_segmentation <- biogeme_database_segmentation(
    database,
    "MALE",
    c(`0` = "female", `1` = "male")
  )
  ga_segmentation <- biogeme_database_segmentation(
    database,
    "GA",
    c(`1` = "GA", `0` = "noGA"),
    reference = "noGA"
  )
  income_segmentation <- biogeme_database_segmentation(
    database,
    "INCOME",
    c(
      `0` = "inc-zero",
      `1` = "inc-under50",
      `2` = "inc-50-100",
      `3` = "inc-100+",
      `4` = "inc-unknown"
    )
  )

  # segmentation_catalogs() creates the same native catalog choices as b21b.
  asc_catalogs <- segmentation_catalogs(
    "asc",
    list(asc_car, asc_train),
    list(gender_segmentation, ga_segmentation),
    maximum_number = 2
  )
  b_cost_catalog <- segmentation_catalogs(
    "b_cost",
    list(b_cost),
    list(ga_segmentation, income_segmentation),
    maximum_number = 1
  )[[1L]]

  # All three travel-time catalogs share one native controller.
  lambda_time <- biogeme_beta("lambda_time", start = 1, lower = -10, upper = 10)
  time_controller <- catalog_controller("train_tt", c("linear", "log", "boxcox"))
  train_tt_catalog <- catalog(
    "train_tt",
    list(
      linear = variable("TRAIN_TT_SCALED"),
      log = logzero(variable("TRAIN_TT_SCALED")),
      boxcox = boxcox(variable("TRAIN_TT_SCALED"), lambda_time)
    ),
    time_controller
  )
  sm_tt_catalog <- catalog(
    "sm_tt",
    list(
      linear = variable("SM_TT_SCALED"),
      log = logzero(variable("SM_TT_SCALED")),
      boxcox = boxcox(variable("SM_TT_SCALED"), lambda_time)
    ),
    time_controller
  )
  car_tt_catalog <- catalog(
    "car_tt",
    list(
      linear = variable("CAR_TT_SCALED"),
      log = logzero(variable("CAR_TT_SCALED")),
      boxcox = boxcox(variable("CAR_TT_SCALED"), lambda_time)
    ),
    time_controller
  )

  utilities <- list(
    `1` = asc_catalogs[[2L]] + b_time * train_tt_catalog +
      b_cost_catalog * variable("TRAIN_COST_SCALED"),
    `2` = b_time * sm_tt_catalog + b_cost_catalog * variable("SM_COST_SCALED"),
    `3` = asc_catalogs[[1L]] + b_time * car_tt_catalog +
      b_cost_catalog * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )

  biogeme_model(
    database = database,
    formula = logit_log_probability(
      utilities = utilities,
      availability = availability,
      alternative = variable("CHOICE")
    ),
    control = biogeme_control(
    output_directory = prepared$output,
      model_name = "b21_multiple_models",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "b21c_process_pareto"
)
database <- swissmetro_data(prepared$data)
model <- build_b21c_multiple_models_model(database)

make_plot <- example_flag(prepared$options$plot, default = TRUE)
reuse_pareto <- example_flag(prepared$options$reuse_pareto, default = FALSE)
pareto_file <- file.path(prepared$output, "b21_multiple_models.pareto")
csv_file <- file.path(prepared$output, "b21_process_pareto.csv")
plot_file <- file.path(prepared$output, "b21_process_pareto.png")

# The native b21c example expects b21a to have produced this Pareto file. To
# make this counterpart runnable from a clean directory, create that exact
# prerequisite afresh unless --reuse-pareto=true is explicitly requested.
if (!reuse_pareto && file.exists(pareto_file)) unlink(pareto_file, force = TRUE)
if (!file.exists(pareto_file)) {
  cat("Creating the native b21a Pareto prerequisite...\n")
  invisible(assisted_specification(
    model,
    objectives = "loglikelihood_dimension",
    pareto_file_name = pareto_file,
    model_name = "b21_multiple_models",
    control = model$control,
    force = TRUE
  ))
}
if (file.exists(csv_file)) unlink(csv_file, force = TRUE)
if (file.exists(plot_file)) unlink(plot_file, force = TRUE)

# recycle=FALSE matches native b21c and guarantees fresh complete estimates.
fit <- pareto_post_processing(
  model,
  pareto_file_name = pareto_file,
  model_name = "b21_multiple_models",
  control = model$control,
  recycle = FALSE,
  plot_file_name = if (make_plot) plot_file else NULL,
  objective_x = 0L,
  objective_y = 1L,
  label_x = "Negative loglikelihood",
  label_y = "Number of parameters"
)

cat(paste(fit$pareto_statistics, collapse = "\n"), "\n", sep = "")
print(fit$summary)
write.csv(fit$summary, csv_file, row.names = TRUE, quote = TRUE)
cat(sprintf("Summary table available in %s\n", basename(csv_file)))

# Append the same short-name explanations that native b21c writes to CSV.
cat("\n\n", file = csv_file, append = TRUE)
for (name in names(fit$description)) {
  if (!identical(name, unname(fit$description[[name]]))) {
    cat(sprintf("%s: %s\n", name, fit$description[[name]]))
    cat(sprintf("%s,%s\n", name, fit$description[[name]]), file = csv_file, append = TRUE)
  }
}

invisible(fit)