test-montecarlo-group5.R
In rbiogeme: Interface to 'Biogeme'

native_montecarlo_b07 <- function(
    data,
    model_name,
    draw_type,
    number_of_draws,
    seed
) {
  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(
    paste0("native_", model_name),
    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)
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_time_s <- beta("b_time_s", 1, NULL, NULL, 0)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  b_time_random <- b_time + b_time_s *
    expressions$Draws("b_time_rnd", draw_type)
  utilities <- reticulate::dict(
    `1` = asc_train + b_time_random * train_tt_scaled +
      b_cost * train_cost_scaled,
    `2` = b_time_random * sm_tt_scaled + b_cost * sm_cost_scaled,
    `3` = asc_car + b_time_random * car_tt_scaled +
      b_cost * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  probability <- models$logit(utilities, availability, choice)
  log_probability <- expressions$log(
    expressions$MonteCarlo(probability)
  )
  biogeme <- biogeme_module$BIOGEME(
    database,
    log_probability,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    analytical_hessian_mode = "automatic",
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  biogeme$model_name <- model_name
  results <- biogeme$estimate()
  bridge <- rbiogeme:::biogeme_bridge()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

r_montecarlo_b07 <- function(
    data,
    model_name,
    draw_type,
    number_of_draws,
    seed
) {
  database <- swissmetro_data(data)
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  b_time <- biogeme_beta("b_time", start = 0)
  b_time_s <- biogeme_beta("b_time_s", start = 1)
  b_cost <- biogeme_beta("b_cost", start = 0)
  b_time_random <- b_time + b_time_s * draw("b_time_rnd", draw_type)
  utilities <- list(
    "1" = asc_train + b_time_random * variable("TRAIN_TT_SCALED") +
      b_cost * variable("TRAIN_COST_SCALED"),
    "2" = b_time_random * variable("SM_TT_SCALED") +
      b_cost * variable("SM_COST_SCALED"),
    "3" = asc_car + b_time_random * 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")
  )
  probability <- logit_probability(
    utilities = utilities,
    availability = availability,
    alternative = variable("CHOICE")
  )
  model <- biogeme_model(
    database = database,
    formula = log(monte_carlo(probability)),
    draws = biogeme_draws(
      name = "b_time_rnd",
      draw_type = draw_type,
      number_of_draws = number_of_draws,
      seed = seed
    ),
    control = biogeme_control(
      model_name = model_name,
      number_of_draws = number_of_draws,
      seed = seed,
      analytical_hessian_mode = "automatic",
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  estimate(model, model_name = model_name, control = model$control)
}

test_that("Monte Carlo Group 5 files are syntactically valid", {
  files <- rbiogeme_example_path( "montecarlo",
    c(
      "plot_b07estimation_monte_carlo.R",
      "plot_b07estimation_monte_carlo_500.R",
      "plot_b07estimation_monte_carlo_anti.R",
      "plot_b07estimation_monte_carlo_anti_500.R",
      "plot_b07estimation_monte_carlo_halton.R",
      "plot_b07estimation_monte_carlo_halton_500.R",
      "plot_b07estimation_monte_carlo_mlhs.R",
      "plot_b07estimation_monte_carlo_mlhs_500.R",
      "plot_b07estimation_monte_carlo_mlhs_anti.R",
      "plot_b07estimation_monte_carlo_mlhs_anti_500.R"
    )
  )
  expect_true(all(file.exists(files)))
  for (file in files) parse(file)
})

test_that("b07 draw-design estimations match native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run native Monte Carlo equivalence tests"
  )
  skip_if_not(
    rbiogeme_test_configure_python(),
    "Set RBIOGEME_PYTHON to a compatible native Biogeme interpreter"
  )
  data_path <- normalizePath(
    rbiogeme_example_path( "montecarlo", "swissmetro.dat"),
    mustWork = TRUE
  )
  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  variants <- data.frame(
    model_name = c(
      "b07estimation_monte_carlo",
      "b07estimation_monte_carlo_500",
      "b07estimation_monte_carlo_anti",
      "b07estimation_monte_carlo_anti_500",
      "b07estimation_monte_carlo_halton",
      "b07estimation_monte_carlo_halton_500",
      "b07estimation_monte_carlo_mlhs",
      "b07estimation_monte_carlo_mlhs_500",
      "b07estimation_monte_carlo_mlhs_anti",
      "b07estimation_monte_carlo_mlhs_anti_500"
    ),
    draw_type = c(
      "NORMAL", "NORMAL",
      "NORMAL_ANTI", "NORMAL_ANTI",
      "NORMAL_HALTON2", "NORMAL_HALTON2",
      "NORMAL_MLHS", "NORMAL_MLHS",
      "NORMAL_MLHS_ANTI", "NORMAL_MLHS_ANTI"
    ),
    stringsAsFactors = FALSE
  )
  temporary_directory <- tempfile("rbiogeme-montecarlo-group5-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  for (index in seq_len(nrow(variants))) {
    number_of_draws <- 64L
    seed <- 1223L
    variant <- variants[index, , drop = FALSE]
    r_fit <- r_montecarlo_b07(
      data,
      variant$model_name,
      variant$draw_type,
      number_of_draws,
      seed
    )
    native <- native_montecarlo_b07(
      data,
      paste0(variant$model_name, "_native"),
      variant$draw_type,
      number_of_draws,
      seed
    )
    native_results <- native$results

    expect_equal(nobs(r_fit), native$number_of_rows, info = variant$model_name)
    expect_identical(r_fit$beta_names, native_results$beta_names)
    expect_equal(
      unname(coef(r_fit)),
      native_results$beta_values,
      tolerance = 1e-7,
      info = variant$model_name
    )
    expect_equal(
      as.numeric(logLik(r_fit)),
      native_results$final_log_likelihood,
      tolerance = 1e-7,
      info = variant$model_name
    )
    expect_equal(r_fit$number_of_draws, native_results$number_of_draws)
    expect_equal(
      r_fit$number_of_excluded_data,
      native_results$number_of_excluded_data
    )
    expect_identical(isTRUE(r_fit$convergence), isTRUE(native_results$convergence))
  }
})