tests/testthat/test-swissmetro-b11b.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b11b <- 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)
  nests_module <- reticulate::import("biogeme.nests", convert = FALSE)
  bridge <- rbiogeme:::biogeme_bridge()

  database <- database_module$Database(
    "swissmetro_native_b11b",
    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_swissmetro <- beta("b_time_swissmetro", 0, NULL, NULL, 0)
  b_time_train <- beta("b_time_train", 0, NULL, NULL, 0)
  b_time_car <- beta("b_time_car", 0, NULL, NULL, 0)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  b_headway_swissmetro <- beta("b_headway_swissmetro", 0, NULL, NULL, 0)
  b_headway_train <- beta("b_headway_train", 0, NULL, NULL, 0)
  ga_train <- beta("ga_train", 0, NULL, NULL, 0)
  ga_swissmetro <- beta("ga_swissmetro", 0, NULL, NULL, 0)
  existing_nest_parameter <- beta("existing_nest_parameter", 1, 1, 5, 0)
  public_nest_parameter <- beta("public_nest_parameter", 1, 1, 5, 0)
  alpha_existing <- beta("alpha_existing", 0.5, 0, 1, 0)
  alpha_public <- 1 - alpha_existing
  utilities <- reticulate::dict(
    `1` = asc_train + b_time_train * variable("TRAIN_TT") / 100 +
      b_cost * train_cost_scaled + b_headway_train * variable("TRAIN_HE") +
      ga_train * ga,
    `2` = asc_sm + b_time_swissmetro * variable("SM_TT") / 100 +
      b_cost * sm_cost_scaled + b_headway_swissmetro * variable("SM_HE") +
      ga_swissmetro * ga,
    `3` = asc_car + b_time_car * variable("CAR_TT") / 100 +
      b_cost * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  nest_existing <- nests_module$OneNestForCrossNestedLogit(
    nest_param = existing_nest_parameter,
    dict_of_alpha = reticulate::dict(`1` = alpha_existing, `2` = 0.0, `3` = 1.0),
    name = "existing"
  )
  nest_public <- nests_module$OneNestForCrossNestedLogit(
    nest_param = public_nest_parameter,
    dict_of_alpha = reticulate::dict(`1` = alpha_public, `2` = 1.0, `3` = 0.0),
    name = "public"
  )
  nests <- nests_module$NestsForCrossNestedLogit(
    choice_set = reticulate::r_to_py(as.integer(c(1L, 2L, 3L))),
    tuple_of_nests = reticulate::tuple(nest_existing, nest_public)
  )

  log_probability <- models$logcnl(utilities, availability, nests, choice)
  estimator <- biogeme_module$BIOGEME(
    database,
    log_probability,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  estimator$model_name <- "b11a_cnl_native_for_b11b"
  estimation_results <- estimator$estimate()
  betas <- estimation_results$get_beta_values()

  probability_train <- models$cnl(utilities, availability, nests, 1L)
  probability_swissmetro <- models$cnl(utilities, availability, nests, 2L)
  probability_car <- models$cnl(utilities, availability, nests, 3L)
  simulations <- reticulate::dict(
    `Prob. train` = probability_train,
    `Prob. Swissmetro` = probability_swissmetro,
    `Prob. car` = probability_car,
    `Elas. 1` = expressions$Derive(probability_train, "TRAIN_TT") *
      variable("TRAIN_TT") / probability_train,
    `Elas. 2` = expressions$Derive(probability_swissmetro, "SM_TT") *
      variable("SM_TT") / probability_swissmetro,
    `Elas. 3` = expressions$Derive(probability_car, "CAR_TT") *
      variable("CAR_TT") / probability_car
  )
  simulator <- biogeme_module$BIOGEME(
    database,
    simulations,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  simulator$model_name <- "b11b_cnl_simul_native"
  simulated <- simulator$simulate(the_beta_values = betas)
  correlation <- nests$correlation(
    parameters = betas,
    alternatives_names = reticulate::dict(
      `1` = "Train", `2` = "Swissmetro", `3` = "Car"
    )
  )
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(estimation_results)),
    values = reticulate::py_to_r(simulated),
    correlation = reticulate::py_to_r(correlation),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

test_that("b11b Swissmetro CNL 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_time_swissmetro <- biogeme_beta("b_time_swissmetro", start = 0)
  b_time_train <- biogeme_beta("b_time_train", start = 0)
  b_time_car <- biogeme_beta("b_time_car", start = 0)
  b_cost <- biogeme_beta("b_cost", start = 0)
  b_headway_swissmetro <- biogeme_beta("b_headway_swissmetro", start = 0)
  b_headway_train <- biogeme_beta("b_headway_train", start = 0)
  ga_train <- biogeme_beta("ga_train", start = 0)
  ga_swissmetro <- biogeme_beta("ga_swissmetro", start = 0)
  existing_nest_parameter <- biogeme_beta("existing_nest_parameter", start = 1, lower = 1, upper = 5)
  public_nest_parameter <- biogeme_beta("public_nest_parameter", start = 1, lower = 1, upper = 5)
  alpha_existing <- biogeme_beta("alpha_existing", start = 0.5, lower = 0, upper = 1)
  alpha_public <- 1 - alpha_existing
  utilities <- list(
    `1` = asc_train + b_time_train * variable("TRAIN_TT") / 100 +
      b_cost * variable("TRAIN_COST_SCALED") + b_headway_train * variable("TRAIN_HE") +
      ga_train * variable("GA"),
    `2` = asc_sm + b_time_swissmetro * variable("SM_TT") / 100 +
      b_cost * variable("SM_COST_SCALED") + b_headway_swissmetro * variable("SM_HE") +
      ga_swissmetro * variable("GA"),
    `3` = asc_car + b_time_car * variable("CAR_TT") / 100 +
      b_cost * variable("CAR_CO_SCALED")
  )
  availability <- list(
    `1` = variable("TRAIN_AV_SP"),
    `2` = variable("SM_AV"),
    `3` = variable("CAR_AV_SP")
  )
  nests <- cross_nested_nests(
    choice_set = c(1L, 2L, 3L),
    nests = list(
      cross_nested_nest(
        existing_nest_parameter,
        list(`1` = alpha_existing, `2` = 0, `3` = 1),
        name = "existing"
      ),
      cross_nested_nest(
        public_nest_parameter,
        list(`1` = alpha_public, `2` = 1, `3` = 0),
        name = "public"
      )
    )
  )
  probability_train <- cross_nested_probability(utilities, availability, nests, 1)
  probability_swissmetro <- cross_nested_probability(utilities, availability, nests, 2)
  probability_car <- cross_nested_probability(utilities, availability, nests, 3)
  model <- cross_nested_logit_model(
    database = database,
    choice = "CHOICE",
    utilities = utilities,
    availability = availability,
    nests = nests,
    control = biogeme_control(
      model_name = "b11a_cnl",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  model$simulations <- list(
    `Prob. train` = probability_train,
    `Prob. Swissmetro` = probability_swissmetro,
    `Prob. car` = probability_car,
    `Elas. 1` = Derive(probability_train, "TRAIN_TT") * variable("TRAIN_TT") / probability_train,
    `Elas. 2` = Derive(probability_swissmetro, "SM_TT") * variable("SM_TT") / probability_swissmetro,
    `Elas. 3` = Derive(probability_car, "CAR_TT") * variable("CAR_TT") / probability_car
  )

  temporary_directory <- tempfile("rbiogeme-b11b-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_fit <- estimate(model, model_name = "b11a_cnl", control = model$control)
  r_simulation <- simulate(
    model,
    beta = r_fit,
    control = biogeme_control(
      model_name = "b11b_cnl_simul",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  native <- native_swissmetro_b11b(data)
  native_values <- as.data.frame(native$values, check.names = FALSE)
  r_values <- as.data.frame(r_simulation, check.names = FALSE)

  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_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-8,
    ignore_attr = TRUE
  )
  expect_equal(
    unname(as.matrix(cross_nested_logit_correlation(
      model,
      beta_values = coef(r_fit),
      alternatives_names = c(`1` = "Train", `2` = "Swissmetro", `3` = "Car")
    ))),
    unname(as.matrix(native$correlation)),
    tolerance = 1e-10
  )
  expect_equal(
    100 * mean(r_values[["Prob. train"]]),
    100 * mean(native_values[["Prob. train"]]),
    tolerance = 1e-10
  )
})