Nothing
# 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().
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.