R/ndset.R

Defines functions classifyNDSet genNDSet genSample addNDSet2D addNDSet criterionPoints

Documented in addNDSet addNDSet2D classifyNDSet criterionPoints genNDSet genSample

## Functions related to multi-objective optimization


#' Calculate the criterion points of a set of points and ranges to find the set
#' of non-dominated points (Pareto points) and classify them into extreme
#' supported, non-extreme supported, non-supported.
#'
#' @param pts A data frame with a column for each variable in the solution
#'   space (can also be a rangePoints).
#' @param obj A p x n matrix(one row for each criterion).
#' @param crit Either max or min.
#' @param labels If \code{NULL} or "n" don't add any labels (empty string). If
#'   'coord' labels are the solution space coordinates. Otherwise number all
#'   points from one based on the solution space points.
#'
#' @return A data frame with columns x1, ..., xn, z1, ..., zp, lbl (label), nD
#'   (non-dominated), ext (extreme), nonExt (non-extreme supported).
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#' @examples
#' A <- matrix( c(3, -2, 1, 2, 4, -2, -3, 2, 1), nc = 3, byrow = TRUE)
#' b <- c(10,12,3)
#' pts <- integerPoints(A, b)
#' obj <- matrix( c(1,-3,1,-1,1,-1), byrow = TRUE, ncol = 3 )
#' criterionPoints(pts, obj, crit = "max", labels = "numb")
criterionPoints<-function(pts, obj, crit, labels = "coord") {
   n <- ncol(obj)
   p <- nrow(obj)
   pts <- .checkPts(pts, stopUnique = FALSE)
   zVal <- pts %*% t(obj)
   zVal <- round(zVal,10)
   zNames <- paste0("z", 1:p)
   colnames(zVal) <- zNames
   iP <- bind_cols(as.data.frame(pts), as.data.frame(zVal))
   iP$lbl <- as.character(1:nrow(iP))
   zVal <- addNDSet(zVal, crit = crit, keepDom = TRUE, dubND = TRUE)
   iP <- full_join(iP, zVal, by = zNames)

   if (is.null(labels)) {
      iP$lbl <- ""
      return(iP)
   }
   if (labels == "n") iP$lbl <- ""
   if (labels == "coord") iP$lbl <- df2String(iP[,1:n])
   return(iP)


#    iP <- as.data.frame(iP)
#    tol <- 1e-4
#    iP$oldLbl <- 1:length(iP$x1)
#    if (crit=="max") iP <- iP[order(-iP$z2,-iP$z1),]
#    if (crit=="min") iP <- iP[order(iP$z2,iP$z1),]
#    iP$lbl <- 1:length(iP$x1)  # note use in alg!
#    # classify non dom
#    iP$nD <- FALSE
#    iP$nD[1] <- TRUE  # upper left point
#    p1 <- iP$z1[1]; p2 <- iP$z2[1]  # current non dom point (due to sorting will p2 always be larger than or equal to current)
#    for (r in 2:length(iP$x1)) { # current point under consideration
#       if (abs(p2 - iP$z2[r]) < tol &
#           abs(p1 - iP$z1[r]) < tol) {
#          iP$nD[r] <- TRUE
#          p1 <- iP$z1[r]
#          p2 <- iP$z2[r]
#          next
#       }
#       if (crit == "max" &
#           p2 - iP$z2[r] > tol &
#           iP$z1[r] > p1 + tol) {
#          iP$nD[r] <- TRUE
#          p1 <- iP$z1[r]
#          p2 <- iP$z2[r]
#          next
#       }
#       if (crit == "min" &
#           iP$z2[r] - p2 > tol &
#           iP$z1[r] < p1 - tol) {
#          iP$nD[r] <- TRUE
#          p1 <- iP$z1[r]
#          p2 <- iP$z2[r]
#          next
#       }
#    }
#    # classify extreme supported
#    idx <- which(iP$nD & !duplicated(cbind(iP$z1,iP$z2), MARGIN = 1) )  # remove duplicated points
#    iP$ext <- FALSE
#    iP$ext[idx[1]] <- TRUE
#    iP$ext[idx[length(idx)]] <- TRUE
#    if (length(idx)<3) {
#       iP$nonExt <- FALSE
#       #return(iP)   # a single extreme
#    } else {
#       nD <- iP[idx,]
#       ul<-1
#       lr<-length(idx)
#       while (ul<length(idx)) {
#          # for (k in 1:1000) {
#          slope <- (nD$z2[lr]-nD$z2[ul])/(nD$z1[lr]-nD$z1[ul])
#          nD$val <- nD$z2-slope*nD$z1
#          #cat("val:",nD$val[ul],"max:",max(nD$val),"min:",min(nD$val),"\n")
#          if (crit=="max") {
#             i <- which.max(nD$val)
#             if (nD$val[ul]<nD$val[i] - tol) {
#                iP$ext[nD$lbl[i]] <- TRUE
#                lr <- i
#             } else {
#                ul <- lr
#                lr<-length(idx)
#             }
#          }
#          if (crit=="min") {
#             i <- which.min(nD$val)
#             if (nD$val[ul]>nD$val[i] + tol) {
#                iP$ext[nD$lbl[i]] <- TRUE
#                lr <- i
#             } else {
#                ul <- lr
#                lr<-length(idx)
#             }
#          }
#       }
#       # classify nonextreme supported
#       idxExt <- which(iP$ext)
#       iP$nonExt <- FALSE
#       if (length(idxExt)>1) {
#          for (i in 2:length(idxExt)) {
#             slope <- (iP$z2[idxExt[i]]-iP$z2[idxExt[i-1]])/(iP$z1[idxExt[i]]-iP$z1[idxExt[i-1]])
#             nDCand <- iP[idxExt[i-1]:idxExt[i],]
#             nDCand <- nDCand[nDCand$nD & !duplicated(cbind(nDCand$z1,nDCand$z2), MARGIN = 1),]
#             nDCand <- nDCand[c(-1,-length(nDCand$nD)),]
#             if (length(nDCand$nD)==0) next   # no points inbetween
#             for (j in 1:length(nDCand$nD)) {
#                slopeCur = (nDCand$z2[j]-iP$z2[idxExt[i-1]])/(nDCand$z1[j]-iP$z1[idxExt[i-1]])
#                if (abs(slope - slopeCur) < tol) iP$nonExt[nDCand$lbl[j]==iP$lbl] <- TRUE
#             }
#          }
#       }
#       # classify duplicates
#       idx <- which(iP$nD)
#       if (!length(idx)<2) {
#          for (i in 2:length(idx)) {
#             if (iP$ext[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
#                iP$ext[i] = TRUE
#                next
#             }
#             if (iP$nonExt[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
#                iP$nonExt[i] = TRUE
#             }
#          }
#       }
#    }
#    # set correct labels
#    iP$lbl <- iP$oldLbl
#    iP$oldLbl <- NULL
}



