inst/examples/indicators/optima.R

In rbiogeme: Interface to 'Biogeme'

# Exact Optima data preparation for the indicators examples.
#
# This helper contains data operations only. The filters are declared as
# Biogeme expressions and materialized by native Biogeme; the model equations
# remain visible in each runnable example script.

read_optima_database <- function(path, name = "optima") {
  if (!is.character(path) || length(path) != 1L || is.na(path) || !nzchar(path)) {
    stop("path must be one non-empty character string.", call. = FALSE)
  }
  if (!file.exists(path) || dir.exists(path)) {
    stop("The Optima data file does not exist: ", path, call. = FALSE)
  }
  data <- read.delim(path, check.names = FALSE, stringsAsFactors = FALSE)
  if (nrow(data) < 1L) {
    stop("The Optima data file has no observations.", call. = FALSE)
  }
  optima_database(data, name = name)
}

optima_database <- function(data, name = "optima") {
  database <- biogeme_database(name, data)
  choice <- variable("Choice")
  car_availability <- variable("CarAvail")

  # These are exactly the two filters in biogeme.data.optima.read_data():
  # remove unavailable choices and observations choosing an unavailable car.
  database <- biogeme_database_remove(database, choice == -1)
  database <- biogeme_database_remove(
    database,
    (choice == 1) & (car_availability == 3)
  )

  # Materialization applies both filters in native Biogeme and preserves the
  # original row identifiers. The normalized-weight scale is then computed
  # from precisely the retained rows, as in the Python helper.
  database <- biogeme_database_materialize(database)
  if (!biogeme_database_has_column(database, "Weight")) {
    stop("The Optima data must contain a Weight column.", call. = FALSE)
  }
  weight_sum <- sum(database$data$Weight)
  if (!is.finite(weight_sum) || weight_sum == 0) {
    stop("The retained Optima weights must have a finite non-zero sum.", call. = FALSE)
  }
  weight_scale <- nrow(database$data) / weight_sum

  define <- function(database, name, expression) {
    biogeme_database_define_variable(database, name, expression)
  }

  age <- variable("age")
  calculated_income <- variable("CalculatedIncome")
  household <- variable("NbHousehold")
  cars <- variable("NbCar")
  bikes <- variable("NbBicy")
  family <- variable("FamilSitu")
  child_residence <- variable("ResidChild")
  education <- variable("Education")
  socio_professional <- variable("SocioProfCat")
  time_pt <- variable("TimePT")
  time_car <- variable("TimeCar")
  marginal_cost_pt <- variable("MarginalCostPT")
  car_cost <- variable("CostCarCHF")
  distance <- variable("distance_km")
  trajectories <- variable("NbTrajects")

  # These derived columns mirror the complete native optima.py helper. They
  # remain symbolic and are evaluated by Biogeme only after compilation.
  database <- define(database, "normalized_weight", variable("Weight") * weight_scale)
  database <- define(database, "livesInUrbanArea", variable("UrbRur") == 2)
  database <- define(database, "owningHouse", variable("OwnHouse") == 1)
  database <- define(database, "ScaledIncome", calculated_income / 1000)
  database <- define(database, "age_65_more", age >= 65)
  database <- define(database, "age_30_less", age <= 30)
  database <- define(database, "age_category", 2 - (age <= 30) + (age >= 65))
  database <- define(
    database,
    "household_size",
    3 - 2 * (household == 1) - (household == 2)
  )
  database <- define(
    database,
    "income_category",
    1 + (calculated_income >= 3250) +
      (calculated_income >= 7000) + (calculated_income >= 15000)
  )
  database <- define(database, "moreThanOneCar", cars > 1)
  database <- define(database, "moreThanOneBike", bikes > 1)
  database <- define(database, "individualHouse", variable("HouseType") == 1)
  database <- define(database, "male", variable("Gender") == 1)
  database <- define(database, "single", (family == 1) + (family == 4) > 0)
  database <- define(
    database,
    "haveChildren",
    (family == 3) + (family == 4) > 0
  )
  database <- define(database, "haveGA", variable("GenAbST") == 1)
  database <- define(database, "highEducation", education >= 6)
  database <- define(database, "artisans", socio_professional == 5)
  database <- define(database, "employees", socio_professional == 6)
  database <- define(
    database,
    "occupation_status",
    4 - 3 * (variable("OccupStat") == 1) -
      2 * (variable("OccupStat") == 2) -
      (variable("OccupStat") == 9)
  )
  database <- define(
    database,
    "childCenter",
    (child_residence == 1) + (child_residence == 2) > 0
  )
  database <- define(
    database,
    "childSuburb",
    (child_residence == 3) + (child_residence == 4) > 0
  )
  database <- define(database, "TimePT_scaled", time_pt / 200)
  database <- define(database, "TimePT_hour", time_pt / 60)
  database <- define(database, "TimeCar_scaled", time_car / 200)
  database <- define(database, "TimeCar_hour", time_car / 60)
  database <- define(
    database,
    "MarginalCostPT_scaled",
    marginal_cost_pt / 10
  )
  database <- define(database, "CostCarCHF_scaled", car_cost / 10)
  database <- define(database, "distance_km_scaled", distance / 5)
  database <- define(database, "PurpHWH", variable("TripPurpose") == 1)
  database <- define(database, "PurpOther", variable("TripPurpose") != 1)
  database <- define(
    database,
    "number_of_trips",
    (trajectories == 1) + 2 * (trajectories == 2) +
      3 * (trajectories >= 3)
  )
  database <- define(
    database,
    "urbanization",
    2 - (variable("TypeCommune") <= 3) +
      (variable("TypeCommune") >= 8)
  )

  # Defining lazy derived columns resets the provisional count in the generic
  # database helper. Restore the count from the already materialized native
  # filter so callers can inspect it without forcing another data operation.
  database$filtered_row_count <- as.integer(nrow(database$data))
  database
}