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 hypvol grid2 grid1 decomp2sigma sigma2decomp sigma2decomp plot.mclustBIC plot.Mclust mvn2plot mclust2Dplot mclust1Dplot 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 pickBIC print.summary.mclustBIC summaryMclustBICn summaryMclustBIC summary.mclustBIC print.mclustBIC mclustBIC predict.Mclust logLik.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 grid1 grid2 hypvol logLik.Mclust Mclust mclust1Dplot mclust2Dplot mclustBIC 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 mvn2plot mvnX mvnXII mvnXXI mvnXXX nMclustParams nVarParams pickBIC plot.Mclust plot.mclustBIC predict.Mclust print.Mclust print.mclustBIC 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 summaryMclustBIC 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)
  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())
  G <- attr(Bic, "G")
  modelNames <- attr(Bic, "modelNames")
  Sumry <- summary(Bic, data, G = G, modelNames = modelNames)
  if(length(Sumry)==0) 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$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", "bic", "loglik", "df", 
                    "hypvol", "parameters", "z", 
                    "classification", "uncertainty")
  structure(ans[orderedNames], class = "Mclust")  
}

print.Mclust <- function(x, digits = getOption("digits"), ...)
{
  cat("\'", class(x)[1], "\' model object:\n", sep = "")
  G <- x$G
  noise <- !is.null(attr(x$BIC, "Vinv"))
  if(G == 0 & noise)
    { cat(" best model: single noise component\n") }
  else
    { M <- mclustModelNames(x$modelName)$type
     cat(" best model: ", M, " (", x$model, ") with ", G, " components\n", 
         if(noise) " and a noise term\n", sep = "") }
  invisible()
}

summary.Mclust <- function(object, parameters = FALSE, classification = FALSE, ...)
{
  # 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")
  else
    title <- paste("Density estimation via Gaussian finite mixture modeling")
  #
  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 = icl(object),
              pro = pro, mean = mean, variance = sigma,
              noise = noise,
              prior = attr(object$BIC, "prior"), 
              classification = object$classification, 
              printParameters = parameters, 
              printClassification = classification)
  class(obj) <- "summary.Mclust"
  return(obj)
}

print.summary.Mclust <- function(x, digits = getOption("digits"), ...)
{
  cat(rep("-", nchar(x$title)),"\n",sep="")
  cat(x$title, "\n")
  cat(rep("-", nchar(x$title)),"\n",sep="")
  #
  if(x$G == 0)
    { cat("\nMclust model with only a noise component:\n\n") }
  else
    { cat("\nMclust ", x$modelName, " (", 
          mclustModelNames(x$modelName)$type, ") model with ", 
          x$G, ifelse(x$G > 1, " components", " component"),
          if(x$noise) "\nand a noise term", ":\n\n",
          sep = "") 
    }
  #
  if(!is.null(x$prior))
    { cat("Prior: ")
      cat(x$prior$functionName, "(", 
          paste(names(x$prior[-1]), x$prior[-1], sep = " = ", 
                collapse = ", "), ")", sep = "")
      cat("\n\n")
  }
  #
  tab <- data.frame("log-likelihood" = x$loglik, "n" = x$n, 
                    "df" = x$df, "BIC" = x$bic, "ICL" = x$icl, 
                    row.names = "")
  print(tab, digits = digits)
  #
  cat("\nClustering table:")
  print(table(factor(x$classification, 
                     levels = { l <- seq_len(x$G)
                                if(is.numeric(x$noise)) l <- c(l,0) 
                                l })),
        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") }
  }
  if(x$printClassification)
  { cat("\nClassification:\n")
    print(x$classification, digits = digits)
  }
  #
  invisible(x)
}

