Nothing
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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.