R/mesh_set.R

In IRSN/DiceView: Methods for Visualization of Computer Experiments Design and Surrogate

#### Excursion sets ####

#' @title Search excursion set of nD function, sampled by a mesh
#' @param f Function to inverse at 'threshold'
#' @param threshold target value to inverse
#' @param sign focus at conservative for above (sign=1) or below (sign=-1) the threshold
#' @param intervals bounds to inverse in, each column contains min and max of each dimension
#' @param mesh.type "unif" or "seq" (default) or "LHS" to preform interval partition
#' @param mesh.sizes number of parts for mesh (duplicate for each dimension if using "seq")
#' @param maxerror_f maximal tolerance on f precision
#' @param ex_filter.tri boolean function to validate a geometry::tri as considered in excursion : 'any' or 'all'
#' @param ... parameters to forward to roots_mesh(...) call
#' @param vectorized boolean: is f already vectorized ? (default: FALSE) or if function: vectorized version of f.
#' @param tol the desired accuracy (convergence tolerance on f arg).
#' @importFrom geometry delaunayn
#' @export
#' @examples
#' # mesh_exsets(function(x) x, threshold=.51, sign=1, intervals=rbind(0,1),
#' #   maxerror_f=1E-2,tol=1E-2) # for faster testing
#' # mesh_exsets(function(x) x, threshold=.50000001, sign=1, intervals=rbind(0,1),
#' #   maxerror_f=1E-2,tol=1E-2) # for faster testing
#' # mesh_exsets(function(x) sum(x), threshold=.51,sign=1, intervals=cbind(rbind(0,1),rbind(0,1)),
#' #   maxerror_f=1E-2,tol=1E-2) # for faster testing
#' # mesh_exsets(sin,threshold=0,sign="sup",interval=c(pi/2,5*pi/2),
#' #   maxerror_f=1E-2,tol=1E-2) # for faster testing
#'
#' if (identical(Sys.getenv("NOT_CRAN"), "true")) { # too long for CRAN on Windows
#'
#'   e = mesh_exsets(function(x) (0.25+x[1])^2+(0.5+x[2])^2 ,
#'                 threshold =0.25,sign=-1, intervals=matrix(c(-1,1,-1,1),nrow=2),
#'                 maxerror_f=1E-2,tol=1E-2) # for faster testing
#'
#'   plot(e$p,xlim=c(-1,1),ylim=c(-1,1));
#'   apply(e$tri,1,function(tri) polygon(e$p[tri,],col=rgb(.4,.4,.4,.4)))
#'
#'   if (requireNamespace("rgl")) {
#'     e = mesh_exsets(function(x) (0.5+x[1])^2+(-0.5+x[2])^2+(0.+x[3])^2,
#'                   threshold = .25,sign=-1, mesh.type="unif",
#'                   intervals=matrix(c(-1,1,-1,1,-1,1),nrow=2),
#'                   maxerror_f=1E-2,tol=1E-2) # for faster testing
#'
#'     rgl::plot3d(e$p,xlim=c(-1,1),ylim=c(-1,1),zlim=c(-1,1));
#'     apply(e$tri,1,function(tri)rgl::lines3d(e$p[tri,]))
#'   }
#' }
mesh_exsets = function (f, vectorized = FALSE, threshold, sign, intervals,
          mesh.type = "seq", mesh.sizes = 11, maxerror_f = 1e-09, tol = .Machine$double.eps^0.25,
          ex_filter.tri = all, ...) {
    if (sign == "lower" || sign == -1 || sign == "inf" || sign == "<" || isFALSE(sign))
        return(
            mesh_exsets(f = function(...) { -f(...) }, vectorized = vectorized, threshold = -threshold, sign = 1,
            intervals = intervals, mesh.type = mesh.type, mesh.sizes = mesh.sizes,
            maxerror_f = maxerror_f, tol = tol, ...))
    if (sign != "upper" && sign != 1 && sign != "sup" && sign != ">" && !isTRUE(sign))
        stop("unknown sign: '", sign, "'")

    if (is.matrix(mesh.type))
        d = ncol(mesh.type)
    else if (is.matrix(intervals))
        d = ncol(intervals)
    else if (is.array(intervals))
        d = 1
    else
        stop("cannot identify dim of mesh")

    if (isTRUE(vectorized) && is.function(f))
        f_vec = function(x, ...) f(x, ...)
    else if (is.function(vectorized)) {
        f_vec = function(x, ...) vectorized(x, ...)
        if (is.null(f)) f = f_vec
    } else if (isFALSE(vectorized) && !is.null(f))
        f_vec = Vectorize.function(f, dim = d)
    else
        stop("Cannot decide how to vectorize f")

    f_0 <- function(...) return(f(...) - threshold)
    f_vec_0 <- function(...) return(f_vec(...) - threshold)
    r <- roots_mesh(f = f_0, vectorized = f_vec_0, intervals = intervals,
                    mesh.type = mesh.type, mesh.sizes = mesh.sizes, maxerror_f = maxerror_f,
                    tol = tol, ...)

    if (all(is.na(r)))
        all_points = attr(r, "mesh")$p
    else
        all_points = rbind(attr(r, "mesh")$p, r)

    if (maxerror_f != 0) { # try add points on each side of the frontier
        f_inf <- function(...) return(f(...) - (threshold - maxerror_f))
        f_vec_inf <- function(...) return(f_vec(...) - (threshold - maxerror_f))
        r_inf <- roots_mesh(f = f_inf, vectorized = f_vec_inf, intervals = intervals,
                            mesh.type = mesh.type, mesh.sizes = mesh.sizes, maxerror_f = maxerror_f,
                            tol = tol, ...)
        if (!all(is.na(r_inf)))
            all_points = rbind(all_points, r_inf)

        f_sup <- function(...) return(f(...) - (threshold + maxerror_f))
        f_vec_sup <- function(...) return(f_vec(...) - (threshold + maxerror_f))
        r_sup <- roots_mesh(f = f_sup, vectorized = f_vec_sup, intervals = intervals,
                            mesh.type = mesh.type, mesh.sizes = mesh.sizes, maxerror_f = maxerror_f,
                            tol = tol, ...)
        if (!all(is.na(r_sup)))
            all_points = rbind(all_points, r_sup)
    }

    new_mesh <- mesh(intervals, mesh.type = all_points)

    new_mesh$y = f_vec(new_mesh$p, ...)
    I = which(apply(new_mesh$tri, 1,
              function(i) ex_filter.tri(new_mesh$y[i] >= (threshold - 2 * maxerror_f))))
    colnames(new_mesh$p) <- colnames(intervals)
    return(list(p = new_mesh$p, tri = new_mesh$tri[I, ], areas = new_mesh$areas[I],
                neighbours = new_mesh$neighbours[I]))
}

