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

native_montecarlo_b03 <- function(number_of_draws, seed, explicit = FALSE) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  draws_module <- reticulate::import("biogeme.draws", convert = FALSE)
  functools <- reticulate::import("functools", convert = FALSE)

  database <- database_module$Database(
    "native_montecarlo_b03",
    reticulate::r_to_py(data.frame(FakeColumn = 1.0))
  )
  if (isTRUE(explicit)) {
    uniform_draw <- expressions$Draws("U", "UNIFORM")
    integrand <- expressions$exp(uniform_draw) + expressions$exp(1 - uniform_draw)
    simulated_integral <- expressions$MonteCarlo(integrand) / 2.0
    halton13_draw <- expressions$Draws("U_halton13", "HALTON13")
    integrand_halton13 <- expressions$exp(halton13_draw) +
      expressions$exp(1 - halton13_draw)
    simulated_integral_halton13 <- expressions$MonteCarlo(integrand_halton13) / 2.0
    mlhs_draw <- expressions$Draws("U_mlhs", "UNIFORM_MLHS")
    integrand_mlhs <- expressions$exp(mlhs_draw) + expressions$exp(1 - mlhs_draw)
    simulated_integral_mlhs <- expressions$MonteCarlo(integrand_mlhs) / 2.0
    custom_draw_type <- "HALTON13"
    custom_generator <- functools$partial(
      draws_module$get_halton_draws,
      base = 13L,
      skip = 10L
    )
  } else {
    integrand <- expressions$exp(expressions$Draws("U", "UNIFORM_ANTI"))
    simulated_integral <- expressions$MonteCarlo(integrand)
    integrand_halton13 <- expressions$exp(
      expressions$Draws("U_halton13", "HALTON13_ANTI")
    )
    simulated_integral_halton13 <- expressions$MonteCarlo(integrand_halton13)
    integrand_mlhs <- expressions$exp(
      expressions$Draws("U_mlhs", "UNIFORM_MLHS_ANTI")
    )
    simulated_integral_mlhs <- expressions$MonteCarlo(integrand_mlhs)
    custom_draw_type <- "HALTON13_ANTI"
    base_halton13 <- functools$partial(
      draws_module$get_halton_draws,
      base = 13L,
      skip = 10L
    )
    custom_generator <- functools$partial(
      draws_module$get_antithetic,
      base_halton13
    )
  }
  true_integral <- expressions$exp(1.0) - 1.0
  simulations <- reticulate::dict(
    `Analytical Integral` = true_integral,
    `Simulated Integral` = simulated_integral,
    `Error             ` = simulated_integral - true_integral,
    `Simulated Integral (Halton13)` = simulated_integral_halton13,
    `Error (Halton13)             ` = simulated_integral_halton13 - true_integral,
    `Simulated Integral (MLHS)` = simulated_integral_mlhs,
    `Error (MLHS)             ` = simulated_integral_mlhs - true_integral
  )
  custom_generator_tuple <- draws_module$RandomNumberGeneratorTuple(
    generator = custom_generator,
    description = if (isTRUE(explicit)) {
      "Halton draws with base 13, skipping 10"
    } else {
      "Antithetic Halton draws with base 13, skipping 10"
    }
  )
  random_number_generators <- reticulate::dict()
  random_number_generators[[custom_draw_type]] <- custom_generator_tuple
  biogeme <- biogeme_module$BIOGEME(
    database,
    simulations,
    random_number_generators = random_number_generators,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  reticulate::py_to_r(
    biogeme$simulate(the_beta_values = reticulate::dict())
  )
}

r_montecarlo_b03 <- function(number_of_draws, seed, explicit = FALSE) {
  database <- biogeme_database(
    "fake_database",
    data.frame(FakeColumn = 1.0)
  )
  if (isTRUE(explicit)) {
    uniform_draw <- draw("U", "UNIFORM")
    integrand <- exp(uniform_draw) + exp(1 - uniform_draw)
    simulated_integral <- monte_carlo(integrand) / 2.0
    halton13_draw <- draw("U_halton13", "HALTON13")
    integrand_halton13 <- exp(halton13_draw) + exp(1 - halton13_draw)
    simulated_integral_halton13 <- monte_carlo(integrand_halton13) / 2.0
    mlhs_draw <- draw("U_mlhs", "UNIFORM_MLHS")
    integrand_mlhs <- exp(mlhs_draw) + exp(1 - mlhs_draw)
    simulated_integral_mlhs <- monte_carlo(integrand_mlhs) / 2.0
    custom_draws <- biogeme_draws(
      name = "U_halton13",
      draw_type = "HALTON13",
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generator = "HALTON13"
    )
  } else {
    integrand <- exp(draw("U", "UNIFORM_ANTI"))
    simulated_integral <- monte_carlo(integrand)
    integrand_halton13 <- exp(draw("U_halton13", "HALTON13_ANTI"))
    simulated_integral_halton13 <- monte_carlo(integrand_halton13)
    integrand_mlhs <- exp(draw("U_mlhs", "UNIFORM_MLHS_ANTI"))
    simulated_integral_mlhs <- monte_carlo(integrand_mlhs)
    custom_draws <- biogeme_draws(
      name = "U_halton13",
      draw_type = "HALTON13_ANTI",
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generator = "HALTON13_ANTI"
    )
  }
  true_integral <- exp(1.0) - 1.0
  simulations <- list(
    `Analytical Integral` = true_integral,
    `Simulated Integral` = simulated_integral,
    `Error             ` = simulated_integral - true_integral,
    `Simulated Integral (Halton13)` = simulated_integral_halton13,
    `Error (Halton13)             ` = simulated_integral_halton13 - true_integral,
    `Simulated Integral (MLHS)` = simulated_integral_mlhs,
    `Error (MLHS)             ` = simulated_integral_mlhs - true_integral
  )
  model <- biogeme_model(
    database = database,
    simulations = simulations,
    draws = custom_draws
  )
  simulation <- simulate(
    model,
    beta = setNames(numeric(0), character(0)),
    control = biogeme_control(
      model_name = if (isTRUE(explicit)) {
        "r_montecarlo_b03_explicit"
      } else {
        "r_montecarlo_b03"
      },
      number_of_draws = as.integer(number_of_draws),
      seed = as.integer(seed),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
  as.data.frame(simulation, check.names = FALSE)
}

test_that("Monte Carlo Group 2 example files are syntactically valid", {
  files <- c(
    rbiogeme_example_path( "montecarlo", "plot_b03antithetic.R"),
    rbiogeme_example_path( "montecarlo", "plot_b03antithetic_explicit.R")
  )
  expect_true(all(file.exists(files)))
  for (file in files) parse(file)
})

test_that("b03 antithetic draws 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"
  )

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

  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b03(number_of_draws, seed)
  native_values <- as.data.frame(
    native_montecarlo_b03(number_of_draws, seed),
    check.names = FALSE
  )
  expect_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
})

test_that("b03 explicit antithetic pairs 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"
  )

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

  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b03(number_of_draws, seed, explicit = TRUE)
  native_values <- as.data.frame(
    native_montecarlo_b03(number_of_draws, seed, explicit = TRUE),
    check.names = FALSE
  )
  expect_identical(names(r_values), names(native_values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native_values)),
    tolerance = 1e-12
  )
})