R/plots.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

# ---- PLOTS -----

#### for BDataFrequencyHistograms -----------------------------------------------



#' Plot histograms as segments
#'
#' Plot a histograms dataset as a set of "stairs" using matplot.
#'
#' @param x a "BDataHistograms" object, of any class
#' @param probability logical, should the histogram represent density (forced if
#' `x` is of class "BDataDensityHistograms") or frequency?
#' @param xlab label for the x-axis
#' @param ylab label for the y-axis
#' @param ... further arguments to [graphics::matplot()]
#'
#' @return nothing, called for its side effect to produce a plot.
#' Each histogram is represented only as its outline.
#' @export
#' @family BDataHistograms
#' @importFrom methods is
#' @examples
#' x = lapply(rep(1000,4), runif)
#' bdfhx = as.BDataFrequencyHistograms(x, breaks="Sturges",basis=B2Unif01Basis5)
#' plot(bdfhx)
#' bdfhx2 = as.BDataFrequencyHistograms(x, breaks=sqrt(0:10)/sqrt(10),basis=B2Unif01Basis5)
#' plot(bdfhx2, probability=TRUE)
#' # plot(bdfhx2) # this generates a warning
plot.BDataHistograms = function(x, probability = is(x, "BDataDensityHistograms"),
                                         xlab="x", ylab=ifelse(probability, "density", "frequency"),
                                         ...) {
  if(probability | is(x, "BDataDensityHistograms")) {
    xx = clo(unclass(x)) %*% diag(1/diff(getBreaks(x)))
  } else {
    xx = unclass(x)
    diffs = diff(getBreaks(x))
    if(sum((diffs - diffs[1])^2/diffs)>1e-12) warning("plotting frequencies in spite of irregular breaks")
  }
  graphics::matplot(
    x = rep(getBreaks(x), each=2),
    y = rbind(0,apply(xx, 1, rep, each=2),0) ,
    type="l",
    xlab=xlab,
    ylab=ylab,
    ...
  )
  graphics::abline(h=0, col="grey80")
  invisible()
}




#' plot the 'BayesSpaceBasis'
#'
#' @param x 'BayesSpaceBasis' object
#' @param ... further arguments to plot.BDataSplines
#' @param as tells in which version to plot the basis
# ' @param close see [clo], if TRUE closable objects will be normalized
# '     to their closed from.
#' @returns typically some representation of the data plotted
#' @export
#' @method plot BayesSpaceBasis
#' @family BayesSpaceBasis
#' @examples
#' print(B2Unif01Basis5)
#plot.BayesSpaceBasis = function(x,xlab="x",ylab="functional value",...) {
#  plot(x$basis,...)
#}
plot.BayesSpaceBasis = function(x,...,
                                as=c("is","logPDF","PDF","CDF","density","logDensity","clr","coefs")) {
    as=match.arg(as)
    if( as=="is" ) {
        plotFuncs(as.function(x),...,support=getSupport(x))
    } else {
        plot(as.BDataSplines(x),...,as=as)
    }
}


#' Plot a set of spline functions (using visualisation tools from fda)
#'
#' @param x 'BDataSplines' object
#' @param ... further arguments to plot, plotFunc, plotCoefs
#' @param as a string defining how the spline should be plotted. The
#'     default "coef" is a parallel coordinate plot of the
#'     coefficients. Other options include "PDF", "CDF", "density",
#'     "logPDF", "logDensity", and "clr" (centered log-ratio, showing
#'     the additive representation in Bayes space).
#' @param close if true, and as is not plotting the spline, but a
#'     closable BData, this is used as an argument to [clo].
#' @returns nothing or an invisible version of the object actually
#'     plotted
#' @export
#' @family BDataSplines
#' @examples
#' X = BDataSplines(matrix(rnorm(30),ncol=5),basis=B2Unif01Basis5)
#' plot(X, as = "coef")
#' plot(X, as = "PDF")
#' plot(X, as = "clr")  # centered log-ratio representation
#'
plot.BDataSplines = function(x,...,as=c("coef","PDF","CDF","density","logPDF","logDensity","clr"),close=TRUE)
                              {
    switch(match.arg(as),
           coef=plotCoefs(x,...),
           PDF=plot(clo(as.BDataPDFs(x),close=close),...),
           CDF=plot(clo(as.BDataCDFs(x),close=close),...),
           density=plot(clo(as.BDataDensities(x),close=close),...),
           logPDF =plot(clo(as.BDataLogPDFs(x),close=close),...),
           logDensity = plot(clo(as.BDataLogDensities(x),close=close),...),
           clr = plot(clr(x),...)
           )
}

