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

assisted_group2_files <- function() {
  rbiogeme_example_path(
    "assisted",
    c("plot_b01model.R", "plot_b02nonlinear.R", "plot_b03alt_spec.R")
  )
}

test_that("assisted group 2 examples are self-contained", {
  files <- assisted_group2_files()
  expect_true(all(file.exists(files)))
  for (file in files) {
    expect_silent(parse(file = file))
    contents <- paste(readLines(file, warn = FALSE), collapse = "\n")
    expect_match(contents, "biogeme_beta")
    expect_match(contents, "estimate_catalog")
    expect_match(contents, "prepare_swissmetro_example")
    expect_match(contents, "filter_purpose = FALSE")
  }
  expect_match(paste(readLines(files[[1L]], warn = FALSE), collapse = "\n"), "b01model")
  expect_match(paste(readLines(files[[2L]], warn = FALSE), collapse = "\n"), "b02nonlinear")
  expect_match(paste(readLines(files[[3L]], warn = FALSE), collapse = "\n"), "b01alt_spec")
})

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

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

native_assisted_group2_results <- function(data, specification) {
  expressions <- reticulate::import("biogeme.expressions", convert = FALSE)
  biogeme_module <- reticulate::import("biogeme.biogeme", convert = FALSE)
  catalog_module <- reticulate::import("biogeme.catalog", convert = FALSE)
  models <- reticulate::import("biogeme.models", convert = FALSE)
  nests_module <- reticulate::import("biogeme.nests", convert = FALSE)
  results_processing <- reticulate::import(
    "biogeme.results_processing",
    convert = FALSE
  )
  bridge <- rbiogeme:::biogeme_bridge()
  parts <- native_assisted_group2_database(
    data,
    paste0("swissmetro_native_assisted_", specification)
  )
  variable <- parts$variable
  beta <- expressions$Beta
  utilities <- NULL
  model_name <- switch(
    specification,
    b01model = "b01model",
    b02nonlinear = "b02nonlinear",
    b03alt_spec = "b01alt_spec"
  )

  if (identical(specification, "b01model")) {
    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_cost <- beta("b_cost", 0, NULL, NULL, 0)
    utilities <- reticulate::dict(
      `1` = asc_train + b_time * parts$train_tt_scaled + b_cost * parts$train_cost_scaled,
      `2` = b_time * parts$sm_tt_scaled + b_cost * parts$sm_cost_scaled,
      `3` = asc_car + b_time * parts$car_tt_scaled + b_cost * parts$car_co_scaled
    )
    availability <- reticulate::dict(
      `1` = parts$train_av_sp,
      `2` = parts$sm_av,
      `3` = parts$car_av_sp
    )
    logit <- models$loglogit(utilities, availability, parts$choice)
    mu_existing <- beta("mu_existing", 1, 1, 10, 0)
    existing <- nests_module$OneNestForNestedLogit(
      nest_param = mu_existing,
      list_of_alternatives = reticulate::r_to_py(as.integer(c(1L, 3L))),
      name = "Existing"
    )
    nests_existing <- nests_module$NestsForNestedLogit(
      choice_set = reticulate::r_to_py(as.integer(c(1L, 2L, 3L))),
      tuple_of_nests = reticulate::tuple(existing)
    )
    nested_existing <- models$lognested(
      utilities,
      availability,
      nests_existing,
      parts$choice
    )
    mu_public <- beta("mu_public", 1, 1, 10, 0)
    public <- nests_module$OneNestForNestedLogit(
      nest_param = mu_public,
      list_of_alternatives = reticulate::r_to_py(as.integer(c(1L, 2L))),
      name = "Public"
    )
    nests_public <- nests_module$NestsForNestedLogit(
      choice_set = reticulate::r_to_py(as.integer(c(1L, 2L, 3L))),
      tuple_of_nests = reticulate::tuple(public)
    )
    nested_public <- models$lognested(
      utilities,
      availability,
      nests_public,
      parts$choice
    )
    expression_catalog <- reticulate::dict()
    reticulate::py_set_item(expression_catalog, "logit", logit)
    reticulate::py_set_item(expression_catalog, "nested existing", nested_existing)
    reticulate::py_set_item(expression_catalog, "nested public", nested_public)
    formula <- catalog_module$Catalog$from_dict(
      catalog_name = "model_catalog",
      dict_of_expressions = expression_catalog
    )
  } else if (identical(specification, "b02nonlinear")) {
    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, 0, 0)
    b_cost <- beta("b_cost", 0, NULL, 0, 0)
    lambda_travel_time <- beta("lambda_travel_time", 1, -10, 10, 0)
    square_tt_coef <- beta("square_tt_coef", 0, NULL, NULL, 0)
    cube_tt_coef <- beta("cube_tt_coef", 0, NULL, NULL, 0)
    power_series <- function(the_variable) {
      the_variable + square_tt_coef * the_variable^2 +
        cube_tt_coef * the_variable * the_variable^3
    }
    controller <- catalog_module$Controller(
      controller_name = "train_tt_catalog",
      specification_names = reticulate::r_to_py(c("linear", "boxcox", "power"))
    )
    train_catalog <- catalog_module$Catalog$from_dict(
      catalog_name = "train_tt_catalog",
      dict_of_expressions = reticulate::dict(
        linear = parts$train_tt_scaled,
        boxcox = models$boxcox(parts$train_tt_scaled, lambda_travel_time),
        power = power_series(parts$train_tt_scaled)
      ),
      controlled_by = controller
    )
    sm_catalog <- catalog_module$Catalog$from_dict(
      catalog_name = "sm_tt_catalog",
      dict_of_expressions = reticulate::dict(
        linear = parts$sm_tt_scaled,
        boxcox = models$boxcox(parts$sm_tt_scaled, lambda_travel_time),
        power = power_series(parts$sm_tt_scaled)
      ),
      controlled_by = controller
    )
    car_catalog <- catalog_module$Catalog$from_dict(
      catalog_name = "car_tt_catalog",
      dict_of_expressions = reticulate::dict(
        linear = parts$car_tt_scaled,
        boxcox = models$boxcox(parts$car_tt_scaled, lambda_travel_time),
        power = power_series(parts$car_tt_scaled)
      ),
      controlled_by = controller
    )
    utilities <- reticulate::dict(
      `1` = asc_train + b_time * train_catalog + b_cost * parts$train_cost_scaled,
      `2` = b_time * sm_catalog + b_cost * parts$sm_cost_scaled,
      `3` = asc_car + b_time * car_catalog + b_cost * parts$car_co_scaled
    )
    availability <- reticulate::dict(
      `1` = parts$train_av_sp,
      `2` = parts$sm_av,
      `3` = parts$car_av_sp
    )
    formula <- models$loglogit(utilities, availability, parts$choice)
  } else {
    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_cost <- beta("b_cost", 0, NULL, NULL, 0)
    time_catalogs <- reticulate::py_get_item(catalog_module$generic_alt_specific_catalogs(
      generic_name = "b_time",
      beta_parameters = reticulate::r_to_py(list(b_time)),
      alternatives = reticulate::r_to_py(c("train", "swissmetro", "car"))
    ), 0L)
    cost_catalogs <- reticulate::py_get_item(catalog_module$generic_alt_specific_catalogs(
      generic_name = "b_cost",
      beta_parameters = reticulate::r_to_py(list(b_cost)),
      alternatives = reticulate::r_to_py(c("train", "swissmetro", "car"))
    ), 0L)
    utilities <- reticulate::dict(
      `1` = asc_train + time_catalogs$train * parts$train_tt_scaled +
        cost_catalogs$train * parts$train_cost_scaled,
      `2` = time_catalogs$swissmetro * parts$sm_tt_scaled +
        cost_catalogs$swissmetro * parts$sm_cost_scaled,
      `3` = asc_car + time_catalogs$car * parts$car_tt_scaled +
        cost_catalogs$car * parts$car_co_scaled
    )
    availability <- reticulate::dict(
      `1` = parts$train_av_sp,
      `2` = parts$sm_av,
      `3` = parts$car_av_sp
    )
    formula <- models$loglogit(utilities, availability, parts$choice)
  }

  estimator <- biogeme_module$BIOGEME(
    parts$database,
    formula,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  estimator$model_name <- model_name
  native_results <- estimator$estimate_catalog()
  keys <- vapply(
    reticulate::iterate(native_results$keys()),
    as.character,
    character(1)
  )
  serialized <- lapply(keys, function(key) {
    result <- reticulate::py_get_item(native_results, key)
    reticulate::py_to_r(bridge$extract_estimation_results(result))
  })
  names(serialized) <- keys
  non_dominated <- results_processing$pareto_optimal(native_results)
  list(
    results = serialized,
    non_dominated = vapply(
      reticulate::iterate(non_dominated$keys()),
      as.character,
      character(1)
    ),
    number_of_rows = nrow(reticulate::py_to_r(parts$database$dataframe))
  )
}

