## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  eval      = FALSE   # heavy computation — run manually
)

## ----packages-----------------------------------------------------------------
# library(rmdcev)
# library(apollo)
# library(dplyr)
# library(tidyr)
# library(bench)      # for timing; install.packages("bench") if needed

## ----fn-to-wide---------------------------------------------------------------
# to_wide <- function(data_long, nalts) {
#   data_long <- as.data.frame(data_long)
# 
#   # Individual-level covariates (b-cols that don't vary within id)
#   # b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific
#   all_b  <- grep("^b[0-9]", names(data_long), value = TRUE)
#   # detect individual-specific cols: same value for all rows of an id
#   ind_b  <- all_b[sapply(all_b, function(v) {
#     all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1))
#   })]
# 
#   ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b)))
# 
#   alt_vars <- data_long %>%
#     select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>%
#     pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b)))
# 
#   db <- left_join(ind_vars, alt_vars, by = "id")
# 
#   # Numeraire quantity
#   pj_xj <- Reduce("+", lapply(1:nalts, function(j) {
#     db[[paste0("price_", j)]] * db[[paste0("quant_", j)]]
#   }))
#   db$quant_num <- db$income - pj_xj
#   db
# }

## ----fn-apollo-prob-----------------------------------------------------------
# make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"),
#                               npsi_j = 1, ind_col = "b2") {
#   model <- match.arg(model)
#   alt_names <- paste0("alt", 1:nalts)
#   all_names <- c("outside", alt_names)
# 
#   function(apollo_beta, apollo_inputs, functionality = "estimate") {
#     apollo_attach(apollo_beta, apollo_inputs)
#     on.exit(apollo_detach(apollo_beta, apollo_inputs))
#     P <- list()
# 
#     # Utility for each non-numeraire alternative
#     V <- setNames(vector("list", nalts + 1), all_names)
#     V[["outside"]] <- 0
#     for (j in seq_len(nalts)) {
#       vj <- get(ind_col) * psi_ind   # individual-specific component
#       for (k in seq_len(npsi_j)) {
#         vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k))
#       }
#       V[[paste0("alt", j)]] <- vj
#     }
# 
#     # Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others
#     alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0
#     alpha <- setNames(
#       c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)),
#       all_names
#     )
# 
#     gamma <- setNames(
#       c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))),
#       all_names
#     )
# 
#     sigma <- exp(ln_sigma)
# 
#     cChoice <- setNames(
#       c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))),
#       all_names
#     )
#     cost <- setNames(
#       c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))),
#       all_names
#     )
#     avail <- setNames(
#       c(list(1), lapply(seq_len(nalts), function(j) 1)),
#       all_names
#     )
# 
#     mdcev_settings <- list(
#       alternatives     = all_names,
#       budget           = income,
#       V                = V,
#       alpha            = alpha,
#       gamma            = gamma,
#       sigma            = sigma,
#       continuousChoice = cChoice,
#       cost             = cost,
#       avail            = avail
#     )
# 
#     P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
#     P <- apollo_prepareProb(P, apollo_inputs, functionality)
#     return(P)
#   }
# }

## ----fn-apollo-beta-----------------------------------------------------------
# make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) {
#   model <- match.arg(model)
#   betas <- c(psi_ind = 0)
#   for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0
#   for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0
#   if (model == "gamma") betas["ln_alpha_out"] <- 0
#   betas["ln_sigma"] <- 0
#   betas
# }

