## Functions for plotting
#' Plot the convex hull of a set of points in 2D.
#'
#' @param pts A matrix with a point in each row.
#' @param drawPoints Draw the points.
#' @param drawLines Draw lines of the facets.
#' @param drawPolygons Fill the hull.
#' @param addText Add text to the points. Currently `coord` (coordinates), `rownames` (rownames)
#' and `both` supported or a vector with text.
#' @param addRays Add the ray defined by `direction`.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of `pts`
#' plus a value greater than on equal zero (minimize objective $i$). If negative, consider the
#' i'th column of `pts` minus a value greater than on equal zero (maximize objective $i$).
#' @param drawPlot Draw the `ggplot`. Set to FALSE if you want to combine hulls in a single plot.
#' @param drawBBoxHull If `addRays` then draw the hull areas hitting the bounding box also.
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#' @param ... Further arguments passed on the the `ggplot` plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsGeom_point`: A list of arguments for [`ggplot2::geom_point`].
#' * `argsGeom_path`: A list of arguments for [`ggplot2::geom_path`].
#' * `argsGeom_polygon`: A list of arguments for [`ggplot2::geom_polygon`].
#' * `argsGeom_label`: A list of arguments for [`ggplot2::geom_label`].
#'
#' @return The `ggplot` object if `drawPlot = TRUE`; otherwise, a list of `ggplot` components.
#' @export
#' @importFrom rlang .data
#' @import ggplot2
#' @examples
#' library(ggplot2)
#' pts<-matrix(c(1,1), ncol = 2, byrow = TRUE)
#' plotHull2D(pts)
#' pts1<-matrix(c(2,2, 3,3), ncol = 2, byrow = TRUE)
#' plotHull2D(pts1, drawPoints = TRUE)
#' plotHull2D(pts1, drawPoints = TRUE, addRays = TRUE, addText = "coord")
#' plotHull2D(pts1, drawPoints = TRUE, addRays = TRUE, addText = "coord", drawBBoxHull = TRUE)
#' plotHull2D(pts1, drawPoints = TRUE, addRays = TRUE, direction = -1, addText = "coord")
#' pts2<-matrix(c(1,1, 2,2, 0,1), ncol = 2, byrow = TRUE)
#' plotHull2D(pts2, drawPoints = TRUE, addText = "coord")
#' plotHull2D(pts2, drawPoints = TRUE, addRays = TRUE, addText = "coord")
#' plotHull2D(pts2, drawPoints = TRUE, addRays = TRUE, direction = -1, addText = "coord")
#' ## Combine hulls
#' ggplot() +
#' plotHull2D(pts2, drawPoints = TRUE, addText = "coord", drawPlot = FALSE) +
#' plotHull2D(pts1, drawPoints = TRUE, drawPlot = FALSE) +
#' gMOIPTheme() +
#' xlab(expression(x[1])) +
#' ylab(expression(x[2]))
#'
#' # Plotting an LP
#' A <- matrix(c(-3,2,2,4,9,10), ncol = 2, byrow = TRUE)
#' b <- c(3,27,90)
#' obj <- c(7.75, 10)
#' pts3 <- cornerPoints(A, b)
#' plotHull2D(pts3, drawPoints = TRUE, addText = "coord", argsGeom_polygon = list(fill = "red"))
plotHull2D <- function(pts,
drawPoints = FALSE,
drawLines = TRUE,
drawPolygons = TRUE,
addText = FALSE,
addRays = FALSE,
direction = 1,
drawPlot = TRUE,
drawBBoxHull = FALSE,
m = apply(pts,2,min)-5,
M = apply(pts,2,max)+5,
...)
{
args <- list(...)
argsGeom_point <- mergeLists(list(size = 1.5, col = "black"), args$argsGeom_point)
argsGeom_path <- mergeLists(list(size = 0.75, col = "grey40"), args$argsGeom_path)
argsGeom_polygon <- mergeLists(list(fill = "black", alpha = 0.2), args$argsGeom_polygon)
argsGeom_label <- mergeLists(list(nudge_x = 0.2), args$argsGeom_label)
pts <- .checkPts(pts, p = 2)
hull <- convexHull(pts, addRays = addRays, direction = direction, m = m, M = M)
set <- hull$pts
hull <- hull$hull
d <- dimFace(set[,1:2, drop = FALSE])
aES <- aes_string(x = colnames(set)[1], y = colnames(set)[2])
aES1 <- aes_string(x = colnames(set)[1], y = colnames(set)[2], label = 'text')
aES2 <- aes_string(x = colnames(set)[1], y = colnames(set)[2], label = 'addText')
lst <- NULL
if (d > 1) { # a polygon
if (drawPolygons) {
ptsT <- set[hull,]
lst <- c(lst,
do.call(geom_polygon, args = c(list(aES, data = ptsT), argsGeom_polygon)))
}
}
if (d > 0) { # a line or polygon
if (drawLines) {
ptsT <- set[hull,]
for (i in 2:nrow(ptsT)){
if (!(ptsT$pt[i-1] == 0 & ptsT$pt[i] == 0) | drawBBoxHull) {
ptsTT <- ptsT[(i-1):i,]
lst <- c(lst,
do.call(geom_path, args = c(list(aES, data = ptsTT), argsGeom_path)))
}
}
if (!(ptsT$pt[1] == 0 & ptsT$pt[nrow(ptsT)] == 0) | drawBBoxHull) {
ptsTT <- ptsT[c(1,nrow(ptsT)),]
lst <- c(lst,
do.call(geom_path, args = c(list(aES, data = ptsTT), argsGeom_path)))
}
}
}
if (drawPoints | d == 0) {
ptsT <- dplyr::filter(set, .data$pt == 1)
lst <- c(lst,
do.call(geom_point, args = c(list(aES, data = ptsT), argsGeom_point)))
}
ptsT <- dplyr::filter(set, .data$pt == 1)
if (length(addText) > 1) {
if (length(addText) == nrow(ptsT)) {
lst <- c(lst,
do.call(geom_label, args = c(list(aes_string(label = 'addText'), data = ptsT), argsGeom_label)))
}
} else if (length(addText) == 1) {
if (addText == TRUE | addText == "coord") {
ptsT$text <- paste0("(",ptsT[,1],",",ptsT[,2],")")
lst <- c(lst,
do.call(geom_label, args = c(list(aES1, data = ptsT), argsGeom_label)))
} else if (addText == "rownames") {
ptsT$text <- rownames(ptsT)
lst <- c(lst,
do.call(geom_label, args = c(list(aES1, data = ptsT), argsGeom_label)))
} else if (addText == "both") {
text <- paste0("(",ptsT[,1],",",ptsT[,2],")")
ptsT$text <- paste(text,rownames(ptsT))
lst <- c(lst,
do.call(geom_label, args = c(list(aES1, data = ptsT), argsGeom_label)))
} else if (addText != FALSE & length(addText) == nrow(ptsT)) {
lst <- c(lst,
do.call(geom_label, args = c(list(aES2, data = ptsT), argsGeom_label)))
}
}
# if (latex) lst <- c(lst, xlab("$x_1$"), ylab("$x_2$"))
# if (!latex) lst <- c(lst, xlab(expression(x[1])), ylab(expression(x[2])))
if (drawPlot) lst <- ggplot() + lst
return(lst)
}
#' Plot the polytope (bounded convex set) of a linear mathematical program (Ax <= b)
#'
#' This is a wrapper function calling [plotPolytope2D()] (2D graphics) and
#' [plotPolytope3D()] (3D graphics).
#'
#' @param A The constraint matrix.
#' @param b Right hand side.
#' @param obj A vector with objective coefficients.
#' @param type A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' @param nonneg A boolean vector of same length as number of variables. If
#' entry k is TRUE then variable k must be non-negative.
#' @param crit Either max or min (only used if add the iso-profit line)
#' @param faces A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' Useful if e.g. want to show the linear relaxation of an IP.
#' @param plotFaces If \code{True} then plot the faces.
#' @param plotFeasible If \code{True} then plot the feasible points/segments
#' (relevant for IPLP/MILP).
#' @param plotOptimum Show the optimum corner solution point (if alternative solutions
#' only one is shown) and add the iso-profit line.
#' @param latex If \code{True} make latex math labels for TikZ.
#' @param labels If \code{NULL} don't add any labels. If 'n' no labels but show the points. If equal
#' `coord` add coordinates to the points. Otherwise number all points from one.
#' @param ... If 2D, further arguments passed on the the `ggplot` plotting functions. This must be
#' done as lists. Currently the following arguments are supported:
#'
#' * `argsFaces`: A list of arguments for [`plotHull2D`].
#' * `argsFeasible`: A list of arguments for `ggplot2` functions:
#' - `geom_point`: A list of arguments for [`ggplot2::geom_point`].
#' - `geom_line`: A list of arguments for [`ggplot2::geom_line`].
#' * `argsLabels`: A list of arguments for `ggplot2` functions:
#' - `geom_text`: A list of arguments for [`ggplot2::geom_text`].
#' * `argsOptimum`:
#' - `geom_point`: A list of arguments for [`ggplot2::geom_point`].
#' - `geom_abline`: A list of arguments for [`ggplot2::geom_abline`].
#' - `geom_label`: A list of arguments for [`ggplot2::geom_label`].
#' * `argsTheme`: A list of arguments for [`ggplot2::theme`].
#'
#' If 3D further arguments passed on the the RGL plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsAxes3d`: A list of arguments for [`rgl::axes3d`].
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`] to open the RGL window.
#' * `argsTitle3d`: A list of arguments for [`rgl::title3d`].
#' * `argsFaces`: A list of arguments for [`plotHull3D`].
#' * `argsFeasible`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#' - `segments3d`: A list of arguments for [`rgl::segments3d`].
#' - `triangles3d`: A list of arguments for [`rgl::triangles3d`].
#' * `argsLabels`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#' - `text3d`: A list of arguments for [`rgl::text3d`].
#' * `argsOptimum`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#'
#' @note The feasible region defined by the constraints must be bounded (i.e. no extreme rays)
#' otherwise you may see strange results.
#'
#' @return If 2D a `ggplot` object. If 3D a RGL window with the 3D plot.
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @import rgl ggplot2
#' @example inst/examples/ex_polytope.R
plotPolytope <- function(A,
b,
obj = NULL,
type = rep("c", ncol(A)),
nonneg = rep(TRUE, ncol(A)),
crit = "max",
faces = type,
plotFaces = TRUE,
plotFeasible = TRUE,
plotOptimum = FALSE,
latex = FALSE,
labels = NULL, ...)
{
if (length(type)!=ncol(A)) stop("Arg. 'type' must be same length as columns in A!")
if (length(faces)!=ncol(A)) stop("Arg. 'faces' must be same length as columns in A!")
if (!is.null(obj))
if (length(obj)!=ncol(A))
stop("Arg. 'obj' must have the same columns as in A and be a single criterion!")
if (ncol(A)==2) return(plotPolytope2D(A, b, obj, type, nonneg, crit, faces, plotFaces,
plotFeasible, plotOptimum, latex, labels, ...))
if (ncol(A)==3) return(plotPolytope3D(A, b, obj, type, nonneg, crit, faces, plotFaces,
plotFeasible, plotOptimum, latex, labels, ...))
stop("Only 2D or 3D variables supported!")
}
#' Plot the polytope (bounded convex set) of a linear mathematical program
#'
#' @param A The constraint matrix.
#' @param b Right hand side.
#' @param obj A vector with objective coefficients.
#' @param type A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' @param nonneg A boolean vector of same length as number of variables. If
#' entry k is TRUE then variable k must be non-negative.
#' @param crit Either max or min (only used if add the iso-profit line)
#' @param faces A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' Useful if e.g. want to show the linear relaxation of an IP.
#' @param plotFaces If \code{True} then plot the faces.
#' @param plotFeasible If \code{True} then plot the feasible points/segments
#' (relevant for ILP/MILP).
#' @param plotOptimum Show the optimum corner solution point (if alternative solutions
#' only one is shown) and add the iso-profit line.
#' @param latex If \code{True} make latex math labels for TikZ.
#' @param labels If \code{NULL} don't add any labels. If 'n' no labels but show the points. If equal
#' `coord` add coordinates to the points. Otherwise number all points from one.
#' @param ... Further arguments passed on the the `ggplot` plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsFaces`: A list of arguments for [`plotHull2D`].
#' * `argsFeasible`: A list of arguments for `ggplotl2` functions:
#' - `geom_point`: A list of arguments for [`ggplot2::geom_point`].
#' - `geom_line`: A list of arguments for [`ggplot2::geom_line`].
#' * `argsLabels`: A list of arguments for `ggplotl2` functions:
#' - `geom_text`: A list of arguments for [`ggplot2::geom_text`].
#' * `argsOptimum`:
#' - `geom_point`: A list of arguments for [`ggplot2::geom_point`].
#' - `geom_abline`: A list of arguments for [`ggplot2::geom_abline`].
#' - `geom_label`: A list of arguments for [`ggplot2::geom_label`].
#' * `argsTheme`: A list of arguments for [`ggplot2::theme`].
#'
#' @return A `ggplot` object.
#' @note In general use [plotPolytope()] instead of this function. The feasible region defined by the constraints must be bounded otherwise you may see
#' strange results.
#' @seealso [plotPolytope()] for examples.
#' @author Lars Relund \email{lars@@relund.dk}
#' @import ggplot2
plotPolytope2D <-
function(A,
b,
obj = NULL,
type = rep("c", ncol(A)),
nonneg = rep(TRUE, ncol(A)),
crit = "max",
faces = rep("c", ncol(A)),
plotFaces = TRUE,
plotFeasible = TRUE,
plotOptimum = FALSE,
latex = FALSE,
labels = NULL,
...)
{
args <- list(...)
argsFaces <- mergeLists(list(argsGeom_polygon = list(fill = "gray90"), argsGeom_path = list(size = 0.5)), args$argsFaces)
argsFeasible <- list()
argsFeasible$geom_point <- mergeLists(list(), args$argsFeasible$geom_point)
argsFeasible$geom_line <- mergeLists(list(), args$argsFeasible$geom_line)
argsLabels <- list()
argsLabels$geom_text <- mergeLists(list(size=3, color = "gray50", hjust = 1), args$argsLabels$geom_text)
argsOptimum <- list()
argsOptimum$geom_abline <- mergeLists(list(lty="dashed"), args$argsOptimum$geom_abline)
argsOptimum$geom_point <- mergeLists(list(color = "red"), args$argsOptimum$geom_point)
argsOptimum$geom_label <- mergeLists(list(nudge_x = 1, nudge_y = 0.5), args$argsOptimum$geom_label)
argsTheme <- mergeLists(list(), args$argsTheme)
if (!is.null(obj) & (!is.vector(obj) | !length(obj) == ncol(A)))
stop("Arg. obj must be a vector of same length as the number of columns in A.")
# Create solution plot
p <- ggplot() #+ coord_fixed(ratio = 1)
if (latex) p <- p + xlab("$x_1$") + ylab("$x_2$")
if (!latex) p <- p + xlab(expression(x[1])) + ylab(expression(x[2]))
if (plotFaces) {
cPoints <- cornerPoints(A, b, faces, nonneg)
cPoints <- unique(cPoints)
p <- p +
do.call(plotHull2D, args = c(list(cPoints, drawPlot = FALSE), argsFaces))
}
# find feasible points
if (all(type == "c")) {
points <- cornerPoints(A, b, type, nonneg)
} else if (all(type == "i")) {
points <- integerPoints(A, b, nonneg)
} else {
pl <- slices(A, b, type, nonneg)
pl <- lapply(pl, unique)
for (i in 1:length(pl)) pl[[i]] <- cbind(pl[[i]],i)
pl <- lapply(pl, function(x) {
colnames(x) <- c("x1", "x2", "g")
rownames(x) <- NULL
x<-data.frame(x)
})
points <- do.call(rbind, pl)
#points <- points[,1:2]
}
points <- as.data.frame(points)
if (plotFeasible) {
if (all(type == "i")) {
p <- p +
do.call(geom_point, args = c(list(aes_string(x = 'x1', y = 'x2'), data=points), argsFeasible$geom_point))
}
if (length(which(type == "c"))==1) {
p <- p +
do.call(geom_line, args = c(list(aes_string(x = 'x1', y = 'x2', group='g'), data=points), argsFeasible$geom_line)) +
do.call(geom_point, args = c(list(aes_string(x = 'x1', y = 'x2'), data=points), argsFeasible$geom_point))
}
}
if (!is.null(labels)) {
tmp <- points[,1:ncol(A)]
tmp <- as.data.frame(tmp)
if (labels == "coord")
tmp$lbl <- df2String(tmp)
else if (labels == "n")
tmp$lbl <- ""
else
tmp$lbl <- 1:nrow(tmp)
if (length(tmp$lbl)>0) {
#p <- p + do.call(geom_point, args = c(list(aes_string(x = 'x1', y = 'x2'), data=tmp), argsLabels$geom_point))
nudgeS=-(max(tmp$x1)-min(tmp$x1))/100
#if (anyDuplicated(cbind(tmp$x1,tmp$x2), MARGIN = 1) > 0)
p <- p +
do.call(ggrepel::geom_text_repel,
args = c(list(aes_string(x = 'x1', y = 'x2', label = 'lbl'), data=tmp),
argsLabels$geom_text))
}
}
if (plotOptimum) {
if (!is.null(obj)) { # add iso-profit line
tmp <- points
tmp$lbl <- df2String(tmp)
tmp$z <- as.matrix(points[,1:2]) %*% obj
if (crit=="max") i <- which.max(tmp$z)
if (crit=="min") i <- which.min(tmp$z)
if (latex) str <- paste0("$z^* = ", round(tmp$z[i], 2) , "$")
if (!latex) str <- paste0("z* = ", round(tmp$z[i],2) )
if (obj[2]!=0) {
p <- p +
do.call(geom_abline, args = c(list(intercept = tmp$z[i]/obj[2], slope = -obj[1]/obj[2]), argsOptimum$geom_abline))
} else {
p <- p +
do.call(geom_abline, args = c(list(xintercept = tmp$x1[i]), argsOptimum))
}
p <- p +
do.call(geom_point, args = c(list(aes_string(x = 'x1', y = 'x2'), data=tmp[i,]), argsOptimum$geom_point))
p <- p +
do.call(geom_label, args = c(list(aes_string(x = 'x1', y = 'x2', label = 'str'), data = tmp[i,]), argsOptimum$geom_label))
}
}
p <- p +
do.call(gMOIPTheme, args = c(list(), argsTheme))
return(p)
}
#' Plot the polytope (bounded convex set) of a linear mathematical program
#'
#' @param A The constraint matrix.
#' @param b Right hand side.
#' @param obj A vector with objective coefficients.
#' @param type A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' @param nonneg A boolean vector of same length as number of variables. If
#' entry k is TRUE then variable k must be non-negative.
#' @param crit Either max or min (only used if add the iso-profit line)
#' @param faces A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' Useful if e.g. want to show the linear relaxation of an IP.
#' @param plotFaces If \code{True} then plot the faces.
#' @param plotFeasible If \code{True} then plot the feasible points/segments
#' (relevant for ILP/MILP).
#' @param plotOptimum Show the optimum corner solution point (if alternative solutions
#' only one is shown) and add the iso-profit line.
#' @param latex If \code{True} make latex math labels for TikZ.
#' @param labels If \code{NULL} don't add any labels. If 'n' no labels but show the points. If equal
#' `coord` add coordinates to the points. Otherwise number all points from one.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsAxes3d`: A list of arguments for [`rgl::axes3d`].
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`] to open the RGL window.
#' * `argsTitle3d`: A list of arguments for [`rgl::title3d`].
#' * `argsFaces`: A list of arguments for [`plotHull3D`].
#' * `argsFeasible`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#' - `segments3d`: A list of arguments for [`rgl::segments3d`].
#' - `triangles3d`: A list of arguments for [`rgl::triangles3d`].
#' * `argsLabels`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#' - `text3d`: A list of arguments for [`rgl::text3d`].
#' * `argsOptimum`: A list of arguments for RGL functions:
#' - `points3d`: A list of arguments for [`rgl::points3d`].
#'
#' @note In general use [plotPolytope()] instead of this function. The feasible region defined by the constraints must be bounded otherwise you may see
#' strange results.
#' @seealso [plotPolytope()] for examples.
#' @return A RGL window with 3D plot.
#' @author Lars Relund \email{lars@@relund.dk}
#' @import rgl
plotPolytope3D <-
function(A,
b,
obj = NULL,
type = rep("c", ncol(A)),
nonneg = rep(TRUE, ncol(A)),
crit = "max",
faces = rep("c", ncol(A)),
plotFaces = TRUE,
plotFeasible = TRUE,
plotOptimum = FALSE,
latex = FALSE,
labels = NULL, ...)
{
# set plot parameters
args <- list(...)
argsPlot3d <- mergeLists(list(xlab = '', box = FALSE, axes = FALSE), args$argsPlot3d)
argsFaces <- mergeLists(list(drawPoints = FALSE, drawLines = TRUE), args$argsFaces)
argsFeasible <- list()
argsFeasible$points3d <- mergeLists(list(size = 3), args$argsFeasible$points3d)
argsFeasible$segments3d <- mergeLists(list(lwd=1), args$argsFeasible$segments3d)
argsFeasible$triangles3d <- mergeLists(list(col="grey100", alpha=0.4, smooth = FALSE), args$argsFeasible$triangles3d)
argsOptimum <- list()
argsOptimum$points3d <- mergeLists(list(col="red", size = 7), args$argsOptimum$points3d)
argsLabels <- list()
argsLabels$points3d <- mergeLists(list(col="grey50", size = 7), args$argsLabels$points3d)
argsLabels$text3d <- mergeLists(list(), args$argsLabels$text3d)
argsAxes3d <- mergeLists(list(edges=c('x', 'y', 'z')), args$argsAxes3d)
argsTitle3d <- mergeLists(list(xlab = 'x1', ylab = 'x2', zlab = 'x3'), args$argsTitle3d)
do.call(plot3d, args = c(list(x = replicate(2, 1:3), type = 'n'), argsPlot3d))
aspect3d("iso")
if (plotFaces) { # plot faces
vertices <- cornerPoints(A, b, faces, nonneg)
do.call(plotHull3D, args = c(list(pts = vertices), argsFaces))
}
if (plotFeasible) {
if (all(type == "c")) {
# don't plot anything
} else if (all(type == "i")) {
iPoints <- integerPoints(A, b, nonneg)
do.call(points3d, args = c(list(iPoints[,1:3]), argsFeasible$points3d))
} else {
pl <- slices(A, b, type, nonneg)
for (i in 1:length(pl)) {
mat <- pl[[i]]
# print(mat)
if (is.null(mat)) next
if (nrow(mat) == 1) {
do.call(points3d, args = c(list(mat), argsFeasible$points3d))
}
if (nrow(mat) == 2) {
do.call(segments3d, args = c(list(mat), argsFeasible$segments3d))
#cyl <- cylinder3d(mat, radius = 0.015)
#shade3d(cyl, col = "black")
}
if (nrow(mat) == 3) do.call(triangles3d, args = c(list(mat), argsFeasible$triangles3d))
if (nrow(mat) > 3) {
hull <- geometry::convhulln(mat, options = "QJ")
tri <- t(hull)
do.call(triangles3d, args = c(list(mat[tri,1], mat[tri,2], mat[tri,3]), argsFeasible$triangles3d))
}
}
}
}
if (plotOptimum) {
if (is.null(obj)) stop("You need to specify the objective coefficients when using argument plotOption = TRUE.")
vertices <- cornerPoints(A, b, type, nonneg)
val <- vertices[,1:3] %*% as.matrix(obj)
idx <- which.max(val)
val <- vertices[idx,1:3]
do.call(points3d, args = c(list(val[1], val[2], val[3]), argsOptimum$points3d))
}
if (!is.null(labels)) {
if (all(type == "c")) {
points <- cornerPoints(A, b, type, nonneg)
} else if (all(type == "i")) {
points <- integerPoints(A, b, nonneg)
} else {
points <- slices(A, b, type, nonneg, collapse = TRUE)
}
points <- as.data.frame(points)
rownames(points) <- NULL
if ((length(which(type == "c"))<length(type) & length(which(type == "c"))>0) |
(length(which(type == "c"))==length(type))) {
#pch3d(iPoints[,1:3], col="black", cex = 0.1, pch = 16)
do.call(points3d, args = c(list(points[,1:3]), argsLabels$points3d))
}
if (labels=="coord")
points$lbl <- df2String(points)
else if (labels == "n")
points$lbl <- ""
else points$lbl <- 1:nrow(points)
do.call(text3d, args = c(list(points[,1:3], texts = points$lbl), argsLabels$text3d))
}
do.call(axes3d, args = argsAxes3d)
do.call(title3d, args = argsTitle3d)
rgl::highlevel()
}
#' Merge two lists to one
#'
#' @param a First list.
#' @param b Second list.
mergeLists <- function (a,b) {
c(a[setdiff(names(a), names(b))], b)
}
#' Plot the lines of a linear mathematical program (Ax = b)
#'
#' @param A The constraint matrix.
#' @param b Right hand side.
#' @param nonneg A boolean vector of same length as number of variables. If
#' entry k is TRUE then variable k must be non-negative and the line is plotted too.
#' @param latex If \code{True} make latex math labels for TikZ.
#' @param ... Further arguments passed on the the `ggplot` plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsTheme`: A list of arguments for [`ggplot2::theme`].
#'
#' @return A `ggplot` object.
#' @note In general you will properly use [plotPolytope()] instead of this function.
#' @seealso [plotPolytope()].
#' @author Lars Relund \email{lars@@relund.dk}
#' @import ggplot2
plotLines2D <-
function(A,
b,
nonneg = rep(TRUE, ncol(A)),
latex = FALSE,
...)
{
args <- list(...)
argsTheme <- mergeLists(list(), args$argsTheme)
# Create solution plot
p <- ggplot() #+ coord_fixed(ratio = 1)
if (latex) p <- p + xlab("$x_1$") + ylab("$x_2$")
if (!latex) p <- p + xlab(expression(x[1])) + ylab(expression(x[2]))
for (i in 1:nrow(A)) {
if (A[i,1] == 0) {
p <- p + geom_hline(yintercept = b[i])
next
}
if (A[i,2] == 0) {
p <- p + geom_vline(xintercept = b[i])
next
}
p <- p + geom_abline(slope = -A[i,1]/A[i,2], intercept = b[i]/A[i,2])
}
if (nonneg[1]) p <- p + geom_vline(xintercept = 0)
if (nonneg[2]) p <- p + geom_hline(yintercept = 0)
l <- cornerPoints(A, b, nonneg = nonneg)
p <- p + xlim(min(l[,1]), max(l[,1])) +
ylim(min(l[,2]), max(l[,2])) +
do.call(gMOIPTheme, args = c(list(), argsTheme))
return(p)
}
#' Create a plot of the criterion space of a bi-objective problem
#'
#' @param A The constraint matrix.
#' @param b Right hand side.
#' @param obj A p x n matrix(one row for each criterion).
#' @param type A character vector of same length as number of variables. If
#' entry k is 'i' variable \eqn{k} must be integer and if 'c' continuous.
#' @param nonneg A boolean vector of same length as number of variables. If
#' entry k is TRUE then variable k must be non-negative.
#' @param crit Either max or min (only used if add the iso-profit line).
#' @param addTriangles Add search triangles defined by the non-dominated extreme
#' points.
#' @param addHull Add the convex hull and the rays.
#' @param plotFeasible If \code{True} then plot the criterion points/slices.
#' @param latex If true make latex math labels for TikZ.
#' @param labels If \code{NULL} don't add any labels. If 'n' no labels but show the points. If equal
#' `coord` add coordinates to the points. Otherwise number all points from one.
#'
#' @note Currently only points are checked for dominance. That is, for MILP
#' models some nondominated points may in fact be dominated by a segment.
#' @return The `ggplot` object.
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @import ggplot2
#' @example inst/examples/ex_criterion.R
plotCriterion2D <- function(A,
b,
obj,
type = rep("c", ncol(A)),
nonneg = rep(TRUE, ncol(A)),
crit = "max",
addTriangles = FALSE,
addHull = TRUE,
plotFeasible = TRUE,
latex = FALSE,
labels = NULL)
{
# Set Custom theme
myTheme <- theme_bw()
myTheme <-
myTheme + theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
#axis.line = element_blank(),
axis.line = element_line(
colour = "black",
size = 0.5,
arrow = arrow(length = unit(0.3, "cm"))
),
#axis.ticks = element_blank()
#axis.text.x = element_text(margin = margin(r = 30))
# axis.ticks.length = unit(0.5,"mm"),
#aspect.ratio=4/3,
legend.position = "none"
)
# First find all relevant points
if (all(type == "i")) {
points <- integerPoints(A, b, nonneg)
} else if (all(type == "c")) {
points <- cornerPoints(A, b, type, nonneg)
} else {
points <- slices(A, b, type, nonneg, collapse = TRUE)
}
points <- criterionPoints(points, obj, crit, labels)
if (all(type == "c")) { # if cont then no non-ext
points$nd[points$nd & !points$se] <- FALSE
points$se[points$nd & !points$se] <- FALSE
points$sne[points$nd & !points$se] <- FALSE
}
#dat <<- points # hack to get data as data frame
if (crit=="max") points <- points[order(-points$z2,-points$z1),]
if (crit=="min") points <- points[order(points$z2,points$z1),]
# Initialize plot
p <- ggplot(points, aes_q(x = quote(z1), y = quote(z2), col = "grey10") )
if (latex) p <- p + xlab("$z_1$") + ylab("$z_2$")
if (!latex) p <- p + xlab(expression(z[1])) + ylab(expression(z[2]))
# Add hull plus rays
if (addHull) {
tmp<-points[points$se & !duplicated(cbind(points$z1,points$z2), MARGIN = 1),]
delta <- max( (max(points$z1)-min(points$z1))/10, (max(points$z2)-min(points$z2))/10 )
if (crit=="max") {
tmp<-rbind(tmp[1:2,],tmp,tmp[1,]) # add rows
tmp$z1[1] <- min(points$z1) - delta
tmp$z2[1] <- min(points$z2) - delta
tmp$z1[2] <- min(points$z1) - delta
tmp$z2[2] <- max(points$z2)
tmp$z1[length(tmp$z1)] <- max(points$z1)
tmp$z2[length(tmp$z1)] <- min(points$z2)- delta
}
if (crit=="min") {
tmp<-rbind(tmp[1,],tmp,tmp[1:2,]) # add rows
tmp$z1[1] <- max(points$z1) + delta
tmp$z2[1] <- min(points$z2)
tmp$z1[length(tmp$z1)-1] <- min(points$z1)
tmp$z2[length(tmp$z1)-1] <- max(points$z2) + delta
tmp$z1[length(tmp$z1)] <- max(points$z1) + delta
tmp$z2[length(tmp$z1)] <- max(points$z2) + delta
}
p <- p + geom_polygon(fill="gray95", col = NA, data=tmp)
}
if (plotFeasible) {
# plot feasible areas
if (all(type == "c")) {
idx <- grDevices::chull(points[,c("z1","z2")])
p <- p + geom_polygon(data = points[idx,], aes_string(x = 'z1', y = 'z2'),
fill=NA, size = 0.5, linetype = 1, col="grey80", alpha=0.6)
p <- p + geom_point(aes_string(colour = 'nd', shape = 'se'), data = points) +
scale_colour_grey(start = 0.6, end = 0)
} else if (all(type == "i")) {
#iPoints <- integerPoints(A, b, nonneg)
#iPoints <- criterionPoints(iPoints, obj, crit, labels)
p <- p + geom_point(aes_string(colour = 'nd', shape = 'se'), data = points) +
# #coord_fixed(ratio = 1) +
scale_colour_grey(start = 0.6, end = 0)
} else {
pl <- slices(A, b, type, nonneg)
pl <- lapply(pl, function(x) {
x <- x %*% t(obj)
if (is.vector(x)) x <- matrix(x, nrow=1)
x[,1:2, drop = FALSE]
})
idx <- lapply(pl, grDevices::chull)
#idx <- lapply(idx, function(x) c(x, x[1]))
pl <- mapply(function(x,y){
x[y,,drop = FALSE]
},
x = pl, y = idx, SIMPLIFY = FALSE)
for (i in 1:length(pl)) pl[[i]] <- cbind(pl[[i]],i)
# pl <- lapply(pl, function(x) {
# colnames(x) <- c("z1", "z2", "g")
# rownames(x) <- NULL
# x<-data.frame(x)
# })
pl <- do.call(rbind, pl)
pl <- as.data.frame(pl)
colnames(pl) <- c("z1", "z2", "g")
rownames(pl) <- NULL
p <- p + geom_polygon(
data = pl,
alpha = 0.6, fill = "grey",
aes_string(
x = "z1",
y = "z2",
group = "g"
)
) + scale_fill_identity()
points$co <- points$se | points$sne
p <- p + geom_point(aes_string(colour = 'co', shape = 'se'), data = points) +
scale_colour_grey(start = 0.8, end = 0)
}
}
# Add triangles
if (addTriangles) {
tmp<-points[points$se | points$sne,]
if (length(tmp$z1)>1) { # triangles
for (r in 1:(dim(tmp)[1] - 1)) {
p <- p +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r],
xend = tmp$z1[r + 1],
yend = tmp$z2[r + 1],
colour = "gray"
) +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r],
xend = tmp$z1[r],
yend = tmp$z2[r + 1],
colour = "gray"
) +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r + 1],
xend = tmp$z1[r + 1],
yend = tmp$z2[r + 1],
colour = "gray"
)
}
}
}
nudgeC=-(max(points$z1)-min(points$z1))/100
if (!is.null(labels) & anyDuplicated(round(cbind(points$z1,points$z2),10), MARGIN = 1) > 0)
p <- p + ggrepel::geom_text_repel(aes_string(label = 'lbl'), size=3,
colour = "gray50", data=points)
if (!is.null(labels) & anyDuplicated(round(cbind(points$z1,points$z2),10), MARGIN = 1) == 0)
p <- p + geom_text(aes_string(label = 'lbl'), nudge_x = nudgeC, nudge_y = nudgeC,
hjust=1, size=3, colour = "gray50", data=points)
p <- p + myTheme
return(p)
}
#' Help function to save the view angle for the RGL 3D plot
#'
#' @param fname The file name of the view.
#' @param overwrite Overwrite existing file.
#' @param print Print the view so can be copied to R code (no file is saved).
#'
#' @note Only save if the file name don't exists.
#' @return NULL (invisible).
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @examples
#' \donttest{
#' view <- matrix( c(-0.412063330411911, -0.228006735444069, 0.882166087627411, 0,
#' 0.910147845745087, -0.0574885793030262, 0.410274744033813, 0, -0.042830865830183,
#' 0.97196090221405, 0.231208890676498, 0, 0, 0, 0, 1), nc = 4)
#'
#' loadView(v = view)
#' A <- matrix( c(3, 2, 5, 2, 1, 1, 1, 1, 3, 5, 2, 4), nc = 3, byrow = TRUE)
#' b <- c(55, 26, 30, 57)
#' obj <- c(20, 10, 15)
#' plotPolytope(A, b, plotOptimum = TRUE, obj = obj, labels = "coord")
#'
#' # Try to modify the angle in the RGL window
#' saveView(print = TRUE) # get the view angle to insert into R code
#' }
saveView <- function(fname = "view.RData", overwrite = FALSE, print = FALSE) {
if (print) {
view <- rgl::par3d()$userMatrix
cat(paste0("view <- matrix( c(", paste0(view, collapse = ", "), "), nc = 4)"))
} else if (!file.exists(fname) | overwrite) {
view <- rgl::par3d()$userMatrix
save(view, file = fname)
message(paste0("RGL view saved to RData file ", fname, "."))
}
return(invisible(NULL))
}
#' Help function to load the view angle for the RGL 3D plot from a file or matrix
#'
#' @param fname The file name of the view.
#' @param v The view matrix.
#' @param clear Call [rgl::clear3d()].
#' @param close Call [rgl::close3d()].
#' @param zoom Zoom level.
#' @param ... Additional parameters passed to [rgl::view3d()].
#'
#' @return NULL
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @examples
#' \donttest{
#' view <- matrix( c(-0.412063330411911, -0.228006735444069, 0.882166087627411, 0,
#' 0.910147845745087, -0.0574885793030262, 0.410274744033813, 0, -0.042830865830183,
#' 0.97196090221405, 0.231208890676498, 0, 0, 0, 0, 1), nc = 4)
#'
#' loadView(v = view)
#' A <- matrix( c(3, 2, 5, 2, 1, 1, 1, 1, 3, 5, 2, 4), nc = 3, byrow = TRUE)
#' b <- c(55, 26, 30, 57)
#' obj <- c(20, 10, 15)
#' plotPolytope(A, b, plotOptimum = TRUE, obj = obj, labels = "coord")
#'
#' # Try to modify the angle in the RGL window
#' saveView(print = TRUE) # get the view angle to insert into R code
#' }
loadView <- function(fname = "view.RData", v = NULL, clear = TRUE, close = FALSE, zoom = 1, ...) {
if (clear) rgl::clear3d()
if (close) rgl::close3d()
if (!is.null(v)) {
rgl::view3d(userMatrix = v, zoom = zoom, ...)
} else {
if (file.exists(fname)) {
load(fname)
rgl::view3d(userMatrix = v)
} else {
warning(paste0("Can't load view in file ", fname, "!"))
}
}
}
#' Create a plot of a discrete non-dominated set.
#'
#' @param points Data frame with non-dominated points.
#' @param crit Either max or min (only used if add the iso-profit line). A vector is currently not
#' supported.
#' @param addTriangles Add search triangles defined by the non-dominated extreme points.
#' @param addHull Add the convex hull and the rays.
#' @param latex If true make latex math labels for TikZ.
#' @param labels If \code{NULL} don't add any labels. If 'n' no labels but show the points. If equal
#' `coord` add coordinates to the points. Otherwise number all points from one.
#'
#' @note Currently only points are checked for dominance. That is, for MILP models some
#' nondominated points may in fact be dominated by a segment.
#' @return The `ggplot` object.
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @import ggplot2
#' @examples
#' dat <- data.frame(z1=c(12,14,16,18,18,18,14,15,15), z2=c(18,16,12,4,2,6,14,14,16))
#' points <- addNDSet(dat, crit = "min", keepDom = TRUE)
#' plotNDSet2D(points, crit = "min", addTriangles = TRUE)
#' plotNDSet2D(points, crit = "min", addTriangles = FALSE)
#' plotNDSet2D(points, crit = "min", addTriangles = TRUE, addHull = FALSE)
#' points <- addNDSet(dat, crit = "max", keepDom = TRUE)
#' plotNDSet2D(points, crit = "max", addTriangles = TRUE)
#' plotNDSet2D(points, crit = "max", addHull = FALSE)
plotNDSet2D <- function(points,
crit,
addTriangles = FALSE,
addHull = TRUE,
latex = FALSE,
labels = NULL)
{
# Set Custom theme
myTheme <- theme_bw()
myTheme <-
myTheme + theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
#axis.line = element_blank(),
axis.line = element_line(
colour = "black",
size = 0.5,
arrow = arrow(length = unit(0.3, "cm"))
),
#axis.ticks = element_blank()
#axis.text.x = element_text(margin = margin(r = 30))
# axis.ticks.length = unit(0.5,"mm"),
#aspect.ratio=4/3,
legend.position = "none"
)
if (crit=="max") points <- points[order(-points$z2,-points$z1),]
if (crit=="min") points <- points[order(points$z2,points$z1),]
# Initialize plot
p <- ggplot(points, aes_q(x = quote(z1), y = quote(z2)) )
if (latex) p <- p + xlab("$z_1$") + ylab("$z_2$")
if (!latex) p <- p + xlab(expression(z[1])) + ylab(expression(z[2]))
# Add hull plus rays
if (addHull) {
di <- .mToDirection(crit, 2)
p <- p +
plotHull2D(points[points$nd, c("z1", "z2")],
addRays = TRUE, drawLines = TRUE, drawBBoxHull = FALSE,
direction = di, drawPlot = FALSE)
}
# Add triangles
if (addTriangles) {
tmp<-points[points$se | points$sne,]
if (length(tmp$z1)>1) { # triangles
for (r in 1:(dim(tmp)[1] - 1)) {
p <- p +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r],
xend = tmp$z1[r + 1],
yend = tmp$z2[r + 1],
colour = "gray"
) +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r],
xend = tmp$z1[r],
yend = tmp$z2[r + 1],
colour = "gray"
) +
geom_segment(
x = tmp$z1[r],
y = tmp$z2[r + 1],
xend = tmp$z1[r + 1],
yend = tmp$z2[r + 1],
colour = "gray"
)
}
}
}
p <- p + geom_point(aes_string(colour = 'nd', shape = 'se'), data = points) +
# #coord_fixed(ratio = 1) +
scale_colour_grey(start = 0.6, end = 0)
nudgeC=-(max(points$z1)-min(points$z1))/100
if (!is.null(labels) &
anyDuplicated(round(cbind(points$z1, points$z2), 10), MARGIN = 1) > 0)
p <- p + ggrepel::geom_text_repel(
aes_string(label = 'lbl'),
size = 3,
colour = "gray50",
data = points
)
if (!is.null(labels) &
anyDuplicated(round(cbind(points$z1, points$z2), 10), MARGIN = 1) == 0)
p <-
p + geom_text(
aes_string(label = 'lbl'),
nudge_x = nudgeC,
nudge_y = nudgeC,
hjust = 1,
size = 3,
colour = "gray50",
data = points
)
p <- p + myTheme
return(p)
}
#' Plot a rectangle defined by two corner points.
#'
#' The rectangle is defined by \{x|a <= x <= b\} where a is the minimum values and
#' b is the maximum values.
#'
#' @param a A vector of length 3.
#' @param b A vector of length 3.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`].
#' * `argsSegments3d`: A list of arguments for [`rgl::segments3d`][rgl::points3d()].
#' * `argsPolygon3d`: A list of arguments for [`rgl::polygon3d`].
#' * `argsShade3d`: A list of arguments for [`rgl::shade3d`][rgl::mesh3d()].
#'
#' @return Object ids (invisible).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D()
#' plotRectangle3D(c(0,0,0), c(1,1,1))
#' plotRectangle3D(c(1,1,1), c(4,4,3), drawPoints = TRUE, drawLines = FALSE,
#' argsPlot3d = list(size=2, type="s", alpha=0.3))
#' ids <- plotRectangle3D(c(2,2,2), c(3,3,2.5), argsPolygon3d = list(alpha = 1) )
#' finalize3D()
#' # pop3d(id = ids) remove last object
#' }
plotRectangle3D <- function(a, b, ...) {
args <- list(...)
argsSegments3d <- mergeLists(list(), args$argsSegments3d)
argsPlot3d <- mergeLists(list(), args$argsPlot3d)
argsShade3d <- mergeLists(list(), args$argsShade3d)
argsPolygon3d <- mergeLists(list(), args$argsPolygon3d)
if (any(a == b)) stop("Define different min and max values!")
if (is.matrix(a)) a <- as.vector(a)
if (is.matrix(b)) b <- as.vector(b)
x<-matrix(c(a,b), ncol=3, byrow = TRUE)
x<-expand.grid(x=c(x[1,1],x[2,1]), y=c(x[1,2],x[2,2]), z=c(x[1,3],x[2,3]))
argsP <- list()
if (!is.null(args$drawPoints)) argsP$drawPoints <- args$drawPoints
if (!is.null(args$drawLines)) argsP$drawLines <- args$drawLines
lst <- do.call(plotHull3D, args = c(list(x), list(argsSegments3d = argsSegments3d,
argsPlot3d = argsPlot3d, argsShade3d = argsShade3d,
argsPolygon3d = argsPolygon3d), argsP))
return(invisible(lst$ids))
}
#' Plot a polygon.
#'
#' @param pts Vertices.
#' @param useShade Plot shade of the polygon.
#' @param useLines Plot lines inside the polygon.
#' @param usePoints Plot point shapes inside the polygon.
#' @param useFrame Plot a frame around the polygon.
#' @param ... Further arguments passed on the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsShade`: A list of arguments for [`rgl::polygon3d`] (n > 4 vertices),
#' [rgl::triangles3d()] (n = 3 vertices) and [rgl::quads3d()] (n = 4 vertices)
#' if `useShade = TRUE`.
#' * `argsFrame`: A list of arguments for [`rgl::lines3d`] if `useFrame = TRUE`.
#' * `argsPoints`: A list of arguments for [`rgl::shade3d`] if `usePoints = TRUE`. It is important
#' to give a texture using `texture`. A texture can be set using [getTexture()].
#' * `argsLines`: A list of arguments for [rgl::persp3d()] when `useLines = TRUE`. Moreover, the list
#' may contain `lines`: number of lines.
#'
#' @return Object ids (invisible).
#' @importFrom stats lm predict
#' @export
#'
#' @examples
#' \donttest{
#' pts <- data.frame(x = c(1,0,0,0.4), y = c(0,1,0,0.3), z = c(0,0,1,0.3))
#' pts <- data.frame(x = c(1,0,0), y = c(0,1,0), z = c(0,0,1))
#'
#' ini3D()
#' plotPolygon3D(pts)
#' finalize3D()
#'
#' ini3D()
#' plotPolygon3D(pts, argsShade = list(color = "red", alpha = 1))
#' finalize3D()
#'
#' ini3D()
#' plotPolygon3D(pts, useFrame = TRUE, argsShade = list(color = "red", alpha = 0.5),
#' argsFrame = list(color = "green"))
#' finalize3D()
#'
#' ini3D()
#' plotPolygon3D(pts, useFrame = TRUE, useLines = TRUE, useShade = TRUE,
#' argsShade = list(color = "red", alpha = 0.2),
#' argsLines = list(color = "blue"))
#' finalize3D()
#'
#' ini3D()
#' ids <- plotPolygon3D(pts, usePoints = TRUE, useFrame = TRUE,
#' argsPoints = list(texture = getTexture(pch = 16, cex = 20)))
#' finalize3D()
#' # pop3d(id = ids) # remove object again
#'
#' # In general you have to finetune size and numbers when you use textures
#' # Different pch
#' for (i in 0:3) {
#' fname <- getTexture(pch = 15+i, cex = 30)
#' ini3D(TRUE)
#' plotPolygon3D(pts, usePoints = TRUE, argsPoints = list(texture = fname))
#' finalize3D()
#' }
#'
#' # Size of pch
#' for (i in 1:4) {
#' fname <- getTexture(pch = 15+i, cex = 10 * i)
#' ini3D(TRUE)
#' plotPolygon3D(pts, usePoints = TRUE, argsPoints = list(texture = fname))
#' finalize3D()
#' }
#'
#' # Number of pch
#' fname <- getTexture(pch = 16, cex = 20)
#' for (i in 1:4) {
#' ini3D(TRUE)
#' plotPolygon3D(pts, usePoints = TRUE,
#' argsPoints = list(texture = fname, texcoords = rbind(pts$x, pts$y, pts$z)*5*i))
#' finalize3D()
#' }
#' }
plotPolygon3D <- function(pts, useShade = TRUE, useLines = FALSE, usePoints = FALSE,
useFrame = TRUE, ...) {
args <- list(...)
# args <- list(argsPoints=argsPoints)
argsShade <- mergeLists(list(color = "black", col = "grey40",
lwd = 2, alpha = 0.2, fill = TRUE),
args$argsShade)
argsFrame <- mergeLists(list(color = "black", lwd = 1, alpha = 0.8), args$argsFrame)
argsPoints <- mergeLists(list(color = "grey", specular = "black", alpha = 0.5,
texcoords = rbind(pts[,1], pts[,2], pts[,3])*25
),
args$argsPoints)
argsLines <- mergeLists(list(color = "black", lwd = 1, alpha = 0.4, lines = 50, back = 'lines',
front = 'lines', lit = FALSE),
args$argsLines)
ids <- NULL
pts <- .checkPts(pts, p = 3, asDF = TRUE)
if (dimFace(pts) != 2) stop("Vertices don't define a polygon!")
if (usePoints) useShade = FALSE
if (useShade) {
if (nrow(pts) > 3) ids <- c(ids, do.call(rgl::polygon3d, args = c(list(pts), argsShade)))
if (nrow(pts) == 3) ids <- c(ids, do.call(rgl::triangles3d, args = c(list(pts), argsShade)))
}
if (useFrame) {
n <- length(pts[,1])
nas <- n+1
prev <- 0L
loop <- integer()
for (i in seq_along(nas)) {
loop <- c(loop, if (i > 1) NA, (prev + 1L):(nas[i] - 1L), prev + 1L)
prev <- nas[i]
}
res <- cbind(pts[loop,1], pts[loop,2], pts[loop,3])
ids <- c(ids, do.call(rgl::lines3d, args = c(list(res), argsFrame)))
}
if (usePoints) {
if (nrow(pts) > 3)
poly <- do.call(rgl::polygon3d, args = c(list(pts, plot = FALSE), argsPoints))
if (nrow(pts) == 3)
poly <- mesh3d(vertices = rbind(pts[,1], pts[,2], pts[,3], 1), triangles = 1:3)
# poly$texcoords <- argsPoints$texcoords
ids <- c(ids, do.call(rgl::shade3d, args = c(list(poly), argsPoints)))
}
if (useLines) {
if (!rgl::cur3d()) stop("Option useLines need an open RGL window!")
limits <- rgl::par3d()$bbox
m <- c(limits[1], limits[3], limits[5])
M <- c(limits[2], limits[4], limits[6])
if (all(m == M)) stop("The RGL window's bounding box is not valid!")
# do the mesh
x <- seq(m[1], M[1], length.out = argsLines$lines)
y <- seq(m[2], M[2], length.out = argsLines$lines)
xy <- expand.grid(x, y)
colnames(pts) <- c("x", "y", "z")
xy[sp::point.in.polygon(xy[,1], xy[,2], pts$x, pts$y) <= 0,] <- NA
z <- predict(lm(z ~ x + y, data = pts), newdata = data.frame(x=xy[,1], y=xy[,2]))
ids <- c(ids, do.call(rgl::persp3d, args = c(list(x, y, z, add = TRUE), argsLines)))
}
print(rgl::highlevel())
return(invisible(ids))
}
#' Save a point symbol as a temporary file.
#'
#' @param pch Point number/symbol.
#' @param cex Point size
#' @param ... Further arguments passed to `plot`.
#'
#' @return The file name.
#' @export
#'
#' @examples
#' \donttest{
#' # Pch shapes
#' generateRPointShapes<-function(){
#' oldPar<-par()
#' par(font=2, mar=c(0.5,0,0,0))
#' y=rev(c(rep(1,6),rep(2,5), rep(3,5), rep(4,5), rep(5,5)))
#' x=c(rep(1:5,5),6)
#' plot(x, y, pch = 0:25, cex=1.5, ylim=c(1,5.5), xlim=c(1,6.5),
#' axes=FALSE, xlab="", ylab="", bg="blue")
#' text(x, y, labels=0:25, pos=3)
#' par(mar=oldPar$mar,font=oldPar$font )
#' }
#' generateRPointShapes()
#'
#' getTexture()
#' }
getTexture <- function(pch = 16, cex = 10, ...) {
filename <- tempfile(fileext = ".png")
grDevices::png(filename = filename)
graphics::plot(1, 1, pch = pch, cex = cex, axes=FALSE, xlab="", ylab="", ...)
grDevices::dev.off()
return(filename)
}
#' Plot a cone defined by a point in 3D.
#'
#' The cones are defined as the point plus R3+.
#'
#' @param pts A matrix with a point in each row.
#' @param drawPoint Draw the points defining the cone.
#' @param drawLines Draw lines of the cone.
#' @param drawPolygons Draw polygons of the cone.
#' @param rectangle Draw the cone as a rectangle.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of `pts`
#' plus a value greater than on equal zero (minimize objective $i$). If negative, consider the
#' i'th column of `pts` minus a value greater than on equal zero (maximize objective $i$).
#' @param useRGLBBox Use the RGL bounding box as ray limits for the cone.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`].
#' * `argsSegments3d`: A list of arguments for [`rgl::segments3d`][rgl::points3d()].
#' * `argsPolygon3d`: A list of arguments for [`rgl::polygon3d`].
#'
#' @return Object ids (invisible).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D(argsPlot3d = list(xlim = c(0,6), ylim = c(0,6), zlim = c(0,6)))
#' plotCones3D(c(4,4,4), drawLines = FALSE, drawPoint = TRUE,
#' argsPlot3d = list(col = "red", size = 10),
#' argsPolygon3d = list(alpha = 1), rectangle = TRUE)
#' plotCones3D(c(1,1,1), rectangle = FALSE)
#' plotCones3D(matrix(c(3,3,3,2,2,2), ncol = 3, byrow = TRUE))
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(xlim = c(0,6), ylim = c(0,6), zlim = c(0,6)))
#' plotCones3D(c(4,4,4), direction = 1)
#' plotCones3D(c(2,2,2), direction = -1)
#' plotCones3D(c(4,2,2), direction = c(1,-1,-1))
#' ids <- plotCones3D(c(2,2,4), direction = c(-1,-1,1))
#' finalize3D()
#' # pop3d(id = ids) # remove last cone
#' }
plotCones3D <-
function(pts,
drawPoint = TRUE,
drawLines = TRUE,
drawPolygons = TRUE,
direction = 1,
rectangle = FALSE,
useRGLBBox = TRUE,
...) {
args <- list(...)
argsPlot3d <- mergeLists(list(), args$argsPlot3d)
argsSegments3d <- mergeLists(list(lwd = 1, col = "grey40"), args$argsSegments3d)
argsPolygon3d <- mergeLists(list(), args$argsPolygon3d)
pts <- .checkPts(pts, p = 3)
for (i in 1:dim(pts)[1]) {
pt <- as.vector(pts[i, ])
lst <- plotHull3D(pt, drawPoints = drawPoint,
drawLines = drawLines, drawPolygons = drawPolygons,
addRays = TRUE, direction = direction, drawBBoxHull = rectangle,
useRGLBBox = useRGLBBox,
argsPlot3d = argsPlot3d, argsSegments3d = argsSegments3d,
argsPolygon3d = argsPolygon3d
)
}
return(invisible(lst$ids))
}
#' Plot a cone defined by a point in 2D.
#'
#' The cones are defined as the point plus/minus rays of R2.
#'
#' @param pts A matrix with a point in each row.
#' @param drawPoint Draw the points defining the cone.
#' @param drawLines Draw lines of the cone.
#' @param drawPolygons Draw polygons of the cone.
#' @param rectangle Draw the cone as a rectangle.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of `pts`
#' plus a value greater than on equal zero (minimize objective $i$). If negative, consider the
#' i'th column of `pts` minus a value greater than on equal zero (maximize objective $i$).
#' @param drawPlot Draw the `ggplot`. Set to FALSE if you want to combine hulls in a single plot.
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#' @param ... Further arguments passed to [plotHull2D]
#'
#' @return A `ggplot` object
#' @export
#' @import ggplot2
#' @examples
#' library(ggplot2)
#' plotCones2D(c(4,4), drawLines = FALSE, drawPoint = TRUE,
#' argsGeom_point = list(col = "red", size = 10),
#' argsGeom_polygon = list(alpha = 0.5), rectangle = TRUE)
#' plotCones2D(c(1,1), rectangle = FALSE)
#' plotCones2D(matrix(c(3,3,2,2), ncol = 2, byrow = TRUE))
#'
#' ## The Danish flag
#' lst <- list(argsGeom_polygon = list(alpha = 0.85, fill = "red"),
#' drawPlot = FALSE, drawPoint = FALSE, drawLines = FALSE)
#' p1 <- do.call(plotCones2D, args = c(list(c(2,4), direction = 1), lst))
#' p2 <- do.call(plotCones2D, args = c(list(c(1,2), direction = -1), lst))
#' p3 <- do.call(plotCones2D, args = c(list(c(2,2), direction = c(1,-1)), lst))
#' p4 <- do.call(plotCones2D, args = c(list(c(1,4), direction = c(-1,1)), lst))
#' ggplot() + p1 + p2 + p3 + p4 + theme_void()
plotCones2D <-
function(pts,
drawPoint = TRUE,
drawLines = TRUE,
drawPolygons = TRUE,
direction = 1,
rectangle = FALSE,
drawPlot = TRUE,
m = apply(pts,2,min)-5,
M = apply(pts,2,max)+5,
...) {
pts <- .checkPts(pts, p = 2)
lst <- NULL
for (i in 1:dim(pts)[1]) {
pt <- as.vector(pts[i, ])
lst <- c(lst,
plotHull2D(pt, drawPoints = drawPoint,
drawLines = drawLines, drawPolygons = drawPolygons,
addRays = TRUE, direction = direction, drawPlot = FALSE,
drawBBoxHull = rectangle, m = m, M = M, ...))
}
if (drawPlot) lst <- ggplot() + lst
return(lst)
}
#' Plot the convex hull of a set of points in 3D.
#'
#' @param pts A matrix with a point in each row.
#' @param drawPoints Draw the points.
#' @param drawLines Draw lines of the facets.
#' @param drawPolygons Fill the facets.
#' @param addText Add text to the points. Currently `coord` (coordinates), `rownames` (rownames)
#' and `both` supported or a vector with text.
#' @param addRays Add the ray defined by `direction`.
#' @param useRGLBBox Use the RGL bounding box when add rays.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of `pts`
#' plus a value greater than on equal zero (minimize objective $i$). If negative, consider the
#' i'th column of `pts` minus a value greater than on equal zero (maximize objective $i$).
#' @param drawBBoxHull If `addRays` then draw the hull areas hitting the bounding box also.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`].
#' * `argsSegments3d`: A list of arguments for [`rgl::segments3d`][rgl::points3d()].
#' * `argsPolygon3d`: A list of arguments for [`rgl::polygon3d`].
#' * `argsShade3d`: A list of arguments for [`rgl::shade3d`][rgl::mesh3d()].
#' * `argsText3d`: A list of arguments for [`rgl::text3d`].
#'
#' @return A list with hull, `pts` classified and object ids (invisible).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D()
#' pts<-matrix(c(0,0,0), ncol = 3, byrow = TRUE)
#' plotHull3D(pts) # a point
#' pts<-matrix(c(1,1,1,2,2,2,3,3,3), ncol = 3, byrow = TRUE)
#' plotHull3D(pts, drawPoints = TRUE) # a line
#' pts<-matrix(c(1,0,0,1,1,1,1,2,2,3,1,1,3,3,3), ncol = 3, byrow = TRUE)
#' plotHull3D(pts, drawLines = FALSE, argsPolygon3d = list(alpha=0.6)) # a polygon
#' pts<-matrix(c(5,5,5,10,10,5,10,5,5,5,5,10), ncol = 3, byrow = TRUE)
#' lst <- plotHull3D(pts, argsPolygon3d = list(alpha=0.9), argsSegments3d = list(color="red"))
#' finalize3D()
#' # pop3d(id = lst$ids) # remove last hull
#'
#' ## Using addRays
#' pts <- data.frame(x = c(1,3), y = c(1,3), z = c(1,3))
#' ini3D(argsPlot3d = list(xlim = c(0,max(pts$x)+10),
#' ylim = c(0,max(pts$y)+10),
#' zlim = c(0,max(pts$z)+10)))
#' plotHull3D(pts, drawPoints = TRUE, addRays = TRUE, , drawBBoxHull = FALSE)
#' plotHull3D(c(4,4,4), drawPoints = TRUE, addRays = TRUE)
#' finalize3D()
#'
#' pts <- data.frame(x = c(4,2.5,1), y = c(1,2.5,4), z = c(1,2.5,4))
#' ini3D(argsPlot3d = list(xlim = c(0,max(pts$x)+10),
#' ylim = c(0,max(pts$y)+10),
#' zlim = c(0,max(pts$z)+10)))
#' plotHull3D(pts, drawPoints = TRUE, addRays = TRUE)
#' finalize3D()
#'
#' pts <- matrix(c(
#' 0, 4, 8,
#' 0, 8, 4,
#' 8, 4, 0,
#' 4, 8, 0,
#' 4, 0, 8,
#' 8, 0, 4,
#' 4, 4, 4,
#' 6, 6, 6
#' ), ncol = 3, byrow = TRUE)
#' ini3D(FALSE, argsPlot3d = list(xlim = c(min(pts[,1])-2,max(pts[,1])+10),
#' ylim = c(min(pts[,2])-2,max(pts[,2])+10),
#' zlim = c(min(pts[,3])-2,max(pts[,3])+10)))
#' plotHull3D(pts, drawPoints = TRUE, addText = "coord")
#' plotHull3D(pts, addRays = TRUE)
#' finalize3D()
#'
#' pts <- genNDSet(3, 100, dubND = FALSE)
#' pts <- as.data.frame(pts[,1:3])
#'
#' ini3D(argsPlot3d = list(
#' xlim = c(0,max(pts[,1])+10),
#' ylim = c(0,max(pts[,2])+10),
#' zlim = c(0,max(pts[,3])+10)))
#' plotHull3D(pts, drawPoints = TRUE, addRays = TRUE)
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(
#' xlim = c(0,max(pts[,1])+10),
#' ylim = c(0,max(pts[,2])+10),
#' zlim = c(0,max(pts[,3])+10)))
#' plotHull3D(pts, drawPoints = TRUE, drawPolygons = TRUE, addText = "coord", addRays = TRUE)
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(
#' xlim = c(0,max(pts[,1])+10),
#' ylim = c(0,max(pts[,2])+10),
#' zlim = c(0,max(pts[,3])+10)))
#' plotHull3D(pts, drawPoints = TRUE, drawLines = FALSE,
#' argsPolygon3d = list(alpha = 1), addRays = TRUE)
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(
#' xlim = c(0,max(pts[,1])+10),
#' ylim = c(0,max(pts[,2])+10),
#' zlim = c(0,max(pts[,3])+10)))
#' plotHull3D(pts, drawPoints = TRUE, argsPolygon3d = list(color = "red"), addRays = TRUE)
#' plotCones3D(pts, argsPolygon3d = list(alpha = 1), rectangle = TRUE)
#' finalize3D()
#' }
plotHull3D <- function(pts,
drawPoints = FALSE,
drawLines = TRUE,
drawPolygons = TRUE,
addText = FALSE,
addRays = FALSE,
useRGLBBox = TRUE,
direction = 1,
drawBBoxHull = TRUE,
...)
{
args <- list(...)
argsSegments3d <- mergeLists(list(lwd = 1, col = "grey40"), args$argsSegments3d)
argsPlot3d <- mergeLists(list(size = 5, col = "black"), args$argsPlot3d)
argsShade3d <- mergeLists(list(col = "grey100"), args$argsShade3d)
argsPolygon3d <- mergeLists(list(color = "grey100", lwd = 2, alpha = 0.1),
args$argsPolygon3d)
argsText3d <- mergeLists(list(), args$argsText3d)
ids <- NULL
pts <- .checkPts(pts, p = 3)
if (length(direction) != 3) direction = rep(direction[1],3)
hull <- convexHull(pts, addRays = addRays, direction = direction, useRGLBBox = useRGLBBox)
set <- hull$pts
hull <- hull$hull
d <- dimFace(set[,1:3, drop = FALSE])
if (d==3) { # then poly define facets
poly <- hull
v <- 1:3
names(v) <- colnames(set[,1:3])
pN <- purrr::map_dfc(v, function(i) if (sign(direction[i]) > 0) max(set[,i]) else min(set[,i]))
#colnames(pN) <- colnames(set[,1:3])
for (i in 1:dim(poly)[1]) {
tri <- poly[i,!is.na(poly[i,])]
pt <- set[tri,1:4]
tri <- 1:nrow(pt)
if (drawLines) {
if (length(tri)>3) { # then have to find the vertex sequence
tri <- convexHull(pt[,1:3])$hull
}
tri1 <- NULL # use tri1 to include the addRays=TRUE case
for (j in 2:length(tri)){
if (!(pt[tri[j-1],4] == 0 && pt[tri[j],4] == 0)) # && !drawBBoxHull
tri1 <- c(tri1,NA,tri[j-1],tri[j])
}
if (!(pt[tri[1], 4] == 0 && pt[tri[length(tri)], 4] == 0)) { # && !drawBBoxHull
tri1 <- c(tri1, NA, tri[1], tri[length(tri)])
}
ids <- c(ids, do.call(rgl::polygon3d, args = c(list(pt[tri1, 1], pt[tri1, 2], pt[tri1, 3], fill = FALSE),
argsSegments3d)))
}
if (drawPolygons) {
if (drawBBoxHull | (!drawBBoxHull & nrow(dplyr::intersect(pt[,1:3],pN)) == 0)) {
if (length(tri)==3) {
ids <- c(ids, do.call(rgl::triangles3d, args = c(list(x = pt[,1:3]), argsPolygon3d) ))
} else if (length(tri)==4){
tri <- convexHull(pt[,1:3])$hull # then have to find the vertex sequence
obj <- rgl::qmesh3d(t(pt[,1:3]),tri, homogeneous = FALSE)
ids <- c(ids, do.call(rgl::shade3d, args = c(list(obj), argsPolygon3d)))
} else {
# idx <- apply(pt[,1:3], 2, function(x) {return(length(unique(x))==1)})
# idx<-which(!idx)
# if (length(idx)==3) coords <- 1:2 else coords <- idx
#plotHull3D(pt[tri,1:3])
tri <- convexHull(pt[,1:3])$hull
comb <- t(utils::combn(3,2))
for (j in 1:3) {
res <- try(
ids <- c(ids, do.call(rgl::polygon3d, args = c(list(
pt[tri, 1], pt[tri, 2], pt[tri, 3], fill = TRUE, coords = comb[j,]
), argsPolygon3d))), silent = TRUE)
if (!inherits(res, "try-error")) break #else {cat(i, " ", j, " ")}
}
}
}
}
}
}
if (d==0) {
ids <- c(ids, do.call(rgl::plot3d, args = c(list(set[,1], set[,2], set[,3], add=TRUE), argsPlot3d)))
}
if (d==1) { #segments3d(set, col = "black", lwd=10, smooth= TRUE)
cyl <- rgl::cylinder3d(set, radius = 0.01)
ids <- c(ids, do.call(rgl::shade3d, args = c(list(cyl), argsShade3d)))
}
if (d==2) {
idx <- apply(set[,1:3], 2, function(x) {return(length(unique(x))==1)})
tri <- hull
if (length(tri)==3) {
if (drawPolygons) ids <- c(ids, do.call(rgl::triangles3d,
args = c(list(set[tri,1], set[tri,2], set[tri,3]), argsPolygon3d) ))
if (drawLines) ids <- c(ids, do.call(rgl::polygon3d,
args = c(list(set[tri,1], set[tri,2], set[tri,3], coords = which(!idx), fill= FALSE),
argsSegments3d) ))
} else {
if (drawPolygons) ids <- c(ids, do.call(rgl::polygon3d,
args = c(list(set[tri,1], set[tri,2], set[tri,3], coords = which(!idx), fill= TRUE),
argsPolygon3d) ))
if (drawLines) ids <- c(ids, do.call(rgl::polygon3d,
args = c(list(set[tri,1], set[tri,2], set[tri,3], coords = which(!idx), fill= FALSE),
argsSegments3d) ))
}
}
if (drawPoints) ids <- c(ids, do.call(plotPoints3D, args = c(list(pts), list(argsPlot3d = argsPlot3d))))
if (length(addText) > 1) {
if (length(addText) == nrow(pts)) {
do.call(rgl::text3d, args = c(list(x = pts, texts = addText), argsText3d) )
}
} else if (length(addText) == 1) {
if (addText == TRUE | addText == "coord") {
text <- paste0("(",pts[,1],",",pts[,2],",",pts[,3],")")
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = pts, texts = text), argsText3d) ))
} else if (addText == "rownames") {
text <- rownames(as.data.frame(pts))
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = pts, texts = text), argsText3d) ))
} else if (addText == "both") {
text <- paste0("(",pts[,1],",",pts[,2],",",pts[,3],")")
text <- paste(text,rownames(as.data.frame(pts)))
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = pts, texts = text), argsText3d) ))
} else if (addText != FALSE & length(addText) == nrow(pts)) {
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = pts, texts = addText), argsText3d) ))
}
}
return(invisible(list(hull = hull, pts = set, ids = ids)))
stop("Error in plotHull3D")
}
#' Plot points in 3D.
#'
#' @param pts A vector or matrix with the points.
#' @param addText Add text to the points. Currently `coord` (coordinates), `rownames` (rownames)
#' and `both` supported or a vector with the text.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`].
#' * `argsPch3d`: A list of arguments for [`rgl::pch3d`].
#' * `argsText3d`: A list of arguments for [`rgl::text3d`].
#'
#' @return Object ids (invisible).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D()
#' pts<-matrix(c(1,1,1,5,5,5), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' plotPoints3D(c(2,3,3), argsPlot3d = list(col = "red", size = 10))
#' plotPoints3D(c(3,2,3), argsPlot3d = list(col = "blue", size = 10, type="p"))
#' plotPoints3D(c(1.5,1.5,1.5), argsPlot3d = list(col = "blue", size = 10, type="p"))
#' plotPoints3D(c(2,2,2, 1,1,1), addText = "coord")
#' ids <- plotPoints3D(c(3,3,3, 4,4,4), addText = "rownames")
#' finalize3D()
#' rgl::rglwidget()
#' # pop3d(ids) # remove the last again
#' }
plotPoints3D <- function(pts, addText = FALSE, ...) {
args <- list(...)
argsPlot3d <- mergeLists(list(size = 5, col = "black", type="p"), args$argsPlot3d)
argsPch3d <- mergeLists(list(litt = TRUE), args$argsPch3d)
argsText3d <- mergeLists(list(), args$argsText3d)
ids <- NULL
p <- pts
if (is.vector(pts)) p <- matrix(pts, ncol = 3, byrow = TRUE)
p <- as.matrix(p)
if (is.null(args$argsPch3d)) ids <- c(ids, do.call(rgl::plot3d, args = c(list(p, add= TRUE), argsPlot3d) ))
if (!is.null(args$argsPch3d)) ids <- c(ids, do.call(rgl::pch3d, args = c(list(p), argsPch3d) ))
if (length(addText) > 1) {
if (length(addText) == nrow(p)) {
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = p, texts = addText), argsText3d) ))
}
} else if (length(addText) == 1) {
if (addText == TRUE | addText == "coord") {
text <- paste0("(",p[,1],",",p[,2],",",p[,3],")")
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = p, texts = text), argsText3d) ))
} else if (addText == "rownames") {
text <- rownames(as.data.frame(p))
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = p, texts = text), argsText3d) ))
} else if (addText == "both") {
text <- paste0("(",p[,1],",",p[,2],",",p[,3],")")
text <- paste(text,rownames(as.data.frame(p)))
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = p, texts = text), argsText3d) ))
} else if (addText != FALSE & length(addText) == nrow(p)) {
ids <- c(ids, do.call(rgl::text3d, args = c(list(x = p, texts = addText), argsText3d) ))
}
}
return(invisible(ids))
}
#' Plot a plane in 3D.
#'
#' @param normal Normal to the plane.
#' @param point A point on the plane.
#' @param offset The offset of the plane (only used if `point = NULL`).
#' @param useShade Plot shade of the plane.
#' @param useLines Plot lines inside the plane.
#' @param usePoints Plot point shapes inside the plane.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists (see examples). Currently the following arguments are supported:
#'
#' * `argsPlanes3d`: A list of arguments for [rgl::planes3d()] used when `useShade = TRUE`.
#' * `argsLines`: A list of arguments for [rgl::persp3d()] when `useLines = TRUE`. Moreover, the list
#' may contain `lines`: number of lines.
#'
#' @return NULL (invisible)
#' @export
#'
#' @examples
#' \donttest{
#' ini3D(argsPlot3d = list(xlim = c(-1,10), ylim = c(-1,10), zlim = c(-1,10)) )
#' plotPlane3D(c(1,1,1), point = c(1,1,1))
#' plotPoints3D(c(1,1,1))
#' plotPlane3D(c(1,2,1), point = c(2,2,2), argsPlanes3d = list(color="red"))
#' plotPoints3D(c(2,2,2))
#' plotPlane3D(c(2,1,1), offset = -6, argsPlanes3d = list(color="blue"))
#' plotPlane3D(c(2,1,1), argsPlanes3d = list(color="green"))
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(xlim = c(-1,10), ylim = c(-1,10), zlim = c(-1,10)) )
#' plotPlane3D(c(1,1,1), point = c(1,1,1), useLines = TRUE, useShade = TRUE)
#' ids <- plotPlane3D(c(1,2,1), point = c(2,2,2), argsLines = list(col="blue", lines = 100),
#' useLines = TRUE)
#' finalize3D()
#' # pop3d(id = ids) # remove last plane
#' }
plotPlane3D <- function(normal, point = NULL, offset = 0, useShade = TRUE, useLines = FALSE,
usePoints = FALSE, ...) {
args <- list(...)
argsPlanes3d <- mergeLists(list(col = "grey", alpha = 0.2), args$argsPlanes3d)
argsLines <- mergeLists(list(back = 'lines', front = 'lines', add = TRUE, lines = 50),
args$argsLines)
ids <- NULL
if (!is.null(point)) offset <- -sum(normal * point)
if (useShade) {
ids <- c(ids, do.call(rgl::planes3d, args = c(list(a = normal, d = offset), argsPlanes3d) ))
}
# else use points or lines
if (!rgl::cur3d()) stop("Option useLines or usePoints need an open rgl window!")
limits <- rgl::par3d()$bbox
m <- c(limits[1], limits[3], limits[5])
M <- c(limits[2], limits[4], limits[6])
# do the mesh
x <- seq(m[1], M[1], length.out = argsLines$lines)
y <- seq(m[2], M[2], length.out = argsLines$lines)
f <- function(x,y){-(normal[1]/normal[3])*x - (normal[2]/normal[3])*y - offset/normal[3]}
z <- outer(x,y,f)
z[z < m[3] | z > M[3]] <- NA
if (useLines) {
ids <- c(ids, do.call(rgl::persp3d, args = c(list(x, y, z), argsLines) ))
}
if (usePoints) {
return(invisible(ids))
# do.call(rgl::plot3d, args = c(list(x, y, z, add = TRUE), argsPoints) )
# persp3d(x, y, z, back = 'lines', front = 'lines', add = TRUE)
#
#
# if (nrow(pts) > 3) poly <- do.call(rgl::polygon3d, args = c(list(pts, plot = FALSE),
# argsPolygon3d))
# if (nrow(pts) == 3) {
# # poly <- do.call(rgl::triangles3d, args = c(list(pts, plot = FALSE), argsPolygon3d))
# poly <- tmesh3d(rbind(pts[,1], pts[,2], pts[,3], 1), indices = 1:3)
# }
# poly$texcoords <- argsPolygon3d$texcoords
# do.call(rgl::shade3d, args = c(list(poly, col = "white", specular = "black"), argsPolygon3d))
}
return(invisible(ids))
}
#' The `ggplot` theme for the package
#'
#' @param ... Further arguments parsed to [ggplot2::theme()].
#'
#' @return The theme object.
#' @export
#'
#' @examples
#' pts <- matrix(c(1,1), ncol = 2, byrow = TRUE)
#' plotHull2D(pts)
#' pts1 <- matrix(c(2,2, 3,3), ncol = 2, byrow = TRUE)
#' pts2 <- matrix(c(1,1, 2,2, 0,1), ncol = 2, byrow = TRUE)
#' ggplot2::ggplot() +
#' plotHull2D(pts2, drawPoints = TRUE, addText = "coord", drawPlot = FALSE) +
#' plotHull2D(pts1, drawPoints = TRUE, drawPlot = FALSE) +
#' gMOIPTheme() +
#' ggplot2::xlab(expression(x[1])) +
#' ggplot2::ylab(expression(x[2]))
gMOIPTheme <- function(...) {
return(
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(),
#axis.line = element_blank(),
axis.line = element_line(colour = "black", size = 0.5,
arrow = arrow(length = unit(0.3,"cm")) ),
#axis.ticks = element_blank()
#axis.text.x = element_text(margin = margin(r = 30))
# axis.ticks.length = unit(0.5,"mm"),
#aspect.ratio=4/3,
legend.position="none",
...
)
)
}
#' Initialize the RGL window.
#'
#' @param new A new window is opened (otherwise the current is cleared).
#' @param clear Clear the current RGL window.
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsPlot3d`: A list of arguments for [`rgl::plot3d`].
#' * `argsAspect3d`: A list of arguments for [`rgl::aspect3d`].
#'
#' @return NULL (invisible).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D()
#' pts<-matrix(c(1,1,1,5,5,5), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' finalize3D()
#'
#' lim <- c(-1, 7)
#' ini3D(argsPlot3d = list(xlim = lim, ylim = lim, zlim = lim))
#' plotPoints3D(pts)
#' finalize3D()
#' }
ini3D <- function(new = TRUE, clear = FALSE, ...){
args <- list(...)
# args <- list(argsPlot3d = list(xlim = c(0,6), ylim = c(0,6), zlim = c(0,6)))
argsPlot3d <-
mergeLists(list(
xlab = '',
ylab = '',
zlab = '',
box = FALSE,
axes = FALSE
), args$argsPlot3d)
argsAspect3d <- mergeLists(list(x = "iso"), args$argsAspect3d)
if (new) rgl::open3d()
if (clear) rgl::clear3d()
do.call(rgl::plot3d, args = c(list(x = c(1,1), y = c(1,1), z = c(1,1), type = 'n'), argsPlot3d))
do.call(rgl::aspect3d, args = argsAspect3d)
rgl::highlevel()
return(invisible(NULL))
}
#' Finalize the RGL window.
#'
#' @param ... Further arguments passed on the the RGL plotting functions. This must be done as
#' lists. Currently the following arguments are supported:
#'
#' * `argsAxes3d`: A list of arguments for [`rgl::axes3d`].
#' * `argsTitle3d`: A list of arguments for [`rgl::title3d`][rgl::axes3d].
#'
#' @return The RGL object (using [rgl::highlevel()]).
#' @export
#'
#' @examples
#' \donttest{
#' ini3D()
#' pts<-matrix(c(1,1,1,5,5,5), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' finalize3D()
#'
#' ini3D()
#' pts<-matrix(c(1,1,1,5,5,5), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' finalize3D(argsAxes3d = list(edges = "bbox"))
#' }
finalize3D <- function(...){
args <- list(...)
argsAxes3d <- mergeLists(list(edges = c('x', 'y', 'z')), args$argsAxes3d)
argsTitle3d <- mergeLists(list(xlab = expression(italic(z)[1]), ylab = expression(italic(z)[2]),
zlab = expression(italic(z)[3])), args$argsTitle3d)
do.call(rgl::axes3d, args = argsAxes3d)
do.call(rgl::title3d, args = argsTitle3d)
rgl::highlevel()
}
#' Convert LaTeX to a png file
#'
#' @param tex TeX string. Remember to escape backslash with \\.
#' @param viewPng View the result in the plots window.
#' @param width Width of the png.
#' @param height Height of the png (`width` are ignored).
#' @param dpi Dpi of the png. Not used if `width` or `height` are specified.
#' @param fontsize Front size used in the LaTeX document.
#' @param calcM Estimate 1 em in pixels in the resulting png.
#' @param crop Call command line program `pdfcrop` (must be installed).
#'
#' @return The filename of the png or a list if `calcM = TRUE`.
#' @export
#'
#' @examples
#' \dontrun{
#' tex <- "$\\mathbb{R}_{\\geqq}$"
#' texToPng(tex, viewPng = TRUE)
#' texToPng(tex, fontsize = 20, viewPng = TRUE)
#' texToPng(tex, height = 50, fontsize = 10, viewPng = TRUE)
#' texToPng(tex, height = 50, fontsize = 50, viewPng = TRUE)
#' tex <- "MMM"
#' texToPng(tex, dpi=72, calcM = TRUE)
#' texToPng(tex, width = 100, calcM = TRUE)
#' f <- texToPng(tex, dpi=300)
#' pngSize(f)
#' }
texToPng <- function(tex, width = NULL, height = NULL, dpi = 72, viewPng = FALSE, fontsize = 12,
calcM = FALSE, crop = FALSE) {
texFile <- tempfile(fileext=".tex")
pdfFile <- gsub("[.]tex", ".pdf", texFile)
pngFile <- gsub("[.]tex", ".png", texFile)
writeLines(c(paste0("\\documentclass[class=scrreprt, fontsize = ", fontsize, "pt]{standalone}"),
"\\usepackage{amsfonts}",
"\\usepackage{amssymb}",
"\\begin{document}",
tex,
"\\end{document}"),
texFile)
system(paste0("pdflatex -shell-escape -output-directory=", tempdir(), " ", texFile))
if (crop) system(paste("pdfcrop", pdfFile, pdfFile))
tmp <- pdftools::pdf_pagesize(pdfFile)
z <- 1
if (!is.null(width)) z <- dpi/72 * tmp$width/width
if (!is.null(height)) z <- dpi/72 * tmp$height/height
dpi <- dpi * 1/z
if (calcM) em1 <- .sizeM(dpi = dpi, fontsize = fontsize)$w
pdftools::pdf_convert(pdfFile, format = "png", dpi = dpi, antialias = TRUE, filenames = pngFile, verbose = FALSE)
if (viewPng) {
img <- png::readPNG(pngFile)
grid::grid.newpage()
grid::grid.raster(img, just = "center", height = 0.2)
}
if (calcM) {
s <- pngSize(pngFile)
emPng <- s$w/em1
return(list(png=pngFile, emPx = em1, emLength = emPng, w = s$w, h = s$h))
}
return(pngFile)
}
#' Estimate 1 em in pixels in the resulting png.
#'
#' @param ... Arguments parsed on to `texToPng`.
#'
#' @return The width and size of the png.
#' @keywords internal
.sizeM <- function(...) {
f <- texToPng("M", ...)
return(pngSize(f))
}
#' To size of the png file.
#'
#' @param png Png file name.
#'
#' @return A list with width and height.
pngSize <- function(png) {
img <- png::readPNG(png)
w <- dim(img)[2]
h <- dim(img)[1]
return(list(w=w,h=h))
}
#' Draw boxes, axes and other text outside the data using TeX strings.
#'
#' @param main The main title for the plot.
#' @param sub The subtitle for the plot.
#' @param xlab The axis labels for the plot .
#' @param ylab The axis labels for the plot .
#' @param zlab The axis labels for the plot .
#' @param line The ``line'' of the plot margin to draw the label on.
#' @param ... Additional parameters which are passed to \code{\link{plotMTeX3D}}.
#'
#' @details The rectangular prism holding the 3D plot has 12 edges. They are identified
#' using 3 character strings. The first character (`x', `y', or `z') selects
#' the direction of the axis. The next two characters are each `-' or `+',
#' selecting the lower or upper end of one of the other coordinates. If only
#' one or two characters are given, the remaining characters default to `-'.
#' For example \code{edge = 'x+'} draws an x-axis at the high level of y and the
#' low level of z.
#'
#' By default, \code{axes3d} uses the \code{\link{bbox3d}} function to draw the axes.
#' The labels will move so that they do not obscure the data. Alternatively,
#' a vector of arguments as described above may be used, in which case
#' fixed axes are drawn using \code{axis3d}.
#'
#' If \code{pos} is a numeric vector of length 3, \code{edge} determines
#' the direction of the axis and the tick marks, and the values of the
#' other two coordinates in \code{pos} determine the position. See the
#' examples.
#' @return The object IDs of objects added to the scene.
#' @export
#'
#' @examples
#' \dontrun{
#' ini3D(argsPlot3d = list(xlim = c(0, 2), ylim = c(0, 2), zlim = c(0, 2)))
#' plotTitleTeX3D(main = "\\LaTeX", sub = "subtitle $\\alpha$",
#' xlab = "$x^1_2$", ylab = "$\\beta$", zlab = "$x\\cdot y$")
#' finalize3D()
#' }
plotTitleTeX3D <- function (main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
zlab = NULL, line = NA, ...) {
save <- rgl::par3d(skipRedraw = TRUE, ignoreExtent = TRUE)
on.exit(rgl::par3d(save))
result <- numeric(0)
if (!is.null(main)) {
aline <- ifelse(is.na(line), 2, line)
result <- c(result, main = plotMTeX3D(main, "x++",
line = aline, ...))
}
if (!is.null(sub)) {
aline <- ifelse(is.na(line), 3, line)
result <- c(result, sub = plotMTeX3D(sub, "x", line = aline,
...))
}
if (!is.null(xlab)) {
aline <- ifelse(is.na(line), 2, line)
result <- c(result, xlab = plotMTeX3D(xlab, "x", line = aline,
...))
}
if (!is.null(ylab)) {
aline <- ifelse(is.na(line), 2, line)
result <- c(result, ylab = plotMTeX3D(ylab, "y", line = aline,
...))
}
if (!is.null(zlab)) {
aline <- ifelse(is.na(line), 2, line)
result <- c(result, zlab = plotMTeX3D(zlab, "z", line = aline,
...))
}
lowlevel(result)
}
#' Get ranges of the bounding box margins
#'
#' @param expand Expand margins.
#' @param ranges The bounding box.
#'
#' @return A list with ranges.
#' @keywords internal
.getRanges <- function (expand = 1.03, ranges = par3d("bbox"))
{
ranges <- list(xlim = ranges[1:2], ylim = ranges[3:4], zlim = ranges[5:6])
strut <- FALSE
ranges <- lapply(ranges, function(r) {
d <- diff(r)
if (d > 0)
return(r)
strut <<- TRUE
if (d < 0)
return(c(0, 1))
else if (r[1] == 0)
return(c(-1, 1))
else return(r[1] + 0.4 * abs(r[1]) * c(-1, 1))
})
ranges$strut <- strut
ranges$x <- (ranges$xlim - mean(ranges$xlim)) * expand +
mean(ranges$xlim)
ranges$y <- (ranges$ylim - mean(ranges$ylim)) * expand +
mean(ranges$ylim)
ranges$z <- (ranges$zlim - mean(ranges$zlim)) * expand +
mean(ranges$zlim)
ranges
}
#' Plot TeX in the margin
#'
#' @param tex TeX string
#' @param edge The position at which to draw the axis or text.
#' @param line The ``line'' of the plot margin to draw the label on.
#' @param at The value of a coordinate at which to draw the axis.
#' @param pos The position at which to draw the axis or text.
#' @param ... Further arguments passed to [plotTeX3D()].
#'
#' @return The object IDs of objects added to the scene.
#' @export
plotMTeX3D <- function (tex, edge, line = 0, at = NULL, pos = NA, ...) {
save <- rgl::par3d(ignoreExtent = TRUE)
on.exit(rgl::par3d(save))
ranges <- .getRanges()
edge <- c(strsplit(edge, "")[[1]], "-", "-")[1:3]
coord <- match(toupper(edge[1]), c("X", "Y", "Z"))
if (coord == 2)
edge[1] <- edge[2]
else if (coord == 3)
edge[1:2] <- edge[2:3]
range <- ranges[[coord]]
if (is.null(at))
at <- mean(range)
newlen <- max(length(tex), length(line), length(at))
tex <- rep(tex, len = newlen)
line <- rep(line, len = newlen)
at <- rep(at, len = newlen)
if (all(is.na(pos))) {
pos <- matrix(NA, 3, length(at))
if (edge[1] == "+")
pos[1, ] <- ranges$x[2]
else pos[1, ] <- ranges$x[1]
if (edge[2] == "+")
pos[2, ] <- ranges$y[2]
else pos[2, ] <- ranges$y[1]
if (edge[3] == "+")
pos[3, ] <- ranges$z[2]
else pos[3, ] <- ranges$z[1]
}
else pos <- matrix(pos, 3, length(at))
pos[coord, ] <- at
ticksize <- 0.05 * (pos[, 1] - c(mean(ranges$x), mean(ranges$y),
mean(ranges$z)))
ticksize[coord] <- 0
plotTeX3D(pos[1, ] + 3 * ticksize[1] * line,
pos[2, ] + 3 * ticksize[2] * line,
pos[3, ] + 3 * ticksize[3] * line, tex,
...)
}
#' Plot TeX at a position.
#'
#' @param x Coordinate.
#' @param y Coordinate.
#' @param z Coordinate.
#' @param cex Expansion factor (you properly have to fine tune it).
#' @param fixedSize Fix the size of the object (no scaling when zoom).
#' @param tex TeX string.
#' @param size Size of the generated png.
#' @param ... Arguments passed on to [rgl::sprites3d()] and [texToPng()].
#'
#' @return The shape ID of the displayed object is returned.
#' @export
#'
#' @examples
#' \dontrun{
#' tex0 <- "$\\mathbb{R}_{\\geqq}$"
#' tex1 <- "\\LaTeX"
#' tex2 <- "This is a title"
#' ini3D(argsPlot3d = list(xlim = c(0, 2), ylim = c(0, 2), zlim = c(0, 2)))
#' plotTeX3D(0.75,0.75,0.75, tex0)
#' plotTeX3D(0.5,0.5,0.5, tex0, cex = 2)
#' plotTeX3D(1,1,1, tex2)
#' finalize3D()
#' ini3D(new = TRUE, argsPlot3d = list(xlim = c(0, 200), ylim = c(0, 200), zlim = c(0, 200)))
#' plotTeX3D(75,75,75, tex0)
#' plotTeX3D(50,50,50, tex1)
#' plotTeX3D(100,100,100, tex2)
#' finalize3D()
#' }
plotTeX3D <- function (x, y, z, tex, cex = graphics::par("cex"), fixedSize = FALSE, size = 480, ...) {
f <- texToPng(tex, width = size, calcM = TRUE, ...)
# expand png so same width and height
system(paste0("convert ", f$png,
" -background none -gravity center -extent ", max(f$w,f$h), "x", max(f$w,f$h), " ", f$png))
tmp <- rgl::par3d()$bbox
radius <- 1/2 * cex * f$emLength * 20/f$emPx * max(c(tmp[2]-tmp[1], tmp[4]-tmp[3],tmp[6]-tmp[5]) * rgl::par3d()$scale)
rgl::sprites3d(x, y, z, texture = f$png,
textype = "rgba", col = "white", lit = FALSE,
radius = radius, fixedSize = fixedSize, ...)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.