tests/testthat/test-swissmetro-b27.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b27 <- function(data, output_directory) {
  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_b27", 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_cost <- beta("b_cost", 0, NULL, NULL, 0)
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_time_s <- beta("b_time_s", 1, NULL, NULL, 0)
  b_time_rnd <- b_time + b_time_s * expressions$Draws("b_time_rnd", "NORMAL")
  v_train <- asc_train + b_time_rnd * train_tt_scaled + b_cost * train_cost_scaled
  v_swissmetro <- asc_sm + b_time_rnd * sm_tt_scaled + b_cost * sm_cost_scaled
  v_car <- asc_car + b_time_rnd * car_tt_scaled + b_cost * car_co_scaled
  log_probability <- expressions$log(
    expressions$MonteCarlo(models$logit(
      reticulate::dict(`1` = v_train, `2` = v_swissmetro, `3` = v_car),
      reticulate::dict(`1` = train_av_sp, `2` = variable("SM_AV"), `3` = car_av_sp),
      choice
    ))
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    number_of_draws = 2000L,
    seed = 1223L,
    calculating_second_derivatives = "never",
    monte_carlo_diagnostic_auto = FALSE,
    monte_carlo_diagnostic_draw_factors = "0.5,1.0,2.0",
    monte_carlo_diagnostic_replications = 1L,
    monte_carlo_diagnostic_time_budget = 300L,
    monte_carlo_diagnostic_max_draws = 4000L,
    user_notes = paste0(
      "Post-estimation Monte Carlo draw-stability diagnostic for a mixed ",
      "logit model using the Swissmetro data."
    ),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b27_monte_carlo_native"
  results <- biogeme$estimate()
  diagnostic <- biogeme$check_monte_carlo_stability(
    estimation_results = results,
    output_directory = output_directory,
    basename = "b27",
    resume = FALSE
  )
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    diagnostic = reticulate::py_to_r(diagnostic$data),
    execution_status = as.character(diagnostic$execution_status),
    diagnostic_conclusion = as.character(diagnostic$diagnostic_conclusion),
    recommendation = as.character(diagnostic$recommendation)
  )
}

normalise_diagnostic_seed <- function(value) {
  value <- as.numeric(value)
  if (value < 0) value <- value + 4294967296
  value
}

test_that("b27 Swissmetro Monte Carlo diagnostic 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)
  r_model <- r_swissmetro_mixture(data, "b27", number_of_draws = 2000L)$model
  r_model$control <- biogeme_control(
    model_name = "b27_monte_carlo",
    number_of_draws = 2000,
    seed = 1223,
    calculating_second_derivatives = "never",
    monte_carlo_diagnostic_auto = FALSE,
    monte_carlo_diagnostic_draw_factors = "0.5,1.0,2.0",
    monte_carlo_diagnostic_replications = 1,
    monte_carlo_diagnostic_time_budget = 300,
    monte_carlo_diagnostic_max_draws = 4000,
    user_notes = paste0(
      "Post-estimation Monte Carlo draw-stability diagnostic for a mixed ",
      "logit model using the Swissmetro data."
    ),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  temporary_directory <- tempfile("rbiogeme-b27-")
  dir.create(temporary_directory, recursive = TRUE)
  native_directory <- file.path(temporary_directory, "native")
  r_directory <- file.path(temporary_directory, "r")
  dir.create(native_directory)
  dir.create(r_directory)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_fit <- estimate(r_model, model_name = "b27_monte_carlo", control = r_model$control)
  native <- native_swissmetro_b27(data, native_directory)
  r_diagnostic <- check_monte_carlo_stability(
    model = r_model,
    fit = r_fit,
    model_name = "b27_monte_carlo",
    control = r_model$control,
    output_directory = r_directory,
    basename = "b27",
    resume = FALSE
  )

  expect_identical(r_diagnostic$execution_status, native$execution_status)
  expect_identical(r_diagnostic$diagnostic_conclusion, native$diagnostic_conclusion)
  expect_identical(r_diagnostic$recommendation, native$recommendation)
  expect_length(r_diagnostic$data$planned_evaluations, 3L)
  expect_length(native$diagnostic$planned_evaluations, 3L)
  for (index in seq_len(3L)) {
    r_plan <- r_diagnostic$data$planned_evaluations[[index]]
    native_plan <- native$diagnostic$planned_evaluations[[index]]
    expect_equal(r_plan$draw_count, native_plan$draw_count)
    expect_equal(r_plan$draw_factor, native_plan$draw_factor)
    expect_equal(r_plan$replication, native_plan$replication)
    expect_equal(
      normalise_diagnostic_seed(r_plan$seed),
      normalise_diagnostic_seed(native_plan$seed)
    )
  }
  expect_length(r_diagnostic$data$completed_evaluations, 3L)
  expect_length(native$diagnostic$completed_evaluations, 3L)
  for (index in seq_len(3L)) {
    r_evaluation <- r_diagnostic$data$completed_evaluations[[index]]
    native_evaluation <- native$diagnostic$completed_evaluations[[index]]
    expect_equal(r_evaluation$draw_count, native_evaluation$draw_count)
    expect_equal(
      normalise_diagnostic_seed(r_evaluation$seed),
      normalise_diagnostic_seed(native_evaluation$seed)
    )
    expect_equal(r_evaluation$objective, native_evaluation$objective, tolerance = 2e-6)
    expect_equal(r_evaluation$gradient, native_evaluation$gradient, tolerance = 2e-6)
  }
})