R/Plotting_Functions.R

Defines functions .bubble .barcode_panel .barchart_panel .bar_code .violin_panel .scatter_panel .mosaic_panel .hc_reorder .reorder_hc .is_cols .heat_legend .diag_panel .density_panel .contour_panel .boxplot_panel .draw_axis .mvn2D_panel .grid_2 .grid_1 plot.MoEClust MoE_Uncertainty.MoEClust MoE_Uncertainty MoE_Similarity.MoEClust MoE_Similarity MoE_plotLogLik.MoEClust MoE_plotLogLik MoE_plotCrit.MoEClust MoE_plotCrit MoE_plotGate.MoEClust MoE_plotGate MoE_gpairs.MoEClust MoE_gpairs

Documented in MoE_gpairs MoE_plotCrit MoE_plotGate MoE_plotLogLik MoE_Similarity MoE_Uncertainty plot.MoEClust

#' Generalised Pairs Plots for MoEClust Mixture Models
#'
#' Produces a matrix of plots showing pairwise relationships between continuous response variables and continuous/categorical/logical/ordinal associated covariates, as well as the clustering achieved, according to fitted MoEClust mixture models.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param subset A list giving named arguments for producing only a subset of panels:
#' \describe{
#' \item{\code{show.map}}{Logical indicating whether to show panels involving the MAP classification (defaults to \code{TRUE}, unless there is only one component, in which case the MAP classification is never plotted.).}
#' \item{\code{data.ind}}{For subsetting response variables: a vector of column indices corresponding to the variables in the columns of \code{res$data} which should be shown. Defaults to all. Can be \code{0}, in order to suppress plotting the response variables.}
#' \item{\code{cov.ind}}{For subsetting covariates: a vector of column indices corresponding to the covariates in the columns \code{res$net.covs} which should be shown. Defaults to all. Can be \code{0}, in order to suppress plotting the covariates.}
#' }
#' The result of the subsetting must include at least two variables, whether they be the MAP classification, a response variable, or a covariate, in order to be valid for plotting purposes. The arguments \code{data.ind} and \code{cov.ind} can also be used to simply reorder the panels, without actually subsetting.
#' @param response.type The type of plot desired for the scatterplots comparing continuous response variables. Defaults to \code{"points"}. See \code{scatter.pars} below.
#'
#' Points can also be sized according to their associated clustering uncertainty with the option \code{"uncertainty"}. In doing so, the transparency of the points will also be proportional to their clustering uncertainty, provided the device supports transparency. See also \code{\link{MoE_Uncertainty}} for an alternative means of visualising observation-specific cluster uncertainties (especially for univariate data). See \code{scatter.pars} below, and note that models fitted via the \code{"CEM"} algorithm will have no associated clustering uncertainty.
#'
#' Alternatively, the bivariate \code{"density"} contours can be displayed (see \code{density.pars}), provided there is at least one Gaussian component in the model. Caution is advised when producing density plots for models with covariates in the expert network; the required number of evaluations of the (multivariate) Gaussian density for each panel (\code{res$G * prod(density.pars$grid.size)}) increases by a factor of \code{res$n}, thus plotting may be slow (particularly for large data sets). See \code{density.pars} below.
#' @param scatter.type A vector of length 2 (or 1) giving the plot type for the upper and lower triangular portions of the plot, respectively, pertaining to the associated covariates. Defaults to \code{"lm"} for covariate vs. response panels and \code{"points"} otherwise. Only relevant for models with continuous covariates in the gating &/or expert network. \code{"ci"} and \code{"lm"} type plots are only produced for plots pairing covariates with response, and never response vs. response or covariate vs. covariate. Note that lines &/or confidence intervals will only be drawn for continuous covariates included in the expert network; to include covariates included only in the gating network also, the options \code{"lm2"} or \code{"ci2"} can be used but this is not generally advisable. See \code{scatter.pars} below.
#' @param conditional A vector of length 2 (or 1) giving the plot type for the upper and lower triangular portions of the plot, respectively, for plots involving a mix of categorical and continuous variables. Defaults to \code{"stripplot"} in the upper triangle and \code{"boxplot"} in the lower triangle (see \code{\link[lattice]{panel.stripplot}} and \code{\link[lattice]{panel.bwplot}}). \code{"violin"} and \code{"barcode"} plots can also be produced. Only relevant for models with categorical covariates in the gating &/or expert network, unless \code{show.MAP} is \code{TRUE}. Comparisons of two categorical variables (which can only ever be covariates or the MAP classification) are always displayed via mosaic plots (see \code{\link[vcd]{strucplot}}). 
#' 
#' All \code{conditional} panel types can be customised further; see \code{stripplot.pars}, \code{boxplot.pars} (for both \code{"boxplot"} and \code{"violin"} plots), \code{barcode.pars}, and \code{mosaic.pars} below. Note that when \code{conditional} is of length 1, that plot type will be used in \emph{both} the upper and lower triangular portions of the plot, where relevant.
#' @param addEllipses Controls whether to add MVN ellipses with axes corresponding to the within-cluster covariances for the response data. The options \code{"inner"} and \code{"outer"} (the default) will colour the axes or the perimeter of those ellipses, respectively, according to the cluster they represent (according to \code{scatter.pars$eci.col}). The option \code{"both"} will obviously colour both the axes and the perimeter. The \code{"yes"} or \code{"no"} options merely govern whether the ellipses are drawn, i.e. \code{"yes"} draws ellipses without any colouring. Ellipses are only ever drawn for multivariate data, and only when \code{response.type} is \code{"points"} or \code{"uncertainty"}.
#'
#' Ellipses are centered on the posterior mean of the fitted values when there are expert network covariates, otherwise on the posterior mean of the response variables. In the presence of expert network covariates, the component-specific covariance matrices are also (by default, via the argument \code{expert.covar} below) modified for plotting purposes via the function \code{\link{expert_covar}}, in order to account for the extra variability of the means, usually resulting in bigger shapes & sizes for the MVN ellipses.
#' @param expert.covar Logical (defaults to \code{TRUE}) governing whether the extra variability in the component means is added to the MVN ellipses corresponding to the component covariance matrices in the presence of expert network covariates. See the function \code{\link{expert_covar}}. Only relevant when \code{response.type} is \code{"points"} or \code{"uncertainty"} when \code{addEllipses} is invoked accordingly, and/or \code{diag.pars$show.dens=TRUE} (see below), and only relevant for models with expert network covariates.
#' @param border.col A vector of length 5 (or 1) containing \emph{border} colours for plots against the MAP classification, response vs. response, covariate vs. response, response vs. covariate, and covariate vs. covariate panels, respectively.
#'
#' Defaults to \code{c("purple", "black", "brown", "brown", "navy")}.
#' @param bg.col A vector of length 5 (or 1) containing \emph{background} colours for plots against the MAP classification, response vs. response, covariate vs. response, response vs. covariate, and covariate vs. covariate panels, respectively.
#'
#' Defaults to \code{c("cornsilk", "white", "palegoldenrod", "palegoldenrod", "cornsilk")}.
#' @param outer.margins A list of length 4 with units as components named bottom, left, top, and right, giving the outer margins; the defaults uses two lines of text. A vector of length 4 with units (ordered properly) will work, as will a vector of length 4 with numeric variables (interpreted as lines).
#' @param outer.labels The default is \code{NULL}, for alternating labels around the perimeter. If \code{"all"}, all labels are printed, and if \code{"none"}, no labels are printed.
#' @param outer.rot A 2-vector (\code{x}, \code{y}) rotating the top/bottom outer labels \code{x} degrees and the left/right outer labels \code{y} degrees. Only works for categorical labels of boxplot and mosaic panels. Defaults to \code{c(0, 90)}.
#' @param gap The gap between the tiles; defaulting to \code{0.05} of the width of a tile.
#' @param buffer The fraction by which to expand the range of quantitative variables to provide plots that will not truncate plotting symbols. Defaults to \code{0.025}, i.e. 2.5 percent of the range. Particularly useful when ellipses are drawn (see \code{addEllipses}) to ensure ellipses are visible in full.
#' @param uncert.cov A logical indicating whether the expansion factor for points on plots involving covariates should also be modified when \code{response.type="uncertainty"}. Defaults to \code{FALSE}, and only relevant for scatterplot and strip plot panels. When \code{TRUE}, \code{scatter.pars$uncert.pch} is invoked as the plotting symbols for covariate-related scatterplot and strip plot panels, otherwise \code{scatter.pars$scat.pch} and \code{stripplot.pars$strip.pch} is invoked for such panels.
#' @param scatter.pars A list supplying select parameters for the continuous vs. continuous scatterplots.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(scat.pch=res$classification, uncert.pch=19,
#'      scat.col=res$classification, scat.size=unit(0.25, "char"), 
#'      eci.col=1:res$G, noise.size=unit(0.2, "char")),}
#' where \code{scat.pch}, \code{scat.col}, and \code{scat.size} give the plotting symbols, colours, and sizes of the points in scatterplot panels, respectively. Note that \code{eci.col} gives both a) the colour of the fitted lines &/or confidence intervals for expert-related panels when \code{scatter.type} is one of \code{"ci"} or \code{"lm"} and b) the colour of the ellipses (if any) when \code{addEllipses} is one of \code{"outer"}, \code{"inner"}, or \code{"both"} and the response data is multivariate. Note that \code{eci.col} will inherit a suitable default from \code{scat.col} instead if the latter is supplied but the former is not. 
#' 
#' Note also that \code{scat.size} will be modified on an observation-by-observation level when \code{response.type} is \code{"uncertainty"}. Furthermore, note that the behaviour for plotting symbols when \code{response.type="uncertainty"} changes compared to \code{response.type="points"} depending on the value of the \code{uncert.cov} argument above. \code{uncert.pch} gives the plotting symbol used for all scatterplot (and strip plot) panels when \code{response.type="uncertainty"} and \code{uncert.cov} is \code{TRUE}. However, when \code{uncert.cov} is \code{FALSE}, \code{scat.pch} is invoked for scatterplots involving covariates and \code{uncert.pch} is used for panels involving only response variables. Finally, \code{noise.size} can be used to modify \code{scat.size} for observations assigned to the noise component (if any), but only when \code{response.type="points"}.
#' @param density.pars A list supplying select parameters for visualising the bivariate density contours, only when \code{response.type} is \code{"density"}.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(grid.size=c(100, 100), dcol="grey50",
#'      nlevels=11, show.labels=TRUE, label.style="mixed"),}
#' where \code{grid.size} is a vector of length two giving the number of points in the x & y direction of the grid over which the density is evaluated, respectively (though \code{density.pars$grid.size} can also be supplied as a scalar, which will be automatically recycled to a vector of length 2), and \code{dcol} is either a single colour or a vector of length \code{nlevels} colours (although note that \code{dcol}, when \emph{not} specified, will be adjusted for transparency). Finally, \code{label.style} can take the values \code{"mixed"}, \code{"flat"}, or \code{"align"}. Note that \code{density.pars$grid.size[1]} is also relevant when \code{diag.pars$show.dens=TRUE} (see below).
#' @param stripplot.pars A list supplying select parameters for continuous vs. categorical panels when one or both of the entries of \code{conditional} is \code{"stripplot"}.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(strip.pch=res$classification, strip.size=unit(0.5, "char"),
#'      strip.col=res$classification, jitter=TRUE, size.noise=unit(0.4, "char")),}
#' where \code{strip.size} and \code{size.noise} retain the definitions for the similar arguments under \code{scatter.pars} above. However, \code{stripplot.pars$size.noise} is invoked regardless of the \code{response.type} (in contrast to \code{scatter.pars$noise.size}). Notably, \code{strip.col} will inherit a suitable default from \code{scatter.pars$scat.col} if the latter is supplied but the former is not. Note also that the \code{strip.pch} default is modified to \code{scatter.pars$uncert.pch} if \code{uncert.cov} is \code{TRUE}.
#' @param boxplot.pars A list supplying select parameters for continuous vs. categorical panels when one or both of the entries of \code{conditional} is \code{"boxplot"} or \code{"violin"}.
#' 
#' \code{NULL} is equivalent to:
#' \preformatted{list(box.pch="|", box.col="black", varwidth=FALSE,
#'      notch=FALSE, notch.frac=0.5, box.fill=1:res$G).}
#' All of the above are relevant for \code{"boxplot"} panels, are passed to \code{\link[lattice]{panel.bwplot}} when producing boxplots, and retain the same definitions as the similarly named arguments therein. However, only \code{box.col}, \code{varwidth}, and \code{box.fill} are relevant for \code{"violin"} panels, and in both cases \code{box.fill} is only invoked for panels where the categorical variable is the MAP classification (i.e. when \code{isTRUE(subset$show.map)}). See \code{diag.pars$hist.color} for controlling the colours of non-MAP-related boxplot/violin panels. Notably, \code{box.fill} will inherit a suitable default from \code{scatter.pars$scat.col} if the latter is supplied but the former is not.
#' @param barcode.pars A list supplying select parameters for continuous vs. categorical panels when one or both of the entries of \code{conditional} is \code{"barcode"}. See the help file for \code{barcode::barcode}.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(bar.col=res$G:1, nint=0, ptsize=unit(0.25, "char"), 
#'      ptpch=1, bcspace=NULL, use.points=FALSE),}
#' where \code{bar.col} is only invoked for panels where the categorical variable is the MAP classification (i.e. when \code{isTRUE(subset$show.map)}) if it is of length greater than 1, otherwise it is used for all relevant panels. See \code{diag.pars$hist.color} for controlling the colours of non-MAP-related barcode panels. Notably, \code{bar.col} will inherit a suitable default from \code{scatter.pars$scat.col} if the latter is supplied but the former is not.
#' @param mosaic.pars A list supplying select parameters for categorical vs. categorical panels (if any). 
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(shade=NULL, gp_labels=grid::gpar(fontsize=9), 
#'      gp_args=list(), gp=list(), mfill=TRUE, mcol=1:res$G).}
#' The current default arguments and values thereof are passed through to \code{\link[vcd]{strucplot}} for producing mosaic tiles. When \code{shade} is not \code{FALSE}, \code{mfill} is a logical which governs the colouring scheme for panels (if any) involving the MAP classification. When \code{mfill} is \code{TRUE} (the default), \code{gp} is invoked here in such a way that tiles will inherit appropriate interior colours via \code{gp$fill} from \code{mcol} and a \code{"black"} outer colour via \code{gp$col}. When \code{mfill} is \code{FALSE}, or the panel involves two categorical covariates, the outer colours are inherited from \code{mcol} and the interior fill colour is inherited from \code{bg.col}. See \code{diag.pars$hist.color} for controlling the interior fill colour of non-MAP-related mosaic panels. Notably, \code{mcol} will inherit a suitable default from \code{scatter.pars$scat.col} if the latter is supplied but the former is not.
#' @param axis.pars A list supplying select parameters for controlling the axes.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(n.ticks=5, axis.fontsize=9).}
#' The argument \code{n.ticks} will be overwritten for categorical variables with fewer than 5 levels.
#' @param diag.pars A list supplying select parameters for panels along the diagonal.
#'
#' \code{NULL} is equivalent to:
#' \preformatted{list(diag.fontsize=9, show.hist=TRUE, show.dens=FALSE,
#'      diagonal=TRUE, hist.color=hist.color, show.counts=TRUE),}
#' where \code{hist.color} is a vector of length 4, giving the colours for the response variables, gating covariates, expert covariates, and covariates entering both networks, respectively. By default, diagonal panels for response variables are \code{ifelse(diag.pars$show.dens, "white", "black")} and covariates of any kind are \code{"dimgrey"}. \code{hist.color} also governs the outer colour for mosaic panels and the fill colour for boxplot, violin, and barcode panels (except for those involving the MAP classification). However, in the case of response vs. (categorical) covariates boxplots and violin plots, the fill colour is always \code{"white"}. The MAP classification is always coloured by cluster membership, by default. The argument \code{show.counts} is only relevant for categorical variables.
#'
#' The argument \code{show.dens} toggles whether parametric density estimates are drawn over the diagonal panels for each response variable. When \code{show.dens=TRUE}, the component densities are shown via thin lines, with colours given by \code{scatter.pars$scat.col}, while a thick \code{"black"} line is used for the overall mixture density. This argument can be used with or without \code{show.hist} also being \code{TRUE}, though density curves will appear bigger when \code{show.hist=FALSE}. Note that \code{show.dens=TRUE} is also affected by the \code{expert.covar} argument above. Finally, the grid size when \code{show.dens=TRUE} is given by \code{max(res$n, density.pars$grid.size[1])}.
#' 
#' When \code{diagonal=TRUE} (the default), the diagonal from the top left to the bottom right is used for displaying the marginal distributions of variables (via histograms, with or without overlaid density estimates, or barplots, as appropriate). Specifying \code{diagonal=FALSE} will place the diagonal running from the top right down to the bottom left.
#' @param ... Catches unused arguments. Alternatively, named arguments can be passed directly here to any/all of \code{scatter.pars}, \code{stripplot.pars}, \code{boxplot.pars}, \code{barcode.pars}, \code{mosaic.pars}, \code{axis.pars}, and \code{diag.pars}.
#'
#' @importFrom lattice "current.panel.limits" "llines" "panel.abline" "panel.bwplot" "panel.histogram" "panel.lines" "panel.points" "panel.rect" "panel.stripplot" "panel.text" "panel.violin" "trellis.grobname" "trellis.par.get" "trellis.par.set"
#' @importFrom matrixStats "colMeans2" "rowLogSumExps"
#' @importFrom mclust "mclust.options" "sigma2decomp"
#' @importFrom vcd "strucplot"
#'
#' @return A generalised pairs plot showing all pairwise relationships between clustered response variables and associated gating &/or expert network continuous &/or categorical variables, coloured according to the MAP classification, with the marginal distributions of each variable along the diagonal.
#' @note For \code{MoEClust} models with more than one expert network covariate, fitted lines produced in continuous covariate vs. continuous response scatterplots via \code{scatter.type="lm"} or \code{scatter.type="ci"} will \strong{NOT} correspond to the coefficients in the expert network (\code{res$expert}).
#'
#' \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_gpairs}} which accepts the default arguments, and also produces other types of plots. Caution is advised producing generalised pairs plots when the dimension of the data is large.
#' 
#' Finally, note that all colour-related defaults in \code{scatter.pars}, \code{stripplot.pars}, \code{barcode.pars}, and \code{mosaic.pars} above assume a specific colour-palette (see \code{\link[mclust]{mclust.options}("classPlotColors")}). Thus, for instance, specifying \code{scatter.pars$scat.col=res$classification} will produce different results compared to leaving this argument unspecified. This is especially true for models with a noise component, for which the default is handled quite differently (for one thing, \code{res$G} is the number of \emph{non-noise} components). Similarly, all \code{pch}-related defaults in \code{scatter.pars} and \code{stripplot.pars} above assume a specific set of plotting symbols also (see \code{\link[mclust]{mclust.options}("classPlotSymbols")}). Generally, all colour and symbol related arguments are strongly recommended to be left at their default values, unless being supplied as a single character string, e.g. \code{"black"} for colours. To help in this regard, colour-related arguments sensibly inherent their default from \code{scatter.pars$scat.col} if that is supplied and the argument in question is not.
#' @export
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @references Murphy, K. and Murphy, T. B. (2020). Gaussian parsimonious clustering models with covariates and a noise component. \emph{Advances in Data Analysis and Classification}, 14(2): 293-325. <\doi{10.1007/s11634-019-00373-8}>.
#'
#' Emerson, J. W., Green, W. A., Schloerke, B., Crowley, J., Cook, D., Hofmann, H. and Wickham, H. (2013). The generalized pairs plot. \emph{Journal of Computational and Graphical Statistics}, 22(1): 79-91.
#' @seealso \code{\link{MoE_clust}}, \code{\link{MoE_stepwise}}, \code{\link{plot.MoEClust}}, \code{\link{MoE_Uncertainty}}, \code{\link{expert_covar}}, \code{\link[lattice]{panel.stripplot}}, \code{\link[lattice]{panel.bwplot}}, \code{\link[lattice]{panel.violin}}, \code{\link[vcd]{strucplot}}, \code{\link[mclust]{mclust.options}}
#' @keywords plotting
#' @usage
#' MoE_gpairs(res,
#'            response.type = c("points", "uncertainty", "density"),
#'            subset = list(...),
#'            scatter.type = c("lm", "points"),
#'            conditional = c("stripplot", "boxplot"),
#'            addEllipses = c("outer", "yes", "no", "inner", "both"),
#'            expert.covar = TRUE,
#'            border.col = c("purple", "black", "brown", "brown", "navy"),
#'            bg.col = c("cornsilk", "white", "palegoldenrod", "palegoldenrod", "cornsilk"),
#'            outer.margins = list(bottom = grid::unit(2, "lines"),
#'                                 left = grid::unit(2, "lines"),
#'                                 top = grid::unit(2, "lines"),
#'                                 right = grid::unit(2, "lines")),
#'            outer.labels = NULL,
#'            outer.rot = c(0, 90),
#'            gap = 0.05,
#'            buffer = 0.025,
#'            uncert.cov = FALSE,
#'            scatter.pars = list(...),
#'            density.pars = list(...),
#'            stripplot.pars = list(...),
#'            boxplot.pars = list(...),
#'            barcode.pars = list(...),
#'            mosaic.pars = list(...),
#'            axis.pars = list(...),
#'            diag.pars = list(...),
#'            ...)
#' @examples
#' \donttest{data(ais)
#' res   <- MoE_clust(ais[,3:7], G=2, gating= ~ BMI, expert= ~ sex,
#'                    network.data=ais, modelNames="EVE")
#' MoE_gpairs(res)
#'
#' # Produce the same plot, but with a violin plot in the lower triangle.
#' # Colour the outline of the mosaic tiles rather than the interior using mfill.
#' # Size points in the response vs. response panels by their clustering uncertainty.
#' 
#' MoE_gpairs(res, conditional=c("stripplot", "violin"),
#'            mfill=FALSE, response.type="uncertainty")
#'
#' # Instead show the bivariate density contours of the response variables (without labels).
#' # (Plotting may be slow when response.type="density" for models with expert covariates.)
#' # Use different colours for histograms of covariates in the gating/expert/both networks.
#' # Also use different colours for response vs. covariate & covariate vs. response panels.
#' 
#' MoE_gpairs(res, response.type="density", show.labels=FALSE,
#'            hist.color=c("black", "cyan", "hotpink", "chartreuse"),
#'            bg.col=c("whitesmoke", "white", "mintcream", "mintcream", "floralwhite"))
#'            
#' # Examine the effect of the expert.covar argument in conjunction with show.dens
#' MoE_gpairs(res, cov.ind=0, expert.covar=TRUE, 
#'            show.dens=TRUE, show.hist=FALSE, grid.size=1000)
#' MoE_gpairs(res, cov.ind=0, expert.covar=FALSE, 
#'            show.dens=TRUE, show.hist=FALSE, grid.size=1000)
#'            
#' # Produce a generalised pairs plot for a model with a noise component.
#' # Reorder the covariates and omit the variables "Hc" and "Hg".
#' # Use barcode plots for the categorical/continuous pairs.
#' # Magnify the size of scatter points assigned to the noise component.
#' 
#' resN  <- MoE_clust(ais[,3:7], G=2, gating= ~ SSF + Ht, expert= ~ sex,
#'                    network.data=ais, modelNames="EEE", tau0=0.1, noise.gate=FALSE)
#'                    
#' MoE_gpairs(resN, data.ind=c(1,2,5), cov.ind=c(3,1,2), 
#'            conditional="barcode", noise.size=grid::unit(0.5, "char"))}
MoE_gpairs          <- function(res, response.type = c("points", "uncertainty", "density"), subset = list(...), scatter.type = c("lm", "points"), conditional = c("stripplot", "boxplot"),
                                addEllipses = c("outer", "yes", "no", "inner", "both"), expert.covar = TRUE, border.col = c("purple", "black", "brown", "brown", "navy"), bg.col = c("cornsilk", "white", "palegoldenrod", "palegoldenrod", "cornsilk"),
                                outer.margins = list(bottom=grid::unit(2, "lines"), left=grid::unit(2, "lines"), top=grid::unit(2, "lines"), right=grid::unit(2, "lines")), outer.labels = NULL, outer.rot = c(0, 90), gap = 0.05, buffer = 0.025, uncert.cov = FALSE, 
                                scatter.pars = list(...), density.pars = list(...), stripplot.pars = list(...), boxplot.pars = list(...), barcode.pars = list(...), mosaic.pars = list(...), axis.pars = list(...), diag.pars = list(...), ...) {
    UseMethod("MoE_gpairs")
}

#' @method MoE_gpairs MoEClust
#' @export
MoE_gpairs.MoEClust <- function(res, response.type = c("points", "uncertainty", "density"), subset = list(...), scatter.type = c("lm", "points"), conditional = c("stripplot", "boxplot"),
                                addEllipses = c("outer", "yes", "no", "inner", "both"), expert.covar = TRUE, border.col = c("purple", "black", "brown", "brown", "navy"), bg.col = c("cornsilk", "white", "palegoldenrod", "palegoldenrod", "cornsilk"),
                                outer.margins = list(bottom=grid::unit(2, "lines"), left=grid::unit(2, "lines"), top=grid::unit(2, "lines"), right=grid::unit(2, "lines")), outer.labels = NULL, outer.rot = c(0, 90), gap = 0.05, buffer = 0.025, uncert.cov = FALSE, 
                                scatter.pars = list(...), density.pars = list(...), stripplot.pars = list(...), boxplot.pars = list(...), barcode.pars = list(...), mosaic.pars = list(...), axis.pars = list(...), diag.pars = list(...), ...) {

  res   <- if(inherits(res, "MoECompare")) res$optimal else res
  opar  <- suppressWarnings(graphics::par(no.readonly=TRUE))
  opar$new          <- FALSE
  on.exit(suppressWarnings(graphics::par(opar)))
  suppressWarnings(graphics::par(pty="m"))
  dat   <- res$data
  net   <- res$net.covs
  z     <- res$z
  claSS <- res$classification
  G     <- res$G
  Gseq  <- seq_len(G)
  both  <- attr(net, "Both")
  gate  <- setdiff(attr(net, "Gating"), both)
  expx  <- setdiff(attr(net, "Expert"), both)
  subset            <- if(!is.null(subset)          && inherits(subset,         "list")) subset
  scatter.pars      <- if(!is.null(scatter.pars)    && inherits(scatter.pars,   "list")) scatter.pars
  density.pars      <- if(!is.null(density.pars)    && inherits(density.pars,   "list")) density.pars
  stripplot.pars    <- if(!is.null(stripplot.pars)  && inherits(stripplot.pars, "list")) stripplot.pars
  boxplot.pars      <- if(!is.null(boxplot.pars)    && inherits(boxplot.pars,   "list")) boxplot.pars
  barcode.pars      <- if(!is.null(barcode.pars)    && inherits(barcode.pars,   "list")) barcode.pars
  mosaic.pars       <- if(!is.null(mosaic.pars)     && inherits(mosaic.pars,    "list")) mosaic.pars
  axis.pars         <- if(!is.null(axis.pars)       && inherits(axis.pars,      "list")) axis.pars
  diag.pars         <- if(!is.null(diag.pars)       && inherits(diag.pars,      "list")) diag.pars
  if(is.null(subset$show.map)) {
    subset$show.map <- (length(unique(claSS)) != 1) && ((G + !is.na(res$hypvol) > 1))
  } else {
    if(length(subset$show.map) > 1  ||
       !is.logical(subset$show.map))                  stop("'subset$show.map' should be a single logical indicator", call.=FALSE)
    if(length(unique(claSS))  == 1)       {           message("'show.map' set to FALSE as there is only one non-empty component\n")
      subset$show.map         <- FALSE
    }
  }
  if(is.null(subset$data.ind)) {
    subset$data.ind <- seq_len(ncol(dat))
  } else if(length(subset$data.ind) < 1  ||
            !all(is.numeric(subset$data.ind)) ||
 !all(subset$data.ind %in% c(0, seq_len(ncol(dat))))) stop("Invalid 'subset$data.ind'", call.=FALSE)
  if(is.null(subset$cov.ind))  {
    subset$cov.ind <- seq_len(ncol(net))
  } else if(length(subset$cov.ind)  < 1  ||
            !all(is.numeric(subset$cov.ind))  ||
 !all(subset$cov.ind  %in% c(0, seq_len(ncol(net))))) stop("Invalid 'subset$cov.ind'",  call.=FALSE)
  subset$data.ind  <- unique(subset$data.ind)
  subset$cov.ind   <- unique(subset$cov.ind)
  if((length(subset$cov.ind[subset$cov.ind     > 0]) +
     length(subset$data.ind[subset$data.ind    > 0]) +
  subset$show.map)  < 1)                              stop("Invalid subsetting", call.=FALSE)
  if((length(c(subset$data.ind,
     subset$cov.ind)) + subset$show.map)  < 1)        stop("Not enough columns to plot based on arguments supplied to 'subset'!", call.=FALSE)
  if((length(c(subset$data.ind,
     subset$cov.ind)) + subset$show.map) == 1)        warning("Try a different plotting method as there is only one panel to display!", call.=FALSE, immediate.=TRUE)
  dat   <- dat[,subset$data.ind, drop=FALSE]
  net   <- net[,subset$cov.ind,  drop=FALSE]
  dcol  <- ncol(dat)  + subset$show.map
  uni.c <- unique(claSS[claSS > 0])
  claSS <- factor(claSS)
  x     <- if(ncol(net) == 0) as.data.frame(dat)      else cbind(dat, net)
  x     <- if(subset$show.map) cbind(MAP  = claSS, x) else x
  clust <- as.character(claSS)
  zc    <- function(x) length(unique(x)) <= 1L
  saxzc <- vapply(x, zc, logical(1L))
  nrm   <- sum(saxzc, na.rm=TRUE)
  if(any(saxzc, na.rm=TRUE)) {                        warning(paste(nrm, "column", ifelse(nrm > 1, "s", ""), " with less than two distinct values eliminated\n"), call.=FALSE)
   dcol <- sum(which(saxzc) < dcol)
   x    <- x[,!saxzc]
  }
  N     <- ncol(x)
  Nseq  <- seq_len(N)
  both  <- which(names(net) %in% gsub("[[:space:]]", ".", both)) + dcol
  gate  <- which(names(net) %in% gsub("[[:space:]]", ".", gate)) + dcol
  expx  <- which(names(net) %in% gsub("[[:space:]]", ".", expx)) + dcol
  both[is.na(both)] <- 0L
  gate[is.na(gate)] <- 0L
  expx[is.na(expx)] <- 0L
  both  <- if(length(both) == 0L) 0L else both
  gate  <- if(length(gate) == 0L) 0L else gate
  expx  <- if(length(expx) == 0L) 0L else expx
  U     <- sort(unique(clust))
  L     <- length(U)
  noise <- any(clust == 0) || !is.na(res$hypvol)
  if(L  <= length(mclust.options("classPlotSymbols"))) {
    symbols   <- mclust.options("classPlotSymbols")
    if(noise) { 
      symbols <- c(symbols[-which(symbols == 4)][Gseq], 4)
    } else    {
      symbols <- symbols[Gseq]
    }
  } else if(L <= 9)  {
    symbols   <- as.character(seq_len(9L))
  } else if(L <= 26) {
    symbols   <- LETTERS
  } 
  if(L <= length(mclust.options("classPlotColors"))) {
    colors    <- mclust.options("classPlotColors")[Gseq]
    if(noise)  { 
      colors  <- c(colors[Gseq], "darkgrey")
    }
  } else if(length(colors) == 1) colors <- rep(colors, L)
  if(length(symbols) < L) {                           warning("More symbols needed to show classification\n", call.=FALSE)
    symbols <- rep(16, L)
  }
  if(length(colors)  < L) {                           warning("More colors needed to show classification\n",  call.=FALSE)
    colors <- rep("black", L)
  }
  scatter.type <- if(length(scatter.type) == 1L) rep(scatter.type, 2L) else scatter.type
  conditional  <- if(length(conditional)  == 1L) rep(conditional,  2L) else conditional
  bg.col       <- if(length(bg.col)       == 1L) rep(bg.col,       5L) else bg.col
  border.col   <- if(length(border.col)   == 1L) rep(border.col,   5L) else border.col
  if(length(bg.col)     != 5)                         stop("'bg.col' must be a vector of length 1 or 5 containing valid colours",     call.=FALSE)
  if(length(border.col) != 5)                         stop("'border.col' must be a vector of length 1 or 5 containing valid colours", call.=FALSE)
  if(!missing(response.type)    && (length(response.type) > 1 ||
     !is.character(response.type)))                   stop("'response.type' must be a single character string",     call.=FALSE)
  response.type         <- match.arg(response.type)
  if(G == 0    && response.type == "density")         stop("'response.type' cannot be \"density\" for models with only a noise component",       call.=FALSE)
  if(length(scatter.type)  != 2 ||
            !all(is.character(scatter.type)))         stop("'scatter.type' must be a character vector of length 2", call.=FALSE)
  if(length(conditional)   != 2 ||
            !all(is.character(conditional)))          stop("'conditional' must be a character vector of length 2",  call.=FALSE)
  if(!all(scatter.type %in% c("ci", "lm", "points",
                              "ci2", "lm2")))         stop("The entries of 'scatter.type' must be one of 'points', 'ci', 'lm', 'ci2', or 'lm2'", call.=FALSE)
  if(!all(conditional  %in% c("stripplot", "violin",
                             "boxplot", "barcode")))  stop("The entries of 'conditional' must be one of 'stripplot', 'boxplot', 'violin' or 'barcode'", call.=FALSE)
  if(!missing(addEllipses) && (length(addEllipses) > 1 ||
     !is.character(addEllipses)))                     stop("'addEllipses' must be a single character string",       call.=FALSE)
  addEllipses  <- match.arg(addEllipses)
  addEllipses  <- ifelse(G == 0,   "no",    addEllipses)
  drawEllipses <- addEllipses   != "no"
  colEllipses  <- addEllipses   != "yes" && drawEllipses
  
  upr.gate <- grepl("2", scatter.type[1L])
  low.gate <- grepl("2", scatter.type[2L])
  upr.exp  <- ifelse(upr.gate, substr(scatter.type[1L], 1L, nchar(scatter.type[1L]) - 1L), scatter.type[1L])
  low.exp  <- ifelse(low.gate, substr(scatter.type[2L], 1L, nchar(scatter.type[2L]) - 1L), scatter.type[2L])
  upr.cond <- conditional[1L]
  low.cond <- conditional[2L]
  if(!inherits(outer.margins, "list")) {
    if(length(outer.margins) == 4)     {
      outer.margins     <- if(inherits(outer.margins[1L], "units")) list(bottom=outer.margins[1L], left=outer.margins[2L], top=outer.margins[3L], right=outer.margins[4L]) else list(bottom=grid::unit(outer.margins[1L], "lines"), left=grid::unit(outer.margins[2L], "lines"), top=grid::unit(outer.margins[3L], "lines"), right=grid::unit(outer.margins[4L], "lines"))
    } else                                            stop("'outer.margins' are not valid", call.=FALSE)
  }
  if(is.null(outer.labels)) {
    lab1   <- switch(EXPR=names(x)[1L], MAP=1L, as.integer(res$d != 1) + 1L)
    lab2   <- switch(EXPR=names(x)[1L], MAP=2L, 1L)
    outer.labels$top    <- logical(N)
    outer.labels$top[seq(lab1,  N, by=2)] <- TRUE
    outer.labels$left   <- logical(N)
    outer.labels$left[seq(lab2, N, by=2)] <- TRUE
    outer.labels$right  <- !outer.labels$left
    outer.labels$bottom <- !outer.labels$top
  } else {
    if(pmatch(as.character(outer.labels),         "all", nomatch=FALSE)) {
      all.labeling      <- TRUE
    } else if(pmatch(as.character(outer.labels), "none", nomatch=FALSE)) {
      all.labeling      <- FALSE
    } else                                            stop("Invalid 'outer.labels'", call.=FALSE)
    outer.labels        <- NULL
    outer.labels$top    <-
    outer.labels$left   <-
    outer.labels$bottom <-
    outer.labels$right  <- rep(all.labeling, N)
  }
  if(length(outer.rot)  != 2    ||
    !all(is.numeric(outer.rot)) ||
     any(outer.rot       < 0))                        stop("Invalid 'outer.rot': must be a strictly non-negative numeric vector of length 2", call.=FALSE)
  claSS                 <- as.integer(levels(claSS))[claSS]
  if(any(C <- claSS == 0))       {
    claSS0              <- factor(claSS, levels=c(sort(levels(factor(claSS)))[-1L], 0L))
    claSS[which(C)]     <- G  + 1L
  } else claSS0         <- factor(claSS)
  x[,1L]                <- if(names(x)[1L]  == "MAP") claSS0 else x[,1L]
  if(length(gap)        != 1    || (!is.numeric(gap)    ||
     gap    < 0))                                     stop("'gap' must be single strictly non-negative number", call.=FALSE)
  if(length(buffer)     != 1    || (!is.numeric(buffer) ||
     buffer < 0))                                     stop("'buffer' must be single strictly non-negative number", call.=FALSE)
  if(is.null(scatter.pars$scat.pch))  {
    scatter.pars$scat.pch         <- symbols[claSS]
  }
  if(is.null(scatter.pars$uncert.pch))       {
    scatter.pars$uncert.pch       <- 19L
  }
  if(is.null(scatter.pars$scat.size)) {
    scatter.pars$size   <- grid::unit(0.25, "char")
  } else scatter.pars$size        <- scatter.pars$scat.size
  if(length(scatter.pars$size)   > 1 ||
       !inherits(scatter.pars$size,         "unit"))  stop("'scatter.pars$scat.size' must be a single item of class 'unit'", call.=FALSE)
  if(is.null(scatter.pars$noise.size)) {
    scatter.pars$noise.size       <- grid::unit(0.2, "char")
  } else scatter.pars$noise.size  <- scatter.pars$noise.size
  if(length(scatter.pars$noise.size)   > 1 ||
     !inherits(scatter.pars$noise.size,     "unit"))  stop("'scatter.pars$noise.size' must be a single item of class 'unit'", call.=FALSE)
  scat.null             <- is.null(scatter.pars$scat.col)
  if(isTRUE(scat.null))  {
    scatter.pars$col    <- colors[claSS]
  } else scatter.pars$col    <- scatter.pars$scat.col
  if(is.null(scatter.pars$eci.col))   {
    if(isFALSE(scat.null))    {
      scatter.pars$eci.col   <- unique(scatter.pars$col[claSS != G + 1])[stats::na.omit(match(Gseq, uni.c))]
    } else scatter.pars$eci.col      <- colors[Gseq]
  }
  if(length(scatter.pars$eci.col)    != G)  {
    if(length(scatter.pars$eci.col)  == 1)  {
      scatter.pars$eci.col   <- rep(scatter.pars$eci.col, G)
    } else                                            stop("'scatter.pars$eci.col' must be of length 1 or G", call.=FALSE)
  }
  if(response.type == "uncertainty")  {
    if(attr(res, "Algo")     == "CEM")                message("Model was fitted by CEM and has no clustering uncertainty\n")
    uncertainty         <- res$uncertainty
    uncertainty         <- res$uncertainty <- (uncertainty - min(uncertainty))/(diff(range(uncertainty)) + .Machine$double.eps)
    bubbleX             <- .bubble(uncertainty, cex=c(0.3, 2.8), alpha=c(0.3, 0.8))
    uncertainty         <- list()
    uncertainty$cex     <- bubbleX$cex
    uncertainty$col     <- if(grDevices::dev.capabilities()$semiTransparency) unname(mapply(grDevices::adjustcolor, col=scatter.pars$col, alpha.f=bubbleX$alpha)) else scatter.pars$col
  }
  if(is.null(density.pars$grid.size)) {
    density.pars$grid.size   <- c(100L, 100L)
  } else if(length(density.pars$grid.size)  == 1) {
    density.pars$grid.size   <- rep(density.pars$grid.size, 2L)
  }
  if(length(density.pars$grid.size)  != 2 || !all(is.numeric(density.pars$grid.size)) ||
            any(density.pars$grid.size       < 10))   stop("Invalid 'density.pars$grid.size'", call.=FALSE)
  if(is.null(density.pars$nlevels))   {
    density.pars$nlevels <- 11
  } else if(length(density.pars$nlevels)   > 1   ||
            !is.numeric(density.pars$nlevels)    ||
            density.pars$nlevels            <= 1)     stop("Invalid 'density.pars$nlevels'", call.=FALSE)
  if(is.null(density.pars$dcol))     {
    density.pars$dcol    <- if(grDevices::dev.capabilities()$semiTransparency) unname(mapply(grDevices::adjustcolor, col="grey50", alpha.f=seq(0.5, 1, length=density.pars$nlevels))) else "grey50"
  } else if(!is.element(length(density.pars$dcol),
            c(1, density.pars$nlevels))     ||
            !is.character(density.pars$dcol))         stop("Invalid 'density.pars$dcol", call.=FALSE)
  density.pars$dcol      <- if(length(density.pars$dcol) == 1) rep(density.pars$dcol, density.pars$nlevels) else density.pars$dcol
  if(is.null(density.pars$show.labels))  {
    density.pars$show.labels <- TRUE
  } else if(length(density.pars$show.labels) > 1 ||
            !is.logical(density.pars$show.labels))    stop("Invalid 'density.pars$show.labels", call.=FALSE)
  if(is.null(density.pars$label.style))  {
    density.pars$label.style <- "mixed"
  } else if(density.pars$show.labels &&
           (length(density.pars$label.style) > 1 ||
            !is.character(density.pars$label.style) ||
            !is.element(density.pars$label.style,
            c("mixed", "flat", "align"))))            stop("Invalid 'density.pars$label.style", call.=FALSE)
  if(is.null(axis.pars$n.ticks))     {
    axis.pars$n.ticks   <- 5
  }
  if(is.null(axis.pars$axis.fontsize))   {
    axis.pars$fontsize  <- 9
  } else axis.pars$fontsize       <- axis.pars$axis.fontsize
  if(axis.pars$n.ticks   < 3)        {                warning("Fewer than 3 axis ticks might cause problems\n", call.=FALSE)
    axis.pars$n.ticks   <- 3
  }
  diagonal              <- diag.pars$diagonal
  if(is.null(diagonal))    {
    diagonal            <- TRUE
  } else if(length(diagonal)  > 1 ||
            !is.logical(diagonal))                    stop("'diag.pars$diagonal' must be a single logical indicator", call.=FALSE)
  if(is.null(diag.pars$diag.fontsize))   {
    diag.pars$fontsize  <- 9
  } else diag.pars$fontsize       <- diag.pars$diag.fontsize
  if(is.null(diag.pars$show.hist))   {
    diag.pars$show.hist <- TRUE
  }
  if(is.null(diag.pars$show.dens))   {
    show.D <- FALSE
  } else {
    show.D <- diag.pars$show.dens
  }
  if(response.type == "density"   || show.D) {
    res$parameters$fits <- array(unlist(lapply(res$expert, "[[", "fitted.values")), dim=c(res$n, res$d, G))
    res$parameters$lpro <- log(if(isTRUE(attr(res, "Gating"))) colMeans2(res$parameters$pro) else res$parameters$pro)
  }
  if(isTRUE(show.D)) {
    diag.pars$grid.size <- density.pars$grid.size[1L]
  }
  res$parameters$varianceX        <- if(isTRUE(expert.covar)  &&
                                        attr(res, "Expert")   &&
                                        (isTRUE(drawEllipses) || 
                                        isTRUE(show.D))) suppressMessages(expert_covar(res, ...)) else res$parameters$variance
  if(is.null(diag.pars$hist.color))  {
   diag.pars$hist.color <- c("black", "dimgrey", "dimgrey", "dimgrey")
  } else {
   diag.pars$hist.color <- if(length(diag.pars$hist.color) == 1) rep(diag.pars$hist.color, 4L) else diag.pars$hist.color
   if(length(diag.pars$hist.color) != 4)              stop("'diag.pars$hist.color' must be a vector of length 1 or 4", call.=FALSE)
  }
  hist.col              <- diag.pars$hist.color
  diag.pars$hist.color  <- replace(Nseq,                 Nseq  <=  dcol, hist.col[1L])
  diag.pars$hist.color  <- replace(diag.pars$hist.color, Nseq %in% gate, hist.col[2L])
  diag.pars$hist.color  <- replace(diag.pars$hist.color, Nseq %in% expx, hist.col[3L])
  diag.pars$hist.color  <- replace(diag.pars$hist.color, Nseq %in% both, hist.col[4L])
  if(is.null(diag.pars$show.counts)) {
    diag.pars$show.counts  <- TRUE
  } else if(length(diag.pars$show.counts) > 1 ||
            !is.logical(diag.pars$show.counts))       stop("'diag.pars$show.counts' must be a single logical indicator", call.=FALSE)
  if(is.null(stripplot.pars$strip.pch))  {
    stripplot.pars$pch            <- if(response.type == "uncertainty" && uncert.cov) scatter.pars$uncert.pch else symbols[claSS]
  } else stripplot.pars$pch       <- if(response.type == "uncertainty" && uncert.cov) scatter.pars$uncert.pch else stripplot.pars$strip.pch
  if(length(uncert.cov)  > 1                  ||
     !is.logical(uncert.cov))                         stop("'uncert.cov' must be a single logical indicator", call.=FALSE)
  uncert.cov            <- uncert.cov   && response.type == "uncertainty"
  if(is.null(stripplot.pars$strip.size) ||
     any(names(list(...)) == "size")) {
    stripplot.pars$size <- if(response.type == "uncertainty" && uncert.cov) .bubble(res$uncertainty, cex=c(0.15, 1.4), alpha=c(0.3, 0.8))$cex else grid::unit(0.5, "char")
  } else stripplot.pars$size      <- if(response.type  == "uncertainty" && uncert.cov) .bubble(res$uncertainty, cex=c(0.15, 1.4), alpha=c(0.3, 0.8))$cex else stripplot.pars$strip.size
  if(any(response.type  != "uncertainty", !uncert.cov) &&
    (length(stripplot.pars$size) > 1    ||
            !inherits(stripplot.pars$size, "unit")))  stop("'stripplot.pars$strip.size' must be a single item of class 'unit'", call.=FALSE)
  if(is.null(stripplot.pars$size.noise)) {
    stripplot.pars$size.noise     <- grid::unit(0.4, "char")
  } else stripplot.pars$size.noise      <- stripplot.pars$size.noise
  if(length(stripplot.pars$size.noise)   > 1 ||
     !inherits(stripplot.pars$size.noise,  "unit"))   stop("'stripplot.pars$size.noise' must be a single item of class 'unit'", call.=FALSE)
  if(noise) {
    if(response.type    == "points")     {
      scatter.pars$size <- replace(rep(scatter.pars$size,   nrow(dat)), clust == 0, scatter.pars$noise.size)  
    }
    stripplot.pars$size <- replace(rep(stripplot.pars$size, nrow(dat)), clust == 0, stripplot.pars$size.noise)
  }
  if(is.null(stripplot.pars$strip.col))   {
    stripplot.pars$col  <- if(response.type  == "uncertainty" && uncert.cov)           uncertainty$col else if(isFALSE(scat.null)) scatter.pars$col else colors[claSS]
  } else stripplot.pars$col       <- if(response.type  == "uncertainty" && uncert.cov) uncertainty$col else stripplot.pars$strip.col
  if(is.null(stripplot.pars$jitter)) {
  stripplot.pars$jitter <- TRUE
  }
  noise.cols <- scatter.pars$col
  noise.cols <- unique(noise.cols)[stats::na.omit(match(if(noise) c(Gseq, G + 1) else Gseq, unique(claSS)))]
  if(is.null(boxplot.pars$box.pch))  {
   boxplot.pars$box.pch <- "|"
  }
  if(is.null(boxplot.pars$box.col))  {
   boxplot.pars$box.col <- "black"
  }
  boxplot.pars$varwidth <- !is.null(boxplot.pars$varwidth) && isTRUE(boxplot.pars$varwidth)
  boxplot.pars$notch    <- !is.null(boxplot.pars$notch)    && isTRUE(boxplot.pars$notch)
  if(is.null(boxplot.pars$notch.frac))    {
    boxplot.pars$notch.frac         <- 0.5
  }
  if(is.null(boxplot.pars$box.fill))      {
   boxplot.pars$box.fill            <- noise.cols
  } else {
   if(length(boxplot.pars$box.fill) == 1) {
    boxplot.pars$box.fill           <- rep(boxplot.pars$box.fill, L)
   }
   if(!(length(boxplot.pars$box.fill) %in% c(1L, L))) stop("'boxplot.pars$box.fill' must be a scalar or a vector with length given by the number of components (incl. the noise component, if any)", call.=FALSE)
  }
  bar.col               <- FALSE
  if(is.null(barcode.pars$bar.col))  {
    barcode.pars$col    <- if(isFALSE(scat.null)) rev(noise.cols) else colors[rev(seq_along(unique(claSS)))]
  } else {
    if(length(barcode.pars$bar.col) == 1) {
      bar.col           <- TRUE
      barcode.pars$col              <- rep(barcode.pars$bar.col, L)
    } else barcode.pars$col         <- barcode.pars$bar.col
    if(!(length(barcode.pars$bar.col) %in% c(1L, L))) stop("'barcode.pars$bar.col' must be a scalar or a vector with length given by the number of components (incl. the noise component, if any)", call.=FALSE)
  }
  if(is.null(barcode.pars$nint))     {
    barcode.pars$nint   <- 0
  }
  if(is.null(barcode.pars$ptsize))   {
    barcode.pars$ptsize <- grid::unit(0.25, "char")
  }
  if(is.null(barcode.pars$ptpch))    {
    barcode.pars$ptpch  <- 1
  }
  if(is.null(barcode.pars$use.points)) {
    barcode.pars$use.points <- FALSE
  }
  if(is.null(mosaic.pars$gp_labels)) {
    mosaic.pars$gp_labels   <- grid::gpar(fontsize=9)
  }
  if(is.null(mosaic.pars$gp_args))   {
    mosaic.pars$gp_args <- list()
  }
  if(is.null(mosaic.pars$shade))     {
    mosaic.pars$shade   <- NULL
  }
  if(is.null(mosaic.pars$mcol))      {
    mosaic.pars$mcol    <- noise.cols
  } else {
    if(length(mosaic.pars$mcol)     == 1) {
      mosaic.pars$mcol  <- rep(mosaic.pars$mcol, L)
    }
    if(!(length(mosaic.pars$mcol) %in% c(1L, L)))     stop("'mosaicpars$mcol' must be a scalar or a vector with length given by the number of components (incl. the noise component, if any)", call.=FALSE)
  }
  mosaic.pars$gp$col    <- if(is.null(mosaic.pars$gp$col)) "black" else mosaic.pars$gp$col
  mosaic.pars$mfill     <- is.null(mosaic.pars$mfill) || isTRUE(mosaic.pars$mfill)
  grid::grid.newpage()
  vp.main    <- grid::viewport(x=outer.margins$bottom, y=outer.margins$left,
                               width=grid::unit(1,  "npc") - outer.margins$right - outer.margins$left,
                               height=grid::unit(1, "npc") - outer.margins$top   - outer.margins$bottom,
                               just=c("left", "bottom"), name="main", clip="off")
  grid::pushViewport(vp.main)
  for(i   in Nseq) {
    for(j in Nseq) {
      MAPpan <- (i == 1L || j == 1L) && names(x)[1L] == "MAP"
      bg     <- if(isTRUE(MAPpan)) bg.col[1L]     else if(i <= dcol && j <= dcol) bg.col[2L]     else if(i > dcol && j > dcol) bg.col[5L]     else if(j > dcol && i <= dcol) bg.col[3L]     else bg.col[4L]
      border <- if(isTRUE(MAPpan)) border.col[1L] else if(i <= dcol && j <= dcol) border.col[2L] else if(i > dcol && j > dcol) border.col[5L] else if(j > dcol && i <= dcol) border.col[3L] else border.col[4L]
      labelj <- ifelse(diagonal, j, N - j + 1L)
      x[is.infinite(x[,i]), i] <- NA
      x[is.infinite(x[,j]), j] <- NA
      vp     <- grid::viewport(x=(labelj - 1L)/N, y=1 - i/N, width=1/N, height=1/N, just=c("left", "bottom"), name=as.character(i * N + j))
      grid::pushViewport(vp)
      vp.in  <- grid::viewport(x=0.5, y=0.5, width=1 - gap, height=1 - gap, just=c("center", "center"), name=paste("IN", as.character(i * N + j)))
      grid::pushViewport(vp.in)
      xpos   <- NULL
      if(i == 1 && outer.labels$top[j])    {
        xpos <- FALSE
      }
      if(i == N && outer.labels$bottom[j]) {
        xpos <- TRUE
      }
      ypos   <- NULL
      if(j == N && outer.labels$right[i])  {
        ypos <- FALSE
      }
      if(j == 1 && outer.labels$left[i])   {
        ypos <- TRUE
      }
      if(!is.null(ypos) && !diagonal)      {
        ypos <- !ypos
      }
      if(i == j) {
        diag.pars$show.dens <- show.D && i <= dcol && !MAPpan
        .diag_panel(x=x[,i], varname=names(x)[i], diag.pars=diag.pars, hist.col=if(i == 1 && names(x)[i] == "MAP") list(noise.cols) else diag.pars$hist.color, axis.pars=axis.pars, xpos=xpos, ypos=ypos, buffer=buffer, index=i, i=i - isTRUE(subset$show.map), res=res, col=noise.cols, outer.rot=outer.rot)
      } else  {
        if(xor(is.factor(x[,i]), is.factor(x[,j]))) {
          if(i < j & upr.cond != "barcode") .boxplot_panel(x=x[,j], y=x[,i], type=upr.cond, axis.pars=axis.pars, xpos=xpos, ypos=ypos, buffer=buffer, stripplot.pars=stripplot.pars, boxplot.pars=boxplot.pars, outer.rot=outer.rot, bg=bg, box.fill=if(isTRUE(MAPpan)) boxplot.pars$box.fill else if(i > dcol && j > dcol) diag.pars$hist.color[j] else "white", border=border)
          if(i > j & low.cond != "barcode") .boxplot_panel(x=x[,j], y=x[,i], type=low.cond, axis.pars=axis.pars, xpos=xpos, ypos=ypos, buffer=buffer, stripplot.pars=stripplot.pars, boxplot.pars=boxplot.pars, outer.rot=outer.rot, bg=bg, box.fill=if(isTRUE(MAPpan)) boxplot.pars$box.fill else if(i > dcol && j > dcol) diag.pars$hist.color[j] else "white", border=border)
          if(i < j & upr.cond == "barcode") {
            grid::pushViewport(grid::viewport(gp=grid::gpar(fill=bg)))
            if(is.factor(x[,i])) {
              .bar_code(x=split(x[,j], x[,i])[length(levels(x[,i])):1L], horizontal=TRUE, xlim=NULL, labelloc=ypos, axisloc=xpos, labelouter=TRUE, 
                        newpage=FALSE, fontsize=axis.pars$fontsize, buffer=buffer, nint=barcode.pars$nint, ptsize=barcode.pars$ptsize, 
                        ptpch=barcode.pars$ptpch, bcspace=barcode.pars$bcspace, use.points=barcode.pars$use.points, outerbox=border, 
                        col=if(isTRUE(MAPpan) || bar.col) barcode.pars$col else rep(diag.pars$hist.color[ifelse(j > dcol && i > dcol, j, i)], nlevels(x[,i])))
            } else {
              if(!is.null(ypos)) ypos <- !ypos
              .bar_code(x=split(x[,i], x[,j])[length(levels(x[,j])):1L], horizontal=FALSE, xlim=NULL, labelloc=xpos, axisloc=ypos, labelouter=TRUE, 
                        newpage=FALSE, fontsize=axis.pars$fontsize, buffer=buffer, nint=barcode.pars$nint, ptsize=barcode.pars$ptsize, 
                        ptpch=barcode.pars$ptpch, bcspace=barcode.pars$bcspace, use.points=barcode.pars$use.points, outerbox=border, 
                        col=if(isTRUE(MAPpan) || bar.col) barcode.pars$col else rep(diag.pars$hist.color[ifelse(j > dcol && i > dcol, j, i)], nlevels(x[,j])))
            }
            grid::popViewport()
          }
          if(i > j & low.cond == "barcode") {
            grid::pushViewport(grid::viewport(gp=grid::gpar(fill=bg)))
            if(is.factor(x[,i])) {
              .bar_code(x=split(x[,j], x[,i])[length(levels(x[,i])):1L], horizontal=TRUE, xlim=NULL, labelloc=ypos, axisloc=xpos, labelouter=TRUE, 
                        newpage=FALSE, fontsize=axis.pars$fontsize, buffer=buffer, nint=barcode.pars$nint, ptsize=barcode.pars$ptsize, 
                        ptpch=barcode.pars$ptpch, bcspace=barcode.pars$bcspace, use.points=barcode.pars$use.points, outerbox=border, 
                        col=if(isTRUE(MAPpan) || bar.col) barcode.pars$col else rep(diag.pars$hist.color[ifelse(j > dcol && i > dcol, i, j)], nlevels(x[,i])))
            } else {
              if(!is.null(ypos)) ypos <- !ypos
              .bar_code(x=split(x[,i], x[,j])[length(levels(x[,j])):1L], horizontal=FALSE, xlim=NULL, labelloc=xpos, axisloc=ypos, labelouter=TRUE, 
                        newpage=FALSE, fontsize=axis.pars$fontsize, buffer=buffer, nint=barcode.pars$nint, ptsize=barcode.pars$ptsize, 
                        ptpch=barcode.pars$ptpch, bcspace=barcode.pars$bcspace, use.points=barcode.pars$use.points, outerbox=border, 
                        col=if(isTRUE(MAPpan) || bar.col) barcode.pars$col else rep(diag.pars$hist.color[ifelse(j > dcol && i > dcol, i, j)], nlevels(x[,j])))
            }
            grid::popViewport()
          }
        }
        if(!any(is.factor(x[,i]), is.factor(x[,j])))  {
          if(response.type == "density" && all(j <= dcol, i <= dcol)) {
            .density_panel(cbind(x[,j], x[,i]), dimens=c(subset$data.ind[j - subset$show.map], subset$data.ind[i - subset$show.map]), res, density.pars, axis.pars, xpos, ypos, buffer, outer.rot, bg, border)
          } else {
            scatter.pars$pch   <- if(response.type == "uncertainty" && (all(j <= dcol, i <= dcol) || isTRUE(uncert.cov))) scatter.pars$uncert.pch else scatter.pars$scat.pch
            .scatter_panel(x=x[,j], y=x[,i], type=ifelse(j > dcol   && i <= dcol, ifelse(upr.gate || (j %in% c(expx, both)), upr.exp, "points"), ifelse(j <= dcol && i <= dcol, ifelse(drawEllipses, "ellipses", "points"), ifelse(j <= dcol && (low.gate || (i %in% c(expx, both))), low.exp, "points"))),
                           scatter.pars=scatter.pars, axis.pars=axis.pars, xpos=xpos, ypos=ypos, buffer=buffer, z=z, G=G, res=res, dimens=c(subset$data.ind[j - subset$show.map], subset$data.ind[i - subset$show.map]), outer.rot=outer.rot, bg=bg, mvn.type=addEllipses, border=border,
                           uncertainty=if(response.type == "uncertainty" && (uncert.cov  || (i <= dcol && j <= dcol))) uncertainty else NA, mvn.col=if(colEllipses) scatter.pars$eci.col)
          }
        }
        if(all(is.factor(x[,i]),  is.factor(x[,j]))) {
          .mosaic_panel(x=x[,j], y=x[,i], mosaic.pars=mosaic.pars, mosaic.inner=if(j == 1) rep(mosaic.pars$mcol, each=nlevels(x[,i]))   else if(i == 1) rep(mosaic.pars$mcol, nlevels(x[,j])) else hist.col[4L], axis.pars=axis.pars, xpos=xpos, ypos=ypos, outer.rot=outer.rot, bg=ifelse(i == 1 || j == 1, ifelse(isFALSE(mosaic.pars$mfill), bg, mosaic.pars$gp$col), bg), mfill=(i == 1 || j == 1) && isTRUE(mosaic.pars$mfill))
        }
      }
      grid::popViewport(1)
      grid::upViewport()
    }
  }
  grid::popViewport()
    invisible()
}

#' Plot MoEClust Gating Network
#'
#' Plots the gating network for fitted MoEClust models, i.e. the observation index against the mixing proportions for that observation, coloured by cluster.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param x.axis Optional argument for the x-axis against which the mixing proportions are plotted. Defaults to \code{1:res$n} if missing. Supplying \code{x.axis} changes the defaults for the \code{type} and \code{xlab} arguments. Users are advised to only use quantities related to the gating network of the fitted model here. Furthermore, use of the \code{x.axis} argument is not recommended for models with more than one gating network covariate.
#' @param type,pch,xlab,ylab,ylim,col These graphical parameters retain their definitions from \code{\link[graphics]{matplot}}. \code{col} defaults to the settings in \code{\link[mclust]{mclust.options}}. Note that the default value of \code{type} changes depending on whether \code{x.axis} is supplied and whether the gating network contains multiple covariates &/or categorical covariates.
#' @param ... Catches unused arguments, or additional arguments to be passed to \code{\link[graphics]{matplot}}.
#'
#' @return A plot of the gating network of the fitted MoEClust model. The parameters of the gating network can also be returned invisibly.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @note \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_plotGate}} which accepts the default arguments, and also produces other types of plots.
#'
#' By default, the noise component (if any) will be coloured \code{"darkgrey"}.
#' @seealso \code{\link{MoE_clust}}, \code{\link{plot.MoEClust}}, \code{\link[graphics]{matplot}}
#' @importFrom mclust "mclust.options"
#' @keywords plotting
#' @export
#' @usage
#' MoE_plotGate(res,
#'              x.axis = NULL,
#'              type = "b",
#'              pch = 1,
#'              xlab = "Observation",
#'              ylab = expression(widehat(tau)[g]),
#'              ylim = c(0, 1),
#'              col = NULL,
#'              ...)
#' @examples
#' data(ais)
#' res   <- MoE_clust(ais[,3:7], gating= ~ BMI, G=3, modelNames="EEE",
#'                    network.data=ais, noise.gate=FALSE, tau0=0.1)
#'
#' # Plot against the observation index and examine the gating network coefficients
#' (gate <- MoE_plotGate(res))
#'
#' # Plot against BMI
#' MoE_plotGate(res, x.axis=ais$BMI, xlab="BMI")
#' 
#' # Plot against a categorical covariate
#' res2  <- MoE_clust(ais[,3:7], gating= ~ sex, G=3, modelNames="EVE", network.data=ais)
#' MoE_plotGate(res2, x.axis=ais$sex, xlab="sex")
MoE_plotGate  <- function(res, x.axis = NULL, type = "b", pch = 1, xlab = "Observation", ylab = expression(widehat(tau)[g]), ylim = c(0, 1), col = NULL, ...) {
    UseMethod("MoE_plotGate")
}

#' @method MoE_plotGate MoEClust
#' @export
MoE_plotGate.MoEClust   <- function(res, x.axis = NULL, type = "b", pch = 1, xlab = "Observation", ylab = expression(widehat(tau)[g]), ylim = c(0, 1), col = NULL, ...) {
  res        <- if(inherits(res, "MoECompare")) res$optimal                else res
  oldpar     <- suppressWarnings(graphics::par(no.readonly=TRUE))
  oldpar$new <- FALSE
  on.exit(suppressWarnings(graphics::par(oldpar)))
  suppressWarnings(graphics::par(pty="m"))
  N          <- res$n
  G          <- res$G
  if(G == 1)                                          message("No clustering has taken place!\n")
  Tau        <- .mat_byrow(res$parameters$pro, nrow=N, ncol=ncol(res$z))
  vars       <- all.vars(stats::as.formula(attr(res$gating, "Formula")))
  ncovs      <- length(vars)  > 1
  axmiss     <- is.null(x.axis)
  tmiss      <- missing(type)
  xlab       <- if(!axmiss   && missing(xlab)) deparse(substitute(x.axis)) else xlab
  x.axis     <- if(axmiss)      seq_len(N) else x.axis
  type       <- if(tmiss     && 
                   !axmiss   && ncovs) "p" else type 
  if(length(x.axis) != N)                             stop("'x.axis' must be of length N", call.=FALSE)
  if(x.fac   <- is.factor(x.axis)) {
    xlev     <- levels(x.axis)
    x.axis   <- as.integer(x.axis)
    xaxt     <- "n"
  } else      {
    type     <- ifelse(tmiss && any(vars %in% names(res$gating$xlevels)), "p", type)
    xaxt     <- "s"
  }
  if(isFALSE(axmiss))   {
    if(isTRUE(ncovs))                                 warning("Function may produce undesirable plot when 'x.axis' is supplied for a model with multiple gating network covariates\n", call.=FALSE, immediate.=TRUE)
    o.axis   <- order(x.axis, decreasing=FALSE)
    x.axis   <- x.axis[o.axis]
    Tau      <- Tau[o.axis,,        drop=FALSE]
  }
  cX         <- missing(col)
  nX         <- is.na(res$hypvol)
  col        <- if(cX) mclust.options("classPlotColors")[seq_len(G)]        else col
  col        <- if(cX  && !nX) c(rep(col, length.out=G), "darkgrey")        else col
  if(missing(ylab))     {
    graphics::matplot(x=x.axis, y=Tau, type=type, pch=pch, xlab="",   ylab="",   xaxt=xaxt, ylim=ylim, col=col, ...)
    dots     <- list(...)
    dots     <- dots[unique(names(dots))]
    dots["cex.axis"]   <- dots["cex.lab"]
    dots["cex.lab"]    <- NULL
    graphics::mtext(side=1, xlab, las=1, line=2,   cex=dots[["cex.axis"]])
    graphics::mtext(side=2, ylab, las=2, line=2.5, cex=dots[["cex.axis"]])
  } else      {
    graphics::matplot(x=x.axis, y=Tau, type=type, pch=pch, xlab=xlab, ylab=ylab, xaxt=xaxt, ylim=ylim, col=col, ...)
  }
  if(x.fac)   {
    graphics::axis(1, at=unique(x.axis), labels=xlev)
  }
    invisible(res$gating)
}

#' Model Selection Criteria Plot for MoEClust Mixture Models
#'
#' Plots the BIC, ICL, AIC, or log-likelihood values of a fitted \code{MoEClust} object.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param criterion The criterion to be plotted. Defaults to \code{"bic"}.
#' @param ... Catches other arguments, or additional arguments to be passed to \code{\link[mclust]{plot.mclustBIC}} (or equivalent functions for the other \code{criterion} arguments). In particular, the argument \code{legendArgs} to \code{\link[mclust]{plot.mclustBIC}} can be passed.
#'
#' @importFrom mclust "plot.mclustBIC" "plot.mclustICL"
#' @return A plot of the values of the chosen \code{criterion}. The values themselves can also be returned invisibly.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @keywords plotting
#' @export
#' @note \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_plotCrit}} which accepts the default arguments, and also produces other types of plots.
#' @seealso \code{\link{MoE_clust}}, \code{\link{plot.MoEClust}}, \code{\link[mclust]{plot.mclustBIC}}
#' @usage
#' MoE_plotCrit(res,
#'              criterion = c("bic", "icl", "aic", "loglik"),
#'              ...)
#' @examples
#' \donttest{# data(ais)
#' # res   <- MoE_clust(ais[,3:7], expert= ~ sex, network.data=ais)
#' # (crit <- MoE_plotCrit(res))}
MoE_plotCrit <- function(res, criterion = c("bic", "icl", "aic", "loglik"), ...) {
    UseMethod("MoE_plotCrit")
}

#' @method MoE_plotCrit MoEClust
#' @export
MoE_plotCrit.MoEClust   <- function(res, criterion = c("bic", "icl", "aic", "loglik"), ...) {
  res        <- if(inherits(res, "MoECompare")) res$optimal else res
  oldpar     <- suppressWarnings(graphics::par(no.readonly=TRUE))
  oldpar$new <- FALSE
  on.exit(suppressWarnings(graphics::par(oldpar)))
  suppressWarnings(graphics::par(pty="m"))
  if(!missing(criterion)      &&
     (length(criterion)  > 1  ||
      !is.character(criterion)))                      stop("'criterion' must be a single character string", call.=FALSE)
  criterion  <- match.arg(criterion)
  crit       <- switch(EXPR=criterion, bic=res$BIC, icl=res$ICL, res$AIC)
  crit2      <- replace(crit, !is.finite(crit), NA)
    switch(EXPR=criterion, bic=plot.mclustBIC(crit2, ...), icl=plot.mclustICL(crit2, ...),
                     aic=plot.mclustAIC(crit2, ...), loglik=plot.mclustLoglik(crit2, ...))
    invisible(crit)
}

#' Plot the Log-Likelihood of a MoEClust Mixture Model
#'
#' Plots the log-likelihood at every iteration of the EM/CEM algorithm used to fit a MoEClust mixture model.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param type,xlab,ylab,xaxt These graphical parameters retain their usual definitions from \code{\link{plot}}.
#' @param ... Catches unused arguments, or additional arguments to be passed to \code{\link{plot}}.
#'
#' @return A plot of the log-likelihood versus the number EM iterations. A list with the vector of log-likelihood values and the final value at convergence can also be returned invisibly.
#' @note \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_plotLogLik}} which accepts the default arguments, and also produces other types of plots.
#'
#' \code{res$LOGLIK} can also be plotted, to compare maximal log-likelihood values for all fitted models.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @keywords plotting
#' @export
#' @seealso \code{\link{MoE_clust}}, \code{\link{plot.MoEClust}},
#' @usage
#' MoE_plotLogLik(res,
#'                type = "l",
#'                xlab = "Iteration",
#'                ylab = "Log-Likelihood",
#'                xaxt = "n",
#'                ...)
#' @examples
#' data(ais)
#' res <- MoE_clust(ais[,3:7], gating= ~ BMI, expert= ~ sex, tau0=0.1,
#'                  G=2, modelNames="EVE", network.data=ais)
#' (ll <- MoE_plotLogLik(res))
MoE_plotLogLik           <- function(res, type = "l", xlab = "Iteration", ylab = "Log-Likelihood", xaxt = "n", ...) {
    UseMethod("MoE_plotLogLik")
}

#' @method MoE_plotLogLik MoEClust
#' @export
MoE_plotLogLik.MoEClust  <- function(res, type = "l", xlab = "Iteration", ylab = "Log-Likelihood", xaxt = "n", ...) {
  res        <- if(inherits(res, "MoECompare")) res$optimal else res
  oldpar     <- suppressWarnings(graphics::par(no.readonly=TRUE))
  oldpar$new <- FALSE
  on.exit(suppressWarnings(graphics::par(oldpar)))
  suppressWarnings(graphics::par(pty="m"))
  xll        <- res$loglik
  if(res$G   == 1)                                    message("EM algorithm not used; no clustering has taken place!\n")
  if(all(xll != cummax(xll)))                         warning("Log-likelihoods are not strictly increasing\n", call.=FALSE)
  base::plot(xll, type = ifelse(length(xll) == 1, "p", type), xlab = xlab, ylab = ylab, xaxt = xaxt, ...)
  if(length(xaxt) == 1   && is.character(xaxt))  {
    seqll    <- seq_along(xll)
    llseq    <- pretty(seqll)
    llseq    <- if(any(llseq  != floor(llseq))) seqll else  llseq
    if(xaxt  == "n") graphics::axis(1, at = llseq, labels = llseq)
  } else                                              stop("'xaxt' must be a single character string", call.=FALSE)
    invisible(list(ll = xll, converge = xll[length(xll)]))
}

#' Plot the Similarity Matrix of a MoEClust Mixture Model
#'
#' Produces a heatmap of the similarity matrix constructed from the \code{res$z} matrix at convergence of a MoEClust mixture model.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param col A vector of colours as per \code{\link[graphics]{image}}. Will be checked for validity.
#' @param reorder A logical (defaults to \code{TRUE}) indicating whether observations should be reordered for visual clarity.
#' @param legend A logical (defaults to \code{TRUE}) indicating whether to append a colour key legend.
#' @param ... Catches unused arguments, or arguments to be passed to \code{\link[stats]{hclust}} when \code{reorder=TRUE}.
#' @note \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_Similarity}} which accepts the default arguments, and also produces other types of plots.
#' @return The similarity matrix in the form of a heatmap is plotted; the matrix itself can also be returned invisibly. The invisibly returned matrix will also be reordered if \code{reordered=TRUE}.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @keywords plotting
#' @export
#' @seealso \code{\link{MoE_clust}}, \code{\link{plot.MoEClust}},
#' @usage
#' MoE_Similarity(res,
#'                col = grDevices::heat.colors(30L, rev=TRUE), 
#'                reorder = TRUE, 
#'                legend = TRUE,
#'                ...)
#' @examples
#' data(ais)
#' mod <- MoE_clust(ais[,3:7], G=2, modelNames="EEE", gating= ~ SSF + Ht,
#'                  expert= ~ sex, network.data=ais, tau0=0.1, noise.gate=FALSE)
#' sim <- MoE_Similarity(mod)
MoE_Similarity <- function(res, col = grDevices::heat.colors(30L, rev=TRUE), reorder = TRUE, legend = TRUE, ...) {
    UseMethod("MoE_Similarity")
}

#' @method MoE_Similarity MoEClust
#' @export
MoE_Similarity.MoEClust       <- function(res, col = grDevices::heat.colors(30L, rev=TRUE), reorder = TRUE, legend = TRUE, ...) {
  sim          <- if(res$G > 1) tcrossprod(res$z) else matrix(1L, nrow=res$n, ncol=res$n)
  if(isTRUE(reorder)) {
    dots       <- list(...)
    args       <- c(list(sim=sim), dots[names(dots) %in% c("method", "members")])
    perm       <- do.call(.reorder_hc, args[unique(names(args))])
    sim        <- if((reorder <- !identical(perm, res$classification))) sim[perm,perm] else sim
  }
  if(!all(.is_cols(col)))                             stop("Invalid 'cols' colour palette supplied", call.=FALSE)
  if(isTRUE(legend))  {
    oldpar     <- suppressWarnings(graphics::par(no.readonly=TRUE))
    oldpar$new <- FALSE
    on.exit(suppressWarnings(graphics::par(oldpar)))
    graphics::par(mar=c(5.1, 4.1, 4.1, 3.1))
  }
  graphics::image(t(sim[seq(from=ncol(sim), to=1L, by=-1L),]), col=col, main="Similarity Matrix",
                  xlab=paste0("Observation 1:N", if(isTRUE(reorder)) " (Reordered)"),
                  ylab=paste0("Observation 1:N", if(isTRUE(reorder)) " (Reordered)"))
  graphics::box(lwd=2)
  if(isTRUE(legend)) suppressWarnings(.heat_legend(data=sim, col=col, cex.lab=0.8))
    invisible(if(isTRUE(reorder)) provideDimnames(sim, base=list(as.character(perm), as.character(perm))) else sim) # fix
}

#' Plot Clustering Uncertainties
#'
#' Plots the clustering uncertainty for every observation from a fitted \code{"MoEClust"} model, including models with a noise component.
#' @param res An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param type The type of plot to be produced (defaults to \code{"barplot"}). The \code{"profile"} option instead displays uncertainties in increasing/decreasing order of magnitude (see \code{decreasing}).
#' @param truth An optional argument giving the true classification of the data. When \code{truth} is supplied and \code{type="barplot"}, misclassified observations are highlighted in a different colour, otherwise observations with uncertainty greater than \code{1/res$G} are given in a different colour. When \code{truth} is supplied and \code{type="profile"}, the uncertainty of misclassified observations is marked by vertical lines on the plot.
#' @param decreasing Logical indicating whether uncertainties should be ordered in decreasing order (defaults to \code{FALSE}). Only relevant when \code{type="profile"}.
#' @param ... Catches unused arguments.
#'
#' @details The y-axis of this plot runs from \code{0} to \code{1 - 1/res$G}, with a horizontal line also drawn at \code{1/res$G}. When \code{type="barplot"}, uncertainties greater than this value are given a different colour when \code{truth} is not supplied, otherwise misclassified observations are given a different colour. Note, however, that \eqn{G^{(0)}}{G^(0)} = \code{res$G + 1} is used in place of \code{res$G} for models with a noise component.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @note \code{\link{plot.MoEClust}} is a wrapper to \code{\link{MoE_Uncertainty}} which accepts the default arguments, and also produces other types of plots.
#'
#' An alternative means of visualising clustering uncertainties (at least for multivariate data) is provided by the functions \code{MoE_gpairs} and \code{plot.MoEClust}, specifically when their argument \code{response.type} is given as \code{"uncertainty"}.
#' @return A plot showing the clustering uncertainty of each observation (sorted in increasing/decreasing order when \code{type="profile"}). The (unsorted) vector of uncertainties can also be returned invisibly. When \code{truth} is supplied, the indices of the misclassified observations are also invisibly returned.
#' @seealso \code{\link{MoE_clust}}, \code{\link{MoE_gpairs}}, \code{\link{plot.MoEClust}}
#' @export
#'
#' @importFrom mclust "classError" "mclust.options"
#' @keywords plotting
#' @usage
#' MoE_Uncertainty(res,
#'                 type = c("barplot", "profile"),
#'                 truth = NULL,
#'                 decreasing = FALSE,
#'                 ...)
#' @examples
#' data(ais)
#' res <- MoE_clust(ais[,3:7], gating= ~ sex, G=3, modelNames="EEE", network.data=ais)
#'
#' # Produce an uncertainty barplot
#' MoE_Uncertainty(res)
#'
#' # Produce an uncertainty profile plot
#' MoE_Uncertainty(res, type="profile")
#'
#' # Let's assume the true clusters correspond to sex
#' (ub <- MoE_Uncertainty(res, truth=ais$sex))
#' (up <- MoE_Uncertainty(res, type="profile", truth=ais$sex))
MoE_Uncertainty          <- function(res, type = c("barplot", "profile"), truth = NULL, decreasing = FALSE, ...) {
    UseMethod("MoE_Uncertainty")
}

#' @method MoE_Uncertainty MoEClust
#' @importFrom mclust "classError" "mclust.options"
#' @export
MoE_Uncertainty.MoEClust <- function(res, type = c("barplot", "profile"), truth = NULL, decreasing = FALSE, ...) {
  res        <- if(inherits(res, "MoECompare")) res$optimal else res
  noise      <- !is.na(res$hypvol)
  oldpar     <- suppressWarnings(graphics::par(no.readonly=TRUE))
  oldpar$new <- FALSE
  on.exit(suppressWarnings(graphics::par(oldpar)))
  suppressWarnings(graphics::par(pty="m", mar=replace(graphics::par()$mar, c(2, 4), c(ifelse(noise, 5.1, 4.1), 3.1))))
  if(!missing(type)      && (length(type) > 1     ||
     !is.character(type)))                            stop("'type' must be a single character string",        call.=FALSE)
  if(length(decreasing)   > 1                     ||
     !is.logical(decreasing))                         stop("'decreasing' must be a single logical indicator", call.=FALSE)
  type       <- match.arg(type)
  n.obs      <- res$n
  ucert      <- res$uncertainty
  if(tmiss   <- !missing(truth)) {
    if(all(!is.factor(truth), !is.numeric(truth)) ||
       n.obs != length(truth))                        stop(paste0("'truth' must be a numeric or factor vector of length N=",  n.obs), call.=FALSE)
    mC       <- classError(classification=res$classification, class=as.numeric(as.factor(truth)))$misclassified
  }
  G          <- res$G + noise
  if(G == 1)                                          message("No clustering has taken place!\n")
  if(attr(res, "Algo") == "CEM")                      message("Model was fitted by CEM: no uncertainties to plot\n")
  oneG       <- 1/G
  min1G      <- 1   - oneG
  yx         <- unique(c(0, pretty(c(0, min1G))))
  YX         <- which.min(abs(yx - min1G))
  yx[YX]     <- min1G
  yx         <- abs(yx[yx < 1])
  cm         <- mclust.options("classPlotColors")
  if(type    == "barplot") {
    cu       <- if(tmiss) replace(rep(cm[1L], n.obs), mC, cm[2L]) else cm[seq_len(2L)][(ucert >= oneG) + 1L]
    cu[ucert == 0] <- NA
    base::plot(ucert, type="h", ylim=range(yx), col=cu, yaxt="n", ylab="", xlab="Observations", lend=1)
    graphics::lines(x=c(0, n.obs), y=c(oneG, oneG), lty=2, col=cm[3L])
    graphics::axis(2, at=yx, labels=replace(yx, YX, ifelse(noise, expression(1 - frac(1, widehat(G^{'(0)'}))), expression(1 - frac(1, hat(G))))), las=2, xpd=TRUE)
    graphics::axis(2, at=oneG, labels=ifelse(noise, expression(frac(1, widehat(G^{'(0)'}))), expression(frac(1, hat(G)))), las=2, xpd=TRUE, side=4)
  } else      {
    ord      <- order(ucert, decreasing=decreasing)
    ucord    <- ucert[ord]
    if(tmiss)  mcO <- which(ord %in% mC)
    base::plot(ucord, type="n", ylim=c(-max(ucert)/32, max(yx)), ylab="", xaxt="n", yaxt="n", xlab=paste0("Observations in order of ", ifelse(decreasing, "decreasing", "increasing"), " uncertainty"))
    graphics::lines(x=c(0, n.obs), y=c(0, 0), lty=3)
    graphics::lines(ucord)
    graphics::points(ucord, pch=15, cex=if(tmiss) replace(rep(0.5, n.obs), mcO, 0.75) else 0.5, col=if(tmiss) replace(rep(1, n.obs), mcO, cm[2L]) else 1)
    graphics::lines(x=c(0, n.obs), y=c(oneG, oneG), lty=2, col=cm[3L])
    graphics::axis(2, at=yx, labels=replace(yx, YX, ifelse(noise, expression(1 - frac(1, widehat(G^{'(0)'}))), expression(1 - frac(1, hat(G))))), las=2, xpd=TRUE)
    graphics::axis(2, at=oneG, labels=ifelse(noise, expression(frac(1, widehat(G^{'(0)'}))), expression(frac(1, hat(G)))), las=2, xpd=TRUE, side=4)
    if(tmiss) {
      Nseq   <- (seq_len(n.obs))
      for(i in mC)  {
        x    <- Nseq[ord == i]
        graphics::lines(c(x, x), c(-max(ucert)/32, ucert[i]), lty=1, lend=1)
      }
    }
  }
  graphics::mtext(paste0("Uncertainty", ifelse(noise, "\n(with a noise component)", "")), side=2, line=3)
  if(tmiss)   {
    U        <- list(Uncertainties = ucert, Misclassified = mC)
    class(U) <- "listof"
  } else U   <- ucert
    invisible(U)
}

#' Plot MoEClust Results
#'
#' Plot results for fitted MoE_clust mixture models with gating &/or expert network covariates: generalised pairs plots, model selection criteria, the log-likelihood vs. the EM iterations, and the gating network are all currently visualisable.
#' @param x An object of class \code{"MoEClust"} generated by \code{\link{MoE_clust}}, or an object of class \code{"MoECompare"} generated by \code{\link{MoE_compare}}. Models with a noise component are facilitated here too.
#' @param what The type of graph requested:
#' \describe{
#' \item{\code{gpairs}}{A generalised pairs plot. To further customise this plot, arguments to \code{\link{MoE_gpairs}} can be supplied.}
#' \item{\code{gating}}{The gating network. To further customise this plot, arguments to \code{\link{MoE_plotGate}} and \code{\link[graphics]{matplot}} can be supplied.}
#' \item{\code{criterion}}{The model selection criteria. To further customise this plot, arguments to \code{\link{MoE_plotCrit}} and \code{\link[mclust]{plot.mclustBIC}} can be supplied.}
#' \item{\code{loglik}}{The log-likelihood vs. the iterations of the EM algorithm. To further customise this plot, arguments to \code{\link{MoE_plotLogLik}} and \code{\link{plot}} can be supplied.}
#' \item{\code{similarity}}{The similarity matrix constructed from \code{x$z} at convergence, in the form of a heatmap. To further customise this plot, arguments to \code{\link{MoE_Similarity}} can be supplied.}
#' \item{\code{uncertainty}}{The clustering uncertainty for every observation. To further customise this plot, arguments to \code{\link{MoE_Uncertainty}} can be supplied.}
#' }
#' By default, all of the above graphs are produced.
#' @param ... Optional arguments to be passed to \code{\link{MoE_gpairs}}, \code{\link{MoE_plotGate}}, \code{\link{MoE_plotCrit}}, \code{\link{MoE_plotLogLik}}, \code{\link{MoE_Similarity}}, \code{\link{MoE_Uncertainty}}, \code{\link[graphics]{matplot}}, \code{\link[mclust]{plot.mclustBIC}} and \code{\link{plot}}. In particular, the argument \code{legendArgs} to \code{\link[mclust]{plot.mclustBIC}} can be passed to \code{\link{MoE_plotCrit}}.
#'
#' @details For more flexibility in plotting, use \code{\link{MoE_gpairs}}, \code{\link{MoE_plotGate}}, \code{\link{MoE_plotCrit}}, \code{\link{MoE_plotLogLik}}, \code{\link{MoE_Similarity}}, and \code{\link{MoE_Uncertainty}} directly.
#' @importFrom lattice "current.panel.limits" "llines" "panel.abline" "panel.bwplot" "panel.histogram" "panel.lines" "panel.points" "panel.rect" "panel.stripplot" "panel.text" "panel.violin" "trellis.grobname" "trellis.par.get" "trellis.par.set"
#' @importFrom matrixStats "rowLogSumExps"
#' @importFrom mclust "mclust.options" "plot.mclustBIC" "plot.mclustICL" "sigma2decomp"
#' @importFrom vcd "strucplot"
#' @note Caution is advised producing generalised pairs plots when the dimension of the data is large.
#'
#' Other types of plots are available by first calling \code{\link{as.Mclust}} on the fitted object, and then calling \code{\link[mclust]{plot.Mclust}} on the results. These can be especially useful for univariate data.
#' @return The visualisation according to \code{what} of the results of a fitted \code{MoEClust} model.
#' @seealso \code{\link{MoE_clust}}, \code{\link{MoE_stepwise}}, \code{\link{MoE_gpairs}}, \code{\link{MoE_plotGate}}, \code{\link{MoE_plotCrit}}, \code{\link{MoE_plotLogLik}}, \code{\link{MoE_Similarity}}, \code{\link{MoE_Uncertainty}}, \code{\link{as.Mclust}}, \code{\link[mclust]{plot.Mclust}}
#' @references Murphy, K. and Murphy, T. B. (2020). Gaussian parsimonious clustering models with covariates and a noise component. \emph{Advances in Data Analysis and Classification}, 14(2): 293-325. <\doi{10.1007/s11634-019-00373-8}>.
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @export
#' @method plot MoEClust
#' @keywords plotting main
#' @usage
#' \method{plot}{MoEClust}(x,
#'      what = c("gpairs", "gating", "criterion", "loglik", "similarity", "uncertainty"),
#'      ...)
#' @examples
#' \donttest{data(ais)
#' res <- MoE_clust(ais[,3:7], gating= ~ BMI, expert= ~ sex,
#'                  G=2, modelNames="EVE", network.data=ais)
#'
#' # Plot the gating network
#' plot(res, what="gating", x.axis=ais$BMI, xlab="BMI")
#'
#' # Plot the log-likelihood
#' plot(res, what="loglik", col="blue")
#'
#' # Plot the uncertainty profile
#' plot(res, what="uncertainty", type="profile")
#'
#' # Produce a generalised pairs plot
#' plot(res, what="gpairs")
#' 
#' # Produce a heatmap of the similarity matrix
#' plot(res, what="similarity")
#' 
#' # Modify the gpairs plot by passing arguments to MoE_gpairs()
#' plot(res, what="gpairs", response.type="density", varwidth=TRUE,
#'      data.ind=c(5,3,4,1,2), jitter=FALSE, show.counts=FALSE)}
plot.MoEClust <- function(x, what=c("gpairs", "gating", "criterion", "loglik", "similarity", "uncertainty"), ...) {
  if(!missing(what) && !all(is.character(what)))      stop("'what' must be a character string", call.=FALSE)
  what        <- match.arg(what, several.ok=TRUE)
  if(interactive()  && length(what) > 1) {
    title     <- c("MoEClust Plots\n\n<Press 0 to exit>")
    what.tmp  <- c(gpairs="Generalised Pairs Plot",
                   gating="Gating Network",
                   criterion="Model Selection Criteria",
                   loglik="Log-Likelihood",
                   similarity="Similarity Matrix",
                   uncertainty="Uncertainty")
    choice    <- utils::menu(what.tmp, graphics=FALSE, title=title)
    while(choice    != 0) {
      switch(EXPR=what[choice], 
             gpairs=MoE_gpairs(x, ...), 
             gating=MoE_plotGate(x, ...),
             criterion=MoE_plotCrit(x, ...), 
             loglik=MoE_plotLogLik(x, ...), 
             similarity=MoE_Similarity(x, ...),
             uncertainty=MoE_Uncertainty(x, ...))
      choice  <- utils::menu(what.tmp, graphics=FALSE, title=title)
    }
  } else if(length(what)  > 1) {                      stop("'what' must be a single character string for non-interactive sessions", call.=FALSE)
  } else {
    switch(EXPR=what, 
           gpairs=MoE_gpairs(x, ...), 
           gating=MoE_plotGate(x, ...),
           criterion=MoE_plotCrit(x, ...), 
           loglik=MoE_plotLogLik(x, ...), 
           similarity=MoE_Similarity(x, ...), 
           uncertainty=MoE_Uncertainty(x, ...))
  }
    invisible()
}

.grid_1   <- function(n, range = c(0, 1), edge = TRUE) {
  if(any(n < 0 | round(n) != n))                      stop("'n' must be nonpositive and integer", call.=FALSE)
  G       <- numeric(n)
  if(edge) {
    G     <- seq(from=min(range), to=max(range), by=abs(diff(range))/(n - 1L))
  } else   {
    lj    <- abs(diff(range))
    incr  <- lj/(2L * n)
    G     <- seq(from=min(range) + incr, to=max(range) - incr, by=2 * incr)
  }
    G
}

.grid_2   <- function(x, y) {
  lx      <- length(x)
  ly      <- length(y)
  xy      <- matrix(0L, nrow=lx * ly, ncol=2)
  l       <- 0L
  for(j   in seq_len(ly))   {
    for(i in seq_len(lx))   {
      l   <- l + 1L
      xy[l,]  <- c(x[i], y[j])
    }
  }
    xy
}

#' @importFrom lattice "panel.lines" "panel.points"
.mvn2D_panel <- function(mu, sigma, k = 15L, alone = FALSE, col = rep("grey30", 3), pch = 8, lty = c(1, 2), lwd = c(2, 1)) {
  panel.points(mu[1L], mu[2L], col=col[3L], pch=pch)
  ev         <- eigen(sigma, symmetric = TRUE)
  s          <- tryCatch(sqrt(rev(sort(ev$values))), warning=function(w) stop("Negative eigenvalues"))
  V          <- ev$vectors[,rev(order(ev$values))]
  theta      <- (0L:k)  * (pi/(2L * k))
  x          <- s[1L]   * cos(theta)
  y          <- s[2L]   * sin(theta)
  xy         <- sweep(tcrossprod(cbind(c(x, -x, -x, x), c(y, y, -y, -y)), V), MARGIN=2L, STATS=mu, FUN="+", check.margin=FALSE)
  l          <- length(x)
  i          <- seq_len(l)
  for(k in seq_len(4L)) {
    panel.lines(xy[i,], col=col[1L], lty=lty[1L], lwd=lwd[1L])
    i        <- i + l
  }
  x          <- s[1L]
  y          <- s[2L]
  xy         <- sweep(tcrossprod(cbind(c(x, -x, 0, 0), c(0, 0, y, -y)),   V), MARGIN=2L, STATS=mu, FUN="+", check.margin=FALSE)
  panel.lines(xy[1L:2L,], col=col[2L], lty=lty[2L], lwd=lwd[2L])
  panel.lines(xy[3L:4L,], col=col[2L], lty=lty[2L], lwd=lwd[2L])
    invisible()
}

.draw_axis <- function(x, y, axis.pars, xpos, ypos, cat.labels = NULL, horiz = NULL, xlim = NULL, ylim = NULL, outer.rot) {
  x.fac    <- is.factor(x)
  y.fac    <- is.factor(y)
  x        <- if(x.fac) x         else as.numeric(x)
  y        <- if(y.fac) y         else as.numeric(y)
  if(is.null(xlim)) {
    px     <- if(x.fac) levels(x) else pretty(x,    axis.pars$n.ticks)
    px     <- if(x.fac) px        else px[px > min(x, na.rm=TRUE)    & px < max(x, na.rm=TRUE)]
  } else {
    px     <- if(x.fac) levels(x) else pretty(xlim, axis.pars$n.ticks)
    px     <- if(x.fac) px        else px[px > min(xlim, na.rm=TRUE) & px < max(xlim, na.rm=TRUE)]
  }
  if(is.null(ylim)) {
    py     <- if(y.fac) levels(y) else pretty(y,    axis.pars$n.ticks)
    py     <- if(y.fac) py        else py[py > min(y, na.rm=TRUE)    & py < max(y, na.rm=TRUE)]
  } else {
    py     <- if(y.fac) levels(y) else pretty(ylim, axis.pars$n.ticks)
    py     <- if(y.fac) py        else py[py > min(ylim, na.rm=TRUE) & py < max(ylim, na.rm=TRUE)]
  }
  k        <- length(cat.labels)
  if(!is.null(xpos)) {
    if(!is.null(cat.labels) && isFALSE(horiz)) {
      grid::grid.text(cat.labels, x=grid::unit(seq_len(k), "native"), y=grid::unit(rep(1 * (1 - xpos), k), "npc") + grid::unit(rep(-1 * xpos + 1 * (1 - xpos), k), "lines"), rot=outer.rot[1L], gp=grid::gpar(fontsize=axis.pars$fontsize))
    } else  grid::grid.xaxis(at=if(x.fac) seq_along(px)      else px, gp=grid::gpar(fontsize=axis.pars$fontsize), main=xpos, label=px)
  }
  if(!is.null(ypos)) {
    if(!is.null(cat.labels) && isTRUE(horiz))  {
      grid::grid.text(cat.labels, y=grid::unit(seq_len(k), "native"), x=grid::unit(rep(1 * (1 - ypos), k), "npc") + grid::unit(rep(-1 * ypos + 1 * (1 - ypos), k), "lines"), rot=outer.rot[2L], gp=grid::gpar(fontsize=axis.pars$fontsize))
    } else  grid::grid.yaxis(at=if(y.fac) rev(seq_along(py)) else py, gp=grid::gpar(fontsize=axis.pars$fontsize), main=ypos, label=py)
  }
}

#' @importFrom lattice "panel.bwplot" "panel.stripplot" "trellis.par.get" "trellis.par.set"
.boxplot_panel <- function(x, y, type, axis.pars, xpos, ypos, buffer, stripplot.pars, boxplot.pars, outer.rot, bg, box.fill, border) {
  xlim         <- NULL
  ylim         <- NULL
  old.color    <- trellis.par.get("box.rectangle")$col
  trellis.par.set(name="box.rectangle", value=list(col="black"))
  trellis.par.set(name="box.umbrella",  value=list(col="black"))
  trellis.par.set(name="box.dot",       value=list(col="black"))
  trellis.par.set(name="plot.symbol",   value=list(col="black"))
  if(is.factor(x)) {
    cat.labels <- levels(x)
    k          <- length(levels(x))
    cat.var    <- as.numeric(x)
    cont.var   <- y
    horiz      <- FALSE
  } else {
    cat.labels <- levels(y)
    k          <- length(levels(y))
    cat.labels <- cat.labels[k:1L]
    cat.var    <- k + 1L - as.numeric(y)
    cont.var   <- x
    horiz      <- TRUE
  }
  grid::grid.rect(gp=grid::gpar(fill=bg, col=border))
  cont.range   <- range(cont.var, na.rm=TRUE)
  if(horiz) {
    xlim       <- cont.range + c(-buffer * (diff(cont.range)), buffer * (diff(cont.range)))
    grid::pushViewport(grid::viewport(xscale=xlim, yscale=c(0.5, max(cat.var, na.rm=TRUE) + 0.5)))
    if(is.null(ypos)) cat.labels <- NULL
    .draw_axis(x=cont.var, y=cat.var, axis.pars=axis.pars, xpos=xpos, ypos=ypos, cat.labels=cat.labels, horiz=horiz, xlim=xlim, ylim=ylim, outer.rot=outer.rot)
    grid::popViewport(1)
    grid::pushViewport(grid::viewport(xscale=xlim, yscale=c(0.5, max(cat.var, na.rm=TRUE) + 0.5), clip=TRUE))
    if(type == "boxplot")   panel.bwplot(cont.var, cat.var, horizontal=horiz, col=boxplot.pars$box.col, pch=boxplot.pars$box.pch, gp=grid::gpar(box.umbrella=list(col=boxplot.pars$box.col)), 
                                         varwidth=boxplot.pars$varwidth, notch=boxplot.pars$notch, notch.frac=boxplot.pars$notch.frac, fill=rev(box.fill))
    if(type == "stripplot") panel.stripplot(cont.var, cat.var, horizontal=horiz, jitter.data=stripplot.pars$jitter, col=stripplot.pars$col, cex=stripplot.pars$size, pch=stripplot.pars$pch)
    if(type == "violin")    .violin_panel(cont.var, cat.var, horizontal=horiz, varwidth=boxplot.pars$varwidth, border=boxplot.pars$box.col, col=rev(box.fill))
  } else    {
    ylim       <- cont.range + c(-buffer * (diff(cont.range)), buffer * (diff(cont.range)))
    grid::pushViewport(grid::viewport(yscale=ylim, xscale=c(0.5, max(cat.var, na.rm=TRUE) + 0.5)))
    if(is.null(xpos)) cat.labels <- NULL
    .draw_axis(x=cat.var, y=cont.var, axis.pars=axis.pars, xpos=xpos, ypos=ypos, cat.labels=cat.labels, horiz=horiz, xlim=xlim, ylim=ylim, outer.rot=outer.rot)
    grid::popViewport(1)
    grid::pushViewport(grid::viewport(yscale=ylim, xscale=c(0.5, max(cat.var, na.rm=TRUE) + 0.5), clip=TRUE))
    if(type == "boxplot")   panel.bwplot(cat.var, cont.var, horizontal=horiz, col=boxplot.pars$box.col, pch=boxplot.pars$box.pch, gp=grid::gpar(box.umbrella=list(col=boxplot.pars$box.col)), 
                                         varwidth=boxplot.pars$varwidth, notch=boxplot.pars$notch, notch.frac=boxplot.pars$notch.frac, fill=box.fill)
    if(type == "stripplot") panel.stripplot(cat.var, cont.var, horizontal=horiz, jitter.data=stripplot.pars$jitter, col=stripplot.pars$col, cex=stripplot.pars$size, pch=stripplot.pars$pch)
    if(type == "violin")    .violin_panel(cat.var, cont.var, horizontal=horiz, varwidth=boxplot.pars$varwidth, border=boxplot.pars$box.col, col=box.fill)
  }
  grid::popViewport(1)
  trellis.par.set(name="box.rectangle", value=list(col=old.color))
  trellis.par.set(name="box.umbrella",  value=list(col=old.color))
  trellis.par.set(name="box.dot",       value=list(col=old.color))
  trellis.par.set(name="plot.symbol",   value=list(col=old.color))
}

#' @importFrom lattice "current.panel.limits" "llines" "ltext" "trellis.par.get"
.contour_panel <- function(x, y, zz, density.pars) {
  levs         <- pretty(zz, density.pars$nlevels)
  cLines       <- grDevices::contourLines(x, y, zz, levels=levs)
  text         <- trellis.par.get("add.text")
  tmp          <- list(col=density.pars$dcol, alpha=text$alpha, cex=text$cex/2, fontfamily=text$fontfamily, fontface=text$fontface, font=text$font, style=density.pars$label.style)
  labels       <- c(tmp, list(labels=format(levs, trim=TRUE)))
  ux           <- sort(unique(x[!is.na(x)]))
  uy           <- sort(unique(y[!is.na(y)]))
  cpl          <- current.panel.limits(unit = "cm")
  asp          <- diff(cpl$ylim)/diff(cpl$xlim)
  ccount       <- 0L
  for(val in cLines) {
    ccount     <- ccount  + 1L
    llines(val, col=density.pars$dcol[ccount], identifier=paste("levelplot", "line", ccount, sep="."))
    if(length(val$x) > 5 && isTRUE(density.pars$show.labels)) {
      slopes   <- diff(val$y)/diff(val$x)
      switch(EXPR=labels$style,
             flat=   {
        txtloc <- which.min(abs(slopes))
        rotang <- 0
      },             {
        rx     <- range(ux)
        ry     <- range(uy)
        depth  <- pmin(pmin(val$x - rx[1L], rx[2L] - val$x)/diff(rx), pmin(val$y - ry[1L], ry[2L] - val$y)/diff(ry))
        if(labels$style  == "align" | 
           depth[txtloc  <- which.min(abs(slopes))] < 0.5) {
          txtloc         <- min(which.max(depth), length(slopes))
          rotang         <- atan(asp * slopes[txtloc] * diff(rx)/diff(ry)) * 180/base::pi  
        } else rotang    <- 0
      })
      ltext(labels$labels[match(val$level, levs)], srt=rotang, adj=c(0.5, 0), col=labels$col[ccount], alpha=labels$alpha, cex=labels$cex, font=labels$font, fontfamily=labels$fontfamily,
            fontface=labels$fontface, x=0.5 * (val$x[txtloc] + val$x[txtloc + 1L]), y=0.5 * (val$y[txtloc] + val$y[txtloc + 1L]), identifier=paste("levelplot", "label", ccount, sep="."))
    }
  }
}

#' @importFrom matrixStats "rowLogSumExps"
.density_panel <- function(dat, dimens, res, density.pars, axis.pars, xpos, ypos, buffer, outer.rot, bg, border) {
  pars         <- res$parameters
  modelName    <- res$modelName
  G            <- res$G
  n            <- res$n
  sigma        <- list(modelName=switch(EXPR=modelName, EEE=, VVV="EEV", modelName), d=2, G=G)
  if(is.element(modelName, c("EEI", "EII", "VII")))           {
    switch(EXPR=modelName, EII=, VII=      {
      sigma    <- c(sigma, list(sigmasq=pars$variance$sigmasq))
    }, EEI=      {
      sigma    <- c(sigma, list(Sigma=pars$variance$Sigma[dimens,dimens]))
    } )
  } else         {
    sigma      <- c(sigma, list(sigma=array(dim=c(2, 2, G))))
    for(k in seq_len(G)) sigma$sigma[,,k] <- pars$variance$sigma[dimens,dimens,k]
  }
  range1       <- range(dat[,1L], na.rm=TRUE)
  range2       <- range(dat[,2L], na.rm=TRUE)
  xlim         <- range1 + c(-buffer * (diff(range1)), buffer * (diff(range1)))
  ylim         <- range2 + c(-buffer * (diff(range2)), buffer * (diff(range2)))
  lx           <- density.pars$grid.size[1L]
  ly           <- density.pars$grid.size[2L]
  x            <- .grid_1(n=lx, range=xlim, edge=TRUE)
  y            <- .grid_1(n=ly, range=ylim, edge=TRUE)
  xy           <- .grid_2(x, y)
  xyn          <- nrow(xy)
  Vinv         <- pars$Vinv
  noise        <- !is.null(Vinv)
  GN           <- G + noise
  gate         <- attr(res, "Gating")
  expx         <- attr(res, "Expert")
  ltau         <- .mat_byrow(pars$lpro, nrow=xyn, ncol=GN)
  mu           <- if(isTRUE(expx)) pars$fits[,dimens,,drop=FALSE] else pars$mean[dimens,, drop=FALSE]
  if(expx)        {
    if(noise)     {
      den      <- cbind(Reduce("+", lapply(seq_len(n), function(i) MoE_dens(data=xy, mus=mu[i,,], sigs=sigma)))/n, log(Vinv)) + ltau
    } else den <- Reduce("+",       lapply(seq_len(n), function(i) MoE_dens(data=xy, mus=mu[i,,], sigs=sigma, Vinv=Vinv)))/n  + ltau
  } else if(gate) {
    den        <- MoE_dens(data=xy, mus=mu, sigs=sigma, Vinv=Vinv) + ltau
  }   else den <- MoE_dens(data=xy, mus=mu, sigs=sigma, log.tau=ltau, Vinv=Vinv)
  zz           <- matrix(exp(rowLogSumExps(den)), nrow=lx, ncol=ly, byrow=FALSE)
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim))
  .draw_axis(x=x, y=y, axis.pars=axis.pars, xpos=xpos, ypos=ypos, cat.labels=NULL, horiz=NULL, xlim=xlim, ylim=ylim, outer.rot=outer.rot)
  grid::popViewport(1)
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim, clip=TRUE))
  grid::grid.rect(gp=grid::gpar(fill=bg, col=border))
  .contour_panel(x, y, zz, density.pars)
  grid::popViewport(1)
}

