R/RcppExports.R

In choicer: Discrete Choice Models for Economic Applications

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Radical inverse (van der Corput) for testing halton.h
#'
#' @param n Sequence index (coerced to uint64_t).
#' @param base Prime base (coerced to uint32_t).
#' @return Radical inverse value in [0, 1).
#' @noRd
halton_radical_inverse <- function(n, base) {
    .Call(`_choicer_halton_radical_inverse`, n, base)
}

#' Wichura AS241 inverse normal CDF for testing halton.h
#'
#' @param p Probability in (0, 1).
#' @return Quantile value.
#' @noRd
halton_inv_normal_cdf <- function(p) {
    .Call(`_choicer_halton_inv_normal_cdf`, p)
}

#' Generate an n x dim matrix of uniform digit-permuted Halton draws for testing
#'
#' Returns an n x dim matrix using global indices 1..n (one row per index,
#' one column per dimension). For scramble=0 (compat mode) the result
#' reproduces randtoolbox::halton(n, dim, normal=FALSE) exactly.
#'
#' @param n Number of Halton points (rows).
#' @param dim Number of dimensions (columns).
#' @param seed Master seed for position-wise digit permutations (coerced to uint64_t).
#'   Ignored when scramble=0.
#' @param scramble 0 = identity (compat), 1 = position-wise digit permutation.
#' @return n x dim arma::mat of uniform [0,1) values.
#' @noRd
halton_generate_uniform <- function(n, dim, seed, scramble) {
    .Call(`_choicer_halton_generate_uniform`, n, dim, seed, scramble)
}

#' Generate a K_w x (S*N) matrix of normal draws for testing halton.h
#'
#' Built via HaltonGen::fill_eta_i for i=1..N.
#'
#' Layout: columns `[(i-1)*S, i*S)` hold eta_i (K_w x S) for individual i.
#' Within individual i, column s holds the K_w variates for draw s (0-based),
#' so `out(k, (i-1)*S + s) = inv_normal_cdf(phi_{PRIMES[k]}((i-1)*S + s + 1))`.
#'
#' @param S   Number of draws per individual.
#' @param N   Number of individuals.
#' @param K_w Number of random-coefficient dimensions.
#' @param seed Master seed for position-wise digit permutations (coerced to uint64_t).
#' @param scramble 0 = identity (compat), 1 = position-wise digit permutation.
#' @return K_w x (S*N) arma::mat of standard-normal draws.
#' @noRd
halton_generate_normal <- function(S, N, K_w, seed, scramble) {
    .Call(`_choicer_halton_generate_normal`, S, N, K_w, seed, scramble)
}

#' Hand-rolled lower Cholesky for testing hb_internal.h
#'
#' @param A Symmetric matrix.
#' @return List with `ok` (FALSE on a non-positive pivot) and `L` (lower
#'   Cholesky factor; zero-filled when `ok` is FALSE).
#' @noRd
hb_test_chol <- function(A) {
    .Call(`_choicer_hb_test_chol`, A)
}

#' Hand-rolled triangular solve for testing hb_internal.h
#'
#' Solves L x = b (`transpose = FALSE`, forward substitution) or L' x = b
#' (`transpose = TRUE`, back substitution) for a lower-triangular L.
#'
#' @param L Lower-triangular matrix.
#' @param b Right-hand side vector.
#' @param transpose Solve against L' instead of L.
#' @return Solution vector.
#' @noRd
hb_test_trisolve <- function(L, b, transpose) {
    .Call(`_choicer_hb_test_trisolve`, L, b, transpose)
}

#' SPD solve via hand-rolled Cholesky + two trisolves for testing
#'
#' @param A Symmetric positive-definite matrix.
#' @param b Right-hand side vector.
#' @return List with `ok` (FALSE when `A` is not SPD) and `x` (solution;
#'   zero-filled when `ok` is FALSE).
#' @noRd
hb_test_spd_solve <- function(A, b) {
    .Call(`_choicer_hb_test_spd_solve`, A, b)
}

#' Fixed-order log-sum-exp with optional implicit outside term for testing
#'
#' @param v Vector of utilities.
#' @param include_outside Add the implicit `exp(0)` outside-option term.
#' @return `log(sum(exp(v)))`, plus 1 inside the sum when
#'   `include_outside = TRUE`.
#' @noRd
hb_test_logsumexp <- function(v, include_outside) {
    .Call(`_choicer_hb_test_logsumexp`, v, include_outside)
}

#' Two-block Gibbs for sigma_d2 under the half-Cauchy scale mixture
#'
#' Runs `n_iter` sweeps of draw_sigma_d2_conditional() on fixed residuals
#' `xi`, starting from sigma_d2 = a_d = 1, one Xoshiro stream per iteration
#' (tag 0). Exercises both the half-Cauchy Makalic-Schmidt mixture
#' (`half_cauchy = TRUE`, scale `s_d`) and the plain IG(c0, d0) fallback.
#'
#' @param xi Vector of alternative-level residuals.
#' @param n_iter Number of Gibbs sweeps.
#' @param seed Master seed (coerced to uint64_t).
#' @param half_cauchy Use the half-Cauchy scale mixture (else IG fallback).
#' @param s_d Half-Cauchy scale.
#' @param c0 IG fallback shape.
#' @param d0 IG fallback scale.
#' @return Numeric vector of `n_iter` sigma_d2 draws.
#' @noRd
hb_test_sigma_d2_gibbs <- function(xi, n_iter, seed, half_cauchy, s_d, c0, d0) {
    .Call(`_choicer_hb_test_sigma_d2_gibbs`, xi, n_iter, seed, half_cauchy, s_d, c0, d0)
}