## Plotting of functions on complex supports
### We want to see a matrix with N columns, but we don't know what the sapply
### returned
rearrangeFromSapply3 <- function(x,N) {
    d = dim(x)
    n = length(d)
    reportIfNot(n<4,n,"To many dimensions in function evaluation")
    if( length(d)==3 ) {
        reportIfNot(d[3]==1||d[2]==1,d,"Functioned returned vectorial value")
        dim(x)=c(d[1],d[2]*d[3])
        t(x)
    }
    if( n<2 ) {
        matrix(x,nrow=N)
    } else {
        reportIfNot(nrow(x)==N,nrow(x),N=N,"The number of rows did not fit")
        x
    }
}

#' Helper function to plot multiple functions
#'
#' This function uses a loop over lines and points to plot multiple
#' functions with different parameters into one plot.
#'
#' @param x a list of functions to be plotted, naming to be compatible
#'     with plot. Alternatively it can be a single function or a list
#'     with a single function which gives a matrix of results, where
#'     the columns correspond to the arguments and the rows to the
#'     different results.
#' @param support a Support object telling what to plot
#' @param col a (vector of) plotting colors recycled over functions
#' @param pch a (vector of) plotting chars recycled over functions
#'     used for the discrete support
#' @param lty a (vector of) line types recycled over the functions
#'     used for the lines of the continuous support
#' @param lwd a (vector of) line widths recycled over the function
#'     used for the lines of the continuous support
#' @param nEval the number of discretization points for the continuous
#'     support
#' @param xs the discretization points for the continuous support
#' @param xL candidates for the discretization points for the
#'     continuous part, unallowed points will be removed.
#' @param X the points for the discrete part of the support
#' @param add a boolean, if TRUE functions are added to an existing
#'     plot, if FALSE a new plot is generated.
#' @param xlim the xlim of a new plot
#' @param ylim the ylim of a new plot, if NULL it will be computed
#'     from the range of the functions.
#' @param xlab the label of the x-axis
#' @param ylab the label of the y-axis
#' @param ... further arguments for plot
#' @param close if TRUE the data is plotted in closed from
#'     (i.e. representing a measure of integral 1), if a numeric
#'     value, it is closed to this numeric value. (not yet
#'     implemented: if 0 the centered version is plotted), if FALSE,
#'     the functions are plotted as is. see [clo].
#' @export
#' @importFrom graphics lines
#' @importFrom graphics par
#' @importFrom graphics points
#' @family Support
#' @returns invisible the results of the computation \itemize{ \item x
#'     - the evaluation points for continuous evaluation \item y - the
#'     evaluated values for the continuous evaluation, each function
#'     gets a column \item X - the evaluation points for the discrete
#'     support \item Y - the evaluated values for the discrete
#'     support, each function gets a column }
#'
#' @examples
#' plotFuncs(list(fx(sin(10*x)),fx(cos(10*x))),support=Support(0,pi,0:3),pch=2:3,col=2:3)
#' X = BDataPDFs(list(fx(dnorm(x,0,1)),fx(dnorm(x,0.3,0.5)),fx(dnorm(x,0,1.5))),
#'               basis=B2N01Basis5)
#' plot(X,col=2:4)
#' plot(as.BDataLogPDFs(X),col=2:4)
#' plot(as.BDataDensities(X),col=2:4)
#' plot(as.BDataLogDensities(X),col=2:4)
#'
plotFuncs = function(x,...,
                             support=getSupport(x),
                             col=par("col"),
                             pch=par("pch"),
                             lty=par("lty"),
                             lwd=par("lwd"),
                             nEval=1500,
                             add=FALSE,
                             xlim=getSupport1DPlottingLimits(support),
                             ylim=NULL,
                             xs=xL[! xL %in% X],
                             xL=if(xlim[1]!=xlim[2])
                                    seq(xlim[1],xlim[2],length.out=nEval),
                             X=getSupportAtoms(support),
                             xlab="x",
                     ylab="f(x)",close=FALSE) {
    if( isFALSE(close) )
        funcs=x
    else
        funcs=clo(x,close=close)
    if( is.function(funcs) )
        funcs=list(funcs)
    N = 0
    if( length(xs)> 0 ) {
        y = rearrangeFromSapply3(sapply(funcs,function(f) f(xs),simplify="array"),length(xs))
        N = ncol(y)
    } else
        y = c()
    if( length(X)>0  ) {
        Y = rearrangeFromSapply3(sapply(funcs,function(f) f(X),simplify="array"),length(X))
        if( N == 0 ) {
            N=ncol(Y)
        } else {
            reportIfNot(N==ncol(Y),N=N,dim(Y),"number of values do not match")
        }
    } else
        Y = c()
    if(!add) {
        if( missing(ylim) || is.null(ylim) || any(!is.finite(ylim))) {
            ylim=boundedRange(y,Y,ylim=ylim)
        }
        plot(0,0,type="n",xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,...)
    }
    for(i in 1:N) {
        if( length(xs) > 0 )
            lines(xs,y[,i],lty=warp(lty,i),lwd=warp(lwd,i),col=warp(col,i))
        if( length(X) > 0 )
            points(X,Y[,i],pch=warp(pch,i),col=warp(col,i))
    }
    return(invisible(list(x=xs,y=y,X=X,Y=Y)))
}