## ----fn-fit-both--------------------------------------------------------------
# fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") {
#   nalts <- sim$data %>% distinct(alt) %>% nrow()
# 
#   # ---- rmdcev ----
#   rmdcev_formula <- as.formula(paste0(
#     "~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ")
#   ))
#   t_rmdcev <- system.time(
#     fit_r <- mdcev(
#       formula          = rmdcev_formula,
#       data             = sim$data,
#       psi_ascs         = 0,
#       model            = model,
#       algorithm        = "MLE",
#       print_iterations = FALSE,
#       backend          = "rstan"
#     )
#   )
# 
#   # ---- Apollo ----
#   database <- to_wide(sim$data, nalts)
#   db_col   <- paste0("b", npsi_j + 1)   # individual-specific column
# 
#   apollo_control <- list(
#     modelName       = paste0("apollo_", model),
#     modelDescr      = paste0("Apollo ", model),
#     indivID         = "id",
#     outputDirectory = tempdir(),
#     nCores          = 1
#   )
#   apollo_beta  <- make_apollo_beta(nalts, model, npsi_j)
#   apollo_fixed <- character(0)
# 
#   apollo_inputs <- apollo_validateInputs(
#     apollo_beta    = apollo_beta,
#     apollo_fixed   = apollo_fixed,
#     database       = database,
#     apollo_control = apollo_control
#   )
# 
#   apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col)
# 
#   t_apollo <- system.time(
#     fit_a <- apollo_estimate(
#       apollo_beta, apollo_fixed, apollo_prob, apollo_inputs,
#       estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
#                                printLevel = 0L)
#     )
#   )
# 
#   # ---- Extract estimates ----
#   r_psi   <- as.numeric(fit_r$stan_fit$par$psi)
#   r_gamma <- as.numeric(fit_r$stan_fit$par$gamma)
#   r_scale <- as.numeric(fit_r$stan_fit$par$scale)
# 
#   a_est   <- fit_a$estimate
#   a_gamma <- exp(a_est[grep("^ln_g", names(a_est))])
#   a_scale <- exp(a_est["ln_sigma"])
# 
#   # Map column names: psi_j1 = alt-varying, psi_ind = individual-varying
#   a_psi_j   <- a_est[grep("^psi_j", names(a_est))]
#   a_psi_ind <- a_est["psi_ind"]
#   a_psi     <- c(a_psi_j, a_psi_ind)
# 
#   list(
#     ll      = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout),
#     psi     = rbind(rmdcev = r_psi, apollo = a_psi),
#     gamma   = rbind(rmdcev = r_gamma[1:min(5, nalts)],
#                     apollo = a_gamma[1:min(5, nalts)]),
#     scale   = c(rmdcev = r_scale, apollo = as.numeric(a_scale)),
#     time_s  = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"])
#   )
# }

## ----param-recovery-----------------------------------------------------------
# set.seed(42)
# true_psi_j   <- 0.5
# true_psi_ind <- -1.5
# true_gamma   <- rep(5, 5)
# true_scale   <- 1.0
# 
# sim_h0 <- GenerateMDCEVData(
#   model       = "hybrid0",
#   nobs        = 500,
#   nalts       = 5,
#   psi_j_parms = true_psi_j,
#   psi_i_parms = true_psi_ind,
#   gamma_parms = true_gamma,
#   scale_parms = true_scale
# )
# 
# cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1)
# 
# cat("=== hybrid0 parameter comparison ===\n")
# cat(sprintf("LL:  rmdcev = %.2f   apollo = %.2f\n",
#             cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"]))
# cat(sprintf("     LL difference = %.6f\n",
#             abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"])))
# 
# cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n")
# print(round(cmp_h0$psi, 4))
# 
# cat("\nGamma (true: all 5)\n")
# print(round(cmp_h0$gamma, 4))
# 
# cat(sprintf("\nScale (true: 1.0):  rmdcev=%.4f  apollo=%.4f\n",
#             cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"]))