#' Gibbs sampler for the hierarchical Bayesian multinomial logit model
#'
#' Runs the adaptive random-walk Metropolis-within-Gibbs sampler for the
#' hierarchical (random-coefficients, panel) multinomial logit with a
#' BLP-style alternative-level random effect: inside utilities
#' \eqn{U_{ijt} = x_{ijt}'\gamma_i + \delta_j + EV1} against an implicit
#' outside option with systematic utility 0, \eqn{\beta_i \sim N(b, W)}
#' (\eqn{\gamma_{ik} = \beta_{ik}} or \eqn{\exp(\beta_{ik})} per
#' \code{rc_dist}), and \eqn{\delta_j = z_j'\theta + \xi_j},
#' \eqn{\xi_j \sim N(0, \sigma_d^2)}.
#'
#' The per-respondent \eqn{\beta_i} updates are parallelized with OpenMP;
#' the \eqn{\delta_j} updates run as a strictly serial sweep (their
#' conditionals are coupled through the shared softmax denominators). Each
#' (iteration, unit) pair uses its own RNG stream, so draws are
#' reproducible given the seed and a fixed thread count; across different
#' thread counts they are invariant only up to floating-point
#' reduction-order round-off (~1e-15), not bitwise (see
#' \code{set_num_threads()}). This is the low-level engine behind
#' \code{\link{run_hmnlogit}}, which handles initialization and
#' post-processing.
#'
#' @param X total_rows x K_struct structural design matrix (inside rows
#'   only, no ASC columns), rows sorted by (person, task, alternative).
#' @param Z J x P alternative-level mean-function design (intercept first).
#' @param M Integer vector: inside alternatives per choice situation.
#' @param choice_pos Integer vector: 1-based within-task position of the
#'   chosen row; 0 = outside option chosen.
#' @param include_outside_option Must be \code{TRUE} (the implicit outside
#'   good anchors the location of delta; a no-outside mode is roadmapped).
#' @param alt_of_row Integer vector: 1-based alternative code per row of X.
#' @param Ti Integer vector: choice situations per respondent.
#' @param rc_dist Integer vector (length K_struct): 0 = normal coordinate,
#'   1 = log-normal (enters utility as \code{exp(beta_ik)}).
#' @param beta_pooled Pooled MNL MLE on the chain scale (log scale for
#'   log-normal coordinates); centers the H_i proposal information.
#' @param delta_init Initial delta (length J).
#' @param theta_init Initial theta (length P).
#' @param b_bar K vector, prior mean of b.
#' @param A K x K prior precision matrix of b.
#' @param nu Inverse-Wishart prior degrees of freedom for W (>= K).
#' @param V K x K inverse-Wishart prior scale matrix for W.
#' @param theta_bar P vector, prior mean of theta.
#' @param A_theta P x P prior precision matrix of theta.
#' @param sd_prior List with elements \code{half_cauchy} (logical),
#'   \code{s_d} (half-Cauchy scale), \code{c0}, \code{d0} (IG fallback).
#' @param R Total number of Gibbs iterations.
#' @param burn Number of initial iterations discarded (0 <= burn < R);
#'   proposal-scale adaptation happens during burn-in only.
#' @param thin Keep every thin-th post-burn-in draw.
#' @param seed Master RNG seed (non-negative; all streams derive from it).
#' @param keep_beta_i 0 = no beta_i output, 1 = online means/SDs,
#'   2 = means/SDs plus the full (K, N, R_keep) draw cube.
#' @param s_init Initial per-respondent proposal scale.
#' @param accept_target Robbins-Monro acceptance target for the beta_i
#'   updates (the delta_j target is fixed at 0.44).
#' @param trace Print progress every \code{trace} iterations (0 = silent).
#' @returns List with \code{bdraw} (R_keep x K), \code{wdraw} (R_keep x
#'   K(K+1)/2, lower triangle of W in row-major order), \code{deltadraw}
#'   (R_keep x J), \code{thetadraw} (R_keep x P), \code{sigma_d2draw},
#'   \code{loglik_trace}, acceptance rates and final proposal scales
#'   (\code{accept_rate_beta}, \code{accept_rate_delta}, \code{s_final},
#'   \code{s_delta_final}), posterior summaries \code{beta_i_mean} /
#'   \code{beta_i_sd} (K x N, \code{NULL} when \code{keep_beta_i = 0}),
#'   \code{beta_i_draws} (K x N x R_keep cube when \code{keep_beta_i = 2}),
#'   \code{delta_mean} / \code{delta_sd} / \code{xi_mean} / \code{xi_sd}
#'   (J x 1), and \code{R_keep}.
#' @examples
#' \donttest{
#' sim <- simulate_hmnl_data(N = 20, T = 2, J = 3, seed = 42)
#' d <- prepare_hmnl_data(sim$data, "task", "alt", "choice",
#'                        c("x1", "x2"), person_col = "pid")
#' out <- choicer:::hmnl_gibbs(d$X, d$Z, d$M, d$choice_pos, TRUE, d$alt_of_row, d$Ti,
#'   rc_dist = d$rc_dist, beta_pooled = rep(0, d$K_struct),
#'   delta_init = rep(0, d$J), theta_init = rep(0, d$P),
#'   b_bar = rep(0, d$K_struct), A = 0.01 * diag(d$K_struct),
#'   nu = d$K_struct + 3, V = (d$K_struct + 3) * diag(d$K_struct),
#'   theta_bar = rep(0, d$P), A_theta = 0.01 * diag(d$P),
#'   sd_prior = list(half_cauchy = TRUE, s_d = 1, c0 = 3, d0 = 3),
#'   R = 300, burn = 100, thin = 1, seed = 7, keep_beta_i = 1,
#'   s_init = 2.38 / sqrt(d$K_struct), accept_target = 0.234)
#' colMeans(out$bdraw)
#' }
#' @keywords internal
hmnl_gibbs <- function(X, Z, M, choice_pos, include_outside_option, alt_of_row, Ti, rc_dist, beta_pooled, delta_init, theta_init, b_bar, A, nu, V, theta_bar, A_theta, sd_prior, R, burn, thin, seed, keep_beta_i, s_init, accept_target, trace = 0L) {
    .Call(`_choicer_hmnl_gibbs`, X, Z, M, choice_pos, include_outside_option, alt_of_row, Ti, rc_dist, beta_pooled, delta_init, theta_init, b_bar, A, nu, V, theta_bar, A_theta, sd_prior, R, burn, thin, seed, keep_beta_i, s_init, accept_target, trace)
}

