tests/testthat/test-montecarlo-group3.R

In rbiogeme: Interface to 'Biogeme'

source(
  rbiogeme_example_path( "montecarlo", "swissmetro_one.R")
)

native_montecarlo_one_inputs <- function(data, database_name) {
  database_module <- reticulate::import("biogeme.database", convert = FALSE)
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)

  database <- database_module$Database(
    database_name,
    reticulate::r_to_py(data[1L, , drop = FALSE])
  )
  variable <- expressions$Variable
  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)
  list(
    database = database,
    choice = variable("CHOICE"),
    sm_av = variable("SM_AV"),
    car_av_sp = car_av_sp,
    train_av_sp = train_av_sp,
    sm_tt_scaled = sm_tt_scaled,
    sm_cost_scaled = sm_cost_scaled,
    train_tt_scaled = train_tt_scaled,
    train_cost_scaled = train_cost_scaled,
    car_tt_scaled = car_tt_scaled,
    car_co_scaled = car_co_scaled,
    expressions = expressions
  )
}

native_montecarlo_b04 <- function(data, number_of_draws, seed) {
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)
  inputs <- native_montecarlo_one_inputs(data, "native_montecarlo_b04")
  expressions <- inputs$expressions
  omega <- expressions$RandomVariable("omega")
  b_time_random <- -2.26 + 1.66 * omega
  utilities <- reticulate::dict(
    `1` = -0.402 + b_time_random * inputs$train_tt_scaled -
      1.29 * inputs$train_cost_scaled,
    `2` = b_time_random * inputs$sm_tt_scaled -
      1.29 * inputs$sm_cost_scaled,
    `3` = 0.137 + b_time_random * inputs$car_tt_scaled -
      1.29 * inputs$car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = inputs$train_av_sp,
    `2` = inputs$sm_av,
    `3` = inputs$car_av_sp
  )
  probability <- models$logit(utilities, availability, inputs$choice)
  simulations <- reticulate::dict(
    Numerical = expressions$IntegrateNormal(probability, "omega")
  )
  biogeme <- biogeme_module$BIOGEME(
    inputs$database,
    simulations,
    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())
  )
}