## ----param-recovery-gamma-----------------------------------------------------
# set.seed(42)
# true_alpha <- 0.5
# 
# sim_g <- GenerateMDCEVData(
#   model        = "gamma",
#   nobs         = 500,
#   nalts        = 5,
#   psi_j_parms  = true_psi_j,
#   psi_i_parms  = true_psi_ind,
#   gamma_parms  = true_gamma,
#   alpha_parms  = true_alpha,
#   scale_parms  = true_scale
# )
# 
# cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1)
# 
# cat("=== gamma parameter comparison ===\n")
# cat(sprintf("LL:  rmdcev = %.2f   apollo = %.2f\n",
#             cmp_g$ll["rmdcev"], cmp_g$ll["apollo"]))
# cat(sprintf("     LL difference = %.6f\n",
#             abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"])))
# 
# cat("\nPsi:\n"); print(round(cmp_g$psi, 4))
# cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4))
# cat(sprintf("\nScale: rmdcev=%.4f  apollo=%.4f\n",
#             cmp_g$scale["rmdcev"], cmp_g$scale["apollo"]))

## ----timing-setup-------------------------------------------------------------
# grid <- expand.grid(
#   n    = c(500, 1000, 5000),
#   J    = c(50, 100),
#   rep  = 1:3,
#   stringsAsFactors = FALSE
# ) %>%
#   filter(!(n == 1000 & J == 100))  # keep n=1000 for J=50 only

## ----timing-mle---------------------------------------------------------------
# timing_mle <- lapply(seq_len(nrow(grid)), function(i) {
#   n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i]
#   set.seed(r * 100 + i)
# 
#   sim <- GenerateMDCEVData(
#     model       = "hybrid0",
#     nobs        = n,
#     nalts       = J,
#     psi_j_parms = 0.5,
#     psi_i_parms = -1.5,
#     gamma_parms = rep(5, J),
#     scale_parms = 1.0
#   )
# 
#   res <- fit_both(sim, model = "hybrid0", npsi_j = 1)
#   data.frame(n=n, J=J, rep=r,
#              rmdcev_s = res$time_s["rmdcev"],
#              apollo_s = res$time_s["apollo"],
#              ll_diff  = abs(res$ll["rmdcev"] - res$ll["apollo"]))
# })
# 
# timing_mle <- do.call(rbind, timing_mle)

## ----timing-mle-summary-------------------------------------------------------
# summary_mle <- timing_mle %>%
#   group_by(n, J) %>%
#   summarise(
#     rmdcev_mean_s = round(mean(rmdcev_s), 1),
#     apollo_mean_s = round(mean(apollo_s), 1),
#     speedup       = round(mean(apollo_s) / mean(rmdcev_s), 1),
#     ll_diff_mean  = signif(mean(ll_diff), 3),
#     .groups = "drop"
#   )
# 
# cat("=== MLE estimation timing (seconds, 3 replicates) ===\n")
# print(summary_mle, n = Inf)

## ----sim-apollo-prob5---------------------------------------------------------
# # Concrete apollo_probabilities for J=5 (used only in this section)
# apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
#   apollo_attach(apollo_beta, apollo_inputs)
#   on.exit(apollo_detach(apollo_beta, apollo_inputs))
#   P <- list()
# 
#   V <- list(
#     outside = 0,
#     alt1    = psi_ind * b2 + psi_j1 * b1_1,
#     alt2    = psi_ind * b2 + psi_j1 * b1_2,
#     alt3    = psi_ind * b2 + psi_j1 * b1_3,
#     alt4    = psi_ind * b2 + psi_j1 * b1_4,
#     alt5    = psi_ind * b2 + psi_j1 * b1_5
#   )
#   mdcev_settings <- list(
#     alternatives     = c("outside","alt1","alt2","alt3","alt4","alt5"),
#     budget           = income, V = V,
#     alpha            = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0),
#     gamma            = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3),
#                             alt4=exp(ln_g4), alt5=exp(ln_g5)),
#     sigma            = exp(ln_sigma),
#     continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2,
#                             alt3=quant_3, alt4=quant_4, alt5=quant_5),
#     cost             = list(outside=1, alt1=price_1, alt2=price_2,
#                             alt3=price_3, alt4=price_4, alt5=price_5),
#     avail            = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1)
#   )
#   P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
#   P <- apollo_prepareProb(P, apollo_inputs, functionality)
#   return(P)
# }

