test-swissmetro-b01d.R
In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b01d <- 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)

  database <- database_module$Database(
    "swissmetro_native_b01d",
    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_cost_scaled <- database$define_variable("TRAIN_COST_SCALED", train_cost / 100)
  sm_cost_scaled <- database$define_variable("SM_COST_SCALED", sm_cost / 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)
  train_tt <- variable("TRAIN_TT")
  sm_tt <- variable("SM_TT")
  car_tt <- variable("CAR_TT")
  utilities <- reticulate::dict(
    `1` = asc_train + b_time * train_tt / 100 + b_cost * train_cost_scaled,
    `2` = asc_sm + b_time * sm_tt / 100 + b_cost * sm_cost_scaled,
    `3` = asc_car + b_time * car_tt / 100 + b_cost * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  log_probability <- models$loglogit(utilities, availability, choice)
  estimator <- biogeme_module$BIOGEME(
    database,
    log_probability,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  estimator$model_name <- "b01a_native_for_b01d"
  estimator$calculate_null_loglikelihood(availability)
  estimation_results <- estimator$estimate()
  betas <- estimation_results$get_beta_values()

  prob_train <- models$logit(utilities, availability, 1L)
  prob_swissmetro <- models$logit(utilities, availability, 2L)
  prob_car <- models$logit(utilities, availability, 3L)
  simulation <- reticulate::dict(
    `Prob. train` = prob_train,
    `Prob. Swissmetro` = prob_swissmetro,
    `Prob. car` = prob_car,
    `logit elas. 1` = train_av_sp * (1 - prob_train) * train_tt * b_time / 100,
    `generic elas. 1` = expressions$Derive(prob_train, "TRAIN_TT") * train_tt / prob_train,
    `logit elas. 2` = variable("SM_AV") * (1 - prob_swissmetro) * sm_tt * b_time / 100,
    `generic elas. 2` = expressions$Derive(prob_swissmetro, "SM_TT") * sm_tt / prob_swissmetro,
    `logit elas. 3` = car_av_sp * (1 - prob_car) * car_tt * b_time / 100,
    `generic elas. 3` = expressions$Derive(prob_car, "CAR_TT") * car_tt / prob_car
  )
  simulator <- biogeme_module$BIOGEME(
    database,
    simulation,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  list(
    estimates = unlist(reticulate::py_to_r(betas), use.names = TRUE),
    values = reticulate::py_to_r(simulator$simulate(the_beta_values = betas))
  )
}

test_that("b01d 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)
  model <- local({
    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 <- biogeme_beta("b_time", start = 0)
    b_cost <- biogeme_beta("b_cost", start = 0)
    utilities <- list(
      `1` = asc_train + b_time * variable("TRAIN_TT") / 100 +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm + b_time * variable("SM_TT") / 100 +
        b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * 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")
    )
    probabilities <- list(
      `Prob. train` = logit_probability(utilities, availability, 1),
      `Prob. Swissmetro` = logit_probability(utilities, availability, 2),
      `Prob. car` = logit_probability(utilities, availability, 3)
    )
    result <- logit_model(database, "CHOICE", utilities, availability)
    result$simulations <- probabilities
    result
  })
  temporary_directory <- tempfile("rbiogeme-b01d-")
  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 = "b01a_r_for_b01d",
    control = biogeme_control(
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )

  model$simulations <- list(
    `Prob. train` = logit_probability(
      model$utilities, model$availability, 1
    ),
    `Prob. Swissmetro` = logit_probability(
      model$utilities, model$availability, 2
    ),
    `Prob. car` = logit_probability(
      model$utilities, model$availability, 3
    ),
    `logit elas. 1` = variable("TRAIN_AV_SP") * (1 - model$simulations[[1]]) *
      variable("TRAIN_TT") * model$parameters$b_time / 100,
    `generic elas. 1` = Derive(model$simulations[[1]], "TRAIN_TT") *
      variable("TRAIN_TT") / model$simulations[[1]],
    `logit elas. 2` = variable("SM_AV") * (1 - model$simulations[[2]]) *
      variable("SM_TT") * model$parameters$b_time / 100,
    `generic elas. 2` = Derive(model$simulations[[2]], "SM_TT") *
      variable("SM_TT") / model$simulations[[2]],
    `logit elas. 3` = variable("CAR_AV_SP") * (1 - model$simulations[[3]]) *
      variable("CAR_TT") * model$parameters$b_time / 100,
    `generic elas. 3` = Derive(model$simulations[[3]], "CAR_TT") *
      variable("CAR_TT") / model$simulations[[3]]
  )
  r_simulation <- simulate(model, beta = r_fit)
  native <- native_swissmetro_b01d(data)

  expect_equal(unname(coef(r_fit)), unname(native$estimates), tolerance = 1e-8)
  expect_equal(nrow(r_simulation$values), 6768L)
  expect_identical(names(r_simulation$values), names(native$values))
  expect_equal(
    unname(as.matrix(r_simulation$values)),
    unname(as.matrix(native$values)),
    tolerance = 1e-8,
    ignore_attr = TRUE
  )
  expect_equal(
    r_simulation$values$`logit elas. 1`,
    r_simulation$values$`generic elas. 1`,
    tolerance = 1e-10
  )
  expect_equal(
    r_simulation$values$`logit elas. 2`,
    r_simulation$values$`generic elas. 2`,
    tolerance = 1e-10
  )
  finite_car <- is.finite(r_simulation$values$`logit elas. 3`) &
    is.finite(r_simulation$values$`generic elas. 3`)
  expect_equal(
    r_simulation$values$`logit elas. 3`[finite_car],
    r_simulation$values$`generic elas. 3`[finite_car],
    tolerance = 1e-10
  )
})