#' Add discrete points to a non-dominated set and classify them into extreme supported, non-extreme
#' supported, non-supported.
#'
#' @param pts A data frame with points to add (a column for each objective).
#' @param nDSet A data frame with current non-dominated set (NULL is none yet). Column names of the
#'   p objectives must be z1, ..., zp.
#' @param crit A max or min vector. If length one assume all objectives are optimized in the same
#'   direction.
#' @param keepDom Keep dominated points in output.
#' @param dubND Duplicated non-dominated points are classified as non-dominated.
#' @param classify Non-dominated points are classified into supported extreme (se), supported
#'   non-extreme (sne) and unsupported (us)
#'
#' @return A data frame with a column for each objective (`z` columns) and nd (non-dominated).
#'   Moreover if `classify` then columns `se`, `sne`, `us` and `cls`.
#'
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#'
#' @examples
#' nDSet <- data.frame(z1=c(12,14,16,18), z2=c(18,16,12,4))
#' pts <- data.frame(z1 = c(18,18,14,15,15), z2=c(2,6,14,14,16))
#' addNDSet(pts, nDSet, crit = "max")
#' addNDSet(pts, nDSet, crit = "max", keepDom = TRUE)
#' addNDSet(pts, nDSet, crit = "min")
#' addNDSet(c(2,2), nDSet, crit = "max")
#' addNDSet(c(2,2), nDSet, crit = "min")
#'
#' addNDSet(c(2,2), crit = "min")
#'
#' \donttest{
#' nDSet <- data.frame(z1=c(12,14,16,18), z2=c(18,16,12,4), z3 = c(1,7,0,6))
#' pts <- data.frame(z1=c(12,14,16,18), z2=c(18,16,12,4), z3 = c(2,2,2,6))
#' crit = c("min", "min", "max")
#' di <- c(1,1,-1)
#' li <- c(-1,20)
#' ini3D(argsPlot3d = list(xlim = li, ylim = li, zlim = li))
#' plotCones3D(nDSet, direction = di, argsPolygon3d = list(color = "green", alpha = 1),
#'             drawPoint = FALSE)
#' plotHull3D(nDSet, addRays = TRUE, direction = di)
#' plotPoints3D(nDSet, argsPlot3d = list(col = "red"), addText = "coord")
#' plotPoints3D(pts, addText = "coord")
#' finalize3D()
#' addNDSet(pts, nDSet, crit, dubND = FALSE)
#' addNDSet(pts, nDSet, crit, dubND = TRUE)
#' addNDSet(pts, nDSet, crit, dubND = TRUE, keepDom = TRUE)
#' addNDSet(pts, nDSet, crit, dubND = TRUE, keepDom = TRUE, classify = FALSE)
#' }
addNDSet<-function(pts, nDSet = NULL, crit = "max", keepDom = FALSE, dubND = FALSE,
                   classify = TRUE)
{
   if (is.data.frame(nDSet))
      nDSet$nd <- NULL; nDSet$se <- NULL; nDSet$sne <- NULL; nDSet$us <- NULL
   if (!is.data.frame(pts)) pts <- as.data.frame(.checkPts(pts, stopUnique = FALSE))
   p <- ncol(pts)
   direction <- .mToDirection(crit, p)
   colnames(pts) <- paste0("z",1:p)
   if (!is.null(nDSet)) {
      if (ncol(pts)!= ncol(nDSet))
         stop("Number of columns (not including classification colunms) must be the same!")
   } else {
      if (nrow(pts) == 1)
         return(mutate(pts, nd = TRUE, se = TRUE, sne = FALSE, us = FALSE, cls = "se"))
      nDSet <- pts[1,]
      pts <- pts[-1,]
   }
   set <- bind_rows(nDSet, pts)

   # new implem
   pf <- dplyr::if_else(direction == 1, "quo(", "quo(desc(")
   sf <- dplyr::if_else(direction == 1, ")", "))")
   args <- paste0("list(", paste0(pf, colnames(nDSet), sf, collapse = ", "), ")")
   args <- eval(parse(text=args))
   set <- set %>% dplyr::arrange(!!!args)  # sort based on direction, i.e. i<j => j cannot dom i
   set$nd <- TRUE
   if (dubND) {
      for (i in 1:(nrow(set)-1)) {
         for (j in (i+1):nrow(set)) {
            if (set$nd[j]) { # must check if i dom j and not eq
               if (all(set[i, 1:p]*direction <= set[j, 1:p]*direction)) { # i may dom j
                  if (!(all(set[i, 1:p]*direction == set[j, 1:p]*direction))) set$nd[j] <- FALSE
               }
            }
         }
      }
   } else {
      for (i in 1:(nrow(set)-1)) {
         for (j in (i+1):nrow(set)) {
            if (set$nd[j]) { # must check if i dom j
               if (all(set[i, 1:p]*direction <= set[j, 1:p]*direction)) set$nd[j] <- FALSE
            }
         }
      }
   }

   # old implem
   # idx <- (nrow(nDSet)+1):nrow(set)
   # set$nd <- TRUE
   # if (dubND) {
   #    for (i in idx) {
   #       for (j in 1:(i-1)) {
   #          if (!set$nd[j]) next  # if j dom don't check against i
   #          less <- set[j,1:p]*direction < set[i,1:p]*direction
   #          eq <- set[j,1:p]*direction == set[i,1:p]*direction
   #          if (all(less | eq) & !all(eq)) {set$nd[i] <- FALSE; break} # j strict dom i
   #          if (all(!less | eq) & !all(eq)) set$nd[j] <- FALSE # i strict dom j
   #       }
   #    }
   # } else {
   #    for (i in idx) {
   #       for (j in 1:(i-1)) {
   #          if (!set$nd[j]) next  # if j dom don't check against i
   #          less <- set[j,1:p]*direction < set[i,1:p]*direction
   #          eq <- set[j,1:p]*direction == set[i,1:p]*direction
   #          if (all(less | eq)) {set$nd[i] <- FALSE;  break} # j dom i
   #          if (all(!less | eq) & !all(eq)) set$nd[j] <- FALSE # i strict dom j
   #       }
   #    }
   # }

   if (!keepDom) set <- set %>% filter(.data$nd)
   if (classify) {
      set1 <- classifyNDSet(set[set$nd,1:p], direction)
      set <- set %>% tibble::rownames_to_column(var = "id")
      ids <- set$id[set$nd]
      set1 <- set1 %>% tibble::add_column(id = set$id[set$nd], .before = 1)
      set <- left_join(x = set, y = set1, by = c("id", colnames(pts)))
      if (keepDom) set <- set %>%
         mutate_if(is.factor, as.character) %>%
         tidyr::replace_na(list(se = FALSE, sne = FALSE, us = TRUE, cls = "d"))
      set <- set %>% select(-id)
   }
   return(set)

   # iP = pts
   # colnames(iP)[1:2] <- paste0("z", 1:2)
   # #iP <- round(iP,10)
   # rownames(iP) <- NULL
   # iP <- rbind(iP, nDSet)
   # tol <- 1e-4
   # iP$oldRowIdx <- 1:length(iP$z1)
   # if (crit=="max") iP <- iP[order(-iP$z2,-iP$z1),]
   # if (crit=="min") iP <- iP[order(iP$z2,iP$z1),]
   #
   # # classify non dom
   # iP$nD <- FALSE
   # iP$nD[1] <- TRUE  # upper left point
   # p1 <- iP$z1[1]; p2 <- iP$z2[1]  # current non dom point (due to sorting will p2 always be larger than or equal to current)
   # for (r in 2:length(iP$z1)) { # current point under consideration
   #    if (abs(p2 - iP$z2[r]) < tol &
   #        abs(p1 - iP$z1[r]) < tol) {
   #       iP$nD[r] <- TRUE
   #       p1 <- iP$z1[r]
   #       p2 <- iP$z2[r]
   #       next
   #    }
   #    if (crit == "max" &
   #        p2 - iP$z2[r] > tol &
   #        iP$z1[r] > p1 + tol) {
   #       iP$nD[r] <- TRUE
   #       p1 <- iP$z1[r]
   #       p2 <- iP$z2[r]
   #       next
   #    }
   #    if (crit == "min" &
   #        iP$z2[r] - p2 > tol &
   #        iP$z1[r] < p1 - tol) {
   #       iP$nD[r] <- TRUE
   #       p1 <- iP$z1[r]
   #       p2 <- iP$z2[r]
   #       next
   #    }
   # }
   # # iP$nD <- TRUE
   # # for (i in 2:nrow(iP)) { # remove non-dom z2
   # #    if (iP$z2[i-1]==iP$z2[i]) iP$nD[i] <- FALSE
   # # }
   # # if (!keepDom) iP <- iP[iP$nD,]
   # # for (i in 2:nrow(iP)) { # remove non-dom z1
   # #    if (iP$z1[i-1]==iP$z1[i]) iP$nD[i] <- FALSE
   # # }
   # if (!keepDom) iP <- iP[iP$nD,]
   # iP$rowIdx <- 1:nrow(iP)
   # # classify extreme supported
   # idx <- which(iP$nD & !duplicated(cbind(iP$z1,iP$z2), MARGIN = 1) )  # remove duplicated points
   # iP$ext <- FALSE
   # iP$ext[idx[1]] <- TRUE
   # iP$ext[idx[length(idx)]] <- TRUE
   # if (length(idx)<3) {
   #    iP$nonExt <- FALSE
   #    #return(iP)   # a single extreme
   # } else {
   #    nD <- iP[idx,]
   #    ul<-1
   #    lr<-length(idx)
   #    while (ul<length(idx)) {
   #       # for (k in 1:1000) {
   #       slope <- (nD$z2[lr]-nD$z2[ul])/(nD$z1[lr]-nD$z1[ul])
   #       nD$val <- nD$z2-slope*nD$z1
   #       #cat("val:",nD$val[ul],"max:",max(nD$val),"min:",min(nD$val),"\n")
   #       if (crit=="max") {
   #          i <- which.max(nD$val)
   #          if (nD$val[ul]<nD$val[i] - tol) {
   #             iP$ext[nD$rowIdx[i]] <- TRUE
   #             lr <- i
   #          } else {
   #             ul <- lr
   #             lr<-length(idx)
   #          }
   #       }
   #       if (crit=="min") {
   #          i <- which.min(nD$val)
   #          if (nD$val[ul]>nD$val[i] + tol) {
   #             iP$ext[nD$rowIdx[i]] <- TRUE
   #             lr <- i
   #          } else {
   #             ul <- lr
   #             lr<-length(idx)
   #          }
   #       }
   #    }
   #    # classify nonextreme supported
   #    idxExt <- which(iP$ext)
   #    iP$nonExt <- FALSE
   #    if (length(idxExt)>1) {
   #       for (i in 2:length(idxExt)) {
   #          slope <- (iP$z2[idxExt[i]]-iP$z2[idxExt[i-1]])/(iP$z1[idxExt[i]]-iP$z1[idxExt[i-1]])
   #          nDCand <- iP[idxExt[i-1]:idxExt[i],]
   #          nDCand <- nDCand[nDCand$nD & !duplicated(cbind(nDCand$z1,nDCand$z2), MARGIN = 1),]
   #          nDCand <- nDCand[c(-1,-length(nDCand$nD)),]
   #          if (length(nDCand$nD)==0) next   # no points inbetween
   #          for (j in 1:length(nDCand$nD)) {
   #             slopeCur = (nDCand$z2[j]-iP$z2[idxExt[i-1]])/(nDCand$z1[j]-iP$z1[idxExt[i-1]])
   #             if (abs(slope - slopeCur) < tol) iP$nonExt[nDCand$rowIdx[j]==iP$rowIdx] <- TRUE
   #          }
   #       }
   #    }
   #    # classify duplicates
   #    idx <- which(iP$nD)
   #    if (!length(idx)<2) {
   #       for (i in 2:length(idx)) {
   #          if (iP$ext[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
   #             iP$ext[i] = TRUE
   #             next
   #          }
   #          if (iP$nonExt[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
   #             iP$nonExt[i] = TRUE
   #          }
   #       }
   #    }
   # }
   # iP$rowIdx <- NULL
   # iP$oldRowIdx <- NULL
   # return(iP)
}