#' @importFrom lattice "llines" "panel.histogram"
.diag_panel <- function(x, varname, diag.pars, hist.col, axis.pars, xpos, ypos, buffer, index, i, res, col, outer.rot) {
  x         <- x[!is.na(x)]
  drange    <- range(as.numeric(x), na.rm=TRUE)
  xlim      <- drange + c(-buffer * (diff(drange)), buffer * (diff(drange)))
  ylim      <- xlim
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim))
  what      <- if(is.factor(x)) x else as.numeric(x)
  .draw_axis(x=what, y=what, axis.pars=axis.pars, xpos=xpos, ypos=ypos, cat.labels=NULL, horiz=NULL, xlim=xlim, ylim=ylim, outer.rot=outer.rot)
  grid::popViewport(1)
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim, clip=TRUE))
  if(!any(diag.pars$show.hist, diag.pars$show.dens && !is.factor(x))) {
    grid::grid.rect()
    grid::grid.text(varname, 0.5, 0.5, gp=grid::gpar(fontsize=diag.pars$fontsize, fontface=2))
  }
  grid::popViewport(1)
  if(diag.pars$show.hist)  {
    if(!is.factor(x)) {
      grid::pushViewport(grid::viewport(xscale=xlim, yscale=c(0, 100), clip=TRUE))
      panel.histogram(as.numeric(x), breaks=NULL, type="percent", col=if(diag.pars$show.dens) "white" else if(index == 1)  hist.col[[index]] else hist.col[index])
    } else {
      grid::pushViewport(grid::viewport(xscale=c(min(as.numeric(x), na.rm=TRUE) - 1, max(as.numeric(x), na.rm=TRUE)  + 1), yscale=c(0, 100), clip=TRUE))
      tabx <- table(x)
      show.counts    <- if(isTRUE(diag.pars$show.counts)) as.numeric(tabx) else FALSE
      .barchart_panel(seq_along(tabx), 100 * tabx/sum(tabx), horizontal=FALSE, col=if(index == 1) hist.col[[index]] else hist.col[index], show.counts=show.counts, fontsize=diag.pars$fontsize)
    }
  }
  if(diag.pars$show.dens) {
    if(!is.factor(x)) {
      if(!diag.pars$show.hist) grid::pushViewport(grid::viewport(xscale=xlim, yscale=c(0, 1), clip=TRUE))
      pars         <- res$parameters
      modelName    <- pars$varianceX$modelName
      G            <- res$G
      n            <- res$n
      sigma        <- list(modelName=modelName, d=1, G=G)
      if(is.element(modelName, c("EEI", "EII", "EEE", "VII", "E", "V"))) {
        switch(EXPR=modelName, E=, V=, EII=, VII=      {
          sigma    <- c(sigma, list(sigmasq=pars$varianceX$sigmasq))
        }, EEI=, EEE=     {
          sigma    <- c(sigma, list(sigmasq=pars$varianceX$Sigma[i,i]))
        } )
        sigma$modelName  <- switch(EXPR=modelName, V=, VII="V", "E")
      } else         {
        sigma      <- c(sigma, list(sigmasq=array(dim=c(1, 1, G))))
        for(k in seq_len(G)) sigma$sigmasq[,,k] <- pars$varianceX$sigma[i,i,k]
        sigma$modelName  <- "V"
      }
      if((lx       <- diag.pars$grid.size) > n) {
        rx         <- range(x, na.rm=TRUE)
        xlim       <- rx + c(-buffer * (diff(rx)), buffer * (diff(rx)))
        xd         <- .grid_1(n=lx, range=xlim, edge=TRUE)
        xn         <- length(xd)  
      } else {
        xd         <- x
        xn         <- n
      }
      Vinv         <- pars$Vinv
      noise        <- !is.null(Vinv)
      GN           <- G + noise
      gate         <- attr(res, "Gating")
      expx         <- attr(res, "Expert")
      ltau         <- .mat_byrow(pars$lpro, nrow=xn, ncol=GN)
      mu           <- if(isTRUE(expx)) pars$fits[,i,,drop=FALSE] else if(res$d == 1) pars$mean else pars$mean[i,, drop=FALSE]
      if(expx)        {
        if(noise)     {
          den      <- cbind(Reduce("+", lapply(seq_len(n), function(i) MoE_dens(data=xd, mus=mu[i,,], sigs=sigma)))/n, log(Vinv)) + ltau
        } else den <- Reduce("+",       lapply(seq_len(n), function(i) MoE_dens(data=xd, mus=mu[i,,], sigs=sigma, Vinv=Vinv)))/n  + ltau
      } else if(gate) {
        den        <- MoE_dens(data=xd, mus=mu, sigs=sigma, Vinv=Vinv) + ltau
      }   else den <- MoE_dens(data=xd, mus=mu, sigs=sigma, log.tau=ltau, Vinv=Vinv)
      zz           <- matrix(exp(rowLogSumExps(den)), nrow=xn, ncol=1L)
      den          <- exp(den)
      oo           <- order(xd)
      if(diag.pars$show.hist) {
        endpoints    <- range(x, finite=TRUE)
        nint         <- round(log2(length(x)) + 1)
        counts       <- graphics::hist(x, breaks=endpoints[1L] + diff(endpoints) * (0L:nint)/nint, plot=FALSE)$counts
        cmax         <- 100 * max(counts)/n
      } else {
        cmax         <- 0.9 - buffer
      }
      dmax           <- max(den, zz)
      for(g in seq_len(GN)) {
        llines(xd[oo], den[oo,g] * cmax/dmax, col=col[g])
      }
      llines(xd[oo], zz[oo] * cmax/dmax, col=hist.col[1L], lwd=2)
    }
  }
  if(any(diag.pars$show.hist, diag.pars$show.dens && !is.factor(x))) {
    grid::grid.text(varname, 0.5, 0.9, gp=grid::gpar(fontsize=diag.pars$fontsize))
    grid::popViewport(1)
  }
}

