R/rhull.R

In blackbox: Black Box Optimization and Exploration of Parameter Space

rExpandedHull <- function(pointnumber, knotshull, hrepr=NULL, vT,extrapol, ## no default, to catch problems
                        focalPoint=NULL) { ## random sampling in a hull generated by 'knots'; but see truer random sampling rhullByvT or simply rvolTriangulation
  ## knotshull should already not be redundant (or hardly so)
  hullmin <- matrixStats::colMins(knotshull, useNames=TRUE) # apply(knotshull, 2, min)
  hullmax <- matrixStats::colMaxs(knotshull, useNames=TRUE) # apply(knotshull, 2, max)
  parnbr <- ncol(knotshull)
  ## rosglobal$par would seem an obvious center but if on an edge, sampled points cannot expand far from it (+ pb with code that samples refpoints from assumedly two points)
  ## use of non-null focalPoint (DISTINCT from bary) allows to sample preferentially a direction from the barycenter   
  if (is.null(focalPoint)) focalPoint <- do.call(blackbox.getOption("barycenterFn"), list(vertices=knotshull)) 
  colNames <- names(focalPoint)
  resu <- as.data.frame(matrix(nrow=0, ncol=parnbr))
  if(is.null(hrepr)) { ## which ideally is not the case
    stop.redef("(!) from rhull: is.null(hrepr). Feasible in principle but probably ill-conceived call to rhull")
    ## otherwise continued with:   (but call to Hreprrational is not very efficient here)
    #hrepr <- Hreprrational(knotshull)
    #hrepr <- q2d(hrepr)
  }
  Lmatrix <- t(chol(cov(knotshull))) ## version ellipse, 2014/02
  for (it in seq(pointnumber)) {
    randomInternalPoint <- rvolTriangulation(n=1L,volTriangulationObj=vT) ## use of randomInternalPoint will allow more genuine extrapol
    if(parnbr>1) {
      ## approximation of sampling in hull by sampling in an ellipse:
      ## random direction:
      randomdangle <- runif(ceiling(parnbr/2))*2*pi
      randomdir <- as.vector(rbind(cos(randomdangle), sin(randomdangle)))[1:parnbr] ## cos(angle1), sin(angle1), cos(angle2...)
      randomdir <- Lmatrix %*% randomdir
      ## finds the intersection of the chull and of the line defined by (focalPoint, direction)
      ## and keeps a defining point (...)[1, 2+seq(parnbr)] :  ## this must be the slow step of the game...
      tmp <- scdd.addHline(fHrepr=hrepr, direction=randomdir, origin=focalPoint)
      # 2015/10/24 non uniform sampling of the ref points introduced to avoid half of the points = focalPoint when focalPoint is at boundary
      if (nrow(tmp)>1) {
        reldist <- tmp[,3]-randomInternalPoint[1]## [,3] to remove cols "a","b"; [1]: any param will give the same result
        prob <- abs(reldist)^parnbr ## sample uniformly wrt to volumes on each side.
        prob <- prob/sum(prob)
        ranpt <- sample(seq(nrow(tmp)), size=1,prob=prob)
      } else ranpt <- 1L
      refpoint <- tmp[ranpt, 2+seq(parnbr)]/tmp[ranpt, 2] ## a point on the envelope... col 2 is not necessarily 1 in the Vrepr!
    } else refpoint <- sample(knotshull, 1)
    ## sampling a distance from center in a sphere uniformly wrt volume:
    dist <- extrapol*(runif(1))^(1/parnbr) ## F^{-1}(u) for sphere of radius 'extrapol'
    ## one problem is that uniformly sampling an (extrapolated) sphere preferentially samples the layers most distant from center
    ## and given how the 'sphere' is defined, this samples low-likelihood layers, missing as inner extrapolated maximum
    resu <- rbind(resu, randomInternalPoint+dist*(refpoint-randomInternalPoint))
  }
  colnames(resu) <- colNames
  return(resu)
} ## end def rhull


