R/mclust.R

Defines functions emXXX cdensXXX mvnXXX emXXI cdensXXI mvnXXI emXII cdensXII mvnXII emX cdensX mvnX simVVV mstepVVV meVVV estepVVV emVVV cdensVVV simVVI mstepVVI meVVI estepVVI emVVI cdensVVI simV mstepV meV simVII mstepVII meVVI meVII estepVII emVII cdensVII simVEV mstepVEV meVEV estepVEV emVEV cdensVEV estepV emV cdensV simVEI mstepVEI meVEI estepVEI emVEI cdensVEI sim nMclustParams nVarParams mvn mstep me mclustVariance estep em checkModelName bic decomp2sigma sigma2decomp sigma2decomp simEVI mstepEVI meEVI estepEVI emEVI cdensEVI simE mstepE meE simEII mstepEII meEII estepEII emEII simEEV mstepEEV meEEV estepEEV emEEV cdensEEV estepE emE cdensE simEEI mstepEEI meEEI estepEEI emEEI cdensEII cdensEEI simEEE mstepEEE meEEE estepEEE emEEE cdensEEE priorControl emControl defaultPrior mclustModelNames mclustModel print.mclustLoglik mclustLoglik mclustBICupdate pickBIC plot.mclustBIC print.summary.mclustBIC summaryMclustBICn summary.mclustBIC print.mclustBIC mclustBIC predict.Mclust logLik.Mclust plot.Mclust print.summary.Mclust summary.Mclust print.Mclust Mclust

Documented in bic cdensE cdensEEE cdensEEI cdensEEV cdensEII cdensEVI cdensV cdensVEI cdensVEV cdensVII cdensVVI cdensVVV cdensX cdensXII cdensXXI cdensXXX checkModelName decomp2sigma defaultPrior em emControl emE emEEE emEEI emEEV emEII emEVI emV emVEI emVEV emVII emVVI emVVV emX emXII emXXI emXXX estep estepE estepEEE estepEEI estepEEV estepEII estepEVI estepV estepVEI estepVEV estepVII estepVVI estepVVV logLik.Mclust Mclust mclustBIC mclustBICupdate mclustLoglik mclustModel mclustModelNames mclustVariance me meE meEEE meEEI meEEV meEII meEVI meV meVEI meVEV meVII meVVI meVVV mstep mstepE mstepEEE mstepEEI mstepEEV mstepEII mstepEVI mstepV mstepVEI mstepVEV mstepVII mstepVVI mstepVVV mvn mvnX mvnXII mvnXXI mvnXXX nMclustParams nVarParams pickBIC plot.Mclust plot.mclustBIC predict.Mclust print.Mclust print.mclustBIC print.mclustLoglik print.summary.Mclust print.summary.mclustBIC priorControl sigma2decomp sim simE simEEE simEEI simEEV simEII simEVI simV simVEI simVEV simVII simVVI simVVV summary.Mclust summary.mclustBIC summary.mclustBIC summaryMclustBICn

Mclust <- function(data, G = NULL, modelNames = NULL, prior = NULL, 
                   control = emControl(), initialization = NULL, 
                   warn = mclust.options("warn"), x = NULL, 
                   verbose = interactive(), ...) 
{
  call <- match.call()
  data <- data.matrix(data)
  d <- ncol(data)
  
  if(!is.null(x))
  {
    if(!inherits(x, "mclustBIC"))
      stop("If provided, argument x must be an object of class 'mclustBIC'")
  }

  mc <- match.call(expand.dots = TRUE)
  mc[[1]] <- as.name("mclustBIC")
  mc[[2]] <- data

  BIC <- eval(mc, parent.frame())
  # get the best model from BIC table
  G <- attr(BIC, "G")
  modelNames <- attr(BIC, "modelNames")
  Sumry <- summary(BIC, data, G = G, modelNames = modelNames)
  
  if(length(Sumry)==0) 
  {
    if(warn) 
      warning("no model(s) could be fitted. Try adjusting G and modelNames arguments")
    return()
  }
  
  if(!(length(G) == 1)) 
    { bestG <- length(tabulate(Sumry$cl))
      if(warn) 
        { if(bestG == max(G) & warn) 
             warning("optimal number of clusters occurs at max choice")
          else if(bestG == min(G) & warn) 
             warning("optimal number of clusters occurs at min choice")
        }
  }
  oldClass(Sumry) <- NULL
  
  Sumry$bic <- Sumry$bic[1]
  Sumry$icl <- icl.Mclust(Sumry)
  Sumry$hypvol <- if(is.null(attr(BIC, "Vinv"))) 
                    as.double(NA) else 1/attr(BIC, "Vinv")
  # df <- (2*Sumry$loglik - Sumry$bic)/log(Sumry$n)
  df <- if(is.null(Sumry$modelName)) NULL 
        else with(Sumry, nMclustParams(modelName, d, G, 
                                       noise = (!is.na(hypvol)),
                                       equalPro = attr(Sumry, "control")$equalPro))

  ans <- c(list(call = call, data = data, BIC = BIC, df = df), Sumry)
  orderedNames <- c("call", "data", "modelName", 
                    "n", "d", "G", 
                    "BIC", "loglik", "df", "bic", "icl",
                    "hypvol", "parameters", "z", 
                    "classification", "uncertainty")
  structure(ans[orderedNames], class = "Mclust")  
}

print.Mclust <- function(x, digits = getOption("digits"), ...)
{
  txt <- paste0("\'", class(x)[1], "\' model object: ")
  noise <- !is.null(attr(x$BIC, "Vinv"))
  if(x$G == 0 & noise)
    { txt <- paste0(txt, "single noise component") }
  else
    { txt <- paste0(txt, "(", x$model, ",", x$G, ")",
                    if(noise) " + noise component") }
  catwrap(txt)
  cat("\n")
  catwrap("\nAvailable components:\n")
  print(names(x))
  invisible(x)
}

summary.Mclust <- function(object, classification = TRUE, parameters = FALSE, ...)
{
  classification <- as.logical(classification)
  parameters <- as.logical(parameters)
  # collect info
  G  <- object$G
  noise <- if(is.na(object$hypvol)) FALSE else object$hypvol
  pro <- object$parameters$pro
  if(is.null(pro)) pro <- 1
  names(pro) <- if(noise) c(seq_len(G),0) else seq(G)
  mean <- object$parameters$mean
  if(object$d > 1)
    { sigma <- object$parameters$variance$sigma }
  else
    { sigma <- rep(object$parameters$variance$sigmasq, object$G)[1:object$G]
      names(sigma) <- names(mean) }
  if(is.null(object$density))
  {
    title <- paste("Gaussian finite mixture model fitted by EM algorithm")
    printClassification <- classification
    classification <- if(printClassification)
    {
      factor(object$classification, 
             levels = { l <- seq_len(object$G)
                        if(is.numeric(noise)) l <- c(l,0) 
                        l })
    } else NULL
  } else
  {
    title <- paste("Density estimation via Gaussian finite mixture modeling")
    printClassification <- FALSE
    classification <- NULL
  }           
  #
  obj <- list(title = title, n = object$n, d = object$d, 
              G = G, modelName = object$modelName, 
              loglik = object$loglik, df = object$df, 
              bic = object$bic, icl = object$icl,
              pro = pro, mean = mean, variance = sigma,
              noise = noise,
              prior = attr(object$BIC, "prior"), 
              printParameters = parameters,
              printClassification = printClassification,
              classification = classification)
  class(obj) <- "summary.Mclust"
  return(obj)
}

print.summary.Mclust <- function(x, digits = getOption("digits"), ...)
{
  txt <- paste(rep("-", min(nchar(x$title), getOption("width"))), collapse = "")
  catwrap(txt)
  catwrap(x$title)
  catwrap(txt)
  #
  cat("\n")
  if(x$G == 0)
  { 
    catwrap("Mclust model with only a noise component:")
  } else
  { 
    catwrap(paste0("Mclust ", x$modelName, " (", 
        mclustModelNames(x$modelName)$type, ") model with ", 
        x$G, ifelse(x$G > 1, " components", " component"),
        if(x$noise) " and a noise term", ":"))
  }
  cat("\n")
  #
  if(!is.null(x$prior))
  { 
    catwrap(paste0("Prior: ", x$prior$functionName, "(", 
                   paste(names(x$prior[-1]), x$prior[-1], sep = " = ", 
                         collapse = ", "), ")", sep = ""))
    cat("\n")
  }
  #
  tab <- data.frame("log-likelihood" = x$loglik, "n" = x$n, 
                    "df" = x$df, "BIC" = x$bic, "ICL" = x$icl, 
                    row.names = "", check.names = FALSE)
  print(tab, digits = digits)
  #
  if(x$printClassification)
  {
    cat("\nClustering table:")
    print(table(x$classification), digits = digits)
  }
  #
  if(x$printParameters)
  { 
    cat("\nMixing probabilities:\n")
    print(x$pro, digits = digits)
    cat("\nMeans:\n")
    print(x$mean, digits = digits)
    cat("\nVariances:\n")
    if(x$d > 1) 
      { for(g in 1:x$G)
           { cat("[,,", g, "]\n", sep = "")
             print(x$variance[,,g], digits = digits) }
      }
    else print(x$variance, digits = digits)
    if(x$noise)
      { cat("\nHypervolume of noise component:\n")
        cat(signif(x$noise, digits = digits), "\n") }
  }
  #
  invisible(x)
}