# old version, wrong df, data not needed, inefficient because compute dens on original scale
# logLik.Mclust <- function(object, data, ...)
# {
#   if(!missing(data))
#     object$data <- data.matrix(data)
#   par <- object$parameters
#   df <- with(object, (G-1) + G*d + nVarParams(modelName, d = d, G = G))
# 
#   #  l <- matrix(as.double(NA), object$n, object$G)
#   #  for(k in seq(object$G))
#   #     { l[,k] <- par$pro[k] * dmvnorm(data, par$mean[,k],
#   #                                           par$variance$sigma[,,k]) }
#   #  l <- sum(log(rowSums(l)))
#   
#   l <- sum(log(do.call("dens", object)))
#   attr(l, "nobs") <- object$n
#   attr(l, "df") <- df
#   class(l) <- "logLik"
#   return(l)
# }

logLik.Mclust <- function(object, ...)
{
  par <- object$parameters
  # df <- with(object, (G-1) + G*d + nVarParams(modelName, d = d, G = G))
  df <- with(object, nMclustParams(modelName, d, G, 
                                   noise = (!is.na(hypvol)), 
                                   equalPro = attr(BIC, "control")$equalPro))
  l <- sum(do.call("dens", c(object, logarithm = TRUE)))
  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
  prior <- object$parameters$pro
  noise <- (!is.na(object$hypvol))
  # old
  # z <- do.call("cdens", object)
  # z <- sweep(z, MARGIN = 1, FUN = "/", STATS = apply(z, 1, max))
  # z <- sweep(z, MARGIN = 2, FUN = "*", STATS = prior/sum(prior))
  # z <- sweep(z, MARGIN = 1, STATS = apply(z, 1, sum), FUN = "/")
  # new: more efficient and accurate
  z <- do.call("cdens", c(object, list(logarithm = TRUE)))
  if(noise)
    z <- cbind(z, log(object$parameters$Vinv))
  z <- sweep(z, MARGIN = 2, FUN = "+", STATS = log(prior/sum(prior)))
  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
          # TODO: remove after check
          # if(!is.logical(noise)) 
          #   { if(any(match(noise, 1:n, nomatch = 0) == 0))
          #       stop("numeric noise must correspond to row indexes of data")
          #     noise <- as.logical(match(1:n, noise, nomatch = 0))
          # }
          # initialization$noise <- noise
          # nnoise <- sum(as.numeric(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 (if not provided) for initialization when data size is
  # larger than the value specified in mclust.options()
  if(n > .mclust$subset & is.null(initialization$subset))
    { initialization$subset <- sample(seq.int(n), 
                                      size = .mclust$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("hcModelNames")[1])
          }
          else {
            hcPairs <- hc(data = data, modelName = "EII")
          } 
        }
        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(modelName = modelName, data = data, 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
      # TODO: remove after check 
      # if (is.logical(subset)) subset <- which(subset)
      if (is.null(initialization$hcPairs)) {
        if (d != 1) {
          if (n > d) {
            hcPairs <- hc(modelName = mclust.options("hcModelNames")[1], 
                          data = data[subset,])
          }
          else {
            hcPairs <- hc(modelName = "EII", 
                          data = data[subset,])
          }
        }
        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(modelName = modelName, z = z, 
                      data = as.matrix(data)[initialization$subset,],
                      prior = prior, control = control, warn = warn)
          #
          #  ctrl <- control
          #  ctrl$itmax[1] <- 1
          #  ms <- me(modelName = modelName, data = as.matrix(data)[
          #           initialization$subset,  ], z = z, prior = prior, control = ctrl)
          #
          es <- do.call("estep", c(list(data = data, warn = warn), ms))
          out <- me(modelName = modelName, data = data, 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(modelName = mclust.options("hcModelNames")[1], 
                          data = data[-noise,  ])
          }
          else {
            hcPairs <- hc(modelName = "EII", data = data[-noise,  ])
          }
        }
        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(modelName = modelName, data = data, 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(modelName = mclust.options("hcModelNames")[1], 
                          data = data[subset,])
          }
          else {
            hcPairs <- hc(modelName = "EII", 
                          data = data[subset,])
          }
        }
        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(modelName = modelName, z = z, 
                      data = as.matrix(data)[subset,],
                      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(modelName = modelName, data = data, 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
  
  cat("Bayesian Information Criterion (BIC):\n")
  NextMethod("print")
  cat("\n")
  cat("Top", pick, "models based on the BIC criterion:\n")
  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(modelName = bestModel, data = data, 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(modelName = bestModel, data = data,
                    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(modelName = bestModel, prior = prior, z = z, 
                data = as.matrix(data)[subset,], control = control, 
                warn = warn)
    es <- do.call("estep", c(list(data = data), ms))
    out <- me(modelName = bestModel, data = data, z = es$z, 
              prior = prior, control = control, warn = warn)
    # perform extra M-step and update parameters
    ms <- mstep(modelName = bestModel, data = data, 
                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
  # TODO: remove after check
  # if(!is.logical(noise)) 
  #   noise <- as.logical(match(1:n, noise, nomatch = 0))
  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(modelName = bestModel, data = data, 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(modelName = bestModel, data = as.matrix(data)[subset,], 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(modelName = bestModel, data = data, 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("\nClassification table for model (", colnames(bic)[1], "):", sep = "")
  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()
}

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)) 
  { 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
}

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, retrunCode = 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) || length(dimdat[dimdat > 1]) == 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) || length(dimdat[dimdat > 1]) == 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) || length(dimdat[dimdat > 1]) == 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) || length(dimdat[dimdat > 1]) == 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) || length(dimdat[dimdat > 1]) == 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) || length(dimdat[dimdat > 1]) == 1
  if(!oneD)
    stop("data must be one-dimensional")
  data <- drop(data)
  n <- length(data)
  pro <- parameters$pro
  pro <- pro/sum(pro)
  l <- length(pro)
  
  mu <- drop(parameters$mean)
  G <- length(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(names(data), NULL)
    return(structure(list(modelName = "E", n=n, d=1, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     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 specified")
  if(sigmasq < 0)
    stop("sigma-squared is negative")
  if(!sigmasq) {
    WARNING <- "sigma-squared vanishes"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,K)
    dimnames(z) <- list(names(data), NULL)
    return(structure(list(modelName = "E", n=n, d=1, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     WARNING = WARNING, returnCode = -1))
  }
  temp <- .Fortran("es1e",
                   as.double(data),
                   as.double(mu),
                   as.double(sigmasq),
                   as.double(pro),
                   as.integer(n),
                   as.integer(G),
                   as.double(if (is.null(Vinv)) -1 else Vinv),
                   double(1),
                   double(n * K),
                   PACKAGE = "mclust")[8:9]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, K)
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "cannot compute E-step"
    if(warn) warning(WARNING)
    z[] <- loglik <- NA
    ret <- -1
  }
  else ret <- 0
  dimnames(z) <- list(names(data),NULL) 
  structure(list(modelName = "E", n = n, d = 1, G = G, 
                 z = z, parameters = parameters, loglik = loglik),
            WARNING = WARNING, returnCode = ret)
}

cdensEEV <- function(data, logarithm = FALSE, parameters, warn = NULL, ...)
{
  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 = "EEV", 
                     WARNING = WARNING, returnCode = 9))
  }
  if (is.null(parameters$variance$scale) ||
        is.null(parameters$variance$shape) ||
        is.null(parameters$variance$orientation)) 
    stop("variance parameters are missing")
  temp <- .Fortran("eseev",
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$scale),
                   as.double(parameters$variance$shape),
                   as.double(aperm(parameters$variance$orientation,c(2,1,3))),
                   as.double(-1),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(-1),
                   double(p),
                   double(p),
                   double(1),
                   double(n * G),
                   PACKAGE = "mclust")[13:14]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, G)
  WARNING <- NULL
  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 = "EEV",
            WARNING = WARNING, returnCode = ret)
}

