tests/testthat/test-swissmetro-b13-panel-simul.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_b13_panel_simul <- function(data, number_of_draws = 128L, seed = 1223L) {
  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)
  single_formula <- reticulate::import(
    "biogeme.jax_calculator.single_formula",
    convert = FALSE
  )
  second_derivatives <- reticulate::import(
    "biogeme.second_derivatives",
    convert = FALSE
  )
  numpy <- reticulate::import("numpy", convert = FALSE)
  bridge <- rbiogeme:::biogeme_bridge()

  database <- database_module$Database(
    "swissmetro_native_b13_panel_simul",
    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)
  database$panel("ID")

  beta <- expressions$Beta
  b_cost <- beta("b_cost", 0, NULL, 0, 0)
  b_time <- beta("b_time", 0, NULL, 0, 0)
  b_time_s <- beta("b_time_s", 1, 1.0e-5, NULL, 0)
  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_car_s <- beta("asc_car_s", 1, 1.0e-5, NULL, 0)
  asc_train <- beta("asc_train", 0, NULL, NULL, 0)
  asc_train_s <- beta("asc_train_s", 1, 1.0e-5, NULL, 0)
  asc_sm <- beta("asc_sm", 0, NULL, NULL, 0)
  asc_sm_s <- beta("asc_sm_s", 1, 1.0e-5, NULL, 0)
  b_time_rnd <- b_time + b_time_s * expressions$Draws("b_time_rnd", "NORMAL_ANTI")
  asc_car_rnd <- asc_car + asc_car_s * expressions$Draws("asc_car_rnd", "NORMAL_ANTI")
  asc_train_rnd <- asc_train + asc_train_s * expressions$Draws("asc_train_rnd", "NORMAL_ANTI")
  asc_sm_rnd <- asc_sm + asc_sm_s * expressions$Draws("asc_sm_rnd", "NORMAL_ANTI")

  utilities <- reticulate::dict(
    `1` = asc_train_rnd + b_time_rnd * train_tt_scaled + b_cost * train_cost_scaled,
    `2` = asc_sm_rnd + b_time_rnd * sm_tt_scaled + b_cost * sm_cost_scaled,
    `3` = asc_car_rnd + b_time_rnd * car_tt_scaled + b_cost * car_co_scaled
  )
  availability <- reticulate::dict(
    `1` = train_av_sp,
    `2` = variable("SM_AV"),
    `3` = car_av_sp
  )
  kernel <- models$logit(utilities, availability, choice)
  trajectory <- expressions$PanelLikelihoodTrajectory(kernel)
  log_probability <- expressions$log(expressions$MonteCarlo(trajectory))

  estimator <- biogeme_module$BIOGEME(
    database,
    log_probability,
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  estimator$model_name <- "b12_panel"
  estimation_results <- estimator$estimate()
  betas <- estimation_results$get_beta_values()

  numerator <- expressions$MonteCarlo(b_time_rnd * trajectory)
  denominator <- expressions$MonteCarlo(trajectory)
  numpy_state <- numpy$random$get_state()
  numpy$random$seed(as.integer(seed))
  direct_loglike <- tryCatch(
    single_formula$calculate_single_formula_from_expression(
      expression = log_probability,
      database = database,
      number_of_draws = as.integer(number_of_draws),
      the_betas = betas,
      second_derivatives_mode = second_derivatives$SecondDerivativesMode$NEVER,
      numerically_safe = FALSE,
      use_jit = TRUE
    ),
    finally = numpy$random$set_state(numpy_state)
  )

  simulator <- biogeme_module$BIOGEME(
    database,
    reticulate::dict(Numerator = numerator, Denominator = denominator),
    number_of_draws = as.integer(number_of_draws),
    seed = as.integer(seed),
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  simulator$model_name <- "b13_panel_simul_native"
  simulator$use_flatten_database <- TRUE
  simulated <- reticulate::py_to_r(simulator$simulate(the_beta_values = betas))
  simulated[["Individual-level parameters"]] <-
    simulated[["Numerator"]] / simulated[["Denominator"]]

  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(estimation_results)),
    direct_loglike = as.numeric(reticulate::py_to_r(direct_loglike)),
    values = simulated,
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe)),
    number_of_individuals = length(unique(reticulate::py_to_r(database$dataframe)$ID))
  )
}