plot.Mclust <- function(x, 
                        what = c("BIC", "classification", "uncertainty", "density"), 
                        dimens = NULL, xlab = NULL, ylab = NULL,  
                        addEllipses = TRUE, main = FALSE,
                        ...) 
{
  
  object <- x # Argh.  Really want to use object anyway
  if(!inherits(object, "Mclust")) 
    stop("object not of class 'Mclust'")
  
  data <- object$data
  p <- ncol(data)
  if(p == 1) 
    colnames(data) <- deparse(object$call$data)
  dimens <- if(is.null(dimens)) seq(p) else dimens[dimens <= p]
  d <- length(dimens)
  
  main <- if(is.null(main) || is.character(main)) FALSE else as.logical(main)
  
  what <- match.arg(what, several.ok = TRUE)
  oldpar <- par(no.readonly = TRUE)

  plot.Mclust.bic <- function(...)
    plot.mclustBIC(object$BIC, xlab = xlab, ...)

  plot.Mclust.classification <- function(...)
  {  
    if(d == 1)
    { 
      mclust1Dplot(data = data[,dimens,drop=FALSE], 
                   what = "classification",
                   classification = object$classification,
                   z = object$z, 
                   xlab = if(is.null(xlab)) colnames(data)[dimens] else xlab, 
                   main = main, ...) 
    }
    if(d == 2) 
    { 
      pars <- object$parameters
      pars$mean <- pars$mean[dimens,,drop=FALSE]
      pars$variance$d <- length(dimens)
      pars$variance$sigma <- pars$variance$sigma[dimens,dimens,,drop=FALSE]
      mclust2Dplot(data = data[,dimens,drop=FALSE], 
                   what = "classification", 
                   classification = object$classification, 
                   parameters = if(addEllipses) pars else NULL,
                   xlab = if(is.null(xlab)) colnames(data)[dimens][1] else xlab, 
                   ylab = if(is.null(ylab)) colnames(data)[dimens][2] else ylab,
                   main = main, ...) 
    }
    if(d > 2)
    { 
      pars <- object$parameters
      pars$mean <- pars$mean[dimens,,drop=FALSE]
      pars$variance$d <- length(dimens)
      pars$variance$sigma <- pars$variance$sigma[dimens,dimens,,drop=FALSE]
      on.exit(par(oldpar))
      par(mfrow = c(d, d), 
          mar = rep(0.2/2,4), 
          oma = rep(3,4))
      for(i in seq(d))
      {
        for(j in seq(d))
        {
          if(i == j)
          {
            plot(data[, dimens[c(j, i)]],
                 type = "n", xlab = "", ylab = "", axes = FALSE)
            text(mean(par("usr")[1:2]), mean(par("usr")[3:4]),
                 labels = colnames(data[, dimens])[i],
                 cex = 1.5, adj = 0.5)
            box()
          } else 
          { 
            coordProj(data = data, 
                      dimens = dimens[c(j,i)], 
                      what = "classification", 
                      classification = object$classification,
                      parameters = object$parameters,
                      addEllipses = addEllipses,
                      main = FALSE, xaxt = "n", yaxt = "n", ...)
          }
          if(i == 1 && (!(j%%2))) axis(3)
          if(i == d && (j%%2))    axis(1)
          if(j == 1 && (!(i%%2))) axis(2)
          if(j == d && (i%%2))    axis(4)
        }
      }
    }
  }

  plot.Mclust.uncertainty <- function(...) 
  {
    pars <- object$parameters
    if(d > 1)
    {
      pars$mean <- pars$mean[dimens,,drop=FALSE]
      pars$variance$d <- length(dimens)
      pars$variance$sigma <- pars$variance$sigma[dimens,dimens,,drop=FALSE]
    }
    #
    if(p == 1 || d == 1)
    { 
      mclust1Dplot(data = data[,dimens,drop=FALSE], 
                   what = "uncertainty", 
                   parameters = pars, z = object$z, 
                   xlab = if(is.null(xlab)) colnames(data)[dimens] else xlab, 
                   main = main, ...) 
    }
    if(p == 2 || d == 2)
    { 
      mclust2Dplot(data = data[,dimens,drop=FALSE], 
                   what = "uncertainty", 
                   parameters = pars,
                   # uncertainty = object$uncertainty,
                   z = object$z,
                   classification = object$classification,
                   xlab = if(is.null(xlab)) colnames(data)[dimens][1] else xlab, 
                   ylab = if(is.null(ylab)) colnames(data)[dimens][2] else ylab,
                   addEllipses = addEllipses, main = main, ...)
    }
    if(p > 2 && d > 2)
    { 
      on.exit(par(oldpar))
      par(mfrow = c(d, d), 
          mar = rep(0,4),
          mar = rep(0.2/2,4), 
          oma = rep(3,4))
      for(i in seq(d))
      { 
        for(j in seq(d)) 
        { 
          if(i == j) 
          { 
            plot(data[, dimens[c(j, i)]], type="n",
                 xlab = "", ylab = "", axes = FALSE)
            text(mean(par("usr")[1:2]), mean(par("usr")[3:4]),
                 labels = colnames(data[,dimens])[i], 
                 cex = 1.5, adj = 0.5)
            box()
          } else 
          { 
            coordProj(data = data, 
                      what = "uncertainty", 
                      parameters = object$parameters,
                      # uncertainty = object$uncertainty,
                      z = object$z,
                      classification = object$classification,
                      dimens = dimens[c(j,i)], 
                      main = FALSE, 
                      addEllipses = addEllipses,
                      xaxt = "n", yaxt = "n", ...)
          }
          if(i == 1 && (!(j%%2))) axis(3)
          if(i == d && (j%%2))    axis(1)
          if(j == 1 && (!(i%%2))) axis(2)
          if(j == d && (i%%2))    axis(4)
        }
      }
    }
  }

  plot.Mclust.density <- function(...)
  {
    if(p == 1)
      { 
        objdens <- as.densityMclust(object)
        plotDensityMclust1(objdens, 
                           xlab = if(is.null(xlab)) colnames(data)[dimens] else xlab, 
                           main = if(main) main else NULL, ...) 
        # mclust1Dplot(data = data,
        #              parameters = object$parameters,
        #              # z = object$z, 
        #              what = "density", 
        #              xlab = if(is.null(xlab)) colnames(data)[dimens] else xlab, 
        #              main = main, ...) 
    }
    if(p == 2) 
      { surfacePlot(data = data, 
                    parameters = object$parameters,
                    what = "density", 
                    xlab = if(is.null(xlab)) colnames(data)[1] else xlab, 
                    ylab = if(is.null(ylab)) colnames(data)[2] else ylab,
                    main = main, ...) 
    }
    if(p > 2) 
      { 
        objdens <- as.densityMclust(object)
        objdens$data <- objdens$data[,dimens,drop=FALSE]
        objdens$varname <- colnames(data)[dimens]
        objdens$range <- apply(data, 2, range)
        objdens$d <- d
        objdens$parameters$mean <- objdens$parameters$mean[dimens,,drop=FALSE]
        objdens$parameters$variance$d <- d
        objdens$parameters$variance$sigma <- 
          objdens$parameters$variance$sigma[dimens,dimens,,drop=FALSE]
        # 
        if (d == 1)
          plotDensityMclust1(objdens, ...)
        else if (d == 2)
          plotDensityMclust2(objdens, ...)
        else
          plotDensityMclustd(objdens, ...)
      }
  }
  
  if(interactive() & length(what) > 1)
    { title <- "Model-based clustering plots:"
      # present menu waiting user choice
      choice <- menu(what, graphics = FALSE, title = title)
      while(choice != 0)
           { if(what[choice] == "BIC")            plot.Mclust.bic(...)
             if(what[choice] == "classification") plot.Mclust.classification(...)
             if(what[choice] == "uncertainty")    plot.Mclust.uncertainty(...)
             if(what[choice] == "density")        plot.Mclust.density(...)
             # re-present menu waiting user choice
             choice <- menu(what, graphics = FALSE, title = title)
           }
  } 
  else 
    { if(any(what == "BIC"))            plot.Mclust.bic(...)
      if(any(what == "classification")) plot.Mclust.classification(...) 
      if(any(what == "uncertainty"))    plot.Mclust.uncertainty(...) 
      if(any(what == "density"))        plot.Mclust.density(...) 
  }
    
  invisible()
}

logLik.Mclust <- function(object, ...)
{
  if(is.null(object$loglik)) 
    l <- sum(do.call("dens", c(object, logarithm = TRUE)))
  else
    l <- object$loglik
  if(is.null(object$df)) 
    {
      noise <- if(is.null(object$hypvol)) FALSE else (!is.na(object$hypvol))
      equalPro <- if(is.null(object$BIC)) FALSE else attr(object$BIC, "control")$equalPro
      df <- with(object, nMclustParams(modelName, d, G, 
                                       noise = noise, equalPro = equalPro))
  } else df <- object$df
  attr(l, "nobs") <- object$n
  attr(l, "df") <- df
  class(l) <- "logLik"
  return(l)
}

predict.Mclust <- function(object, newdata, ...)
{
  if(!inherits(object, "Mclust")) 
    stop("object not of class 'Mclust'")
  if(missing(newdata))
    { newdata <- object$data }
  newdata <- as.matrix(newdata)
  if(ncol(object$data) != ncol(newdata))
    { stop("newdata must match ncol of object data") }
  #
	object$data <- newdata
  z <- do.call("cdens", c(object, list(logarithm = TRUE)))
  pro <- object$parameters$pro
	pro <- pro/sum(pro)
  noise <- (!is.na(object$hypvol))
  z <- if(noise) cbind(z, log(object$parameters$Vinv))
       else      cbind(z) # drop redundant attributes
  z <- sweep(z, MARGIN = 2, FUN = "+", STATS = log(pro))
  z <- sweep(z, MARGIN = 1, FUN = "-", STATS = apply(z, 1, logsumexp))
  z <- exp(z)
  cl <- c(seq(object$G), if(noise) 0)
  colnames(z) <- cl
  cl <- cl[apply(z, 1, which.max)]
  out <- list(classification = cl, z = z)
  return(out) 
}