#' @describeIn plotFuncs plot method for BDataDensities
#' @export
#' @method plot BDataDensities
plot.BDataDensities = plotFuncs
formals(plot.BDataDensities)["ylab"]=expression(expression(d*X/d*P))

#' @describeIn plotFuncs plot method for BDataLogDensities
#' @export
#' @method plot BDataLogDensities
plot.BDataLogDensities = plotFuncs
formals(plot.BDataLogDensities)["ylab"]=expression(expression(ln(dX/d*P)))

#' @describeIn plotFuncs plot method for BDataPDFs
#' @export
#' @method plot BDataPDFs
plot.BDataPDFs = plotFuncs
formals(plot.BDataPDFs)["ylab"]=expression(expression(ln(pdf(x))==ln(d*X/d*lambda)))


#' @describeIn plotFuncs plot method for BDataLogPDFs
#' @export
#' @method plot BDataLogPDFs
plot.BDataLogPDFs = plotFuncs
formals(plot.BDataPDFs)["ylab"]=expression(expression(pdf(x)==d*X/d*lambda))

# Fix only for the raw routine
formals(plotFuncs)$close=FALSE

#' Plot multiple coefficients in a parallel coordinate plot
#'
#' The function currently uses the same axes,
#'
#' @param x a matrix of coefficients
#' @param type character plotting type: "p" points, "l" lines, "b"
#'     both, "n" nothing
#' @param commonAxis logical, if TRUE all coefficients are plotted on
#'     the same axis. if FALSE individual axes are plotted. FALSE is
#'     not yet implemented.
#' @param sequence gives the sequence of the coefficients to be
#'     plotted.
#' @param col a (vector of) plotting colors recycled over functions
#' @param pch a (vector of) plotting chars recycled over functions and
#'     used for the points
#' @param cex a (vector of) character extends recycled over the data
#'     rows and used for points
#' @param lty a (vector of) line types recycled over the functions
#'     used for the lines of the continuous support
#' @param lwd a (vector of) line widths recycled over the function
#'     used for the lines of the continuous support
#' @param nEval the number of discretization points for the continuous
#'     support
#' @param add a boolean, if TRUE functions are added to an existing
#'     plot, if FALSE a new plot is generated.
#' @param xlim the xlim of a new plot
#' @param ylim the ylim of a new plot, if NULL it will be computed
#'     from the range of the functions.
#' @param xlab the label of the x-axis
#' @param ylab the label of the y-axis
#' @param axes a list of scalings for the individual parallel axes
#'     (not yet implemented)
#' @param xaxt a character specifying the type of xaxis, please see
#'     par or plot for a description. Use "n" to remove it and "s" to plot.
#' @param ... further arguments for plot
#'
#' @returns The coefficient matrix `x` (invisibly), called for the side
#'     effect of producing the parallel coordinate plot.
#' @export
#' @importFrom graphics points
#' @importFrom graphics plot
#' @importFrom graphics axis
#' @examples
#' plotCoefs(matrix(rnorm(20),nrow=5))
#'
#'
#'
plotCoefs = function(x,type="b",...,
                     commonAxis=TRUE,
                     sequence=1:ncol(x),
                     col=par("col"),
                     pch=par("pch"),
                     lty=par("lty"),
                     lwd=par("lwd"),
                     cex=par("cex"),
                     nEval=NULL,
                     add=FALSE,
                     xlim=c(0.5,length(sequence)+0.5),
                     ylim=range(c(x)),
                     xlab="j",
                     axes=rep(list(c(0,1)),ncol(x)),
                     ylab=expression(c["."*j]),xaxt="s") {
    if( !is.matrix(x) )
        x=matrix(x,nrow=1)
    t = 1:length(sequence)
    if( !commonAxis )
        warning("plotCoefs: commonAxis=FALSE not yet implemented.")
    if(!add) {
        if( missing(ylim) || is.null(ylim) ) {
            ylim=range(x)
        }
        plot(0,0,type="n",xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,...,xaxt="n")
        if( !isTRUE(xaxt=="n") )
            graphics::axis(1,t,paste(sequence))
    }
    for(i in 1:nrow(x)) {
        points(t,x[i,sequence],type=type,
              lty=warp(lty,i),lwd=warp(lwd,i),col=warp(col,i),
              pch=warp(pch,i),cex=warp(pch,i))
    }
    return(invisible(x))
}