emEEV <- function(data, parameters, prior = NULL, control = emControl(), 
                  warn = NULL, ...)
{
  z <- estepEEV(data, parameters = parameters, warn = warn)$z  
  meEEV(data, z = z, prior = prior, control = control, 
        Vinv = parameters$Vinv, warn = warn)
}

estepEEV <- 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 = "EEV", 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) ||
        is.null(parameters$variance$orientation)) 
    stop("variance parameters are missing")
  temp <- .Fortran("eseev",
                   as.double(data),
                   as.double(mu),
                   as.double(parameters$variance$scale),
                   as.double(parameters$variance$shape),
                   as.double(aperm(parameters$variance$orientation,c(2,1,3))),
                   as.double(pro),
                   as.integer(n),
                   as.integer(p),
                   as.integer(G),
                   as.double(if (is.null(Vinv)) -1 else Vinv),
                   double(p),
                   double(p),
                   double(1),
                   double(n * K),
                   PACKAGE = "mclust")[13:14]
  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 = "EEV", n = n, d = p, G = G, 
                 z = z, parameters = parameters, loglik = loglik),
            WARNING = WARNING, returnCode = ret)
}

meEEV <- 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) || length(dimdat[dimdat > 1]) == 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 = "EEV", d = p, G = G, 
                     scale = NA, shape = rep(NA,p), orientation = array(NA,c(p,p,G)))
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance, Vinv=Vinv)
    return(structure(list(modelName="EEV", 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")
  lwork <- max(3 * min(n, p) + max(n, p), 5 * min(n, p))
  storage.mode(z) <- "double"
  if(is.null(prior)) {
    temp <- .Fortran("meeev",
                     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),
                     as.integer(lwork),
                     double(p * G),
                     double(1),
                     double(p),
                     double(p * p * G),
                     double(K),
                     double(lwork),
                     double(p),
                     PACKAGE = "mclust")[7:16]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEV"), 
                                                 prior[names(prior) !="functionName"]))
    temp <- .Fortran("meeevp",
                     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),
                     as.integer(lwork),
                     double(p * G),
                     double(1),
                     double(p),
                     double(p * p * G),
                     double(K),
                     double(lwork),
                     double(p),
                     PACKAGE = "mclust")[11:20]
  }
  z <- temp[[1]]
  its <- temp[[2]]
  err <- temp[[3]]
  loglik <- temp[[4]]
  lapackSVDinfo <- temp[[5]]
  mu <- matrix(temp[[6]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  scale <- temp[[7]]
  shape <- temp[[8]]
  O <- aperm(array(temp[[9]], c(p, p, G)),c(2,1,3))
  pro <- temp[[10]]
  WARNING <- NULL
  if(lapackSVDinfo) {
    if(lapackSVDinfo > 0) {
      WARNING <- "LAPACK DGESVD fails to converge"
    }
    else {
      WARNING <- "input error for LAPACK DGESVD"
    }
    z[] <- O[] <- shape[] <- NA
    scale <- loglik <- NA
    sigma <- array(NA, c(p, p, G))
    ret <- -9
  }
  else if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "singular covariance"
    if(warn) warning(WARNING)
    shape[] <- NA
    mu[] <- pro[] <- z[] <- loglik <- NA
    sigma <- array(NA, c(p, p, G))
    ret <- -1
  }
  else if(loglik <  - signif(.Machine$double.xmax, 6)) {
    if(control$equalPro) {
      WARNING <- "a 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 <- NA
    sigma <- array(NA, c(p, p, G))
    ret <- if(control$equalPro) -2 else -3
  }
  else {
    sigma <- scale * shapeO(shape, O, transpose = FALSE)
    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(O) <- list(dimnames(data)[[2]], dimnames(data)[[2]], 
                      NULL)
  ## Sigma = scale * O %*% diag(shape) %*% t(O)
  variance <- list(modelName = "EEV", d = p, G = G, sigma = sigma,
                   scale = scale, shape = shape, orientation = O) 
  parameters <- list(pro=pro, mean=mu, variance=variance, Vinv=Vinv) 
  structure(list(modelName = "EEV", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters, control = control,
                 loglik = loglik),
            info = info, WARNING = WARNING, returnCode = ret)
}

mstepEEV <- function(data, z, prior = NULL, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options("warn")
  dimdat <- dim(data)
  oneD <- is.null(dimdat) || length(dimdat[dimdat > 1]) == 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 = "EEV", d = p, G = G, 
                     scale = NA, shape = rep(NA,p), orientation=array(NA,c(p,p,G))) 
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance)
    return(structure(list(modelName="EEV", prior=prior, n=n, d=p, 
                          G=G, z=z, parameters=parameters), 
                     WARNING = WARNING, returnCode = 9))
    
  }
  #  shape <- sqrt(rev(sort(shape/exp(sum(log(shape))/p))))
  if(any(is.na(z)) || any(z < 0) || any(z > 1)) stop(
    "improper specification of z")
  lwork <- max(3 * min(n, p) + max(n, p), 5 * min(n, p), G)
  if(is.null(prior)) {
    temp <- .Fortran("mseev",
                     as.double(data),
                     as.double(z),
                     as.integer(n),
                     as.integer(p),
                     as.integer(G),
                     double(lwork),
                     as.integer(lwork),
                     double(p * G),
                     double(1),
                     double(p),
                     double(p * p * G),
                     double(G),
                     PACKAGE = "mclust")[7:12]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EEV"), 
                                                 prior[names(prior) != "functionName"]))
    temp <- .Fortran("mseevp",
                     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(lwork),
                     as.integer(lwork),
                     double(p * G),
                     double(1),
                     double(p),
                     double(p * p * G),
                     double(G),
                     PACKAGE = "mclust")[11:16]
  }
  lapackSVDinfo <- temp[[1]]
  mu <- matrix(temp[[2]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  scale <- temp[[3]]
  shape <- temp[[4]]
  O <- aperm( array(temp[[5]], c(p, p, G)), c(2,1,3))
  pro <- temp[[6]]
  WARNING <- NULL
  if(lapackSVDinfo) {
    if(lapackSVDinfo > 0) {
      WARNING <- "LAPACK DGESVD fails to converge"
      if(warn) warning(WARNING)
      ret <- -4
    }
    else {
      WARNING <- "input error for LAPACK DGESVD"
      if(warn) warning(WARNING)
      ret <- -5
    }
    O[] <- shape[] <- scale <- NA
    sigma <- array(NA, c(p, p, G))
  }
  else if(any(c(abs(scale), shape) > signif(.Machine$double.xmax, 6))) {
    WARNING <- "cannot compute M-step"
    if(warn) warning(WARNING)
    mu[] <- pro[] <- scale <- O[] <- shape[] <- NA
    sigma <- array(NA, c(p, p, G))
    ret <- -1
  }
  else {
    sigma <- scale * shapeO(shape, O, transpose = FALSE)
    ret <- 0
  }
  dimnames(z) <- list(dimnames(data)[[1]], NULL)
  dimnames(mu) <- list(dimnames(data)[[2]], NULL)
  dimnames(sigma) <- dimnames(O) <- 
    list(dimnames(data)[[2]], dimnames(data)[[2]], NULL)
  variance <- list(modelName = "EEV", d = p, G = G, sigma = sigma,
                   scale = scale, shape = shape, orientation = O)
  parameters <- list(pro=pro, mean=mu, variance=variance)
  structure(list(modelName = "EEV", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters), 
            WARNING = WARNING, returnCode = ret)
}

simEEV <- 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 = "EEV"))
  }
  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")
  sss <- sqrt(parameters$variance$scale * shape)
  for(k in 1:G) {
    m <- ctabel[k]
    cholSigma <- t(parameters$variance$orientation[,  , k]) * sss
    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 = "EEV")
}