.heat_legend <- function(data, col = NULL, cex.lab = 1, ...) {
  if(length(cex.lab) > 1    || 
     (!is.numeric(cex.lab)  ||
     cex.lab <= 0))                                   stop("Invalid 'cex.lab' supplied", call.=FALSE)
  if(!is.numeric(data))                               stop("'data' must be numeric",     call.=FALSE)
  if(missing(col))   {
    col      <- grDevices::heat.colors(30L, rev=TRUE)
  } else if(!all(.is_cols(col)))                      stop("Invalid 'col' colour palette supplied", call.=FALSE)
  bx         <- graphics::par("usr")
  xpd        <- graphics::par()$xpd
  box.cx     <- c(bx[2L] + (bx[2L]  - bx[1L])/1000, bx[2L] + (bx[2L] - bx[1L])/1000 + (bx[2L] - bx[1L])/50)
  box.cy     <- c(bx[3L],   bx[3L])
  box.sy     <- (bx[4L]  -  bx[3L]) / length(col)
  xx         <- rep(box.cx, each    = 2L)
  graphics::par(xpd = TRUE)
  
  for(i in seq_along(col)) {
    yy       <- c(box.cy[1L] + (box.sy * (i - 1L)),
                  box.cy[1L] + (box.sy * (i)),
                  box.cy[1L] + (box.sy * (i)),
                  box.cy[1L] + (box.sy * (i - 1L)))
    graphics::polygon(xx, yy, col   =  col[i], border = col[i])
  }
  
  graphics::par(new = TRUE)
  yrange     <- range(data, na.rm = TRUE)
  base::plot(0, 0, type = "n",  ylim = yrange, yaxt = "n", ylab = "", xaxt = "n", xlab = "", frame.plot = FALSE)
  graphics::axis(side = 4, las = 2, tick = FALSE, line = 0.1, cex.axis = cex.lab)
  suppressWarnings(graphics::par(xpd = xpd))
    invisible()
}