mclustBIC <- function(data, G = NULL, modelNames = NULL, 
                      prior = NULL, control = emControl(), 
                      initialization = list(hcPairs = NULL, 
                                            subset = NULL, 
                                            noise = NULL),  
                      Vinv = NULL, warn = mclust.options("warn"), 
                      x = NULL, verbose = interactive(), ...)
{
  dimData <- dim(data)
  oneD <- (is.null(dimData) || length(dimData[dimData > 1]) == 1)
  if(!oneD && length(dimData) != 2)
    stop("data must be a vector or a matrix")
  if(oneD) 
    { data <- drop(as.matrix(data))
      n <- length(data)
      d <- 1
  } else {
      data <- as.matrix(data)
      n <- nrow(data)
      d <- ncol(data)
  }

  if(is.null(x)) 
    { if(is.null(modelNames)) 
        { if(d == 1) 
            { modelNames <- c("E", "V") 
          } else {
              modelNames <- mclust.options("emModelNames")
              if(n <= d) 
                { # select only spherical and diagonal models
                  m <- match(modelNames, c("EII", "VII", "EEI", 
                                           "VEI", "EVI", "VVI"),
                             nomatch = 0)
                  modelNames <- modelNames[m]
                }
          }
      }
    
      if(!is.null(prior))
        { # remove models not available with prior
          modelNames <- setdiff(modelNames, c("EVE","VEE","VVE","EVV"))
      }
    
      if(is.null(G)) 
        { G <- if (is.null(initialization$noise)) 1:9 else 0:9 }
      else 
        { G <- sort(as.integer(unique(G))) }
    
      if(is.null(initialization$noise)) 
        { if (any(G > n)) G <- G[G <= n] }
      else 
        {
          noise <- initialization$noise
          if(is.logical(noise)) noise <- which(noise)
          if(any(match(noise, 1:n, nomatch = 0) == 0))
            stop("numeric or logical vector for noise must correspond to row indexes of data")
          initialization$noise <- noise
          nnoise <- length(noise)
          if(any(G > (n-nnoise))) G <- G[G <= n-nnoise]
        }
    
      if(!is.null(initialization$subset))
        { subset <- initialization$subset
          if(is.logical(subset)) subset <- which(subset)
          initialization$subset <- subset
          if(any(G > n)) G <- G[G <= n]
      }
      Gall <- G
      Mall <- modelNames
    }
  else 
    {
      if(!missing(prior) || !missing(control) || 
         !missing(initialization) || !missing(Vinv))
         stop("only G and modelNames may be specified as arguments when x is supplied")
      prior <- attr(x,"prior") 
      control <- attr(x,"control")
      initialization <- attr(x,"initialization")
      Vinv <- attr(x,"Vinv")
      warn <- attr(x,"warn")
      Glabels <- dimnames(x)[[1]]
      Mlabels <- dimnames(x)[[2]]
      if(is.null(G)) G <- Glabels
      if(is.null(modelNames)) modelNames <- Mlabels
      Gmatch <- match(as.character(G), Glabels, nomatch = 0)
      Mmatch <- match(modelNames, Mlabels, nomatch = 0)
      if(all(Gmatch) && all(Mmatch)) 
        { out <- x[as.character(G),modelNames,drop=FALSE]
          mostattributes(out) <- attributes(x)
          attr(out, "dim") <- c(length(G), length(modelNames))
          attr(out, "dimnames") <- list(G, modelNames)
          attr(out, "G") <- as.numeric(G)
          attr(out, "modelNames") <- modelNames
          attr(out, "returnCodes") <- 
          attr(x, "returnCodes")[as.character(G),modelNames,drop=FALSE]
          return(out)
      }
      Gall <- sort(as.numeric(unique(c(as.character(G), Glabels))))
      Mall <- unique(c(modelNames, Mlabels))
  }
  
  if(any(as.logical(as.numeric(G))) < 0) 
    { if(is.null(initialization$noise)) 
        { stop("G must be positive") }
    else {
          stop("G must be nonnegative")
    }
  }
  
  if(d == 1 && any(nchar(modelNames) > 1)) 
    { Emodel <- any(sapply(modelNames, function(x)
                    charmatch("E", x, nomatch = 0)[1]) == 1)
      Vmodel <- any(sapply(modelNames, function(x)
                    charmatch("V", x, nomatch = 0)[1]) == 1)
      modelNames <- c("E", "V")[c(Emodel, Vmodel)]
  }
  
  # set subset for initialization when subset is not, no hcPairs is provided, and
  # data size is larger than the value specified in mclust.options()
  if(is.null(initialization$subset) & 
     is.null(initialization$hcPairs) & 
     n > mclust.options("subset"))
  { 
    initialization$subset <- sample(seq.int(n), 
                                    size = mclust.options("subset"),
                                    replace = FALSE) 
  }
  
  l <- length(Gall)
  m <- length(Mall)
  if(verbose) 
    { cat("fitting ...\n")
      flush.console()
      pbar <- txtProgressBar(min = 0, max = l*m+1, style = 3)
      on.exit(close(pbar))
      ipbar <- 0
  }
  EMPTY <- -.Machine$double.xmax
  BIC <- RET <- matrix(EMPTY, nrow = l, ncol = m, 
                       dimnames = list(as.character(Gall), as.character(Mall)))
  if(!is.null(x)) 
    { BIC[dimnames(x)[[1]],dimnames(x)[[2]]] <- x
      RET[dimnames(x)[[1]],dimnames(x)[[2]]] <- attr(x, "returnCodes")
      BIC <- BIC[as.character(G),modelNames,drop=FALSE]
      RET <- RET[as.character(G),modelNames,drop=FALSE]
  }
  G <- as.numeric(G)
  Glabels <- as.character(G)
  Gout <- G
  if(is.null(initialization$noise)) 
  {
    ## standard case ----
    if (G[1] == 1) 
    {
      for(mdl in modelNames[BIC["1",] == EMPTY]) 
      {
        out <- mvn(modelName = mdl, data = data, prior = prior)
        BIC["1", mdl] <- bic(modelName = mdl, loglik = out$loglik, 
                             n = n, d = d, G = 1, equalPro = FALSE)
        RET["1", mdl] <- attr(out, "returnCode")
        if(verbose) 
          { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
      }
      if (l == 1) 
      {
        BIC[BIC == EMPTY] <- NA
        if(verbose) 
          { ipbar <- l*m+1; setTxtProgressBar(pbar, ipbar) }
        return(structure(BIC, G = G, modelNames = modelNames, prior = prior, 
                         control = control, initialization = initialization, 
                         warn = warn, n = n, d = d, oneD = oneD,
                         returnCodes = RET, class =  "mclustBIC"))
      }
      G <- G[-1]
      Glabels <- Glabels[-1]
    }
    if (is.null(initialization$subset)) 
    {
      ## all data in initial hierarchical clustering phase (no subset) ----
      if (is.null(initialization$hcPairs)) 
      { 
        if (d != 1) 
        {
          if (n > d) 
          {
            hcPairs <- hc(data = data, 
                          modelName = mclust.options("hcModelName"),
                          use = mclust.options("hcUse"))
          } else 
          {
            hcPairs <- hc(data = data, 
                          modelName = "EII",
                          use = mclust.options("hcUse"))
          } 
        }
        else {
          hcPairs <- NULL 
          # hcPairs <- hc(data = data, modelName = "E")
        }
      }
      else hcPairs <- initialization$hcPairs
      if (d > 1 || !is.null(hcPairs))  clss <- hclass(hcPairs, G)
      for (g in Glabels) 
      {
        if (d > 1 || !is.null(hcPairs)) 
        {
          cl <- clss[,g]
        } else 
        {
          cl <- qclass(data, as.numeric(g))
        }
        if(verbose) 
          { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        z <- unmap(cl, groups = 1:max(cl))
        if(any(apply( z, 2, max) == 0) & warn) 
        { #  missing groups
          if(warn) warning("there are missing groups")
          small <- sqrt(.Machine$double.neg.eps)
          z[z < small] <- small
          z <-  t(apply( z, 1, function(x) x/sum(x)))
        }
        for(modelName in na.omit(modelNames[BIC[g,] == EMPTY])) 
        {  
          out <- me(data = data, modelName = modelName, z = z,
                    prior = prior, control = control, warn = warn)
          BIC[g, modelName] <- bic(modelName = modelName, 
                                   loglik = out$loglik,
                                   n = n, d = d, G = as.numeric(g), 
                                   equalPro = control$equalPro)
          RET[g, modelName] <- attr(out, "returnCode")
          if(verbose) 
            { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        }
      }
    } else 
    {
      ## initial hierarchical clustering phase on a subset ----
      subset <- initialization$subset
      if (is.null(initialization$hcPairs)) {
        if (d != 1) {
          if (n > d) {
            hcPairs <- hc(data = data[subset,],
                          modelName = mclust.options("hcModelName"), 
                          use = mclust.options("hcUse"))
          }
          else {
            hcPairs <- hc(data = data[subset,],
                          modelName = "EII", 
                          use = mclust.options("hcUse"))
          }
        }
        else {
            hcPairs <- NULL
            # hcPairs <- hc(modelName = "E", data = data[subset])
        }
      }
      else hcPairs <- initialization$hcPairs
      if (d > 1 || !is.null(hcPairs)) clss <- hclass(hcPairs, G)
      for (g in Glabels) {
        if (d > 1 || !is.null(hcPairs)) {
          cl <- clss[, g]
        }
        else {
          cl <- qclass(data[subset], as.numeric(g))
        }
        if(verbose) 
          { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        z <- unmap(cl, groups = 1:max(cl))
        if(any(apply( z, 2, max) == 0) & warn) 
          { #  missing groups
            if(warn) warning("there are missing groups")
            small <- sqrt(.Machine$double.neg.eps)
            z[z < small] <- small
            z <-  t(apply( z, 1, function(x) x/sum(x)))
        }
        for (modelName in modelNames[!is.na(BIC[g,])]) {
          ms <- mstep(data = as.matrix(data)[initialization$subset,],
	              modelName = modelName, z = z, 
                      prior = prior, control = control, warn = warn)
          #
          #  ctrl <- control
          #  ctrl$itmax[1] <- 1
          #  ms <- me( data = as.matrix(data)[initialization$subset,  ],
	  #      modelName = modelName, z = z, prior = prior, control = ctrl)
          #
          es <- do.call("estep", c(list(data = data, warn = warn), ms))
          out <- me(data = data, modelName = modelName, z = es$z, 
                    prior = prior, control = control, warn = warn)
          BIC[g, modelName] <- bic(modelName = modelName, 
                                   loglik = out$loglik,
                                   n = n, d = d, G = as.numeric(g), 
                                   equalPro = control$equalPro)
          RET[g, modelName] <- attr(out, "returnCode")
          if(verbose) 
            { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        }
      }
    }
  }
  else 
  { 
    ## noise case ----
    noise <- initialization$noise
    if (is.null(Vinv) || Vinv <= 0)
      Vinv <- hypvol(data, reciprocal = TRUE)
    
    if (is.null(initialization$subset)) 
    {
      ## all data in initial hierarchical clustering phase (no subset) ----
      if(nnoise == n) 
        stop("All observations cannot be initialised as noise!")
      if (!G[1]) 
      {
        hood <- n * log(Vinv)
        BIC["0",] <- 2 * hood - log(n)
        if (l == 1) 
        { return(structure(BIC, G = G, modelNames = modelNames, 
                           prior = prior, control = control, 
                           initialization = list(hcPairs = hcPairs,
                                                 noise = initialization$noise), 
                           warn = warn, n = n, d = d, oneD = oneD,
                           returnCodes = RET, class =  "mclustBIC"))
        }
        G <- G[-1]
        Glabels <- Glabels[-1]
      }
      if (is.null(initialization$hcPairs)) 
      {
        if (d != 1) {
          if (n > d) {
            hcPairs <- hc(data = data[-noise,],
                          modelName = mclust.options("hcModelName"), 
                          use = mclust.options("hcUse"))
          }
          else {
            hcPairs <- hc(data = data[-noise,],
                          modelName = "EII",
                          use = mclust.options("hcUse"))
          }
        }
        else {
          hcPairs <- NULL 
          #    hcPairs <- hc(modelName = "E", data = data[-noise])
        }
      }
      else hcPairs <- initialization$hcPairs
      if (d > 1 || !is.null(hcPairs)) clss <- hclass(hcPairs, G)
      if(verbose) 
        { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }

      z <- matrix(0, n, max(G) + 1)
      for (g in Glabels) 
      {
        z[] <- 0
        k <- as.numeric(g)
        if(d > 1 || !is.null(hcPairs)) { 
          cl <- clss[,g]
        }
        else {
          cl <- qclass(data[-noise], k = k)
        }
        z[-noise,1:k] <- unmap(cl, groups = 1:max(cl))
        if(any(apply(z[-noise,1:k,drop=FALSE], 2, max) == 0) & warn) 
        { # missing groups
          if(warn) warning("there are missing groups")
          # todo: should be pmax(...) qui sotto??
          z[-noise,1:k] <- max(z[-noise,1:k], sqrt(.Machine$double.neg.eps))
          # todo: should be t(...) qui sotto??
          z[-noise,1:k] <- apply(z[-noise,1:k,drop=FALSE], 1, 
                                 function(z) z/sum(z))
        }
        z[noise, k+1] <- 1
        K <- 1:(k+1) 
        for (modelName in na.omit(modelNames[BIC[g,] == EMPTY])) 
        {  
          out <- me(data = data, modelName = modelName, z = z[, K], 
                    prior = prior, Vinv = Vinv, 
                    control = control, warn = warn)
          BIC[g, modelName] <- bic(modelName = modelName, 
                                   loglik = out$loglik, 
                                   n = n, d = d, G = k, 
                                   noise = TRUE, 
                                   equalPro = control$equalPro)
          RET[g, modelName] <- attr(out, "returnCode")
          if(verbose) 
            { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        }
      }
    } 
    else
    { 
      ## initial hierarchical clustering phase on a subset ----
      subset <- initialization$subset
      subset <- setdiff(subset, noise) # remove from subset noise obs
      initialization$subset <- subset
      if(length(subset) == 0) 
        stop("No observations in the initial subset after removing the noise!")
      
      if (!G[1]) 
      {
        hood <- n * log(Vinv)
        BIC["0",] <- 2 * hood - log(n)
        if (l == 1) 
        { return(structure(BIC, G = G, modelNames = modelNames, 
                           prior = prior, control = control, 
                           initialization = list(hcPairs = hcPairs, 
                                                 subset = initialization$subset), 
                           warn = warn, n = n, d = d, oneD = oneD,
                           returnCodes = RET, class =  "mclustBIC"))
        }
        G <- G[-1]
        Glabels <- Glabels[-1]
      }
      
      if (is.null(initialization$hcPairs)) 
      {
        if (d != 1) {
          if (n > d) {
            hcPairs <- hc(data = data[subset,],
                          modelName = mclust.options("hcModelName"), 
                          use = mclust.options("hcUse"))
          }
          else {
            hcPairs <- hc(data = data[subset,],
                          modelName = "EII", 
                          use = mclust.options("hcUse"))
          }
        }
        else {
            hcPairs <- NULL
            # hcPairs <- hc(modelName = "E",  data = data[subset])
        }
      }
      else hcPairs <- initialization$hcPairs
      if (d > 1 || !is.null(hcPairs)) clss <- hclass(hcPairs, G)
      if(verbose) 
        { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
      for (g in Glabels) 
      {
        k <- as.numeric(g)
        if (d > 1 || !is.null(hcPairs)) {
          cl <- clss[, g]
        }
        else {
          cl <- qclass(data[subset], k = k)
        }
        z <- unmap(cl, groups = 1:max(cl))
        if(any(apply(z, 2, max) == 0) & warn)
          { #  missing groups
            if(warn) warning("there are missing groups")
            small <- sqrt(.Machine$double.neg.eps)
            z[z < small] <- small
            z <- t(apply( z, 1, function(x) x/sum(x)))
        }
        for (modelName in na.omit(modelNames[BIC[g,] == EMPTY]))
        {
          ms <- mstep(data = as.matrix(data)[subset,],
                      modelName = modelName, z = z, 
                      prior = prior, control = control, warn = warn)
          es <- do.call("estep", c(list(data = data, warn = warn), ms))
          
          if(is.na(es$loglik))
          {  BIC[g, modelName] <- NA
             RET[g, modelName] <- attr(es, "returnCode") }
          else 
          {
            es$z <- cbind(es$z, 0)
            es$z[noise,] <- matrix(c(rep(0,k),1), byrow = TRUE,
                                   nrow = length(noise), ncol = k+1)
            out <- me(data = data, modelName = modelName, z = es$z, 
                      prior = prior, Vinv = Vinv, 
                      control = control, warn = warn)
            BIC[g, modelName] <- bic(modelName = modelName, 
                                     loglik = out$loglik,
                                     n = n, d = d, G = k, 
                                     noise = TRUE,
                                     equalPro = control$equalPro)
            RET[g, modelName] <- attr(out, "returnCode")
          }
          if(verbose) 
            { ipbar <- ipbar+1; setTxtProgressBar(pbar, ipbar) }
        }
      }
    }
  }
  
  if(verbose) 
    { ipbar <- l*m+1; setTxtProgressBar(pbar, ipbar) }

  if(!is.null(prior) & any(is.na(BIC)))
    warning("The presence of BIC values equal to NA is likely due to one or more of the mixture proportions being estimated as zero, so that the model estimated reduces to one with a smaller number of components.")  
  
  structure(BIC, G = Gout, modelNames = modelNames, 
            prior = prior, Vinv = Vinv, control = control, 
            initialization = list(hcPairs = hcPairs, 
                                  subset = initialization$subset,
                                  noise = initialization$noise), 
            warn = warn, n = n, d = d, oneD = oneD,
            criterion = "BIC", returnCodes = RET, 
            class = "mclustBIC")
}

print.mclustBIC <- function(x, pick = 3, ...)
{
  subset <- !is.null(attr(x, "subset"))
  oldClass(x) <- attr(x, "args") <- NULL
  attr(x, "criterion") <- NULL
  attr(x, "control") <- attr(x, "initialization") <- NULL
  attr(x, "oneD") <- attr(x, "warn") <- attr(x, "Vinv") <- NULL
  attr(x, "prior") <- attr(x, "G") <- attr(x, "modelNames") <- NULL
  ret <- attr(x, "returnCodes") == -3
  n <- attr(x, "n")
  d <- attr(x, "d")
  attr(x, "returnCodes") <- attr(x, "n") <- attr(x, "d") <- NULL
  
  catwrap("Bayesian Information Criterion (BIC):")
  NextMethod("print")
  cat("\n")
  catwrap(paste("Top", pick, "models based on the BIC criterion:"))
  print(pickBIC(x, pick), ...)
  invisible()
}

summary.mclustBIC <- function(object, data, G, modelNames, ...)
{
  mc <- match.call(expand.dots = FALSE)
  if(missing(data)) 
    { if(!missing(G)) 
        object <- object[rownames(object) %in% G,,drop=FALSE]
      if(!missing(modelNames)) 
        object <- object[,colnames(object) %in% modelNames,drop=FALSE]
      ans <- pickBIC(object, ...)
      class(ans) <- "summary.mclustBIC" 
  } else
    { if(is.null(attr(object,"initialization")$noise)) 
        { mc[[1]] <- as.name("summaryMclustBIC") }
      else 
      { mc[[1]] <- as.name("summaryMclustBICn") }
      warn <- attr(object, "warn")
      ans <- eval(mc, parent.frame())
      if(length(ans) == 0) return(ans) 
      Glabels <- dimnames(object)[[1]]
      if(length(Glabels) != 1 && (!missing(G) && length(G) > 1)) 
        { Grange <- range(as.numeric(Glabels))
          if(match(ans$G, Grange, nomatch = 0) & warn)
            warning("best model occurs at the min or max of number of components considered!")
      }
  }
  ans
}

summaryMclustBIC <- function (object, data, G = NULL, modelNames = NULL, ...) 
{
  dimData <- dim(data)
  oneD <- (is.null(dimData) || length(dimData[dimData > 1]) == 1)
  if (!oneD && length(dimData) != 2) 
    stop("data must be a vector or a matrix")
  if (oneD) {
    data <- drop(as.matrix(data))
    n <- length(data)
    d <- 1
  }
  else {
    data <- as.matrix(data)
    n <- nrow(data)
    d <- ncol(data)
  }
  initialization <- attr(object, "initialization")
  hcPairs <- initialization$hcPairs
  subset <- initialization$subset
  prior <- attr(object, "prior")
  control <- attr(object, "control")
  warn <- attr(object, "warn")
  oldClass(object) <- NULL
  attr(object, "prior") <- attr(object, "warn") <- NULL
  attr(object, "modelNames") <- attr(object, "oneD") <- NULL
  attr(object, "initialization") <- attr(object, "control") <- NULL
  d <- if (is.null(dim(data))) 1 else ncol(data)
  if(is.null(G)) 
    G <- dimnames(object)[[1]]
  if(is.null(modelNames)) 
    modelNames <- dimnames(object)[[2]]
  bestBICs <- pickBIC(object[as.character(G), modelNames, drop = FALSE], k = 3)
  if(all(is.na(bestBICs))) 
  { return(structure(list(), bestBICvalues = bestBICs, prior = prior, 
                     control = control, initialization = initialization, 
                     class = "summary.mclustBIC")) 
  }
  temp <- unlist(strsplit(names(bestBICs)[1], ","))
  bestModel <- temp[1]
  G <- as.numeric(temp[2])
  if(G == 1) 
  { 
    out <- mvn(modelName = bestModel, data = data, prior = prior)
    ans <- c(list(bic = bestBICs, 
                  z = unmap(rep(1,n)),
                  classification = rep(1, n), 
                  uncertainty = rep(0, n)), 
             out)
    orderedNames <- c("modelName", "n", "d", "G", "bic", "loglik", 
                      "parameters", "z", "classification", "uncertainty")
    return(structure(ans[orderedNames], bestBICvalues = bestBICs, 
                     prior = prior, control = control, 
                     initialization = initialization, 
                     class = "summary.mclustBIC"))
  }
  
  if(is.null(subset)) 
  {
    if(d > 1 || !is.null(hcPairs))
      { z <- unmap(hclass(hcPairs, G)) }
    else 
      { z <- unmap(qclass(data, G), groups = 1:G) }
    out <- me(data = data, modelName = bestModel, z = z, 
              prior = prior, control = control, warn = warn)
    if(sum((out$parameters$pro - colMeans(out$z))^2) > 
       sqrt(.Machine$double.eps))
      { # perform extra M-step and update parameters
        ms <- mstep(data = data, modelName = bestModel, z = out$z,
	            prior = prior, warn = warn)
        if(attr(ms, "returnCode") == 0)
          out$parameters <- ms$parameters
    }
  }
  else 
  {
    if(d > 1 || !is.null(hcPairs)) 
      { z <- unmap(hclass(hcPairs, G)) }
    else 
      { z <- unmap(qclass(data[subset], G)) }
    ms <- mstep(data = as.matrix(data)[subset,], modelName = bestModel,
                prior = prior, z = z, control = control, warn = warn)
    es <- do.call("estep", c(list(data = data), ms))
    out <- me(data = data, modelName = bestModel, z = es$z, 
              prior = prior, control = control, warn = warn)
    # perform extra M-step and update parameters
    ms <- mstep(data = data, modelName = bestModel, 
                z = out$z, prior = prior, warn = warn)
    if(attr(ms, "returnCode") == 0)
      out$parameters <- ms$parameters
  }
  obsNames <- if (is.null(dim(data))) names(data) else dimnames(data)[[1]]
  classification <- map(out$z, warn = warn)
  uncertainty <- 1 - apply(out$z, 1, max)
  names(classification) <- names(uncertainty) <- obsNames
  
  ans <- c(list(bic = bic(bestModel, out$loglik, out$n, out$d, out$G, 
                          noise = FALSE, equalPro = control$equalPro),
                # bic = as.vector(bestBICs[1]), 
                classification = classification, 
                uncertainty = uncertainty), 
           out)
  orderedNames <- c("modelName", "n", "d", "G", "bic", "loglik", 
                    "parameters", "z", "classification", "uncertainty")
  structure(ans[orderedNames], bestBICvalues = bestBICs, prior = prior, 
            control = control, initialization = initialization, 
            class = "summary.mclustBIC")
}

summaryMclustBICn <- function(object, data, G = NULL, modelNames = NULL, ...)
{
  dimData <- dim(data)
  oneD <- is.null(dimData) || length(dimData[dimData > 1]) == 1
  if(!oneD && length(dimData) != 2)
    stop("data must be a vector or a matrix")
  if(oneD) {
    data <- drop(as.matrix(data))
    n <- length(data)
    d <- 1
  }
  else {
    data <- as.matrix(data)
    n <- nrow(data)
    d <- ncol(data)
  }
  initialization <- attr(object, "initialization")
  hcPairs <- initialization$hcPairs
  subset <- initialization$subset
  noise <- initialization$noise
  if(is.logical(noise)) 
    noise <- which(noise)
  prior <- attr(object, "prior")
  control <- attr(object, "control")
  warn <- attr(object, "warn")
  Vinv <- attr(object, "Vinv")
  oldClass(object) <- NULL
  attr(object, "control") <- attr(object, "initialization") <- NULL
  attr(object, "prior") <- attr(object, "Vinv") <- NULL
  attr(object, "warn") <- NULL
  ##
  if (is.null(G)) G <- dimnames(object)[[1]]
  if (is.null(modelNames)) modelNames <- dimnames(object)[[2]]
  bestBICs <- pickBIC(object[as.character(G), modelNames, drop = FALSE], k = 3)
  if(all(is.na(bestBICs))) 
  { return(structure(list(), bestBICvalues = bestBICs, prior = prior, 
                     control = control, initialization = initialization, 
                     class = "summary.mclustBIC"))
  }
  temp <- unlist(strsplit(names(bestBICs)[1], ","))
  bestModel <- temp[1] 
  G <- as.numeric(temp[2])
  if(G == 0) 
  { ans <- list(bic = bestBICs[1], 
                z = unmap(rep(0,n)),
                classification = rep(0, n), 
                uncertainty = rep(0, n), 
                n = n, d = ncol(data), 
                modelName = bestModel, G = 0, 
                loglik = n * log(Vinv), 
                Vinv = Vinv,
                parameters = NULL)
  orderedNames <- c("modelName", "n", "d", "G", "bic", "loglik", "Vinv", 
                    "parameters", "z", "classification", "uncertainty")
  return(structure(ans[orderedNames], bestBICvalues = bestBICs, 
                   prior = prior, control = control, 
                   initialization = initialization, 
                   class = "summary.mclustBIC"))
  }
  G1 <- G + 1

  if(is.null(subset)) 
  {
    z <- matrix(0, n, G1)
    if(d > 1 || !is.null(hcPairs))
    { z[-noise, 1:G] <- unmap(hclass(hcPairs, G)) }
    else 
    { z[-noise, 1:G] <- unmap(qclass(data[-noise], G)) }
    z[noise, G1] <- 1
    out <- me(data = data, modelName = bestModel, z = z, 
              prior = prior, Vinv = Vinv,
              control = control, warn = warn)
  }
  else
  { subset <- setdiff(subset, noise) # set subset among those obs not noise
  if(d > 1 || !is.null(hcPairs))
  { z <- unmap(hclass(hcPairs, G)) }
  else 
  { z <- unmap(qclass(data[subset], G)) }
  ms <- mstep(data = as.matrix(data)[subset,], modelName = bestModel, z = z, 
              prior = prior, control = control, warn = warn)
  es <- do.call("estep", c(list(data = data, warn = warn), ms))
  es$z <- cbind(es$z, 0)
  es$z[noise,] <- matrix(c(rep(0,G),1), byrow = TRUE,
                         nrow = length(noise), ncol = G+1)
  out <- me(data = data, modelName = bestModel, z = es$z,
            prior = prior, Vinv = Vinv, 
            control = control, warn = warn)
  }
  
  obsNames <- if(is.null(dim(data))) 
    names(data) else dimnames(data)[[1]]
  classification <- map(out$z, warn = warn)
  classification[classification == G1] <- 0
  uncertainty <- 1 - apply(out$z, 1, max)
  names(classification) <- names(uncertainty) <- obsNames
  ans <- c(list(bic = as.vector(bestBICs[1]), classification = classification, 
                uncertainty = uncertainty, Vinv = Vinv), out)
  orderedNames <- c("modelName", "n", "d", "G", "bic", "loglik", "parameters", 
                    "Vinv", "z", "classification", "uncertainty")
  structure(ans[orderedNames], 
            bestBICvalues = bestBICs, 
            prior = prior, control = control, 
            initialization = initialization, 
            class = "summary.mclustBIC")
}

print.summary.mclustBIC <- function(x, digits = getOption("digits"), ...)
{
  if("classification" %in% names(x))
  { 
    bic <- attr(x,"bestBICvalues")
    l <- length(bic)
    if(l == 1) 
    { 
      cat("BIC value:\n")
      print(bic, digits = digits)
    } else 
    {
      cat("Best BIC values:\n")
      bic <- drop(as.matrix(bic))
      bic <- rbind(BIC = bic, "BIC diff" = bic - max(bic))
      print(bic, digits = digits)
    }
    cat("\n")
    catwrap(paste0("Classification table for model (", 
                   if(l == 1) names(bic)[1] else colnames(bic)[1],
                   "):"))
    print(table(x$classification), digits = digits, ...)
  } else 
  {
    cat("Best BIC values:\n")
    x <- if(length(x) == 0) attr(x,"bestBICvalues") else drop(as.matrix(x))
    x <- rbind(BIC = x, "BIC diff" = x - max(x))
    print(x, digits = digits)
  }
  invisible()
}


plot.mclustBIC <- function(x, G = NULL, modelNames = NULL, 
                           symbols = NULL, colors = NULL, 
                           xlab = NULL, ylab = "BIC",
                           legendArgs = list(x = "bottomright", ncol = 2, cex = 1, inset = 0.01), 
                           ...)
{
  args <- list(...)
  if(is.null(xlab)) xlab <- "Number of components"
  subset <- !is.null(attr(x, "initialization")$subset)
  noise <- !is.null(attr(x, "initialization")$noise)
  ret <- attr(x, "returnCodes") == -3
  # legendArgsDefault <- list(x = "bottomright", ncol = 2, cex = 1, inset = 0.01)
  legendArgsDefault <- eval(formals(plot.mclustBIC)$legendArgs)
  legendArgs <- append(as.list(legendArgs), legendArgsDefault)
  legendArgs <- legendArgs[!duplicated(names(legendArgs))]

  n <- ncol(x)
  dnx <- dimnames(x)
  x <- matrix(as.vector(x), ncol = n)
  dimnames(x) <- dnx
  if(is.null(modelNames))
    modelNames <- dimnames(x)[[2]]
  if(is.null(G))
    G <- as.numeric(dimnames(x)[[1]])
  # BIC <- x[as.character(G), modelNames, drop = FALSE]
  # X <- is.na(BIC)
  # nrowBIC <- nrow(BIC)
  # ncolBIC <- ncol(BIC)
  if(is.null(symbols)) 
  { colNames <- dimnames(x)[[2]]
    m <- length(modelNames)
    if(is.null(colNames)) 
    { symbols <- if(m > 9) LETTERS[1:m] else as.character(1:m)
      names(symbols) <- modelNames
    }
    else 
      { symbols <- mclust.options("bicPlotSymbols")[modelNames] }
  }
  if(is.null(colors)) 
  { colNames <- dimnames(x)[[2]]
    if(is.null(colNames)) 
    { colors <- 1:m
      names(colors) <- modelNames
    }
    else { 
      # colors <- mclust.options("bicPlotColors")[modelNames] 
      colors <- mclust.options("bicPlotColors") 
      if(!is.null(names(colors)) & 
         !any(names(colors) == ""))
        colors <- colors[modelNames]
    }
  }
  x <- x[,modelNames, drop = FALSE]
  ylim <- if(is.null(args$ylim)) 
          range(as.vector(x[!is.na(x)])) else args$ylim
  
  matplot(as.numeric(dnx[[1]]), x, type = "b", 
          xaxt = "n", xlim = range(G), ylim = ylim,
          pch = symbols, col = colors, lty = 1,
          xlab = xlab, ylab = ylab, main = "")
  axis(side = 1, at = as.numeric(dnx[[1]]))
  if(!is.null(legendArgs))
    { do.call("legend", c(list(legend = modelNames, col = colors, pch = symbols),
              legendArgs)) }
  invisible(symbols)
}

pickBIC <- function(x, k = 3, ...)
{
  if(!is.matrix(x))
  { warning("sorry, the pickBIC function cannot be applied to the provided argument!")
    return() }
  Glabels <- dimnames(x)[[1]]
  modelNames <- dimnames(x)[[2]]
  mis <- is.na(x)
  if(all(mis) & mclust.options("warn")) 
  { warning("none of the selected models could be fitted")
    return(rep(NA,k))
  }
  x[mis] <-  - .Machine$double.xmax
  x <- data.frame(as.vector(x), Glabels[as.vector(row(x))], 
                  modelNames[as.vector(col(x))])
  # x <- x[rev(order(x[,1])),] 
  # order by including first simpler models if ties are present
  x <- x[order(-x[, 1], x[,2], x[,3]),]
  namesx <- apply(x[,-1,drop = FALSE], 1, function(z) 
    paste(as.character(z[2]), as.character(z[1]), sep = ","))
  k <- min(k, nrow(x))
  x <- x[1:k,1]
  x[x ==  - .Machine$double.xmax] <- NA
  namesx <- namesx[1:k]
  namesx[is.na(x)] <- " "
  names(x) <- namesx
  x
}

mclustBICupdate <- function(BIC, ...)
{
  args  <- list(...)
  nargs <- length(args)
  BIC1 <- BIC
  if(length(args) > 1)
  { # recursively call the function when multiple arguments
    BIC2 <- mclustBICupdate(args[[1]], args[[-1]])
  } else {
    BIC2 <- args[[1]] }
  if(is.null(BIC1)) return(BIC2)
  if(is.null(BIC2)) return(BIC1)
  
  stopifnot(inherits(BIC1, c("mclustBIC", "mclustSBIC", "mclustICL")) & 
            inherits(BIC2, c("mclustBIC", "mclustSBIC", "mclustICL")))
  stopifnot(all.equal(attributes(BIC1)[c("n", "d")],
                      attributes(BIC2)[c("n", "d")]))

  G <- unique(c(rownames(BIC1), rownames(BIC2)))
  modelNames <- unique(c(colnames(BIC1), colnames(BIC2)))
  BIC <- matrix(as.double(NA), nrow = length(G), ncol = length(modelNames),
                dimnames = list(G, modelNames))
  BIC[rownames(BIC1),colnames(BIC1)] <- BIC1[rownames(BIC1),colnames(BIC1)]
  BIC[rownames(BIC2),colnames(BIC2)] <- BIC2[rownames(BIC2),colnames(BIC2)]
  r <- intersect(rownames(BIC1), rownames(BIC2))
  c <- intersect(colnames(BIC1), colnames(BIC2))
  BIC[r,c] <- pmax(BIC1[r,c], BIC2[r,c], na.rm = TRUE)
  
  attr <- if(pickBIC(BIC2,1) > pickBIC(BIC1,1)) 
          attributes(BIC2) else attributes(BIC1) 
  attr$dim <- dim(BIC)
  attr$dimnames <- dimnames(BIC)
  attr$G <- as.numeric(G)
  attr$modelNames <- modelNames
  attr$returnCodes <- NULL
  attributes(BIC) <- attr
  return(BIC)
}

mclustLoglik <- function(object, ...)
{
  stopifnot(inherits(object, "mclustBIC"))
  BIC <- object
  G <- as.numeric(rownames(BIC))
  modelNames <- colnames(BIC)
  n <- attr(BIC, "n")
  d <- attr(BIC, "d")
  noise <- if(is.null(attr(BIC, "noise"))) FALSE else TRUE
  
  loglik <- matrix(as.double(NA), nrow = length(G), ncol = length(modelNames),
                   dimnames = list(G, modelNames))
  for(i in seq_along(G))
    for(j in seq_along(modelNames))
    {
      npar <- nMclustParams(G = G[i], modelName = modelNames[j], 
                            d = d, noise = noise)
      loglik[i,j] <- 0.5*(BIC[i,j] + npar*log(n))
    }
  
  mostattributes(loglik) <- attributes(BIC)
  attr(loglik, "criterion") <- "loglik"
  class(loglik) <- "mclustLoglik"
  return(loglik)
}

print.mclustLoglik <- function(x, ...)
{
  oldClass(x) <- attr(x, "args") <- NULL
  attr(x, "criterion") <- NULL
  attr(x, "control") <- attr(x, "initialization") <- NULL
  attr(x, "oneD") <- attr(x, "warn") <- attr(x, "Vinv") <- NULL
  attr(x, "prior") <- attr(x, "G") <- attr(x, "modelNames") <- NULL
  attr(x, "returnCodes") <- attr(x, "n") <- attr(x, "d") <- NULL
  
  catwrap("Log-likelihood:")
  NextMethod("print")
  invisible()
}

mclustModel <- function(data, BICvalues, G=NULL, modelNames=NULL, ...)
{
  mc <- match.call(expand.dots = FALSE)
  if (is.null(attr(BICvalues,"initialization")$noise)) {
    mc[[1]] <- as.name("summaryMclustBIC")
  }
  else {
    mc[[1]] <- as.name("summaryMclustBICn")
  }
  nm <- names(mc)
  mc[1:3] <- mc[c(1,3,2)]
  nm[1:3] <- nm[c(1,3,2)]
  nm[nm == "BICvalues"] <- "object" 
  names(mc) <- nm
  ans <- eval(mc, parent.frame())
  ans$classification <- ans$uncertainty <- NULL
  attr( ans, "bestBICvalues") <- NULL
  attr( ans, "prior") <- NULL
  attr( ans, "control") <- NULL
  attr( ans, "initialization") <- NULL
  oldClass(ans) <- "mclustModel"
  ans
}

mclustModelNames <- function(model)
{
  type <- switch(EXPR = as.character(model),
                 "E" = "univariate, equal variance",
                 "V" = "univariate, unequal variance",
                 "EII" = "spherical, equal volume",
                 "VII" = "spherical, varying volume",
                 "EEI" = "diagonal, equal volume and shape",
                 "VEI" = "diagonal, equal shape",
                 "EVI" = "diagonal, equal volume, varying shape",
                 "VVI" = "diagonal, varying volume and shape",
                 "EEE" = "ellipsoidal, equal volume, shape and orientation",
                 "EVE" = "ellipsoidal, equal volume and orientation",
                 "VEE" = "ellipsoidal, equal shape and orientation",
                 "VVE" = "ellipsoidal, equal orientation",
                 "EEV" = "ellipsoidal, equal volume and shape",
                 "VEV" = "ellipsoidal, equal shape",
                 "EVV" = "ellipsoidal, equal volume",
                 "VVV" = "ellipsoidal, varying volume, shape, and orientation",
                 "X"   = "univariate normal",
                 "XII" = "spherical multivariate normal",
                 "XXI" = "diagonal multivariate normal",
                 "XXX" = "ellipsoidal multivariate normal",
                 warning("invalid model"))
  return(list(model = model, type = type))
}

defaultPrior <- function(data, G, modelName, ...)
{
  aux <- list(...)
  if(is.null(aux$shrinkage)) {
    shrinkage <- 0.01
  }
  else if(is.na(aux$shrinkage) || !aux$shrinkage) {
    shrinkage <- 0
  }
  else if(aux$shrinkage < 0) {
    stop("negative value given for shrinkage")
  }
  else {
    shrinkage <- aux$shrinkage
  }
  if(is.null(aux$mean)) {
    mean <- if (is.null(dim(data))) 
      mean(data) else colMeans(data)
  }
  else if(any(is.na(aux$mean))) {
    if(shrinkage)
      stop("positive shrinkage with no prior mean specified")
    mean <- if (is.null(dim(data))) 
      mean(data) else colMeans(data)
  }
  else {
    if(!shrinkage)
      stop("prior mean specified but not shrinkage")
    mean <- aux$mean
  }
  switch(EXPR = modelName,
         E = ,
         V = ,
         X = {
           dof <- 3
           if(is.null(aux$scale)) {
             scale <- var(data)/G^2
           }
           else {
             scale <- aux$scale
           }
           list(shrinkage = shrinkage, mean = mean, dof = dof,
                scale = scale)
         },
         ##
         EII = ,
         VII = ,
         XII = ,
         EEI = ,
         EVI = ,
         VEI = ,
         VVI = ,
         XXI = {
           n <- nrow(data)
           p <- ncol(data)
           dof <- p + 2
           if(is.null(aux$scale)) {
             fac <- (1/G)^(2/p)
             scale <- (fac * sum(apply(data, 2, var)))/
               p
           }
           else {
             scale <- aux$scale
           }
           list(shrinkage = shrinkage, mean = mean, dof = dof,
                scale = scale)
         },
         ##
         EEE = ,
         EVE = ,
         VEE = ,
         VVE = ,
         EEV = ,
         VEV = ,
         EVV = ,
         VVV = ,
         XXX = {
           n <- nrow(data)
           p <- ncol(data)
           dof <- p + 2
           if(is.null(aux$scale)) {
             fac <- (1/G)^(2/p)
             if(n > p) {
               scale <- fac * var(data)
             }
             else {
               scale <- fac * diag(apply(data,
                                         2, var))
             }
           }
           else {
             scale <- aux$scale
           }
           list(shrinkage = shrinkage, mean = mean, dof = dof,
                scale = scale)
         },
         stop("no default prior for this model"))
}

emControl <- function(eps = .Machine$double.eps, 
                      tol = c(1.0e-05, sqrt(.Machine$double.eps)), 
                      itmax = c(.Machine$integer.max, .Machine$integer.max), 
                      equalPro = FALSE)
{
  if(any(eps < 0)) stop("eps is negative")
  if(any(eps >= 1))
    stop("eps is not less than 1")
  if(any(tol < 0))
    stop("tol is negative")
  if(any(tol >= 1))
    stop("tol is not less than 1")
  if(any(itmax < 0))
    stop("itmax is negative")
  if(length(tol) == 1)
    tol <- rep(tol, 2)
  if(length(itmax) == 1)
    itmax <- c(itmax, .Machine$integer.max)
  i <- is.infinite(itmax)
  if(any(i))
    itmax[i] <- .Machine$integer.max
  list(eps = eps, tol = tol, itmax = itmax, equalPro = equalPro)
}

priorControl <- function(functionName = "defaultPrior", ...)
{
  c(list(functionName = functionName), list(...))
}

cdensEEE <- function(data, logarithm = FALSE, parameters, warn = NULL, ...)
{
  dimdat <- dim(data)
  if(is.null(dimdat) || length(dimdat) > 2)
    stop("data must be a matrix or a vector")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  mu <- as.matrix(parameters$mean)
  G <- ncol(mu)
  if(any(is.na(unlist(parameters[c("pro", "mean", "variance")]))) ||
       any(is.null(parameters[c("pro", "mean", "variance")]))) 
  { WARNING <- "parameters are missing"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,G)
    dimnames(z) <- list(dimnames(data)[[1]], NULL)
    return(structure(z, logarithm = logarithm, modelName = "EEE", 
                     WARNING = WARNING, returnCode = 9))
  }
  if(is.null(parameters$variance$cholSigma))
    stop("variance parameters are missing")
  temp <- .Fortran("eseee",
                   as.logical(1),
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$cholSigma),
                   as.double(-1),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(-1),
                   double(p),
                   double(1),
                   double(n * G),
                   PACKAGE = "mclust")[10:12]
  lapackCholInfo <- temp[[1]][1]
  loglik <- temp[[2]]
  z <- matrix(temp[[3]], n, G)
  WARNING <- NULL
  if(lapackCholInfo) {
    if(lapackCholInfo > 0) {
      WARNING <- "sigma is not positive definite"
      if(warn) warning(WARNING)
    }
    else {
      WARNING <- "input error for LAPACK DPOTRF"
      if(warn) warning(WARNING)
    }
    z[] <- NA
    ret <- -9 
  }
  else if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    z[] <- NA
    ret <- -1
  }
  else {
    if (!logarithm) z <- exp(z)
    ret <- 0
  }
  dimnames(z) <- list(dimnames(data)[[1]],NULL)
  structure(z, logarithm = logarithm, modelName = "EEE",
            WARNING = WARNING, returnCode = ret)
}

emEEE <- function(data, parameters, prior = NULL, control = emControl(), 
                  warn = NULL, ...)
{
  z <- estepEEE(data, parameters = parameters, warn = warn)$z  
  meEEE(data, z = z, prior = prior, control = control, 
        Vinv = parameters$Vinv, warn = warn)
}

estepEEE <- function(data, parameters, warn = NULL, ...)
{
  if (is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  if(is.null(dimdat) || length(dimdat) > 2)
    stop("data must be a matrix or a vector")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  pro <- parameters$pro
  pro <- pro/sum(pro)
  l <- length(pro)
  mu <- as.matrix(parameters$mean)
  G <- ncol(mu)
  noise <- (l == G + 1)
  if(!noise) {
    if(l != G)
      stop("pro improperly specified")
    K <- G
    Vinv <- NULL
  }
  else {
    K <- G + 1
    Vinv <- parameters$Vinv
    if(is.null(Vinv) || Vinv <= 0)
      Vinv <- hypvol(data, reciprocal = TRUE)
  }
  if(any(is.na(unlist(parameters[c("pro", "mean", "variance")]))) ||
       any(is.null(parameters[c("pro", "mean", "variance")]))) {
    WARNING <- "parameters are missing"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,K)
    dimnames(z) <- list(dimnames(data)[[1]], NULL)
    return(structure(list(modelName = "EEE", n=n, d=p, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  if (is.null(parameters$variance$cholSigma))
    stop("variance parameters are missing")
  temp <- .Fortran("eseee",
                   as.logical(1),
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$cholSigma),
                   as.double(pro),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(if (is.null(Vinv)) -1 else Vinv),
                   double(p),
                   double(1),
                   double(n * K),
                   PACKAGE = "mclust")[10:12]
  lapackCholInfo <- temp[[1]][1]
  loglik <- temp[[2]]
  z <- matrix(temp[[3]], n, K)
  WARNING <- NULL
  if(lapackCholInfo) {
    if(lapackCholInfo > 0) {
      WARNING <- "sigma is not positive definite"
      warning(WARNING)
      ret <- -4 
    }
    else {
      WARNING <- "input error for LAPACK DPOTRF"
      warning(WARNING)
      ret <- -5
    }
    z[] <- loglik <- NA
  }
  else if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    z[] <- loglik <- NA
    ret <- -1
  }
  else ret <- 0
  dimnames(z) <- list(dimnames(data)[[1]],NULL)
  structure(list(modelName = "EEE", n = n, d = p, G = G, 
                 z = z, parameters = parameters, loglik = loglik),
            WARNING = WARNING, returnCode = ret)
}

meEEE <- function(data, z, prior = NULL, control = emControl(), 
                  Vinv = NULL, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)
  if(oneD || length(dimdat) != 2)
    stop("data should in the form of a matrix")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  z <- as.matrix(z)
  dimz <- dim(z)
  if(dimz[1] != n)
    stop("data and z should have the same row dimension")
  K <- dimz[2]
  if (!is.null(Vinv)) {
    G <- K - 1
    if(Vinv <= 0) Vinv <- hypvol(data, reciprocal = TRUE)
  }
  else G <- K
  if(all(is.na(z))) {
    WARNING <- "z is missing"
    if(warn) warning(WARNING)
    variance <- list(modelName = "EEE", d = p, G = G, 
                     Sigma = matrix(as.double(NA), p, p), cholSigma = matrix(as.double(NA), p, p)) 
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance, Vinv=Vinv)
    return(structure(list(modelName="EEE", prior=prior, n=n, d=p, 
                          G=G, z=z, parameters=parameters, 
                          control=control, loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  if(any(is.na(z)) || any(z < 0) || any(z > 1))
    stop("improper specification of z")
  storage.mode(z) <- "double"
  if(is.null(prior)) {
    temp <- .Fortran("meeee",
                     as.logical(control$equalPro),
                     as.double(data),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(if (is.null(Vinv)) -1 else Vinv),
                     z,
                     as.integer(control$itmax[1]),
                     as.double(control$tol[1]),
                     as.double(control$eps),
                     double(p * G),
                     double(p * p),
                     double(K),
                     double(p),
                     PACKAGE = "mclust")[7:13]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEE"), 
                                                 prior[names(prior) != "functionName"]))
    temp <- .Fortran("meeeep",
                     as.logical(control$equalPro),
                     as.double(data),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(if (is.null(Vinv)) -1 else Vinv),
                     as.double(priorParams$shrinkage),
                     as.double(priorParams$mean),
                     as.double(if(any(priorParams$scale != 0)) chol(priorParams$
                                                                      scale) else priorParams$scale),
                     as.double(priorParams$dof),
                     z,
                     as.integer(control$itmax[1]),
                     as.double(control$tol[1]),
                     as.double(control$eps),
                     double(p * G),
                     double(p * p),
                     double(K),
                     double(p),
                     PACKAGE = "mclust")[c(11:17, 10)]
  }
  z <- temp[[1]]
  its <- temp[[2]]
  err <- temp[[3]]
  loglik <- temp[[4]]
  mu <- matrix(temp[[5]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  cholSigma <- matrix(temp[[6]], p, p)
  pro <- temp[[7]]
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    mu[] <- pro[] <- z[] <- loglik <- NA
    sigma <- array(NA, c(p, p, G))
    Sigma <- matrix( NA, p, p)
    ret <- -1
  }
  else if(loglik <  - signif(.Machine$double.xmax, 6)) {
    if(control$equalPro) {
      WARNING <- "z column sum fell below threshold"
      if(warn) warning(WARNING)
    }
    else {
      WARNING <- "mixing proportion fell below threshold"
      if(warn) warning(WARNING)
    }
    mu[] <- pro[] <- z[] <- loglik <- logprior <- NA
    sigma <- array(NA, c(p, p, G))
    Sigma <- matrix(as.double(NA), p, p)
    ret <- if(control$equalPro) -2 else -3
  }
  else {
    Sigma <- unchol(cholSigma, upper = TRUE)
    sigma <- array(0, c(p, p, G))
    for(k in 1:G)
      sigma[,  , k] <- Sigma
    if(its >= control$itmax[1]) {
      WARNING <- "iteration limit reached"
      if(warn) warning(WARNING)
      its <-  - its
      ret <- 1
    }
    else ret <- 0
  }
  info <- c(iterations = its, error = err)
  dimnames(z) <- list(dimnames(data)[[1]], NULL)
  dimnames(mu) <- list(dimnames(data)[[2]], NULL)
  dimnames(Sigma) <- dimnames(cholSigma) <- 
    list(dimnames(data)[[2]], dimnames(data)[[2]])
  dimnames(sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]],
                          NULL)
  variance <- list(modelName = "EEE", d = p, G = G,
                   sigma = sigma, Sigma = Sigma, cholSigma = cholSigma) 
  parameters <- list(pro=pro, mean=mu, variance=variance, Vinv=Vinv)
  structure(list(modelName = "EEE", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters, control = control,
                 loglik = loglik), 
            info = info, WARNING = WARNING, returnCode = ret)
}

mstepEEE <- function(data, z, prior = NULL,  warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)
  if(oneD || length(dimdat) != 2)
    stop("data should be a matrix or a vector")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  ##
  z <- as.matrix(z)
  dimz <- dim(z)
  if(dimz[1] != n)
    stop("row dimension of z should equal data length")
  G <- dimz[2]
  if(all(is.na(z))) {
    WARNING <- "z is missing"
    if(warn) warning(WARNING)
    variance <- list(modelName = "EEE", d = p, G = G, 
                     sigma <- array(NA, c(p,p, G)), 
                     Sigma = matrix(as.double(NA), p, p), cholSigma = matrix(as.double(NA), p, p)) 
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance)
    return(structure(list(modelName="EEE", prior=prior, n=n, d=p, 
                          G=G, z=z, parameters=parameters, 
                          loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  if(any(is.na(z)) || any(z < 0) || any(z > 1))
    stop("improper specification of z")
  if(is.null(prior)) {
    temp <- .Fortran("mseee",
                     as.double(data),
                     as.double(z),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     double(p),
                     double(p * G),
                     double(p * p),
                     double(G),
                     PACKAGE = "mclust")[7:9]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEE"), 
                                                 prior[names(prior) != "functionName"]))
    temp <- .Fortran("mseeep",
                     as.double(data),
                     as.double(z),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(priorParams$shrinkage),
                     as.double(priorParams$mean),
                     as.double(if(any(priorParams$scale != 0)) chol(priorParams$scale) else priorParams$scale),
                     as.double(priorParams$dof),
                     double(p),
                     double(p * G),
                     double(p * p),
                     double(G),
                     PACKAGE = "mclust")[11:13]
  }
  mu <- matrix(temp[[1]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  cholSigma <- matrix(temp[[2]], p, p)
  pro <- temp[[3]]
  sigma <- array(0, c(p, p, G))
  Sigma <- unchol(cholSigma, upper = TRUE)
  for(k in 1:G)
    sigma[,  , k] <- Sigma
  WARNING <- NULL
  if(any(mu > signif(.Machine$double.xmax, 6))) {
    WARNING <- "cannot compute M-step"
    if(warn) warning(WARNING)
    mu[] <- sigma[] <- Sigma[] <- cholSigma[] <- NA
    ret <- -1
  }
  else ret <- 0
  dimnames(z) <- list(dimnames(data)[[1]], NULL)
  dimnames(mu) <- list(dimnames(data)[[2]], NULL)
  dimnames(Sigma) <- dimnames(cholSigma) <- 
    list(dimnames(data)[[2]], dimnames(data)[[2]])
  dimnames(sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]],
                          NULL)
  variance <- list(modelName = "EEE", d = p, G = G, 
                   sigma = sigma, Sigma = Sigma, cholSigma= cholSigma)
  parameters <- list(pro=pro, mean=mu, variance=variance)
  structure(list(modelName = "EEE", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters), 
            WARNING = WARNING, returnCode = ret)
}

simEEE <- function(parameters, n, seed = NULL, ...)
{
  if(!is.null(seed)) set.seed(seed)
  mu <- as.matrix(parameters$mean)
  d <- nrow(mu)
  G <- ncol(mu)
  if(any(is.na(parameters[c("mean", "variance")])) || 
     any(is.null(parameters[c("mean", "variance")]))) 
    { warning("parameters are missing")
      return(structure(matrix(as.double(NA), n, d + 1), modelName = "EEE"))
  }
  pro <- parameters$pro
  if(is.null(pro))
    pro <- rep(1/G, G)
  clabels <- sample(1:G, size = n, replace = TRUE, prob = pro)
  ctabel <- tabulate(clabels, nbins=G)
  x <- matrix(0, n, d)
  if(is.null(cholSigma <- parameters$variance$cholSigma)) {
    if(is.null(Sigma <- parameters$variance$Sigma)) {
      stop("variance parameters must inlcude either Sigma or cholSigma"
      )
    }
    cholSigma <- chol(Sigma)
  }
  for(k in 1:G) {
    m <- ctabel[k]
    x[clabels == k,] <- sweep(matrix(rnorm(m * d), nrow = m, ncol = d) %*% 
                              cholSigma, MARGIN = 2, STATS = mu[,k], FUN = "+")
  }
  dimnames(x) <- list(NULL, paste0("x", 1:d))
  structure(cbind(group = clabels, x), modelName = "EEE")
}

cdensEEI <- function(data, logarithm = FALSE, parameters, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  if(is.null(dimdat) || length(dimdat) != 2)
    stop("data must be a matrix")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  mu <- as.matrix(parameters$mean)
  G <- ncol(mu)
  if(any(is.na(unlist(parameters[c("pro", "mean", "variance")]))) ||
       any(is.null(parameters[c("pro", "mean", "variance")]))) {
    WARNING <- "parameters are missing"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,G)
    dimnames(z) <- list(dimnames(data)[[1]], NULL)
    return(structure(z, logarithm = logarithm, modelName = "EEI", 
                     WARNING = WARNING, returnCode = 9))
  }
  if (is.null(parameters$variance$scale) ||
        is.null(parameters$variance$shape)) 
    stop("variance parameters are missing")
  temp <- .Fortran("eseei",
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$scale),
                   as.double(parameters$variance$shape),
                   as.double(-1),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(-1),
                   double(1),
                   double(n * G),
                   PACKAGE = "mclust")[10:11]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, G)
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "sigma-squared falls below threshold"
    if(warn) warning(WARNING)
    z[] <- NA
    ret <- -1 
  }
  else {
    if (!logarithm) z <- exp(z)
    ret <- 0
  }
  dimnames(z) <- list(dimnames(data)[[1]],NULL)
  structure(z, logarithm = logarithm, modelName = "EEI",
            WARNING = WARNING, returnCode = ret)
}

cdensEII <- function(data, logarithm = FALSE, parameters, warn = NULL, ...)
  {
    if(is.null(warn)) warn <- mclust.options("warn")
    dimdat <- dim(data)
    if(is.null(dimdat) || length(dimdat) != 2)
      stop("data must be a matrix")
    data <- as.matrix(data)
    n <- nrow(data)
    p <- ncol(data)
    mu <- as.matrix(parameters$mean)
    G <- ncol(mu)
    if(any(is.na(unlist(parameters[c("pro", "mean", "variance")]))) ||
         any(is.null(parameters[c("pro", "mean", "variance")]))) {
      WARNING <- "parameters are missing"
      if(warn) warning(WARNING)
      z <- matrix(as.double(NA),n,G)
      dimnames(z) <- list(dimnames(data)[[1]], NULL)
      return(structure(z, logarithm = logarithm, modelName = "EII", 
                       WARNING = WARNING, returnCode = 9))
    }
    sigmasq <- parameters$variance$sigmasq
    if(sigmasq < 0)
      stop("sigma-squared is negative")
    if(!sigmasq) {
      WARNING <- "sigma-squared vanishes"
      if(warn) warning(WARNING)
      z <- matrix(as.double(NA),n,G)
      dimnames(z) <- list(dimnames(data)[[1]], NULL)
      return(structure(z, logarithm = logarithm, modelName = "EII", 
                       WARNING = WARNING, returnCode = 9))
    }
    temp <- .Fortran("eseii",
                     as.double(data),
                     as.double(mu),
                     as.double(sigmasq),
                     as.double(-1),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(-1),
                     double(1),
                     double(n * G),
                     PACKAGE = "mclust")[9:10]
    loglik <- temp[[1]]
    z <- matrix(temp[[2]], n, G)
    WARNING <- NULL
    if(loglik > signif(.Machine$double.xmax, 6)) {
      WARNING <- "sigma-squared falls below threshold"
      if(warn) warning(WARNING)
      z[] <- NA
      ret <- -1
    }
    else {
      if (!logarithm) z <- exp(z)
      ret <- 0
    }
    dimnames(z) <- list(dimnames(data)[[1]],NULL)
    structure(z, logarithm = logarithm, modelName = "EII",
              WARNING = WARNING, returnCode = ret)
  }

emEEI <- function(data, parameters, prior = NULL, control = emControl(), 
                  warn = NULL, ...)
{
  z <- estepEEI(data, parameters = parameters, warn = warn)$z  
  meEEI(data, z = z, prior = prior, control = control, 
        Vinv = parameters$Vinv, warn = warn)
}

estepEEI <- function(data, parameters, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  if(is.null(dimdat) || length(dimdat) != 2)
    stop("data must be a matrix")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  pro <- parameters$pro
  pro <- pro/sum(pro)
  l <- length(pro)
  mu <- as.matrix(parameters$mean)
  G <- ncol(mu)
  noise <- (l == G + 1)
  if(!noise) {
    if(l != G)
      stop("pro improperly specified")
    K <- G
    Vinv <- NULL
  }
  else {
    K <- G + 1
    Vinv <- parameters$Vinv
    if(is.null(Vinv) || Vinv <= 0)
      Vinv <- hypvol(data, reciprocal = TRUE)
  } 
  if(any(is.na(unlist(parameters[c("pro", "mean", "variance")]))) ||
       any(is.null(parameters[c("pro", "mean", "variance")]))) {
    WARNING <- "parameters are missing"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,K)
    dimnames(z) <- list(dimnames(data)[[1]], NULL)
    return(structure(list(modelName = "EEI", n=n, d=p, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  if (is.null(parameters$variance$scale) ||
        is.null(parameters$variance$shape)) 
    stop("variance parameters are missing")
  temp <- .Fortran("eseei",
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$scale),
                   as.double(parameters$variance$shape),
                   as.double(pro),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(if (is.null(Vinv)) -1 else Vinv),
                   double(1),
                   double(n * K),
                   PACKAGE = "mclust")[10:11]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, K)
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    z[] <- loglik <- NA
    ret <- -1
  }
  else ret <- 0
  dimnames(z) <- list(dimnames(data)[[1]],NULL)
  structure(list(modelName = "EEI", n = n, d = p, G = G, 
                 z = z, parameters = parameters, loglik = loglik),
            WARNING = WARNING, returnCode = ret)
}

meEEI <- function(data, z, prior = NULL, control = emControl(), 
                  Vinv = NULL, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)
  if(oneD || length(dimdat) > 2)
    stop("data  should be in the form of a matrix")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  z <- as.matrix(z)
  dimz <- dim(z)
  if(dimz[1] != n)
    stop("data and z should have the same row dimension")
  K <- dimz[2]
  if (!is.null(Vinv)) {
    G <- K - 1
    if(Vinv <= 0) Vinv <- hypvol(data, reciprocal = TRUE)
  }
  else G <- K
  if(all(is.na(z))) {
    WARNING <- "z is missing"
    if(warn) warning(WARNING)
    variance <- list(modelName = "EEI", d = p, G = G, 
                     scale = NA, shape = rep(NA,p)) 
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance, Vinv=Vinv)
    return(structure(list(modelName="EEI", prior=prior, n=n, d=p, 
                          G=G, z=z, parameters=parameters, 
                          control=control, loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  if(any(is.na(z)) || any(z < 0) || any(z > 1))
    stop("improper specification of z")
  storage.mode(z) <- "double"
  if(is.null(prior)) {
    temp <- .Fortran("meeei",
                     as.logical(control$equalPro),
                     as.double(data),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(if (is.null(Vinv)) -1 else Vinv),
                     z,
                     as.integer(control$itmax[1]),
                     as.double(control$tol[1]),
                     as.double(control$eps),
                     double(p * G),
                     double(1),
                     double(p),
                     double(K),
                     PACKAGE = "mclust")[7:14]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEI"), 
                                                 prior[names(prior) != "functionName"]))
    temp <- .Fortran("meeeip",
                     as.logical(control$equalPro),
                     as.double(data),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(if (is.null(Vinv)) -1 else Vinv),
                     as.double(priorParams$shrinkage),
                     as.double(priorParams$mean),
                     as.double(priorParams$scale),
                     as.double(priorParams$dof),
                     z,
                     as.integer(control$itmax[1]),
                     as.double(control$tol[1]),
                     as.double(control$eps),
                     double(p * G),
                     double(1),
                     double(p),
                     double(K),
                     PACKAGE = "mclust")[11:18]
  }
  z <- temp[[1]]
  its <- temp[[2]]
  err <- temp[[3]]
  loglik <- temp[[4]]
  mu <- matrix(temp[[5]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  scale <- temp[[6]]
  shape <- temp[[7]]
  pro <- temp[[8]]
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    sigma <- array(NA, c(p, p, G))
    Sigma <- matrix(as.double(NA), p, p)
    mu[] <- pro[] <- z[] <- loglik <- shape[] <- NA
    ret <- -1
  }
  else if(loglik <  - signif(.Machine$double.xmax, 6)) {
    if(control$equalPro) {
      WARNING <- "z column sum fell below threshold"
      if(warn) warning(WARNING)
    }
    else {
      WARNING <- "mixing proportion fell below threshold"
      if(warn) warning(WARNING)
    }
    sigma <- array(NA, c(p, p, G))
    Sigma <- matrix(as.double(NA), p, p)
    mu[] <- pro[] <- z[] <- loglik <- shape[] <- NA
    ret <- if(control$equalPro) -2 else -3
  }
  else {
    sigma <- array(0, c(p, p, G))
    Sigma <- diag(scale * shape)
    for(k in 1:G)
      sigma[,  , k] <- Sigma
    if(its >= control$itmax[1]) {
      WARNING <- "iteration limit reached"
      if(warn) warning(WARNING)
      its <-  - its
      ret <- 1
    }
    else ret <- 0
  }
  info <- c(iterations = its, error = err)
  dimnames(z) <- list(dimnames(data)[[1]], NULL)
  dimnames(mu) <- list(dimnames(data)[[2]], NULL)
  dimnames(Sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]])
  dimnames(sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]],
                          NULL)
  variance <- list(modelName = "EEI", d = p, G = G, 
                   sigma = sigma, Sigma = Sigma, 
                   scale = scale, shape = shape)
  parameters <- list(pro=pro, mean=mu, variance=variance, Vinv=Vinv)
  structure(list(modelName = "EEI", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters, control = control,
                 loglik = loglik), 
            info = info, WARNING = WARNING, returnCode = ret)
}