emEII <- function(data, parameters, prior = NULL, control = emControl(), 
                  warn = NULL, ...)
{
  z <- estepEII(data, parameters = parameters, warn = warn)$z
  meEII(data, z = z, prior = prior, control = control, 
        Vinv = parameters$Vinv, warn = warn)
}

estepEII <- 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)
  p <- ncol(data)
  n <- nrow(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 = "EII", n=n, d=p, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     WARNING = WARNING, returnCode = 9))
  }
  sigmasq <- parameters$variance$sigmasq
  if(is.null(sigmasq))
    if(warn) warning("variance parameters are missing")
  if(sigmasq < 0)
    stop("sigma-squared is negative")
  if(!sigmasq) {
    WARNING <- "sigma-squared vanishes"
    if(warn) warning(WARNING)
    z <- matrix(as.double(NA),n,K)
    dimnames(z) <- list(dimnames(data)[[1]], NULL)
    return(structure(list(modelName = "EII", n=n, d=p, G=G, z=z,
                          parameters=parameters, loglik=NA), 
                     WARNING = WARNING, returnCode = -1))
  }
  temp <- .Fortran("eseii",
                   as.double(data),
                   as.double(mu),
                   as.double(sigmasq),
                   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")[9:10]
  loglik <- temp[[1]]
  z <- matrix(temp[[2]], n, K)
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6)) {
    WARNING <- "sigma-squared falls below threshold"
    if(warn) warning(WARNING)
    z[] <- loglik <- NA
    ret <- -1
  }
  else ret <- 0
  dimnames(z) <- list(dimnames(data)[[1]],NULL)
  structure(list(modelName = "EII", n = n, d = p, G = G, 
                 z = z, parameters = parameters, loglik = loglik),
            WARNING = WARNING, returnCode = ret)
}

