R/RcppExports.R

In hmcdm: Hidden Markov Cognitive Diagnosis Models for Learning

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

bijectionvector <- function(K) {
    .Call(`_hmcdm_bijectionvector`, K)
}

#' @title Convert integer to attribute pattern
#' @description Based on the bijective relationship between natural numbers and sum of powers of two,
#'  convert integer between 0 and 2^K-1 to K-dimensional attribute pattern.
#' @param K An \code{int} for the number of attributes
#' @param CL An \code{int} between 0 and 2^K-1
#' @return A \code{vec} of the K-dimensional attribute pattern corresponding to CL.
#' @examples
#' inv_bijectionvector(4,0)
#' @export
inv_bijectionvector <- function(K, CL) {
    .Call(`_hmcdm_inv_bijectionvector`, K, CL)
}

rwishart <- function(df, S) {
    .Call(`_hmcdm_rwishart`, df, S)
}

rinvwish <- function(df, Sig) {
    .Call(`_hmcdm_rinvwish`, df, Sig)
}

rmultinomial <- function(ps) {
    .Call(`_hmcdm_rmultinomial`, ps)
}

rDirichlet <- function(deltas) {
    .Call(`_hmcdm_rDirichlet`, deltas)
}

dmvnrm <- function(x, mean, sigma, logd = FALSE) {
    .Call(`_hmcdm_dmvnrm`, x, mean, sigma, logd)
}

rmvnrm <- function(mu, sigma) {
    .Call(`_hmcdm_rmvnrm`, mu, sigma)
}

#' @title Generate random Q matrix
#' @description Creates a random Q matrix containing three identity matrices after row permutation
#' @param J An \code{int} that represents the number of items
#' @param K An \code{int} that represents the number of attributes/skills
#' @return A dichotomous \code{matrix} for Q.
#' @examples 
#' random_Q(15,4)
#' @export
random_Q <- function(J, K) {
    .Call(`_hmcdm_random_Q`, J, K)
}

#' @title Generate ideal response matrix
#' @description Based on the Q matrix and the latent attribute space, generate the ideal response matrix for each skill pattern
#' @param K An \code{int} of the number of attributes
#' @param J An \code{int} of the number of items
#' @param Q A J-by-K Q \code{matrix}
#' @return A J-by-2^K ideal response \code{matrix}
#' @examples 
#' Q = random_Q(15,4)
#' ETA = ETAmat(4,15,Q)
#' @export
ETAmat <- function(K, J, Q) {
    .Call(`_hmcdm_ETAmat`, K, J, Q)
}

#' @title Generate monotonicity matrix
#' @description Based on the latent attribute space, generate a matrix indicating whether it is possible to
#' transition from pattern cc to cc' under the monotonicity learning assumption.
#' @param K An \code{int} of the number of attributes.
#' @return A 2^K-by-2^K dichotomous \code{matrix} of whether it is possible to transition between two patterns 
#' @examples
#' TP = TPmat(4)
#' @export
TPmat <- function(K) {
    .Call(`_hmcdm_TPmat`, K)
}

crosstab <- function(V1, V2, TP, nClass, col_dim) {
    .Call(`_hmcdm_crosstab`, V1, V2, TP, nClass, col_dim)
}

resp_miss <- function(Responses, Test_order, Test_versions) {
    .Call(`_hmcdm_resp_miss`, Responses, Test_order, Test_versions)
}

#' @title Compute item pairwise odds ratio
#' @description Based on a response matrix, calculate the item pairwise odds-ratio according do (n11*n00)/(n10*n01), where nij is the
#' number of people answering both item i and item j correctly
#' @param N An \code{int} of the sample size
#' @param J An \code{int} of the number of items
#' @param Yt An N-by-J response \code{matrix}
#' @return A J-by-J upper-triangular \code{matrix} of the item pairwise odds ratios
#' @examples 
#' \donttest{
#' N = dim(Y_real_array)[1]
#' J = nrow(Q_matrix)
#' OddsRatio(N,J,Y_real_array[,,1])}
#' @export
OddsRatio <- function(N, J, Yt) {
    .Call(`_hmcdm_OddsRatio`, N, J, Yt)
}