.is_cols      <- function(cols) {
    vapply(cols, function(x)    { tryCatch(is.matrix(grDevices::col2rgb(x)), error = function(e) FALSE) }, logical(1L))
}

.reorder_hc   <- function(sim, ...) {
  dis         <- 1 - sim
  if(is.matrix(dis)) {
    disd      <- stats::as.dist(dis)
  } else  {
    disd      <- dis
    dis       <- as.matrix(dis)
  }
  n           <- nrow(dis)
    if(n <= 2) seq_len(n) else .hc_reorder(stats::hclust(disd, ...), dis, ...)$order
}

.hc_reorder   <- function(x, dis, ...) {
  dis         <- as.matrix(dis)
  merges      <- x$merge
  n           <- nrow(merges)
  endpoints   <- matrix(0L, n, 2L)
  dir         <- matrix(1L, n, 2L)
  for(i in seq_len(n)) {
    j    <- merges[i,1L]
    k    <- merges[i,2L]
    if((j < 0) && (k   < 0)) {
      endpoints[i,1L] <- -j
      endpoints[i,2L] <- -k
    } else if(j  < 0)  {
      j  <- -j
      endpoints[i,1L] <- j
      e1 <- endpoints[k,1L]
      e2 <- endpoints[k,2L]
      if(dis[j,e1] < dis[j,e2]) {
        endpoints[i,2L] <- e2
      } else     {
        endpoints[i,2L] <- e1
        dir[i,2L] <- -1L
      }
    }
    else if(k  < 0)    {
      k  <- -k
      endpoints[i,2L]   <- k
      e1 <- endpoints[j,1L]
      e2 <- endpoints[j,2L]
      if(dis[k,e1] < dis[k,e2]) {
        endpoints[i,1L] <- e2
        dir[i,1L] <- -1L
      } else    {
        endpoints[i,1L] <- e1
      }
    } else {
      ek1     <- endpoints[k,1L]
      ek2     <- endpoints[k,2L]
      ej1     <- endpoints[j,1L]
      ej2     <- endpoints[j,2L]
      d11     <- dis[ej1,ek1]
      d12     <- dis[ej1,ek2]
      d21     <- dis[ej2,ek1]
      d22     <- dis[ej2,ek2]
      dmin    <- min(d11, d12, d21, d22)
      if(dmin == d21)    {
        endpoints[i,1L] <- ej1
        endpoints[i,2L] <- ek2
      } else if(dmin    == d11) {
        endpoints[i,1L] <- ej2
        endpoints[i,2L] <- ek2
        dir[i,1L] <- -1L
      } else if(dmin    == d12) {
        endpoints[i,1L] <- ej2
        endpoints[i,2L] <- ek1
        dir[i,1L] <-
        dir[i,2L] <- -1L
      }
      else {
        endpoints[i,1L] <- ej1
        endpoints[i,2L] <- ek1
        dir[i,2L] <- -1L
      }
    }
  }
  for(i in n:2L)   {
    if(dir[i,1L]  == -1L)       {
      m  <- merges[i,1L]
      if(m > 0)    {
        m1    <- merges[m,1L]
        merges[m,1L] <- merges[m,2L]
        merges[m,2L] <- m1
        if(dir[m,1L] == dir[m,2L]) {
          dir[m,]    <- -dir[m,]
        }
      }
    }
    if(dir[i,2L] == -1L)        {
      m  <- merges[i,2L]
      if(m > 0)   {
        m1    <- merges[m,1L]
        merges[m,1L] <- merges[m,2L]
        merges[m,2L] <- m1
        if(dir[m,1L] == dir[m,2L]) {
          dir[m,]    <- -dir[m,]
        }
      }
    }
  }
  clusters <- as.list(seq_len(n))
  for(i in seq_len(n))           {
    j    <- merges[[i,1L]]
    k    <- merges[[i,2L]]
    clusters[[i]]    <- if((j < 0) && (k < 0)) c(-j,-k) else if(j < 0) c(-j, clusters[[k]]) else if(k < 0) c(clusters[[j]], -k) else c(clusters[[j]], clusters[[k]])
  }
  x1     <- x
  x1$merge    <- merges
  x1$order    <- clusters[[n]]
    return(x1)
}