meEII <- 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) || length(dimdat[dimdat > 1]) == 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]
  # number of groups
  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 = "EII", d = p, G = G, sigmasq = NA)
    parameters <- list(pro=rep(NA,G), mean=matrix(as.double(NA),p,G), 
                       variance=variance, Vinv=Vinv)
    return(structure(list(modelName="EII", 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("meeii",
                     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(K),
                     PACKAGE = "mclust")[7:13]
  }
  else {
    priorParams <- do.call(prior$functionName, c(list(data = 
                                                        data, G = G, modelName = "EII"), prior[names(prior) !=
                                                                                                 "functionName"]))
    temp <- .Fortran("meeiip",
                     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(K),
                     PACKAGE = "mclust")[c(11:17, 10)]
  }
  mu <- matrix(temp[[5]], p, G)
  dimnames(mu) <- list(NULL, as.character(1:G))
  z <- temp[[1]]
  its <- temp[[2]]
  err <- temp[[3]]
  loglik <- temp[[4]]
  sigmasq <- temp[[6]]
  Sigma <- diag(rep(sigmasq, p))
  pro <- temp[[7]]
  WARNING <- NULL
  if(loglik > signif(.Machine$double.xmax, 6) || 
       sigmasq <= max(control$eps,0)) {
    WARNING <- "sigma-squared falls below threshold"
    if(warn) warning(WARNING)
    mu[] <- pro[] <- sigmasq <- z[] <- loglik <- NA
    sigma <- array(NA, c(p, p, G))
    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[] <- sigmasq <- z[] <- loglik <- NA
    sigma <- array(NA, c(p, p, G))
    ret <- if(control$equalPro) -2 else -3
  }
  else {
    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) <- list(dimnames(data)[[2]], dimnames(data)[[2]])
  dimnames(sigma) <- list(dimnames(data)[[2]], dimnames(data)[[2]],
                          NULL) 
  variance <- list(modelName = "EII", d = p, G = G, sigma = sigma, 
                   Sigma = Sigma, sigmasq = sigmasq, scale = sigmasq)
  parameters <- list(pro=pro, mean=mu, variance = variance, Vinv=Vinv)
  structure(list(modelName = "EII", prior = prior, n = n, d = p, G = G, 
                 z = z, parameters = parameters, control = control, 
                 loglik = loglik),
            info = info, WARNING = WARNING, returnCode = ret)
}

mstepEII <- function(data, z, prior = NULL, warn = NULL, ...)
{
  if(is.null(warn)) warn <- mclust.options