getMode <- function(sorted_vec, size) {
    .Call(`_hmcdm_getMode`, sorted_vec, size)
}

Sparse2Dense <- function(Y_real_array, Test_order, Test_versions) {
    .Call(`_hmcdm_Sparse2Dense`, Y_real_array, Test_order, Test_versions)
}

Dense2Sparse <- function(Y_sim, Test_order, Test_versions) {
    .Call(`_hmcdm_Dense2Sparse`, Y_sim, Test_order, Test_versions)
}

Mat2Array <- function(Q_matrix, T) {
    .Call(`_hmcdm_Mat2Array`, Q_matrix, T)
}

Array2Mat <- function(r_stars) {
    .Call(`_hmcdm_Array2Mat`, r_stars)
}

#' @title Generate a list of Q-matrices for each examinee.
#' @description Generate a list of length N. Each element of the list is a JxK Q_matrix of all items
#' administered across all time points to the examinee, in the order of administration.
#' @param Q_matrix A J-by-K matrix, indicating the item-skill relationship.
#' @param Design_array An N-by-J-by-L array indicating whether examinee n has taken item j at l time point.
#' @return A list length of N. Each element of the list is a JxK Q_matrix for each examinee.
#' @examples 
#' \donttest{
#' Q_examinee = Q_list_g(Q_matrix, Design_array)}
#' @export
Q_list_g <- function(Q_matrix, Design_array) {
    .Call(`_hmcdm_Q_list_g`, Q_matrix, Design_array)
}

design_array <- function(Test_order, Test_versions, Jt) {
    .Call(`_hmcdm_design_array`, Test_order, Test_versions, Jt)
}

point_estimates_learning <- function(output, model, N, K, T, alpha_EAP = TRUE) {
    .Call(`_hmcdm_point_estimates_learning`, output, model, N, K, T, alpha_EAP)
}

Learning_fit_g <- function(output, model, Y_real_array, Q_matrix, Design_array, Q_examinee = NULL, Latency_array = NULL, G_version = NA_integer_, R = NULL) {
    .Call(`_hmcdm_Learning_fit_g`, output, model, Y_real_array, Q_matrix, Design_array, Q_examinee, Latency_array, G_version, R)
}

parm_update_HO_g <- function(Design_array, alphas, pi, lambdas, thetas, response, itempars, Q_examinee, theta_propose, deltas_propose, Q_matrix) {
    .Call(`_hmcdm_parm_update_HO_g`, Design_array, alphas, pi, lambdas, thetas, response, itempars, Q_examinee, theta_propose, deltas_propose, Q_matrix)
}

Gibbs_DINA_HO_g <- function(Response, Q_matrix, Design_array, theta_propose, deltas_propose, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_DINA_HO_g`, Response, Q_matrix, Design_array, theta_propose, deltas_propose, chain_length, burn_in)
}

parm_update_HO_RT_sep_g <- function(Design_array, alphas, pi, lambdas, thetas, latency, RT_itempars, taus, phi_vec, tauvar, response, itempars, Q_matrix, Q_examinee, G_version, theta_propose, a_sigma_tau0, rate_sigma_tau0, deltas_propose, a_alpha0, rate_alpha0) {
    .Call(`_hmcdm_parm_update_HO_RT_sep_g`, Design_array, alphas, pi, lambdas, thetas, latency, RT_itempars, taus, phi_vec, tauvar, response, itempars, Q_matrix, Q_examinee, G_version, theta_propose, a_sigma_tau0, rate_sigma_tau0, deltas_propose, a_alpha0, rate_alpha0)
}

Gibbs_DINA_HO_RT_sep_g <- function(Response, Latency, Q_matrix, Design_array, G_version, theta_propose, deltas_propose, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_DINA_HO_RT_sep_g`, Response, Latency, Q_matrix, Design_array, G_version, theta_propose, deltas_propose, chain_length, burn_in)
}