r_assisted_group2_model <- function(database, specification) {
  if (identical(specification, "b01model")) {
    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_cost <- biogeme_beta("b_cost", start = 0)
    utilities <- list(
      `1` = asc_train + b_time * variable("TRAIN_TT_SCALED") +
        b_cost * variable("TRAIN_COST_SCALED"),
      `2` = b_time * variable("SM_TT_SCALED") + b_cost * variable("SM_COST_SCALED"),
      `3` = asc_car + b_time * 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")
    )
    choice <- variable("CHOICE")
    mu_existing <- biogeme_beta("mu_existing", start = 1, lower = 1, upper = 10)
    mu_public <- biogeme_beta("mu_public", start = 1, lower = 1, upper = 10)
    log_probability <- catalog(
      "model_catalog",
      list(
        logit = logit_log_probability(utilities, availability, choice),
        `nested existing` = nested_log_probability(
          utilities, availability,
          nested_nests(c(1, 2, 3), list(nested_nest(mu_existing, c(1, 3), "Existing"))),
          choice
        ),
        `nested public` = nested_log_probability(
          utilities, availability,
          nested_nests(c(1, 2, 3), list(nested_nest(mu_public, c(1, 2), "Public"))),
          choice
        )
      )
    )
  } else if (identical(specification, "b02nonlinear")) {
    asc_car <- biogeme_beta("asc_car", start = 0)
    asc_train <- biogeme_beta("asc_train", start = 0)
    b_time <- biogeme_beta("b_time", start = 0, upper = 0)
    b_cost <- biogeme_beta("b_cost", start = 0, upper = 0)
    lambda_travel_time <- biogeme_beta(
      "lambda_travel_time", start = 1, lower = -10, upper = 10
    )
    square_tt_coef <- biogeme_beta("square_tt_coef", start = 0)
    cube_tt_coef <- biogeme_beta("cube_tt_coef", start = 0)
    power_series <- function(the_variable) {
      the_variable + square_tt_coef * the_variable^2 +
        cube_tt_coef * the_variable * the_variable^3
    }
    controller <- catalog_controller(
      "train_tt_catalog", c("linear", "boxcox", "power")
    )
    train_catalog <- catalog(
      "train_tt_catalog",
      list(
        linear = variable("TRAIN_TT_SCALED"),
        boxcox = boxcox(variable("TRAIN_TT_SCALED"), lambda_travel_time),
        power = power_series(variable("TRAIN_TT_SCALED"))
      ),
      controller
    )
    sm_catalog <- catalog(
      "sm_tt_catalog",
      list(
        linear = variable("SM_TT_SCALED"),
        boxcox = boxcox(variable("SM_TT_SCALED"), lambda_travel_time),
        power = power_series(variable("SM_TT_SCALED"))
      ),
      controller
    )
    car_catalog <- catalog(
      "car_tt_catalog",
      list(
        linear = variable("CAR_TT_SCALED"),
        boxcox = boxcox(variable("CAR_TT_SCALED"), lambda_travel_time),
        power = power_series(variable("CAR_TT_SCALED"))
      ),
      controller
    )
    log_probability <- logit_log_probability(
      list(
        `1` = asc_train + b_time * train_catalog + b_cost * variable("TRAIN_COST_SCALED"),
        `2` = b_time * sm_catalog + b_cost * variable("SM_COST_SCALED"),
        `3` = asc_car + b_time * car_catalog + b_cost * variable("CAR_CO_SCALED")
      ),
      list(`1` = variable("TRAIN_AV_SP"), `2` = variable("SM_AV"), `3` = variable("CAR_AV_SP")),
      variable("CHOICE")
    )
  } else {
    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_cost <- biogeme_beta("b_cost", start = 0)
    time_catalogs <- generic_alt_specific_catalogs(
      "b_time", list(b_time), c("train", "swissmetro", "car")
    )[[1L]]
    cost_catalogs <- generic_alt_specific_catalogs(
      "b_cost", list(b_cost), c("train", "swissmetro", "car")
    )[[1L]]
    log_probability <- logit_log_probability(
      list(
        `1` = asc_train + time_catalogs$train * variable("TRAIN_TT_SCALED") +
          cost_catalogs$train * variable("TRAIN_COST_SCALED"),
        `2` = time_catalogs$swissmetro * variable("SM_TT_SCALED") +
          cost_catalogs$swissmetro * variable("SM_COST_SCALED"),
        `3` = asc_car + time_catalogs$car * variable("CAR_TT_SCALED") +
          cost_catalogs$car * variable("CAR_CO_SCALED")
      ),
      list(`1` = variable("TRAIN_AV_SP"), `2` = variable("SM_AV"), `3` = variable("CAR_AV_SP")),
      variable("CHOICE")
    )
  }
  biogeme_model(
    database,
    formula = log_probability,
    control = biogeme_control(
      model_name = switch(
        specification,
        b01model = "b01model",
        b02nonlinear = "b02nonlinear",
        b03alt_spec = "b01alt_spec"
      ),
      generate_html = FALSE,
      generate_yaml = FALSE,
      save_iterations = FALSE
    )
  )
}

