tests/testthat/test-swissmetro-b04.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b04 <- function(data, seed = 73129L, slices = 5L) {
  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)
  numpy <- reticulate::import("numpy", convert = FALSE)

  database <- database_module$Database(
    "swissmetro_native_b04",
    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)

  beta <- expressions$Beta
  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_train <- beta("asc_train", 0, NULL, NULL, 0)
  asc_sm <- beta("asc_sm", 0, NULL, NULL, 1)
  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` = asc_sm + 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
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    models$loglogit(utilities, availability, choice),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b04_native"
  biogeme$calculate_null_loglikelihood(availability)
  results <- biogeme$estimate()

  # The native split uses NumPy's global random state. Seed it immediately
  # before validate(), matching the bridge's temporary seeded validation call.
  numpy$random$seed(as.integer(seed))
  validation <- biogeme$validate(results, slices = as.integer(slices))
  folds <- lapply(seq_len(reticulate::py_len(validation)), function(index) {
    fold <- reticulate::py_get_item(validation, index - 1L)
    simulated_values <- reticulate::py_to_r(fold$simulated_values)
    simulated_values <- as.data.frame(simulated_values, check.names = FALSE)
    estimation_data <- fold$estimation_modeling_elements$database$dataframe
    validation_data <- fold$validation_modeling_elements$database$dataframe
    list(
      simulated_values = simulated_values,
      estimation_sample_size = reticulate::py_len(estimation_data),
      validation_sample_size = reticulate::py_len(validation_data),
      estimation_indices = as.integer(reticulate::py_to_r(estimation_data$index$tolist())),
      validation_indices = as.integer(reticulate::py_to_r(validation_data$index$tolist()))
    )
  })
  bridge <- rbiogeme:::biogeme_bridge()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe)),
    folds = folds
  )
}

test_that("b04 Swissmetro validation 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)
  database <- swissmetro_data(data)
  model <- swissmetro_mnl_model(database)
  temporary_directory <- tempfile("rbiogeme-b04-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  control <- biogeme_control(
    model_name = "b04_r",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  r_fit <- estimate(model, model_name = "b04_r", control = control)
  seed <- 73129L
  r_validation <- validate(
    model = model,
    fit = r_fit,
    folds = 5L,
    seed = seed,
    control = control
  )
  native <- native_swissmetro_b04(data, seed = seed, slices = 5L)
  native_results <- native$results

  expect_equal(nobs(r_fit), native$number_of_rows)
  expect_identical(r_fit$beta_names, native_results$beta_names)
  expect_equal(unname(coef(r_fit)), native_results$beta_values, tolerance = 1e-8)
  expect_equal(
    as.numeric(logLik(r_fit)),
    native_results$final_log_likelihood,
    tolerance = 1e-8
  )
  expect_equal(length(r_validation), 5L)
  expect_equal(length(r_validation), length(native$folds))
  for (index in seq_along(r_validation)) {
    r_fold <- r_validation[[index]]
    native_fold <- native$folds[[index]]
    expect_equal(r_fold$fold, index)
    expect_equal(r_fold$estimation_sample_size, native_fold$estimation_sample_size)
    expect_equal(r_fold$validation_sample_size, native_fold$validation_sample_size)
    expect_equal(r_fold$estimation_indices, native_fold$estimation_indices)
    expect_equal(r_fold$validation_indices, native_fold$validation_indices)
    r_values <- as.data.frame(r_fold$simulated_values, check.names = FALSE)
    native_values <- native_fold$simulated_values
    expect_identical(names(r_values), names(native_values))
    expect_equal(
      unname(as.matrix(r_values)),
      unname(as.matrix(native_values)),
      tolerance = 1e-8
    )
  }
})