native_montecarlo_b05 <- function(data, number_of_draws, seed) {
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)
  inputs <- native_montecarlo_one_inputs(data, "native_montecarlo_b05")
  expressions <- inputs$expressions
  omega <- expressions$RandomVariable("omega")
  b_time_random <- -2.26 + 1.66 * omega
  b_time_random_normal <- -2.26 + 1.66 * expressions$Draws("b_normal", "NORMAL")
  b_time_random_anti <- -2.26 + 1.66 * expressions$Draws("b_anti", "NORMAL_ANTI")
  b_time_random_halton <- -2.26 + 1.66 * expressions$Draws("b_halton", "NORMAL_HALTON2")
  b_time_random_mlhs <- -2.26 + 1.66 * expressions$Draws("b_mlhs", "NORMAL_MLHS")
  b_time_random_antimlhs <- -2.26 + 1.66 *
    expressions$Draws("b_antimlhs", "NORMAL_MLHS_ANTI")
  conditional_logit <- function(random_coefficient) {
    utilities <- reticulate::dict(
      `1` = -0.402 + random_coefficient * inputs$train_tt_scaled -
        1.29 * inputs$train_cost_scaled,
      `2` = random_coefficient * inputs$sm_tt_scaled -
        1.29 * inputs$sm_cost_scaled,
      `3` = 0.137 + random_coefficient * inputs$car_tt_scaled -
        1.29 * inputs$car_co_scaled
    )
    availability <- reticulate::dict(
      `1` = inputs$train_av_sp,
      `2` = inputs$sm_av,
      `3` = inputs$car_av_sp
    )
    models$logit(utilities, availability, inputs$choice)
  }
  simulations <- reticulate::dict(
    Numerical = expressions$IntegrateNormal(
      conditional_logit(b_time_random),
      "omega"
    ),
    MonteCarlo = expressions$MonteCarlo(conditional_logit(b_time_random_normal)),
    Antithetic = expressions$MonteCarlo(conditional_logit(b_time_random_anti)),
    Halton = expressions$MonteCarlo(conditional_logit(b_time_random_halton)),
    MLHS = expressions$MonteCarlo(conditional_logit(b_time_random_mlhs)),
    `Antithetic MLHS` = expressions$MonteCarlo(
      conditional_logit(b_time_random_antimlhs)
    )
  )
  biogeme <- biogeme_module$BIOGEME(
    inputs$database,
    simulations,
    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_b04 <- function(data, number_of_draws, seed) {
  inputs <- prepare_swissmetro_one_database(data)
  omega <- random_variable("omega")
  b_time_random <- -2.26 + 1.66 * omega
  utilities <- list(
    `1` = -0.402 + b_time_random * inputs$train_tt_scaled -
      1.29 * inputs$train_cost_scaled,
    `2` = b_time_random * inputs$sm_tt_scaled -
      1.29 * inputs$sm_cost_scaled,
    `3` = 0.137 + b_time_random * inputs$car_tt_scaled -
      1.29 * inputs$car_co_scaled
  )
  availability <- list(
    `1` = inputs$train_av_sp,
    `2` = inputs$sm_av,
    `3` = inputs$car_av_sp
  )
  probability <- logit_probability(utilities, availability, inputs$choice)
  model <- biogeme_model(
    database = inputs$database,
    simulations = list(Numerical = integrate_normal(probability, "omega"))
  )
  simulation <- simulate(
    model,
    beta = setNames(numeric(0), character(0)),
    control = biogeme_control(
      model_name = "r_montecarlo_b04",
      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)
}

r_montecarlo_b05 <- function(data, number_of_draws, seed) {
  inputs <- prepare_swissmetro_one_database(data)
  omega <- random_variable("omega")
  b_time_random <- -2.26 + 1.66 * omega
  b_time_random_normal <- -2.26 + 1.66 * draw("b_normal", "NORMAL")
  b_time_random_anti <- -2.26 + 1.66 * draw("b_anti", "NORMAL_ANTI")
  b_time_random_halton <- -2.26 + 1.66 * draw("b_halton", "NORMAL_HALTON2")
  b_time_random_mlhs <- -2.26 + 1.66 * draw("b_mlhs", "NORMAL_MLHS")
  b_time_random_antimlhs <- -2.26 + 1.66 *
    draw("b_antimlhs", "NORMAL_MLHS_ANTI")
  conditional_logit <- function(random_coefficient) {
    logit_probability(
      utilities = list(
        `1` = -0.402 + random_coefficient * inputs$train_tt_scaled -
          1.29 * inputs$train_cost_scaled,
        `2` = random_coefficient * inputs$sm_tt_scaled -
          1.29 * inputs$sm_cost_scaled,
        `3` = 0.137 + random_coefficient * inputs$car_tt_scaled -
          1.29 * inputs$car_co_scaled
      ),
      availability = list(
        `1` = inputs$train_av_sp,
        `2` = inputs$sm_av,
        `3` = inputs$car_av_sp
      ),
      alternative = inputs$choice
    )
  }
  model <- biogeme_model(
    database = inputs$database,
    simulations = list(
      Numerical = integrate_normal(conditional_logit(b_time_random), "omega"),
      MonteCarlo = monte_carlo(conditional_logit(b_time_random_normal)),
      Antithetic = monte_carlo(conditional_logit(b_time_random_anti)),
      Halton = monte_carlo(conditional_logit(b_time_random_halton)),
      MLHS = monte_carlo(conditional_logit(b_time_random_mlhs)),
      `Antithetic MLHS` = monte_carlo(conditional_logit(b_time_random_antimlhs))
    )
  )
  simulation <- simulate(
    model,
    beta = setNames(numeric(0), character(0)),
    control = biogeme_control(
      model_name = "r_montecarlo_b05",
      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 3 files are syntactically valid", {
  files <- c(
    rbiogeme_example_path( "montecarlo", "swissmetro_one.R"),
    rbiogeme_example_path( "montecarlo", "plot_b04normal_mixture_numerical.R"),
    rbiogeme_example_path( "montecarlo", "plot_b05normal_mixture_monte_carlo.R")
  )
  expect_true(all(file.exists(files)))
  for (file in files) parse(file)
})

test_that("b04 numerical mixture simulation matches 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
  )
  temporary_directory <- tempfile("rbiogeme-montecarlo-group3-b04-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b04(data, number_of_draws, seed)
  native_values <- as.data.frame(
    native_montecarlo_b04(data, 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("b05 mixture integration methods 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
  )
  temporary_directory <- tempfile("rbiogeme-montecarlo-group3-b05-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  data <- read.delim(data_path, check.names = FALSE, stringsAsFactors = FALSE)
  number_of_draws <- 32L
  seed <- 1223L
  r_values <- r_montecarlo_b05(data, number_of_draws, seed)
  native_values <- as.data.frame(
    native_montecarlo_b05(data, 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
  )
})