#' @importFrom vcd "strucplot"
.mosaic_panel <- function(x, y, mosaic.pars, mosaic.inner, axis.pars, xpos, ypos, outer.rot, bg, mfill) {
  if(isTRUE(mfill))        {
    mosaic.pars$gp$fill   <- mosaic.inner
    mosaic.pars$gp$col    <- bg
  } else       {
    mosaic.pars$gp$fill   <- bg
    mosaic.pars$gp$col    <- mosaic.inner
  }
  if(!is.null(xpos)       && !is.null(ypos)) {
    suppressWarnings(strucplot(table(y, x), margins=c(0, 0, 0, 0), newpage=FALSE, pop=FALSE, keep_aspect_ratio=FALSE, shade=mosaic.pars$shade, legend=FALSE, gp=mosaic.pars$gp,
                     gp_args=mosaic.pars$gp_args, labeling_args=list(tl_labels=c(xpos, !ypos), gp_labels=mosaic.pars$gp_labels, varnames=c(FALSE, FALSE), rot_labels=rep(outer.rot, 2))))
  } else if(is.null(xpos) && is.null(ypos))  {
      suppressWarnings(strucplot(table(y, x), margins=c(0, 0, 0, 0), shade=mosaic.pars$shade, legend=FALSE, gp=mosaic.pars$gp, gp_args=mosaic.pars$gp_args, newpage=FALSE, pop=FALSE, keep_aspect_ratio=FALSE, labeling=NULL))
  } else if(is.null(xpos)) {
      suppressWarnings(strucplot(table(y, x), margins=c(0, 0, 0, 0), newpage=FALSE, pop=FALSE, keep_aspect_ratio=FALSE, shade=mosaic.pars$shade, legend=FALSE, gp=mosaic.pars$gp, gp_args=mosaic.pars$gp_args,
                       labeling_args=list(labels=c(TRUE, FALSE), tl_labels=c(ypos, FALSE), gp_labels=mosaic.pars$gp_labels, varnames=c(FALSE, FALSE), rot_labels=rep(outer.rot, 2))))
  } else suppressWarnings(strucplot(table(y, x), margins=c(0, 0, 0, 0), newpage=FALSE, pop=FALSE, keep_aspect_ratio=FALSE, shade=mosaic.pars$shade, legend=FALSE, gp=mosaic.pars$gp, gp_args=mosaic.pars$gp_args,
                          labeling_args=list(labels=c(FALSE, TRUE), tl_labels=c(FALSE, !xpos), gp_labels=mosaic.pars$gp_labels, varnames=c(FALSE, FALSE), rot_labels=rep(outer.rot, 2))))
}

