R/mesh_set.R

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

#### Generalization of root finding ####

#' @title One Dimensional Root (Zero) Finding
#' @description Search one root with given precision (on y). Iterate over uniroot as long as necessary.
#' @param f the function for which the root is sought.
#' @param lower the lower end point of the interval to be searched.
#' @param upper the upper end point of the interval to be searched.
#' @param maxerror_f the maximum error on f evaluation (iterates over uniroot to converge).
#' @param f_lower the same as f(lower).
#' @param f_upper the same as f(upper).
#' @param tol the desired accuracy (convergence tolerance on f arg).
#' @param convexity the learned convexity factor of the function, used to reduce the boundaries for uniroot.
#' @param ... additional named or unnamed arguments to be passed to f.
#' @import stats
#' @export
#' @author Yann Richet, IRSN
#' @examples
#' f=function(x) {cat("f");1-exp(x)}; f(root(f,lower=-1,upper=2))
#' f=function(x) {cat("f");exp(x)-1}; f(root(f,lower=-1,upper=2))
#'
#' .f = function(x) 1-exp(1*x)
#' f=function(x) {cat("f");y=.f(x);points(x,y,pch=20,col=rgb(0,0,0,.2));y}
#' plot(.f,xlim=c(-1,2)); f(root(f,lower=-1,upper=2))
#'
#' .f = function(x) exp(10*x)-1
#' f=function(x) {cat("f");y=.f(x);points(x,y,pch=20);y}
#' plot(.f,xlim=c(-1,2)); f(root(f,lower=-1,upper=2))
#'
#' .f = function(x) exp(100*x)-1
#' f=function(x) {cat("f");y=.f(x);points(x,y,pch=20);y}
#' plot(.f,xlim=c(-1,2)); f(root(f,lower=-1,upper=2))
#'
#' f=function(x) {cat("f");exp(100*x)-1}; f(root(f,lower=-1,upper=2))
root = function(f,lower,upper,maxerror_f=1E-7, f_lower = f(lower, ...),f_upper = f(upper, ...), tol = .Machine$double.eps^0.25, convexity=0, ...) {
    if (upper < lower)
        return(root(f,lower=upper,upper=lower,maxerror_f=maxerror_f, f_lower = f(upper, ...), f_upper=f(lower, ...), tol = tol, convexity=convexity, ...))

    r = NULL
    try(
        r <- uniroot(f,lower=lower,upper=upper,f.lower = f_lower,f.upper = f_upper,tol=tol,...)
        ,silent=T)
    if (is.null(r)) {
        warning(paste0("No root found in [",lower,",",upper,"] -> [",f_lower,",",f_upper,"]"))
        return(NULL)
    }

    r_root = r$root
    f_root = f(r_root, ...)
    err = abs(f_root)/maxerror_f

    # cat(paste0("[ ",lower," | ",r_root," | ",upper," ] -> ^",exp(convexity)," -> [ ",f_lower," | ",f_root," | ",f_upper," ]\n"))

    if (err>1) {
        if (f_lower*f_root<0) {
            # cat("<")
            x0 = r_root - (r_root-lower) * max(1E-3,(f_root/(f_root-f_lower))^exp(convexity))
            f_x0 = f(x0, ...)
            # cat(paste0(" x0:",x0," -> ",f_x0))
            if (f_x0*f_root<0){
                # cat("+")
                return(root(f, lower=x0,upper=r_root,maxerror_f=maxerror_f,f_lower=f_x0,f_upper=f_root, convexity=convexity-f_x0/(f_root-f_lower) ,tol=tol,...))
            } else {
                # cat("-")
                return(root(f, lower=lower,upper=x0,maxerror_f=maxerror_f,f_lower=f_lower,f_upper=f_x0, convexity=convexity-f_x0/(f_root-f_lower) ,tol=tol,...))
            }
        } else { #(f_upper*f_root<0)
            # cat(">")
            x0 = r_root + (upper-r_root) * max(1E-3,(f_root/(f_root-f_upper))^exp(convexity))
            f_x0 = f(x0, ...)
            # cat(paste0(" x0:",x0," -> ",f_x0))
            if (f_x0*f_root<0){
                # cat("-")
                return(root(f, lower=r_root,upper=x0,maxerror_f=maxerror_f,f_lower=f_root,f_upper=f_x0, convexity=convexity-f_x0/(f_upper-f_root) ,tol=tol,...))
            }else{
                # cat("+")
                return(root(f, lower=x0,upper=upper,maxerror_f=maxerror_f,f_lower=f_x0,f_upper=f_upper, convexity=convexity-f_x0/(f_upper-f_root) ,tol=tol,...))
            }
        }
    } else
        r$root
}

