Nothing
## **** 1D profile CI bounds by NLOPT ****
calcBounds1D <- function(CIdlr,## <0
string1, string2, lowval, hival, maxparval, CIvar) {
INFO <- list(fittedNames=blackbox.getOption("fittedNames"),
fittedparamnbr=blackbox.getOption("fittedparamnbr"),
rosglobal=blackbox.getOption("rosglobal"), ## FR->FR dangereux
optimizers=blackbox.getOption("optimizers"),
fitobject=blackbox.getOption("fitobject"),
scalefactor=blackbox.getOption("scalefactor")
)
notinKgspace <- CIvar %w/o% INFO$fittedNames ## e.g. latt2Ns2 %w/o% canonicals is latt2Ns2
locchull <- providefullhull(notinKgspace)[[1]] ## 'full' dimensional, vertices constraints and affectedConstraints
#
if (length(notinKgspace)==0L){ ## no extra composite hull
if(is.na(maxparval)) message.redef("(!) From calcBounds1D(): maxparval is NA were a value is expected")
maxpt <- INFO$rosglobal$par
maxlogL <- INFO$rosglobal$value
} else {
## ML in composite space hull. rosglobal$par may not be in the alternative hull, and a higher max may be found in this hull
if( ! is.na(maxparval)) message.redef("(!) From calcBounds1D(): maxparval has a value were NA is expected")
initpt <- INFO$fitobject$x[which.max(INFO$fitobject$fitted.values), ]
initpt <- fromFONKtoanyspace(initpt, colnames(locchull$vertices))
ui <- -locchull$a;ci <- -locchull$b ## within-hullness occurs for ui %*% x -ci => 0
candidates <- generateInitpts(bestProfiledOut=initpt, vertices=locchull$vertices, ui=ui, ci=ci, hrepr=locchull$Hrep, fixedlist=NA, precision="default")
###################
bestresu <- list(value=- Inf)
for (candidate in candidates) {
resu <- optimWrapper( ##purefn,
initval=candidate$par, gr=NULL, ## not initpt, 14/10/2011
chullformats=locchull,
control=list( ## parscale is provided within optimWrapper
fnscale=-1/INFO$scalefactor, trace=FALSE, maxit=10000))
if (resu$value>bestresu$value) bestresu <- resu
} ## end loop over candidate starting points
maxpt <- bestresu$par
maxpt <- fromFONKtoanyspace(maxpt, colnames(locchull$vertices))
maxparval <- maxpt[CIvar]
maxlogL <- tofKpredict.nohull(maxpt,fixedlist = NULL)
}
#
hull_a <- locchull$a
hull_b <- locchull$b ## within-hullness occurs for ui %*% x -ci => 0 ie a %*% x -b <= 0
xnames <- colnames(locchull$vertices)
parPos <- which(xnames==CIvar)
seuillogL <- maxlogL+CIdlr
eval_g0 <- function(x,a,b) {
names(x) <- xnames
logL <- tofKpredict.nohull(x,fixedlist=NULL)
c(seuillogL-logL, a %*% x-b) ## using max() appears fatal...
} ## must be <0
if ( "NLOPT_LD_MMA_for_CI" %in% INFO$optimizers) {
SignedParPosIndic <- as.numeric(seq(length(xnames))==parPos) ## for eval_grad_f0
eval_jac_g0 <- function( x, a, b ) { ## (jacobian:) gradients of constraints
# one row for each constraint, one col for each dim of x
# first row is -dlogL/dx, other rows are trivial.
locfn <- function(x) {
names(x) <- xnames
return(- tofKpredict.nohull(x,fixedlist=NULL))
}
firstrow <- grad(func=locfn, x=x)
rbind(firstrow,a)
}
locOptimizer <- "NLOPT_LD_MMA"
xtol_rel <- 1e-4 ## -4 is default value; e-12 is never attained;
} else{
# NLOPT_LN_COBYLA seems to behave oddly. Even given a good starting point, it may look near the (lower) bound after a few iterations, optimize from there, and terminate at a worse solution than the initial point.
eval_grad_f0 <- eval_jac_g0 <- NULL
locOptimizer <- "NLOPT_LN_COBYLA"
xtol_rel <- 1.0e-12
}
lb <- matrixStats::colMins(locchull$vertices, useNames=TRUE) # apply(locchull$vertices,2,min)
ub <- matrixStats::colMaxs(locchull$vertices, useNames=TRUE) # apply(locchull$vertices,2,max)
CIpoints <- matrix(nrow=0, ncol=INFO$fittedparamnbr)
CIlo <- NA
if (abs(lowval-maxparval)<abs(maxparval*1e-08)) { ## not an error ## lowval==maxparval if relative diff is <1e-16 in absolute value
# CIlo <- NA ## already so
} else {
errmesbeg <- paste(string1, " for ", CIvar, " could not be computed", sep="")
objfn <- function(x,a,b) {x[parPos]} ## pour la borne inf
if ( locOptimizer == "NLOPT_LD_MMA") {
eval_grad_f0 <- function(x,a,b) { SignedParPosIndic }
}
optr <- try(nloptr(x0=maxpt, eval_f=objfn, lb=lb, ub=ub,
eval_g_ineq=eval_g0, a=hull_a, b=hull_b,
eval_grad_f=eval_grad_f0, eval_jac_g_ineq=eval_jac_g0,
opts=list(algorithm=locOptimizer, xtol_rel=xtol_rel, maxeval=-1)))
if(inherits(optr,"try-error")) {
errmsg <- paste(errmesbeg," (maybe out of sampled range)", sep="")
message.redef(errmsg)
} else {
solution <- optr$solution
names(solution) <- xnames
check_tol <- xtol_rel*max(abs(hull_a %*% abs(solution))) ## cannot be more stringent ?
if(any(hull_a %*% solution -hull_b>check_tol)) {
errmsg <- paste(errmesbeg," (optimization result does not satisfy constraints)", sep="")
warning(errmsg)
} else {
CIpoint <- solution
if (length(notinKgspace)>0L){ ## extra composite hull
CIpoint <- tofullKrigingspace(CIpoint) ## will be used by predict.HLfit
}
logl <- purefn(CIpoint,testhull = FALSE)
if( logl<seuillogL-1e-4) {
errmsg <- paste(errmesbeg," (optimization result has suspiciously low likelihood)", sep="")
warning(errmsg)
} else if( logl>seuillogL+1e-4) {
errmsg <- paste(errmesbeg," (maybe out of sampled range: optimization result has suspiciously high likelihood)", sep="")
message.redef(errmsg)
} else {
CIlo <- solution[parPos]
CIpoints <- rbind(CIpoints, CIpoint)
}
}
}
}
CIup <- NA ## default values...
if (abs(hival-maxparval)<abs(maxparval*1e-08)) { ## not an error
} else {
errmesbeg <- paste(string2, " for ", CIvar, " could not be computed", sep="")
if ( locOptimizer == "NLOPT_LD_MMA") {
SignedParPosIndic <- - SignedParPosIndic
}
objfn <- function(x,a,b) {-x[parPos]} ## pour la borne sup
optr <- try(nloptr(x0=maxpt, eval_f=objfn, lb=lb, ub=ub,
eval_g_ineq=eval_g0, a=hull_a, b=hull_b,
eval_grad_f=eval_grad_f0, eval_jac_g_ineq=eval_jac_g0,
opts=list(algorithm=locOptimizer, xtol_rel=xtol_rel, maxeval=-1)))
if(inherits(optr,"try-error")) {
errmsg <- paste(errmesbeg," (maybe out of sampled range)", sep="")
message.redef(errmsg)
} else {
solution <- optr$solution
names(solution) <- xnames
check_tol <- xtol_rel*max(abs(hull_a %*% abs(solution))) ## cannot be more stringent ?
if(any(hull_a %*% solution -hull_b>check_tol)) {
errmsg <- paste(errmesbeg," (optimization result does not satisfy constraints)", sep="")
warning(errmsg)
} else {
CIpoint <- solution
if (length(notinKgspace)>0L){ ## extra composite hull
CIpoint <- tofullKrigingspace(CIpoint) ## will be used by predict.HLfit
}
logl <- purefn(CIpoint,testhull = FALSE)
if( logl<seuillogL-1e-4) {
errmsg <- paste(errmesbeg," (optimization result has suspiciously low likelihood)", sep="")
warning(errmsg)
} else if( logl>seuillogL+1e-4) {
errmsg <- paste(errmesbeg," (maybe out of sampled range: optimization result has suspiciously high likelihood)", sep="")
message.redef(errmsg)
} else {
CIup <- solution[parPos]
CIpoints <- rbind(CIpoints, CIpoint)
}
}
}
}
return(list(CIlo=CIlo, CIup=CIup, CIpoints=CIpoints))
} ## end def bounds1D
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.