#### Plot meshes ####

#' @title Plot a one dimensional mesh
#' @param mesh 1-dimensional mesh to draw
#' @param y ordinate value where to draw the mesh
#' @param color.nodes color of the mesh nodes
#' @param color.mesh color of the mesh elements
#' @param alpha transparency of the mesh elements & nodes
#' @param ... optional arguments passed to plot function
#' @import graphics
#' @export
#' @examples
#' plot_mesh(mesh_exsets(function(x) x, threshold=.51, sign=1, intervals=rbind(0,1)))
#' plot_mesh(mesh_exsets(function(x) (x-.5)^2, threshold=.1, sign=-1, intervals=rbind(0,1)))
plot_mesh = function(mesh,y=0,color.nodes='black',color.mesh='darkgray',alpha=0.4, ...){
    nodes.rgb=col2rgb(color.nodes)/255
    p = plot(x=mesh$p,y=rep(y,length(mesh$p)),ylab="",col=rgb(nodes.rgb[1,],nodes.rgb[2,],nodes.rgb[3,],alpha),...)
    mesh.rgb=col2rgb(color.mesh)/255
    apply(mesh$tri,1,
          function(tri) {
              lines(mesh$p[tri,],y=rep(0,length(tri)),col=rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,],alpha))
              polygon(x=c(mesh$p[tri,],rev(mesh$p[tri,])),y=c(-1,-1,1,1),border = NA,col=rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,],alpha))
          })
    invisible(p)
}