#' @title One Dimensional Multiple Roots (Zero) Finding
#' @description Search multiple roots of 1D function, sampled/splitted by a (1D) mesh
#' @param f Function to find roots
#' @param interval bounds to inverse in
#' @param split.size number of parts to perform uniroot inside
#' @param split function or "unif" or "seq" (default) to preform interval partition
#' @param maxerror_f the maximum error on f evaluation (iterates over uniroot to converge).
#' @param tol the desired accuracy (convergence tolerance on f arg).
#' @param ... additional named or unnamed arguments to be passed to f.
#' @return array of x, so f(x)=target
#' @import stats
#' @export
#' @examples
#' roots(sin,interval=c(pi/2,5*pi/2))
#' roots(sin,interval=c(pi/2,1.5*pi/2))
#'
#' f=function(x)exp(x)-1;
#' f(roots(f,interval=c(-1,2)))
#'
#' f=function(x)exp(1000*x)-1;
#' f(roots(f,interval=c(-1,2)))
roots = function(f,interval, maxerror_f=1E-7,split="seq",split.size=11,tol=.Machine$double.eps^0.25,...) {
    lower = min(interval)
    upper = max(interval)

    if (is.function(split))
        intervals = split(min=lower,max=upper,n=split.size)
    else if (is.numeric(split))
        intervals = split
    else if (split=="seq")
        intervals = seq(from=lower,to=upper,length.out=split.size)
    else if (split=="unif")
        intervals = runif(min=lower,max=upper,n=split.size)
    else stop("unsupported mesh: function(min,max,n)/array/'seq'/'fun': type is ",typeof(split)," : ",paste0(split,collapse=";",sep=","))

    # roots = foreach (i = 1:(length(intervals)-1)) %dopar% {
    #     lower.i = intervals[i]
    #     upper.i = intervals[i+1]
    #     r = NULL
    #     try({r = root(f,lower=lower.i,upper=upper.i,maxerror_f = maxerror_f,...)},silent = T)
    #     r
    # }
    I.list = lapply(seq_len(length(intervals)-1), function(i) c(intervals[i],intervals[i+1]))
    roots = parallel::mclapply(I.list,
                               function(I) {
                                   r = NULL;
                                   try({r = root(f,lower=I[1],upper=I[2],maxerror_f = maxerror_f,tol=tol,...)},silent = F);
                                   r
                                })
    roots = unlist(Filter(Negate(is.null), roots))

    if (length(roots)>0)
        r = matrix(unlist(roots),ncol=1)
    else
        r = NA
    attr(r,"mesh") <- list(p=matrix(intervals,ncol=1))
    return(r)
}

#' @title Minimal distance between one point to many points
#' @param x one point
#' @param X matrix of points (same number of columns than x)
#' @param norm normalization vecor of distance (same number of columns than x)
#' @return minimal distance
#' @export
#' @examples
#' min_dist(runif(3),matrix(runif(30),ncol=3))
min_dist <- function (x, X, norm=rep(1,ncol(X))){
    return(min(sqrt(rowSums(((X - matrix(x, nrow = nrow(X), ncol = ncol(X), byrow = TRUE))/norm)^2))))
}

