tests/testthat/test-assisted-group1.R

In rbiogeme: Interface to 'Biogeme'

test_that("assisted group 1 examples are self-contained", {
  scripts <- c(
    "plot_simple_example.R",
    "plot_b00logit.R"
  )
  paths <- file.path(rbiogeme_example_path( "assisted"), scripts)
  expect_true(all(file.exists(paths)))
  for (path in paths) {
    expect_error(parse(file = path), NA)
    source_text <- paste(readLines(path, warn = FALSE), collapse = "\n")
    expect_match(source_text, "biogeme_beta\\(")
    expect_match(source_text, "prepare_swissmetro_example")
    expect_match(source_text, "#")
  }

  b00 <- paste(
    readLines(paths[[2L]], warn = FALSE),
    collapse = "\n"
  )
  expect_match(b00, "build_b00logit_model")
  expect_match(b00, "filter_purpose = FALSE")
  expect_match(b00, "model_name = \"b00logit\"")
})

native_assisted_b00 <- function(data) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)

  database <- database_module$Database(
    "swissmetro_native_assisted_b00",
    reticulate::r_to_py(data)
  )
  variable <- expressions$Variable
  choice <- variable("CHOICE")
  database$remove(choice == 0)
  ga <- variable("GA")
  sp <- variable("SP")
  train_cost <- database$define_variable(
    "TRAIN_COST",
    variable("TRAIN_CO") * (ga == 0)
  )
  sm_cost <- database$define_variable(
    "SM_COST",
    variable("SM_CO") * (ga == 0)
  )
  train_av_sp <- database$define_variable(
    "TRAIN_AV_SP",
    variable("TRAIN_AV") * (sp != 0)
  )
  car_av_sp <- database$define_variable(
    "CAR_AV_SP",
    variable("CAR_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
  )
  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)
  utilities <- reticulate::dict(
    `1` = asc_train + b_time * train_tt_scaled + b_cost * train_cost_scaled,
    `2` = b_time * sm_tt_scaled + b_cost * sm_cost_scaled,
    `3` = asc_car + b_time * car_tt_scaled + b_cost * 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 <- "b00logit_native_assisted"
  estimator$calculate_null_loglikelihood(availability)
  result <- estimator$estimate()
  bridge <- rbiogeme:::biogeme_bridge()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(result)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("assisted b00logit matches native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run full assisted 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)
  database <- swissmetro_data(data, filter_purpose = FALSE)
  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)
  model <- logit_model(
    database = database,
    choice = "CHOICE",
    utilities = list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = 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")
    )
  )

  temporary_directory <- tempfile("rbiogeme-assisted-b00-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)
  fit <- estimate(
    model,
    model_name = "b00logit_r",
    control = biogeme_control(
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  native <- native_assisted_b00(data)

  expect_equal(nobs(fit), native$number_of_rows)
  expect_identical(fit$beta_names, native$results$beta_names)
  expect_equal(unname(coef(fit)), native$results$beta_values, tolerance = 1e-8)
  expect_equal(as.numeric(logLik(fit)), native$results$final_log_likelihood, tolerance = 1e-8)
  expect_equal(fit$number_of_excluded_data, native$results$number_of_excluded_data)
  expect_identical(isTRUE(fit$convergence), isTRUE(native$results$convergence))
})