## ----sim-convergence----------------------------------------------------------
# set.seed(42)
# J5   <- 5
# n500 <- 500
# 
# sim5 <- GenerateMDCEVData(
#   model       = "hybrid0",
#   nobs        = n500,
#   nalts       = J5,
#   psi_j_parms = 0.5,
#   psi_i_parms = -1.5,
#   gamma_parms = rep(5, J5),
#   scale_parms = 1.0
# )
# 
# # ── Fit both packages ──────────────────────────────────────────────────────
# fit_r5 <- mdcev(
#   formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0,
#   model = "hybrid0", algorithm = "MLE",
#   print_iterations = FALSE, backend = "rstan"
# )
# 
# database5 <- to_wide(sim5$data, J5)
# ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id",
#             outputDirectory = tempdir(), nCores = 1)
# ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1)
# ai5 <- suppressMessages(
#   apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0),
#                         database = database5, apollo_control = ac5))
# fit_a5 <- suppressMessages(suppressWarnings(
#   apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5,
#     estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
#                              printLevel = 0L))
# ))
# 
# # ── 1% price-change policy ─────────────────────────────────────────────────
# baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500)
# price_change    <- 0.01 * colMeans(baseline_prices)   # 1% of mean price per alt
# 
# # rmdcev: additive policy vector (0 for numeraire, then per-alt changes)
# policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE)
# policies5$price_p[[2]] <- c(0, price_change)
# 
# df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1)
# 
# t_rmdcev5 <- system.time(
#   demand_r5 <- mdcev.sim(
#     df_sim5$df_indiv, df_common = df_sim5$df_common,
#     sim_options = df_sim5$sim_options,
#     cond_err = 0, nerrs = 100, sim_type = "demand"
#   )
# )
# 
# # Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index]
# mean_demand <- function(demand_list, policy_idx, alt_idx) {
#   mean(sapply(demand_list, function(ind)
#     mean(sapply(ind, function(s) s[policy_idx, alt_idx]))))
# }
# d_base_r   <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a))
# d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a))
# delta_r    <- d_policy_r - d_base_r
# 
# # Apollo: predict at baseline and policy databases
# # Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt)
# database5_pol <- database5
# for (j in 1:J5) {
#   database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j]
# }
# # Recompute quant_num at observed quantities under new prices
# pj_xj_pol <- Reduce("+", lapply(1:J5, function(j)
#   database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]]))
# database5_pol$quant_num <- database5_pol$income - pj_xj_pol
# 
# ai5_pol <- suppressMessages(
#   apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0),
#                         database = database5_pol, apollo_control = ac5))
# 
# t_apollo5_base <- system.time(
#   pred_base <- suppressMessages(suppressWarnings(
#     apollo_prediction(fit_a5, apollo_probabilities_j5, ai5,
#                       prediction_settings = list(nRep = 100, silent = TRUE))
#   ))
# )
# t_apollo5_pol <- system.time(
#   pred_pol  <- suppressMessages(suppressWarnings(
#     apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol,
#                       prediction_settings = list(nRep = 100, silent = TRUE))
#   ))
# )
# t_apollo5 <- t_apollo5_base + t_apollo5_pol
# 
# delta_a <- sapply(1:J5, function(j) {
#   mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) -
#   mean(pred_base[[paste0("alt", j, "_cont_mean")]])
# })
# 
# cat("=== Demand change: 1% price increase (n=500, J=5) ===\n")
# cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff"))
# for (j in 1:J5) {
#   cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j],
#               abs(delta_r[j] - delta_a[j])))
# }
# cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a))))
# cat(sprintf("rmdcev %% change: %s%%\n",
#             paste(round(100 * delta_r / d_base_r, 2), collapse = ", ")))
# 
# cat(sprintf("\nTiming:  rmdcev %.2f s  |  Apollo (2 predictions) %.2f s\n",
#             t_rmdcev5["elapsed"], t_apollo5["elapsed"]))

