tests/testthat/test-swissmetro-b05c.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b05c <- function(
    data,
    beta_values,
    number_of_draws = 128L,
    seed = 1223L
) {
  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_b05c",
    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")
  utilities <- reticulate::dict(
    `1` = asc_train + b_time_rnd * train_tt_scaled + b_cost * train_cost_scaled,
    `2` = asc_sm + b_time_rnd * sm_tt_scaled + b_cost * sm_cost_scaled,
    `3` = asc_car + b_time_rnd * car_tt_scaled + b_cost * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  conditional_probability <- models$logit(utilities, availability, choice)
  integral <- expressions$MonteCarlo(conditional_probability)
  integral_square <- expressions$MonteCarlo(
    conditional_probability * conditional_probability
  )
  integration_error <- (
    (integral_square - integral * integral) / 2.0
  )^0.5
  simulations <- reticulate::dict(
    Numerator = expressions$MonteCarlo(b_time_rnd * conditional_probability),
    Denominator = integral,
    Integral = integral,
    `Integration error` = integration_error
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    simulations,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- "b05normal_mixture_simul"
  simulated <- biogeme$simulate(
    the_beta_values = reticulate::r_to_py(as.list(beta_values))
  )
  list(
    values = reticulate::py_to_r(simulated),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe)),
    number_of_draws = as.integer(number_of_draws)
  )
}

test_that("b05c Swissmetro simulation 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)
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  b_cost <- biogeme_beta("b_cost", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1)
  b_time_rnd <- b_time + b_time_s * draw("b_time_rnd", "NORMAL")
  utilities <- list(
    `1` = asc_train + b_time_rnd * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    `2` = asc_sm + b_time_rnd * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    `3` = asc_car + b_time_rnd * 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")
  )
  conditional_probability <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )
  integral <- monte_carlo(conditional_probability)
  integral_square <- monte_carlo(conditional_probability * conditional_probability)
  simulations <- list(
    Numerator = monte_carlo(b_time_rnd * conditional_probability),
    Denominator = integral,
    Integral = integral,
    `Integration error` = sqrt((integral_square - integral * integral) / 2.0)
  )
  draws <- biogeme_draws(
    name = "b_time_rnd",
    draw_type = "NORMAL",
    number_of_draws = 128L,
    seed = 1223L
  )
  model <- biogeme_model(database = database, formula = log(integral), draws = draws)
  beta_values <- c(
    asc_train = -0.4,
    b_time = -2.0,
    b_time_s = 1.5,
    b_cost = -1.2,
    asc_car = 0.1
  )
  temporary_directory <- tempfile("rbiogeme-b05c-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_simulation <- simulate(
    model,
    expressions = simulations,
    beta = beta_values,
    control = biogeme_control(
      model_name = "b05normal_mixture_simul_r",
      number_of_draws = 128L,
      seed = 1223L,
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  native <- native_swissmetro_b05c(
    data,
    beta_values = beta_values,
    number_of_draws = 128L,
    seed = 1223L
  )
  native_values <- as.data.frame(native$values, check.names = FALSE)
  r_values <- as.data.frame(r_simulation, check.names = FALSE)

  expect_equal(nrow(r_values), native$number_of_rows)
  expect_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
  expect_equal(r_simulation$number_of_draws, native$number_of_draws)
})