#' @title Plot a two dimensional mesh
#' @param mesh 2-dimensional mesh to draw
#' @param color.nodes color of the mesh nodes
#' @param color.mesh color of the mesh elements
#' @param alpha transparency of the mesh elements & nodes
#' @param ... optional arguments passed to plot function
#' @import graphics
#' @export
#' @examples
#' plot2d_mesh(mesh_exsets(f = function(x) sin(pi*x[1])*sin(pi*x[2]),
#'                         threshold=0,sign=1, mesh.type="unif",mesh.size=11,
#'                         intervals = matrix(c(1/2,5/2,1/2,5/2),nrow=2)))
plot2d_mesh = function(mesh,color.nodes='black',color.mesh='darkgray',alpha=0.4,...){
    nodes.rgb=col2rgb(color.nodes)/255
    p = plot(mesh$p,col=rgb(nodes.rgb[1,],nodes.rgb[2,],nodes.rgb[3,],alpha),...)
    mesh.rgb=col2rgb(color.mesh)/255
    apply(mesh$tri,1,
          function(tri)
              polygon(mesh$p[c(tri,tri[1]),],border = rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,],alpha), col = rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,],alpha)))
    # sapply(1:nrow(mesh$tri),function(i) {xy=colMeans(mesh$p[mesh$tri[i,],]);text(x = xy[1],y=xy[2],paste0(i))})
    # apply(mesh$tri,1,function(tri) points(mesh$p[tri,], col = rgb(0,0,0,1),pch=20))
    invisible(p)
}

#' @title Plot a three dimensional mesh
#' @param mesh 3-dimensional mesh to draw
#' @param engine3d 3d framework to use: 'rgl' if installed or 'scatterplot3d' (default)
#' @param color.nodes color of the mesh nodes
#' @param color.mesh color of the mesh elements
#' @param alpha transparency of the mesh elements & nodes
#' @param ... optional arguments passed to plot function
#' @export
#' @examples
#' if (identical(Sys.getenv("NOT_CRAN"), "true")) { # too long for CRAN on Windows
#'
#'   plot3d_mesh(mesh_exsets(function(x) (0.5+x[1])^2+(-0.5+x[2])^2+(0.+x[3])^2,
#'                           threshold = .25,sign=-1, mesh.type="unif",
#'                           maxerror_f=1E-2,tol=1E-2, # faster display
#'                           intervals=matrix(c(-1,1,-1,1,-1,1),nrow=2)),
#'                           engine3d='scatterplot3d')
#'
#'   if (requireNamespace("rgl")) {
#'     plot3d_mesh(mesh_exsets(function(x) (0.5+x[1])^2+(-0.5+x[2])^2+(0.+x[3])^2,
#'                             threshold = .25,sign=-1, mesh.type="unif",
#'                             maxerror_f=1E-2,tol=1E-2, # faster display
#'                             intervals=matrix(c(-1,1,-1,1,-1,1),nrow=2)),engine3d='rgl')
#'   }
#' }
plot3d_mesh = function(mesh,engine3d=NULL,color.nodes='black',color.mesh='darkgray',alpha=0.4,...){
    nodes.rgb=col2rgb(color.nodes)/255
    package = load3d(engine3d)
    if (is.null(package)) return()
    p3d = plot3d(mesh$p,col=rgb(nodes.rgb[1,],nodes.rgb[2,],nodes.rgb[3,],alpha), package=package,...)
    mesh.rgb=col2rgb(color.mesh)/255
    apply(mesh$tri,1,function(tri) {
        quads3d(mesh$p[tri,],col=rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,]),alpha=alpha/10, package=package)
        quads3d(mesh$p[tri,][c(4,3,2,1),],col=rgb(mesh.rgb[1,],mesh.rgb[2,],mesh.rgb[3,]),alpha=alpha/10, package=package)
        # triangles3d(mesh$p[tri,][-1,],col=color,alpha=0.05, package=package)
        # triangles3d(mesh$p[tri,][-2,],col=color,alpha=0.05, package=package)
        # triangles3d(mesh$p[tri,][-3,],col=color,alpha=0.05, package=package)
        # triangles3d(mesh$p[tri,][-4,],col=color,alpha=0.05, package=package)
    }) #rgl::lines3d(mesh$p[t(combn(tri,2)),],col=color))
    invisible(p3d)
}


#### Utility functions ####

#' @title Generalize expand.grid() for multi-columns data. Build all combinations of lines from X1 and X2. Each line may hold multiple columns.
#' @param X1 variable values, possibly with many columns
#' @param X2 variable values, possibly with many columns
#' combn.design(matrix(c(10,20),ncol=1),matrix(c(1,2,3,4,5,6),ncol=2))
#' combn.design(matrix(c(10,20,30,40),ncol=2),matrix(c(1,2,3,4,5,6),ncol=2))
combn.design <- function(X1,X2) {
    n1 = nrow(X1)
    n2 = nrow(X2)
    n = expand.grid(1:n1,1:n2)
    C = matrix(NA,ncol=ncol(X1)+ncol(X2),nrow=n1*n2)
    for (i in 1:(n1*n2)) {
        C[i,1:ncol(X1)] = X1[n[i,1],]
        C[i,ncol(X1)+(1:ncol(X2))] = X2[n[i,2],]
    }
    C
}

