# R/hull.R In gMOIP: Tools for 2D and 3D Plots of Single and Multi-Objective Linear/Integer Programming Models

```## Functions related to the convex hull

#' Return the dimension of the convex hull of a set of points.
#'
#' @param pts A matrix/data frame/vector that can be converted to a matrix with a row for each
#'   point.
#' @param dim The dimension of the points, i.e. assume that column 1-dim specify the points. If
#'   NULL assume that the dimension are the number of columns.
#'
#' @return The dimension of the object.
#' @export
#'
#' @examples
#' ## In 1D
#' pts <- matrix(c(3), ncol = 1, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(1,3,4), ncol = 1, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 2D
#' pts <- matrix(c(3,3,6,3,3,6), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(1,1,2,2,3,3), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 3D
#' pts <- c(3,3,3,6,3,3,3,6,3,6,6,3)
#' dimFace(pts, dim = 3)
#' pts <- matrix( c(1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix( c(1,1,1,2,2,2), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(2,2,2,3,2,2), ncol=3, byrow= TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0,0,0,1,1,0,2,2,0,5,2,0,6,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0,0,0,1,1,0,2,2,0,0,2,1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 4D
#' pts <- matrix(c(2,2,2,3,2,2,3,4,1,2,3,4), ncol=4, byrow= TRUE)
#' dimFace(pts,)
dimFace<-function(pts, dim = NULL) {
if (is.null(dim)) dim <- ncol(pts)
if (is.vector(pts)) pts <- matrix(pts, ncol = dim, byrow = TRUE)
x <- as.matrix(unique(pts[,1:dim,drop=FALSE]))
if (nrow(x)==1) return(0)
m1 <- x[2:dim(x),]
m2 <- matrix(rep(x[1,], times = dim(x)-1), ncol = dim(x), byrow = TRUE)
x <- m1 - m2
return(Matrix::rankMatrix(x))
}

#' Efficient test for points inside a convex hull in p dimensions.
#'
#' @param pts A \eqn{nxp} array to test, \eqn{n} data points, in dimension \eqn{p}. If you have
#'   many points to test, it is most efficient to call this function once with the entire set.
#' @param vertices A \eqn{mxp} array of vertices of the convex hull. May contain redundant
#'   (non-vertex) points.
#' @param hull Tessellation (or triangulation) generated by convhulln (only works if the dimension
#'   of the hull is \eqn{p}). If hull is \code{NULL}, then it will be generated.
#' @param tol Tolerance on the tests for inclusion in the convex hull. You can think of `tol` as the
#'   difference a point value may be different from the values of the hull, and still be perceived
#'   as on the surface of the hull. Because of numerical slop nothing can ever be done exactly here.
#'   In higher dimensions, the numerical issues of floating point arithmetic will probably suggest a
#'   larger value of `tol`. `tol` is not used if the dimension of the hull is larger than one and
#'   not equal \eqn{p}.
#'
#' @note Some of the code are inspired by the [Matlab
#'   code](https://www.mathworks.com/matlabcentral/fileexchange/10226-inhull) by
#'   John D'Errico and
#'   [how to find a point inside a hull](https://stat.ethz.ch/pipermail/r-help/2009-December/415377.html).
#'   If the dimension of the hull is below \eqn{p} then PCA may be used to check (a
#'   warning will be given).
#'
#' @return An integer vector of length \eqn{n} with values 1 (inside hull), -1 (outside hull) or 0
#'   (on hull to precision indicated by tol).
#'
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#'
#' @examples
#' ## In 1D
#' vertices <- matrix(4, ncol = 1)
#' pt <- matrix(c(2,4), ncol = 1, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(1,4), ncol = 1)
#' pt <- matrix(c(1,3,4,5), ncol = 1, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' ## In 2D
#' vertices <- matrix(c(2,4), ncol = 2)
#' pt <- matrix(c(2,4, 1,1), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(0,0, 3,3), ncol = 2, byrow = TRUE)
#' pt <- matrix(c(0,0, 1,1, 2,2, 3,3, 4,4), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(0,0, 0,3, 3,0), ncol = 2, byrow = TRUE)
#' pt <- matrix(c(0,0, 1,1, 4,4), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' ## in 3D
#' vertices <- matrix(c(2,2,2), ncol = 3, byrow = TRUE)
#' pt <- matrix(c(1,1,1, 3,3,3, 2,2,2, 3,3,2), ncol = 3, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(2,2,2, 4,4,4), ncol = 3, byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,1,1, 2,2,2, 3,3,3, 4,4,4, 3,3,2), ncol = 3, byrow = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(1,0,0, 1,1,0, 1,0,1), ncol = 3, byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,0.1,0.2, 3,3,2), ncol = 3, byrow = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(2,2,2, 2,4,4, 2,2,4, 4,4,2, 4,2,2, 2,4,2, 4,2,4, 4,4,4), ncol = 3,
#'             byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,1,1, 3,3,3, 2,2,2, 3,3,2), ncol = 3, byrow = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' ## In 5D
#' vertices <- matrix(c(4,0,0,0,0, 0,4,0,0,0, 0,0,4,0,0, 0,0,0,4,0, 0,0,0,0,4, 0,0,0,0,0),
#'             ncol = 5, byrow = TRUE)
#' pt <- matrix(c(0.1,0.1,0.1,0.1,0.1, 3,3,3,3,3, 2,0,0,0,0), ncol = 5, byrow = TRUE)
#' inHull(pt, vertices)
inHull <- function(pts, vertices, hull=NULL,
tol=mean(mean(abs(as.matrix(vertices))))*sqrt(.Machine\$double.eps)) {
# ensure arguments are matrices (not data frames) and get sizes
pts <- .checkPts(pts)
vertices <- .checkPts(vertices)

p <- ncol(vertices)   # dimension of points
cx <- nrow(pts)   # points to test
d <- dimFace(vertices, dim = p)
if (d == 0) return(-as.integer(apply(pts, 1, FUN = function(x) abs(x - vertices[1,]) > tol)))
if (p == 1) { #1D space - check if between max and min
m <- min(vertices)
M <- max(vertices)
val <- rep(0, cx)  # first set all on hull
val[pts > M + tol | pts < m - tol] <- -1
val[pts < M - tol & pts > m + tol] <- 1
return(val)
}
if (d == 1) { # a line (i.e. no interior points)
denom <- vertices[2,]-vertices[1,]  # the direction d of the line p = d * t + p0
# case: if d_i = 0 => p_i = p0_i otherwise not on line
zeros <- denom == 0
if (any(zeros)) {
val <- apply(pts[, zeros, drop=F], 1, FUN = function(x) {any(x != vertices[1, zeros])})  # are some points not on the line
# reduce dimension where d_i = 0
vertices <- vertices[, !zeros, drop = F]
pts <- pts[, !zeros, drop=F]
res <- inHull(pts, vertices, tol = tol)
res[res == 1] <- 0   # set interior pts to true value
res[res == -1] <- 1
res <- as.logical(res)
return(-as.integer(res | val))
} else {
l <- apply(pts, 1, FUN = function(x) {(x-vertices[1,])/denom})  # calc t values t = d_i/(p_i-p0_i)
tRng <- range(apply(vertices, 1, FUN = function(x) {(x - vertices[1,1])/denom})) # t range for vertices
eq <- apply(l, 2, FUN = function(x) abs(max(x) - min(x)) < tol) # check if values in each col are equal (i.e. on line)
l <- apply(l, 2, FUN = function(x) max(x) < tRng + tol & min(x) > tRng - tol) # check if t values in range
l <- eq & l
return(-as.integer(!l))
}
}
if (p == d) {
if (is.null(hull)) hull <- geometry::convhulln(vertices)
nt <- nrow(hull)    # number of simplexes in hull
# find normal vectors to each simplex
nrmls <- matrix(NA, nt, p)         # predefine each nrml as NA, degenerate
degenflag <- matrix(TRUE, nt, 1)
for (i in  1:nt) {
nullsp <-
t(MASS::Null(t(vertices[hull[i,-1],] - matrix(vertices[hull[i, 1],], p - 1, p, byrow = TRUE))))
if (dim(nullsp) == 1) {
nrmls[i, ] <- nullsp
degenflag[i] <- FALSE
}
}
# Warn of degenerate faces, and remove corresponding normals
if (sum(degenflag) > 0)
warning(sum(degenflag), " degenerate faces in convex hull")
nrmls <- nrmls[!degenflag, ]
nt <- nrow(nrmls)
# find center point in hull, and any (1st) point in the plane of each simplex
center = colMeans(vertices)
a <- vertices[hull[!degenflag, 1], ]
# scale normal vectors to unit length and ensure pointing inwards
nrmls <- nrmls / matrix(apply(nrmls, 1, function(x) sqrt(sum(x ^ 2))), nt, p)
dp <- sign(apply((matrix(center, nt, p, byrow = TRUE) - a) * nrmls, 1, sum))
nrmls <- nrmls * matrix(dp, nt, p)
# if  min across all faces of dot((x - a),nrml) is
#      +ve then x is inside hull
#      0   then x is on hull
#      -ve then x is outside hull
# Instead of dot((x - a),nrml)  use dot(x,nrml) - dot(a, nrml)
aN <- diag(a %*% t(nrmls))
val <- apply(pts %*% t(nrmls) - matrix(aN, cx, nt, byrow = TRUE), 1, min)
# code  values inside 'tol' to zero, return sign as integer
val[abs(val) < tol] <- 0
return(as.integer(sign(val)))
}
if (p != d) {  # note p>2
outside <- rep(2,cx)
row.names(pts) <- 1:nrow(pts)
nrp <- nrow(vertices)
done <- FALSE
while(!done){
pt <- pts[which(outside != -1),,drop=F]
set <- rbind(vertices,pt)
dS <- dimFace(set, dim = p)
if (p == dS) {
val <- geometry::convhulln(set)
idx <- as.integer(row.names(pt[unique(val[val>nrp])-nrp,,drop=F]))
outside[idx] <- -1
done <- all(outside != 2)
}
if (dS != p) {
idx <- as.integer(row.names(pt))
warning(
"Using PCA to transform. Results for pts index ",
paste(idx, collapse = ","),
" might not be accurate!"
)
set <- stats::princomp(set)
idx <- which(set\$sdev > tol)
set <- stats::predict(set)[,idx]
# set <- prcomp(set, scale. = T, retx = T, rank = dS)\$x
val <- inHull(set[(nrp+1):nrow(set),], set[1:nrp,], tol = 5*tol)
val[val == 1] <- 0
outside[as.integer(row.names(pt))] <- val
return(outside)
}
}
return(outside)
}
}

#' Find segments (lines) of a face.
#'
#' @param vertices A mxp array of vertices of the convex hull, as used by
#'   convhulln.
#' @param hull Tessellation (or triangulation) generated by convhulln If hull is
#'   left empty or not supplied, then it will be generated.
#' @param tol Tolerance on the tests for inclusion in the convex hull. You can
#'   think of tol as the distance a point may possibly lie outside the hull, and
#'   still be perceived as on the surface of the hull. Because of numerical slop
#'   nothing can ever be done exactly here. I might guess a semi-intelligent
#'   value of tol to be
#'
#'   tol = 1.e-13*mean(abs(vertices(:)))
#'
#'   In higher dimensions, the numerical issues of floating point arithmetic
#'   will probably suggest a larger value of tol.
#'
#' @return A matrix with segments.
#' @author Lars Relund \email{lars@@relund.dk}
hullSegment <- function(vertices, hull=geometry::convhulln(vertices),
tol=mean(mean(abs(vertices)))*sqrt(.Machine\$double.eps)) {
vertices <- as.matrix(vertices)
p <- ncol(vertices)   # columns in vertices
nt <- nrow(hull)    # number of simplexes in hull
# find normal vectors to each simplex
nrmls <- matrix(NA, nt, p)         # predefine each nrml as NA, degenerate
degenflag <- matrix(TRUE, nt, 1)
for (i in  1:nt) {
nullsp <-
t(MASS::Null(t(
vertices[hull[i,-1],] - matrix(vertices[hull[i, 1],], p - 1, p, byrow =
TRUE)
)))
if (dim(nullsp) == 1) {
nrmls[i, ] <- nullsp
degenflag[i] <- FALSE
}
}
nt <- nrow(nrmls)
# find center point in hull, and any (1st) point in the plane of each simplex
center = colMeans(vertices)
a <- vertices[hull[, 1], ]
# scale normal vectors to unit length and ensure pointing inwards
nrmls <- nrmls / matrix(apply(nrmls, 1, function(x) sqrt(sum(x ^ 2))), nt, p)
dp <- sign(apply((matrix(center, nt, p, byrow = TRUE) - a) * nrmls, 1, sum))
nrmls <- nrmls * matrix(dp, nt, p)
idx <- !duplicated(round(nrmls,10))
nrmls <- nrmls[idx,]
a <- a[idx,]
aN <- diag(a %*% t(nrmls))
res <- NULL
m <- nrow(vertices)
for( i in 1:m )
for( j in 1:m )
if( i < j ) {
# Middle of the segment
p <- .5 * vertices[i,] + .5 * vertices[j,]
# Check if it is at the intersection of two planes
tmp <- p %*% t(nrmls) - aN #matrix(aN, 1, nrow(nrmls), byrow = TRUE)
if(sum( abs( tmp ) < 1e-6 ) >= 2) res <- rbind(res, c(i,j))
}
return(as.matrix(res))
}

#' Add all points on the bounding box hit by the rays.
#'
#' @param pts A data frame with all points
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of the `pts`
#'   plus a value greater than on equal zero. If negative, consider the i'th column of the `pts`
#'   minus a value greater than on equal zero.
#'
#' @note Assume that `pts` has been checked using [.checkPts].
#' @return The points merged with the points on the bounding box. The column `pt` equals 1 if
#'   points from pts and zero otherwise.
#' @export
#'
#' @examples
#' pts <- genNDSet(3,10)[,1:3]
#' addRays(pts, dir = c(-1,-1,1), m = c(0,0,0), M = c(100,100,100))
#' pts <- genSample(5,20)[,1:5]
addRays <- function(pts, m = apply(pts,2,min)-5, M = apply(pts,2,max)+5, direction = 1) {
pts <- as.data.frame(pts)
p <- ncol(pts)
if (length(direction) != p) direction = rep(direction,p)
v <- purrr::map_dbl(1:p, function(i) if (sign(direction[i]) > 0) M[i] else m[i])
names(v) <- colnames(pts)
set <- cbind(pts, pt = 1) # point in (1=original set, 0=artificial)

for (i in 1:nrow(pts)) {
x <- dplyr::bind_rows(pts[i,,drop=FALSE],v)
x <- as.list(x)
x <- expand.grid(x)
x\$pt <- 0
set <- rbind(set,x)
}
# colnames(set)[1:p] <- paste0("z",1:p)
# dplyr::distinct(set, V1, .keep_all = TRUE)
# rownames(set) <- 1:dim(set)
set <- set[rownames(unique(set[,1:p,drop=FALSE])),, drop=FALSE]
# rownames(set)[set\$pt==0] <- (nrow(pts)+1):nrow(set)  # may give error if number in first rows
return(set)
}

#' Find the convex hull of a set of points.
#'
#' @param pts A matrix with a point in each row.
#' @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 tol Tolerance on std. dev. if using PCA.
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#'
#' @return A list with \code{hull} equal a matrix with row indices of the vertices defining each
#'   facet in the hull and \code{pts} equal the input points (and dummy points) and columns:
#'   \code{pt}, true if a point in the original input; false if a dummy point (a point on a ray).
#'   \code{vtx}, TRUE if a vertex in the hull.
#' @export
#'
#' @examples
#' ## 1D
#' pts<-matrix(c(1,2,3), ncol = 1, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#'
#' ## 2D
#' pts<-matrix(c(1,1, 2,2), ncol = 2, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#' plotHull2D(pts, drawPoints = TRUE)
#' plotHull2D(pts, addRays = TRUE, drawPoints = TRUE)
#' pts<-matrix(c(1,1, 2,2, 0,1), ncol = 2, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts)
#' plotHull2D(pts, drawPoints = TRUE)
#' convexHull(pts, addRays = TRUE, direction = c(-1,1))
#'
#' ## 3D
#' pts<-matrix(c(1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a point
#' convexHull(pts)
#' pts<-matrix(c(0,0,0,1,1,1,2,2,2,3,3,3), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#' pts<-matrix(c(0,0,0,0,1,1,0,2,2,0,0,2), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts)
#' pts<-matrix(c(1,0,0,1,1,1,1,2,2,3,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts) # a polyhedron
#' pts<-matrix(c(1,1,1,2,2,1,2,1,1,1,1,2), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polytope (polyhedron)
#' convexHull(pts)
#'
#' ini3D(argsPlot3d = list(xlim = c(0,3), ylim = c(0,3), zlim = c(0,3)))
#' pts<-matrix(c(1,1,1,2,2,1,2,1,1,1,1,2), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' plotHull3D(pts, argsPolygon3d = list(color = "red"))
#' convexHull(pts)
#' finalize3D()
convexHull <- function(pts, addRays = FALSE, useRGLBBox = FALSE, direction = 1,
tol = mean(mean(abs(pts)))*sqrt(.Machine\$double.eps)*2,
m = apply(pts,2,min)-5, M = apply(pts,2,max)+5) {
pts <- .checkPts(pts)
p <- ncol(pts)
if (length(direction) != p) direction = rep(direction,p)
# print(set)
if (rgl::rgl.cur() > 0 & useRGLBBox) {
limits <- rgl::par3d()\$bbox
for (i in 1:dim(pts)) {
pt <- as.vector(pts[i,])
if (!(limits < pt && pt < limits &&
limits < pt && pt < limits &&
limits < pt && pt < limits)){
stop("The point is not in the interior of the current bounding box. Resize your axes.")
}
}
m <- c(limits, limits, limits)
M <- c(limits, limits, limits)
set <- addRays(pts, m, M, direction)
} else set <- addRays(pts, m, M, direction = direction)
} else {
set <- cbind(pts, pt = 1) # point in (1=original set, 0=artificial)
}
d <- dimFace(set[,1:p, drop = FALSE])
# cat("Object of dimension",d,"\n")

## note p can be > 3 !!
hull <- NULL
if (d==0) { # a point
hull <- 1
} else if (d==p) { # same dimension as space
if (p == 1) hull <- grDevices::chull(cbind(set[,1:p],0))
if (p == 2) hull <- grDevices::chull(set[,1:p])
if (p > 2) hull <- geometry::convhulln(set[,1:p], return.non.triangulated.facets = TRUE)
} else if (d==1) { # a line (d<p)
l <- nrow(set)
n <- p
comb <- t(utils::combn(n,2))
for (i in 1:dim(comb)) {  # simple projection down on each axis
pt <- unique(set[,comb[i,]])
if (l == dim(pt)) {
hull <- grDevices::chull(pt)
if (length(hull)==2) {
hull <- as.numeric(rownames(pt[hull,]))
break
}
}
}
} else if (d==2) { # a plane (d<p)
l <- nrow(set)
tmp <- cbind(set, id = 1:l)
comb <- t(utils::combn(p,2))
for (i in 1:dim(comb)) {  # simple projection down on each axis (keep the projection with highest number of vertices)
pt <- unique(set[,comb[i,]])
iL <- 0
if (l == dim(pt)) {
hull <- grDevices::chull(pt)
if (length(hull)>iL) {
iL <- length(hull)
hB <- rownames(pt[hull,])
}
}
}
hull <- tmp[hB,'id'] # return idx of set
} else { # d > 2 & d < p
warning(
"Using PCA to transform. Results might not be accurate!"
)
set1 <- stats::prcomp(set[,1:p])
# set1 <- stats::princomp(set[,1:p])
idx <- which(set1\$sdev > tol)
set1 <- stats::predict(set1)[,idx, drop = FALSE]
if (ncol(set1) == 1) hull <- grDevices::chull(cbind(set1,0))
if (ncol(set1) == 2) hull <- grDevices::chull(set1)
if (ncol(set1) > 2) hull <- geometry::convhulln(set1, return.non.triangulated.facets = TRUE)
}
if (is.vector(hull)) hull <- matrix(hull, nrow = 1)
## classify
idx <- unique(as.vector(hull))
idx <- idx[!is.na(idx)]
idx <- sort(idx)
set <- as.data.frame(set)
set\$vtx <- FALSE
set\$vtx[idx] <- TRUE
# remove unwarnted dummy points
set\$oldId <- 1:nrow(set)
tmp <- dplyr::filter(set, .data\$vtx | .data\$pt, .preserve = TRUE)
if (nrow(set) != nrow(tmp)) {
tmp\$newId <- 1:nrow(tmp)
mch <- dplyr::left_join(set, tmp, by = 'oldId') %>% dplyr::select(.data\$newId)
hull <- apply(hull, c(1,2), function(id) {if (is.na(id)) return(NA) else return(mch[id,1])} )
set <- dplyr::select(tmp, -.data\$oldId, -.data\$newId)
} else set <- dplyr::select(tmp, -.data\$oldId)
return(list(hull = hull, pts = set))
stop("Cannot find the vertices!")
}
```

## Try the gMOIP package in your browser

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

gMOIP documentation built on Aug. 23, 2021, 5:09 p.m.