parm_update_HO_RT_joint_g <- function(Design_array, alphas, pi, lambdas, thetas, latency, RT_itempars, taus, phi_vec, Sig, response, itempars, Q_matrix, Q_examinee, G_version, sig_theta_propose, S, p, deltas_propose, a_alpha0, rate_alpha0) {
    .Call(`_hmcdm_parm_update_HO_RT_joint_g`, Design_array, alphas, pi, lambdas, thetas, latency, RT_itempars, taus, phi_vec, Sig, response, itempars, Q_matrix, Q_examinee, G_version, sig_theta_propose, S, p, deltas_propose, a_alpha0, rate_alpha0)
}

Gibbs_DINA_HO_RT_joint_g <- function(Response, Latency, Q_matrix, Design_array, G_version, sig_theta_propose, deltas_propose, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_DINA_HO_RT_joint_g`, Response, Latency, Q_matrix, Design_array, G_version, sig_theta_propose, deltas_propose, chain_length, burn_in)
}

parm_update_rRUM_g <- function(Design_array, alphas, pi, taus, R, r_stars, pi_stars, Q_matrix, responses, X_ijk, Smats, Gmats, dirich_prior) {
    invisible(.Call(`_hmcdm_parm_update_rRUM_g`, Design_array, alphas, pi, taus, R, r_stars, pi_stars, Q_matrix, responses, X_ijk, Smats, Gmats, dirich_prior))
}

Gibbs_rRUM_indept_g <- function(Response, Q_matrix, R, Design_array, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_rRUM_indept_g`, Response, Q_matrix, R, Design_array, chain_length, burn_in)
}

parm_update_NIDA_indept_g <- function(Design_array, alphas, pi, taus, R, Q_matrix, responses, X_ijk, Smats, Gmats, dirich_prior) {
    invisible(.Call(`_hmcdm_parm_update_NIDA_indept_g`, Design_array, alphas, pi, taus, R, Q_matrix, responses, X_ijk, Smats, Gmats, dirich_prior))
}