#' @title Test if points are in a hull
#' @param x points to test
#' @param p points defining the hull
#' @param h hull itself (built from p if given as NULL (default))
#' @export
#' @importFrom geometry convhulln
#' @importFrom geometry inhulln
#' @examples
#' is_in.p(x=-0.5,p=matrix(c(0,1),ncol=1))
#' is_in.p(x=0.5,p=matrix(c(0,1),ncol=1))
#' is_in.p(x=matrix(-.5,ncol=2,nrow=1),p=matrix(c(0,0,1,1,0,0),ncol=2))
#' is_in.p(x=matrix(.25,ncol=2,nrow=1),p=matrix(c(0,0,1,1,0,0),ncol=2))
#' is_in.p(x=matrix(-.5,ncol=3,nrow=1),p=matrix(c(0,0,0,1,0,0,0,1,0,0,0,1),ncol=3,byrow = TRUE))
#' is_in.p(x=matrix(.25,ncol=3,nrow=1),p=matrix(c(0,0,0,1,0,0,0,1,0,0,0,1),ncol=3,byrow = TRUE))
is_in.p = function(x,p,h=NULL) {
    if (is.null(h)) { # use p points, build h hull
        if (ncol(p)==1) return(x>min(p) && x<max(p))
        if (!is.matrix(x)) x = matrix(x,ncol=ncol(p))
        h = geometry::convhulln(p)
    }
    geometry::inhulln(ch = h,p=x)
}

#' Checks if some point belongs to a given mesh
#' @param x point to check
#' @param mesh mesh identifying the set which X may belong
#' @import stats
#' @export
#' @examples
#' is_in.mesh(-0.5,mesh=geometry::delaunayn(matrix(c(0,1),ncol=1),output.options =TRUE))
#' is_in.mesh(0.5,mesh=geometry::delaunayn(matrix(c(0,1),ncol=1),output.options =TRUE))
#'
#' x =matrix(-.5,ncol=2,nrow=1)
#' is_in.mesh(x,mesh=geometry::delaunayn(matrix(c(0,0,1,1,0,0),ncol=2),output.options =TRUE))
#'
#' x =matrix(.5,ncol=2,nrow=1)
#' is_in.mesh(x,mesh=geometry::delaunayn(matrix(c(0,0,1,1,0,0),ncol=2),output.options =TRUE))
is_in.mesh = function(x,mesh) {
    for (i in 1:nrow(mesh$tri)) {
        # if (mesh$areas[i]>.Machine$double.eps) # otherwise, flat triangle NOT RELIABLE
            if (all(dist(mesh$p[mesh$tri[i,],, drop = FALSE])>.Machine$double.eps)) {# not too close points
                #if (qr(mesh$p[mesh$tri[i,],, drop = FALSE])$rank == ncol(mesh$p)) # check this is not a degenerated (flat) triangle
                isin = FALSE
                try(isin <- is_in.p(x, mesh$p[mesh$tri[i,],, drop = FALSE]),silent = T) # to avoid coplanar errors in qhull
                if (isin) return(TRUE)
            } # else warning("Too flat element (",i,"): ",capture.output(print(mesh$p[mesh$tri[i,],, drop = FALSE])))
        # else warning("Too small element (",i,"): area=",mesh$areas[i])
        }
    return(FALSE)
}

#' @title Checks if some points belong to a given mesh
#' @param X points to check
#' @param mesh mesh identifying the set which X may belong
#' @export
#' @examples
#' X = matrix(runif(100),ncol=2);
#' inside = are_in.mesh(X,mesh=geometry::delaunayn(matrix(c(0,0,1,1,0,0),ncol=2),output.options =TRUE))
#' print(inside)
#' plot(X,col=rgb(1-inside,0,0+inside))
are_in.mesh = function(X,mesh) {
    #apply(X,1,is_in.mesh,mesh)
    #unlist(foreach(i = 1:nrow(X)) %dopar% {is_in.mesh(X[i,],mesh)})
    X.list = lapply(seq_len(nrow(X)), function(i) X[i,])
    array(unlist(parallel::mclapply(X.list,is_in.mesh,mesh)))
}

