inst/examples/profiling/plot_cnl_profiling.R

In rbiogeme: Interface to 'Biogeme'

#!/usr/bin/env Rscript

# Profiling the cross-nested logit expression.
#
# This is the R counterpart of Biogeme's plot_cnl_profiling.py. The model
# specification is deliberately written out here so the example is
# self-contained. `profile_jax()` sends the compiled expression graph to the
# native Biogeme engine; R does not evaluate a callback during profiling.

library(rbiogeme)

# prepare_swissmetro_example() is defined in ../swissmetro/example_utils.R.
# It handles command-line data/output setup only. The database transformations
# and the complete CNL specification remain visible in this script.
script_path <- commandArgs(trailingOnly = FALSE)
script_path <- sub("^--file=", "", script_path[startsWith(script_path, "--file=")][[1L]])
source(file.path(dirname(normalizePath(script_path)), "..", "swissmetro", "example_utils.R"))

prepared <- prepare_swissmetro_example(
  commandArgs(trailingOnly = TRUE),
  default_model = "cnl_profiling"
)

# swissmetro_data() performs the same native database operations as the
# profiling example: it keeps PURPOSE 1 and 3, removes CHOICE == 0, creates
# the cost and availability variables, and scales time/cost by 100.
database <- swissmetro_data(prepared$data, filter_purpose = TRUE)

# Parameters retain the exact names, starts, bounds, and fixed status of the
# Python example. The fixed asc_sm parameter is written with fixed = TRUE.
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_time_swissmetro <- biogeme_beta("b_time_swissmetro", start = 0)
b_time_train <- biogeme_beta("b_time_train", start = 0)
b_time_car <- biogeme_beta("b_time_car", start = 0)
b_cost <- biogeme_beta("b_cost", start = 0)
b_headway_swissmetro <- biogeme_beta("b_headway_swissmetro", start = 0)
b_headway_train <- biogeme_beta("b_headway_train", start = 0)
ga_train <- biogeme_beta("ga_train", start = 0)
ga_swissmetro <- biogeme_beta("ga_swissmetro", start = 0)

# Nest parameters and allocation parameters use the native CNL bounds.
existing_nest_parameter <- biogeme_beta(
  "existing_nest_parameter", start = 1, lower = 1, upper = 5
)
public_nest_parameter <- biogeme_beta(
  "public_nest_parameter", start = 1, lower = 1, upper = 5
)
alpha_existing <- biogeme_beta(
  "alpha_existing", start = 0.5, lower = 0, upper = 1
)
alpha_public <- 1 - alpha_existing

# Utilities are ordinary neutral Biogeme expressions. Multiplication by a
# variable means symbolic multiplication, not local R evaluation.
v_train <- asc_train +
  b_time_train * variable("TRAIN_TT_SCALED") +
  b_cost * variable("TRAIN_COST_SCALED") +
  b_headway_train * variable("TRAIN_HE") +
  ga_train * variable("GA")
v_swissmetro <- asc_sm +
  b_time_swissmetro * variable("SM_TT_SCALED") +
  b_cost * variable("SM_COST_SCALED") +
  b_headway_swissmetro * variable("SM_HE") +
  ga_swissmetro * variable("GA")
v_car <- asc_car +
  b_time_car * variable("CAR_TT_SCALED") +
  b_cost * variable("CAR_CO_SCALED")

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")
)

# Each CNL allocation is a native expression. Numeric zero/one allocations
# are retained in the compiled tree, as in OneNestForCrossNestedLogit.
existing_nest <- cross_nested_nest(
  existing_nest_parameter,
  list(`1` = alpha_existing, `2` = 0, `3` = 1),
  name = "existing"
)
public_nest <- cross_nested_nest(
  public_nest_parameter,
  list(`1` = alpha_public, `2` = 1, `3` = 0),
  name = "public"
)
nests <- cross_nested_nests(
  choice_set = c(1, 2, 3),
  nests = list(existing_nest, public_nest)
)

# cross_nested_logit_model() is a declarative wrapper. The likelihood itself
# is constructed and evaluated by native Biogeme after compilation.
control <- biogeme_control(
    output_directory = prepared$output,
  generate_html = FALSE,
  generate_yaml = FALSE,
  save_iterations = FALSE
)
model <- cross_nested_logit_model(
  database = database,
  choice = "CHOICE",
  utilities = utilities,
  availability = availability,
  nests = nests,
  control = control
)

# Use the same starting values as the native profiling script. The fixed
# asc_sm entry is included because the native evaluator accepts the complete
# beta mapping, including fixed parameters.
the_betas <- c(
  asc_car = 0,
  asc_train = 0,
  asc_sm = 0,
  b_time_swissmetro = 0,
  b_time_train = 0,
  b_time_car = 0,
  b_cost = 0,
  b_headway_swissmetro = 0,
  b_headway_train = 0,
  ga_train = 0,
  ga_swissmetro = 0,
  existing_nest_parameter = 1,
  public_nest_parameter = 1,
  alpha_existing = 0.5
)

# The four cases are the native example's function-only, gradient, Hessian,
# and BHHH profiles. Each case is evaluated twice by profile_jax().
cases <- list(
  list(label = "Function only", gradient = FALSE, hessian = FALSE, bhhh = FALSE),
  list(label = "Function + gradient", gradient = TRUE, hessian = FALSE, bhhh = FALSE),
  list(label = "Function + gradient + Hessian", gradient = TRUE, hessian = TRUE, bhhh = FALSE),
  list(label = "Function + gradient + BHHH", gradient = TRUE, hessian = FALSE, bhhh = TRUE)
)

result <- profile_jax(
  model = model,
  beta_values = the_betas,
  cases = cases,
  model_name = "cnl_profiling",
  control = control,
  numerically_safe = FALSE
)
print(result)

invisible(result)