#' Plotting BDataSample Objects
#'
#' @param x the data
#' @param ... further generic graphic parameters
#' @param as the from of plotting
#' @param pch the plotting character
#' @param method the method argument of the stripchart method.
#' @param add a boolean, if FALSE a new plot is generated.
#' @returns `NULL` invisibly, called for the side effect of producing a plot.
#' @importFrom graphics stripchart
#' @export
#' @method plot BDataSamples
#'
#' @examples
#' s <- BDataSamples(list(rnorm(200), rnorm(200,1)), basis=B2N01Basis5)
#' plot(s)
#'
plot.BDataSamples <- function(x,...,pch=20,method="stack",as=c("points"),
                              add=FALSE) {

    as = match.arg(as)
    switch(as,
           points=graphics::stripchart(x,...,pch=pch,method=method,add=add),
           boxplots=graphics::boxplot(x,...,pch=pch,method=method,add=add)
           )
}


#' Plotting varKernel objects
#'
#' Plots a varKernel object as an image or contour plot.
#'
#' @param x an function of two arguments to be plotted. The function
#'     is expected to be in an outer vectorized way, i.e. that it
#'     returns a matrix, if both arguments are vectors.
#' @param support the support of the first argument of the function
#' @param SD a support discretization object used to discretized the function
#' @param col a (vector of) plotting colors recycled over functions
#' @param pch a (vector of) plotting chars recycled over functions
#'     used for the discrete support
#' @param lty a (vector of) line types recycled over the functions
#'     used for the lines of the continuous support
#' @param cex see par: size of plotting symbols
#' @param type either "both", "image", "contour" defining whether it
#'     should be represented by an image or a contour or both.
#' @param nEval gives the number of point to discretized the continuous
#'     intervall of the suppurt.
#' @param add if FALSE a new plot is generated and axes are drawn
#' @param C the variance kernel evaluated at SD$all
#' @param V the variation kernel evaluated at SD$all
#' @param lwd line width for the countour plot
#' @param nEval the number of discretization points for the continuous
#'     support
#' @param add a boolean, if TRUE functions are added to an existing
#'     plot, if FALSE a new plot is generated.
#' @param xlim the xlim of a new plot
#' @param ylim the ylim of a new plot, if NULL it will be computed
#'     from the range of the functions.
#' @param zlim the limits for plotting the covariance values
#' @param breaks the values for which to change color and where to draw the countour lines.
#' @param ncol the number of colours, it is one less then the number of breaks
#' @param cols the colours to represent the values between the breaks
#' @param xlab the label of the x-axis
#' @param ylab the label of the y-axis y-axis
#' @param zlab the label of the covariance values, not used yet
#' @param ... further arguments for plot
#' @param dWidth discrete Width, the width to draw the discrete covariance values, which would theoretically take only an infinitely think line.
#' @param notColour the colour to plot outside the range of values
#'
#' @returns `NULL` invisibly, called for the side effect of producing an
#'     image and/or contour plot of the kernel.
#' @export
#'
#' @examples
#' C = varKernel(function(x,y) exp(-outer(x,y,"-")^2),basis=B2N01Basis5)
#' plot(C)
#' V = variationKernel(function(x) 1-exp(-outer(x,x,"-")^2),basis=B2N01Basis5)
#' plot(V)