#' Gibbs sampler for the hierarchical Bayesian multinomial probit model
#'
#' Runs the fully conjugate Albert-Chib Gibbs sampler for the hierarchical
#' multinomial probit with iid \eqn{N(0, \sigma^2)} utility shocks in
#' un-differenced utility space: inside utilities
#' \eqn{U_{ijt} = x_{ijt}'\beta_i + \delta_j + \epsilon} against a
#' stochastic implicit outside option \eqn{U_{iot} = \epsilon},
#' \eqn{\beta_i \sim N(b, W)}, and \eqn{\delta_j = z_j'\theta + \xi_j},
#' \eqn{\xi_j \sim N(0, \sigma_d^2)}. The chain runs on the non-identified
#' parameterization (free \eqn{\sigma^2}); identified quantities are
#' obtained by normalizing each draw by the matching power of \eqn{\sigma}
#' (handled by \code{\link{run_hmnprobit}}).
#'
#' The latent sweep and the \eqn{\beta_i} draws are parallelized with
#' OpenMP across respondents; the \eqn{\delta_j} draws are parallelized
#' across alternatives (conditionally independent given the augmented
#' utilities — unlike the HMNL, whose delta sweep must be serial). Each
#' (iteration, unit) pair uses its own RNG stream, so draws are
#' reproducible given the seed and a fixed thread count; across different
#' thread counts they are invariant only up to floating-point
#' reduction-order round-off (~1e-15), not bitwise.
#'
#' @param X total_rows x K_struct structural design matrix (inside rows
#'   only), rows sorted by (person, task, alternative).
#' @param Z J x P alternative-level mean-function design (intercept first).
#' @param M Integer vector: inside alternatives per choice situation.
#' @param choice_pos Integer vector: 1-based within-task position of the
#'   chosen row; 0 = outside option chosen.
#' @param include_outside_option Must be \code{TRUE} (the outside good
#'   anchors the location of delta; a no-outside mode is roadmapped).
#' @param alt_of_row Integer vector: 1-based alternative code per row of X.
#' @param Ti Integer vector: choice situations per respondent.
#' @param delta_init Initial delta (length J), raw scale.
#' @param theta_init Initial theta (length P), raw scale.
#' @param b_bar K vector, prior mean of b.
#' @param A K x K prior precision matrix of b.
#' @param nu Inverse-Wishart prior degrees of freedom for W (>= K).
#' @param V K x K inverse-Wishart prior scale matrix for W.
#' @param theta_bar P vector, prior mean of theta.
#' @param A_theta P x P prior precision matrix of theta.
#' @param sd_prior List with elements \code{half_cauchy} (logical),
#'   \code{s_d}, \code{c0}, \code{d0} — the sigma_d prior (see
#'   \code{\link{run_hmnlogit}}).
#' @param a0,s0 Inverse-gamma prior shape/scale for the (non-identified)
#'   shock variance \eqn{\sigma^2}.
#' @param R Total number of Gibbs iterations.
#' @param burn Number of initial iterations discarded (0 <= burn < R).
#' @param thin Keep every thin-th post-burn-in draw.
#' @param seed Master RNG seed (non-negative; all streams derive from it).
#' @param keep_beta_i 0 = no beta_i output, 1 = online means/SDs
#'   (identified scale), 2 = means/SDs plus the full (K, N, R_keep) cube of
#'   per-draw-normalized \eqn{\beta_i / \sigma} draws.
#' @param trace Print progress every \code{trace} iterations (0 = silent).
#' @returns List with RAW draw matrices \code{bdraw}, \code{wdraw} (lower
#'   triangle, row-major), \code{deltadraw}, \code{thetadraw},
#'   \code{sigma_d2draw}, \code{sigma2draw}, identified-scale summaries
#'   \code{beta_i_mean} / \code{beta_i_sd} / \code{beta_i_draws} /
#'   \code{delta_mean} / \code{delta_sd} / \code{xi_mean} / \code{xi_sd},
#'   and \code{R_keep}.
#' @examples
#' \donttest{
#' sim <- simulate_hmnp_data(N = 30, T = 2, J = 3, seed = 42)
#' d <- prepare_hmnp_data(sim$data, "task", "alt", "choice",
#'                        c("x1", "x2"), person_col = "pid")
#' out <- choicer:::hmnp_gibbs(d$X, d$Z, d$M, d$choice_pos, TRUE, d$alt_of_row, d$Ti,
#'   delta_init = rep(0, d$J), theta_init = rep(0, d$P),
#'   b_bar = rep(0, d$K_struct), A = 0.01 * diag(d$K_struct),
#'   nu = d$K_struct + 3, V = (d$K_struct + 3) * diag(d$K_struct),
#'   theta_bar = rep(0, d$P), A_theta = 0.01 * diag(d$P),
#'   sd_prior = list(half_cauchy = TRUE, s_d = 1, c0 = 3, d0 = 3),
#'   a0 = 3, s0 = 3, R = 300, burn = 100, thin = 1, seed = 7,
#'   keep_beta_i = 1)
#' colMeans(out$bdraw / sqrt(as.numeric(out$sigma2draw)))
#' }
#' @keywords internal
hmnp_gibbs <- function(X, Z, M, choice_pos, include_outside_option, alt_of_row, Ti, delta_init, theta_init, b_bar, A, nu, V, theta_bar, A_theta, sd_prior, a0, s0, R, burn, thin, seed, keep_beta_i, trace = 0L) {
    .Call(`_choicer_hmnp_gibbs`, X, Z, M, choice_pos, include_outside_option, alt_of_row, Ti, delta_init, theta_init, b_bar, A, nu, V, theta_bar, A_theta, sd_prior, a0, s0, R, burn, thin, seed, keep_beta_i, trace)
}

#' Log-likelihood and gradient for multinomial logit model
#'
#' Computes the log-likelihood and its gradient for the Multinomial Logit model using OpenMP for parallelization.
#' Allows for inclusion of alternative-specific constants, outside option, and observation weights.
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates. Stacks M\[i] x K matrices for individual i.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @returns List with loglikelihood and gradient evaluated at input arguments
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mnl_data(dt, "id", "alt", "choice", c("x1", "x2"))
#' theta <- rep(0, ncol(d$X) + nrow(d$alt_mapping) - 1)
#' result <- choicer:::mnl_loglik_gradient_parallel(theta, d$X, d$alt_idx,
#'   d$choice_idx, d$M, d$weights)
#' result$objective  # negative log-likelihood
#' }
#' @keywords internal
mnl_loglik_gradient_parallel <- function(theta, X, alt_idx, choice_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_loglik_gradient_parallel`, theta, X, alt_idx, choice_idx, M, weights, use_asc, include_outside_option)
}

#' BHHH/OPG information matrix for multinomial logit model
#'
#' Computes the weighted outer product of per-individual scores
#' \eqn{\sum_i w_i\, s_i s_i^\top} for the Multinomial Logit model. The
#' per-individual score \eqn{s_i} is the (positive) gradient of individual
#' \eqn{i}'s log-likelihood contribution and is weight-free; the supplied
#' \code{weights} enter only as the leading multiplier. Passing
#' \code{weights = w} yields the ordinary weighted BHHH/OPG information; passing
#' \code{weights = w^2} yields the sandwich \emph{meat}
#' \eqn{B = \sum_i w_i^2 s_i s_i^\top} used for robust (WESML) inference.
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates. Stacks M\[i] x K matrices for individual i.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @returns A symmetric positive-semidefinite information matrix
#'   \eqn{\sum_i w_i\, s_i s_i^\top} (same sign convention as the negated Hessian).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' B <- choicer:::mnl_bhhh_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$choice_idx, fit$data$M, fit$data$weights)
#' dim(B)
#' }
#' @keywords internal
mnl_bhhh_parallel <- function(theta, X, alt_idx, choice_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_bhhh_parallel`, theta, X, alt_idx, choice_idx, M, weights, use_asc, include_outside_option)
}

mnl_scores_parallel <- function(theta, X, alt_idx, choice_idx, M, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_scores_parallel`, theta, X, alt_idx, choice_idx, M, use_asc, include_outside_option)
}

#' Prediction of choice probabilities and utilities based on fitted model
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates. Stacks M\[i] x K matrices for individual i.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @returns List with choice probability and utility for each choice situation evaluated at input arguments
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' pred <- choicer:::mnl_predict(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$M, use_asc = TRUE)
#' head(pred$choice_prob)
#' }
#' @keywords internal
mnl_predict <- function(theta, X, alt_idx, M, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_predict`, theta, X, alt_idx, M, use_asc, include_outside_option)
}

#' Prediction of market shares based on fitted model
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates. Stacks M\[i] x K matrices for individual i.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @returns vector with predicted market shares for each alternative
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' shares <- choicer:::mnl_predict_shares(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$M, fit$data$weights, use_asc = TRUE)
#' shares
#' }
#' @keywords internal
mnl_predict_shares <- function(theta, X, alt_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_predict_shares`, theta, X, alt_idx, M, weights, use_asc, include_outside_option)
}

#' BLP95 contraction mapping to find delta given target shares
#'
#' @param delta J x 1 vector with initial guess for deltas (ASCs)
#' @param target_shares J x 1 vector with target shares for each alternative
#' @param X sum(M) x K design matrix with covariates. M\[i] x K matrix for individual i
#' @param beta K x 1 vector with model parameters
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @param tol convergence tolerance
#' @param max_iter maximum number of iterations
#' @returns vector with contraction's delta (ASCs) output
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' beta <- coef(fit)[fit$param_map$beta]
#' delta <- blp_contraction(rep(0, J), rep(1/J, J), fit$data$X,
#'   beta, fit$data$alt_idx, fit$data$M, fit$data$weights)
#' delta
#' }
#' @export
blp_contraction <- function(delta, target_shares, X, beta, alt_idx, M, weights, include_outside_option = FALSE, tol = 1e-8, max_iter = 1000L) {
    .Call(`_choicer_blp_contraction`, delta, target_shares, X, beta, alt_idx, M, weights, include_outside_option, tol, max_iter)
}

#' Hessian matrix for multinomial logit model
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates. Stacks M\[i] x K matrices for individual i.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @returns Hessian matrix of the negative log-likelihood
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' H <- choicer:::mnl_loglik_hessian_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$choice_idx, fit$data$M, fit$data$weights)
#' dim(H)
#' }
#' @keywords internal
mnl_loglik_hessian_parallel <- function(theta, X, alt_idx, choice_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_loglik_hessian_parallel`, theta, X, alt_idx, choice_idx, M, weights, use_asc, include_outside_option)
}