#' @importFrom lattice "panel.abline"
.scatter_panel <- function(x, y, type, scatter.pars, axis.pars, xpos, ypos, buffer, z, G, res, dimens, outer.rot, bg, uncertainty, mvn.col, mvn.type, border) {
  xrange  <- range(x, na.rm=TRUE)
  yrange  <- range(y, na.rm=TRUE)
  xlim    <- xrange + c(-buffer * (diff(xrange)), buffer * (diff(xrange)))
  ylim    <- yrange + c(-buffer * (diff(yrange)), buffer * (diff(yrange)))
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim))
  .draw_axis(x=x, y=y, axis.pars=axis.pars, xpos=xpos, ypos=ypos, cat.labels=NULL, horiz=NULL, xlim=xlim, ylim=ylim, outer.rot=outer.rot)
  grid::popViewport(1)
  grid::pushViewport(grid::viewport(xscale=xlim, yscale=ylim, clip=TRUE))
  grid::grid.rect(gp=grid::gpar(fill=bg, col=border))
  uncert  <- !all(is.na(uncertainty))
  grid::grid.points(x, y, pch=scatter.pars$pch, size=scatter.pars$size, gp=grid::gpar(col=if(uncert) uncertainty$col else scatter.pars$col, cex=if(uncert) uncertainty$cex else 1))
  switch(EXPR=type, ellipses= {
    mu    <- array(res$parameters$mean[dimens,], c(2, G))
    sigma <- array(res$parameters$varianceX$sigma[dimens, dimens,], c(2, 2, G))
    for(g in rev(seq_len(G))) try(.mvn2D_panel(mu=mu[,g], sigma=sigma[,,g], k=15L, col=if(mvn.type == "inner") c("grey30", mvn.col[g], mvn.col[g]) else if(mvn.type == "outer") c(mvn.col[g], "grey30", "grey30") else if(mvn.type == "both") rep(mvn.col[g], 3L)), silent=TRUE)
  }, lm=   {
    for(g in seq_len(G)) {
      xy.lm <- stats::lm(y ~ x, weights=z[,g])
      panel.abline(xy.lm$coef[1L], xy.lm$coef[2L], col=scatter.pars$eci.col[g], lwd=1)
    }
  }, ci=   {
    for(g in seq_len(G)) {
      xy.lm <- stats::lm(y ~ x, weights=z[,g])
      xy    <- data.frame(x = seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE), length.out=25L))
      yhat  <- stats::predict(xy.lm, newdata=xy, interval="confidence")
      ci    <- data.frame(lower=yhat[, "lwr"], upper=yhat[, "upr"])
      panel.abline(xy.lm$coef[1L], xy.lm$coef[2L], col=scatter.pars$eci.col[g], lwd=1)
      grid::grid.lines(x=xy$x, y=ci$lower, default.units="native", gp=grid::gpar(col=scatter.pars$eci.col[g], lwd=1, lty=2))
      grid::grid.lines(x=xy$x, y=ci$upper, default.units="native", gp=grid::gpar(col=scatter.pars$eci.col[g], lwd=1, lty=2))
    }
  })
  grid::popViewport(1)
}

