R/vennplot.R

In vennplot: Venn Diagrams in 2D and 3D

#' Draw Venn and Euler diagram in 2D or 3D
#'
#' @param disjoint.combinations Named numeric vector or data.frame where each column should be factor.  See Details.
#' @param vars Extract specific variables of data.frame as \code{disjoint.combinations}. If \code{vars = NULL}, all the information of data.frame will be extracted.
#' @param Delta The length of step for method "lineSearch" or the initial interval of test points for method "NelderMead".
#' @param ThreeD Draw Venn diagram in 3D. See Examples.
#' @param lambda It can be \code{NULL} or a numeric vector. If \code{lambda = NULL}, the loss function optimize lambda, else, based on the given lambda, loss function will calculate stress respectively then return the minimum one and corresponding lambdas.
#' @param stressWay If data set can be separated into a few groups, there will be two ways to express stress: one is to sum up all the stress (named "sum"; default), the other is to use total TSS divide by total RSS (named "combine").
#' @param delta Closeness between groups.
#' @param weight The weight of \code{disjoint.combinations}. It should have the same length with \code{disjoint.combinations}.
#' @param expand If some balls should not intersect and the code fails to detect it. It is possible to be fixed manually but sacrificing stress.
#' @param twoWayGenerate Boolean factor, if false, any missing intersections are set as zero.
#' @param scaleSearch Provide multiple methods to optimize scale lambda. The default method is "NelderMead". See Details.
#' @param twoWaySearch If two way intersections are missing, multiple methods are available to generate two way intersections. The default method is "lineSearch". See Details.
#' @param scaleSeachTolerance A list with tolerance value and boolean factor " proportional". The loop of NelderMead and lineSearch in scaleSearch will end when the difference or  proportional difference matches the tolerance value.
#' @param distanceTolerance A list with tolerance value and boolean factor " proportional". The Newton method of finding distance will end when the difference or  proportional difference matches the tolerance value.
#' @param lossTolerance A list with ToleranceofLoss, maximumStep, ALPHA, ToleranceofStepsize and boolean factor "proportional". If ALPHA is null, the step size will be searched through Newton method and it will stop when step reaches the maximum step or the difference matches ToleranceofStepsize; else step size will be fixed with ALPHA . The loss will end when the difference or proportional difference or the total loss value matches the "ToleranceofLoss".
#' @param stressBound The loop of method NelderMead will stop when stress is beyond the stressBound.
#' @param maximumStep The maximum searching step for method NelderMead and Newton method of calculating distance.
#' @param planeSize The plane size of calculating disjoint intersections numerically.
#' @param lower The lower bound of the interval to be searched for the "goldenSectionSearch" and "L-BFGS-B". See Details.
#' @param upper The upper bound of the interval to be searched for the "goldenSectionSearch" and "Brent". See Details.
#' @param control A list of control parameters. See Details
#' @param hessian Logical. A numerically differentiated Hessian matrix be returned or not. See Details.
#' @param mar Plot margins.
#' @param cols Color of balls. If \code{NULL}, rainbow color will be set.
#' @param alpha Color darkness.
#' @param smooth For 3D plot, if true, the balls will be much more smoother. However, based on the high resolution, if the number of balls is too much, when rotating, the new window stumbles.
#' @param ... Any further graphical parameters to be passed to the \code{plot} function.
#'
#' @details
#' 1. One way sets must be given in \code{disjoint.combination}. e.g.\code{disjoint.combination = c( B=2, AB=0.5)} is not allowed.  \code{disjoint.combination = c(A = 0, B=2, AB=0.5)} works.
#' 2. Except "NelderMead" and "lineSearch", "goldenSectionSearch" in \code{scaleSearch} and \code{twoWaySearch} is based on \code{\link{optimize}} and the rest methods are based on \code{\link{optim}}.
#' 3. \code{lower}, \code{upper}, \code{control} and \code{hessian} share the same parameters with \code{\link{optim}}, and \code{lower}, \code{upper} can also be used in \code{\link{optimize}}
#'
#' @author Zehao Xu and Wayne Oldford
#' @return An object of the class \code{vennplot} with following components:
#' \describe{
#'   \item{xy}{centres of the balls (columns are (\code{x}, \code{y}) or (\code{x}, \code{y}, \code{z}) coordinates).}
#'   \item{radius}{radii of the balls.}
#'   \item{loss}{total loss of \code{vennplot}.}
#'   \item{stress}{stress value for solution.}
#' }
#' @examples
#' # 3D Venn plot with arbitray sets
#' disjoint.combinations = c(A=80, B=50,C=100, D = 100,E = 100,
#'                           "A&C"=30, "A&D"= 30,"B&E" = 30, "A&E" = 40, h = 40, "B&h" = 10)
#' ve = vennplot(disjoint.combinations, ThreeD = TRUE)
#'
#' # data frame
#' vennplot(disjoint.combinations = sharks, vars = c("Au","USA","Fa","Sex"),
#'          scaleSearch = "lineSearch", expand = 1.1)
#'
#' @export
vennplot <- function(disjoint.combinations = NULL, vars = NULL, Delta = 0.1,
                     ThreeD = FALSE, lambda = NULL, stressWay = c("sum","combine"),
                     delta = 0.01, weight = NULL, expand = NULL, twoWayGenerate = FALSE,
                     scaleSearch = c("NelderMead", "lineSearch", "goldenSectionSearch",
                                     "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
                     twoWaySearch = c("lineSearch", "NelderMead", "goldenSectionSearch",
                                      "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
                     scaleSeachTolerance = list(value = 1e-5,  proportional = FALSE),
                     distanceTolerance = list(value = 1e-5,  proportional = FALSE),
                     lossTolerance = list(ToleranceofLoss = 1e-10, maximumStep = 10, ALPHA = 1e-2,
                                          ToleranceofStepsize = 1e-5,  proportional = FALSE),
                     stressBound = 1e-3, maximumStep = 50, planeSize = 50,
                     lower=  -Inf, upper = Inf, control = list(), hessian = FALSE,
                     mar = rep(1,4), cols = NULL, alpha = 0.3, smooth = FALSE, ...){

  if(is.null(disjoint.combinations)) stop("combinations should not be empty")

  if (is.data.frame(disjoint.combinations)){
    disjoint.combinations <- extractCombinations(disjoint.combinations, vars = vars)
  }
  combinations <- disjoint2combinations(disjoint.combinations)
  # reorder disjoint.combinations with ways
  reOrder <- reorderDisjointCombinations(combinations, disjoint.combinations)
  disjoint.combinations <- reOrder$newDisjointCombinations
  # scale proportionately
  combProp <- combinations/sum(disjoint.combinations)
  disProp <- disjoint.combinations/sum(disjoint.combinations)
  # Check weight vector
  if(is.null(weight)){
    weight <- rep(1,length(disjoint.combinations))
  } else {
    if(length(weight)!=length(disjoint.combinations)){
      stop("weight must be NULL or the same length as the number of disjoint combinations")
    }
  }
  weight <- weight[reOrder$newOrder]
  disjointSetNames <- names(disProp)
  names(weight) <- disjointSetNames

  # Get the number of intersections
  numWays <- str_count(disjointSetNames, pattern = "&") + 1

  if(max(numWays)==2){
    nonEmptyTwoWays <- which(combProp!=0)
    weight <- weight[nonEmptyTwoWays]
    disjointSetNames <- disjointSetNames[nonEmptyTwoWays]
    numWays <- numWays[nonEmptyTwoWays]
    combProp <- combProp[nonEmptyTwoWays]
    disProp <- disProp[nonEmptyTwoWays]
  }

  oneWays <- which(numWays == 1)
  # Check colours
  if (is.null(cols)) {
    cols <- rainbow(length(oneWays), alpha = alpha)
  }
  oneWaySetName <- disjointSetNames[oneWays]
  oneWaySet <- combProp[oneWays]
  m <- length(oneWays)

  # and for those that are larger than one way set
  if (length(disjointSetNames) == length(oneWays)) {
    largerThanOneWaySetName <- NULL
  } else {
    largerThanOneWaySetName <- disjointSetNames[-oneWays]
    # All names appearing in two and higher order ways must appear
    # as individual input sets too
    if (!all(unique(unlist(str_split(largerThanOneWaySetName,"&"))) %in% oneWaySetName)){
      stop("Some intersection sets contain sets which do not appear individually.")
    }
  }
  # Detect disjoint groups of sets and return groups in a list
  groups <- groupDetection(largerThanOneWaySetName = largerThanOneWaySetName,
                           oneWaySetName= oneWaySetName)
  # Order groups from largest to smallest
  groupOrder <- order(sapply(groups, length), decreasing = T)
  groups <- groups[groupOrder]
  ngroups <- length(groups)

  # get the combinations and proportions separated by groups

  combPropGroup <- list()
  disPropGroup <- list()
  for(i in 1:ngroups){
    groupMember <- groups[[i]]
    groupWay <- oneWaySetName[groupMember]
    for(j in 1:length(disjointSetNames)){
      if(any(groupWay %in% str_split(disjointSetNames[j],"&")[[1]])){
        groupWay <- c(groupWay, disjointSetNames[j])}
    }
    groupWay <- unique(groupWay)
    combPropGroup[[i]] <- combProp[which(disjointSetNames %in% groupWay)]
    disPropGroup[[i]] <- disProp[which(disjointSetNames %in% groupWay)]
  }

  # Calculate radius
  radius <- if (ThreeD){
    (3*oneWaySet/(4*pi))^(1/3)
  } else {
    sqrt(oneWaySet/pi)
  }
  # Get radii within each group
  radiusGroup <- lapply(groups,
                        function(grp){radius[grp]})

  weightGroup <- lapply(disPropGroup,
                        function(prop){
                          weight[which(names(weight) %in% names(prop)==TRUE)]})
  # lossFunction
  groupCentre <- list()
  if(length(lower) != length(upper)){
    stop("the length of lower vector must be equal to the length of upper vector")
  }
  if(length(lower) == 1 && length(upper) == 1){
    lower = rep(lower,2)
    upper = rep(upper,2)
  }
  if(length(hessian) == 1){hessian = rep(hessian, 2)}

  scaleSearch <- match.arg(scaleSearch)
  twoWaySearch <- match.arg(twoWaySearch)
  stressWay <- match.arg(stressWay)
  stressGroup <- rep(0, ngroups)
  RSSGroup <- rep(0, ngroups)
  TSSGroup <- rep(0, ngroups)
  loss <- 0
  for (gn in 1:ngroups){
    lossFunctionOutput <- lossFunction(gn = gn, combProp = combPropGroup[[gn]], Delta = Delta,
                                       radius = radiusGroup[[gn]], disProp = disPropGroup[[gn]], lambda = lambda,
                                       weight = weightGroup[[gn]], ThreeD = ThreeD,
                                       method = c(scaleSearch, twoWaySearch), twoWayGenerate = twoWayGenerate,
                                       lower=  lower, upper = upper, control = control, hessian = hessian,
                                       expand = expand, scaleSeachTolerance = scaleSeachTolerance,
                                       distanceTolerance = distanceTolerance, lossTolerance = lossTolerance,
                                       stressBound = stressBound, maximumStep = maximumStep, planeSize = planeSize)
    groupCentre[[gn]] <- lossFunctionOutput$centre
    stressGroup[gn] <- lossFunctionOutput$stress
    RSSGroup[gn] <- lossFunctionOutput$RSS
    TSSGroup[gn] <- lossFunctionOutput$TSS
    loss <- loss + lossFunctionOutput$loss
  }
  if(stressWay == "sum"){
    stress <- sum(stressGroup)
  } else {
    #combine stress
    stress <- sum(RSSGroup)/sum(TSSGroup)
  }
  #combine all groups
  if(ngroups > 1){
    combinedGroups <- combineGroups(groupCentre = groupCentre, radiusGroup = radiusGroup,delta = delta)
    xy <- combinedGroups$xy
    radius <- combinedGroups$radius
    oneWaySetName <- combinedGroups$namexy
  } else {
    xy <- groupCentre[[1]]
    radius <- radiusGroup[[1]]
    oneWaySetName <- names(radius)
  }
  if(ThreeD){
    sphere(xy = xy, radius = radius, cols = cols, alpha = alpha, oneWaySetName = oneWaySetName, smooth = smooth)
  } else {
    plot.new()
    plotCircle(xy, radius, oneWaySetName, cols, mar = mar)
  }
  list(xy = xy, radius = radius, loss = loss, stress = stress)
}
#--- helper functions -----------------------------------------------------

#library(rgl)
#library(stringr)
#library(Rcpp)
#sourceCpp(file = "your path/vennCpp2.cpp")

#lossFunction, the main function of vennplot; return centres of each group
lossFunction <- function(gn, combProp, radius, Delta, lambda,disProp,
                         weight, ThreeD, method, twoWayGenerate,
                         lower, upper, control, hessian, expand, scaleSeachTolerance,
                         distanceTolerance, lossTolerance, stressBound, maximumStep, planeSize){
  disjointSetNames <- names(combProp)
  numWays <- str_count(disjointSetNames,pattern = "&")+1
  #one way set
  oneWaySetName <- disjointSetNames[which(numWays==1)]
  m <- length(oneWaySetName)
  if(m == 1){
    if(ThreeD){centre <- matrix(c(0,0,0),nrow = 1)}else{centre <- matrix(c(0,0),nrow=1)}
    rownames(centre) <- oneWaySetName
    loss <- 0
    oneWayStress <- calculateStress(centre, radius, disProp, weight, ThreeD, planeSize,twoWayGenerate)
    RSS <- oneWayStress$RSS
    TSS <- oneWayStress$TSS
    stress <- oneWayStress$stress
  }else{
    complete <- TRUE
    if(max(numWays)>=3){
      firstGenerate <- twoWayGeneration(combProp = combProp, mu = 1)
      #To determine which example we use: Example 1 or Example 2
      if(firstGenerate$resam != 0){complete <- FALSE}
      newTworWaySet <- firstGenerate$newTworWaySet
    }else{newTworWaySet <- combProp[which(numWays == 2)]}
    EDandInitialLocation <- EuclideanDistance(newTworWaySet = newTworWaySet, oneWaySetName = oneWaySetName,
                                              radius = radius, ThreeD = ThreeD, initial = TRUE,
                                              expand = expand, distanceTolerance = distanceTolerance,
                                              maximumStep = maximumStep)
    ED <- EDandInitialLocation$ED
    xy <- EDandInitialLocation$xy
    if(twoWayGenerate == TRUE &&  complete == FALSE){
      out <- solveWithMu(ED = ED, xy = xy, combProp = combProp, ThreeD = ThreeD,
                         Delta = Delta,  radius = radius, disProp = disProp, weight = weight, method = method,
                         firstGenerate = firstGenerate,
                         lower=  lower, upper = upper, control = control, hessian = hessian,
                         expand = expand, scaleSeachTolerance = scaleSeachTolerance,
                         distanceTolerance = distanceTolerance, lossTolerance = lossTolerance,
                         stressBound = stressBound, maximumStep = maximumStep, planeSize = planeSize,
                         twoWayGenerate = twoWayGenerate)
      ED <- out$ED
      xy <- out$xy
    }
    result <- solveWithLambda(ED = ED, xy = xy, ThreeD = ThreeD, lambda = lambda, Delta = Delta,
                              radius = radius, disProp = disProp, weight = weight, method = method,
                              lower=  lower, upper = upper, control = control, hessian = hessian,
                              scaleSeachTolerance = scaleSeachTolerance, lossTolerance = lossTolerance,
                              stressBound = stressBound, maximumStep = maximumStep, planeSize = planeSize,
                              twoWayGenerate = twoWayGenerate)
    centre <- result$xy
    rownames(centre) <- oneWaySetName
    stress <- result$stress
    loss <- result$loss
    RSS <- result$RSS
    TSS <- result$TSS
  }
  list(centre = centre, loss = loss, stress = stress, RSS = RSS, TSS = TSS)
}
#optimize lambda to get centre
solveWithLambda <- function(ED , xy , ThreeD , lambda, Delta, radius ,
                            disProp, weight, method, lower, upper, control,
                            hessian, scaleSeachTolerance, lossTolerance,
                            stressBound, maximumStep, planeSize, twoWayGenerate) {
  stress_0_Calculation <- calculateStress(xy = xy,radius = radius,
                                          disProp = disProp,
                                          weight = weight,ThreeD = ThreeD, planeSize = planeSize,
                                          twoWayGenerate = twoWayGenerate)
  stress_0 <- stress_0_Calculation$stress
  lambda_1_Stress <- findOptimalStress(lambda = 1, xy = xy,weight = weight,
                                       radius = radius, disProp = disProp, ED = ED,
                                       ThreeD = ThreeD, lossTolerance = lossTolerance, planeSize = planeSize,
                                       twoWayGenerate = twoWayGenerate)
  stress_1 <- lambda_1_Stress$stress
  if(ThreeD){offset <- 3}else{offset <- 2}
  if(min(stress_0,stress_1) == stress_0 && dim(xy)[1]<=offset){
    stress <- stress_0
    loss <- 0
    RSS <- stress_0_Calculation$RSS
    TSS <- stress_0_Calculation$TSS
  }else{
    if(is.null(lambda)){
      if(method[1] == "NelderMead"){
        Centre <- list();
        lambda_2_Stress <- findOptimalStress(lambda = 1+Delta, xy = xy, weight = weight,
                                             radius = radius, disProp = disProp,
                                             ED = ED, ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        stress_2 <- lambda_2_Stress$stress
        lambda_3_Stress <- findOptimalStress(lambda = 1-Delta, xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        stress_3 <- lambda_3_Stress$stress
        Centre[[1]] <- lambda_1_Stress$xy
        Centre[[2]] <- lambda_2_Stress$xy
        Centre[[3]] <- lambda_3_Stress$xy
        STRESS <- c(stress_1,stress_2,stress_3)
        xy <- Centre[[which(STRESS == min(STRESS))[1]]]
        loss <- c(lambda_1_Stress$loss, lambda_2_Stress$loss, lambda_3_Stress$loss)[which(STRESS == min(STRESS))[1]]
        LAMBDA <- c(1,1+Delta,1-Delta)[order(STRESS)]
        STRESS <- sort(STRESS)
        count <- 0
        while (count< maximumStep && min(STRESS) > stressBound){
          count <- count+1
          lambda_0 <- mean(LAMBDA[1:2])
          lambda_R <- lambda_0 + (lambda_0 - LAMBDA[3])
          lambda_R_stress <- findOptimalStress(lambda = lambda_R, xy = xy, weight = weight,
                                               radius = radius, disProp = disProp, ED = ED,
                                               ThreeD = ThreeD, lossTolerance = lossTolerance,
                                               planeSize = planeSize, twoWayGenerate = twoWayGenerate)
          stress_R <- lambda_R_stress$stress
          #Reflection:
          if(STRESS[1] <= stress_R &&  stress_R< STRESS[2]){
            STRESS <- c(STRESS[1] , stress_R , STRESS[2])
            LAMBDA <- c(LAMBDA[1], lambda_R, LAMBDA[2])
          } else if (stress_R < STRESS[1]){
            #Expansion
            xy <- lambda_R_stress$xy
            loss <- lambda_R_stress$loss
            lambda_E <- lambda_0 + 2*(lambda_R - lambda_0)
            lambda_E_stress <- findOptimalStress(lambda = lambda_E, xy = xy, weight = weight,
                                                 radius = radius, disProp = disProp, ED = ED,
                                                 ThreeD = ThreeD, lossTolerance = lossTolerance,
                                                 planeSize = planeSize, twoWayGenerate = twoWayGenerate)
            stress_E = lambda_E_stress$stress
            if(stress_E<stress_R){
              STRESS <- c(stress_E , stress_R , STRESS[1])
              LAMBDA <- c(lambda_E, lambda_R, LAMBDA[1])
              xy <- lambda_E_stress$xy
              loss <- lambda_E_stress$loss
            }else if(stress_R<=stress_E && stress_E <STRESS[1]){
              STRESS <- c(stress_R, stress_E, STRESS[1])
              LAMBDA <- c(lambda_R, lambda_E, LAMBDA[1])
            }
            else if(STRESS[1]<=stress_E && stress_E <STRESS[2]){
              STRESS <- c(stress_R, STRESS[1],stress_E)
              LAMBDA <- c(lambda_R, LAMBDA[1], lambda_E)
            }else{
              STRESS <- c(stress_R, STRESS[1],STRESS[2])
              LAMBDA <- c(lambda_R, LAMBDA[1],LAMBDA[2])
            }
          } else if(stress_R >= STRESS[2]) {
            #Contraction
            lambda_C <-  lambda_0 + 0.5*(LAMBDA[3] - lambda_0)
            lambda_C_stress <- findOptimalStress(lambda = lambda_C, xy = xy, weight = weight,
                                                 radius = radius, disProp = disProp, ED = ED,
                                                 ThreeD = ThreeD, lossTolerance = lossTolerance,
                                                 planeSize = planeSize, twoWayGenerate = twoWayGenerate)
            stress_C <- lambda_C_stress$stress
            if (stress_C < STRESS[3]){
              STRESS <- c(STRESS[1],STRESS[2],stress_C)
              LAMBDA <- c(LAMBDA[1],LAMBDA[2],lambda_C)[order(STRESS)]
              STRESS <- sort(STRESS)
              if(min(STRESS) == stress_C){
                xy <- lambda_C_stress$xy
                loss <- lambda_C_stress$loss
              }
            } else {
              #shrink
              LAMBDA[2] <- (LAMBDA[1]+LAMBDA[2])/2
              STRESS[2] <- findOptimalStress(lambda = LAMBDA[2], xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize,twoWayGenerate = twoWayGenerate)$stress
              LAMBDA[3] <- (LAMBDA[1]+LAMBDA[3])/2
              STRESS[3] <- findOptimalStress(lambda = LAMBDA[3], xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize,twoWayGenerate = twoWayGenerate)$stress
            }
          }
          if(allConnectedCpp(xy = xy, radius = radius, ThreeD = ThreeD) == FALSE){break}
          if(BoolScaleNMCpp( proportional = scaleSeachTolerance$ proportional,
                             value = scaleSeachTolerance$value,
                             LAMBDA = LAMBDA, STRESS = STRESS) == FALSE){break}
        }
        if(allConnectedCpp(xy = xy, radius = radius, ThreeD = ThreeD) == FALSE) {
          lambdaStress <-  findOptimalStress(lambda = LAMBDA[2], xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
          xy <- lambdaStress$xy
          if(allConnectedCpp(xy = xy, radius = radius, ThreeD = ThreeD) == FALSE) {
            lambdaStress <-  findOptimalStress(lambda = LAMBDA[3], xy = xy, weight = weight,
                                               radius = radius, disProp = disProp, ED = ED,
                                               ThreeD = ThreeD, lossTolerance = lossTolerance,
                                               planeSize = planeSize, twoWayGenerate = twoWayGenerate)
            xy <- lambdaStress$xy
          }
          stress <- lambdaStress$stress
          loss <- lambdaStress$loss
          RSS <- lambdaStress$RSS
          TSS <- lambdaStress$TSS
        } else {
          lambdaStress <-  findOptimalStress(lambda = LAMBDA[1], xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
          xy <- lambdaStress$xy
          stress <- lambdaStress$stress
          loss <- lambdaStress$loss
          RSS <- lambdaStress$RSS
          TSS <- lambdaStress$TSS
        }
      } else if (method[1] == "lineSearch"){
        lambda_2_Stress <- findOptimalStress(lambda = 1+Delta, xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        stress_2 <- lambda_2_Stress$stress
        lambda_3_Stress <- findOptimalStress(lambda = 1-Delta, xy = xy, weight = weight,
                                             radius = radius, disProp = disProp, ED = ED,
                                             ThreeD = ThreeD, lossTolerance = lossTolerance,
                                             planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        stress_3 <- lambda_3_Stress$stress
        if(min(stress_1,stress_2,stress_3) == stress_1){
          xy <- lambda_1_Stress$xy
          loss <- lambda_1_Stress$loss
          stress <- stress_1
          RSS <- lambda_1_Stress$RSS
          TSS <- lambda_1_Stress$TSS
        } else {
          Center <- list()
          Loss <- c()
          if(min(stress_1,stress_2,stress_3) == stress_2){
            stress_n <- stress_2
            stress <- stress_2
            Center[[1]] <- lambda_2_Stress$xy
            Loss[1] <- lambda_2_Stress$loss
            lambda <- 1+Delta
            shrinkage <- TRUE
          } else {
            stress_n <- stress_3
            stress <- stress_3
            Center[[1]] <- lambda_3_Stress$xy
            Loss[1] <- lambda_3_Stress$loss
            lambda <- 1-Delta
            shrinkage <- FALSE
          }
          RSSVec <- c()
          TSSVec <- c()
          n <- 1
          while(stress_n <= stress){
            stress <- stress_n
            if(shrinkage){
              lambda <- lambda + Delta
            } else {
              lambda <- lambda - Delta
            }
            n <- n+1
            lambdaStress <- findOptimalStress(lambda = lambda, xy = Center[[n-1]], weight = weight,
                                              radius = radius, disProp = disProp, ED = ED,
                                              ThreeD = ThreeD, lossTolerance = lossTolerance,
                                              planeSize = planeSize, twoWayGenerate = twoWayGenerate)
            Center[[n]] <- lambdaStress$xy
            Loss[n] <- lambdaStress$loss
            stress_n <- lambdaStress$stress
            RSSVec[n] <- lambdaStress$RSS
            TSSVec[n] <- lambdaStress$TSS
            if(BoolScaleLCpp( proportional = scaleSeachTolerance$ proportional,
                              value = scaleSeachTolerance$value, stress_n = stress_n, stress = stress)){break}
          }
          xy <- Center[[n-1]]
          loss <- Loss[n-1]
          RSS <- RSSVec[n-1]
          TSS <- TSSVec[n-1]
        }
      }else if(method[1] == "goldenSectionSearch"){

        Optimization <- optimize(f = findOptimalStressLambda, lower=  lower[1], upper = upper[1],
                                 xy = lambda_1_Stress$xy, weight = weight,
                                 radius = radius, disProp = disProp, ED = ED,
                                 ThreeD = ThreeD, lossTolerance = lossTolerance,
                                 planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        if(is.null(Optimization$minimum)||is.infinite(Optimization$minimum)){
          stop("This method does not converge")
        }
        lambdaStress <- findOptimalStress(lambda = Optimization$minimum, xy = lambda_1_Stress$xy,weight = weight,
                                          radius = radius, disProp = disProp, ED = ED,
                                          ThreeD = ThreeD, lossTolerance = lossTolerance,
                                          planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        xy <- lambdaStress$xy
        loss <- lambdaStress$loss
        stress <- Optimization$objective
        RSS <- lambdaStress$RSS
        TSS <- lambdaStress$TSS
      } else {
        ## method[1] must be an method appropriate for optim(...)
        Optimization <- optim(par = 1, fn = findOptimalStressLambda,method = method[1],
                              lower = lower[1], upper = upper[1],
                              control = control, hessian = hessian[1],
                              xy = lambda_1_Stress$xy, weight = weight,
                              radius = radius, disProp = disProp, ED = ED,
                              ThreeD = ThreeD, lossTolerance = lossTolerance,
                              planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        if(is.null(Optimization$par)||is.infinite(Optimization$par)){
          stop("This method does not converge")
        }
        lambdaStress <- findOptimalStress(lambda = Optimization$par, xy = lambda_1_Stress$xy,weight = weight,
                                          radius = radius, disProp = disProp, ED = ED,
                                          ThreeD = ThreeD, lossTolerance = lossTolerance,
                                          planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        xy <- lambdaStress$xy
        loss <- lambdaStress$loss
        stress <- Optimization$value
        RSS <- lambdaStress$RSS
        TSS <- lambdaStress$TSS
      }
    } else {
      Centre <- list()
      RSSVec <- c()
      TSSVec <- c()
      STRESS <- rep(0,length(lambda))
      Loss <- rep(1,length(lambda))
      for(i in 1:length(lambda)){
        out <- findOptimalStress(lambda = lambda[i], xy = xy, weight = weight,
                                 radius = radius, disProp = disProp, ED = ED,
                                 ThreeD = ThreeD, lossTolerance = lossTolerance,
                                 planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        Centre[[i]] <- out$xy
        Loss[i] <- out$loss
        STRESS[i] <- out$stress
        RSSVec[i] <- out$RSS
        TSSVec[i] <- out$TSS
      }
      index <- which(STRESS == min(STRESS))[1]
      stress <- STRESS[index]
      xy <- Centre[[index]]
      loss <- Loss[index]
      RSS <- RSSVec[index]
      TSS <- TSSVec[index]
    }
  }
  list(xy = xy, loss = loss, stress = stress, RSS = RSS, TSS = TSS)
}
# if two way intersections missing, optimize mu to get Euclidean ditance
solveWithMu <- function(ED, xy, combProp, ThreeD, radius, disProp,
                        weight, Delta, method, firstGenerate,
                        lower, upper, control, hessian, expand, scaleSeachTolerance,
                        distanceTolerance, lossTolerance,
                        stressBound, maximumStep, planeSize, twoWayGenerate) {
  stress_0 <- calculateStress(xy = xy,radius = radius,disProp = disProp,weight = weight,ThreeD = ThreeD,
                              planeSize = planeSize, twoWayGenerate = twoWayGenerate)$stress
  mu_1_Stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight,
                                   radius = radius, disProp = disProp, ED = ED,
                                   ThreeD = ThreeD, lossTolerance = lossTolerance,
                                   planeSize = planeSize, twoWayGenerate = twoWayGenerate)
  stress_1 <- mu_1_Stress$stress
  if(min(stress_0,stress_1) != stress_0){
    if(method[2] == "NelderMead"){
      EDhat <- list()
      Centre <- list()
      EDplus <- newEuclideanDistance(firstGenerate = firstGenerate, mu = 1+Delta,
                                     combProp = combProp, radius = radius, ThreeD = ThreeD,
                                     expand = expand, distanceTolerance = distanceTolerance,
                                     maximumStep = maximumStep)
      mu_2_Stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight,radius = radius, disProp = disProp,
                                       ED = EDplus, ThreeD = ThreeD, lossTolerance = lossTolerance,
                                       planeSize = planeSize, twoWayGenerate = twoWayGenerate)
      EDplusplus <- newEuclideanDistance(firstGenerate = firstGenerate, mu = 1+2*Delta,
                                         combProp = combProp, radius = radius, ThreeD = ThreeD,
                                         expand = expand, distanceTolerance = distanceTolerance,
                                         maximumStep = maximumStep)
      mu_3_Stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight,radius = radius, disProp = disProp,
                                       ED = EDplusplus, ThreeD = ThreeD, lossTolerance = lossTolerance,
                                       planeSize = planeSize, twoWayGenerate = twoWayGenerate)
      Centre[[1]] <- mu_1_Stress$xy
      Centre[[2]] <- mu_2_Stress$xy
      Centre[[3]] <- mu_3_Stress$xy
      EDhat[[1]] <- ED
      EDhat[[2]] <- EDplus
      EDhat[[3]] <- EDplusplus
      STRESS <- c(stress_1, mu_2_Stress$stress, mu_3_Stress$stress)
      ED <- EDhat[[which(STRESS == min(STRESS))[1]]]
      centre <- Centre[[which(STRESS == min(STRESS))[1]]]
      MU <- c(1,1+Delta,1+2*Delta)[order(STRESS)]
      STRESS <- sort(STRESS)
      count <- 0
      while(count<maximumStep && min(STRESS)> stressBound) {
        #reflection
        count <- count+1
        mu_0 <- mean(MU[1:2])
        mu_R <- mu_0 + (mu_0 - MU[3])
        EDreflection <- newEuclideanDistance(firstGenerate = firstGenerate, mu = mu_R,
                                             combProp = combProp, radius = radius, ThreeD = ThreeD,
                                             expand = expand,distanceTolerance = distanceTolerance,
                                             maximumStep = maximumStep)
        mu_R_stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight, radius = radius, disProp = disProp,
                                         ED = EDreflection,ThreeD = ThreeD, lossTolerance = lossTolerance,
                                         planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        stress_R <- mu_R_stress$stress
        #Reflection:
        if(STRESS[1] <= stress_R &&  stress_R< STRESS[2]){
          STRESS <- c(STRESS[1] , stress_R , STRESS[2])
          MU <- c(MU[1], mu_R, MU[2])
        }else if(stress_R < STRESS[1]){
          #Expansion
          centre <- mu_R_stress$xy
          ED <- EDreflection
          mu_E <- mu_0 + 2*(mu_R - mu_0)
          EDexpansion <- newEuclideanDistance(firstGenerate = firstGenerate, mu = mu_E,
                                              combProp = combProp, radius = radius, ThreeD = ThreeD,
                                              expand = expand, distanceTolerance = distanceTolerance,
                                              maximumStep = maximumStep)
          mu_E_stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight,
                                           radius = radius, disProp = disProp, ED = EDexpansion,
                                           ThreeD = ThreeD, lossTolerance = lossTolerance,
                                           planeSize = planeSize, twoWayGenerate = twoWayGenerate)
          stress_E <- mu_E_stress$stress
          if(stress_E<stress_R){
            STRESS <- c(stress_E , stress_R , STRESS[1])
            MU <- c(mu_E, mu_R, MU[1])
            ED <- EDexpansion
            centre <- mu_E_stress$xy
          }else if(stress_R<=stress_E && stress_E <STRESS[1]){
            STRESS <- c(stress_R, stress_E, STRESS[1])
            MU <- c(mu_R, mu_E, MU[1])
          }
          else if(STRESS[1]<=stress_E && stress_E <STRESS[2]){
            STRESS <- c(stress_R, STRESS[1],stress_E)
            MU <- c(mu_R, MU[1], mu_E)
          }else{
            STRESS <- c(stress_R, STRESS[1],STRESS[2])
            MU <- c(mu_R, MU[1],MU[2])
          }
        }else if(stress_R >= STRESS[2]){
          #Contraction
          mu_C <-  mu_0 + 0.5*(MU[3] - mu_0)
          EDcontraction <- newEuclideanDistance(firstGenerate = firstGenerate, mu = mu_C,
                                                combProp = combProp,
                                                radius = radius, ThreeD = ThreeD,
                                                expand = expand, distanceTolerance = distanceTolerance,
                                                maximumStep = maximumStep)
          mu_C_stress <- findOptimalStress(lambda = 1, xy = xy, weight = weight,radius = radius, disProp = disProp,
                                           ED = EDcontraction, ThreeD = ThreeD, lossTolerance = lossTolerance,
                                           planeSize = planeSize, twoWayGenerate = twoWayGenerate)
          stress_C <- mu_C_stress$stress
          if (stress_C < STRESS[3]){
            STRESS <- c(STRESS[1],STRESS[2],stress_C)
            MU <- c(MU[1],MU[2],mu_E)[order(STRESS)]
            STRESS <- sort(STRESS)
            if(min(STRESS) == stress_C){
              ED <- EDcontraction
              centre <- mu_C_stress$xy}
          }else{
            #shrink
            MU[2] <- (MU[1] + MU[2])/2
            MU[3] <- (MU[1] + MU[3])/2
          }
        }
        if(BoolScaleNMCpp( proportional = scaleSeachTolerance$ proportional,
                           value = scaleSeachTolerance$value,
                           LAMBDA = MU, STRESS = STRESS) == FALSE){break}
      }
      xy <- centre
    } else if(method[2] == "lineSearch") {
      stress_n <- stress_1
      Center <- list()
      Center[[1]] <- mu_1_Stress$xy
      n <- 1
      mu <- 1
      EDhat <- list()
      EDhat[[1]] <- ED
      while (stress_n <= stress_1){
        stress_1 <- stress_n
        mu <- mu + Delta
        ED <- newEuclideanDistance(firstGenerate = firstGenerate, mu = mu,
                                   combProp = combProp, radius = radius, ThreeD = ThreeD,
                                   expand = expand, distanceTolerance = distanceTolerance,
                                   maximumStep = maximumStep)
        L <- findOptimalStress(lambda = 1, xy = Center[[n]], weight = weight,
                               radius = radius, disProp = disProp, ED = ED,
                               ThreeD = ThreeD, lossTolerance = lossTolerance,
                               planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        n <- n+1
        stress_n <- L$stress
        Center[[n]] <- L$xy
        EDhat[[n]] <- ED
      }
      xy <- Center[[n-1]]
      ED <- EDhat[[n-1]]
    } else if(method[2] == "goldenSectionSearch"){
      Optimization <- optimize(f = findOptimalStressMu, lower = lower[2], upper = upper[2],
                               xy = mu_1_Stress$xy, weight = weight,
                               firstGenerate = firstGenerate,combProp = combProp,
                               radius = radius, disProp = disProp, ED = ED,ThreeD = ThreeD,
                               expand = expand, lossTolerance = lossTolerance,
                               distanceTolerance = distanceTolerance,
                               maximumStep = maximumStep,
                               planeSize = planeSize, twoWayGenerate = twoWayGenerate)
      if(is.null(Optimization$minimum) || is.infinite(Optimization$minimum)){
        stop("This method does not converge")
      }else{
        ED <- newEuclideanDistance(firstGenerate = firstGenerate, mu = Optimization$minimum,
                                   combProp = combProp, radius = radius, ThreeD = ThreeD,
                                   expand = expand, distanceTolerance = distanceTolerance,
                                   maximumStep = maximumStep)
        muStress <- findOptimalStress(lambda = 1, xy = mu_1_Stress$xy,weight = weight,
                                      radius = radius, disProp = disProp, ED = ED,
                                      ThreeD = ThreeD, lossTolerance = lossTolerance,
                                      planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        xy <- muStress$xy
      }
    }else{
      Optimization <- optim(par = 1, fn = findOptimalStressMu,method = method[2],
                            lower = lower[2], upper = upper[2],
                            control = control, hessian = hessian[2],
                            xy = mu_1_Stress$xy, weight = weight,
                            firstGenerate = firstGenerate,combProp = combProp,
                            radius = radius, disProp = disProp, ED = ED,ThreeD = ThreeD,
                            expand = expand, lossTolerance = lossTolerance,
                            distanceTolerance = distanceTolerance,
                            maximumStep = maximumStep,
                            planeSize = planeSize, twoWayGenerate = twoWayGenerate)
      if(is.null(Optimization$apr) || is.infinite(Optimization$apr)){
        stop("This method does not converge")
      }else{
        ED <- newEuclideanDistance(firstGenerate = firstGenerate, mu = Optimization$apr,
                                   combProp = combProp, radius = radius, ThreeD = ThreeD,
                                   expand = expand, distanceTolerance = distanceTolerance,
                                   maximumStep = maximumStep)
        L <- findOptimalStress(lambda = 1, xy = mu_1_Stress$xy,weight = weight,
                               radius = radius, disProp = disProp, ED = ED,
                               ThreeD = ThreeD, lossTolerance = lossTolerance,
                               planeSize = planeSize, twoWayGenerate = twoWayGenerate)
        xy <- L$xy
      }
    }
  }
  list(xy = xy, ED = ED)
}
# draw sphere
sphere = function(xy, radius, cols, alpha, oneWaySetName, smooth){
  open3d()
  if(smooth){
    n <- dim(xy)[1]
    f <- function(s, t) cbind(r * cos(s) * cos(t) + x0,
                              r * sin(s) * cos(t) + y0,
                              r * sin(t) + z0)
    for(i in 1:n){
      x0 = xy[i,1]
      y0 = xy[i,2]
      z0 = xy[i,3]
      r = radius[i]
      persp3d(f, slim = c(0, pi), tlim = c(0, 2*pi), n = 101, add = T, col = cols[i], alpha = alpha)
    }
  } else {
    spheres3d(xy[,1], xy[,2], xy[,3], radius = radius,
              color = cols, alpha = alpha)
  }
  text3d(xy[,1], xy[,2], xy[,3], oneWaySetName)
}
# if two way intersections missing, given mu to get a new Euclidean distance
newEuclideanDistance <- function(firstGenerate, mu, combProp, radius, ThreeD, expand,
                                 distanceTolerance, maximumStep){
  disjointSetNames <- names(combProp)
  numWays <- str_count(disjointSetNames,pattern = "&")+1
  #one way set
  oneWaySetName <- disjointSetNames[which(numWays==1)]
  oneWaySet <- combProp[which(numWays==1)]
  nextGenerate <- lapply(firstGenerate$New,
                         function(a, oneWaySet, mu = mu){
                           higherWay <- a[1];newGenerate <- a[-1]
                           newGenerateName <- names(newGenerate)
                           for(i in 1:length(newGenerate)){
                             newGenerate[i] <- higherWay*mu^(str_count(names(higherWay),pattern = "&")+1-2)
                             newGenerateNameSeparate <- str_split(newGenerateName,"&")[[i]]
                             minOneWay <- min(oneWaySet[which((names(oneWaySet) %in% newGenerateNameSeparate) == TRUE)])
                             if(newGenerate[i]> minOneWay){
                               newGenerate[i] <- runif(1,min = higherWay,max = minOneWay)
                             }
                           }
                           newGenerate
                         }, oneWaySet = oneWaySet, mu = mu)
  newTworWaySet <- c(combProp[which(numWays == 2)], unlist(nextGenerate))
  EuclideanDistance(newTworWaySet = newTworWaySet,
                    oneWaySetName = oneWaySetName, radius = radius,
                    ThreeD = ThreeD, initial = FALSE, expand = expand,
                    distanceTolerance = distanceTolerance, maximumStep = maximumStep)
}

# given lambda, optimize stress to get centres and corresponding values
findOptimalStress <- function(lambda, xy, weight, radius, disProp, ED, ThreeD, lossTolerance, planeSize, twoWayGenerate){
  if(is.null(lossTolerance$ALPHA)){
    bool <- TRUE
    #Just satisfy ``double'' input in Cpp; ALPHA will be generated through Newton method
    ALPHA <- 0.01
  } else{
    bool <- FALSE
    ALPHA <- lossTolerance$ALPHA
  }
  L <- lossCpp(xy = xy,radius = radius, lambda = lambda, ED = ED,
               ThreeD = ThreeD, ToleranceofLoss = lossTolerance$ToleranceofLoss, maximumStep = lossTolerance$maximumStep,
               ToleranceofStepsize = lossTolerance$ToleranceofStepsize,  proportional = lossTolerance$ proportional,
               ALPHA = ALPHA, Bool = bool)
  stressCalculation <- calculateStress(xy = L$xy,radius = radius,
                                       disProp = disProp, weight = weight, ThreeD = ThreeD,
                                       planeSize = planeSize, twoWayGenerate = twoWayGenerate)
  stress <- stressCalculation$stress
  RSS <- stressCalculation$RSS
  TSS <- stressCalculation$TSS
  ## Return a list
  list(xy = L$xy, stress = stress, loss = L$loss, RSS = RSS, TSS = TSS)
}
# used in optim(...) or optimization(...) function, optimize lambda with minimum stress
findOptimalStressLambda<- function(lambda, xy, weight, radius, disProp, ED, ThreeD, lossTolerance, planeSize, twoWayGenerate){
  stressValues <- findOptimalStress(lambda = lambda, xy = xy,weight = weight,
                                    radius = radius, disProp = disProp,
                                    ED = ED,ThreeD = ThreeD, lossTolerance = lossTolerance,
                                    planeSize = planeSize, twoWayGenerate = twoWayGenerate)
  ## Return the stress only
  stressValues$stress

}
# used in optim(...) or optimization(...) function, optimize mu with minimum stress
findOptimalStressMu<- function(mu, xy, weight, firstGenerate, combProp, radius, disProp, ED,
                               ThreeD, expand, lossTolerance, distanceTolerance, maximumStep, planeSize, twoWayGenerate){
  ED <- newEuclideanDistance(firstGenerate = firstGenerate, mu = mu,
                             combProp = combProp, radius = radius, ThreeD = ThreeD,
                             expand = expand, distanceTolerance = distanceTolerance,
                             maximumStep = maximumStep)
  stressValues <- findOptimalStress(lambda = 1, xy = xy, weight = weight,
                                    radius = radius, disProp = disProp, ED = ED,
                                    ThreeD = ThreeD, lossTolerance = lossTolerance,
                                    planeSize = planeSize, twoWayGenerate = twoWayGenerate)
  ## Return the stress only
  stressValues$stress
}

#transform disjointcombinations to not disjoint combinations
disjoint2combinations <- function(combinations){
  disjointSetNames <- names(combinations)
  #numWays <- str_length(name)
  numWays <- str_count(disjointSetNames,pattern = "&")+1
  combinations <- combinations[order(numWays)]
  numWays <- sort(numWays)
  if(max(numWays)!=1){
    for(i in 1:length(numWays)){
      if(numWays[i] == max(numWays)){break} else {
        Index <- sapply(str_split(names(combinations[-i]),pattern = "&"),
                        function(a){bool <- 0;
                        if(all(str_split(names(combinations[i]),pattern = "&")[[1]]%in%a))
                        {bool <- 1};
                        bool})
        if(length(which(Index == 1))!=0){combinations[i] <- combinations[i] + sum(combinations[-i][which(Index == 1)])}
      }
    }
  }

  combinations
}

#reorder
reorderDisjointCombinations <- function(combinations, disjoint.combinations){
  m <- length(combinations)
  newDisjointCombinations <- rep(0, m)
  combinationsName <- names(combinations)
  newOrder <- rep(0, m)
  disjoint.combinationsName <- names(disjoint.combinations)
  for(i in 1:m){
    wayOrder <- which(disjoint.combinationsName%in%combinationsName[i] == TRUE)
    newDisjointCombinations[i] <- disjoint.combinations[wayOrder]
    newOrder[i] <- wayOrder
  }
  names(newDisjointCombinations) <- combinationsName
  list(newDisjointCombinations = newDisjointCombinations, newOrder = newOrder)
}

# If input is a data list(frame), extract combinations from it
extractCombinations <- function(data, vars) {
  data <- as.data.frame(data)
  #turn character to numeric
  if(is.null(vars)==FALSE){
    allColName <- colnames(data)
    vars <- match.arg(vars , allColName,several.ok = T)
    data <- data[, which(allColName%in%vars) ]
  }
  if(is.null(dim(data))){
    data <- as.data.frame(data)
    colnames(data) <- vars
  }
  if(dim(data)[2] == 1){
    warning("Meaningless factor data frame")
  }
  colReduceForSure <- rep(0,dim(data)[2])
  for(i in 1:dim(data)[2]){
    if(is.numeric(data[,i]) && all(unique(data[,i])%in%c(0,1)) == FALSE){
      colReduceForSure[i] <- i
    } else if(length(unique(data[,i])) == 1){
      colReduceForSure[i] <- i
    }
  }
  # get rid of some numeric columns
  if(all(colReduceForSure==0)==FALSE){
    colReduceForSure <- colReduceForSure[which(colReduceForSure!=0)]
    data <- data[,-colReduceForSure]
    warning(cat(paste("Non-factor column(s)",toString(colReduceForSure),
                      "has(have) been ignored"), "\n"))
  }
  colName <- colnames(data)
  n <- dim(data)[1]
  p <- dim(data)[2]
  colAdd <- list()
  colReduce <- rep(0,p)
  for(i in 1:p){
    uniqueName <- unique(data[,i])
    if(all(uniqueName%in%c(0,1))) {next
    } else if (mode(levels(data[,i])) == "character"){
      colReduce[i] <- i
      newDataSet <- matrix(rep(0,n*length(uniqueName)),nrow = n)
      for(j in 1:length(uniqueName)){
        newDataSet[which(data[,i] == uniqueName[j]),j] <- 1
      }
      colnames(newDataSet) <- uniqueName
      colAdd[[i]] <- newDataSet
    }
  }
  colReduce <- colReduce[which(colReduce!=0)]
  if(length(colReduce)!=0) {
    newdata <- as.matrix(data[,-colReduce])
    if(dim(newdata)[2] == 1){
      colnames(newdata) <- colName[-colReduce]
    }
    for(i in 1:length(colAdd)){
      colAddName <- colnames(colAdd[[i]])
      sumAll <- apply(colAdd[[i]], 2, "sum")
      deleteCol <- which(sumAll == min(sumAll))[1]
      newInput <- as.matrix(colAdd[[i]][,-deleteCol])
      if(dim(newInput)[2] == 1){
        colnames(newInput) <- colAddName[-deleteCol]
      }
      newdata <- cbind(newdata, newInput)
    }
  } else {newdata <- data}
  rowSumZero <- which(apply(newdata,1,"sum")==0)
  if(length(rowSumZero)!=0){
    newdata <- as.matrix(newdata[-which(apply(newdata,1,"sum")==0),])
  }
  if(dim(newdata)[2] == 1) {
    disjoint.combinations <- c(sum(newdata))
    names(disjoint.combinations) <- colnames(newdata)
  } else {
    #OUTCOME is not disjoint
    G <- list()
    for(i in 1:dim(newdata)[2]){
      G[[i]] <- t(combn(dim(newdata)[2],i))
    }
    #OUTCOME is disjoint
    newColName <- colnames(newdata)
    nameList <- apply(newdata,1, function(a){paste(newColName[which(a==1)],collapse = "&")})
    disjointOutput <- aggregate(data.frame(count = nameList), list(name = nameList), length)
    disjoint.combinations <- disjointOutput[,2]
    names(disjoint.combinations) <- disjointOutput[,1]
    numWays <- str_count(disjointOutput[,1], pattern = "&")+1
    oneWays <- which(numWays == 1)
    if(length(oneWays) != length(newColName)){
      notIn <- which(newColName %in% disjointOutput[,1][oneWays] == FALSE)
      newIn <- rep(0, length(notIn))
      names(newIn) <- newColName[notIn]
      disjoint.combinations <- c(disjoint.combinations, newIn)
    }
  }
  disjoint.combinations
}

# Calculates distance between two circles
Distance <- function(r1,r2,S,ThreeD, expand, distanceTolerance, maximumStep) {
  if (ThreeD){
    if(S == 0){
      if(is.null(expand)) {
        d <- r1+r2
      } else {
        d <- (r1+r2)*expand
      }
    } else if(abs(S - min(4*pi/3*r1^3,4*pi/3*r2^3)) < distanceTolerance$value) {
      d <- max(r1,r2) - min(r1,r2)
    } else {
      theta <- matrix(c(0,0),nrow = 2)
      thetanew <- theta+1
      f1 <- pi/3*r1^3*(1-cos(theta[1]))^2*(2+cos(theta[1])) +
        pi/3*r2^3*(1-cos(theta[2]))^2*(2+cos(theta[2])) - S
      f2 <- r1*sin(theta[1]) - r2*sin(theta[2])
      k <- 0
      while(k < maximumStep){
        theta <- thetanew
        f1 <- pi/3*r1^3*(1-cos(theta[1]))^2*(2+cos(theta[1])) +
          pi/3*r2^3*(1-cos(theta[2]))^2*(2+cos(theta[2])) - S
        f2 <- r1*sin(theta[1]) - r2*sin(theta[2])
        f <- matrix(c(f1,f2),nrow = 2)
        g <- matrix(c(pi*r1^3*sin(theta[1])^3, pi*r2^3*sin(theta[2])^3,
                      r1*cos(theta[1]),-r2*cos(theta[2])),nrow = 2, byrow = T)
        g <- solve(g)
        thetanew <- theta - g%*%f
        k <- k+1
        if(BoolDistanceCpp( proportional = distanceTolerance$ proportional,
                            value = distanceTolerance$value,
                            f1 = f1, f2 = f2,
                            thetanew = thetanew, theta = theta) == FALSE){break}
      }
      if(any(thetanew>pi) || any(thetanew<0) || k == maximumStep){
        theta1 <- seq(0,pi, length = 200)
        theta2 <- seq(0,pi, length = 200)
        searchMatrix <- distanceCpp(r1, r2,theta1, theta2, S,ThreeD)
        index <- which(searchMatrix== min(searchMatrix),arr.ind=T)[1,]
        d <- r1*cos( theta1[index[1]]) + r2*cos(theta2[index[2]])
      } else {
        d <- r1*cos(thetanew[1]) + r2*cos(thetanew[2])
      }
    }
  }
  else{
    if(S==0){
      if(is.null(expand)){
        d <- r1+r2
      } else {
        d <- (r1+r2)*expand
      }
    }else if( abs(S - min(pi*r1^2,pi*r2^2)) < distanceTolerance$value ){
      d <- max(r1,r2) - min(r1,r2)
    }else{
      theta <- matrix(c(0,0),nrow = 2)
      thetanew <- theta+1
      f1 <- theta[1]*r1^2 - sin(2*theta[1])*r1^2/2 +theta[2]*r2^2 - sin(2*theta[2])*r2^2/2  - S
      f2 <- r1*sin(theta[1]) - r2*sin(theta[2])
      k <- 0
      while(k < maximumStep){
        theta <- thetanew
        f1 <- theta[2]*r2^2 - sin(2*theta[2])*r2^2/2 + theta[1]*r1^2 - sin(2*theta[1])*r1^2/2 - S
        f2 <- r1*sin(theta[1]) - r2*sin(theta[2])
        f <- matrix(c(f1,f2),nrow = 2)
        g <- matrix(c(r1^2-r1^2*cos(2*theta[1]), r2^2-r2^2*cos(2*theta[2]),
                      r1*cos(theta[1]),-r2*cos(theta[2])),nrow = 2, byrow = T)
        g <- solve(g)
        thetanew <- theta - g%*%f
        k <- k+1
        if(BoolDistanceCpp( proportional = distanceTolerance$ proportional,
                            value = distanceTolerance$value,
                            f1 = f1, f2 = f2,
                            thetanew = thetanew, theta = theta) == FALSE){break}
      }
      if(any(thetanew>pi) || any(thetanew<0) || k == maximumStep){
        theta1 <- seq(0,pi, length = 200)
        theta2 <- seq(0,pi, length = 200)
        searchMatrix <- distanceCpp(r1, r2,theta1, theta2, S,ThreeD)
        index <- which(searchMatrix== min(searchMatrix),arr.ind=T)[1,]
        d <- r1*cos( theta1[index[1]]) + r2*cos(theta2[index[2]])
      } else {
        d <- r1*cos(thetanew[1]) + r2*cos(thetanew[2])
      }
    }
  }
  return(d)
}
# Plots circles of Venn diagram
plotCircle <-  function(xy, radius,name, col, mar, ...) {
  par(mar = mar)
  a1 <-  range(c((xy + radius)[,1], (xy - radius)[,1]))
  a2 <- range(c((xy + radius)[,2], (xy - radius)[,2]))
  plot.window(a1, a2, "", asp = 1)
  theta <- seq.int(360)/360*2*pi
  for (i in 1:length(radius)){
    polygon(xy[i,1] +  radius[i]*cos(theta), xy[i,2] + radius[i]*sin(theta), col = col[i],border = col[i])
  }
  text(xy, name)
}

#rotate centres until two groups totally separated
rotateCentres <- function(xy,transxy,radius1,radius2, delta){
  for(i in 1:35){
    theta <- i/18*pi
    if(dim(xy)[2] == 3){
      rotation1 <- matrix(c(1,0,0,0,cos(theta),-sin(theta),0,sin(theta),cos(theta)),nrow = 3,byrow = T)
      rotation2 <- matrix(c(cos(theta),0,sin(theta),0,1,0,-sin(theta),0,cos(theta)),nrow = 3,byrow = T)
      rotation3 <- matrix(c(cos(theta),-sin(theta),0,sin(theta),cos(theta),0,0,0,1),nrow = 3,byrow = T)
      rotation <- rotation1%*%rotation2%*%rotation3}
    else{
      rotation <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),nrow = 2,byrow = T)
    }
    newxy <- t(rotation%*%t(transxy))
    newxy <- t(transxy[1,] - newxy[1,] + t(newxy))
    Judgement <- allDisjointCpp(xy,newxy,radius1,radius2, delta)
    if (Judgement!=0 ){break}
  }
  list(Judgement =Judgement, newxy = newxy)
}

# after finding centres, combine them with reasonable distance
combineGroups <- function(groupCentre, radiusGroup, delta){
  groupNum <- length(groupCentre)
  radius <- unlist(radiusGroup)
  namexy <- lapply(radiusGroup,
                   function(a){
                     names(a)
                   })
  groupxy <- groupCentre
  radiusxy <- radiusGroup
  for (gn in 1:(groupNum-1)){
    #Randomly find a minimum x, minimum y, maximum x, or maximum y of the second group
    cornerrandom = ceiling(runif(1,0,4))
    if(cornerrandom == 1){
      corner <- which(groupCentre[[gn+1]][,1] == min(groupCentre[[gn+1]][,1]))[1]
    }else if(cornerrandom == 2){
      corner <- which(groupCentre[[gn+1]][,1] == max(groupCentre[[gn+1]][,1]))[1]
    }else if(cornerrandom == 3){
      corner <- which(groupCentre[[gn+1]][,2] == min(groupCentre[[gn+1]][,2]))[1]
    }else{
      corner <- which(groupCentre[[gn+1]][,2] == max(groupCentre[[gn+1]][,2]))[1]
    }
    #corresponding radius
    radiusvec <- radiusGroup[[gn+1]][corner]
    #generate a centre to check which can make this centre is totally disjoint with the former group
    xyvec <- transCpp(xy = groupCentre[[gn]], radius = radiusGroup[[gn]], radiusvec = radiusvec, radiusall = radius)
    #move the second group
    transxy <- t(t(groupCentre[[gn+1]]) + xyvec - groupCentre[[gn+1]][corner,])
    #If two groups are totally seperated
    if(allDisjointCpp(groupCentre[[gn]],transxy,radiusGroup[[gn]],radiusGroup[[gn+1]],delta) == 0){
      out <- rotateCentres(groupCentre[[gn]],transxy,radiusGroup[[gn]],radiusGroup[[gn+1]],delta)
      Judgement <- out$Judgement
      newxy <- out$newxy
      while(Judgement == 0){
        xyvec <- transCpp(xy = groupCentre[[gn]], radius = radiusGroup[[gn]],
                          radiusvec = radiusvec, radiusall = radius)
        transxy <- t(t(groupCentre[[gn+1]]) + xyvec - groupCentre[[gn+1]][corner,])
        out <- rotateCentres(groupCentre[[gn]],transxy,radius[[gn]],radius[[gn+1]],delta)
        Judgement <- out$Judgement
        newxy <- out$newxy
      }
      rownames(newxy) <- namexy[[gn+1]]
    }else{newxy <- transxy}
    #renew the second group centre
    groupxy[[gn+1]] <- newxy
  }
  #Move the last group to the former one with reasonable distance
  for(gn in 1:(groupNum-1)){
    center1 <- rep(0,dim(groupxy[[gn]])[2])
    center2 <- rep(0,dim(groupxy[[gn+1]])[2])
    for(j in 1:length(center1)){
      center1[j] <- (max(groupxy[[gn]][,j]+radiusxy[[gn]])+
                       min(groupxy[[gn]][,j]-radiusxy[[gn]]))/2
      center2[j] <- (max(groupxy[[gn+1]][,j]+radiusxy[[gn+1]])+
                       min(groupxy[[gn+1]][,j]-radiusxy[[gn+1]]))/2
    }
    direction <- center1-center2
    if(delta==0){delta <- (1e-4)}
    # scale it
    direction <- direction/sqrt(sum(direction^2))*delta/10
    newxy <- closeCpp(groupxy[[gn]],groupxy[[gn+1]],radiusxy[[gn]],
                      radiusxy[[gn+1]],delta = delta,direction)$xy
    groupxy[[gn+1]] <- rbind(groupxy[[gn]],newxy)
    radiusxy[[gn+1]] <- c(radiusxy[[gn]],radiusxy[[gn+1]])
    namexy[[gn+1]] <- c(namexy[[gn]],namexy[[gn+1]])
  }
  list(xy = groupxy[[groupNum]], radius = radiusxy[[groupNum]],namexy = namexy[[groupNum]])
}

#given centres, numerically calculate each disjoint part's area and return stress
calculateStress <- function(xy, radius, disProp, weight, ThreeD, planeSize, twoWayGenerate){
  m <- length(radius)
  l <- list()
  oneWaySetName <- names(radius)
  if(ThreeD){
    xuan <- (max(xy[,1]) - min(xy[,1])+2*max(radius))/planeSize
    yuan <- (max(xy[,2]) - min(xy[,2])+2*max(radius))/planeSize
    zuan <- (max(xy[,3]) - min(xy[,3])+2*max(radius))/planeSize
    myList <- list()
    for(k in 1:m){
      for(i in 1:planeSize){
        myList[[i]] <- matrix(rep(0,planeSize^2),nrow =planeSize)
      }
      l[[k]] <- binaryIndexThreeDCpp(myList, xy, radius,k ,yuan, xuan, zuan, planeSize)
    }
    numericArea <- goThroughPixelThreeDCpp(l, m, planeSize)
  }else{
    xuan <- (max(xy[,1]) - min(xy[,1])+2*max(radius))/planeSize
    yuan <- (max(xy[,2]) - min(xy[,2])+2*max(radius))/planeSize
    for (k in 1:m){
      l[[k]] <- binaryIndexCpp(matrix(rep(0,planeSize^2),nrow =planeSize),xy,radius,k,yuan,xuan, planeSize)
    }
    numericArea <- goThroughPixelCpp(myList =  l, m = m, num = planeSize)
  }
  if(m == 1){
    TSS <- (sum(numericArea)/(planeSize*planeSize))^2*weight
    RSS <- 0
    stress <- 0
  } else {
    numericArea <- numericArea[which(getRidofZeroCpp(numericArea)!=0),]
    numericAreaMatrix <- matrix(unlist(unique(as.data.frame(numericArea))),ncol = m)
    numericAreaVector <- countCpp(numericArea,numericAreaMatrix)
    numericAreaVectorName <- apply(numericAreaMatrix,1,
                                   function(a, oneWaySetName)
                                   {paste(oneWaySetName[which(a!=0)],collapse="&")
                                   }, oneWaySetName = oneWaySetName)
    names(numericAreaVector) <- numericAreaVectorName
    lengthDisProp <- length(disProp)
    cstar <- rep(0, lengthDisProp)
    astar <- disProp
    wstar <- weight
    lengthNumericAreaVector <- length(numericAreaVector)
    numericAreaVectorIndex <- rep(0,lengthDisProp)
    numericAreaVectorNameSplit <- str_split(numericAreaVectorName,pattern = "&")
    for(i in 1:lengthDisProp){
      splitDisPropName <- str_split(names(disProp),pattern = "&")[[i]]
      index <- 0
      for (j in 1:lengthNumericAreaVector){
        if(j %in% numericAreaVectorIndex) {next}
        if(all(numericAreaVectorNameSplit[[j]] %in% splitDisPropName,
               length(numericAreaVectorNameSplit[[j]]) == length(splitDisPropName))){
          index <- j
          numericAreaVectorIndex[i] <- index
          break
        }
      }
      if(index != 0){
        cstar[i] <- numericAreaVector[index]
      }
    }
    unnecessaryOverlayIndex <- which(c(1:lengthNumericAreaVector) %in% numericAreaVectorIndex == FALSE)
    lengthofThisIndex <- length(unnecessaryOverlayIndex)
    if(lengthofThisIndex != 0 && twoWayGenerate == FALSE){
      cstar <- c(cstar, numericAreaVector[unnecessaryOverlayIndex])
      astar <- c(disProp, rep(0, lengthofThisIndex))
      wstar <- c(weight, rep(1, lengthofThisIndex))
    }
    fit <- lm(cstar~astar-1, weights = wstar)
    RSS <- sum(fit$residuals^2*wstar)
    TSS <- sum(cstar^2*wstar)
    stress <- RSS/TSS
  }
  list(stress = stress, RSS = RSS, TSS = TSS)
}

#given combinations, detect groups
groupDetection <- function(largerThanOneWaySetName, oneWaySetName){
  vertex <- list()
  if(is.null(largerThanOneWaySetName)){
    for(i in 1:length(oneWaySetName)){
      vertex[[i]] <- i
    }
  } else {
    largerThanOneWaySetLength <- length(largerThanOneWaySetName)
    for(i in 1:largerThanOneWaySetLength){
      Boolean1 <- rep(FALSE, largerThanOneWaySetLength)
      for(j in i:largerThanOneWaySetLength){
        Boolean1[j] <- any((str_split(largerThanOneWaySetName[i],"&")[[1]] %in% str_split(largerThanOneWaySetName[j],"&")[[1]])==TRUE)
      }
      Boolean1True <- which(Boolean1 == TRUE)
      vertex[[i]] <- unique(unlist(str_split(largerThanOneWaySetName[Boolean1True],"&")))
      if(i >1){
        Boolean2 <- rep(FALSE, i-1)
        for(j in 1:(i-1)){
          Boolean2[j] <- any((vertex[[i]]%in%vertex[[j]])==TRUE)
        }
        Boolean2True <- which(Boolean2 == TRUE)
        if(length(Boolean2True) == 1){
          vertex[[Boolean2True]] <- unique(c(vertex[[Boolean2True]],vertex[[i]]))
          vertex[[i]] <- NULL
        }
      }
    }
    if(length(which(sapply(vertex, is.null)))!=0){  vertex <- vertex[-which(sapply(vertex, is.null))]}
    vertexLength <- length(vertex)
    if(length(unlist(vertex)) != length(oneWaySetName)){
      difference <- length(oneWaySetName) - length(unlist(vertex))
      rest <- oneWaySetName[which((oneWaySetName%in%unlist(vertex))== FALSE)]
      for(i in 1:difference){
        vertex[[i+vertexLength]] <- rest[i]
      }
    }
    vertex <- lapply(vertex,
                     function(a,oneWaySetName){
                       which((oneWaySetName%in%a)==TRUE)
                     }, oneWaySetName = oneWaySetName)
  }
  vertex
}
# if two way intersections are missing, generate them
twoWayGeneration <- function(combProp, mu){

  #numWays <- str_length(names(combProp))
  numWays <- str_count(names(combProp),pattern = "&")+1
  oneWaySet <- combProp[which(numWays==1)]
  twoWaySet <- combProp[which(numWays == 2)]
  twoWaySetName <- names(twoWaySet)
  #highWay gives way larger than 2
  highWaySet <- sort(combProp[which(numWays > 2)],decreasing = T)
  highWaySetName <- names(highWaySet)

  highWays <- str_count(highWaySetName,pattern = "&")+1
  resam <- 0
  New <- list()
  newName <- c()
  for (i in 1:length(highWays)){
    com <- combn(highWays[i],2)
    L <- dim(com)[2]
    newTwoWaySet <- rep(0,1e5)
    newTwoWaySetName <- rep(0,1e5)
    for (j in 1:L){
      highWaySetNameSplit <- str_split(highWaySetName[i],"&")[[1]][com[,j]]
      value_a <- oneWaySet[which(names(oneWaySet)==highWaySetNameSplit[1])]
      value_b <- oneWaySet[which(names(oneWaySet)==highWaySetNameSplit[2])]
      valuemin <- min(value_a,value_b)
      highWaySetNamePaste <- paste(highWaySetNameSplit,collapse="&")
      highWaySetNamePasteReorder <- paste(c(highWaySetNameSplit[2], highWaySetNameSplit[1]),collapse="&")
      bool1 <- any( highWaySetNamePaste %in% twoWaySetName, highWaySetNamePasteReorder%in%twoWaySetName)
      if(bool1 == FALSE){
        bool2 <- any( highWaySetNamePaste %in% newName, highWaySetNamePasteReorder %in% newName )
        if(bool2 == FALSE){
          resam <- resam+1
          newName[2*resam - 1] <- highWaySetNamePaste
          newName[2*resam] <- highWaySetNamePasteReorder

          if(highWaySet[i] == 0){
            highWaySet[i] <- min(highWaySet[which(highWaySet!=0)])*0.01
          }
          if((highWaySet[i]*mu^(highWays[i]-2)) < valuemin){
            newTwoWaySet[resam] <- (highWaySet[i]*mu^(highWays[i]-2))
          }else{
            newTwoWaySet[resam] <- runif(1,highWaySet[i],valuemin)
          }
          newTwoWaySetName[resam] <- highWaySetNamePaste
        }
      }
    }
    newTwoWayIndex <- which(newTwoWaySet!=0)
    newTwoWaySet <- newTwoWaySet[newTwoWayIndex]
    if(length(newTwoWaySet) != 0){
      newTwoWaySetName <- newTwoWaySetName[newTwoWayIndex]
      names(newTwoWaySet) <- newTwoWaySetName
      New[[i]] <- c(highWaySet[i], newTwoWaySet)
    }
  }
  if(length(New)== 0){
    newTwoWaySet <- twoWaySet
  } else {
    anyNullList <- which(sapply(New, is.null))
    if(length(anyNullList) !=0 ){
      New <- New[-anyNullList]
    }

    newHighWaySet <- unlist(New)
    if(any(is.na(newHighWaySet))){newHighWaySet <- newHighWaySet[-which(is.na(newHighWaySet))]}
    numWays <- str_count(names(newHighWaySet),pattern = "&")+1
    newTwoWaySet <- c(twoWaySet, newHighWaySet[which(numWays == 2)])
  }
  list(newTworWaySet = newTwoWaySet, resam = resam, New = New)
}

#calculate Euclidean distance matrix and intial location (based on Gram matrix)
EuclideanDistance <- function(newTworWaySet, oneWaySetName, radius, ThreeD, initial, expand, distanceTolerance,
                              maximumStep){
  m <- length(oneWaySetName)
  newTworWaySetName <- names(newTworWaySet)
  ED <- matrix(rep(0, m^2),ncol = m)
  for(i in 1:m){
    for(j in i:m){
      if(j == i){next}else{
        sij1 <- newTworWaySet[which(newTworWaySetName == paste(c(oneWaySetName[i],oneWaySetName[j]),collapse = "&"))]
        sij2 <- newTworWaySet[which(newTworWaySetName == paste(c(oneWaySetName[j],oneWaySetName[i]),collapse = "&"))]
        if(length(sij1) == 0 && length(sij2) == 0 ){
          sij <- 0
        } else { sij <-  c(sij1,sij2) }
        ED[i,j] <- Distance(radius[i], radius[j], sij, ThreeD = ThreeD, expand = expand,
                            distanceTolerance = distanceTolerance, maximumStep = maximumStep)
      }
    }
  }
  ED <- ED + t(ED)
  rownames(ED) <- oneWaySetName
  colnames(ED) <- oneWaySetName
  if(initial == TRUE){
    D2 <- ED^2
    J <- diag(m)-1/m
    G <- -0.5*J%*%D2%*%J
    Em <- svd(G)
    U <- (Em$u) %*% sqrt(diag(Em$d))
    if(ThreeD){
      if (dim(U)[1]>=3){xy <- U[,1:3]
      }else{xy <- cbind(U,rep(0,2))}
    }else{
      xy <- U[,1:2]
    }
    step <- 1
    while((allConnectedCpp(xy,radius,ThreeD) == F)){
      meanvec <- apply(xy,2,"mean")
      deriction <- t(t(xy) - meanvec)
      xy <- xy - 0.1*step*deriction
      step <- step + 1
    }
    out <- list(ED = ED, xy = xy)
  }else{out <- ED}
  out
}