#' Compute aggregate elasticities for MNL model
#'
#' Computes the aggregate elasticity matrix (weighted average of individual
#' elasticities) for the Multinomial Logit model.
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param choice_idx N x 1 vector (kept for API consistency, but not used)
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param elast_var_idx 1-based index of the column in X for which to compute the elasticity
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option
#' @returns J x J matrix of aggregate elasticities
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' elas <- choicer:::mnl_elasticities_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$choice_idx, fit$data$M, fit$data$weights, elast_var_idx = 1L)
#' elas
#' }
#' @keywords internal
mnl_elasticities_parallel <- function(theta, X, alt_idx, choice_idx, M, weights, elast_var_idx, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_elasticities_parallel`, theta, X, alt_idx, choice_idx, M, weights, elast_var_idx, use_asc, include_outside_option)
}

#' Compute MNL diversion ratios (parallelized over individuals)
#'
#' Computes the diversion ratio matrix DR(j->k), which measures the fraction
#' of demand lost by alternative j that is captured by alternative k.
#' For MNL: DR(j->k) = sum_n(w_n * P_nj * P_nk) / sum_n(w_n * P_nj * (1 - P_nj))
#'
#' @param theta K + J - 1 or K + J vector with model parameters
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param use_asc whether to use alternative-specific constants
#' @param include_outside_option whether to include outside option
#' @returns J x J matrix where entry (k, j) = DR(j->k). Diagonal is 0.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_mnlogit(dt, "id", "alt", "choice", c("x1", "x2"))
#' dr <- choicer:::mnl_diversion_ratios_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$M, fit$data$weights)
#' dr
#' }
#' @keywords internal
mnl_diversion_ratios_parallel <- function(theta, X, alt_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_mnl_diversion_ratios_parallel`, theta, X, alt_idx, M, weights, use_asc, include_outside_option)
}

#' Gibbs sampler for the Bayesian multinomial probit model
#'
#' Runs the McCulloch-Rossi (1994) Gibbs sampler with Albert-Chib data
#' augmentation for the multinomial probit model in utility differences
#' against a base alternative. The chain operates on the non-identified
#' parameterization (unrestricted \code{Sigma}); identified quantities are
#' obtained by normalizing each draw by \code{sigma_11} (handled by
#' \code{\link{run_mnprobit}}).
#'
#' The latent-utility sweep is parallelized with OpenMP across choice
#' situations (they are conditionally independent given \code{beta} and
#' \code{Sigma}). Each (iteration, observation) pair uses its own RNG
#' stream, so draws are reproducible given the seed and a fixed thread
#' count; across different thread counts they are invariant only up to
#' floating-point reduction-order round-off (~1e-15), not bitwise
#' (see \code{set_num_threads()}).
#'
#' @param X (N*p) x K stacked design matrix of utility differences. Rows are
#'   grouped by choice situation, with the p = J - 1 difference rows of
#'   situation i ordered by alternative.
#' @param y N vector of choices: 0 for the base alternative, j in 1..p for
#'   the j-th non-base alternative.
#' @param p Number of utility differences (J - 1).
#' @param beta_bar K vector, prior mean of beta.
#' @param A K x K prior precision matrix of beta.
#' @param nu Inverse-Wishart prior degrees of freedom (>= p).
#' @param V p x p inverse-Wishart prior scale matrix.
#' @param R Total number of Gibbs iterations.
#' @param burn Number of initial iterations to discard (0 <= burn < R).
#' @param thin Keep every thin-th post-burn-in draw.
#' @param seed Master RNG seed (non-negative; all streams derive from it).
#' @param trace Print progress every \code{trace} iterations (0 = silent).
#' @returns List with \code{betadraw} (R_keep x K), \code{sigmadraw}
#'   (R_keep x p(p+1)/2, lower triangle of Sigma in row-major order), and
#'   \code{R_keep}.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 100; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mnp_data(dt, "id", "alt", "choice", c("x1", "x2"))
#' out <- choicer:::mnp_gibbs(d$X, d$y, d$p,
#'   beta_bar = rep(0, d$K), A = 0.01 * diag(d$K),
#'   nu = d$p + 3, V = (d$p + 3) * diag(d$p),
#'   R = 500, burn = 100, thin = 1, seed = 42)
#' colMeans(out$betadraw)
#' }
#' @keywords internal
mnp_gibbs <- function(X, y, p, beta_bar, A, nu, V, R, burn, thin, seed, trace = 0L) {
    .Call(`_choicer_mnp_gibbs`, X, y, p, beta_bar, A, nu, V, R, burn, thin, seed, trace)
}

rnorm_cpp <- function(n, seed) {
    .Call(`_choicer_rnorm_cpp`, n, seed)
}

rgamma_cpp <- function(n, a, seed) {
    .Call(`_choicer_rgamma_cpp`, n, a, seed)
}

rtruncnorm_cpp <- function(n, mu, sigma, a, b, seed) {
    .Call(`_choicer_rtruncnorm_cpp`, n, mu, sigma, a, b, seed)
}

rmvnorm_cpp <- function(n, mu, Sigma, seed) {
    .Call(`_choicer_rmvnorm_cpp`, n, mu, Sigma, seed)
}

rwishart_cpp <- function(df, S, seed) {
    .Call(`_choicer_rwishart_cpp`, df, S, seed)
}

riwishart_cpp <- function(df, V, seed) {
    .Call(`_choicer_riwishart_cpp`, df, V, seed)
}

build_L_mat <- function(L_params, K_w, rc_correlation) {
    .Call(`_choicer_build_L_mat`, L_params, K_w, rc_correlation)
}

#' Reconstruct variance matrix L from L_params
#'
#' @param L_params flattened choleski decomposition version of the random coefficient parameters matrix
#' @param K_w dimension of the random coefficient parameter (symmetric) matrix
#' @param rc_correlation whether random coefficients are correlated
#' @returns matrix equal to LL', where L is the choleski decomposition of random coefficient matrix
#' @examples
#' L_params <- c(log(1.0), 0.3, log(0.5))
#' Sigma <- choicer:::build_var_mat(L_params, K_w = 2, rc_correlation = TRUE)
#' Sigma  # 2x2 covariance matrix
#' @keywords internal
build_var_mat <- function(L_params, K_w, rc_correlation) {
    .Call(`_choicer_build_var_mat`, L_params, K_w, rc_correlation)
}