plot.varKernel <- function(x,...,
                           support=getSupport(x),
                           SD=supportDiscretization(support,nEval),
                           C = x(SD$all,SD$all),
                           col=par("col"),
                           pch=par("pch"),
                           lty=par("lty"),
                           lwd=par("lwd"),
                           cex=par("lwd"),
                           type=c("both","image","contour"),
                           nEval=513,
                           add=FALSE,
                           xlim=getSupport1DPlottingLimits(support),
                           ylim=rev(xlim),
                           zlim=range(breaks),
                           breaks=pretty(c(C),ncol),
                           ncol=10,
                           cols=grDevices::hcl.colors(length(breaks)-1, "YlOrRd", rev = TRUE),
                           dWidth=(xlim[2]-xlim[1])/200,
                           xlab="x",
                           ylab=xlab,zlab="c",
                           notColour="#FFFFFF") {
    type <- match.arg(type)
    Ccc  <- C[SD$ix,SD$ix]
    Ccd  <- C[SD$ix,SD$iX]
    Cdd  <- C[SD$iX,SD$iX]
    if( length(breaks)!=length(cols)+1) {
      message(list(ncol=ncol,breaks=breaks,cols=cols,range=range(C)))
      warning("There are ",length(breaks)," breaks and ",length(cols), "colors\n")
    }
         if( is.null(Ccc) ) { # pure discrete
             N   = length(SD$X)
             doAdd=add
             if( type %in% c("both","image") ) {
                 graphics::image(x=SD$Xplot,y=SD$Xplot,
                       z=Cdd,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=cols,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( type %in% c("both","contour") ) {
                 graphics::contour(x=SD$Xplot,y=SD$Xplot,
                         z=Cdd,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=col,lwd=lwd,lty=lty,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( !add ) {
                 graphics::axis(side=1,at=as.numeric(SD$Xplot),labels=SD$Xlab)
                 graphics::axis(side=2,at=as.numeric(SD$Xplot),labels=SD$Xlab)
             }
         } else { # continues
             doAdd=add
             if( type %in% c("both","image") ) {
                 graphics::image(x=SD$x,y=SD$x,
                       z=Ccc,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=cols,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( type %in% c("both","contour") ) {
                 graphics::contour(x=SD$x,y=SD$x,
                         z=Ccc,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=col,lwd=lwd,lty=lty,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( type %in% c("both","image") ) {
                 if(length(SD$Xplot)>0) {
                     for(i in 1:length(SD$Xplot)) {
                         graphics::image(SD$Xplot[i]+c(-dWidth/2,dWidth/2),SD$x,t(Ccd[,c(i,i)]),breaks=breaks,col=cols,add=TRUE)
                         graphics::image(SD$x,SD$Xplot[i]+c(-dWidth/2,dWidth/2),Ccd[,c(i,i)],breaks=breaks,col=cols,add=TRUE)
                     }
                     zcol=c(notColour,cols,notColour)[outer(c(Cdd),breaks,">=")+1]
                     graphics::points(c(outer(SD$Xplot,rep(1,length(SD$Xplot)))),
                            c(outer(rep(1,length(SD$Xplot))),SD$Xplot),
                            pch=pch,col=zcol,cex=cex)
                 }
             } else {
                 graphics::abline(h=SD$Xplot)
                 graphics::abline(v=SD$Xplot)
             }
             if( !add ) {
                 graphics::axis(side=1)
                 graphics::axis(side=2)
                 if( length(SD$Xplot)>0 ) {
                     graphics::axis(side=3,at=as.numeric(SD$Xplot),labels=as.character(SD$Xlab))
                     graphics::axis(side=4,at=as.numeric(SD$Xplot),labels=as.character(SD$Xlab))
                 }
             }
         }
}

#' @describeIn plot.varKernel method for variationKernel
#' @export
#' @method plot variationKernel
plot.variationKernel <- function(x,...,
                                 support=getSupport(x),
                                 SD=supportDiscretization(support,nEval),
                                 V = x(SD$all),
                                 col=par("col"),
                                 pch=par("pch"),
                                 lty=par("lty"),
                                 lwd=par("lwd"),
                                 cex=par("lwd"),
                                 type=c("both","image","contour"),
                                 nEval=513,
                                 add=FALSE,
                                 xlim=getSupport1DPlottingLimits(support),
                                 ylim=rev(xlim),
                                 zlim=range(breaks),
                                 breaks=pretty(V,ncol),
                                 ncol=10,
                                 cols=grDevices::hcl.colors(length(breaks)-1, "YlOrRd", rev = TRUE),
                                 dWidth=(xlim[2]-xlim[1])/200,
                                 xlab="x",
                                 ylab="x",zlab="c",
                                 notColour="#FFFFFF")
{
    Ccc=V[SD$ix,SD$ix]
    Ccd=V[SD$ix,SD$iX]
    Cdd=V[SD$iX,SD$iX]
    xs=SD$x
    type=match.arg(type)
    if( length(breaks)!=length(cols)+1) {
      message(list(ncol=ncol,breaks=breaks,cols=cols,range=range(c(Ccc,Ccd,Cdd))))
      warning("There are ",length(breaks)," breaks and ",length(cols), "colors\n")
    }
    if( is.null(Ccc) ) { # pure discrete
             N   = length(SD$Xplot)
             doAdd=add
             if( type %in% c("both","image") ) {
                 graphics::image(x=SD$Xplot,y=SD$Xplot,
                       z=Cdd,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=cols,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( type %in% c("both","contour") ) {
                 graphics::contour(x=SD$Xplot,y=SD$Xplot,
                         z=Cdd,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=col,lwd=lwd,lty=lty,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( !add ) {
                 graphics::axis(side=1,at=as.numeric(SD$Xplot),labels=as.character(SD$Xlab))
                 graphics::axis(side=2,at=as.numeric(SD$Xplot),labels=as.character(SD$Xlab))
             }
         } else { # continues
             doAdd=add
             if( type %in% c("both","image") ) {
                 graphics::image(x=xs,y=xs,
                       z=Ccc,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=cols,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
             }
             if( type %in% c("both","contour") ) {
                 graphics::contour(x=xs,y=xs,
                         z=Ccc,zlim=zlim,xlim=xlim,ylim=ylim,breaks=breaks,col=col,lwd=lwd,lty=lty,xlab=xlab,ylab=ylab,add=doAdd,xaxt="n",yaxt="n")
                 doAdd=TRUE
                 abline(b=1)
             }
             if( type %in% c("both","image") ) {
                 if(length(SD$Xplot)>0) {
                     for(i in 1:length(SD$Xplot)) {
                         graphics::image(SD$Xplot[i]+c(-dWidth/2,dWidth/2),xs,t(Ccd[,c(i,i)]),breaks=breaks,col=cols,add=TRUE)
                         graphics::image(xs,SD$Xplot[i]+c(-dWidth/2,dWidth/2),Ccd[,c(i,i)],breaks=breaks,col=cols,add=TRUE)
                     }
                     zcol=c(notColour,cols,notColour)[outer(c(Cdd),breaks,">=")+1]
                     graphics::points(c(outer(SD$Xplot,rep(1,length(SD$Xplot)))),
                            c(outer(rep(1,length(SD$Xplot))),SD$Xplot),
                            pch=pch,col=zcol,cex=cex)
                 }
             } else {
                 graphics::abline(h=SD$Xplot)
                 graphics::abline(v=SD$Xplot)
             }
             if( !add ) {
                 graphics::axis(side=1)
                 graphics::axis(side=2)
                 if( length(SD$Xplot)>0 ) {
                     graphics::axis(side=3,at=SD$Xplot,labels=as.character(SD$Xlab))
                     graphics::axis(side=4,at=SD$Xplot,labels=as.character(SD$Xlab))
                 }
             }
         }
}

#' Functional Boxplot for BData objects
#'
#' Currently this function used the [fda::fbplot] function on the
#' evaluated functions in compbination with [depthRank] as depth measure.
#'
#' @param x a BData object
#' @param ... further parameters for [graphics::boxplot] or
#'     [fda::fbplot]
#' @param depth the depth used to identify the middle region,
#'     increasing inwards towards the median.
#' @param depthType if depth is missing, an algorithm to compute
#'     it. The default is "depthquantile", which uses the spatial
#'     depth, but ensures a 50% quantile to be used, "spatialdepth"
#'     used the spatialdepth directly, the others use the
#'     1-[centralityQuantiles] with the given method.
#' @param support the support on which the functions should be plotted
#' @param SD support discretization of the support (see
#'     [supportDiscretization]
#' @param nEval the number of datapoints
#' @param xlim the plotting limits in x direction.
#' @seealso fda::fbplot
#' @returns `NULL` invisibly, called for the side effect of producing a
#'     functional boxplot.
#' @importFrom graphics boxplot
#' @export
#'
#' @examples
#' X = BDataSplines(matrix(rnorm(10*5),ncol=5),basis=B2Unif01Basis5)
#' boxplot(X)
#' boxplot(clr(X),depthType="meandistance")
#' boxplot(clr(X))
#'
boxplot.BData <- function(x,...,depth,
                          depthType=c("depthquantile", "spatialdepth", "meandistance", "median"),
                          support=getSupport(x),
                          SD=supportDiscretization(support,nEval),
                          nEval=513,
                          xlim=getSupport1DPlottingLimits(support)
                          ) {

    .=NULL # just here to avoid warnings
    if( missing(depth) ) {
        depth=switch(match.arg(depthType),
                     spatialdepth=spatialDepth(x),
                     meandistance=1-centralityQuantile(x,distType="meandistance"),
                     median=1-centralityQuantile(x,distType="median"),
                     depthquantile=1-centralityQuantile(x,distType="spatialdepth"),
                     default=stop("Unknown depth type used"))
    }
    xx= SD$x
    XX= SD$X
    fda::fbplot(x%..%c(.(xx)),x=xx,...,xlim=xlim,depth=depth)
    if( length(XX) > 0 )
        graphics::boxplot(t(x %..% c(.(SD$X))),at=as.numeric(SD$Xplot),...,add=TRUE)
}
#' @describeIn boxplot.BData method for BDataSplines uses classical boxplot for the coefficients.
#' @export
boxplot.BDataSplines <- function(x,...
                                 ) {
    graphics::boxplot(unclass(x),...)
}


#' @describeIn boxplot.BData method for BDataSamples uses classical boxplot for the the datasets.
#' @export
boxplot.BDataSamples <- function(x,...
                                 ) {
    graphics::boxplot(unclass(x),...)
}