#' Add 2D discrete points to a non-dominated set and classify them into extreme
#' supported, non-extreme supported, non-supported.
#'
#' @param pts A data frame. It is assumed that z1 and z2 are in the two first columns.
#' @param nDSet A data frame with current non-dominated set (NULL is none yet).
#' @param crit Either max or min.
#' @param keepDom Keep dominated points.
#'
#' @return A data frame with columns z1 and z2, nD (non-dominated),
#'         ext (extreme), nonExt (non-extreme supported).
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#'
#' @examples
#' nDSet <- data.frame(z1=c(12,14,16,18), z2=c(18,16,12,4))
#' pts <- data.frame(z1 = c(18,18,14,15,15), z2=c(2,6,14,14,16))
#' addNDSet2D(pts, nDSet, crit = "max")
#' addNDSet2D(pts, nDSet, crit = "max", keepDom = TRUE)
#' addNDSet2D(pts, nDSet, crit = "min")
addNDSet2D<-function(pts, nDSet = NULL, crit = "max", keepDom = FALSE) {
   nDSet$nD <- NULL; nDSet$ext <- NULL; nDSet$nonExt <- NULL
   if (!is.null(nDSet))
      if (ncol(pts)!=ncol(nDSet))
         stop("Number of columns minus classification colunms must be the same!")
   iP = pts
   colnames(iP)[1:2] <- paste0("z", 1:2)
   #iP <- round(iP,10)
   rownames(iP) <- NULL
   iP <- rbind(iP, nDSet)
   tol <- 1e-4
   iP$oldRowIdx <- 1:length(iP$z1)
   if (crit=="max") iP <- iP[order(-iP$z2,-iP$z1),]
   if (crit=="min") iP <- iP[order(iP$z2,iP$z1),]

   # classify non dom
   iP$nD <- FALSE
   iP$nD[1] <- TRUE  # upper left point
   p1 <- iP$z1[1]; p2 <- iP$z2[1]  # current non dom point (due to sorting will p2 always be larger than or equal to current)
   for (r in 2:length(iP$z1)) { # current point under consideration
      if (abs(p2 - iP$z2[r]) < tol &
          abs(p1 - iP$z1[r]) < tol) {
         iP$nD[r] <- TRUE
         p1 <- iP$z1[r]
         p2 <- iP$z2[r]
         next
      }
      if (crit == "max" &
          p2 - iP$z2[r] > tol &
          iP$z1[r] > p1 + tol) {
         iP$nD[r] <- TRUE
         p1 <- iP$z1[r]
         p2 <- iP$z2[r]
         next
      }
      if (crit == "min" &
          iP$z2[r] - p2 > tol &
          iP$z1[r] < p1 - tol) {
         iP$nD[r] <- TRUE
         p1 <- iP$z1[r]
         p2 <- iP$z2[r]
         next
      }
   }
   # iP$nD <- TRUE
   # for (i in 2:nrow(iP)) { # remove non-dom z2
   #    if (iP$z2[i-1]==iP$z2[i]) iP$nD[i] <- FALSE
   # }
   # if (!keepDom) iP <- iP[iP$nD,]
   # for (i in 2:nrow(iP)) { # remove non-dom z1
   #    if (iP$z1[i-1]==iP$z1[i]) iP$nD[i] <- FALSE
   # }
   if (!keepDom) iP <- iP[iP$nD,]
   iP$rowIdx <- 1:nrow(iP)
   # classify extreme supported
   idx <- which(iP$nD & !duplicated(cbind(iP$z1,iP$z2), MARGIN = 1) )  # remove duplicated points
   iP$ext <- FALSE
   iP$ext[idx[1]] <- TRUE
   iP$ext[idx[length(idx)]] <- TRUE
   if (length(idx)<3) {
      iP$nonExt <- FALSE
      #return(iP)   # a single extreme
   } else {
      nD <- iP[idx,]
      ul<-1
      lr<-length(idx)
      while (ul<length(idx)) {
         # for (k in 1:1000) {
         slope <- (nD$z2[lr]-nD$z2[ul])/(nD$z1[lr]-nD$z1[ul])
         nD$val <- nD$z2-slope*nD$z1
         #cat("val:",nD$val[ul],"max:",max(nD$val),"min:",min(nD$val),"\n")
         if (crit=="max") {
            i <- which.max(nD$val)
            if (nD$val[ul]<nD$val[i] - tol) {
               iP$ext[nD$rowIdx[i]] <- TRUE
               lr <- i
            } else {
               ul <- lr
               lr<-length(idx)
            }
         }
         if (crit=="min") {
            i <- which.min(nD$val)
            if (nD$val[ul]>nD$val[i] + tol) {
               iP$ext[nD$rowIdx[i]] <- TRUE
               lr <- i
            } else {
               ul <- lr
               lr<-length(idx)
            }
         }
      }
      # classify nonextreme supported
      idxExt <- which(iP$ext)
      iP$nonExt <- FALSE
      if (length(idxExt)>1) {
         for (i in 2:length(idxExt)) {
            slope <- (iP$z2[idxExt[i]]-iP$z2[idxExt[i-1]])/(iP$z1[idxExt[i]]-iP$z1[idxExt[i-1]])
            nDCand <- iP[idxExt[i-1]:idxExt[i],]
            nDCand <- nDCand[nDCand$nD & !duplicated(cbind(nDCand$z1,nDCand$z2), MARGIN = 1),]
            nDCand <- nDCand[c(-1,-length(nDCand$nD)),]
            if (length(nDCand$nD)==0) next   # no points inbetween
            for (j in 1:length(nDCand$nD)) {
               slopeCur = (nDCand$z2[j]-iP$z2[idxExt[i-1]])/(nDCand$z1[j]-iP$z1[idxExt[i-1]])
               if (abs(slope - slopeCur) < tol) iP$nonExt[nDCand$rowIdx[j]==iP$rowIdx] <- TRUE
            }
         }
      }
      # classify duplicates
      idx <- which(iP$nD)
      if (!length(idx)<2) {
         for (i in 2:length(idx)) {
            if (iP$ext[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
               iP$ext[i] = TRUE
               next
            }
            if (iP$nonExt[i-1] & abs(iP$z2[i-1]-iP$z2[i])<tol & abs(iP$z1[i-1]-iP$z1[i])<tol) {
               iP$nonExt[i] = TRUE
            }
         }
      }
   }
   iP$rowIdx <- NULL
   iP$oldRowIdx <- NULL
   return(iP)
}


#' Generate a sample of points in dimension $p$.
#'
#' @param p Dimension of the points.
#' @param n Number of samples generated.
#' @param range The range of the points in each dimension (a vector or matrix with `p` rows).
#' @param random Random sampling.
#' @param sphere Generate points on a sphere.
#' @param box Generate points in boxes.
#' @param ... Further arguments passed on to the method for generating points. This must be done as
#'   lists (see examples). Currently the following arguments are supported:
#'
#'   * `argsSphere`: A list of arguments for generating points on a sphere:
#'      - `radius`: The radius of the sphere.
#'      - `center`: The center of the sphere.
#'      - `plane`: The plane used.
#'      - `below`: Either true (generate points below the plane), false (generate points above the
#'                 plane) or `NULL` (generated on the whole sphere).
#'      - `factor`: If using a plane. Then the factor to multiply `n` with, so generate enough points
#'                  below/above the plane.
#'      - `closeToPlane`: If TRUE only return points close to the plane.
#'   * `argsBox`: A list of arguments for generating points inside boxes:
#'      - `intervals`: Number of intervals to split the length of the range into. That is, each
#'                     range is divided into `intervals` (sub)intervals and only the lowest/highest
#'                     subrange is used.
#'      - `cor`: How to correlate indices. If `'idxAlt'` then alternate the intervals (high/low)
#'               for each dimension. For instance if `p = 3` and the first dimension is in the high
#'               interval range then the second will be in the low interval range and third in the
#'               high interval range again. If `idxRand` then choose the low/high interval range
#'               for each dimension based on `prHigh`. If `idxSplit` then select
#'               `floor(p/2):ceiling(p/2)` dimensions for the high interval range and the other for
#'               the low interval range.
#'      - `prHigh`: Probability for choosing the high interval range in each dimension.
#'
#' @details Note having ranges with different length when using the sphere method, doesn't make
#'   sense. The best option is properly to use a center and radius here. Moreover, as for higher
#'   `p` you may have to use a larger radius than half of the desired interval range.
#'
#' @return A data frame with `p` columns
#' @export
#'
#' @examples
#' ### Using random
#' ## p = 2
#' range <- matrix(c(1,100, 50,100), ncol = 2, byrow = TRUE )
#' pts <- genSample(2, 1000, range = range, random = TRUE)
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plot(pts)
#'
#' \donttest{
#' ## p = 3
#' range <- matrix(c(1,100, 50,100, 10,50), ncol = 2, byrow = TRUE )
#' ini3D()
#' pts <- genSample(3, 1000, range = range, random = TRUE)
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plotPoints3D(pts)
#' finalize3D()
#' }
#'
#' ## other p
#' p <- 10
#' range <- c(1,100)
#' pts <- genSample(p, 1000, range = range, random = TRUE)
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#'
#' ### Using sphere
#' ## p = 2
#' range <- c(1,100)
#' cent <- rep(range[1] + (range[2]-range[1])/2, 2)
#' pts <- genSample(2, 1000, range = range)
#' dim(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plot(pts, asp=1)
#' abline(sum(cent^2)/cent[1], -cent[2]/cent[1])
#'
#' cent <- c(100,100)
#' r <- 75
#' planeC <- c(cent+r/3)
#' planeC <- c(planeC, -sum(planeC^2))
#' pts <- genSample(2, 100,
#'   argsSphere = list(center = cent, radius = r, below = FALSE, plane = planeC, factor = 6))
#' dim(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plot(pts, asp=1)
#' abline(-planeC[3]/planeC[1], -planeC[2]/planeC[1])
#'
#' pts <- genSample(2, 100, argsSphere = list(center = cent, radius = r, below = NULL))
#' dim(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plot(pts, asp=1)
#'
#' \donttest{
#' ## p = 3
#' ini3D()
#' range <- c(1,100)
#' cent <- rep(range[1] + (range[2]-range[1])/2, 3)
#' pts <- genSample(3, 1000, range = range)
#' dim(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' rgl::spheres3d(cent, radius=49.5, color = "grey100", alpha=0.1)
#' plotPoints3D(pts)
#' rgl::planes3d(cent[1],cent[2],cent[3],-sum(cent^2), alpha = 0.5, col = "red")
#' finalize3D()
#'
#' ini3D()
#' cent <- c(100,100,100)
#' r <- 75
#' planeC <- c(cent+r/3)
#' planeC <- c(planeC, -sum(planeC^2))
#' pts <- genSample(3, 100,
#'   argsSphere = list(center = cent, radius = r, below = FALSE, plane = planeC, factor = 6))
#' rgl::spheres3d(cent, radius=r, color = "grey100", alpha=0.1)
#' plotPoints3D(pts)
#' rgl::planes3d(planeC[1],planeC[2],planeC[3],planeC[4], alpha = 0.5, col = "red")
#' finalize3D()
#'
#' ini3D()
#' pts <- genSample(3, 10000, argsSphere = list(center = cent, radius = r, below = NULL))
#' Rfast::colMinsMaxs(as.matrix(pts))
#' rgl::spheres3d(cent, radius=r, color = "grey100", alpha=0.1)
#' plotPoints3D(pts)
#' finalize3D()
#' }
#'
#' ## Other p
#' p <- 10
#' cent <- rep(0,p)
#' r <- 100
#' pts <- genSample(p, 100000, argsSphere = list(center = cent, radius = r, below = NULL))
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' apply(pts,1, function(x){sqrt(sum((x-cent)^2))}) # test should be approx. equal to radius
#'
#'
#' ### Using box
#' ## p = 2
#' range <- matrix(c(1,100, 50,100), ncol = 2, byrow = TRUE )
#' pts <- genSample(2, 1000, range = range, box = TRUE, argsBox = list(cor = "idxAlt"))
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plot(pts)
#'
#' pts <- genSample(2, 1000, range = range, box = TRUE, argsBox = list(cor = "idxAlt",
#'                  intervals = 6))
#' plot(pts)
#'
#' pts <- genSample(2, 1000, range = range, box = TRUE, argsBox = list(cor = "idxRand"))
#' plot(pts)
#' pts <- genSample(2, 1000, range = range, box = TRUE,
#'                  argsBox = list(cor = "idxRand", prHigh = c(0.1,0.6)))
#' points(pts, pch = 3, col = "red")
#' pts <- genSample(2, 1000, range = range, box = TRUE,
#'                  argsBox = list(cor = "idxRand", prHigh = c(0,0)))
#' points(pts, pch = 4, col = "blue")
#'
#' pts <- genSample(2, 1000, range = range, box = TRUE, argsBox = list(cor = "idxSplit"))
#' plot(pts)
#'
#' \donttest{
#' ## p = 3
#' range <- matrix(c(1,100, 1,200, 1,50), ncol = 2, byrow = TRUE )
#' ini3D(argsPlot3d = list(box = TRUE, axes = TRUE))
#' pts <- genSample(3, 1000, range = range, box = TRUE, , argsBox = list(cor = "idxAlt"))
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' plotPoints3D(pts)
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(box = TRUE, axes = TRUE))
#' pts <- genSample(3, 1000, range = range, box = TRUE, ,
#'                  argsBox = list(cor = "idxAlt", intervals = 6))
#' plotPoints3D(pts)
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(box = TRUE, axes = TRUE))
#' pts <- genSample(3, 1000, range = range, box = TRUE, , argsBox = list(cor = "idxRand"))
#' plotPoints3D(pts)
#' pts <- genSample(3, 1000, range = range, box = TRUE, ,
#'                  argsBox = list(cor = "idxRand", prHigh = c(0.1,0.6,0.1)))
#' plotPoints3D(pts, argsPlot3d = list(col="red"))
#' finalize3D()
#'
#' ini3D(argsPlot3d = list(box = TRUE, axes = TRUE))
#' pts <- genSample(3, 1000, range = range, box = TRUE, , argsBox = list(cor = "idxSplit"))
#' plotPoints3D(pts)
#' finalize3D()
#' }
#'
#' ## other p
#' p <- 10
#' range <- c(1,100)
#' pts <- genSample(p, 1000, range = range, box = TRUE, argsBox = list(cor = "idxSplit"))
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
genSample <- function(p, n, range = c(1,100), random = FALSE, sphere = TRUE, box = FALSE, ...) {
   if (!is.matrix(range)) range <- matrix(range, ncol = 2)
   if (nrow(range) == 1) range <-  matrix(rep(range, each=p), nrow=p)
   args <- list(...)

   if (random) {
      sphere = FALSE
      argsRandom <- mergeLists(list(), args$argsRandom)
      set <- apply(range, 1, function(x) sample(x[1]:x[2], n, replace = TRUE) )
   }

   if (sphere) {
      argsSphere <-
         mergeLists(list(
            radius = min(Rfast::rowMins((range[,2] - range[,1])/2, value = TRUE)),
            center = range[,1] + (range[,2]-range[,1])/2,
            plane = NULL,
            below = TRUE,
            factor = 2 + 1/log(n),
            closeToPlane = FALSE
         ), args$argsSphere)
      if (is.null(argsSphere$plane))
         argsSphere$plane <- c(argsSphere$center, -sum(argsSphere$center ^ 2))
      if (!is.null(argsSphere$below)) n1 <- floor(n * argsSphere$factor) else n1 <- n # increase so have n points below/above
      if (!is.null(argsSphere$closeToPlane))
         if (argsSphere$closeToPlane) n1 <- floor(n * n * (2 + 1/n * 1000))

      # generate points https://math.stackexchange.com/questions/87230/
      #                 picking-random-points-in-the-volume-of-sphere-with-uniform-probability
      pts <- matrix(stats::rnorm(p*n1), nrow = p)  # sample
      pts <- t(apply(pts, 2, function(x){x * argsSphere$radius/sqrt(sum(x^2))}))
      # old method
      # sp <- sphereplot::pointsphere(n,c(0,360),c(-90,90),c(argsSphere$radius,argsSphere$radius))
      # pts <- sphereplot::sph2car(sp, deg = TRUE)

      pts <- pts + matrix(rep(argsSphere$center, dim(pts)[1]), byrow = TRUE, ncol=p)  # shift center
      if (!is.null(argsSphere$below)) {
         set <- matrix(rep(0,p), nrow = 1)
         if(argsSphere$below) {
            for (i in 1:nrow(pts)) { # remove all points below/above the plane
               if (sum(argsSphere$plane[1:p] *  pts[i,]) <= -argsSphere$plane[p+1]) {
                  set <- rbind(set, pts[i,])
               }
               if (nrow(set) == n + 1) break
            }
         }
         if (!argsSphere$below) {
            for (i in 1:nrow(pts)) { # remove all points below/above the plane
               if (sum(argsSphere$plane[1:p] *  pts[i,]) >= -argsSphere$plane[p+1]) {
                  set <- rbind(set, pts[i,])
               }
               if (nrow(set) == n + 1) break
            }
         }
         set <- set[-1,]  # remove first dummy row
      } else if (argsSphere$closeToPlane) {
         set <- matrix(rep(0,p), nrow = 1)
         delta <- argsSphere$radius/25
         for (i in 1:nrow(pts)) { # remove all points not close to the plane
            if (sum(argsSphere$plane[1:p] *  pts[i,]) >= -argsSphere$plane[p+1] - delta &
                sum(argsSphere$plane[1:p] *  pts[i,]) <= -argsSphere$plane[p+1] + delta) {
               set <- rbind(set, pts[i,])
            }
            if (nrow(set) == n + 1) break
         }
         set <- set[-1,]
      } else set <- pts
      set <- round(set)
   }

   if (box) {
      argsBox <-
         mergeLists(list(
            intervals = 3,
            cor = "idxAlt",
            prHigh = rep(0.5, p)
         ), args$argsBox)
      invLength <- ceiling((range[,2]-range[,1])/argsBox$intervals)
      rngL <- cbind(range[,1], range[,1] + invLength)
      rngH <- cbind(range[,2] - invLength, range[,2])
      ptsH <- apply(rngH, 1, function(x) sample(x[1]:x[2], n, replace = TRUE) )
      ptsL <- apply(rngL, 1, function(x) sample(x[1]:x[2], n, replace = TRUE) )
      ptsL <- split(ptsL, seq(nrow(ptsL)))
      ptsH <- split(ptsH, seq(nrow(ptsH)))
      if (argsBox$cor == "idxAlt") {
         high <- stats::rbinom(n, 1, 0.5)
         high <- lapply(high, function(x) rep(c(x,!x),times = p)[1:p])
      }
      if (argsBox$cor == "idxRand") {
         high <-
            sapply(argsBox$prHigh, function(pr)
               sample(0:1, n, replace = TRUE, prob = c(1 - pr, pr)))
         high <- split(high, seq(nrow(high)))
      }
      if (argsBox$cor == "idxSplit") {
         if (p %% 2 == 0) r <- p/2  # even number
         if (p %% 2 == 1) r <- floor(p/2):ceiling(p/2)
         s <- expand.grid(as.data.frame(matrix(rep(0:1,p), ncol = p)))
         s <- s[apply(s, 1, function(x) sum(x) %in% r),]
         rownames(s) <- NULL
         high <- sample(1:nrow(s), n, replace = TRUE)
         high <- lapply(high, function(i) s[i,])
      }
      lst <- lapply(1:length(ptsL),function(i) list(pts = rbind(high[[i]], ptsL[[i]], ptsH[[i]]) ))
      set <- t(sapply(lst, function(x) apply(x$pts,2, function(x) x[x[1]+2])))
   }
   if (nrow(set) != n) warning("Only ", nrow(set), " samples generated!")
   return(set)
}




