tests/testthat/test-swissmetro-b01e.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b01e <- function(data, settings) {
  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)
  bridge <- rbiogeme:::biogeme_bridge()

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

  lapply(seq_len(nrow(settings)), function(index) {
    row <- settings[index, , drop = FALSE]
    native <- biogeme_module$BIOGEME(
      database,
      models$loglogit(utilities, availability, choice),
      infeasible_cg = isTRUE(row$infeasible_cg),
      initial_radius = as.numeric(row$initial_radius),
      second_derivatives = as.numeric(row$second_derivatives),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
    native$model_name <- paste0("b05normal_mixture_algo_", row$name)
    native$calculate_null_loglikelihood(availability)
    result <- native$estimate()
    reticulate::py_to_r(bridge$extract_estimation_results(result))
  })
}

test_that("b01e algorithm settings match 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)
  settings <- data.frame(
    infeasible_cg = c(TRUE, FALSE),
    initial_radius = c(0.1, 10.0),
    second_derivatives = c(0.0, 1.0),
    name = c(
      "cg_TRUE_radius_0.1_second_deriv_0.0",
      "cg_FALSE_radius_10.0_second_deriv_1.0"
    ),
    stringsAsFactors = FALSE
  )
  database <- swissmetro_data(data)
  model <- swissmetro_mnl_model(database)
  original_directory <- getwd()
  temporary_directory <- tempfile("rbiogeme-b01e-")
  dir.create(temporary_directory, recursive = TRUE)
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_results <- lapply(seq_len(nrow(settings)), function(index) {
    row <- settings[index, , drop = FALSE]
    estimate(
      model,
      model_name = paste0("b05normal_mixture_algo_", row$name),
      control = biogeme_control(
        infeasible_cg = isTRUE(row$infeasible_cg),
        initial_radius = as.numeric(row$initial_radius),
        second_derivatives_percentage = as.numeric(row$second_derivatives),
        generate_html = FALSE,
        generate_yaml = FALSE,
        save_iterations = FALSE
      )
    )
  })
  native_results <- native_swissmetro_b01e(data, settings)

  expect_length(r_results, 2L)
  for (index in seq_along(r_results)) {
    r_fit <- r_results[[index]]
    native <- native_results[[index]]
    expect_equal(unname(coef(r_fit)), native$beta_values, tolerance = 1e-6)
    expect_equal(
      as.numeric(logLik(r_fit)),
      native$final_log_likelihood,
      tolerance = 1e-8
    )
    expect_equal(r_fit$number_of_excluded_data, native$number_of_excluded_data)
    expect_identical(isTRUE(r_fit$convergence), isTRUE(native$convergence))
    expect_equal(r_fit$gradient_norm, native$gradient_norm, tolerance = 1e-3)
    expect_true(nzchar(r_fit$termination_reason))
    expect_true(nzchar(native$termination_reason))
    expect_match(r_fit$termination_reason, "Relative gradient")
    expect_match(native$termination_reason, "Relative gradient")
  }
})