tests/testthat/helper-swissmetro-mixtures.R

In rbiogeme: Interface to 'Biogeme'

native_swissmetro_mixture <- function(data, kind, number_of_draws = 256L) {
  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)
  bridge <- rbiogeme:::biogeme_bridge()

  database <- database_module$Database(
    paste0("swissmetro_native_", kind),
    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)
  if (identical(kind, "b26")) database$panel("ID")

  beta <- expressions$Beta
  asc_car <- beta("asc_car", 0, NULL, NULL, 0)
  asc_train <- beta("asc_train", 0, NULL, NULL, 0)
  asc_sm <- beta("asc_sm", 0, NULL, NULL, 1)
  b_cost <- beta("b_cost", 0, NULL, NULL, 0)
  b_time <- beta("b_time", 0, NULL, NULL, 0)
  b_time_s <- beta("b_time_s", 1, NULL, NULL, 0)

  random_number_generators <- NULL
  if (identical(kind, "b24")) {
    b_time_draws <- expressions$Draws("b_time_rnd", "NORMAL_HALTON5")
  } else {
    b_time_draws <- expressions$Draws("b_time_rnd", "TRIANGULAR")
    if (identical(kind, "b25")) {
      generator <- reticulate::py_eval(
        "lambda sample_size, number_of_draws: __import__('numpy').random.triangular(-1, 0, 1, (sample_size, number_of_draws))",
        convert = FALSE
      )
      tuple <- reticulate::import("biogeme.draws", convert = FALSE)$RandomNumberGeneratorTuple
      random_number_generators <- reticulate::dict(
        TRIANGULAR = tuple(generator, "Draws from a triangular distribution")
      )
    }
  }

  if (identical(kind, "b26")) {
    asc_car_s <- beta("asc_car_s", 1, NULL, NULL, 0)
    asc_train_s <- beta("asc_train_s", 1, NULL, NULL, 0)
    asc_sm_s <- beta("asc_sm_s", 1, NULL, NULL, 0)
    asc_car_draws <- expressions$Draws("asc_car_rnd", "TRIANGULAR")
    asc_train_draws <- expressions$Draws("asc_train_rnd", "TRIANGULAR")
    asc_sm_draws <- expressions$Draws("asc_sm_rnd", "TRIANGULAR")
    asc_car <- asc_car + asc_car_s * asc_car_draws
    asc_train <- asc_train + asc_train_s * asc_train_draws
    asc_sm <- asc_sm + asc_sm_s * asc_sm_draws
    random_number_generators <- reticulate::dict(
      TRIANGULAR = {
        generator <- reticulate::py_eval(
          "lambda sample_size, number_of_draws: __import__('numpy').random.triangular(-1, 0, 1, (sample_size, number_of_draws))",
          convert = FALSE
        )
        tuple <- reticulate::import("biogeme.draws", convert = FALSE)$RandomNumberGeneratorTuple
        tuple(generator, "Draws from a triangular distribution")
      }
    )
  }

  v_train <- asc_train + (b_time + b_time_s * b_time_draws) * train_tt_scaled +
    b_cost * train_cost_scaled
  v_swissmetro <- asc_sm + (b_time + b_time_s * b_time_draws) * sm_tt_scaled + b_cost * sm_cost_scaled
  v_car <- asc_car + (b_time + b_time_s * b_time_draws) * car_tt_scaled + b_cost * car_co_scaled
  conditional_probability <- models$logit(
    reticulate::dict(`1` = v_train, `2` = v_swissmetro, `3` = v_car),
    reticulate::dict(`1` = train_av_sp, `2` = variable("SM_AV"), `3` = car_av_sp),
    choice
  )
  if (identical(kind, "b26")) {
    conditional_probability <- expressions$PanelLikelihoodTrajectory(conditional_probability)
  }
  log_probability <- expressions$log(expressions$MonteCarlo(conditional_probability))
  arguments <- list(
    database = database,
    formulas = log_probability,
    number_of_draws = as.integer(number_of_draws),
    seed = 1223L,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  if (identical(kind, "b24")) arguments$analytical_hessian_mode <- "automatic"
  if (identical(kind, "b25")) {
    arguments$analytical_hessian_mode <- "automatic"
    arguments$random_number_generators <- random_number_generators
  }
  if (identical(kind, "b26")) {
    arguments$calculating_second_derivatives <- "never"
    arguments$random_number_generators <- random_number_generators
  }
  biogeme <- do.call(biogeme_module$BIOGEME, arguments)
  biogeme$model_name <- paste0(kind, "_native")
  results <- biogeme$estimate()
  list(
    results = reticulate::py_to_r(bridge$extract_estimation_results(results)),
    number_of_rows = nrow(reticulate::py_to_r(database$dataframe))
  )
}