#' @title Multi Dimensional Multiple Roots (Zero) Finding, sampled by a mesh
#' @param f Function (one or more dimensions) to find roots of
#' @param intervals bounds to inverse in, each column contains min and max of each dimension
#' @param mesh function or "unif" or "seq" (default) to preform interval partition
#' @param mesh.sizes number of parts for mesh (duplicate for each dimension if using "seq")
#' @param vectorized is f already vectorized ? (default: no)
#' @param maxerror_f the maximum error on f evaluation (iterates over uniroot to converge).
#' @param tol the desired accuracy (convergence tolerance on f arg).
#' @param ... Other args for f
#' @import stats
#' @importFrom DiceDesign lhsDesign
#' @importFrom parallel mclapply
#' @export
#' @return matrix of x, so f(x)=0
#' @examples
#' mesh_roots(function(x) x-.51, intervals=rbind(0,1))
#' mesh_roots(function(x) sum(x)-.51, intervals=cbind(rbind(0,1),rbind(0,1)))
#' mesh_roots(sin,intervals=c(pi/2,5*pi/2))
#' mesh_roots(f = function(x) sin(pi*x[1])*sin(pi*x[2]),
#'            intervals = matrix(c(1/2,5/2,1/2,5/2),nrow=2))
#'
#' r = mesh_roots(function(x) (0.25+x[1])^2+(0.5+x[2])^2-.25,
#'                intervals=matrix(c(-1,1,-1,1),nrow=2))
#' plot(r,xlim=c(-1,1),ylim=c(-1,1))
#'
#' r = mesh_roots(function(x) (0.5+x[1])^2+(-0.5+x[2])^2+(0.+x[3])^2- .25,
#'                mesh.sizes = 11,
#'                intervals=matrix(c(-1,1,-1,1,-1,1),nrow=2))
#' scatterplot3d::scatterplot3d(r,xlim=c(-1,1),ylim=c(-1,1),zlim=c(-1,1))
#'
#' mesh_roots(function(x)exp(x)-1,intervals=c(-1,2))
#' mesh_roots(function(x)exp(1000*x)-1,intervals=c(-1,2))
mesh_roots = function(f,vectorized=FALSE,intervals, mesh="seq",mesh.sizes=11,maxerror_f=1E-7,tol=.Machine$double.eps^0.25,...) {
    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)) {
        lowers = apply(mesh,2,min)
        uppers = apply(mesh,2,max)
        d = ncol(mesh)
    } else { d=1 }

    if (d==1)
        return(roots(f,interval = intervals,split = mesh,split.size = mesh.sizes,maxerror_f=maxerror_f,tol=tol,...))

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

    # Create mesh points matrix
    if (is.matrix(mesh)) {
        ridge.points = mesh
    } else if (is.function(mesh)) {
        ridge.points = mesh(prod(mesh.size), lowers, uppers)
    } else if (mesh == "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 == "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 == "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 : ",mesh)

    # Add bounds if needed
    b = c(lowers[1], uppers[1])
    for (id in 1:(d-1)) {
        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,])

    # y = apply(ridge.points,1,f,...)
    # plot(ridge.points,col=rgb(y>0,0,y<0))

    simplexes <- geometry::delaunayn(ridge.points,output.options = TRUE)

    # Now find ridges where target is bounded, and then root it
    if (isTRUE(vectorized))
        f_vec=function(x,...) f(x,...)
    else if (is.function(vectorized))
        f_vec=function(x,...) vectorized(x,...)
    else
        f_vec=Vectorize.function(f,dim=ncol(ridge.points))
    ridge.y = f_vec(ridge.points,...)

    ############### Benchmark: find ridges that crosses zero ###################
    # library(microbenchmark)
    #
    # f = function(x) exp(-x^2)-0.5
    # Npoints = 1000
    # ridge.points = matrix(runif(2*Npoints,-1,1),ncol=2)
    # simplexes <- geometry::delaunayn(ridge.points,output.options = TRUE)
    # ridge.y = f(simplexes$p)
    #
    # tcombn2 = function(x) (utils::combn(x,2))
    # tcombn2_arrangements = function(x) t(arrangements::combinations(x,2))
    # cbind.list = function(l)
    #     matrix(unlist(l),ncol=length(l),byrow=F) #do.call(cbind,l)
    # tcombn2_my = function(x){
    #     l=length(x)
    #     #C = matrix(NA,nrow=2,factorial(l-1))
    #     #i = 1
    #     #for (j in 1:(l-1)) {
    #     #    for (k in (j+1):l) {
    #     #        C[,i] = x[c(j,k)]
    #     #        i = i+1
    #     #    }
    #     #}
    #     #C
    #     do.call(cbind,lapply(1:(l-1),
    #                function(j) cbind.list(lapply((j+1):l,
    #                                       function(k) x[c(j,k)]))))
    # }
    #
    # microbenchmark({tcombn2(1:5)})
    # microbenchmark({tcombn2_arrangements(1:5)})
    # microbenchmark({tcombn2_my(1:5)})
    #
    # tcombn2 = function(x) t(arrangements::combinations(x,2)) # faster !
    #
    # # for x for (baseline)
    # microbenchmark({
    #   ridges = NULL
    #   for (i in 1:nrow(simplexes$tri)) {
    #       #if (!all(ridge.y[simplexes$tri[i,]]<0)) {
    #       more_ridges = tcombn2(simplexes$tri[i,])
    #       for (j in 1:nrow(more_ridges))
    #           if (ridge.y[more_ridges[j,1]]*ridge.y[more_ridges[j,2]]<0)
    #               ridges = rbind(ridges,more_ridges[j,])
    #      #}
    #   }
    # },times=10)
    #
    # # lapply x for
    # microbenchmark({
    #     ridges = lapply(X=as.list(as.data.frame(t(simplexes$tri))),FUN=function(tri_i) {
    #         #if (all(ridge.y[tri_i]<0)) return()
    #         all_ridges = tcombn2(simplexes$tri[i,])
    #         more_ridges=NULL
    #         for (j in 1:nrow(all_ridges))
    #             if (ridge.y[all_ridges[j,1]]*ridge.y[all_ridges[j,2]]<0) # not same sign
    #                 more_ridges = rbind(more_ridges,all_ridges[j,])
    #         return(more_ridges)
    #     })
    #     ridges = do.call(rbind,ridges)
    # },times=10)
    #
    # # lapply x apply
    # microbenchmark({
    #   ridges = lapply(X=as.list(as.data.frame(t(simplexes$tri))),FUN=function(tri_i) {
    #       #if (all(ridge.y[tri_i]<0)) return()
    #       all_ridges = tcombn2(simplexes$tri[i,])
    #       #yprod_all_ridges = apply(all_ridges,1,function(r) ridge.y[r[1]] * ridge.y[r[2]])
    #       #yprod_all_ridges = apply(all_ridges,1,function(r) prod(ridge.y[r]))
    #       yprod_all_ridges = all_ridges[
    #           unlist(lapply(1:nrow(all_ridges),
    #                  function(i)
    #                      if (ridge.y[all_ridges[i,1]] * ridge.y[all_ridges[i,2]]<0)
    #                          return(i)
    #                         else return(NULL))),]
    #       return(all_ridges[which(yprod_all_ridges<0),])
    #   })
    #   ridges = do.call(rbind,ridges)
    # },times=10)
    #
    ############################################################################

    if (requireNamespace("arrangements"))
        tcombn2 = function(x) arrangements::combinations(x,2)
    else
        tcombn2 = function(x) t(utils::combn(x,2))

    #Get all ridges where we will do uniroot later
    ridges = NULL
    for (i in 1:nrow(simplexes$tri)) {
        #if (!all(ridge.y[simplexes$tri[i,]]<0)) {
        more_ridges = tcombn2(simplexes$tri[i,])
        for (j in 1:nrow(more_ridges))
            if (ridge.y[more_ridges[j,1]]*ridge.y[more_ridges[j,2]]<0)
                ridges = rbind(ridges,more_ridges[j,])
       #}
    }

    if (is.null(ridges)) {
        r = NA
        if (!is.null(ridge.points)) {
            attr(r,"mesh") <- simplexes
        }
        return(r)
    }

    ridges = unique(as.matrix(ridges))

    # for (i in 1:nrow(ridges))
    #     lines(x=rbind(ridge.points[ridges[i,1],],ridge.points[ridges[i,2],]))


    # r = foreach(i.r = 1:nrow(ridges),.combine = rbind,.verbose = FALSE,.errorhandling='stop') %dopar% { # No error expected, so stop if occurs
    #     X = ridge.points[ridges[i.r,],,drop=FALSE]
    #     #y=f_vec(X,...)
    #     y = ridge.y[ridges[i.r,]]
    #     if (isFALSE((all(y>0) || all(y<0)))) {
    #         f.r = function(alpha,...)
    #             as.numeric(f(alpha*X[1,] + (1-alpha)*X[2,],...))
    #         alpha_i = root(f.r,lower=0,upper=1, maxerror_f=maxerror_f,...)
    #         return(alpha_i*X[1,] + (1-alpha_i)*X[2,])
    #     } else return(NULL)
    # }

    # r.list = lapply(seq_len(nrow(ridges)), function(i.r) ridges[i.r,])
    # r = parallel::mclapply(r.list,
    #                        function(ridge_i) {
    #                             X = ridge.points[ridge_i,,drop=FALSE]
    #                             y = ridge.y[ridge_i]
    #                             if (isFALSE((all(y>0) || all(y<0)))) {
    #                                 f.r = function(alpha,...)
    #                                     as.numeric(f(alpha*X[1,] + (1-alpha)*X[2,],...))
    #                                 alpha_i = root(f.r,lower=0,upper=1, maxerror_f=maxerror_f,tol=tol,...)
    #                                 return(alpha_i*X[1,] + (1-alpha_i)*X[2,])
    #                             } else return(NULL)
    #                         })
    # r = do.call("rbind",r)

    r = Apply.function(X=ridges,MARGIN=1,FUN=function(ridge_i) {
        X = ridge.points[ridge_i,,drop=FALSE]
        y = ridge.y[ridge_i]
        if (isFALSE((all(y>0) || all(y<0)))) {
            f.r = function(alpha,...)
                as.numeric(f(alpha*X[1,] + (1-alpha)*X[2,],...))
            alpha_i = root(f.r,lower=0,upper=1, maxerror_f=maxerror_f,tol=tol,...)
            return(alpha_i*X[1,] + (1-alpha_i)*X[2,])
        } else return(NULL)
    },.combine = rbind)

    if (is.null(r) || length(r)==0)
        r = NA
    else {
        colnames(r) <- colnames(intervals)
        rownames(r) <- NULL
    }
    if (!is.null(ridge.points)) {
        attr(r,"mesh") <- simplexes
    }
    # TODO attr(r,"values") <- y

    # points(r,pch=20)
    return(r)
}