#' Generate a sample of nondominated points.
#'
#' @param p Dimension of the points.
#' @param n Number of samples generated (note only a subset of these will be non-dominated).
#' @param range The range of the points in each dimension (a vector or matrix with `p` rows).
#' @param random Random sampling.
#' @param sphere Generate points on a sphere.
#' @param box Generate points in boxes.
#' @param keep Keep dominated points also.
#' @param crit Criteria used (a vector of min/max).
#' @param dubND Should duplicated non-dominated points be considered as non-dominated.
#' @param ... Further arguments passed on to [`genSample`].
#'
#' @return A data frame with `p+1` columns (last one indicate if dominated or not).
#' @export
#'
#' @examples
#' \donttest{
#' range <- matrix(c(1,100, 50,100, 10,50), ncol = 2, byrow = TRUE )
#' pts <- genNDSet(3, 50, range = range, random = TRUE, keep = TRUE)
#' head(pts)
#' Rfast::colMinsMaxs(as.matrix(pts))
#' 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)))
#' plotPoints3D(pts[,1:3])
#' plotPoints3D(pts[pts$nd,1:3], argsPlot3d = list(col = "red", size = 10))
#' plotCones3D(pts[pts$nd,1:3], argsPolygon3d = list(alpha = 1))
#' finalize3D()
#'
#' ini3D()
#' range <- c(1,100)
#' cent <- rep(range[1] + (range[2]-range[1])/2, 3)
#' pts <- genNDSet(3, 100, range = range, sphere = TRUE, keep = TRUE,
#'        argsSphere = list(center = cent))
#' rgl::spheres3d(cent, radius=49.5, color = "grey100", alpha=0.1)
#' plotPoints3D(pts)
#' plotPoints3D(pts[pts$nd,], argsPlot3d = list(col = "red", size = 10))
#' rgl::planes3d(cent[1],cent[2],cent[3],-sum(cent^2), alpha = 0.5, col = "red")
#' finalize3D()
#'
#' ini3D()
#' cent <- c(100,100,100)
#' r <- 75
#' planeC <- c(cent+r/3)
#' planeC <- c(planeC, -sum(planeC^2))
#' pts <- genNDSet(3, 100, keep = TRUE,
#'   argsSphere = list(center = cent, radius = r, below = FALSE, plane = planeC, factor = 6))
#' rgl::spheres3d(cent, radius=r, color = "grey100", alpha=0.1)
#' plotPoints3D(pts)
#' plotPoints3D(pts[pts$nd,], argsPlot3d = list(col = "red", size = 10))
#' rgl::planes3d(planeC[1],planeC[2],planeC[3],planeC[4], alpha = 0.5, col = "red")
#' finalize3D()
#' }
genNDSet <-
   function(p,
            n,
            range = c(1, 100),
            random = FALSE,
            sphere = TRUE,
            box = FALSE,
            keep = FALSE,
            crit = "min",
            dubND = TRUE,
            ...) {

   # if (p!=3) stop("Currently only works for p = 3!")
   set <- genSample(p, n, range = range, random = random, sphere = sphere, box = box, ...)
   set <- addNDSet(set, crit = crit, keepDom = keep, dubND = dubND, classify = FALSE)
   # set <- set[order(set[,1],set[,2],set[,3]),]
   # set <- as.data.frame(set)
   # colnames(set) <- c("x","y","z")
   # set$dom <- FALSE
   # for (i in 1:(dim(set)[1]-1)) {
   #    for (j in (i+1):dim(set)[1]) {
   #       if (!set$dom[j]) {
   #          if (set$x[i]<=set$x[j] && set$y[i]<=set$y[j] && set$z[i]<=set$z[j]) set$dom[j] <- TRUE
   #       }
   #    }
   # }
   # if (!keep) set <- set[!set$dom,1:3]
   # row.names(set) <- NULL
   return(set)
}



