tests/testthat/test-swissmetro-b21c.R

In rbiogeme: Interface to 'Biogeme'

native_b21c_specification <- function(data_path, pareto_file_name) {
  native_examples <- "/Users/bierlair/MyFiles/github/biogeme/docs/source/examples/swissmetro"
  sys <- reticulate::import("sys", convert = FALSE)
  sys$path$insert(0L, native_examples)
  old_directory <- getwd()
  setwd(dirname(pareto_file_name))
  on.exit(setwd(old_directory), add = TRUE)
  file.copy(data_path, file.path(getwd(), "swissmetro.dat"), overwrite = TRUE)

  native_environment <- new.env(parent = emptyenv())
  reticulate::source_python(
    file.path(native_examples, "plot_b21b_multiple_models_spec.py"),
    envir = native_environment
  )
  native_biogeme <- native_environment$the_biogeme
  assisted_module <- reticulate::import("biogeme.assisted", convert = FALSE)
  objectives_module <- reticulate::import("biogeme.multiobjectives", convert = FALSE)
  results_processing <- reticulate::import(
    "biogeme.results_processing",
    convert = FALSE
  )
  specification <- reticulate::import(
    "biogeme.catalog.specification",
    convert = FALSE
  )$Specification
  specification$all_results <- reticulate::dict()
  specification$model_names <- NULL

  assisted <- assisted_module$AssistedSpecification(
    biogeme_object = native_biogeme,
    multi_objectives = objectives_module$loglikelihood_dimension,
    pareto_file_name = pareto_file_name
  )
  assisted$run()
  post_processing <- assisted_module$ParetoPostProcessing(
    biogeme_object = native_biogeme,
    pareto_file_name = pareto_file_name
  )
  native_results <- post_processing$reestimate(recycle = FALSE)
  compiled <- results_processing$compile_estimation_results(
    native_results,
    use_short_names = TRUE
  )
  native_summary <- reticulate::py_to_r(reticulate::py_get_item(compiled, 0L))
  bridge <- rbiogeme:::biogeme_bridge()
  keys <- vapply(
    reticulate::iterate(native_results$keys()),
    as.character,
    character(1)
  )
  serialized <- lapply(keys, function(key) {
    reticulate::py_to_r(
      bridge$extract_estimation_results(reticulate::py_get_item(native_results, key))
    )
  })
  names(serialized) <- keys
  list(
    results = serialized,
    summary_rows = nrow(native_summary),
    statistics = vapply(
      reticulate::iterate(post_processing$pareto$statistics()),
      as.character,
      character(1)
    )
  )
}

build_b21c_test_model <- function(database) {
  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)
  gender <- biogeme_database_segmentation(
    database,
    "MALE",
    c(`0` = "female", `1` = "male")
  )
  ga <- biogeme_database_segmentation(
    database,
    "GA",
    c(`1` = "GA", `0` = "noGA"),
    reference = "noGA"
  )
  income <- biogeme_database_segmentation(
    database,
    "INCOME",
    c(
      `0` = "inc-zero",
      `1` = "inc-under50",
      `2` = "inc-50-100",
      `3` = "inc-100+",
      `4` = "inc-unknown"
    )
  )
  asc_catalogs <- segmentation_catalogs(
    "asc",
    list(asc_car, asc_train),
    list(gender, ga),
    maximum_number = 2
  )
  b_cost_catalog <- segmentation_catalogs(
    "b_cost",
    list(b_cost),
    list(ga, income),
    maximum_number = 1
  )[[1L]]
  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(
    "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(
    "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(
    "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 +
      b_cost_catalog * variable("TRAIN_COST_SCALED"),
    `2` = b_time * sm_tt + b_cost_catalog * variable("SM_COST_SCALED"),
    `3` = asc_catalogs[[1L]] + b_time * car_tt +
      b_cost_catalog * variable("CAR_CO_SCALED")
  )
  biogeme_model(
    database = database,
    formula = logit_log_probability(
      utilities = utilities,
      availability = list(
        `1` = variable("TRAIN_AV_SP"),
        `2` = variable("SM_AV"),
        `3` = variable("CAR_AV_SP")
      ),
      alternative = variable("CHOICE")
    ),
    control = biogeme_control(
      model_name = "b21_multiple_models",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

test_that("b21c Swissmetro Pareto post-processing matches native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run full Swissmetro equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )
  data_path <- rbiogeme_test_swissmetro_path()
  skip_if(!nzchar(data_path), "Set RBIOGEME_SWISSMETRO_DATA to the Swissmetro .dat file")

  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  temporary_directory <- tempfile("rbiogeme-b21c-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  pareto_file <- file.path(getwd(), "b21_multiple_models.pareto")
  native <- native_b21c_specification(data_path, pareto_file)
  database <- swissmetro_data(data)
  model <- build_b21c_test_model(database)
  plot_file <- file.path(getwd(), "b21_process_pareto.png")
  r_fit <- pareto_post_processing(
    model,
    pareto_file_name = pareto_file,
    model_name = "b21_multiple_models",
    control = model$control,
    recycle = FALSE,
    plot_file_name = plot_file,
    objective_x = 0L,
    objective_y = 1L,
    label_x = "Negative loglikelihood",
    label_y = "Number of parameters"
  )

  expect_equal(length(r_fit$results), length(native$results))
  expect_equal(nrow(r_fit$summary), native$summary_rows)
  expect_true(file.exists(plot_file))
  expect_equal(r_fit$pareto_statistics, native$statistics)
  expect_setequal(names(r_fit$results), names(native$results))
  for (configuration in names(native$results)) {
    r_result <- r_fit$results[[configuration]]
    native_result <- native$results[[configuration]]
    # The GA-only cost/Box-Cox specification has a nearly flat optimizer
    # direction. Native repeated re-estimation can stop at points differing
    # by a few 1e-3 while retaining the same likelihood to floating precision.
    coefficient_tolerance <- if (identical(
      configuration,
      "asc:MALE-GA;b_cost:GA;train_tt:boxcox"
    )) 5e-3 else 1e-7
    expect_identical(r_result$beta_names, native_result$beta_names, info = configuration)
    expect_equal(
      unname(coef(r_result)),
      native_result$beta_values,
      tolerance = coefficient_tolerance,
      info = configuration
    )
    expect_equal(
      as.numeric(logLik(r_result)),
      native_result$final_log_likelihood,
      tolerance = 1e-7,
      info = configuration
    )
  }
})