## all outputVars must be included in fittedNames (or add fromFONK code )
rhullByEI <- function(n, tryn=100*n ,vT, object, fixed=NULL, outputVars=blackbox.getOption("fittedNames"),blockSize=100L) {
  ## so that oldXnew matrices will contain less than blackbox.getOption("max_mat_size")~1e7 values,
  maxn <- floor(blackbox.getOption("max_mat_size")/nrow(object$data))
  if (maxn <= n) {
    locmess <- paste("From rhullByEI(): 'maxn': ",maxn,"<=",n," ('n'). 'n' reduced to")
    n <- ceiling(maxn/10)
    message.redef(paste(locmess,n))
  }
  if (tryn > maxn) {
    locmess <- paste("From rhullByEI(): 'tryn' reduced from",tryn,"to",maxn)
    message.redef(locmess)
    tryn <- maxn
  }
  trypoints <- data.frame(rvolTriangulation(tryn, vT))
  colnames(trypoints) <- colnames(vT$vertices) ## supposeque non null...
  if (! is.null(fixed)) {
    trypoints <- cbind(trypoints, fixed)
  }
  #if ("Migraine" %in% blackbox.getOption("usedBy")) { ## condition makes sense... rethink later
  trypoints <- apply(trypoints, 1, tofullKrigingspace)
  if (length(outputVars)>1L) {
    trypoints <- t(trypoints)
  } else trypoints <- matrix(trypoints,ncol=1L)
  #}
  colnames(trypoints) <- outputVars ## 'apply' feature
  ## F I X M E j'ai decoupe calcPredVar mais pas point predict(), qui peut aussi être limitant.
  ## il faudrait decouper predict ~ comme ca:
  if ( tryn > blockSize) {
    fulldata <- as.data.frame(trypoints)
    slices <- unique(c(seq(0L,tryn,blockSize),tryn))
    slicefn <- function(it) {
      slice <- (slices[it]+1L):slices[it+1L]
      predict(object, newdata=fulldata[slice,,drop=FALSE],variances=list(predVar=TRUE))
    }
    resus <- lapply(seq_len(length(slices)-1L),slicefn)
    trypred <- do.call(rbind,resus)
    trySE <- do.call(c,lapply(resus,attr,which="predVar"))
  } else {
    trypred <- predict(object, newdata=as.data.frame(trypoints), variances=list(predVar=TRUE))
    trySE <- attr(trypred, "predVar")
  }
  trySE[trySE<0] <- 0
  trySE <- sqrt(trySE)
  ## '-' pcq -ln(L) est prédit
  tryQ <- - trypred + 1.96*trySE ## improvement function for candidate points
  Qmax <- object$Qmax
  expectedImprovement <- trySE*dnorm((Qmax-tryQ)/trySE)+(tryQ-Qmax)*pnorm((tryQ-Qmax)/trySE) ## 7.5 p. 121
  trypoints <- trypoints[order(expectedImprovement, decreasing=TRUE)[seq_len(n)], outputVars, drop=FALSE]
  return(trypoints) ## with names of the restricted space represented by vT
}

rhullByvT <- function(n, vT, fixed=NULL, outputVars=blackbox.getOption("fittedNames")) {
  trypoints <- data.frame(rvolTriangulation(n, vT))
  colnames(trypoints) <- colnames(vT$vertices) ## suppose que non null...
  if (! is.null(fixed)) {
    trypoints <- cbind(trypoints, fixed)
  }
  trypoints <- apply(trypoints, 1, tofullKrigingspace)
  if (length(outputVars)>1L) {
    trypoints <- t(trypoints)
  } else trypoints <- matrix(trypoints,ncol=1L)
  colnames(trypoints) <- outputVars ## 'apply' feature
  return(trypoints) ## with names of the restricted space represented by vT
}