Gibbs_NIDA_indept_g <- function(Response, Q_matrix, R, Design_array, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_NIDA_indept_g`, Response, Q_matrix, R, Design_array, chain_length, burn_in)
}

parm_update_DINA_FOHM <- function(N, J, K, nClass, nT, Y, TP, ETA, ss, gs, CLASS, pi, Omega) {
    invisible(.Call(`_hmcdm_parm_update_DINA_FOHM`, N, J, K, nClass, nT, Y, TP, ETA, ss, gs, CLASS, pi, Omega))
}

Gibbs_DINA_FOHM_g <- function(Response, Q_matrix, Design_array, chain_length, burn_in) {
    .Call(`_hmcdm_Gibbs_DINA_FOHM_g`, Response, Q_matrix, Design_array, chain_length, burn_in)
}

#' @title Gibbs sampler for learning models
#' @description Runs MCMC to estimate parameters of any of the listed learning models. 
#' @param Response An \code{array} of dichotomous item responses. t-th slice is an N-by-J matrix of responses at time t.
#' @param Q_matrix A J-by-K Q-matrix. 
#' @param model A \code{charactor} of the type of model fitted with the MCMC sampler, possible selections are 
#' "DINA_HO": Higher-Order Hidden Markov Diagnostic Classification Model with DINA responses;
#' "DINA_HO_RT_joint": Higher-Order Hidden Markov DCM with DINA responses, log-Normal response times, and joint modeling of latent
#' speed and learning ability; 
#' "DINA_HO_RT_sep": Higher-Order Hidden Markov DCM with DINA responses, log-Normal response times, and separate modeling of latent
#' speed and learning ability; 
#' "rRUM_indept": Simple independent transition probability model with rRUM responses
#' "NIDA_indept": Simple independent transition probability model with NIDA responses
#' "DINA_FOHM": First Order Hidden Markov model with DINA responses
#' @param Design_array An \code{array} of dimension N-by-J-by-L indicating the items assigned (1/0) to each subject at each time point. 
#' Either 'Design_array' or both 'Test_order' & 'Test_versions' need to be provided to run HMCDM. 
#' @param Test_order Optional. A \code{matrix} of the order of item blocks for each test version.
#' @param Test_versions Optional. A \code{vector} of the test version of each learner.
#' @param chain_length An \code{int} of the MCMC chain length.
#' @param burn_in An \code{int} of the MCMC burn-in chain length.
#' @param Latency_array Optional. A \code{array} of the response times. t-th slice is an N-by-J matrix of response times at time t.
#' @param G_version Optional. An \code{int} of the type of covariate for increased fluency (1: G is dichotomous depending on whether all skills required for
#' current item are mastered; 2: G cumulates practice effect on previous items using mastered skills; 3: G is a time block effect invariant across 
#' subjects with different attribute trajectories)
#' @param theta_propose Optional. A \code{scalar} for the standard deviation of theta's proposal distribution in the MH sampling step.
#' @param deltas_propose Optional. A \code{vector} for the band widths of each lambda's proposal distribution in the MH sampling step.
#' @param R Optional. A reachability \code{matrix} for the hierarchical relationship between attributes. 
#' @return A \code{list} of parameter samples and Metropolis-Hastings acceptance rates (if applicable).
#' @author Susu Zhang
#' @examples
#' \donttest{
#' output_FOHM = hmcdm(Y_real_array, Q_matrix, "DINA_FOHM", Design_array, 100, 30)
#' }
#' @export
hmcdm <- function(Response, Q_matrix, model, Design_array = NULL, Test_order = NULL, Test_versions = NULL, chain_length = 100L, burn_in = 50L, G_version = NA_integer_, theta_propose = 0., Latency_array = NULL, deltas_propose = NULL, R = NULL) {
    .Call(`_hmcdm_hmcdm`, Response, Q_matrix, model, Design_array, Test_order, Test_versions, chain_length, burn_in, G_version, theta_propose, Latency_array, deltas_propose, R)
}

sim_resp_DINA <- function(J, K, ETA, Svec, Gvec, alpha) {
    .Call(`_hmcdm_sim_resp_DINA`, J, K, ETA, Svec, Gvec, alpha)
}

simDINA_g <- function(alphas, itempars, Q_matrix, Design_array) {
    .Call(`_hmcdm_simDINA_g`, alphas, itempars, Q_matrix, Design_array)
}

pYit_DINA <- function(ETA_it, Y_it, itempars) {
    .Call(`_hmcdm_pYit_DINA`, ETA_it, Y_it, itempars)
}

sim_resp_rRUM <- function(J, K, Q, rstar, pistar, alpha) {
    .Call(`_hmcdm_sim_resp_rRUM`, J, K, Q, rstar, pistar, alpha)
}

simrRUM_g <- function(alphas, r_stars_mat, pi_stars, Q_matrix, Design_array) {
    .Call(`_hmcdm_simrRUM_g`, alphas, r_stars_mat, pi_stars, Q_matrix, Design_array)
}

pYit_rRUM <- function(alpha_it, Y_it, pi_star_it, r_star_it, Q_it) {
    .Call(`_hmcdm_pYit_rRUM`, alpha_it, Y_it, pi_star_it, r_star_it, Q_it)
}

sim_resp_NIDA <- function(J, K, Q, Svec, Gvec, alpha) {
    .Call(`_hmcdm_sim_resp_NIDA`, J, K, Q, Svec, Gvec, alpha)
}

simNIDA_g <- function(alphas, Svec, Gvec, Q_matrix, Design_array) {
    .Call(`_hmcdm_simNIDA_g`, alphas, Svec, Gvec, Q_matrix, Design_array)
}

pYit_NIDA <- function(alpha_it, Y_it, Svec, Gvec, Q_it) {
    .Call(`_hmcdm_pYit_NIDA`, alpha_it, Y_it, Svec, Gvec, Q_it)
}

#' @title Simulate responses from the specified model (entire cube)
#' @description Simulate a cube of responses from the specified model for all persons on items across all time points. 
#' Currently available models are `DINA`, `rRUM`, and `NIDA`.
#' @param model The cognitive diagnostic model under which the item responses are generated
#' @param alphas An N-by-K-by-L \code{array} of attribute patterns of all persons across L time points 
#' @param Q_matrix A J-by-K of Q-matrix
#' @param Design_array A N-by-J-by-L array indicating whether item j is administered to examinee i at l time point.
#' @param itempars A J-by-2 \code{mat} of item parameters (slipping: 1st col, guessing: 2nd col).
#' @param r_stars A J-by-K \code{mat} of item penalty parameters for missing skills.
#' @param pi_stars A length J \code{vector} of item correct response probability with all requisite skills.
#' @param Svec A length K \code{vector} of slipping probability in applying mastered skills
#' @param Gvec A length K \code{vector} of guessing probability in applying mastered skills
#' @return An \code{array} of item responses from the specified model of examinees across all time points.
#' @examples
#' \donttest{
#' ## DINA ##
#' N = nrow(Design_array)
#' J = nrow(Q_matrix)
#' thetas_true = rnorm(N, 0, 1.8)
#' lambdas_true <- c(-2, .4, .055)
#' Alphas <- sim_alphas(model="HO_joint", 
#'                     lambdas=lambdas_true, 
#'                     thetas=thetas_true, 
#'                     Q_matrix=Q_matrix, 
#'                     Design_array=Design_array)
#' itempars_true <- matrix(runif(J*2,.1,.2), ncol=2)
#' 
#' Y_sim <- sim_hmcdm(model="DINA",Alphas,Q_matrix,Design_array,
#'                    itempars=itempars_true)
#'                    
#' ## rRUM ##
#' J = nrow(Q_matrix)
#' K = ncol(Q_matrix)
#' Smats <- matrix(runif(J*K,.1,.3),c(J,K))
#' Gmats <- matrix(runif(J*K,.1,.3),c(J,K))
#' r_stars <- Gmats / (1-Smats)
#' pi_stars <- apply((1-Smats)^Q_matrix, 1, prod)
#' 
#' Y_sim <- sim_hmcdm(model="rRUM",Alphas,Q_matrix,Design_array,
#'                    r_stars=r_stars,pi_stars=pi_stars)
#' 
#' ## NIDA ##
#' K = ncol(Q_matrix)
#' Svec <- runif(K,.1,.3)
#' Gvec <- runif(K,.1,.3)
#' 
#' Y_sim <- sim_hmcdm(model="NIDA",Alphas,Q_matrix,Design_array,
#'                    Svec=Svec,Gvec=Gvec)
#' }
#' @export
sim_hmcdm <- function(model, alphas, Q_matrix, Design_array, itempars = NULL, r_stars = NULL, pi_stars = NULL, Svec = NULL, Gvec = NULL) {
    .Call(`_hmcdm_sim_hmcdm`, model, alphas, Q_matrix, Design_array, itempars, r_stars, pi_stars, Svec, Gvec)
}

J_incidence_cube_g <- function(Q_matrix, Design_array) {
    .Call(`_hmcdm_J_incidence_cube_g`, Q_matrix, Design_array)
}

G2vec_efficient_g <- function(ETA, J_incidence, alphas_i, t, Q_matrix, Design_array, i) {
    .Call(`_hmcdm_G2vec_efficient_g`, ETA, J_incidence, alphas_i, t, Q_matrix, Design_array, i)
}

#' @title Simulate item response times based on Wang et al.'s (2018) joint model of response times and accuracy in learning
#' @description Simulate a cube of subjects' response times across time points according to a variant of the logNormal model
#' @param alphas An N-by-K-by-T \code{array} of attribute patterns of all persons across T time points 
#' @param Q_matrix A J-by-K  Q-matrix for the test
#' @param Design_array A N-by-J-by-L array indicating whether item j is administered to examinee i at l time point.
#' @param RT_itempars A J-by-2 \code{matrix} of item time discrimination and time intensity parameters
#' @param taus A length N \code{vector} of latent speed of each person
#' @param phi A \code{scalar} of slope of increase in fluency over time due to covariates (G)
#' @param G_version An \code{int} of the type of covariate for increased fluency (1: G is dichotomous depending on whether all skills required for
#' current item are mastered; 2: G cumulates practice effect on previous items using mastered skills; 3: G is a time block effect invariant across 
#' subjects with different attribute trajectories)
#' @return A \code{cube} of response times of subjects on each item across time
#' @examples
#' N = dim(Design_array)[1]
#' J = nrow(Q_matrix)
#' K = ncol(Q_matrix)
#' L = dim(Design_array)[3]
#' class_0 <- sample(1:2^K, N, replace = TRUE)
#' Alphas_0 <- matrix(0,N,K)
#' mu_thetatau = c(0,0)
#' Sig_thetatau = rbind(c(1.8^2,.4*.5*1.8),c(.4*.5*1.8,.25))
#' Z = matrix(rnorm(N*2),N,2)
#' thetatau_true = Z%*%chol(Sig_thetatau)
#' thetas_true = thetatau_true[,1]
#' taus_true = thetatau_true[,2]
#' G_version = 3
#' phi_true = 0.8
#' for(i in 1:N){
#'   Alphas_0[i,] <- inv_bijectionvector(K,(class_0[i]-1))
#' }
#' lambdas_true <- c(-2, .4, .055)     
#' Alphas <- sim_alphas(model="HO_joint", 
#'                      lambdas=lambdas_true, 
#'                      thetas=thetas_true, 
#'                      Q_matrix=Q_matrix, 
#'                      Design_array=Design_array)
#' RT_itempars_true <- matrix(NA, nrow=J, ncol=2)
#' RT_itempars_true[,2] <- rnorm(J,3.45,.5)
#' RT_itempars_true[,1] <- runif(J,1.5,2)
#' ETAs <- ETAmat(K,J,Q_matrix)
#' L_sim <- sim_RT(Alphas,Q_matrix,Design_array,RT_itempars_true,taus_true,phi_true,G_version)
#' @export
sim_RT <- function(alphas, Q_matrix, Design_array, RT_itempars, taus, phi, G_version) {
    .Call(`_hmcdm_sim_RT`, alphas, Q_matrix, Design_array, RT_itempars, taus, phi, G_version)
}

dLit <- function(G_it, L_it, RT_itempars_it, tau_i, phi) {
    .Call(`_hmcdm_dLit`, G_it, L_it, RT_itempars_it, tau_i, phi)
}

simulate_alphas_HO_sep_g <- function(lambdas, thetas, Q_matrix, Design_array, alpha0) {
    .Call(`_hmcdm_simulate_alphas_HO_sep_g`, lambdas, thetas, Q_matrix, Design_array, alpha0)
}

pTran_HO_sep_g <- function(alpha_prev, alpha_post, lambdas, theta_i, Q_i, Design_array, t, i) {
    .Call(`_hmcdm_pTran_HO_sep_g`, alpha_prev, alpha_post, lambdas, theta_i, Q_i, Design_array, t, i)
}

simulate_alphas_HO_joint_g <- function(lambdas, thetas, Q_matrix, Design_array, alpha0) {
    .Call(`_hmcdm_simulate_alphas_HO_joint_g`, lambdas, thetas, Q_matrix, Design_array, alpha0)
}

pTran_HO_joint_g <- function(alpha_prev, alpha_post, lambdas, theta_i, Q_i, Design_array, t, i) {
    .Call(`_hmcdm_pTran_HO_joint_g`, alpha_prev, alpha_post, lambdas, theta_i, Q_i, Design_array, t, i)
}

simulate_alphas_indept <- function(taus, alpha0s, L, R) {
    .Call(`_hmcdm_simulate_alphas_indept`, taus, alpha0s, L, R)
}

simulate_alphas_indept_g <- function(taus, N, L, R, alpha0) {
    .Call(`_hmcdm_simulate_alphas_indept_g`, taus, N, L, R, alpha0)
}

pTran_indept <- function(alpha_prev, alpha_post, taus, R) {
    .Call(`_hmcdm_pTran_indept`, alpha_prev, alpha_post, taus, R)
}

simulate_alphas_FOHM <- function(Omega, N, L, alpha0) {
    .Call(`_hmcdm_simulate_alphas_FOHM`, Omega, N, L, alpha0)
}

rAlpha <- function(Omega, N, L, alpha1) {
    .Call(`_hmcdm_rAlpha`, Omega, N, L, alpha1)
}

#' @title Generate a random transition matrix for the first order hidden Markov model
#' @description Generate a random transition matrix under nondecreasing learning trajectory assumption
#' @param TP A 2^K-by-2^K dichotomous matrix of indicating possible transitions under the monotonicity assumption, created with
#' the TPmat function 
#' @return A 2^K-by-2^K transition matrix, the (i,j)th element indicating the transition probability of transitioning from i-th class to j-th class.
#' @examples
#' K = ncol(Q_matrix)
#' TP = TPmat(K)
#' Omega_sim = rOmega(TP)
#' @export
rOmega <- function(TP) {
    .Call(`_hmcdm_rOmega`, TP)
}

#' @title Generate attribute trajectories under the specified hidden Markov models
#' @description Based on the learning model parameters, create cube of attribute patterns
#' of all subjects across time. 
#' Currently available learning models are Higher-order hidden Markov DCM('HO_sep'), 
#' Higher-order hidden Markov DCM with learning ability as a random effect('HO_joint'), 
#' the simple independent-attribute learning model('indept'), 
#' and the first order hidden Markov model('FOHM'). 
#' @param model The learning model under which the attribute trajectories are generated. Available options are: 'HO_joint', 'HO_sep', 'indept', 'FOHM'.
#' @param lambdas A \code{vector} of transition model coefficients. With 'HO_sep' model specification, `lambdas` should be a length 4 \code{vector}. First entry is intercept of the logistic transition
#' model, second entry is the slope of general learning ability, third entry is the slope for number of other mastered skills,
#' fourth entry is the slope for amount of practice.
#' With 'HO_joint' model specification, `lambdas` should be a length 3 \code{vector}. First entry is intercept of the logistic transition
#' model, second entry is the slope for number of other mastered skills, third entry is the slope for amount of practice.
#' @param thetas A length N \code{vector} of learning abilities of each subject.
#' @param Q_matrix A J-by-K Q-matrix
#' @param Design_array A N-by-J-by-L array indicating items administered to examinee n at time point l.
#' @param taus A length K \code{vector} of transition probabilities from 0 to 1 on each skill
#' @param Omega A 2^K-by-2^K \code{matrix} of transition probabilities from row pattern to column pattern
#' @param N An \code{int} of number of examinees.
#' @param L An \code{int} of number of time points.
#' @param R A K-by-K dichotomous reachability \code{matrix} indicating the attribute hierarchies. The k,k'th entry of R is 1 if k' is prereq to k.
#' @param alpha0 Optional. An N-by-K \code{matrix} of subjects' initial attribute patterns.
#' @return An N-by-K-by-L \code{array} of attribute patterns of subjects at each time point.
#' @examples
#' \donttest{
#' ## HO_joint ##
#' N = nrow(Design_array)
#' J = nrow(Q_matrix)
#' K = ncol(Q_matrix)
#' L = dim(Design_array)[3]
#' class_0 <- sample(1:2^K, N, replace = TRUE)
#' Alphas_0 <- matrix(0,N,K)
#' for(i in 1:N){
#'   Alphas_0[i,] <- inv_bijectionvector(K,(class_0[i]-1))
#' }
#' thetas_true = rnorm(N, 0, 1.8)
#' lambdas_true <- c(-2, .4, .055)
#' Alphas <- sim_alphas(model="HO_joint", 
#'                     lambdas=lambdas_true, 
#'                     thetas=thetas_true, 
#'                     Q_matrix=Q_matrix, 
#'                     Design_array=Design_array)
#' 
#' ## HO_sep ##
#' N = dim(Design_array)[1]
#' J = nrow(Q_matrix)
#' K = ncol(Q_matrix)
#' L = dim(Design_array)[3]
#' class_0 <- sample(1:2^K, N, replace = L)
#' Alphas_0 <- matrix(0,N,K)
#' for(i in 1:N){
#'   Alphas_0[i,] <- inv_bijectionvector(K,(class_0[i]-1))
#' }
#' thetas_true = rnorm(N)
#' lambdas_true = c(-1, 1.8, .277, .055)
#' Alphas <- sim_alphas(model="HO_sep", 
#'                      lambdas=lambdas_true, 
#'                      thetas=thetas_true, 
#'                      Q_matrix=Q_matrix, 
#'                      Design_array=Design_array)
#' 
#' ## indept ##
#' N = dim(Design_array)[1]
#' K = dim(Q_matrix)[2]
#' L = dim(Design_array)[3]
#' tau <- numeric(K)
#' for(k in 1:K){
#'   tau[k] <- runif(1,.2,.6)
#' }
#' R = matrix(0,K,K)
#' p_mastery <- c(.5,.5,.4,.4)
#' Alphas_0 <- matrix(0,N,K)
#' for(i in 1:N){
#'   for(k in 1:K){
#'     prereqs <- which(R[k,]==1)
#'     if(length(prereqs)==0){
#'       Alphas_0[i,k] <- rbinom(1,1,p_mastery[k])
#'     }
#'     if(length(prereqs)>0){
#'       Alphas_0[i,k] <- prod(Alphas_0[i,prereqs])*rbinom(1,1,p_mastery)
#'     }
#'   }
#' }
#' Alphas <- sim_alphas(model="indept", taus=tau, N=N, L=L, R=R)
#' 
#' ## FOHM ##
#' N = dim(Design_array)[1]
#' K = ncol(Q_matrix)
#' L = dim(Design_array)[3]
#' TP <- TPmat(K)
#' Omega_true <- rOmega(TP)
#' class_0 <- sample(1:2^K, N, replace = L)
#' Alphas_0 <- matrix(0,N,K)
#' for(i in 1:N){
#'   Alphas_0[i,] <- inv_bijectionvector(K,(class_0[i]-1))
#' }
#' Alphas <- sim_alphas(model="FOHM", Omega = Omega_true, N=N, L=L)
#' }
#' @export
sim_alphas <- function(model, lambdas = NULL, thetas = NULL, Q_matrix = NULL, Design_array = NULL, taus = NULL, Omega = NULL, N = NA_integer_, L = NA_integer_, R = NULL, alpha0 = NULL) {
    .Call(`_hmcdm_sim_alphas`, model, lambdas, thetas, Q_matrix, Design_array, taus, Omega, N, L, R, alpha0)
}