mstepEEI <- function(data, z, prior = NULL, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)
  if(oneD || length(dimdat) != 2)
    stop("data should be a matrix or a vector")
  data <- as.matrix(data)
  n <- nrow(data)
  p <- ncol(data)
  z <- as.matrix(z)
  dimz <- dim(z)
  if(dimz[1] != n)
    stop("row dimension of z should equal data length")
  G <- dimz[2]
  if(all(is.na(z))) {
    WARNING <- "z is missing"
    if(warn) warning(WARNING)
    variance <- list(modelName = "EEI", d = p, G = G, 
                     scale = NA, shape = rep(NA,p)) 
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance)
    return(structure(list(modelName="EEI", prior=prior, n=n, d=p, 
                          G=G, z=z, parameters=parameters), 
                     WARNING = WARNING, returnCode = 9))
    
  }
  if(any(is.na(z)) || any(z < 0) || any(z > 1))
    stop("improper specification of z")
  if(is.null(prior)) {
    temp <- .Fortran("mseei",
                     as.double(data),
                     as.double(z),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     double(p * G),
                     double(1),
                     double(p),
                     double(G),
                     PACKAGE = "mclust")[6:9]
  }
  else {
    storage.mode(z) <- "double"
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEI"), prior[names(
                                                          prior) != "functionName"]))
    temp <- .Fortran("mseeip",
                     as.double(data),
                     as.double(z),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     as.double(priorParams$shrinkage),
                     as.double(priorParams$mean),
                     as.double(priorParams$scale),
                     as.double(priorParams$dof),
                     double(p * G),
                     double(1),
                     double(p),
                     double(G),
                     PACKAGE = "mclust")[10:13]
  }
  mu <- matrix(temp[[1]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  scale <- temp[[2]]
  shape <- temp[[3]]
  pro <- temp[[4]]
  WARNING <- NULL
  if(any(c(shape, scale) > signif(.Machine$double.xmax, 6)) || any(!c(
    scale, shape))) {
    WARNING <- "cannot compute M-step"
    if(warn) warning(WARNING)
    mu[] <- pro[] <- scale <- shape[] <- NA
    sigma <- Sigma <- array(NA, c(p, p, G))
    ret <- -1
  }
  else {
    sigma <- array(0, c(p, p, G))
    Sigma <- diag(scale * shape)
    for(k in 1:G)
      sigma[,  , k] <- Sigma
    ret <- 0
  }
  dimnames(z) <- list(dimnames(data)[[1]], NULL)
  dimnames(mu) <- list(dimnames(data)[[2]], NULL)
  dimnames(Sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]])
  dimnames(sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]],
                          NULL)
  variance <- list(modelName = "EEI", d = p, G = G, 
                   sigma = sigma, Sigma = Sigma, 
                   scale = scale, shape = shape)
  parameters <- list(pro=pro, mean=mu, variance=variance)
  structure(list(modelName = "EEI", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters), 
            WARNING = WARNING, returnCode = ret)
}

