# This file contains functions to call phyclust in C.
phyclust.se <- function(X, K, EMC = .EMC, manual.id = NULL, byrow = TRUE){
label <- NULL
if(K <= 0){
stop("K > 0")
}
if(! is.null(manual.id)){
if(! is.null(label)){
stop("The manual.id is only for unsupervised clustering.")
}
if(any(manual.id <= 0 | manual.id > K)){
stop("The manual.id is not correct.")
}
manual.id <- manual.id - 1
} else{
if(EMC$init.method == "manualMu"){
stop("The manual.id is missing.")
}
}
if(byrow){
X <- t(X)
}
N.X.org <- ncol(X)
L <- nrow(X)
EMC <- check.EMC(EMC)
EMC <- translate.EMC.se(EMC)
label <- check.label(label, N.X.org, K, byrow)
ret <- .Call("R_phyclust_se",
as.integer(N.X.org),
as.integer(L),
as.integer(K),
as.integer(X),
EMC,
as.integer(manual.id),
label,
PACKAGE = "phyclust")
if(!is.finite(ret$logL)){
stop("The logL is not finite.\n")
}
ret$class.id <- ret$class.id + 1
ret <- translate.ret.se(ret, EMC)
class(ret) <- "phyclust"
ret
} # End of phyclust().
### For internal used.
translate.ret.se <- function(ret, EMC = NULL){
ret$Mu <- matrix(ret$Mu, nrow = ret$K, byrow = TRUE)
if(! is.null(ret$Z.normalized)){
ret$Z.normalized <- matrix(ret$Z.normalized, nrow = ret$N.X.org, byrow = TRUE)
}
ret$QA$pi <- matrix(ret$QA$pi, nrow = ret$K, byrow = TRUE)
ret$QA$kappa <- matrix(ret$QA$kappa, nrow = ret$K, byrow = TRUE)
ret$QA$Tt <- matrix(ret$QA$Tt, nrow = ret$K, byrow = TRUE)
if(!is.null(EMC)){
ret$init.procedure <- .init.procedure[EMC$init.procedure + 1]
ret$init.method <- .init.method[EMC$init.method + 1]
ret$substitution.model <-
as.character(.substitution.model$model[EMC$substitution.model + 1])
ret$edist.model <- .edist.model[EMC$edist.model + 1]
ret$QA$identifier <- .identifier[EMC$identifier + 1]
ret$code.type <- .code.type[EMC$code.type + 1]
ret$em.method <- .em.method[EMC$em.method + 1]
ret$boundary.method <- .boundary.method[EMC$boundary.method + 1]
}
if(ret$QA$identifier == "EE"){
ret$QA$pi <- matrix(ret$QA$pi[1,], nrow = 1)
ret$QA$kappa <- ret$QA$kappa[1]
ret$QA$Tt <- ret$QA$Tt[1]
} else if(ret$QA$identifier == "EV"){
ret$QA$pi <- matrix(ret$QA$pi[1,], nrow = 1)
ret$QA$kappa <- ret$QA$kappa[1]
} else if(ret$QA$identifier == "VE"){
ret$QA$Tt <- ret$QA$Tt[1]
}
if(ret$code.type == "NUCLEOTIDE"){
colnames(ret$QA$pi) <- as.character(.nucleotide$code[1:ncol(ret$QA$pi)])
} else if(ret$code.type == "SNP"){
colnames(ret$QA$pi) <- as.character(.snp$code[1:ncol(ret$QA$pi)])
}
rownames(ret$QA$pi) <- paste("k=", 1:nrow(ret$QA$pi), sep = "")
if(ret$substitution.model %in% c("JC69", "K80", "SNP_JC69")){
ret$QA$pi <- NULL
}
if(ret$substitution.model %in%
c("JC69", "F81", "SNP_JC69", "SNP_F81", "E_F81", "E_SNP_F81")){
ret$QA$kappa <- NULL
}
# ret$label.method <- .label.method[ret$label.method + 1]
### For se model.
ret$SE$type <- as.logical(ret$SE$type)
ret$SE$model <- .se.model[ret$SE$model + 1]
ret$SE$f.err <- matrix(ret$SE$f.err, nrow = 4, byrow = TRUE)
colnames(ret$SE$f.err) <- as.character(.nucleotide$code[1:ncol(ret$SE$f.err)])
rownames(ret$SE$f.err) <- as.character(.nucleotide$code[1:nrow(ret$SE$f.err)])
ret
} # End of translate.ret.se().
translate.EMC.se <- function(EMC){
EMC$exhaust.iter <- as.integer(EMC$exhaust.iter)
EMC$fixed.iter <- as.integer(EMC$fixed.iter)
EMC$short.iter <- as.integer(EMC$short.iter)
EMC$EM.iter <- as.integer(EMC$EM.iter)
EMC$short.eps <- as.double(EMC$short.eps)
EMC$EM.eps <- as.double(EMC$EM.eps)
EMC$cm.reltol <- as.double(EMC$cm.reltol)
EMC$cm.maxit <- as.integer(EMC$cm.maxit)
EMC$nm.abstol.Mu.given.QA <- as.double(EMC$nm.abstol.Mu.given.QA)
EMC$nm.abstol.QA.given.Mu <- as.double(EMC$nm.abstol.QA.given.Mu)
EMC$nm.reltol.Mu.given.QA <- as.double(EMC$nm.reltol.Mu.given.QA)
EMC$nm.reltol.QA.given.Mu <- as.double(EMC$nm.reltol.QA.given.Mu)
EMC$nm.maxit.Mu.given.QA <- as.integer(EMC$nm.maxit.Mu.given.QA)
EMC$nm.maxit.QA.given.Mu <- as.integer(EMC$nm.maxit.QA.given.Mu)
EMC$est.non.seg.site <- as.integer(EMC$est.non.seg.site)
EMC$min.n.class <- as.integer(EMC$min.n.class)
if(EMC$init.procedure[1] %in% .init.procedure){
EMC$init.procedure <- as.integer(which(EMC$init.procedure ==
.init.procedure) - 1)
} else{
stop("The initial procedure is not found.")
}
if(EMC$init.method[1] %in% .init.method){
EMC$init.method <- as.integer(which(EMC$init.method == .init.method) - 1)
} else{
stop("The initial method is not found.")
}
if(EMC$substitution.model[1] %in% as.character(.substitution.model$model)){
EMC$substitution.model <-
as.integer(which(EMC$substitution.model ==
as.character(.substitution.model$model)) - 1)
} else{
stop("The substitution model is not found.")
}
if(EMC$edist.model[1] %in% .edist.model){
EMC$edist.model <- as.integer(which(EMC$edist.model == .edist.model) - 1)
} else{
stop("The edist model is not found.")
}
if(EMC$identifier[1] %in% .identifier){
EMC$identifier <- as.integer(which(EMC$identifier == .identifier) - 1)
} else{
stop("The identifier is not found.")
}
if(EMC$code.type[1] %in% .code.type){
EMC$code.type <- as.integer(which(EMC$code.type == .code.type) - 1)
} else{
stop("The code type is not found.")
}
if(EMC$em.method[1] %in% .em.method){
EMC$em.method <- as.integer(which(EMC$em.method == .em.method) - 1)
} else{
stop("The em method is not found.")
}
if(EMC$boundary.method[1] %in% .boundary.method){
EMC$boundary.method <- as.integer(which(EMC$boundary.method ==
.boundary.method) - 1)
} else{
stop("The boundary method is not found.")
}
EMC$max.init.iter <- as.integer(EMC$max.init.iter)
### For se model.
if(EMC$code.type == as.integer(which(.code.type == "NUCLEOTIDE") - 1)){
EMC$se.type <- as.integer(1)
EMC$se.model <- as.integer(which(EMC$se.model == .se.model) - 1)
EMC$se.constant <- as.double(EMC$se.constant)
} else{
stop("The sequencing error model only supports NUCLEOTIDE data.")
}
EMC
} # End of translate.EMC.se().
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.