tests/testthat/test-swissmetro-b21a.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b21a <- function(data, pareto_file_name) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  catalog_module <- reticulate::import("biogeme.catalog", convert = FALSE)
  segmentation_module <- reticulate::import("biogeme.segmentation", convert = FALSE)
  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)
  models <- reticulate::import("biogeme.models", convert = FALSE)
  bridge <- rbiogeme:::biogeme_bridge()

  specification <- reticulate::import(
    "biogeme.catalog.specification",
    convert = FALSE
  )$Specification
  specification$all_results <- reticulate::dict()
  specification$model_names <- NULL

  database <- database_module$Database(
    "swissmetro_native_b21a",
    reticulate::r_to_py(data)
  )
  variable <- expressions$Variable
  purpose <- variable("PURPOSE")
  choice <- variable("CHOICE")
  database$remove(((purpose != 1) * (purpose != 3) + (choice == 0)) > 0)
  ga <- variable("GA")
  sp <- variable("SP")
  sm_cost <- database$define_variable("SM_COST", variable("SM_CO") * (ga == 0))
  train_cost <- database$define_variable("TRAIN_COST", variable("TRAIN_CO") * (ga == 0))
  car_av_sp <- database$define_variable("CAR_AV_SP", variable("CAR_AV") * (sp != 0))
  train_av_sp <- database$define_variable("TRAIN_AV_SP", variable("TRAIN_AV") * (sp != 0))
  train_tt_scaled <- database$define_variable("TRAIN_TT_SCALED", variable("TRAIN_TT") / 100)
  train_cost_scaled <- database$define_variable("TRAIN_COST_SCALED", train_cost / 100)
  sm_tt_scaled <- database$define_variable("SM_TT_SCALED", variable("SM_TT") / 100)
  sm_cost_scaled <- database$define_variable("SM_COST_SCALED", sm_cost / 100)
  car_tt_scaled <- database$define_variable("CAR_TT_SCALED", variable("CAR_TT") / 100)
  car_co_scaled <- database$define_variable("CAR_CO_SCALED", variable("CAR_CO") / 100)

  gender <- database$generate_segmentation(
    variable = variable("MALE"),
    mapping = reticulate::dict(`0` = "female", `1` = "male")
  )
  ga_segmentation <- database$generate_segmentation(
    variable = ga,
    mapping = reticulate::dict(`1` = "GA", `0` = "noGA")
  )
  income <- database$generate_segmentation(
    variable = variable("INCOME"),
    mapping = reticulate::dict(
      `0` = "inc-zero",
      `1` = "inc-under50",
      `2` = "inc-50-100",
      `3` = "inc-100+",
      `4` = "inc-unknown"
    )
  )
  beta <- expressions$Beta
  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_train <- beta("asc_train", 0, NULL, NULL, 0)
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  asc_controller <- catalog_module$Controller(
    controller_name = "asc",
    specification_names = reticulate::r_to_py(c("no_seg", "GA", "MALE", "MALE-GA"))
  )
  asc_options <- function(the_beta) reticulate::dict(
    no_seg = the_beta,
    GA = segmentation_module$Segmentation(the_beta, list(ga_segmentation))$segmented_beta(),
    MALE = segmentation_module$Segmentation(the_beta, list(gender))$segmented_beta(),
    `MALE-GA` = segmentation_module$Segmentation(
      the_beta,
      list(gender, ga_segmentation)
    )$segmented_beta()
  )
  asc_car_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "segmented_asc_car",
    dict_of_expressions = asc_options(asc_car),
    controlled_by = asc_controller
  )
  asc_train_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "segmented_asc_train",
    dict_of_expressions = asc_options(asc_train),
    controlled_by = asc_controller
  )
  cost_controller <- catalog_module$Controller(
    controller_name = "b_cost",
    specification_names = reticulate::r_to_py(c("no_seg", "INCOME", "GA"))
  )
  b_cost_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "segmented_b_cost",
    dict_of_expressions = reticulate::dict(
      no_seg = b_cost,
      INCOME = segmentation_module$Segmentation(b_cost, list(income))$segmented_beta(),
      GA = segmentation_module$Segmentation(b_cost, list(ga_segmentation))$segmented_beta()
    ),
    controlled_by = cost_controller
  )
  lambda_time <- beta("lambda_time", 1, -10, 10, 0)
  time_controller <- catalog_module$Controller(
    controller_name = "train_tt",
    specification_names = reticulate::r_to_py(c("linear", "log", "boxcox"))
  )
  time_options <- function(the_time) reticulate::dict(
    linear = the_time,
    log = expressions$logzero(the_time),
    boxcox = models$boxcox(the_time, lambda_time)
  )
  train_tt_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "train_tt",
    dict_of_expressions = time_options(train_tt_scaled),
    controlled_by = time_controller
  )
  sm_tt_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "sm_tt",
    dict_of_expressions = time_options(sm_tt_scaled),
    controlled_by = time_controller
  )
  car_tt_catalog <- catalog_module$Catalog$from_dict(
    catalog_name = "car_tt",
    dict_of_expressions = time_options(car_tt_scaled),
    controlled_by = time_controller
  )
  utilities <- reticulate::dict(
    `1` = asc_train_catalog + b_time * train_tt_catalog + b_cost_catalog * train_cost_scaled,
    `2` = b_time * sm_tt_catalog + b_cost_catalog * sm_cost_scaled,
    `3` = asc_car_catalog + b_time * car_tt_catalog + b_cost_catalog * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  estimator <- biogeme_module$BIOGEME(
    database,
    models$loglogit(utilities, availability, choice),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  estimator$model_name <- "b21_multiple_models"
  assisted <- assisted_module$AssistedSpecification(
    biogeme_object = estimator,
    multi_objectives = objectives_module$loglikelihood_dimension,
    pareto_file_name = pareto_file_name
  )
  results <- assisted$run()
  compiled <- results_processing$compile_estimation_results(
    results,
    use_short_names = TRUE
  )
  native_summary <- reticulate::py_to_r(reticulate::py_get_item(compiled, 0L))
  keys <- vapply(reticulate::iterate(results$keys()), as.character, character(1))
  serialized <- lapply(keys, function(key) {
    reticulate::py_to_r(bridge$extract_estimation_results(reticulate::py_get_item(results, key)))
  })
  names(serialized) <- keys
  list(
    results = serialized,
    summary_rows = nrow(native_summary),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("b21a Swissmetro assisted specification 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-b21a-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  native <- native_swissmetro_b21a(data, file.path(getwd(), "native.pareto"))
  database <- swissmetro_data(data)
  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"
    )
  )
  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)
  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]]
  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
  )
  model <- biogeme_model(
    database,
    formula = logit_log_probability(
      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")
      ),
      alternative = variable("CHOICE")
    )
  )
  control <- biogeme_control(
    model_name = "b21_multiple_models",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- assisted_specification(
    model,
    objectives = "loglikelihood_dimension",
    pareto_file_name = file.path(getwd(), "r.pareto"),
    model_name = "b21_multiple_models",
    control = control,
    force = TRUE
  )

  expect_equal(length(r_fit$results), length(native$results))
  expect_equal(nrow(r_fit$summary), native$summary_rows)
  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]]
    expect_identical(r_result$beta_names, native_result$beta_names)
    expect_equal(unname(coef(r_result)), native_result$beta_values, tolerance = 1e-7)
    expect_equal(
      as.numeric(logLik(r_result)),
      native_result$final_log_likelihood,
      tolerance = 1e-7
    )
  }
})