simEEI <- function(parameters, n, seed = NULL, ...)
{
  if(!is.null(seed)) set.seed(seed)
  mu <- as.matrix(parameters$mean)
  d <- nrow(mu)
  G <- ncol(mu)
  if(any(is.na(parameters[c("mean", "variance")])) || 
     any(is.null(parameters[c("mean", "variance")]))) 
    { warning("parameters are missing")
      return(structure(matrix(as.double(NA), n, d + 1), modelName = "EEI"))
  }
  pro <- parameters$pro
  if(is.null(pro))
    pro <- rep(1/G, G)
  clabels <- sample(1:G, size = n, replace = TRUE, prob = pro)
  ctabel <- tabulate(clabels, nbins=G)
  x <- matrix(0, n, d)
  shape <- parameters$variance$shape
  if(length(shape) != d)
    stop("shape incompatible with mean")
  cholSigma <- diag(sqrt(parameters$variance$scale * shape))
  for(k in 1:G) {
    m <- ctabel[k]
    x[clabels == k,  ] <- sweep(matrix(rnorm(m * d), nrow = m,
                                       ncol = d) %*% cholSigma, MARGIN = 2, STATS = mu[, k],
                                FUN = "+")
  }
  dimnames(x) <- list(NULL, paste0("x", 1:d))
  structure(cbind(group = clabels, x), modelName = "EEI")
}

