Nothing
funHDDC <-
function(data, K=1:10, model="AkjBkQkDk", threshold=0.2, itermax=200, eps=1e-6,init='kmeans', criterion="bic",
algo='EM', d_select="Cattell", init.vector=NULL, show=TRUE, mini.nb=c(5, 10),
min.individuals=2, mc.cores=1, nb.rep=1, keepAllRes=TRUE, kmeans.control = list(), d_max=100){
#Options removed from call
com_dim=NULL
scaling=FALSE
noise.ctrl=1e-8
#nb.rep parameter
if ((init=="random")&(nb.rep<20)){
nb.rep=20
}
#
# CONTROLS
#
call = match.call()
.hddc_control(call)
# Control of match.args:
criterion = .myAlerts(criterion, "criterion", "singleCharacterMatch.arg", "HDDC: ", c("bic", "icl"))
algo = .myAlerts(algo, "algo", "singleCharacterMatch.arg", "HDDC: ", c('EM', 'CEM', 'SEM'))
d_select = .myAlerts(d_select, "d_select", "singleCharacterMatch.arg", "HDDC: ", c("cattell", "bic"))
init = .myAlerts(init, "init", "singleCharacterMatch.arg", "HDDC: ", c('random', 'kmeans', 'mini-em', 'param', "vector"))
# We get the model names, properly ordered
model = .hdc_getTheModel(model, all2models = TRUE)
# kmeans controls
kmeans.control = .default_kmeans_control(kmeans.control)
if (scaling) {
stop('Scaling is not implemented!')
} else scaling <- NULL
BIC <- ICL <- c()
fdobj = data
if (class(fdobj)!='list') {x = t(fdobj$coefs)}
else {x = t(fdobj[[1]]$coefs); for (i in 2:length(fdobj)) x = cbind(x,t(fdobj[[i]]$coefs))}
p <- ncol(x)
#
# Preparing the parallel
#
if(d_select=="bic"){
# If the dimension selection is done with BIC, we don't care of the threshold
threshold = "bic"
}
if(max(table(K))>1) warning("The number of clusters, K, is made unique (repeated values are not tolerated).")
K = sort(unique(K))
if(any(K==1)){
# K=1 => only one model required
K = K[K!=1]
addrows = data.frame(model="AKJBKQKDK", K=1, threshold)
} else {
addrows = c()
}
mkt_Expand = expand.grid(model=model, K=K, threshold=threshold)
mkt_Expand = do.call(rbind, replicate(nb.rep, mkt_Expand, simplify=FALSE))
mkt_Expand = rbind(addrows, mkt_Expand) #no need for several runs for K==1
model = as.character(mkt_Expand$model)
K = mkt_Expand$K
threshold = mkt_Expand$threshold
# We transform it into an univariate form
mkt_univariate = apply(mkt_Expand, 1, paste, collapse= "_")
# Mon 'caller' for LTBM that will be used in multi-cores lapply
hddcWrapper = function(mkt_univariate, ...){
mkt_splitted = strsplit(mkt_univariate, "_")
# on retrouve model, K and threshold
model = sapply(mkt_splitted, function(x) x[1])
K = sapply(mkt_splitted, function(x) as.numeric(x[2]))
threshold = sapply(mkt_splitted, function(x) ifelse(x[3]=="bic","bic",as.numeric(x[3])))
# (::: is needed for windows multicore)
res = "unknown error"
# try(res <- HDclassif:::funhddc_main(model=model, K=K, threshold=threshold, ...))
try(res <- .funhddc_main(model=model, K=K, threshold=threshold, ...),silent=TRUE)
res
}
# We reset the number of cores to use
nRuns = length(mkt_univariate)
if(nRuns<mc.cores) mc.cores = nRuns
# We swicth to the right number of cores + a warning if necessary
max_nb_of_cores = parallel::detectCores()
if(mc.cores>max_nb_of_cores){
warning("The argument mc.cores is greater than its maximun.\nmc.cores was set to ", max_nb_of_cores)
mc.cores = max_nb_of_cores
}
#
# Parallel estimations
#
if(mc.cores == 1){
# If there is no need for parallel, we just use lapply // in order to have the same output
par.output = lapply(mkt_univariate, hddcWrapper, fdobj=fdobj, method=d_select, algo=algo, itermax=itermax, eps=eps, init=init, init.vector=init.vector, mini.nb=mini.nb, min.individuals=min.individuals, noise.ctrl=noise.ctrl, com_dim=com_dim, kmeans.control=kmeans.control, d_max=d_max)
} else if(Sys.info()[['sysname']] == 'Windows'){
# we use parLapply:
## create clusters
cl = parallel::makeCluster(mc.cores)
## load the packages
loadMyPackages = function(x){
# loadMyPackages = function(none, myFuns){
# we load the package
library(funHDDC)
# we add the functions in the global env
# for(i in 1:length(myFuns)){
# funName = names(myFuns)[i]
# fun = myFuns[[i]]
# assign(funName, fun, .GlobalEnv)
# }
}
## Create the functions to export
# myFuns = list(hddc_main = hddc_main, hddc_e_step=hddc_e_step, hddc_m_step=hddc_m_step, hdc_getComplexity=hdc_getComplexity, hdc_myEigen=hdc_myEigen, hdclassif_dim_choice=hdclassif_dim_choice, hdclassif_bic=hdclassif_bic)
# par.setup = parallel::parLapply(cl, 1:length(cl), loadMyPackages, myFuns=myFuns)
par.setup = parallel::parLapply(cl, 1:length(cl), loadMyPackages)
## run the parallel
par.output = NULL
try(par.output <- parallel::parLapply(cl, mkt_univariate, hddcWrapper, DATA=fdobj, method=d_select, algo=algo, itermax=itermax, eps=eps, init=init, init.vector=ifelse(missing(init.vector), NA, init.vector), mini.nb=mini.nb, min.individuals=min.individuals, noise.ctrl=noise.ctrl, com_dim=com_dim, kmeans.control=kmeans.control, d_max=d_max))
## Stop the clusters
parallel::stopCluster(cl)
if(is.null(par.output)) stop("Unknown error in the parallel computing. Try mc.cores=1 to detect the problem.")
} else {
# we use mclapply
par.output = NULL
try(par.output <- parallel::mclapply(mkt_univariate, hddcWrapper, DATA=fdobj, method=d_select, algo=algo, itermax=itermax, eps=eps, init=init, init.vector=init.vector, mini.nb=mini.nb, min.individuals=min.individuals, noise.ctrl=noise.ctrl, com_dim=com_dim, kmeans.control=kmeans.control, d_max=d_max, mc.cores=mc.cores))
if(is.null(par.output)) stop("Unknown error in the parallel computing. Try mc.cores=1 to detect the problem.")
}
#
# The results are retrieved
#
getElement = function(x, what, valueIfNull = -Inf){
# attention si x est le modele nul
if(length(x)==1) return(valueIfNull)
if(!is.list(x) && !what %in% names(x)) return(NA)
x[[what]][length(x[[what]])]
}
getComment = function(x){
# we get the error message
if(length(x)==1) return(x)
return("")
}
# All likelihoods
LL_all = sapply(par.output, getElement, what="loglik")
comment_all = sapply(par.output, getComment)
# If no model is valid => problem
if(all(!is.finite(LL_all))){
warning("All models diverged.")
allCriteria = data.frame(model=model, K=K, threshold=threshold, LL = LL_all, BIC=NA, comment=comment_all)
res = list()
res$allCriteria = allCriteria
return(res)
}
# We select, for each (Q,K), the best run
n = nrow(mkt_Expand)
modelKeep = sapply(unique(mkt_univariate), function(x) (1:n)[mkt_univariate==x][which.max(LL_all[mkt_univariate==x])])
# => we select only the best models
LL_all = LL_all[modelKeep]
comment_all = comment_all[modelKeep]
par.output = par.output[modelKeep]
BIC = sapply(par.output, getElement, what="BIC")
ICL = sapply(par.output, getElement, what="ICL")
comp_all = sapply(par.output, getElement, what="complexity", valueIfNull=NA)
model = model[modelKeep]
threshold = threshold[modelKeep]
K = K[modelKeep]
# We define the criterion of model selection
CRIT = switch(criterion,
bic = BIC,
icl = ICL)
# The order of the results
myOrder = order(CRIT, decreasing = TRUE)
# On sauvegarde le bon modele + creation de l'output
qui = which.max(CRIT)
prms = par.output[[qui]]
prms$criterion = CRIT[qui]
names(prms$criterion) = criterion
# Other output
prms$call = call
# We add the complexity
names(comp_all) = mkt_univariate[modelKeep]
prms$complexity_allModels = comp_all
# Display
if(show){
if(n>1) cat("funHDDC: \n")
model2print = sapply(model, function(x) sprintf("%*s", max(nchar(model)), x))
K2print = as.character(K)
K2print = sapply(K2print, function(x) sprintf("%*s", max(nchar(K2print)), x))
thresh2print = as.character(threshold)
thresh_width = max(nchar(thresh2print))
thresh2print = sapply(thresh2print, function(x) sprintf("%s%s", x, paste0(rep("0", thresh_width - nchar(x)), collapse="") ))
# on cree une data.frame
myResMat = cbind(model2print[myOrder], K2print[myOrder], thresh2print[myOrder],.addCommas(prms$complexity_allModels[myOrder]), .addCommas(CRIT[myOrder]), comment_all[myOrder])
myResMat = as.data.frame(myResMat)
names(myResMat) = c("model", "K", "threshold", "complexity",toupper(criterion), "comment")
row.names(myResMat) = 1:nrow(myResMat)
# if no problem => no comment
if(all(comment_all == "")) myResMat$comment = NULL
print(myResMat)
msg = switch(criterion, bic="BIC", icl="ICL")
cat("\nSELECTED: model ", prms$model, " with ", prms$K, " clusters.\n")
cat("Selection Criterion: ", msg, ".\n", sep="")
}
# We also add the matrix of all criteria
allCriteria = data.frame(model=model[myOrder], K=K[myOrder], threshold=threshold[myOrder], LL=LL_all[myOrder],complexity=prms$complexity_allModels[myOrder], BIC=BIC[myOrder], ICL=ICL[myOrder], rank = 1:length(myOrder))
# we add the comments if necessary
if(any(comment_all != "")) allCriteria$comment = comment_all[myOrder]
prms$allCriteria = allCriteria
# If all the results are kept
if(keepAllRes){
all_results = par.output
names(all_results) = mkt_univariate[modelKeep]
prms$all_results = all_results
}
# Other stuff
prms$scaling <- scaling
prms$threshold <- threshold[qui]
return(prms)
}
.funhddc_main <- function(fdobj, K, model, threshold, method, algo, itermax, eps, init, init.vector, mini.nb, min.individuals, noise.ctrl, com_dim=NULL, kmeans.control, d_max, ...){
ModelNames <- c("AKJBKQKDK", "AKBKQKDK", "ABKQKDK", "AKJBQKDK", "AKBQKDK", "ABQKDK", "AKJBKQKD", "AKBKQKD", "ABKQKD", "AKJBQKD", "AKBQKD", "ABQKD")
if (class(fdobj)!='list') {DATA = t(fdobj$coefs)}
else {DATA = t(fdobj[[1]]$coefs); for (i in 2:length(fdobj)) DATA = cbind(DATA,t(fdobj[[i]]$coefs))}
p <- ncol(DATA)
N <- nrow(DATA)
com_ev <- NULL
# We set d_max to a proper value
d_max = min(N, p, d_max)
# if ( any(model==ModelNames[7:12]) ){
# # Common dimension models
# MU <- colMeans(DATA)
# if (N<p) {
# Y <- (DATA-matrix(MU, N, p, byrow=TRUE))/sqrt(N)
# YYt <- tcrossprod(Y)
# com_ev <- hdc_myEigen(YYt, d_max, only.values = TRUE)$values
# } else{
# S <- crossprod(DATA-matrix(MU, N, p, byrow=TRUE))/N
# com_ev <- hdc_myEigen(S, d_max, only.values = TRUE)$values
# }
# if(is.null(com_dim)) com_dim <- .hdclassif_dim_choice(com_ev, N, method, threshold, FALSE, noise.ctrl)
# }
if (K>1){
t <- matrix(0, N, K)
if(init == "vector"){
init.vector = unclass(init.vector)
name <- unique(init.vector)
for (i in 1:K) t[which(init.vector==name[i]), i] <- 1
} else if (init=='kmeans') {
kmc = kmeans.control
cluster <- kmeans(DATA, K, iter.max=kmc$iter.max, nstart=kmc$nstart, algorithm=kmc$algorithm, trace=kmc$trace)$cluster
for (i in 1:K) t[which(cluster==i), i] <- 1
} else if (init=='mini-em'){
prms_best <- 1
for (i in 1:mini.nb[1]){
prms <- .funhddc_main(DATA, K, model, threshold, method, algo, mini.nb[2], 0, 'random', mini.nb = mini.nb, min.individuals = min.individuals, noise.ctrl = noise.ctrl, com_dim = com_dim, d_max=d_max)
if(length(prms)!=1){
if (length(prms_best)==1) prms_best <- prms
else if (prms_best$loglik[length(prms_best$loglik)]<prms$loglik[length(prms$loglik)]) prms_best <- prms
}
}
if (length(prms_best)==1) return(1)
t <- prms_best$posterior
} else {
t <- t(rmultinom(N, 1, rep(1/K, K)))
compteur=1
while(min(colSums(t))<1 && (compteur <- compteur+1)<5) t <- t(rmultinom(N, 1, rep(1/K, K)))
if(min(colSums(t))<1) return("Random initialization failed (n too small)")
}
} else t <- matrix(1, N, 1)
likely <- c()
I <- 0
test <- Inf
while ((I <- I+1)<=itermax && test>=eps){
if (algo!='EM' && I!=1) t <- t2
# Error catching
if (K>1){
if(any(is.na(t))) return("unknown error: NA in t_ik")
if(any(colSums(t>1/K)<min.individuals)) return("pop<min.individuals")
}
m <- .funhddc_m_step(fdobj, K, t, model, threshold, method, noise.ctrl, com_dim, d_max)
t <- .funhddc_e_step(fdobj, m)
L <- t$L
t <- t$t
if (algo=='CEM') {
t2 <- matrix(0, N, K)
t2[cbind(1:N, max.col(t))] <- 1
} else if(algo=='SEM') {
t2 <- matrix(0, N, K)
for (i in 1:N) t2[i, ] <- t(rmultinom(1, 1, t[i, ]))
}
likely[I] <- L
if (I!=1) test <- abs(likely[I]-likely[I-1])
}
# We retrieve the parameters
# a
if ( model%in%c('AKBKQKDK', 'AKBQKDK', 'AKBKQKD', 'AKBQKD') ) {
a <- matrix(m$a[, 1], 1, m$K, dimnames=list(c("Ak:"), 1:m$K))
} else if(model=='AJBQD') {
a <- matrix(m$a[1, ], 1, m$d[1], dimnames=list(c('Aj:'), paste('a', 1:m$d[1], sep='')))
} else if ( model%in%c('ABKQKDK', 'ABQKDK', 'ABKQKD', 'ABQKD', "ABQD") ) {
a <- matrix(m$a[1], dimnames=list(c('A:'), c('')))
} else a <- matrix(m$a, m$K, max(m$d), dimnames=list('Class'=1:m$K, paste('a', 1:max(m$d), sep='')))
# b
if ( model%in%c('AKJBQKDK', 'AKBQKDK', 'ABQKDK', 'AKJBQKD', 'AKBQKD', 'ABQKD', 'AJBQD', "ABQD") ) {
b <- matrix(m$b[1], dimnames=list(c('B:'), c('')))
} else b <- matrix(m$b, 1, m$K, dimnames=list(c("Bk:"), 1:m$K))
# d, mu, prop
d <- matrix(m$d, 1, m$K, dimnames=list(c('dim:'), "Intrinsic dimensions of the classes:"=1:m$K))
mu <- matrix(m$mu, m$K, p, dimnames=list('Class'=1:m$K, 'Posterior group means:'=paste('V', 1:p, sep='')))
prop <- matrix(m$prop, 1, m$K, dimnames=list(c(''), 'Posterior probabilities of groups'=1:m$K))
# Other elements
complexity <- .hdc_getComplexity(m, p)
class(b) <- class(a) <- class(d) <- class(prop) <- class(mu) <- 'hd'
cls <- max.col(t)
converged = test<eps
params = list(model=model, K=K, d=d, a=a, b=b, mu=mu, prop=prop, ev=m$ev, Q=m$Q,fpca=m$fpcaobj, loglik=likely[length(likely)], loglik_all = likely, posterior=t, class=cls, com_ev=com_ev, N=N, complexity=complexity, threshold=threshold, d_select=method, converged=converged)
# We compute the BIC / ICL
bic_icl = .hdclassif_bic(params, p)
params$BIC = bic_icl$bic
params$ICL = bic_icl$icl
# We set the class
class(params) <- 'funHDDC'
return(params)
}
.funhddc_e_step <- function(fdobj, par){
if (class(fdobj)!='list') {x = t(fdobj$coefs)}
else {x = t(fdobj[[1]]$coefs); for (i in 2:length(fdobj)) x = cbind(x,t(fdobj[[i]]$coefs))}
p <- ncol(x)
N <- nrow(x)
K <- par$K
a <- par$a
b <- par$b
mu <- par$mu
d <- par$d
prop <- par$prop
Q <- par$Q
b[b<1e-6] <- 1e-6
K_pen <- matrix(0,K,N)
for (i in 1:K) {
s <- sum(log(a[i,1:d[i]]))
X <- x - matrix(mu[i,], N, p, byrow=TRUE)
Qi = par$fpcaobj[[i]]$W %*% Q[[i]]
proj <- (X%*%Qi)%*%t(Qi)
A <- (-proj)%*%Qi%*%sqrt(diag(1/a[i,1:d[i]],d[i]))
B <- X-proj
K_pen[i,] <- rowSums(A^2)+1/b[i]*rowSums(B^2)+s+(p-d[i])*log(b[i])-2*log(prop[i])+p*log(2*pi)
}
A <- -1/2*t(K_pen)
L <- sum(log(rowSums(exp(A-apply(A,1,max))))+apply(A,1,max))
t <- matrix(0,N,K)
for (i in 1:K) t[,i] <- 1/rowSums(exp((K_pen[i,]-t(K_pen))/2))
list(t=t,L=L)
}
.funhddc_m_step <- function(fdobj, K, t, model, threshold, method, noise.ctrl, com_dim, d_max){
if (class(fdobj)!='list') { x = t(fdobj$coefs)
} else {x = t(fdobj[[1]]$coefs); for (i in 2:length(fdobj)) x = cbind(x,t(fdobj[[i]]$coefs))}
# Some parameters
N <- nrow(x)
p <- ncol(x)
prop <- c()
n <- colSums(t)
prop <- n/N
mu <- matrix(NA, K, p)
for (i in 1:K) mu[i, ] <- colSums(x*t[, i])/n[i]
ind <- apply(t>0, 2, which)
n_bis <- c()
for(i in 1:K) n_bis[i] <- length(ind[[i]])
#
#Calculation on Var/Covar matrices
#
# we keep track of the trace (== sum of eigenvalues) to compute the b
traceVect = c()
ev <- matrix(0, K, p)
Q <- vector(mode='list', length=K)
fpcaobj = list()
for (i in 1:K){
donnees <- .mypca.fd(fdobj,t[,i])
traceVect[i] = sum(diag(donnees$valeurs_propres))
ev[i, ] <- donnees$valeurs_propres
Q[[i]] <- donnees$U
fpcaobj[[i]] = donnees
}
#Intrinsic dimensions selection
# browser()
if (model%in%c("AJBQD", "ABQD")){
d <- rep(com_dim, length=K)
} else if ( model%in%c("AKJBKQKD", "AKBKQKD", "ABKQKD", "AKJBQKD", "AKBQKD", "ABQKD") ){
dmax <- min(apply((ev>noise.ctrl)*rep(1:ncol(ev), each=K), 1, which.max))-1
if(com_dim>dmax) com_dim <- max(dmax, 1)
d <- rep(com_dim, length=K)
} else {
d <- .hdclassif_dim_choice(ev, n, method, threshold, FALSE, noise.ctrl)
}
#Setup of the Qi matrices
for(i in 1:K) Q[[i]] <- matrix(Q[[i]][, 1:d[i]], p, d[i])
#Calculation of the remaining parameters of the selected model
# PARAMETER a
ai <- matrix(NA, K, max(d))
if ( model%in%c('AKJBKQKDK', 'AKJBQKDK', 'AKJBKQKD', 'AKJBQKD') ){
for (i in 1:K) ai[i, 1:d[i]] <- ev[i, 1:d[i]]
} else if ( model%in%c('AKBKQKDK', 'AKBQKDK' , 'AKBKQKD', 'AKBQKD') ){
for (i in 1:K) ai[i, ] <- rep(sum(ev[i, 1:d[i]])/d[i], length=max(d))
} else if(model=="AJBQD"){
for (i in 1:K) ai[i, ] <- ev[1:d[1]]
} else if(model=="ABQD") {
ai[] <- sum(ev[1:d[1]])/d[1]
} else {
a <- 0
eps <- sum(prop*d)
for (i in 1:K) a <- a + sum(ev[i, 1:d[i]])*prop[i]
ai <- matrix(a/eps, K, max(d))
}
# PARAMETER b
bi <- c()
denom = min(N,p)
if ( model%in%c('AKJBKQKDK', 'AKBKQKDK', 'ABKQKDK', 'AKJBKQKD', 'AKBKQKD', 'ABKQKD') ){
for(i in 1:K){
remainEV = traceVect[i] - sum(ev[i, 1:d[i]])
# bi[i] <- sum(ev[i, (d[i]+1):min(N, p)])/(p-d[i])
bi[i] <- remainEV/(p-d[i]) #pour moi c'est p au lieu de denom
}
} else if ( model%in%c("ABQD", "AJBQD") ){
remainEV = traceVect - sum(ev[1:d[1]])
# bi[1:K] <- sum(ev[(d[1]+1):min(N, p)])/(min(N, p)-d[1])
bi[1:K] <- remainEV/(denom-d[1])
} else {
b <- 0
eps <- sum(prop*d)
for(i in 1:K){
remainEV = traceVect[i] - sum(ev[i, 1:d[i]])
# b <- b + sum(ev[i, (d[i]+1):min(N, p)])*prop[i]
b <- b + remainEV*prop[i]
}
bi[1:K] <- b/(min(N, p)-eps)
}
# We adjust the values of b if they are too low
#bi[bi<noise.ctrl] = noise.ctrl
list(model=model, K=K, d=d, a=ai, b=bi, mu=mu, prop=prop, ev=ev, Q=Q,fpcaobj = fpcaobj)
}
.mypca.fd <- function(fdobj, Ti){
if (class(fdobj)=="list"){
#sauvegarde de la moyenne avant centrage
mean_fd<-list()
for (i in 1:length(fdobj)){
mean_fd[[i]]<-fdobj[[i]]
}
#centrage des objets fonctionnels
for (i in 1:length(fdobj)){
coefmean <- apply(t(as.matrix(Ti) %*% matrix(1,1,nrow(fdobj[[i]]$coefs))) * fdobj[[i]]$coefs, 1, sum) / sum(Ti)
fdobj[[i]]$coefs <- sweep(fdobj[[i]]$coefs, 1, coefmean)
mean_fd[[i]]$coefs = as.matrix(data.frame(mean=coefmean))
}
#Cr?ation des matrices des produits des bases de fonction
for (i in 1:length(fdobj)){
name<-paste('W_var',i,sep='')
W_fdobj<-inprod(fdobj[[i]]$basis,fdobj[[i]]$basis)
assign(name,W_fdobj)
}
#Ajout des 0 ? gauche et ? droite des matrices W avant leur fusion en matrice phi
prow<-dim(W_fdobj)[[1]]
pcol<-length(fdobj)*prow
W1<-cbind(W_fdobj,matrix(0,nrow=prow,ncol=(pcol-ncol(W_fdobj))))
W_list<-list()
for (i in 2:(length(fdobj))){
W2<-cbind(matrix(0,nrow=prow,ncol=(i-1)*ncol(W_fdobj)),get(paste('W_var',i,sep='')),matrix(0,nrow=prow,ncol=(pcol-i*ncol(W_fdobj))))
W_list[[i-1]]<-W2
}
#Cr?ation de la matrice phi
W_tot<-rbind(W1,W_list[[1]])
if (length(fdobj)>2){
for(i in 2:(length(fdobj)-1)){
W_tot<-rbind(W_tot,W_list[[i]])
}
}
W_tot[W_tot<1e-15]=0
#Cr?ation de la matrice de coef
coef<-t(fdobj[[1]]$coefs)
for (i in 2:length(fdobj)){
coef<-cbind(coef,t(fdobj[[i]]$coefs))
}
#mat_interm<-1/sqrt((sum(Ti)-1))*coef%*%(W_tot)^(1/2)
#D?but Partie mise en commentaire
#mat_interm<-1/sqrt((sum(Ti)))*coef%*%chol(W_tot,pivot=TRUE)
#cov<-(.repmat(Ti,n=pcol,p=1)*t(mat_interm))%*%mat_interm
#Fin de partie mise en commentaire
#D?but correction
#Construction matrice triangulaire de Choleski
W_m <- chol(W_tot)
#Matrice de covariance
mat_cov <- crossprod(t(.repmat(sqrt(Ti),n=dim(t(coef))[[1]],p=1)*t(coef)))/sum(Ti)
cov = W_m %*% mat_cov %*% t(W_m)
#Fin correction
valeurs<-Eigen(cov)
valeurs_propres<-valeurs$values
vecteurs_propres<-valeurs$vectors
#bj<-solve((W_tot)^(1/2))%*%vecteurs_propres
bj<-solve(chol(W_tot))%*%vecteurs_propres
fonctionspropres<-fdobj[[1]]
fonctionspropres$coefs<-bj
scores<-coef%*%W_tot%*%bj
varprop<-valeurs_propres/sum(valeurs_propres)
pcafd<-list(valeurs_propres=valeurs_propres,harmonic=fonctionspropres,scores=scores,covariance=cov,U=bj,varprop=varprop,meanfd=mean_fd,W=W_tot)
}else if (class(fdobj)!="list") {
#Calcul de la moyenne par groupe
mean_fd<-fdobj
#Centrer les objets fonctionnels par groupe
coefmean <- apply(t(as.matrix(Ti) %*% matrix(1,1,nrow(fdobj$coefs))) * fdobj$coefs, 1, sum) / sum(Ti)
fdobj$coefs <- sweep(fdobj$coefs, 1, coefmean)
mean_fd$coefs = as.matrix(data.frame(mean=coefmean))
#Calcul de la matrice des produits scalaires
W<-inprod(fdobj$basis,fdobj$basis)
#pour ?viter les soucis num?riques on arrondit ? 0 les produits scalaires tr?s petits
W[W<1e-15]=0
#on centre les coefficients
coef<-t(fdobj$coefs)
#mat_interm<-1/sqrt((sum(Ti)-1))*coef%*%(W)^(1/2)
#D?but partie mise en commentaire
#mat_interm<-1/sqrt((sum(Ti)))*coef%*%chol(W)
#cov<-(.repmat(Ti,n=dim(W)[[1]],p=1)*t(mat_interm))%*%mat_interm
#Fin de la partie mise en commentaire
#D?but correction
#Construction matrice triangulaire de Choleski
W_m <- chol(W)
#Matrice de covariance
mat_cov <- crossprod(t(.repmat(sqrt(Ti),n=dim(t(coef))[[1]],p=1)*t(coef)))/sum(Ti)
cov = W_m %*% mat_cov %*% t(W_m)
#Fin correction
valeurs<-Eigen(cov)
valeurs_propres<-valeurs$values
vecteurs_propres<-valeurs$vectors
#Calcul de U
#U<-chol(W)%*%vecteurs_propres
#Calcul des coefficients des fonctions propres
fonctionspropres<-fdobj
bj<-solve(chol(W))%*%vecteurs_propres
fonctionspropres$coefs<-bj
#calcul des scores selon la formule de pca.fd
scores<-inprod(fdobj,fonctionspropres)
varprop<-valeurs_propres/sum(valeurs_propres)
pcafd <-list(valeurs_propres=valeurs_propres,harmonic=fonctionspropres,scores=scores,covariance=cov,U=bj,meanfd=mean_fd,W=W)
}
class(pcafd) <- "pca.fd"
return(pcafd)
}
.hddc_ari <- function(x,y){
#This function is drawn from the mclust package
x <- as.vector(x)
y <- as.vector(y)
tab <- table(x, y)
if (all(dim(tab) == c(1, 1))) return(1)
a <- sum(choose(tab, 2))
b <- sum(choose(rowSums(tab), 2)) - a
c <- sum(choose(colSums(tab), 2)) - a
d <- choose(sum(tab), 2) - a - b - c
ARI <- (a - (a + b) * (a + c)/(a + b + c + d))/((a + b + a + c)/2 - (a + b) * (a + c)/(a + b + c + d))
return(ARI)
}
####
#### CONTROLS ####
####
.hddc_control = function(call){
prefix = "HDDC: "
.myCallAlerts(call, "data", "list,fd", 3, TRUE, prefix)
.myCallAlerts(call, "K", "integerVector", 3, FALSE, prefix)
.myCallAlerts(call, "model", "vector", 3, FALSE, prefix)
.myCallAlerts(call, "threshold", "numericVectorGE0LE1", 3, FALSE, prefix)
.myCallAlerts(call, "criterion", "character", 3, FALSE, prefix)
.myCallAlerts(call, "com_dim", "singleIntegerGE1", 3, FALSE, prefix)
.myCallAlerts(call, "itermax", "singleIntegerGE0", 3, FALSE, prefix)
.myCallAlerts(call, "eps", "singleNumericGE0", 3, FALSE, prefix)
.myCallAlerts(call, "graph", "singleLogical", 3, FALSE, prefix)
.myCallAlerts(call, "algo", "singleCharacter", 3, FALSE, prefix)
.myCallAlerts(call, "d_select", "singleCharacter", 3, FALSE, prefix)
.myCallAlerts(call, "init", "singleCharacter", 3, FALSE, prefix)
.myCallAlerts(call, "show", "singleLogical", 3, FALSE, prefix)
.myCallAlerts(call, "mini.nb", "integerVectorGE1", 3, FALSE, prefix)
.myCallAlerts(call, "scaling", "singleLogical", 3, FALSE, prefix)
.myCallAlerts(call, "min.individuals", "singleIntegerGE2", 3, FALSE, prefix)
.myCallAlerts(call, "noise.ctrl", "singleNumericGE0", 3, FALSE, prefix)
.myCallAlerts(call, "mc.cores", "singleIntegerGE1", 3, FALSE, prefix)
.myCallAlerts(call, "nb.rep", "singleIntegerGE1", 3, FALSE, prefix)
.myCallAlerts(call, "keepAllRes", "singleLogical", 3, FALSE, prefix)
.myCallAlerts(call, "d_max", "singleIntegerGE1", 3, FALSE, prefix)
####
#### SPECIFIC controls
####
# Getting some elements
data = eval.parent(call[["data"]], 2)
K = eval.parent(call[["K"]], 2)
init = eval.parent(call[["init"]], 2)
criterion = eval.parent(call[["criterion"]], 2)
if (is.fd(data)){
# No NA in the data:
if (any(is.na(data$coefs))) stop("NA values in the data are not supported. Please remove them beforehand.")
# Size of the data
if(any(K>2*NROW(data$coefs))) stop("The number of observations must be at least twice the number of clusters ")
}else{
# No NA in the data:
for (i in 1:length(data)){
if (any(is.na(data[[i]]$coefs))) stop("NA values in the data are not supported. Please remove them beforehand.")
}
# Size of the data
if(any(K>2*NROW(data[[1]]$coefs))) stop("The number of observations must be at least twice the number of clusters ")
}
# Initialization Controls
if(!is.null(init)){
# we get the value of the initialization
init = .myAlerts(init, "init", "singleCharacterMatch.arg", "HDDC: ", c('random', 'kmeans', 'mini-em', 'param', "vector"))
# Custom initialization => controls and setup
if(init == "vector"){
fdobj = data
if (class(fdobj)!='list') {x = t(fdobj$coefs)}
else {x = t(fdobj[[1]]$coefs); for (i in 2:length(fdobj)) x = cbind(x,t(fdobj[[i]]$coefs))}
.myCallAlerts(call, "init.vector", "(integer,factor)Vector", 3, FALSE, prefix)
init.vector = eval.parent(call[["init.vector"]], 2)
if(is.null(init.vector)) stop("HDDC: When init='vector', the argument 'init.vector' should be provided.")
if(length(unique(K))>1) stop("HDDC: Several number of classes K cannot be estimated when init='vector'.")
init.vector <- unclass(init.vector)
if(K!=max(init.vector)) stop("The number of class K, and the number of classes in the initialization vector are different")
if( length(init.vector)!=nrow(x) ) stop("The size of the initialization vector is different of the size of the data")
}
# The param init
if (init=='param' && nrow(data)<ncol(data)){
stop("The 'param' initialization can't be done when N<p")
}
# The mini.em init
if (init=='mini-em'){
mini.nb = eval.parent(call[["mini.nb"]], 2)
if(!is.null(mini.nb) && length(mini.nb)!=2){
stop("The parameter mini.nb must be a vector of length 2 with integers\n")
}
}
}
}
.default_kmeans_control = function(control){
.myAlerts(control,"kmeans.control","list","kmeans controls: ")
#
# Default values of the control parameters
#
myDefault = list()
myDefault$iter.max = 10
myDefault$nstart = 1
myDefault$algorithm = c("Hartigan-Wong", "Lloyd", "Forgy","MacQueen")
myDefault$trace = FALSE
#
# Types of each arg
#
myTypes = c("singleIntegerGE1", "singleIntegerGE1", "match.arg", "singleLogical")
#
# Recreation of the kmeans controls + Alerts
#
control = .matchTypeAndSetDefault(control, myDefault, myTypes, "kmeans list of controls: ")
return(control)
}
#=================================#
# This file contains all the
# "control" functions
#=================================#
# Possible elements of myAlerts:
#
.myCallAlerts = function(call, name, myType, nParents=1, mustBeThere=FALSE, prefix=""){
# This function basically calls the function myAlerts, but the arguments are different
if( name %in% names(call) ){
# we check the element exists => to provide a fine error
what = call[[name]]
val = try(eval.parent(what, nParents), silent = TRUE)
# browser()
if( "try-error" %in% class(val) ){
if( class(what)=="name" ){
# it means the variable was not found
stop(prefix,"For argument '",name,"': object '",what,"' not found.", call. = FALSE)
} else {
stop(prefix,"For argument '",name,"': expression ",as.character(as.expression(what))," could not be evaluated.", call. = FALSE)
}
} else {
a = .myAlerts(val, name, myType, prefix)
return(a)
}
} else if(mustBeThere) {
stop(prefix, "The argument '", name, "' must be provided.", call. = FALSE)
}
}
.myAlerts = function(x, name, myType, prefix="", charVec){
# Format of my types:
# - single => must be of lenght one
# - Vector => must be a vector
# - Matrix => must be a matrix
# - GE/GT/LE/LT: greater/lower than a given value
# - predefinedType => eg: numeric, integer, etc
# - match.arg => very specific => should match the charVec
# If there is a parenthesis => the class must be of specified types:
# ex: "(list, data.frame)" must be a list of a data.frame
ignore.case = TRUE
firstMsg = paste0(prefix,"The argument '",name,"' ")
# simple function to extract a pattern
# ex: if my type is VectorIntegerGE1 => myExtract("GE[[:digit:]]+","VectorIntegerGE1") => 1
myExtract = function(expr, text, trim=2){
start = gregexpr(expr,text)[[1]] + trim
length = attr(start,"match.length") - trim
res = substr(text,start,start+length-1)
as.numeric(res)
}
#
# General types handling
#
loType = tolower(myType)
if(grepl("single",loType)){
if(length(x)!=1) stop(firstMsg,"must be of length one.", call. = FALSE)
}
if(grepl("vector",loType) && !grepl("factor",loType)){
if(!is.vector(x)) stop(firstMsg,"must be a vector.", call. = FALSE)
if(is.list(x)) stop(firstMsg,"must be a vector (and not a list).", call. = FALSE)
}
res = .checkTheTypes(loType, x)
if(!res$OK) stop(firstMsg,res$message, call. = FALSE)
# # INTEGER is a restrictive type that deserves some explanations (not included in getTheTypes)
# if(grepl("integer",loType)){
# if(grepl("single",loType)){
# if(!is.numeric(x)) stop(firstMsg,"must be an integer (right now it is not even numeric).", call. = FALSE)
# if(!(is.integer(x) || x%%1==0)) stop(firstMsg,"must be an integer.", call. = FALSE)
# } else {
# if(!is.numeric(x)) stop(firstMsg,"must be composed of integers (right now it is not even numeric).", call. = FALSE)
# if(!(is.integer(x) || all(x%%1==0))) stop(firstMsg,"must be composed of integers.", call. = FALSE)
# }
# }
# GE: greater or equal // GT: greater than // LE: lower or equal // LT: lower than
if(grepl("ge[[:digit:]]+",loType)){
n = myExtract("ge[[:digit:]]+", loType)
if( !all(x>=n) ) stop(firstMsg,"must be greater than, or equal to, ", n,".", call. = FALSE)
}
if(grepl("gt[[:digit:]]+",loType)){
n = myExtract("gt[[:digit:]]+", loType)
if( !all(x>n) ) stop(firstMsg,"must be strictly greater than ", n,".", call. = FALSE)
}
if(grepl("le[[:digit:]]+",loType)){
n = myExtract("le[[:digit:]]+", loType)
if( !all(x<=n) ) stop(firstMsg,"must be lower than, or equal to, ",n,".", call. = FALSE)
}
if(grepl("lt[[:digit:]]+",loType)){
n = myExtract("lt[[:digit:]]+", loType)
if( !all(x<n) ) stop(firstMsg,"must be strictly lower than ", n,".", call. = FALSE)
}
#
# Specific Types Handling
#
if(grepl("match.arg",loType)){
if(ignore.case){
x = toupper(x)
newCharVec = toupper(charVec)
} else {
newCharVec = charVec
}
if( is.na(pmatch(x, newCharVec)) ){
n = length(charVec)
if(n == 1){
msg = paste0("'",charVec,"'")
} else {
msg = paste0("'", paste0(charVec[1:(n-1)], collapse="', '"), "' or '",charVec[n],"'")
}
stop(firstMsg, "must be one of:\n", msg, ".", call. = FALSE)
} else {
qui = pmatch(x, newCharVec)
return(charVec[qui])
}
}
}
.matchTypeAndSetDefault = function(myList, myDefault, myTypes, prefix){
# Cette fonction:
# i) verifie que tous les elements de la liste sont valides
# ii) mes les valeurs par defauts si elles certaines valeurs sont manquantes
# iii) Envoie des messages d'erreur si les typages ne sont pas bons
# En fait cette fonction "coerce" myList en ce qu'il faudrait etre (donne par myDefault)
# 1) check that the names of the list are valid
if(is.null(myList)) myList = list()
list_names = names(myList)
if(length(list_names)!=length(myList) || any(list_names=="")){
stop(prefix,"The elements of the list should be named.", call. = FALSE)
}
obj_names = names(myDefault)
isHere = pmatch(list_names,obj_names)
if(anyNA(isHere)){
if(sum(is.na(isHere))==1) stop(prefix, "The following argument is not defined: ",paste(list_names[is.na(isHere)],sep=", "), call. = FALSE)
else stop(prefix, "The following arguments are not defined: ",paste(list_names[is.na(isHere)],sep=", "), call. = FALSE)
}
# 2) We set the default values and run Warnings
res = list()
for(i in 1:length(obj_names)){
obj = obj_names[i]
qui = which(isHere==i) # qui vaut le numero de l'objet dans myList
type = myTypes[i] # we extract the type => to control for "match.arg" type
if(length(qui)==0){
# we set to the default if it's missing
if(type == "match.arg") {
res[[obj]] = myDefault[[i]][1]
} else {
res[[obj]] = myDefault[[i]]
}
} else {
# we check the integrity of the value
val = myList[[qui]]
if(type == "match.arg"){
# If the value is to be a match.arg => we use our controls and not
# directly the one of the function match.arg()
charVec = myDefault[[i]]
.myAlerts(val, obj, "singleCharacterMatch.arg", prefix, charVec)
val = match.arg(val, charVec)
} else {
.myAlerts(val, obj, type, prefix)
}
res[[obj]] = val
}
}
return(res)
}
.checkTheTypes = function(str, x){
# This function takes in a character string describing the types of the
# element x => it can be of several types
# types that are controlled for:
allTypes = c("numeric", "integer", "character", "logical", "list", "data.frame", "matrix", "factor")
OK = FALSE
message = c()
for(type in allTypes){
if(grepl(type, str)){
# we add the type of the control
message = c(message, type)
if(type == "numeric"){
if(!OK & is.numeric(x)){
OK = TRUE
}
} else if(type == "integer"){
if(is.numeric(x) && (is.integer(x) || all(x%%1==0))){
OK = TRUE
}
} else if(type == "character"){
if(is.character(x)){
OK = TRUE
}
} else if(type == "logical"){
if(is.logical(x)){
OK = TRUE
}
} else if(type == "list"){
if(is.list(x)){
OK = TRUE
}
} else if(type == "data.frame"){
if(is.data.frame(x)){
OK=TRUE
}
} else if(type == "matrix"){
if(is.matrix(x)){
OK = TRUE
}
} else if(type == "factor"){
if(is.factor(x)){
OK = TRUE
}
}
}
if(OK) break
}
if(length(message) == 0) OK = TRUE #ie there is no type to be searched
else if(length(message) >= 3){
n = length(message)
message = paste0("must be of type: ", paste0(message[1:(n-1)], collapse = ", "), " or ", message[n], ".")
} else {
message = paste0("must be of type: ", paste0(message, collapse = " or "), ".")
}
return(list(OK=OK, message=message))
}
.hdclassif_dim_choice <- function(ev, n, method, threshold, graph, noise.ctrl){
# Selection of the intrinsic dimension
# browser()
N <- sum(n)
prop <- n/N
K = ifelse(is.matrix(ev), nrow(ev), 1)
# browser()
if(is.matrix(ev) && K>1){
p <- ncol(ev)
if(method=="cattell"){
dev <- abs(apply(ev, 1, diff))
max_dev <- apply(dev, 2, max, na.rm=TRUE)
dev <- dev/rep(max_dev, each=p-1)
d <- apply((dev>threshold)*(1:(p-1))*t(ev[, -1]>noise.ctrl), 2, which.max)
if(graph){
op = par(mfrow=c(K*(K<=3)+2*(K==4)+3*(K>4 && K<=9)+4*(K>9), 1+floor(K/4)-1*(K==12)+1*(K==7)))
for(i in 1:K){
sub1 <- paste("Class #", i, ", d", i, "=", d[i], sep="")
Nmax <- max(which(ev[i, ]>noise.ctrl))-1
plot(dev[1:(min(d[i]+5, Nmax)), i], type="l", col="blue", main=paste("Cattell's Scree-Test\n", sub1, sep=""), ylab=paste("threshold =", threshold), xlab="Dimension", ylim=c(0, 1.05))
abline(h=threshold, lty=3)
points(d[i], dev[d[i], i], col='red')
}
par(op)
}
} else if(method=="bic"){
d <- rep(0, K)
if(graph) op = par(mfrow=c(K*(K<=3)+2*(K==4)+3*(K>4 && K<=9)+4*(K>9), 1*(1+floor(K/4)-1*(K==12)+1*(K==7))))
for (i in 1:K) {
B <- c()
Nmax <- max(which(ev[i, ]>noise.ctrl))-1
p2 <- sum(!is.na(ev[i, ]))
Bmax <- -Inf
for (kdim in 1:Nmax){
if ((d[i]!=0 & kdim>d[i]+10)) break
a <- sum(ev[i, 1:kdim])/kdim
b <- sum(ev[i, (kdim+1):p2])/(p2-kdim)
if (b<0 | a<0){
B[kdim] <- -Inf
} else {
L2 <- -1/2*(kdim*log(a)+(p2-kdim)*log(b)-2*log(prop[i])+p2*(1+1/2*log(2*pi))) * n[i]
B[kdim] <- 2*L2 - (p2+kdim*(p2-(kdim+1)/2)+1) * log(n[i])
}
if ( B[kdim]>Bmax ){
Bmax <- B[kdim]
d[i] <- kdim
}
}
if(graph){
plot(B, type='l', col=4, main=paste("class #", i, ", d=", d[i], sep=''), ylab='BIC', xlab="Dimension")
points(d[i], B[d[i]], col=2)
}
}
if(graph) par(op)
}
} else{
ev <- as.vector(ev)
p <- length(ev)
if(method=="cattell"){
dvp <- abs(diff(ev))
Nmax <- max(which(ev>noise.ctrl))-1
if (p==2) d <- 1
else d <- max(which(dvp[1:Nmax]>=threshold*max(dvp[1:Nmax])))
diff_max <- max(dvp[1:Nmax])
if(graph){
plot(dvp[1:(min(d+5, p-1))]/diff_max, type="l", col="blue", main=paste("Cattell's Scree-Test\nd=", d, sep=''), ylab=paste("threshold =", threshold, sep=' '), xlab='Dimension', ylim=c(0, 1.05))
abline(h=threshold, lty=3)
points(d, dvp[d]/diff_max, col='red')
}
} else if(method=="bic"){
d <- 0
Nmax <- max(which(ev>noise.ctrl))-1
B <- c()
Bmax <- -Inf
for (kdim in 1:Nmax){
if (d!=0 && kdim>d+10) break
a <- sum(ev[1:kdim])/kdim
b <- sum(ev[(kdim+1):p])/(p-kdim)
if (b<=0 | a<=0) B[kdim] <- -Inf
else{
L2 <- -1/2*(kdim*log(a)+(p-kdim)*log(b)+p*(1+1/2*log(2*pi)))*N
B[kdim] <- 2*L2 - (p+kdim*(p-(kdim+1)/2)+1)*log(N)
}
if ( B[kdim]>Bmax ){
Bmax <- B[kdim]
d <- kdim
}
}
if(graph){
plot(B, type='l', col=4, main=paste("BIC criterion\nd=", d, sep=''), ylab='BIC', xlab="Dimension")
points(d, B[d], col=2)
}
}
}
return(d)
}
.hdclassif_bic <- function(par, p, data=NULL){
model <- par$model
K <- par$K
d <- par$d
b <- par$b
a <- par$a
mu <- par$mu
N <- par$N
prop <- par$prop
if(length(b)==1){
#update of b to set it as variable dimension models
eps <- sum(prop*d)
n_max <- if(model%in%c("ABQD", "AJBQD")) length(par$ev) else ncol(par$ev)
b <- b*(n_max-eps)/(p-eps)
b <- rep(b, length=K)
}
if (length(a)==1) a <- matrix(a, K, max(d))
else if (length(a)==K) a <- matrix(a, K, max(d))
else if (model=='AJBQD') a <- matrix(a, K, d[1], byrow=TRUE)
if(min(a, na.rm=TRUE)<=0 | any(b<0)) return(-Inf)
if (is.null(par$loglik)){
som_a <- c()
for (i in 1:K) som_a[i] <- sum(log(a[i, 1:d[i]]))
L <- -1/2*sum(prop * (som_a + (p-d)*log(b) - 2*log(prop) + p*(1+log(2*pi))))*N
} else if (model%in%c("ABQD", "AJBQD")){
Q <- rep(list(par$Q), K)
K_pen <- matrix(0, K, N)
for (i in 1:K) {
s <- sum(log(a[i, 1:d[i]]))
X <- data-matrix(mu[i, ], N, p, byrow=TRUE)
proj <- (X%*%Q[[i]])%*%t(Q[[i]])
A <- (-proj)%*%Q[[i]]%*%sqrt(diag(1/a[i, 1:d[i]], d[i]))
B <- X-proj
K_pen[i, ] <- rowSums(A^2)+1/b[i]*rowSums(B^2)+s+(p-d[i])*log(b[i])-2*log(prop[i])+p*log(2*pi)
}
A <- -1/2*t(K_pen)
L <- sum(log(rowSums(exp(A-apply(A, 1, max))))+apply(A, 1, max))
} else L <- par$loglik[length(par$loglik)]
ro <- K*p+K-1
tot <- sum(d*(p-(d+1)/2))
D <- sum(d)
d <- d[1]
to <- d*(p-(d+1)/2)
if (model=='AKJBKQKDK') m <- ro+tot+D+K
else if (model=='AKBKQKDK') m <- ro+tot+2*K
else if (model=='ABKQKDK') m <- ro+tot+K+1
else if (model=='AKJBQKDK') m <- ro+tot+D+1
else if (model=='AKBQKDK') m <- ro+tot+K+1
else if (model=='ABQKDK') m <- ro+tot+2
bic <- -(-2*L+m*log(N))
#calcul ICL
t = par$posterior
#if(!is.null(t)){
# means we are in HDDC
Z = ((t - apply(t, 1, max))==0) + 0
icl = bic - 2*sum(Z*log(t+1e-15))
# } else {
# # Si HDDA, entropie est nulle => car classes pures
# icl = bic
#}
return(list(bic = bic, icl = icl))
}
.hdc_getComplexity = function(par, p){
model <- par$model
K <- par$K
d <- par$d
b <- par$b
a <- par$a
mu <- par$mu
N <- par$N
prop <- par$prop
ro <- K*p+K-1
tot <- sum(d*(p-(d+1)/2))
D <- sum(d)
d <- d[1]
to <- d*(p-(d+1)/2)
if (model=='AKJBKQKDK') m <- ro+tot+D+K
else if (model=='AKBKQKDK') m <- ro+tot+2*K
else if (model=='ABKQKDK') m <- ro+tot+K+1
else if (model=='AKJBQKDK') m <- ro+tot+D+1
else if (model=='AKBQKDK') m <- ro+tot+K+1
else if (model=='ABQKDK') m <- ro+tot+2
return(m)
}
.hdc_getTheModel = function(model, all2models = FALSE){
# Function used to get the models from number or names
model_in = model
if(!is.vector(model)) stop("The argument 'model' must be a vector.")
if(anyNA(model)) stop("The argument 'model' must not contain any NA.")
ModelNames <- c("AKJBKQKDK", "AKBKQKDK", "ABKQKDK", "AKJBQKDK", "AKBQKDK", "ABQKDK", "AKJBKQKD", "AKBKQKD", "ABKQKD", "AKJBQKD", "AKBQKD", "ABQKD", "AJBQD", "ABQD")
model = toupper(model)
if(length(model)==1 && model=="ALL"){
if(all2models) model <- 1:14
else return("ALL")
}
qui = which(model %in% 1:14)
model[qui] = ModelNames[as.numeric(model[qui])]
# We order the models properly
qui = which(!model%in%ModelNames)
if (length(qui)>0){
if(length(qui)==1){
msg = paste0("(e.g. ", model_in[qui], " is incorrect.)")
} else {
msg = paste0("(e.g. ", paste0(model_in[qui[1:2]], collapse=", or "), " are incorrect.)")
}
stop("Invalid model name ", msg)
}
# warning:
if(max(table(model))>1) warning("The model vector, argument 'model', is made unique (repeated values are not tolerated).")
mod_num <- c()
for(i in 1:length(model)) mod_num[i] <- which(model[i]==ModelNames)
mod_num <- sort(unique(mod_num))
model <- ModelNames[mod_num]
return(model)
}
####
#### Utilities ####
####
.addCommas = function(x) sapply(x, .addCommas_single )
.addCommas_single = function(x){
# Cette fonction ajoute des virgules pour plus de
# visibilite pour les (tres longues) valeurs de vraisemblance
if(!is.finite(x)) return(as.character(x))
s = sign(x)
x = abs(x)
decimal = x - floor(x)
if(decimal>0) dec_string = substr(decimal, 2, 4)
else dec_string = ""
entier = as.character(floor(x))
quoi = rev(strsplit(entier, "")[[1]])
n = length(quoi)
sol = c()
for(i in 1:n){
sol = c(sol, quoi[i])
if(i%%3 == 0 && i!=n) sol = c(sol, ",")
}
res = paste0(ifelse(s==-1, "-", ""), paste0(rev(sol), collapse=""), dec_string)
res
}
.repmat <- function(v,n,p){
if (p==1){M = cbind(rep(1,n)) %*% v}
else { cat('!'); M = matrix(rep(v,n),n,(length(v)*p),byrow=T)}
M
}
.diago <- function(v){
if (length(v)==1){ res = v }
else { res = diag(v)}
res
}
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.