#' @importFrom lattice "current.panel.limits" "trellis.grobname" "trellis.par.get"
.violin_panel  <- function(x, y, box.ratio = 1, box.width = box.ratio/(1 + box.ratio), horizontal = TRUE, alpha = plot.polygon$alpha,
                           border = plot.polygon$border, lty = plot.polygon$lty, lwd = plot.polygon$lwd, col = plot.polygon$col,
                           varwidth = FALSE, bw = NULL, adjust = NULL, kernel = NULL, window = NULL, width = NULL, n = 50,
                           from = NULL, to = NULL, cut = NULL, na.rm = TRUE, ..., identifier = "violin") {
  if(all(is.na(x) | is.na(y))) return()
  x           <- as.numeric(x)
  y           <- as.numeric(y)
  plot.polygon           <- trellis.par.get("plot.polygon")
  darg        <- list()
  darg$bw     <- bw
  darg$adjust <- adjust
  darg$kernel <- kernel
  darg$window <- window
  darg$width  <- width
  darg$n      <- n
  darg$from   <- from
  darg$to     <- to
  darg$cut    <- cut
  darg$na.rm  <- na.rm
  my.density  <- function(x) {
    answer    <- try(do.call(stats::density, c(list(x=x), darg)), silent=TRUE)
    if(!inherits(answer, "try-error"))     answer              else list(x=rep(x[1L], 3L), y=c(0L, 1L, 0L))
  }
  numeric.list           <- if(horizontal) split(x, factor(y)) else split(y, factor(x))
  levels.fos  <- as.numeric(names(numeric.list))
  d.list      <- lapply(numeric.list, my.density)
  dx.list     <- lapply(d.list, "[[", "x")
  dy.list     <- lapply(d.list, "[[", "y")
  max.d       <- vapply(dy.list, max, numeric(1L))
  if(varwidth)   max.d[] <- max(max.d)
  cur.limits  <- current.panel.limits()
  xscale      <- cur.limits$xlim
  yscale      <- cur.limits$ylim
  height      <- box.width
  col         <- if(length(col) == 1)      rep(col, n)         else col
  if(horizontal) {
    for(i in seq_along(levels.fos)) {
      if(is.finite(max.d[i])) {
        grid::pushViewport(grid::viewport(y=grid::unit(levels.fos[i], "native"), height=grid::unit(height, "native"), yscale=c(max.d[i] * c(-1, 1)), xscale=xscale))
        grid::grid.polygon(x=c(dx.list[[i]], rev(dx.list[[i]])), y=c(dy.list[[i]], -rev(dy.list[[i]])), default.units="native",
                           name=trellis.grobname(identifier, type="panel", group=0), gp=grid::gpar(fill=col[i], col=border, lty=lty, lwd=lwd, alpha=alpha))
        grid::popViewport()
      }
    }
  }
  else {
    for(i in seq_along(levels.fos)) {
      if(is.finite(max.d[i])) {
        grid::pushViewport(grid::viewport(x=grid::unit(levels.fos[i], "native"), width=grid::unit(height, "native"), xscale=c(max.d[i] * c(-1, 1)), yscale=yscale))
        grid::grid.polygon(y=c(dx.list[[i]],  rev(dx.list[[i]])),
                           x=c(dy.list[[i]], -rev(dy.list[[i]])), default.units="native",
                           name=trellis.grobname(identifier, type="panel", group=0), gp=grid::gpar(fill=col[i], col=border, lty=lty, lwd=lwd, alpha=alpha))
        grid::popViewport()
      }
    }
  }
    invisible()
}