r_swissmetro_mixture <- function(data, kind, number_of_draws = 256L) {
  database <- swissmetro_data(data, panel = identical(kind, "b26"))
  asc_car <- biogeme_beta("asc_car", start = 0)
  asc_train <- biogeme_beta("asc_train", start = 0)
  asc_sm <- biogeme_beta("asc_sm", start = 0, fixed = TRUE)
  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)
  if (identical(kind, "b26")) {
    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)
  }
  generator <- if (kind %in% c("b24", "b27")) NULL else "TRIANGULAR"
  draw_type <- if (identical(kind, "b24")) {
    "NORMAL_HALTON5"
  } else if (identical(kind, "b27")) {
    "NORMAL"
  } else {
    "TRIANGULAR"
  }
  b_time_draws <- biogeme_draws(
    "b_time_rnd", draw_type, number_of_draws = number_of_draws,
    seed = 1223, generator = generator
  )
  b_time_rnd <- b_time + b_time_s * b_time_draws
  asc_car_rnd <- asc_car
  asc_train_rnd <- asc_train
  asc_sm_rnd <- asc_sm
  draws <- list(b_time_draws)
  if (identical(kind, "b26")) {
    asc_car_draws <- biogeme_draws("asc_car_rnd", "TRIANGULAR", number_of_draws = number_of_draws, seed = 1223, generator = generator)
    asc_train_draws <- biogeme_draws("asc_train_rnd", "TRIANGULAR", number_of_draws = number_of_draws, seed = 1223, generator = generator)
    asc_sm_draws <- biogeme_draws("asc_sm_rnd", "TRIANGULAR", number_of_draws = number_of_draws, seed = 1223, generator = generator)
    asc_car_rnd <- asc_car + asc_car_s * asc_car_draws
    asc_train_rnd <- asc_train + asc_train_s * asc_train_draws
    asc_sm_rnd <- asc_sm + asc_sm_s * asc_sm_draws
    draws <- c(draws, list(asc_car_draws, asc_train_draws, asc_sm_draws))
  }
  v_train <- asc_train_rnd + b_time_rnd * variable("TRAIN_TT_SCALED") + b_cost * variable("TRAIN_COST_SCALED")
  v_swissmetro <- asc_sm_rnd + b_time_rnd * variable("SM_TT_SCALED") + b_cost * variable("SM_COST_SCALED")
  v_car <- asc_car_rnd + b_time_rnd * variable("CAR_TT_SCALED") + b_cost * variable("CAR_CO_SCALED")
  conditional_probability <- logit_probability(
    utilities = list(`1` = v_train, `2` = v_swissmetro, `3` = v_car),
    availability = list(`1` = variable("TRAIN_AV_SP"), `2` = variable("SM_AV"), `3` = variable("CAR_AV_SP")),
    alternative = variable("CHOICE")
  )
  if (identical(kind, "b26")) conditional_probability <- panel_likelihood_trajectory(conditional_probability)
  control <- biogeme_control(
    model_name = paste0(kind, "_r"),
    number_of_draws = number_of_draws,
    seed = 1223,
    generate_html = FALSE,
    generate_yaml = FALSE,
    save_iterations = FALSE
  )
  if (kind %in% c("b24", "b25")) control$analytical_hessian_mode <- "automatic"
  if (identical(kind, "b26")) control$calculating_second_derivatives <- "never"
  model <- biogeme_model(
    database = database,
    formula = log(monte_carlo(conditional_probability)),
    draws = draws,
    control = control
  )
  list(model = model, database = database)
}