#' Log-likelihood and gradient for Mixed Logit
#'
#' Computes the log-likelihood and its gradient for the Mixed Logit model using
#' OpenMP for parallelization. Allows for inclusion of alternative-specific
#' constants, outside option, observation weights, correlated random coefficients.
#'
#' @param theta vector collecting model parameters (beta, mu, L, delta (ASCs))
#' @param X design matrix for covariates with fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for covariates with random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with choice situation draws; K_w x S x N
#' @param rc_dist K_w x 1 integer vector indicating distribution of random coefficients: 0 = normal, 1 = log-normal
#' @param rc_correlation whether random coefficients should be correlated
#' @param rc_mean whether to estimate means for random coefficients. If so, mean parameters (mu) should be included in theta after beta parameters.
#' @param use_asc whether to use alternative-specific constants. If so, parameters should be included in theta after beta and L (and mu, if applicable).
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns List with loglikelihood and gradient evaluated at input arguments
#' @note For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE,
#'   the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta),
#'   where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2).
#'   This differs from the textbook parameterization exp(mu_k + L_k * eta).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' K_x <- ncol(d$X); K_w <- ncol(d$W); J <- nrow(d$alt_mapping)
#' theta <- rep(0, K_x + K_w + J - 1)
#' result <- choicer:::mxl_loglik_gradient_parallel(theta, d$X, d$W, d$alt_idx,
#'   d$choice_idx, d$M, d$weights, eta, rc_dist = rep(0L, K_w),
#'   rc_correlation = FALSE, rc_mean = FALSE)
#' result$objective
#' }
#' @keywords internal
mxl_loglik_gradient_parallel <- function(theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_loglik_gradient_parallel`, theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Utility to compute analytical Jacobian of random coefficient matrix transformed by vech (dVech(Sigma) / dTheta)
#'
#' @param L_params flattened choleski decomposition version of the random coefficient parameters matrix
#' @param K_w dimension of the random coefficient parameter (symmetric) matrix
#' @param rc_correlation whether random coefficients are correlated
#' @returns Jacobian (dVech(Sigma) / dTheta)
#' @examples
#' L_params <- c(log(0.8), 0.2, log(0.6))
#' J_mat <- choicer:::jacobian_vech_Sigma(L_params, K_w = 2, rc_correlation = TRUE)
#' dim(J_mat)  # 3 x 3 for K_w=2 correlated
#' @keywords internal
jacobian_vech_Sigma <- function(L_params, K_w, rc_correlation = TRUE) {
    .Call(`_choicer_jacobian_vech_Sigma`, L_params, K_w, rc_correlation)
}

#' Analytical Hessian of the log-likelihood v2
#'
#' Computes the Hessian of the log-likelihood for the Mixed Logit model using
#' OpenMP for parallelization. Mirrors the parameters of mxl_loglik_gradient_parallel.
#'
#' @param theta vector collecting model parameters (beta, mu, L, delta (ASCs))
#' @param X design matrix for covariates with fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for covariates with random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with choice situation draws; K_w x S x N
#' @param rc_dist K_w x 1 integer vector indicating distribution of random coefficients: 0 = normal, 1 = log-normal
#' @param rc_correlation whether random coefficients should be correlated
#' @param rc_mean whether to estimate means for random coefficients.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns Hessian evaluated at input arguments
#' @note For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE,
#'   the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta),
#'   where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2).
#'   This differs from the textbook parameterization exp(mu_k + L_k * eta).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' theta <- rep(0, ncol(d$X) + ncol(d$W) + nrow(d$alt_mapping) - 1)
#' H <- choicer:::mxl_hessian_parallel(theta, d$X, d$W, d$alt_idx, d$choice_idx,
#'   d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
#'   rc_correlation = FALSE, rc_mean = FALSE)
#' dim(H)
#' }
#' @keywords internal
mxl_hessian_parallel <- function(theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_hessian_parallel`, theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' BHHH (outer product of gradients) information matrix for Mixed Logit
#'
#' Computes the BHHH approximation to the observed information matrix for the
#' Mixed Logit model: \eqn{H_{BHHH} = \sum_i w_i \cdot s_i s_i^\top}, where
#' \eqn{s_i} is the per-individual score (gradient of \eqn{\log \bar{P}_i}).
#' This outer product of gradients (OPG) estimator provides an alternative to
#' the analytical Hessian for standard error computation that scales to large
#' problems where the analytical Hessian is infeasible (e.g., many alternatives
#' or simulation draws).
#'
#' @param theta vector collecting model parameters (beta, mu, L, delta (ASCs))
#' @param X design matrix for covariates with fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for covariates with random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives within each choice set; 1-based indexing
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 1-based indexing relative to X; 0 is used if include_outside_option=True
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with choice situation draws; K_w x S x N
#' @param rc_dist K_w x 1 integer vector indicating distribution of random coefficients: 0 = normal, 1 = log-normal
#' @param rc_correlation whether random coefficients should be correlated
#' @param rc_mean whether to estimate means for random coefficients.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to 0 (if so, the outside option is not included in the data)
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns n_params x n_params PSD matrix representing the observed information
#'   matrix estimated by the outer product of gradients (same sign convention
#'   as the negated Hessian returned by \code{mxl_hessian_parallel}, so it can
#'   be inverted directly to obtain vcov).
#' @note The BHHH/OPG estimator is only asymptotically equivalent to the
#'   Hessian-based information matrix at the true MLE. In finite samples it can
#'   underestimate standard errors, particularly when the model is mis-specified
#'   or away from the optimum.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' theta <- rep(0, ncol(d$X) + ncol(d$W) + nrow(d$alt_mapping) - 1)
#' H <- choicer:::mxl_bhhh_parallel(theta, d$X, d$W, d$alt_idx, d$choice_idx,
#'   d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
#'   rc_correlation = FALSE, rc_mean = FALSE)
#' dim(H)
#' }
#' @keywords internal
mxl_bhhh_parallel <- function(theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_bhhh_parallel`, theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

mxl_scores_parallel <- function(theta, X, W, alt_idx, choice_idx, M, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_scores_parallel`, theta, X, W, alt_idx, choice_idx, M, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Per-observation simulated choice probabilities for Mixed Logit
#'
#' Returns the simulated choice probability for each (individual, alternative)
#' row of `X`, averaged over the supplied Halton draws. Mirrors `mnl_predict`.
#'
#' @param theta parameter vector (beta, \[mu\], L, delta)
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters are estimated
#' @param use_asc whether ASCs are included
#' @param include_outside_option whether the outside option is present
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns List with `choice_prob` (length sum(M)), `utility` (length sum(M),
#'   simulated mean of the deterministic + W*gamma component), and, when
#'   `include_outside_option = TRUE`, `choice_prob_outside` (length N).
#' @keywords internal
mxl_predict <- function(theta, X, W, alt_idx, M, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_predict`, theta, X, W, alt_idx, M, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Simulated expected logsum (inclusive value) for Mixed Logit
#'
#' Computes the simulated expected logsum (expected maximum utility, up to an
#' additive constant) for each choice situation:
#' \deqn{logsum_i = (1/S) \sum_s \log \sum_j \exp(V_{ij}^s),}
#' where the inner sum runs over individual i's alternatives and includes the
#' outside option's \eqn{\exp(0)} term when `include_outside_option = TRUE`.
#' The log-sum-exp must be averaged *across draws*: applying log-sum-exp to
#' the draw-averaged utilities returned by `mxl_predict` understates the
#' expectation because log-sum-exp is convex (Jensen's inequality).
#'
#' @param theta parameter vector (beta, \[mu\], L, delta)
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters are estimated
#' @param use_asc whether ASCs are included
#' @param include_outside_option whether the outside option is present
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns Vector of length N with the simulated expected logsum per choice
#'   situation.
#' @note For log-normal random coefficients (rc_dist=1) with rc_mean=TRUE,
#'   the distribution is a shifted log-normal: beta_k = exp(mu_k) + exp(L_k * eta),
#'   where exp(mu_k) shifts the location and exp(L_k * eta) ~ LogNormal(0, sigma_k^2).
#'   This differs from the textbook parameterization exp(mu_k + L_k * eta).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' fit <- run_mxlogit(input_data = d, eta_draws = eta)
#' ls <- choicer:::mxl_logsum(coef(fit), d$X, d$W, d$alt_idx, d$M, eta,
#'   rc_dist = rep(0L, ncol(d$W)), rc_correlation = FALSE, rc_mean = FALSE)
#' head(ls)
#' }
#' @keywords internal
mxl_logsum <- function(theta, X, W, alt_idx, M, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_logsum`, theta, X, W, alt_idx, M, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Predicted aggregate market shares for Mixed Logit
#'
#' Exported wrapper around the internal `mxl_predict_shares_internal`. Parses
#' `theta` using the standard parameter ordering and returns the simulated
#' weighted-average market shares.
#'
#' @param theta parameter vector (beta, \[mu\], L, delta)
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters are estimated
#' @param use_asc whether ASCs are included
#' @param include_outside_option whether outside option is included
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns Vector of length J (or J+1 with outside option) of predicted shares.
#' @keywords internal
mxl_predict_shares <- function(theta, X, W, alt_idx, M, weights, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_predict_shares`, theta, X, W, alt_idx, M, weights, eta_draws, rc_dist, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Diversion ratios for Mixed Logit (simulated, derivative-based)
#'
#' Computes the matrix of attribute-based diversion ratios for a fitted
#' Mixed Logit model. DR(k, j) is the fraction of demand lost by alternative
#' `j` that is captured by alternative `k` when a marginal change in
#' alternative j's `elast_var` attribute reduces s_j.
#'
#' In MNL the per-draw realized coefficient is a constant, so it cancels in
#' the ratio and the result is independent of the variable chosen. In MXL,
#' the realized coefficient \eqn{\beta_{ik}^s} varies across individuals
#' and draws, so the diversion ratio depends on which attribute is perturbed.
#' For a variable with a fixed coefficient the dependence again vanishes
#' (the constant cancels); for a random-coefficient variable it does not.
#'
#' @param theta parameter vector (beta, \[mu\], L, delta)
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param elast_var_idx 1-based index of the perturbed variable
#' @param is_random_coef TRUE if the variable is in W (random coef), FALSE if in X (fixed)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters are estimated
#' @param use_asc whether ASCs are included
#' @param include_outside_option whether outside option is included
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns J x J (or (J+1) x (J+1)) matrix of diversion ratios with zero diagonal.
#' @keywords internal
mxl_diversion_ratios_parallel <- function(theta, X, W, alt_idx, M, weights, eta_draws, rc_dist, elast_var_idx, is_random_coef, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_diversion_ratios_parallel`, theta, X, W, alt_idx, M, weights, eta_draws, rc_dist, elast_var_idx, is_random_coef, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' BLP contraction mapping for mixed logit
#'
#' Finds the ASC (delta) parameters such that predicted market shares
#' match target shares, using the contraction mapping of Berry, Levinsohn,
#' and Pakes (1995).
#'
#' @param delta J-1 or J vector with initial guess for deltas (ASCs)
#' @param target_shares J vector with target market shares
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param beta K_x vector with fixed coefficients
#' @param mu K_w vector with mean parameters (raw, will be transformed if log-normal)
#' @param L_params Cholesky parameters vector
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters represent means (TRUE) or are zero (FALSE)
#' @param include_outside_option whether outside option is included
#' @param tol convergence tolerance (default 1e-8)
#' @param max_iter maximum iterations (default 1000)
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations, \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns vector with converged delta (ASC) values
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' fit <- run_mxlogit(input_data = d, eta_draws = eta)
#' pm <- fit$param_map
#' delta <- mxl_blp_contraction(rep(0, J), rep(1/J, J), d$X, d$W,
#'   coef(fit)[pm$beta], rep(0, ncol(d$W)), coef(fit)[pm$sigma],
#'   d$alt_idx, d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
#'   rc_correlation = FALSE, rc_mean = FALSE)
#' delta
#' }
#' @export
mxl_blp_contraction <- function(delta, target_shares, X, W, beta, mu, L_params, alt_idx, M, weights, eta_draws, rc_dist, rc_correlation = TRUE, rc_mean = FALSE, include_outside_option = FALSE, tol = 1e-8, max_iter = 1000L, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_blp_contraction`, delta, target_shares, X, W, beta, mu, L_params, alt_idx, M, weights, eta_draws, rc_dist, rc_correlation, rc_mean, include_outside_option, tol, max_iter, gen_seed, gen_scramble, gen_S)
}

#' Compute aggregate elasticities for mixed logit model
#'
#' Computes the aggregate elasticity matrix (weighted average of individual
#' elasticities) for the Mixed Logit model. The elasticity E(i,j) represents
#' the percentage change in the probability of choosing alternative i when
#' the attribute of alternative j changes by 1%.
#'
#' @param theta parameter vector (beta, \[mu\], L, delta)
#' @param X design matrix for fixed coefficients; sum(M_i) x K_x
#' @param W design matrix for random coefficients; sum(M_i) x K_w or J x K_w
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing
#' @param choice_idx N x 1 vector (kept for API consistency, not used)
#' @param M N x 1 vector with number of alternatives for each individual
#' @param weights N x 1 vector with weights for each observation
#' @param eta_draws Array with draws; K_w x S x N
#' @param rc_dist K_w vector indicating distribution (0=normal, 1=log-normal)
#' @param elast_var_idx 1-based index of the variable for elasticity computation
#' @param is_random_coef TRUE if variable is in W (random coef), FALSE if in X (fixed coef)
#' @param rc_correlation whether random coefficients are correlated
#' @param rc_mean whether mu parameters are estimated
#' @param use_asc whether ASCs are included
#' @param include_outside_option whether outside option is included
#' @param gen_seed Integer master seed for the on-the-fly Halton generator. \code{< 0}
#'   (default) uses the materialized \code{eta_draws} cube; \code{>= 0} generates draws
#'   on the fly from this seed.
#' @param gen_scramble Integer scramble mode for on-the-fly generation: \code{0} =
#'   identity permutations (plain Halton, compat), \code{1} = seeded position-wise digit permutations.
#' @param gen_S Integer number of draws per individual, used only when \code{gen_seed >= 0}.
#' @returns J x J matrix of aggregate elasticities
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 3
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), w1 = rnorm(.N))]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_mxl_data(dt, "id", "alt", "choice", "x1", "w1")
#' eta <- get_halton_normals(50, d$N, ncol(d$W))
#' fit <- run_mxlogit(input_data = d, eta_draws = eta)
#' elas <- choicer:::mxl_elasticities_parallel(coef(fit), d$X, d$W, d$alt_idx,
#'   d$choice_idx, d$M, d$weights, eta, rc_dist = rep(0L, ncol(d$W)),
#'   elast_var_idx = 1L, is_random_coef = FALSE,
#'   rc_correlation = FALSE, rc_mean = FALSE)
#' elas
#' }
#' @keywords internal
mxl_elasticities_parallel <- function(theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, elast_var_idx, is_random_coef, rc_correlation = TRUE, rc_mean = FALSE, use_asc = TRUE, include_outside_option = FALSE, gen_seed = -1L, gen_scramble = 1L, gen_S = 0L) {
    .Call(`_choicer_mxl_elasticities_parallel`, theta, X, W, alt_idx, choice_idx, M, weights, eta_draws, rc_dist, elast_var_idx, is_random_coef, rc_correlation, rc_mean, use_asc, include_outside_option, gen_seed, gen_scramble, gen_S)
}

#' Log-likelihood and gradient for Nested Logit model
#'
#' Computes the log-likelihood and its gradient for the Nested Logit model using OpenMP for parallelization.
#' Especially handles singleton nests by fixing their lambda parameters to 1. Only non-singleton nests have a inclusive value coefficient estimated in theta.
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model parameters.
#'        Order: `[beta (K), lambda (n_non_singleton_nests), delta (n_delta)]`
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 0 for outside option,
#'        1-based index relative to rows in X_i otherwise.
#' @param nest_idx J x 1 vector with indices of nests for each alternative; 1-based indexing (1 to n_nests).
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns List with loglikelihood and gradient evaluated at input arguments
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
#' theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
#' result <- choicer:::nl_loglik_gradient_parallel(theta, d$X, d$alt_idx,
#'   d$choice_idx, d$nest_idx, d$M, d$weights)
#' result$objective
#' }
#' @keywords internal
nl_loglik_gradient_parallel <- function(theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_loglik_gradient_parallel`, theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc, include_outside_option)
}

#' BHHH/OPG information matrix for the Nested Logit model
#'
#' Computes the weighted outer product of per-individual scores
#' \eqn{\sum_i w_i\, s_i s_i^\top} for the Nested Logit model. The
#' per-individual score \eqn{s_i} (over the beta, lambda and delta/ASC blocks)
#' is the (positive) gradient of individual \eqn{i}'s log-likelihood
#' contribution and is weight-free; the supplied \code{weights} enter only as
#' the leading multiplier. Passing \code{weights = w} yields the ordinary
#' weighted BHHH/OPG information; passing \code{weights = w^2} yields the
#' sandwich \emph{meat} \eqn{B = \sum_i w_i^2 s_i s_i^\top} for robust (WESML)
#' inference. Singleton-nest lambdas are fixed to 1 and contribute no score
#' (mirroring the gradient kernel).
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model parameters.
#'        Order: `[beta (K), lambda (n_non_singleton_nests), delta (n_delta)]`
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 0 for outside option,
#'        1-based index relative to rows in X_i otherwise.
#' @param nest_idx J x 1 vector with indices of nests for each alternative; 1-based indexing (1 to n_nests).
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns A symmetric positive-semidefinite information matrix
#'   \eqn{\sum_i w_i\, s_i s_i^\top} (same sign convention as the negated Hessian).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
#' theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
#' B <- choicer:::nl_bhhh_parallel(theta, d$X, d$alt_idx, d$choice_idx,
#'   d$nest_idx, d$M, d$weights)
#' dim(B)
#' }
#' @keywords internal
nl_bhhh_parallel <- function(theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_bhhh_parallel`, theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc, include_outside_option)
}

nl_scores_parallel <- function(theta, X, alt_idx, choice_idx, nest_idx, M, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_scores_parallel`, theta, X, alt_idx, choice_idx, nest_idx, M, use_asc, include_outside_option)
}

#' Numerical Hessian of the log-likelihood via finite differences
#'
#' @param theta (K + n_delta + n_nests) vector with model parameters.
#'        Order: `[beta (K), delta (n_delta), lambda (n_nests)]`
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param choice_idx N x 1 vector with indices of chosen alternatives; 0 for outside option,
#'        1-based index relative to rows in X_i otherwise.
#' @param nest_idx J x 1 vector with indices of nests for each alternative; 1-based indexing (1 to n_nests).
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @param eps finite difference step size
#' @returns Hessian evaluated at input arguments
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' K_x <- ncol(d$X); K_l <- length(unique(d$nest_idx))
#' theta <- c(rep(0, K_x), rep(0.5, K_l), rep(0, J - 1))
#' H <- choicer:::nl_loglik_numeric_hessian(theta, d$X, d$alt_idx, d$choice_idx,
#'   d$nest_idx, d$M, d$weights)
#' dim(H)
#' }
#' @keywords internal
nl_loglik_numeric_hessian <- function(theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE, eps = 1e-6) {
    .Call(`_choicer_nl_loglik_numeric_hessian`, theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc, include_outside_option, eps)
}

#' Analytical Hessian of the negated log-likelihood for the Nested Logit model
#'
#' Computes the exact (analytical) Hessian of the negated log-likelihood for the
#' Nested Logit model using OpenMP parallelisation with thread-local accumulators.
#' Covers all parameter blocks: beta-beta, beta-lambda, beta-delta,
#' lambda-lambda, lambda-delta, and delta-delta. Singleton nests (lambda fixed
#' to 1, not estimated) contribute no rows or columns to the lambda blocks.
#'
#' @param theta (K + n_non_singleton_nests + n_delta) parameter vector.
#'   Order: \code{[beta (K), lambda (n_non_singleton_nests), delta (n_delta)]}.
#'   Same layout as \code{nl_loglik_gradient_parallel}.
#' @param X sum(M) x K design matrix of covariates.
#' @param alt_idx sum(M)-length integer vector of 1-based alternative indices.
#' @param choice_idx N-length integer vector of 1-based chosen alternative
#'   indices; 0 indicates the outside option was chosen.
#' @param nest_idx J-length integer vector of 1-based nest indices for each
#'   inside alternative.
#' @param M N-length integer vector of alternative-set sizes.
#' @param weights N-length numeric vector of individual weights.
#' @param use_asc Logical; whether alternative-specific constants are included.
#' @param include_outside_option Logical; whether an outside option (V=0) is
#'   present.
#' @returns A symmetric (P x P) matrix: the Hessian of the negated
#'   log-likelihood evaluated at \code{theta}. Structurally identical to the
#'   output of \code{nl_loglik_numeric_hessian}; suitable for
#'   \code{invert_hessian()}.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' d <- prepare_nl_data(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' K_x <- ncol(d$X)
#' K_l <- sum(table(d$nest_idx) > 1)   # number of non-singleton nests (= 2)
#' theta <- c(rep(0, K_x), rep(0.8, K_l), rep(0, J - 1))
#' H <- choicer:::nl_loglik_hessian_parallel(theta, d$X, d$alt_idx, d$choice_idx,
#'   d$nest_idx, d$M, d$weights)
#' dim(H)
#' }
#' @keywords internal
nl_loglik_hessian_parallel <- function(theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_loglik_hessian_parallel`, theta, X, alt_idx, choice_idx, nest_idx, M, weights, use_asc, include_outside_option)
}

#' Prediction of choice probabilities and utilities for the Nested Logit model
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model
#'        parameters. Order: `[beta (K), lambda (non-singleton), delta]`.
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param nest_idx J x 1 vector with nest indices for each alternative; 1-based indexing.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns List with `choice_prob` (joint P_ij per stacked row) and `utility` (V_ij).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' pred <- choicer:::nl_predict(coef(fit), fit$data$X, fit$data$alt_idx, fit$data$M,
#'   fit$data$nest_idx, use_asc = TRUE)
#' head(pred$choice_prob)
#' }
#' @keywords internal
nl_predict <- function(theta, X, alt_idx, M, nest_idx, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_predict`, theta, X, alt_idx, M, nest_idx, use_asc, include_outside_option)
}

#' Prediction of market shares for the Nested Logit model
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model
#'        parameters. Order: `[beta (K), lambda (non-singleton), delta]`.
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param nest_idx J x 1 vector with nest indices for each alternative; 1-based indexing.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns vector with predicted market shares (outside-option share first when present).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' shares <- choicer:::nl_predict_shares(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$M, fit$data$weights, fit$data$nest_idx, use_asc = TRUE)
#' shares
#' }
#' @keywords internal
nl_predict_shares <- function(theta, X, alt_idx, M, weights, nest_idx, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_predict_shares`, theta, X, alt_idx, M, weights, nest_idx, use_asc, include_outside_option)
}

#' Compute aggregate elasticities for the Nested Logit model
#'
#' Computes the aggregate (weighted-average) elasticity matrix for the Nested
#' Logit model. Reduces to the MNL elasticities when all lambda = 1.
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model
#'        parameters. Order: `[beta (K), lambda (non-singleton), delta]`.
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param choice_idx N x 1 vector (kept for API consistency, not used).
#' @param nest_idx J x 1 vector with nest indices for each alternative; 1-based indexing.
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param elast_var_idx 1-based index of the column in X for which to compute the elasticity.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns J x J matrix of aggregate elasticities (row = responding alt, col = perturbed alt).
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' elas <- choicer:::nl_elasticities_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$choice_idx, fit$data$nest_idx, fit$data$M, fit$data$weights,
#'   elast_var_idx = 1L)
#' elas
#' }
#' @keywords internal
nl_elasticities_parallel <- function(theta, X, alt_idx, choice_idx, nest_idx, M, weights, elast_var_idx, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_elasticities_parallel`, theta, X, alt_idx, choice_idx, nest_idx, M, weights, elast_var_idx, use_asc, include_outside_option)
}

#' Compute Nested Logit diversion ratios (parallelized over individuals)
#'
#' Computes the diversion ratio matrix DR(j->k) for the Nested Logit model.
#' Entry (k, j) = fraction of demand lost by alternative j captured by k.
#' Reduces to the MNL diversion ratios when all lambda = 1.
#'
#' @param theta (K + n_non_singleton_nests + n_delta) vector with model
#'        parameters. Order: `[beta (K), lambda (non-singleton), delta]`.
#' @param X sum(M) x K design matrix with covariates.
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param nest_idx J x 1 vector with nest indices for each alternative; 1-based indexing.
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param use_asc whether to use alternative-specific constants.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @returns J x J matrix where entry (k, j) = DR(j->k). Diagonal is 0.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' dr <- choicer:::nl_diversion_ratios_parallel(coef(fit), fit$data$X, fit$data$alt_idx,
#'   fit$data$nest_idx, fit$data$M, fit$data$weights)
#' dr
#' }
#' @keywords internal
nl_diversion_ratios_parallel <- function(theta, X, alt_idx, nest_idx, M, weights, use_asc = TRUE, include_outside_option = FALSE) {
    .Call(`_choicer_nl_diversion_ratios_parallel`, theta, X, alt_idx, nest_idx, M, weights, use_asc, include_outside_option)
}

#' BLP95 contraction mapping for the Nested Logit model
#'
#' Damped iterative fixed point recovering delta given target shares, using the
#' NL probability structure. `damping = 1` reproduces the plain BLP update.
#'
#' @param delta J x 1 vector with initial guess for deltas (ASCs).
#' @param target_shares vector with target shares (outside-option share first when present).
#' @param X sum(M) x K design matrix with covariates.
#' @param beta K x 1 vector with fixed coefficients.
#' @param lambda full nest dissimilarity vector of length n_nests (singletons = 1).
#' @param alt_idx sum(M) x 1 vector with indices of alternatives; 1-based indexing.
#' @param nest_idx J x 1 vector with nest indices for each alternative; 1-based indexing.
#' @param M N x 1 vector with number of alternatives for each individual.
#' @param weights N x 1 vector with weights for each observation.
#' @param include_outside_option whether to include outside option normalized to V=0, lambda=1.
#' @param damping damping factor for the update (default 1.0 = plain BLP).
#' @param tol convergence tolerance.
#' @param max_iter maximum number of iterations.
#' @returns vector with contraction's delta (ASCs) output.
#' @examples
#' \donttest{
#' library(data.table)
#' set.seed(42)
#' N <- 50; J <- 4
#' dt <- data.table(id = rep(1:N, each = J), alt = rep(1:J, N))
#' dt[, `:=`(x1 = rnorm(.N), x2 = rnorm(.N))]
#' dt[, nest := ifelse(alt <= 2, "A", "B")]
#' dt[, choice := 0L]
#' dt[, choice := sample(c(1L, rep(0L, J - 1))), by = id]
#' fit <- run_nestlogit(dt, "id", "alt", "choice", c("x1", "x2"), "nest")
#' beta <- coef(fit)[fit$param_map$beta]
#' lambda <- rep(1, length(unique(fit$data$nest_idx)))
#' lambda[as.integer(names(which(table(fit$data$nest_idx) > 1)))] <-
#'   coef(fit)[fit$param_map$lambda]
#' delta <- nl_blp_contraction(rep(0, J), rep(1/J, J), fit$data$X, beta, lambda,
#'   fit$data$alt_idx, fit$data$nest_idx, fit$data$M, fit$data$weights)
#' delta
#' }
#' @export
nl_blp_contraction <- function(delta, target_shares, X, beta, lambda, alt_idx, nest_idx, M, weights, include_outside_option = FALSE, damping = 1.0, tol = 1e-8, max_iter = 1000L) {
    .Call(`_choicer_nl_blp_contraction`, delta, target_shares, X, beta, lambda, alt_idx, nest_idx, M, weights, include_outside_option, damping, tol, max_iter)
}

#' Query choicer OpenMP thread settings
#'
#' @return A list with OpenMP availability, active/max thread settings, CPU
#'   thread capacity reported by OpenMP, thread limits, and relevant
#'   environment variables.
#' @export
thread_info <- function() {
    .Call(`_choicer_thread_info`)
}

get_num_threads <- function() {
    invisible(.Call(`_choicer_get_num_threads`))
}

#' Set the number of OpenMP threads used by choicer
#'
#' @param n_threads Positive integer number of threads.
#' @return Invisibly returns `NULL`.
#' @export
set_num_threads <- function(n_threads) {
    invisible(.Call(`_choicer_set_num_threads`, n_threads))
}