test_that("assisted group 2 estimates match native Biogeme", {
  skip_if_not(
    identical(Sys.getenv("RBIOGEME_RUN_INTEGRATION"), "1"),
    "Set RBIOGEME_RUN_INTEGRATION=1 to run full assisted 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)

  for (specification in c("b01model", "b02nonlinear", "b03alt_spec")) {
    temporary_directory <- tempfile(paste0("rbiogeme-assisted-", specification, "-"))
    r_directory <- file.path(temporary_directory, "r")
    native_directory <- file.path(temporary_directory, "native")
    dir.create(r_directory, recursive = TRUE)
    dir.create(native_directory, recursive = TRUE)
    database <- swissmetro_data(
      data,
      name = paste0("swissmetro_assisted_", specification),
      filter_purpose = FALSE
    )
    original_directory <- getwd()
    on.exit(setwd(original_directory), add = TRUE)
    setwd(r_directory)
    r_fit <- estimate_catalog(
      r_assisted_group2_model(database, specification),
      model_name = switch(
        specification,
        b01model = "b01model",
        b02nonlinear = "b02nonlinear",
        b03alt_spec = "b01alt_spec"
      ),
      control = biogeme_control(
        generate_html = FALSE,
        generate_yaml = FALSE,
        save_iterations = FALSE
      ),
      force = TRUE
    )
    setwd(native_directory)
    native <- native_assisted_group2_results(data, specification)

    expect_equal(length(r_fit$results), length(native$results), info = specification)
    expect_equal(nobs(r_fit$results[[1L]]), native$number_of_rows, info = specification)
    expect_setequal(names(r_fit$results), names(native$results))
    for (configuration in names(native$results)) {
      r_result <- r_fit$results[[configuration]]
      native_result <- native$results[[configuration]]
      expect_identical(r_result$beta_names, native_result$beta_names, info = configuration)
      expect_equal(
        unname(coef(r_result)),
        native_result$beta_values,
        tolerance = 1e-7,
        info = configuration
      )
      expect_equal(
        as.numeric(logLik(r_result)),
        native_result$final_log_likelihood,
        tolerance = 1e-7,
        info = configuration
      )
    }
    expect_setequal(r_fit$non_dominated, native$non_dominated)
  }
})