## ----make-concrete-prob-------------------------------------------------------
# # Dynamically generate a concrete apollo_probabilities function for any J.
# # apollo_prediction() does source-code inspection and cannot use factory closures,
# # so we build and eval() a function body with the right number of alternatives.
# make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") {
#   all_str     <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")")
#   V_items     <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col,
#                               " + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n    ")
#   alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ")
#   gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ")
#   choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ")
#   cost_items  <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ")
#   avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ")
# 
#   eval(parse(text = sprintf(
# 'function(apollo_beta, apollo_inputs, functionality = "estimate") {
#   apollo_attach(apollo_beta, apollo_inputs)
#   on.exit(apollo_detach(apollo_beta, apollo_inputs))
#   P <- list()
#   V <- list(outside = 0, %s)
#   mdcev_settings <- list(
#     alternatives = %s, budget = income, V = V,
#     alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma),
#     continuousChoice = list(%s), cost = list(%s), avail = list(%s)
#   )
#   P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
#   P <- apollo_prepareProb(P, apollo_inputs, functionality)
#   return(P)
# }', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items)))
# }

## ----timing-sim-grid----------------------------------------------------------
# timing_sim <- lapply(c(50, 100), function(J) {
#   n <- 500
#   set.seed(J)
# 
#   sim <- GenerateMDCEVData(
#     model = "hybrid0", nobs = n, nalts = J,
#     psi_j_parms = 0.5, psi_i_parms = -1.5,
#     gamma_parms = rep(5, J), scale_parms = 1.0
#   )
# 
#   # ── rmdcev ──────────────────────────────────────────────────────────────
#   fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0",
#                  algorithm = "MLE", print_iterations = FALSE, backend = "rstan")
#   bp   <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n)
#   pc   <- 0.01 * colMeans(bp)
#   pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE)
#   pols$price_p[[2]] <- c(0, pc)
#   df_s <- PrepareSimulationData(fit_r, pols, nsims = 1)
# 
#   rmdcev_times <- sapply(1:3, function(r)
#     system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common,
#                           sim_options = df_s$sim_options,
#                           cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"])
# 
#   # ── Apollo ──────────────────────────────────────────────────────────────
#   database <- to_wide(sim$data, J)
#   ac  <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id",
#               outputDirectory = tempdir(), nCores = 1)
#   ab  <- make_apollo_beta(J, npsi_j = 1)
#   pfn <- make_apollo_prob_concrete(J)
#   ai  <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0),
#                                database = database, apollo_control = ac)
#   fit_a <- apollo_estimate(ab, character(0), pfn, ai,
#              estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
#                                       printLevel = 0L))
# 
#   # Policy database: 1% additive price increase
#   db_pol <- database
#   for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j]
#   pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j)
#     db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]]))
#   db_pol$quant_num <- db_pol$income - pj_xj_pol
#   ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0),
#                                   database = db_pol, apollo_control = ac)
# 
#   apollo_times <- sapply(1:3, function(r) {
#     t1 <- system.time(apollo_prediction(fit_a, pfn, ai,
#            prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
#     t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol,
#            prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
#     t1 + t2
#   })
# 
#   data.frame(n = n, J = J,
#              rmdcev_mean = round(mean(rmdcev_times), 2),
#              rmdcev_sd   = round(sd(rmdcev_times),   2),
#              apollo_mean = round(mean(apollo_times),  2),
#              apollo_sd   = round(sd(apollo_times),    2),
#              speedup     = round(mean(apollo_times) / mean(rmdcev_times), 1))
# })
# 
# timing_sim <- do.call(rbind, timing_sim)
# print(timing_sim, row.names = FALSE)

## ----session------------------------------------------------------------------
# sessionInfo()