#' @title Extract points of mesh which belong to the mesh triangulation (may not contain all points)
#' @param mesh mesh (list(p,tri,...) from geometry)
#' @return points coordinates inside the mesh triangulation
#' @export
points_in.mesh = function(mesh) {
    if (is.null(mesh)) return(NULL)
    mesh$p[unique(array(mesh$tri)),, drop = FALSE]
}
#' @title Extract points of mesh which do not belong to the mesh triangulation (may not contain all points)
#' @param mesh (list(p,tri,...) from geometry)
#' @return points coordinates outside the mesh triangulation
#' @export
points_out.mesh = function(mesh) {
    if (is.null(mesh)) return(NULL)
    mesh$p[-unique(array(mesh$tri)),, drop = FALSE]
}

#' @title Checks if a mesh is valid
#' @param x mesh to check
#' @export
#' @return TRUE if mesh is valid
is.mesh = function(x) {
    is.list(x) && all(c("p","tri") %in% names(x))
}

#### Build mesh (on hypercube) ####

#' @title Builds a mesh from a design aor set of points
#' @param intervals bounds to inverse in, each column contains min and max of each dimension
#' @param mesh.type function or "unif" or "seq" (default) or "LHS" to preform interval partition
#' @param mesh.sizes number of parts for mesh (duplicate for each dimension if using "seq")
#' @return delaunay mesh (list(p,tri,...) from geometry)
#' @export
#' @examples
#' mesh = mesh(intervals=matrix(c(0,1,0,1),ncol=2),mesh.type="unif",mesh.sizes=10)
#' plot2d_mesh(mesh)
mesh = function(intervals, mesh.type = "seq", mesh.sizes = 11) {
    # setup bounds & dim
    if (is.matrix(intervals)) {
        lowers = apply(intervals, 2, min)
        uppers = apply(intervals, 2, max)
        d = ncol(intervals)
    } else if (is.null(intervals) && is.matrix(mesh.type)) {
        lowers = apply(mesh.type, 2, min)
        uppers = apply(mesh.type, 2, max)
        d = ncol(mesh.type)
    } else
        d = 1

    if (length(mesh.sizes) != d)
        mesh.size = rep(mesh.sizes, d)

    if (is.matrix(mesh.type)) {
        ridge.points = mesh.type
    } else if (is.function(mesh.type)) {
        ridge.points = mesh.type(prod(mesh.size), lowers, uppers)
    } else if (mesh.type == "seq") {
        ridge.points = matrix(seq(f = lowers[1], to = uppers[1],
                                  length.out = mesh.size[1]), ncol = 1)
        for (i in 2:d) {
            ridge.points = combn.design(ridge.points, matrix(seq(f = lowers[i], to = uppers[i], length.out = mesh.size[i]), ncol = 1))
        }
    } else if (mesh.type == "unif") {
        ridge.points = matrix(runif(n = d * prod(mesh.size),
                                    min = 0, max = 1), ncol = d)
        ridge.points = matrix(lowers, ncol = d, nrow = nrow(ridge.points),
                              byrow = TRUE) +
            ridge.points * (matrix(uppers, ncol = d, nrow = nrow(ridge.points), byrow = TRUE) -
                                matrix(lowers, ncol = d, nrow = nrow(ridge.points), byrow = TRUE))

    } else if (mesh.type == "LHS") {
        ridge.points = DiceDesign::lhsDesign(prod(mesh.size), dimension = d)$design
        ridge.points = matrix(lowers, ncol = d, nrow = nrow(ridge.points), byrow = TRUE) +
            ridge.points * (matrix(uppers, ncol = d, nrow = nrow(ridge.points), byrow = TRUE) -
                                matrix(lowers, ncol = d, nrow = nrow(ridge.points), byrow = TRUE))

    } else stop("unsupported mesh setup : ", mesh.type)

    # add bounds
    b = cbind(lowers[1], uppers[1])
    if (d>1) for (id in 1:(d - 1)) { # efficient way for factorial design
        b = rbind(cbind(b, lowers[id + 1]), cbind(b, uppers[id + 1]))
    }
    for (i in 1:nrow(b)) if (min_dist(b[i, ], ridge.points) > .Machine$double.eps)
        ridge.points = rbind(ridge.points, b[i, ])

    return( geometry::delaunayn(ridge.points, output.options = TRUE) )
}