.bar_code       <- function(x, outer.margins = list(bottom = grid::unit(2, "lines"), left = grid::unit(2, "lines"), top = grid::unit(2, "lines"), right = grid::unit(2, "lines")),
                            horizontal = TRUE, xlim = NULL, nint = 0, main = "", xlab = "", labelloc = TRUE, axisloc = TRUE, labelouter = FALSE,  newpage = TRUE,
                            fontsize = 9, ptsize = grid::unit(0.25, "char"), ptpch = 1, bcspace = NULL, use.points = FALSE, buffer = 0.025, log = FALSE, outerbox = "black", col = "black") {
  if(!is.null(labelloc)) {
    if(labelloc == "right" || labelloc == "top")
      labelloc  <- FALSE
    if(labelloc == "left"  || labelloc == "bottom")
      labelloc  <- TRUE
  }
  if(!is.null(axisloc))  {
    if(axisloc  == "right" || axisloc  == "top")
      axisloc   <- FALSE
    if(axisloc  == "left"  || axisloc  == "bottom")
      axisloc   <- TRUE
  }
  x             <- if(is.vector(x) && !inherits(x, "list")) list(x) else x
  names(x)      <- if(is.null(names(x))) as.character(seq_along(x)) else names(x)
  x             <- if(is.matrix(x))                as.data.frame(x) else x
  xlaboffset    <- grid::unit(2.5, "lines")
  if(newpage)      grid::grid.newpage()
  grid::grid.text(main, 0.5, grid::unit(1, "npc") - grid::unit(1, "lines"), gp=grid::gpar(fontface="bold"))
  if(!is.null(axisloc) && !axisloc && main != ""  && !labelouter) outer.margins$top <- outer.margins$top + grid::unit(2, "lines")
  if(!is.null(axisloc) && xlab != ""  && !labelouter)   {
    if(axisloc)  {
      if(horizontal) { outer.margins$bottom       <- outer.margins$bottom + grid::unit(1.5, "lines")
      } else outer.margins$top     <- outer.margins$top + grid::unit(1.5, "lines")
    } else       {
      if(horizontal) { outer.margins$top          <- outer.margins$top    + grid::unit(1.5, "lines")
      } else     {
        outer.margins$left  <- outer.margins$left  + grid::unit(1.5, "lines")
        outer.margins$right <- outer.margins$right - grid::unit(1.5, "lines")
      }
    }
  }
  if(horizontal) {
    thisangle   <- 0
    thisjust    <- c("left", "bottom")
  } else         {
    thisangle   <- 90
    thisjust    <- c("left", "top")
    grid::pushViewport(grid::viewport(x=0, y=0, width=grid::convertHeight(grid::unit(1, "npc"), "inches"), height=grid::convertWidth(grid::unit(1, "npc"), "inches"), just=c("left", "bottom")))
  }
  outer.margins <- if(labelouter) list(bottom=grid::unit(0, "lines"), left=grid::unit(0, "lines"), top=grid::unit(0, "lines"), right=grid::unit(0, "lines")) else outer.margins
  vp.main       <- grid::viewport(x=outer.margins$left, y=outer.margins$bottom,
                                  width=grid::unit(1,  "npc") - outer.margins$right - outer.margins$left,
                                  height=grid::unit(1, "npc") - outer.margins$top   - outer.margins$bottom,
                                  just=thisjust, angle=thisangle, name="main", clip="off")
  grid::pushViewport(vp.main)
  grid::grid.rect(gp=grid::gpar(col=outerbox))
  .barcode_panel(x, horizontal=horizontal, nint=nint, xlim=xlim, labelloc=labelloc, labelouter=labelouter,
                 fontsize=fontsize, ptsize=ptsize, bcspace=bcspace, use.points=use.points, xlab=xlab,
                 xlaboffset=xlaboffset, axisloc=axisloc, buffer=buffer, log=log, col=col)
  grid::popViewport(1)
  if(!horizontal) grid::popViewport(1)
}

#' @importFrom lattice "current.panel.limits" "panel.abline" "panel.rect" "panel.text" "trellis.par.get"
.barchart_panel <- function(x, y, box.ratio = 1, box.width = box.ratio/(1 + box.ratio),
                            horizontal = TRUE, origin = NULL, reference = TRUE, stack = FALSE,
                            groups = NULL, col = if(is.null(groups)) plot.polygon$col else superpose.polygon$col,
                            border = if(is.null(groups)) plot.polygon$border else superpose.polygon$border,
                            lty = if(is.null(groups)) plot.polygon$lty else superpose.polygon$lty,
                            lwd = if(is.null(groups)) plot.polygon$lwd else superpose.polygon$lwd,
                            show.counts = FALSE, ..., fontsize = 9, identifier = "barchart") {
  plot.polygon      <- trellis.par.get("plot.polygon")
  superpose.polygon <- trellis.par.get("superpose.polygon")
  reference.line    <- trellis.par.get("reference.line")
  keep          <- (function(x, y, groups, subscripts, ...) {
    !is.na(x) & !is.na(y) & if(is.null(groups)) TRUE else !is.na(groups[subscripts]) })(x=x, y=y, groups=groups, ...)
  if(!any(keep))   return()
  x             <- as.numeric(x[keep])
  y             <- as.numeric(y[keep])
  if(!is.null(groups))    {
    groupSub    <- function(groups, subscripts, ...) groups[subscripts[keep]]
    if(!is.factor(groups)) groups <- factor(groups)
    nvals       <- nlevels(groups)
    groups      <- as.numeric(groupSub(groups, ...))
  }
  if(horizontal) {
    if(is.null(groups))   {
      if(is.null(origin)) {
        origin  <- current.panel.limits()$xlim[1L]
        reference   <- FALSE
      }
      height    <- box.width
      if(reference)  panel.abline(v=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      panel.rect(x=rep(origin, length(y)), y=y, height=rep(height, length(y)), width=x - origin, border=border, col=col, lty=lty, lwd=lwd, just=c("left", "centre"), identifier=identifier)
    } else if(stack)      {
      if(!is.null(origin) && origin != 0)             warning("'origin' forced to 0 for stacked bars\n", call.=FALSE)
      col       <- rep(col,    length.out=nvals)
      border    <- rep(border, length.out=nvals)
      lty       <- rep(lty,    length.out=nvals)
      lwd       <- rep(lwd,    length.out=nvals)
      height    <- box.width
      if(reference) panel.abline(v=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      for(i in unique(y)) {
        ok      <- y == i
        ord     <- sort.list(groups[ok])
        pos     <- x[ok][ord] > 0
        nok     <- sum(pos, na.rm=TRUE)
        if(nok   > 0) panel.rect(x=cumsum(c(0, x[ok][ord][pos][-nok])), y=rep(i, nok), col=col[groups[ok][ord][pos]], border=border[groups[ok][ord][pos]], lty=lty[groups[ok][ord][pos]],
                                 lwd=lwd[groups[ok][ord][pos]], height=rep(height, nok), width=x[ok][ord][pos], just=c("left", "centre"), identifier=paste(identifier, "pos", i, sep="."))
        neg     <- x[ok][ord] < 0
        nok     <- sum(neg, na.rm=TRUE)
        if(nok   > 0) panel.rect(x=cumsum(c(0, x[ok][ord][neg][-nok])), y=rep(i, nok), col=col[groups[ok][ord][neg]], border=border[groups[ok][ord][neg]], lty=lty[groups[ok][ord][neg]],
                                 lwd=lwd[groups[ok][ord][neg]], height=rep(height, nok), width=x[ok][ord][neg], just=c("left", "centre"), identifier=paste(identifier, "neg", i, sep="."))
      }
    } else {
      if(is.null(origin)) {
        origin  <- current.panel.limits()$xlim[1L]
        reference   <- FALSE
      }
      col       <- rep(col,    length.out=nvals)
      border    <- rep(border, length.out=nvals)
      lty       <- rep(lty,    length.out=nvals)
      lwd       <- rep(lwd,    length.out=nvals)
      height    <- box.width/nvals
      if(reference) panel.abline(v=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      for(i in unique(y)) {
        ok      <- y == i
        nok     <- sum(ok, na.rm=TRUE)
        panel.rect(x=rep(origin, nok), y=(i + height * (groups[ok] - (nvals + 1)/2)), col=col[groups[ok]], border=border[groups[ok]], lty=lty[groups[ok]],
                   lwd=lwd[groups[ok]], height=rep(height, nok), width=x[ok] - origin, just=c("left", "centre"), identifier=paste(identifier, "y", i, sep="."))
      }
    }
  } else {
    if(is.null(groups))   {
      if(is.null(origin)) {
        origin  <- current.panel.limits()$ylim[1L]
        reference   <- FALSE
      }
      width     <- box.width
      y.fix     <- y - origin
      fix.y     <- ifelse(!isFALSE(show.counts), 80, 90)
      y.fix     <- if(max(y.fix) > fix.y) (y.fix * fix.y)/max(y.fix) else y.fix
      if(reference) panel.abline(h=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      panel.rect(x=x, y=rep(origin, length(x)), col=col, border=border, lty=lty, lwd=lwd, width=rep(width, length(x)), height=y.fix, just=c("centre", "bottom"), identifier=identifier)
      if(!identical(show.counts, FALSE)) {
        panel.text(x=x, y=y.fix, label=show.counts, adj=c(0.5, -0.5), identifier=identifier, gp=grid::gpar(fontsize=fontsize), cex=0.8)
      }
    } else if(stack) {
      if(!is.null(origin) && origin != 0)             warning("'origin' forced to 0 for stacked bars\n", call.=FALSE)
      col       <- rep(col,    length.out=nvals)
      border    <- rep(border, length.out=nvals)
      lty       <- rep(lty,    length.out=nvals)
      lwd       <- rep(lwd,    length.out=nvals)
      width     <- box.width
      if(reference) panel.abline(h=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      for(i in unique(x)) {
        ok      <- x == i
        ord     <- sort.list(groups[ok])
        pos     <- y[ok][ord] > 0
        nok     <- sum(pos, na.rm=TRUE)
        if(nok   > 0) panel.rect(x=rep(i, nok), y=cumsum(c(0, y[ok][ord][pos][-nok])), col=col[groups[ok][ord][pos]], border=border[groups[ok][ord][pos]], lty=lty[groups[ok][ord][pos]],
                                 lwd=lwd[groups[ok][ord][pos]], width=rep(width, nok), height=y[ok][ord][pos], just=c("centre", "bottom"), identifier=paste(identifier, "pos", i, sep="."))
        neg     <- y[ok][ord] < 0
        nok     <- sum(neg, na.rm=TRUE)
        if(nok   > 0) panel.rect(x=rep(i, nok), y=cumsum(c(0, y[ok][ord][neg][-nok])), col=col[groups[ok][ord][neg]], border=border[groups[ok][ord][neg]], lty=lty[groups[ok][ord][neg]],
                                 lwd=lwd[groups[ok][ord][neg]], width=rep(width, nok), height=y[ok][ord][neg], just=c("centre", "bottom"), identifier=paste(identifier, "neg", i, sep="."))
      }
    } else {
      if(is.null(origin)) {
        origin  <- current.panel.limits()$ylim[1L]
        reference   <- FALSE
      }
      col       <- rep(col,    length.out=nvals)
      border    <- rep(border, length.out=nvals)
      lty       <- rep(lty,    length.out=nvals)
      lwd       <- rep(lwd,    length.out=nvals)
      width     <- box.width/nvals
      if(reference) panel.abline(h=origin, col=reference.line$col, lty=reference.line$lty, lwd=reference.line$lwd, identifier=paste(identifier, "abline", sep="."))
      for(i in unique(x)) {
        ok      <- x == i
        nok     <- sum(ok, na.rm=TRUE)
        panel.rect(x=(i + width * (groups[ok] - (nvals + 1L)/2)), y=rep(origin, nok), col=col[groups[ok]], border=border[groups[ok]], lty=lty[groups[ok]],
                   lwd=lwd[groups[ok]], width=rep(width, nok), height=y[ok] - origin, just=c("centre", "bottom"), identifier=paste(identifier, "x", i, sep="."))
      }
    }
  }
}

.barcode_panel  <- function(x, horizontal = TRUE, xlim = NULL, labelloc = TRUE, axisloc = TRUE, labelouter = FALSE,
                            nint = 0, fontsize = 9, ptsize = grid::unit(0.25, "char"), ptpch = 1, bcspace = NULL, xlab = "",
                            xlaboffset = grid::unit(2.5, "lines"), use.points = FALSE, buffer = 0.025, log = FALSE, col = "black") {
  if(!inherits(x, "list"))                            stop("x must be a list", call.=FALSE)
  K             <- length(x)
  for(i in seq_len(K))   x[[i]] <- x[[i]][!is.na(x[[i]])]
  maxct         <- 0
  ux     <- unlist(x)
  if(is.null(xlim)) {
    drux <- diff(range(ux))
    minx <- min(ux) - buffer * drux
    maxx <- max(ux) + buffer * drux
  } else  {
    minx <- xlim[1L]
    maxx <- xlim[2L]
  }
  xleft  <- grid::unit(1, "strwidth", names(x)[1L])
  for(i in seq_len(K))   {
    y    <- x[[i]]
    if(length(y)    > 0) {
      z  <- if(nint > 0) graphics::hist(y, breaks=pretty(ux, n=nint), plot=FALSE)$counts else table(y)
      maxct     <- max(maxct, max(z))
      xleft     <- max(xleft, grid::unit(1, "strwidth", names(x)[i]))
    }
  }
  maxct  <- ifelse(log, log(maxct + 3), maxct + 3)
  xleft  <- 1.05 * xleft
  if(is.null(labelloc) || !labelloc)    {
    xright      <- xleft
    xleft       <- grid::unit(0, "npc")
    xtextloc    <- grid::unit(1, "npc") - xright
    xtextalign  <- "left"
  } else  {
    xright      <- grid::unit(0, "npc")
    xtextloc    <- xleft
    xtextalign  <- "right"
  }
  if(labelouter) {
    xleft       <- grid::unit(0, "npc")
    xright      <- grid::unit(0, "npc")
    xtextloc    <- grid::unit(ifelse(is.null(labelloc) || !labelloc, 1.02, -1.02), "npc")
  }
  if(is.null(bcspace)) bcspace <- max(0.2, 1.5/(maxct + 1))
  grid::pushViewport(grid::viewport(x=xleft, y=grid::unit(0, "npc"), width=grid::unit(1, "npc") - xleft - xright,
                     height=grid::unit(1, "npc"), xscale=c(minx, maxx), just=c("left", "bottom")))
  if(!is.null(axisloc)) {
    grid::grid.xaxis(main=axisloc, gp=grid::gpar(fontsize = fontsize))
    if(axisloc)  { grid::grid.text(xlab, x=grid::unit(0.5, "npc"), y = grid::unit(0, "npc") - xlaboffset)
    } else         grid::grid.text(xlab, x=grid::unit(0.5, "npc"), y = grid::unit(1, "npc") + xlaboffset)
  }
  grid::popViewport(1)
  for(i in seq_len(K))  {
    y    <- x[[i]]
    if(!is.null(labelloc)) grid::grid.text(names(x)[i], x=xtextloc, y=grid::unit((i - 1L)/K, "npc") + 0.5 * grid::unit(1/K, "npc"), just=xtextalign, gp=grid::gpar(fontsize=fontsize))
    if(nint > 0) {
     zhist <- graphics::hist(y, breaks=pretty(unlist(x), n=nint), plot=FALSE)
     z     <- zhist$counts
     mids  <- zhist$mids
    } else  {
     z     <- table(y)
     mids  <- as.numeric(names(z))
    }
    if(length(mids) > 0) {
      vp.barcode <- grid::viewport(x=xleft, y=grid::unit((i - 1L)/K, "npc") + grid::unit(0.05/K, "npc"), width=grid::unit(1, "npc") - xleft - xright, height=grid::unit(1/K, "npc") * bcspace - grid::unit(0.05/K, "npc"), xscale=c(minx, maxx), yscale=c(0, 1), just=c("left", "bottom"), name="barcode", clip="off")
      grid::pushViewport(vp.barcode)
      grid::grid.segments(grid::unit(mids[z > 0], "native"), 0, gp=grid::gpar(col=col[i]), grid::unit(mids[z > 0], "native"), 1)
      grid::popViewport(1)
      vp.hist    <- grid::viewport(x=xleft, y=grid::unit((i - 1L)/K, "npc") + grid::unit(1/K, "npc") * bcspace, width=grid::unit(1, "npc") - xright - xleft, height=grid::unit(1/K, "npc") - grid::unit(1/K, "npc") * bcspace, xscale=c(minx, maxx), yscale=c(0, 1), just=c("left", "bottom"), name="hist", clip="off")
      grid::pushViewport(vp.hist)
      vp.buffer  <- grid::viewport(x=0, y=0.05, width=1, height=0.9, just=c("left", "bottom"), xscale=c(minx, maxx), yscale=c(0, 1))
      grid::pushViewport(vp.buffer)
      for(j in seq_along(z)) {
        if(z[j] > 1) {
          xx     <- rep(mids[j], z[j] - 1L)
          yy     <- if(log) log(2L + seq_len(z[j] - 1L))/maxct else seq_len(z[j] - 1L)/maxct
          if(use.points) {
            grid::grid.points(grid::unit(xx, "native"), yy, pch=ptpch, size=ptsize, gp=grid::gpar(col=col[i]))
          } else  {
            yy   <- if(log) c(yy, log(2L + z[j])/maxct)        else c(yy, (z[j])/maxct)
            grid::grid.segments(grid::unit(mids[j], "native"), gp=grid::gpar(col=col[i]), grid::unit(1/maxct, "npc"), grid::unit(mids[j], "native"), grid::unit(max(yy), "npc"))
          }
        }
      }
      grid::popViewport(2)
    }
  }
}

.bubble <- function(x, cex = c(0.2, 3), alpha = c(0.1, 1)) {
  cex   <- cex[!is.na(cex)]
  alpha <- alpha[!is.na(alpha)]
  r     <- sqrt(x/pi)
  r     <- (r - min(r, na.rm=TRUE))/(diff(range(r, na.rm=TRUE)) + sqrt(.Machine$double.eps))
  cex   <- r  * diff(range(cex))   + min(cex)
  alpha <- x  * diff(range(alpha)) + min(alpha)
    return(list(cex = cex, alpha = alpha))
}

#' @method plot mclustAIC
#' @importFrom mclust "plot.mclustBIC"
#' @export
plot.mclustAIC  <- function (x, ylab = "AIC", ...) {
    plot.mclustBIC(x, ylab = ylab, ...)
}

#' @method plot mclustDF
#' @importFrom mclust "plot.mclustBIC"
#' @export
plot.mclustDF   <- function (x, ylab = "DF", ...)  {
    plot.mclustBIC(x, ylab = ylab, ...)
}

#' @method plot mclustITER
#' @importFrom mclust "plot.mclustBIC"
#' @export
plot.mclustITER <- function (x, ylab = "Iterations", ...)  {
    plot.mclustBIC(x, ylab = ylab, ...)
}

#' @method plot mclustLoglik
#' @importFrom mclust "plot.mclustBIC"
#' @export
plot.mclustLoglik <- function (x, ylab = "Log-Likelihood", ...)  {
    plot.mclustBIC(x, ylab = ylab, ...)
}

Try the MoEClust package in your browser

Any scripts or data that you put into this service are public.

MoEClust documentation built on Dec. 28, 2022, 2:24 a.m.