cdensE <- function(data, logarithm = FALSE, parameters, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)
  if(!oneD)
    stop("data must be one-dimensional")
  data <- drop(data)
  n <- length(data)
  mu <- drop(parameters$mean)
  G <- length(mu)
  if(any(is.na(unlist(parameters[c("mean", "variance")]))) ||
       any(is.null(parameters[c("mean", "variance")]))) {
    WARNING <- "parameters are missing"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,G)
    dimnames(z) <- list(names(data), NULL)
    return(structure(z, logarithm = logarithm, modelName = "E",
                     WARNING = WARNING, returnCode = 9))
  }
  sigmasq <- parameters$variance$sigmasq
  if(is.null(sigmasq))
    stop("variance parameters are missing")
  if(length(sigmasq) > 1)
    if(warn) warning("more than one sigma-squared given")
  if(sigmasq < 0)
    stop("sigma-squared is negative")
  if(!sigmasq) {
    WARNING <- "sigma-squared vanishes"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,G)
    dimnames(z) <- list(names(data), NULL)
    return(structure(z, logarithm = logarithm, modelName = "E",
                     WARNING = WARNING, returnCode = 9))
  }
  temp <- .Fortran("es1e",
                   as.double(data),
                   as.double(mu),
                   as.double(sigmasq),
                   as.double(-1),
                   as.integer(n),
                   as.integer(G),
                   as.double(-1),
                   double(1),
                   double(n * G),
                   PACKAGE = "mclust")[8:9]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, G)
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "sigma-squared falls below threshold"
    if(warn) warning(WARNING)
    z[] <- NA
    ret <- -1
  }
  else {
    if (!logarithm) z <- exp(z)
    ret <- 0
  } 
  dimnames(z) <- list(names(data),NULL)
  structure(z, logarithm = logarithm, modelName = "E", 
            WARNING = WARNING, returnCode = ret) 
}

emE <- function(data, parameters, prior = NULL, control = emControl(), 
                warn = NULL, ...)
{
  z <- estepE(data, parameters = parameters, warn = warn)$z
  meE(data, z = z, prior = prior, control = control, 
      Vinv = parameters$Vinv, warn = warn)
}

estepE <- function(data, parameters, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- (is.null(dimdat) || NCOL(data) == 1)