#' Classify a set of nondominated points
#'
#' The classification is supported (true/false), extreme (true/false), supported non-extreme
#' (true/false)
#'
#' @param pts A set of non-dominated points. It is assumed that `ncol(pts)` equals the number of
#'   objectives ($p$).
#' @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 (minimize objective $i$). If negative, consider the
#'   i'th column of the `pts` minus a value greater than on equal zero (maximize objective $i$).
#'
#' @note It is assumed that `pts` are nondominated.
#'
#' @return The ND set with classification columns.
#' @export
#'
#' @examples
#' \donttest{
#' pts <- matrix(c(0,0,1, 0,1,0, 1,0,0, 0.5,0.2,0.5, 0.25,0.5,0.25), ncol = 3, byrow = TRUE)
#' ini3D(argsPlot3d = list(xlim = c(min(pts[,1])-2,max(pts[,1])+2),
#'   ylim = c(min(pts[,2])-2,max(pts[,2])+2),
#'   zlim = c(min(pts[,3])-2,max(pts[,3])+2)))
#' plotHull3D(pts, addRays = TRUE, argsPolygon3d = list(alpha = 0.5), useRGLBBox = TRUE)
#' pts <- classifyNDSet(pts[,1:3])
#' plotPoints3D(pts[pts$se,1:3], argsPlot3d = list(col = "red"))
#' plotPoints3D(pts[!pts$sne,1:3], argsPlot3d = list(col = "black"))
#' plotPoints3D(pts[!pts$us,1:3], argsPlot3d = list(col = "blue"))
#' plotCones3D(pts[,1:3], rectangle = TRUE, argsPolygon3d = list(alpha = 1))
#' finalize3D()
#' pts
#'
#' pts <- matrix(c(0,0,1, 0,1,0, 1,0,0, 0.2,0.1,0.1, 0.1,0.45,0.45), ncol = 3, byrow = TRUE)
#' di <- -1 # maximize
#' ini3D(argsPlot3d = list(xlim = c(min(pts[,1])-1,max(pts[,1])+1),
#'   ylim = c(min(pts[,2])-1,max(pts[,2])+1),
#'   zlim = c(min(pts[,3])-1,max(pts[,3])+1)))
#' plotHull3D(pts, addRays = TRUE, argsPolygon3d = list(alpha = 0.5), direction = di,
#'            addText = "coord")
#' pts <- classifyNDSet(pts[,1:3], direction = di)
#' plotPoints3D(pts[pts$se,1:3], argsPlot3d = list(col = "red"))
#' plotPoints3D(pts[!pts$sne,1:3], argsPlot3d = list(col = "black"))
#' plotPoints3D(pts[!pts$us,1:3], argsPlot3d = list(col = "blue"))
#' plotCones3D(pts[,1:3], rectangle = TRUE, argsPolygon3d = list(alpha = 1), direction = di)
#' finalize3D()
#' pts
#'
#' pts <- matrix(c(0,0,1, 0,0,1, 0,1,0, 0.5,0.2,0.5, 1,0,0, 0.5,0.2,0.5, 0.25,0.5,0.25), ncol = 3,
#'               byrow = TRUE)
#' classifyNDSet(pts)
#'
#' pts <- genNDSet(3,50)[,1:3]
#' ini3D(argsPlot3d = list(xlim = c(0,max(pts$z1)+2),
#'   ylim = c(0,max(pts$z2)+2),
#'   zlim = c(0,max(pts$z3)+2)))
#' plotHull3D(pts, addRays = TRUE, argsPolygon3d = list(alpha = 0.5))
#' pts <- classifyNDSet(pts[,1:3])
#' plotPoints3D(pts[pts$se,1:3], argsPlot3d = list(col = "red"))
#' plotPoints3D(pts[!pts$sne,1:3], argsPlot3d = list(col = "black"))
#' plotPoints3D(pts[!pts$us,1:3], argsPlot3d = list(col = "blue"))
#' finalize3D()
#' pts
#' }
classifyNDSet <- function(pts, direction = 1) {
   pts <- .checkPts(pts, stopUnique = FALSE)
   p <- ncol(pts)
   colnames(pts) <- paste0("z", 1:p)
   idx <- duplicated(pts)
   pts <- pts %>% dplyr::as_tibble() %>% dplyr::mutate(id = 1:nrow(pts)) #%>%  tibble::rownames_to_column(var = "rn")
   if (nrow(pts) == 1) {
      pts <- as.data.frame(pts)
      return(cbind(select(pts,-id), se = TRUE, sne = FALSE, us = FALSE, cls = "se"))
   }
   if (length(direction) != p) direction = rep(direction[1],p)

   # find hull of the unique points and classify
   set <- convexHull(pts[!idx,1:p], addRays = TRUE, direction = direction)
   hull <- set$hull
   set <- set$pts
   d <- dimFace(set[,1:p])
   if (d != p) stop("The points including rays don't seem to define a hull of dimension ", p, "!")
   set <- dplyr::mutate(set, se = dplyr::if_else(.data$vtx,TRUE,FALSE))
   set <- dplyr::mutate(set, sne = FALSE, us = FALSE, id = 1:nrow(set))
   chk <- set %>% dplyr::filter(!.data$vtx)
   if (nrow(chk) != 0) {
      val <- inHull(chk[,1:p], set[set$vtx,1:p])
      set$us[chk$id[which(val == 1)]] <- TRUE
      set$sne[chk$id[which(val == 0)]] <- TRUE
   }
   set <- set %>% # tidy and add old id
      dplyr::filter(.data$pt == 1) %>%
      dplyr::mutate(cls = dplyr::if_else(.data$se, "se", dplyr::if_else(.data$sne, "sne", "us"))) %>%
      dplyr::select(tidyselect::all_of(1:p), c("se", "sne", "us", "cls")) %>%
      dplyr::mutate(id = which(!idx))
   set1 <- set %>% left_join(x = set, y = pts[idx,], by = paste0("z", 1:p)) # match id of duplicates
   set1 <- set1 %>%
      dplyr::filter(!is.na(.data$id.y)) %>%
      dplyr::mutate(id.x = .data$id.y) %>% dplyr::select(.data$z1:.data$id.x)
   set <- bind_rows(set, pts[idx,]) %>% arrange(id) %>% select(-id)
   if (nrow(set1) > 0) for (i in 1:nrow(set1)) set[set1$id.x[i],] <- set1[i, 1:(p+4)]
   return(set)
}

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.