test_that("b13 panel 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, panel = TRUE)
  make_components <- function(b12_bounds) {
    if (b12_bounds) {
      b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
      b_time <- biogeme_beta("b_time", start = 0, upper = 0)
      b_time_s <- biogeme_beta("b_time_s", start = 1, lower = 1.0e-5)
      asc_car_s <- biogeme_beta("asc_car_s", start = 1, lower = 1.0e-5)
      asc_train_s <- biogeme_beta("asc_train_s", start = 1, lower = 1.0e-5)
      asc_sm_s <- biogeme_beta("asc_sm_s", start = 1, lower = 1.0e-5)
    } else {
      b_cost <- biogeme_beta("b_cost", start = 0)
      b_time <- biogeme_beta("b_time", start = 0)
      b_time_s <- biogeme_beta("b_time_s", start = 1)
      asc_car_s <- biogeme_beta("asc_car_s", start = 1)
      asc_train_s <- biogeme_beta("asc_train_s", start = 1)
      asc_sm_s <- biogeme_beta("asc_sm_s", start = 1)
    }
    asc_car <- biogeme_beta("asc_car", start = 0)
    asc_train <- biogeme_beta("asc_train", start = 0)
    asc_sm <- biogeme_beta("asc_sm", start = 0)
    b_time_rnd <- b_time + b_time_s * draw("b_time_rnd", "NORMAL_ANTI")
    asc_car_rnd <- asc_car + asc_car_s * draw("asc_car_rnd", "NORMAL_ANTI")
    asc_train_rnd <- asc_train + asc_train_s * draw("asc_train_rnd", "NORMAL_ANTI")
    asc_sm_rnd <- asc_sm + asc_sm_s * draw("asc_sm_rnd", "NORMAL_ANTI")
    utilities <- list(
      `1` = asc_train_rnd + b_time_rnd * variable("TRAIN_TT_SCALED") + b_cost * variable("TRAIN_COST_SCALED"),
      `2` = asc_sm_rnd + b_time_rnd * variable("SM_TT_SCALED") + b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car_rnd + b_time_rnd * 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")
    )
    kernel <- logit_probability(utilities, availability, variable("CHOICE"))
    trajectory <- panel_likelihood_trajectory(kernel)
    log_probability <- log(monte_carlo(trajectory))
    numerator <- monte_carlo(b_time_rnd * trajectory)
    denominator <- monte_carlo(trajectory)
    list(
      model = biogeme_model(
        database = database,
        formula = log_probability,
        simulations = if (b12_bounds) NULL else list(
          Numerator = numerator,
          Denominator = denominator
        ),
        draws = list(
          biogeme_draws("b_time_rnd", "NORMAL_ANTI", 128L, 1223L),
          biogeme_draws("asc_car_rnd", "NORMAL_ANTI", 128L, 1223L),
          biogeme_draws("asc_train_rnd", "NORMAL_ANTI", 128L, 1223L),
          biogeme_draws("asc_sm_rnd", "NORMAL_ANTI", 128L, 1223L)
        ),
        control = biogeme_control(
          model_name = if (b12_bounds) "b12_panel" else "b13_panel_simul",
          number_of_draws = 128L,
          seed = 1223L,
          generate_html = FALSE,
          generate_yaml = FALSE,
          save_iterations = FALSE
        )
      ),
      log_probability = log_probability,
      numerator = numerator,
      denominator = denominator,
      draws = list(
        biogeme_draws("b_time_rnd", "NORMAL_ANTI", 128L, 1223L),
        biogeme_draws("asc_car_rnd", "NORMAL_ANTI", 128L, 1223L),
        biogeme_draws("asc_train_rnd", "NORMAL_ANTI", 128L, 1223L),
        biogeme_draws("asc_sm_rnd", "NORMAL_ANTI", 128L, 1223L)
      )
    )
  }

  estimation <- make_components(TRUE)
  simulation <- make_components(FALSE)
  temporary_directory <- tempfile("rbiogeme-b13-")
  dir.create(temporary_directory, recursive = TRUE)
  original_directory <- getwd()
  setwd(temporary_directory)
  on.exit(setwd(original_directory), add = TRUE)

  r_fit <- estimate(estimation$model, model_name = "b12_panel", control = estimation$model$control)
  r_loglike <- simulate_single_formula(
    simulation$model,
    simulation$log_probability,
    r_fit,
    number_of_draws = 128L,
    seed = 1223L
  )
  r_values <- as.data.frame(simulate(
    simulation$model,
    beta = r_fit,
    control = biogeme_control(
      model_name = "b13_panel_simul",
      number_of_draws = 128L,
      seed = 1223L,
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  ), check.names = FALSE)
  r_values[["Individual-level parameters"]] <-
    r_values[["Numerator"]] / r_values[["Denominator"]]

  native <- native_swissmetro_b13_panel_simul(data, 128L, 1223L)
  expect_equal(nobs(r_fit), native$number_of_individuals)
  expect_identical(r_fit$beta_names, native$results$beta_names)
  expect_equal(unname(coef(r_fit)), native$results$beta_values, tolerance = 1e-7)
  expect_equal(as.numeric(logLik(r_fit)), native$results$final_log_likelihood, tolerance = 1e-7)
  expect_equal(r_loglike, native$direct_loglike, tolerance = 1e-7)
  expect_identical(names(r_values), names(native$values))
  expect_equal(
    unname(as.matrix(r_values)),
    unname(as.matrix(native$values)),
    tolerance = 1e-7,
    ignore_attr = TRUE
  )
})