#### 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 function or "unif" or "seq" (default) 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 mesh_roots(...) call
#' @param vectorized is f already vectorized ? (default: no)
#' @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="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="seq", mesh.sizes=11, maxerror_f=1E-9,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=mesh, mesh.sizes=mesh.sizes, maxerror_f=maxerror_f,tol=tol,...))

    if (sign!="upper" && sign!=1 && sign!="sup" && sign!=">" && !isTRUE(sign))
        stop("unknown sign: '",sign,"'")

    f_0 <- function(...) return(f(...)-threshold)
    r <- mesh_roots(f=f_0, vectorized=vectorized, intervals=intervals,mesh=mesh,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)
    new_mesh <- geometry::delaunayn(all_points,output.options=TRUE)

    if (isTRUE(vectorized))
        f_vec=function(x,...) f(x,...)
    else if (is.function(vectorized))
        f_vec=function(x,...) vectorized(x,...)
    else
        f_vec=Vectorize.function(f,dim=ncol(new_mesh$p))
    new_mesh$y = f_vec(new_mesh$p,...)
    # I = foreach (i.I = 1:nrow(new_mesh$tri)) %do% {
    #     X = new_mesh$p[new_mesh$tri[i.I,],,drop=FALSE] # NEVER USE "drop=F" => if "F" object exists, it will not crash but do a mess !!!
    #     y = f_vec(X,...) #apply(X,1,f, ...)
    #     if (all(y>=(threshold-2*maxerror_f))) {
    #         return(i.I)
    #     } else {
    #         return(NULL)
    #     }
    # }
    # I = unlist(Filter(Negate(is.null), I))
    I = which(apply(new_mesh$tri,1,function(i) ex_filter.tri(new_mesh$y[i]>=(threshold-10*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 color of the mesh
#' @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='black',...){
    col.rgb=col2rgb(color)/255
    plot(x=mesh$p,y=rep(y,length(mesh$p)),ylab="",col=rgb(col.rgb[1,],col.rgb[2,],col.rgb[3,],0.4),...)
    apply(mesh$tri,1,function(tri) lines(mesh$p[tri,],y=rep(y,length(tri)),col=color))
}

#' @title Plot a two dimensional mesh
#' @param mesh 2-dimensional mesh to draw
#' @param color color of the mesh
#' @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="unif",mesh.size=11,
#'                         intervals = matrix(c(1/2,5/2,1/2,5/2),nrow=2)))
plot2d_mesh = function(mesh,color='black',...){
    col.rgb=col2rgb(color)/255
    plot(mesh$p,col=rgb(col.rgb[1,],col.rgb[2,],col.rgb[3,],0.4),...)
    apply(mesh$tri,1,function(tri) polygon(mesh$p[c(tri,tri[1]),],border = color, col = rgb(0,0,0,.2)))
    # 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))
}

#' @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 color of the mesh
#' @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="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="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='black',...){
    col.rgb=col2rgb(color)/255
    package = load3d(engine3d)
    if (is.null(package)) return()
    plot3d(mesh$p,col=rgb(col.rgb[1,],col.rgb[2,],col.rgb[3,],0.4), package=package,...)
    apply(mesh$tri,1,function(tri) {
        quads3d(mesh$p[tri,],col=color,alpha=0.05, package=package)
        quads3d(mesh$p[tri,][c(4,3,2,1),],col=color,alpha=0.05, 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))
}


#### 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]
}