Nothing
#*******************************************************************************
#
# Local Approximate Gaussian Process Regression
# Copyright (C) 2013, The University of Chicago
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
#
# Questions? Contact Robert B. Gramacy (rbg@vt.edu)
#
#*******************************************************************************
## auto.rho:
##
## initial rho parameter chosen to balance objective values
## with (sum of) squared C values nearby to the constraint
## boundary
auto.rho <- function(obj, C, equal, init=10, ethresh=1e-1)
{
if(init > length(obj)) init <- length(obj)
if(length(obj) != nrow(C)) stop("length(obj) != nrow(C)")
## completely valid runs
v <- apply(C, 1, function(x) { all(x[equal] <= 0) && all(x[!equal] < ethresh) })
## smallest sum of squared invalids
Civ <- C[!v,,drop=FALSE]
valid <- Civ <= 0
Civ[valid] <- 0
Civs2 <- rowSums(Civ^2)
if(length(Civs2) == 0) return(1/2)
else cm2 <- min(Civs2)
## smallest valid objective restricted to first init (10)
f <- obj[1:init]; fv <- f[v[1:init]]
if(length(fv) > 0) fm <- min(fv)
else fm <- median(f)
## or if no valids use the median of first init (10)
rho <- cm2/(2*abs(fm))
return(rho)
}
## obj.alM:
##
## objective function for maximizing EI or (minimizing) EY under
## the augmented Langrangian
obj.alM <- function(x, alM, by, fgpi, fnorm, Cgpi, Cnorm, lambda, alpha, ybest,
slack, equal, N, fnc, Bscale)
{
## print(x)
if(is.null(nrow(x))) x <- matrix(x, nrow=1)
eyei <- alM(x, fgpi, fnorm, Cgpi, Cnorm, lambda, alpha, ybest,
slack, equal, N, fnc, Bscale)
## print(c(cbind(x), eyei$ei, eyei$ey, by))
if(by == "ei") return (-eyei$ei)
else return(eyei$ey)
}
## optim.alM:
##
## wrapper function for using optim to maximize EI or (minimize) EY
## under the augmented Lagrangian. This function works best with nomax = -1,
## which currently encodes the slack variable version, which is deterministic;
## with nomax = 0 the original Monte Carlo EI is used and, since the MC is
## a little noisy, optim might not work well
optim.alM <- function(xstart, B, alM, by, fgpi, fnorm, Cgpi, Cnorm, lambda, alpha,
ybest, slack, equal, N, fn, Bscale)
{
out <- optim(par=xstart, fn=obj.alM, method="L-BFGS-B", lower=B[,1], upper=B[,2],
alM=alM, by=by, fgpi=fgpi, fnorm=fnorm, Cgpi=Cgpi, Cnorm=Cnorm, lambda=lambda,
alpha=alpha, ybest=ybest, slack=slack, equal=equal, N=N, fnc=fn, Bscale=Bscale)
return(out)
}
## optim.auglag:
##
## Optimization of known or estimated objective under unknown constraints via
## the Augmented Lagriangian framework with GP surrogate modeling
## of the constraint functions
optim.auglag <- function(fn, B, fhat=FALSE, equal=FALSE, ethresh=1e-2, slack=FALSE,
cknown=NULL, start=10, end=100, Xstart=NULL, sep=TRUE, ab=c(3/2,8),
lambda=1, rho=NULL, urate=10, ncandf=function(t) { t }, dg.start=c(0.1,1e-6),
dlim=sqrt(ncol(B))*c(1/100,10), Bscale=1, ey.tol=1e-2, N=1000,
plotprog=FALSE, verb=2, ...)
{
## check start
if(start >= end) stop("must have start < end")
## check sep and determine whether to use GP or GPsep commands
if(sep) { newM <- newGPsep; mleM <- mleGPsep; updateM <- updateGPsep;
alM <- alGPsep; deleteM <- deleteGPsep; nd <- nrow(B) }
else { newM <- newGP; mleM <- mleGP; updateM <- updateGP;
alM <- alGP; deleteM <- deleteGP; nd <- 1 }
formals(newM)$dK <- TRUE;
## get initial designwhich.min
X <- dopt.gp(start, Xcand=lhs(10*start, B))$XX
## X <- lhs(start, B)
X <- rbind(Xstart, X)
start <- nrow(X)
## first run to determine dimensionality of the constraint
out <- fn(X[1,]*Bscale, ...)
nc <- length(out$c)
## check equal argument
if(length(equal) == 1) {
if(!(equal %in% c(0,1)))
stop("equal should be a logical scalar or vector of lenghth(fn(x)$c)")
if(equal) equal <- rep(1, nc)
else equal <- rep(0, nc)
} else {
if(length(equal) != nc)
stop("equal should be a logical scalar or vector of lenghth(fn(x)$c)")
if(! any(equal != 0 | equal != 1))
stop("equal should be a logical scalar or vector of lenghth(fn(x)$c)")
}
equal <- as.logical(equal)
## allocate progress objects, and initialize
prog <- obj <- rep(NA, start)
C <- matrix(NA, nrow=start, ncol=nc)
obj[1] <- out$obj; C[1,] <- out$c
if(all(out$c[!equal] <= 0) && all(abs(out$c[equal]) < ethresh)) prog[1] <- out$obj
else prog[1] <- Inf
## remainder of starting run
for(t in 2:start) { ## now that fn is vectorized we can probably remove for
out <- fn(X[t,]*Bscale, ...)
obj[t] <- out$obj; C[t,] <- out$c
## update best so far
if(all(out$c[!equal] <= 0) && all(abs(out$c[equal]) < ethresh) &&
out$obj < prog[t-1]) prog[t] <- out$obj
else prog[t] <- prog[t-1]
}
## handle initial lambda and rho values
if(length(lambda) == 1) lambda <- rep(lambda, nc)
if(is.null(rho)) rho <- auto.rho(obj, C, equal)
else if(length(rho) != 1 || rho <= 0) stop("rho should be a positive scalar")
## calculate AL for data seen so far
if(!slack) { ## Original AL
Ce <- matrix(rep(as.numeric(equal), start), nrow=start, byrow=TRUE)
Ce[Ce != 0] <- -Inf
Cm <- pmax(C, Ce)
} else { ## AL with slack variables
S <- pmax(- C - rho*matrix(rep(lambda, start), nrow=start, byrow=TRUE), 0)
S[, equal] <- 0
Cm <- C+S
}
al <- obj + Cm %*% lambda + Cm^2 %*% rep(1/(2*rho), nc)
## best auglag seen so far
ybest <- min(al)
since <- 0
## best valid so far
m2 <- prog[start]
## initializing constraint surrogates
Cgpi <- rep(NA, nc)
d <- matrix(NA, nrow=nc, ncol=nd)
Cnorm <- rep(NA, nc)
for(j in 1:nc) {
if(j %in% cknown) { Cnorm[j] <- 1; Cgpi[j] <- -1 }
else {
Cnorm[j] <- max(abs(C[,j]))
Cgpi[j] <- newM(X, C[,j]/Cnorm[j], dg.start[1], dg.start[2])
d[j,] <- mleM(Cgpi[j], param="d", tmin=dlim[1], tmax=dlim[2], ab=ab,
verb=verb-1)$d
}
}
ds <- matrix(rowMeans(d, na.rm=TRUE), nrow=1)
## possibly initialize objective surrogate
if(fhat) {
fnorm <- max(abs(obj))
fgpi <- newM(X, obj/fnorm, dg.start[1], dg.start[2])
df <- mleM(fgpi, param="d", tmin=dlim[1], tmax=dlim[2], ab=ab, verb=verb-1)$d
dfs <- matrix(df, nrow=1)
} else { fgpi <- -1; fnorm <- 1 }
## init for loop
mei <- Inf
new.params <- FALSE
## keeping track
meis <- mei
lambdas <- as.numeric(lambda)
rhos <- rho
## iterating over the black box evaluations
for(t in (start+1):end) {
## lambda and rho update
if(!slack) { ## Original AL
Ce <- matrix(rep(as.numeric(equal), nrow(C)), nrow=nrow(C), byrow=TRUE)
Ce[Ce != 0] <- -Inf
Cm <- pmax(C, Ce)
al <- obj + C %*% lambda + Cm^2 %*% rep(1/(2*rho), nc)
if(!is.finite(m2)) { ## update if no valid has been found
lambda.new <- rep(0, length=ncol(C))
rho.new <- auto.rho(obj, C, equal)
if(rho.new >= rho) rho.new <- 9*rho/10
else if(rho.new < rho/2) rho.new <- rho/2
} else { ## update if any valid has been found
ck <- C[which.min(al),]
lambda.new <- pmax(0, lambda + (1/rho) * ck)
if(any(ck[equal] > 0) || any(abs(ck[!equal]) > ethresh)) rho.new <- rho/2
else rho.new <- rho
}
} else { ## Slack Variable AL
S <- pmax(- C - rho*matrix(rep(lambda, nrow(C)), nrow=nrow(C), byrow=TRUE), 0)
S[, equal] <- 0
Cm <- C+S
al <- obj + Cm %*% lambda + Cm^2 %*% rep(1/(2*rho), nc)
cmk <- Cm[which.min(al),]
if(!is.finite(m2)) { ## update if no valid has been found
## same as in non-slack case above
lambda.new <- rep(0, length=ncol(C))
## lambda.new <- lambda + (1/rho) * cmk
rho.new <- auto.rho(obj, C, equal)
## rho.new <- auto.rho(obj, Cm)
if(rho.new >= rho) rho.new <- 9*rho/10
else if(rho.new < rho/2) rho.new <- rho/2
} else { ## update if any valid has been found
## specific to slack variable case
lambda.new <- lambda + (1/rho) * cmk
## if(mean(cmk^2) > quantile(rowMeans(Cm^2), p=0.05)) rho.new <- rho/2
if(any(cmk[!equal] > 0) || any(abs(cmk[equal]) > ethresh)) rho.new <- rho/2
else rho.new <- rho
}
}
## printing progress
if(any(lambda.new != lambda) || rho.new != rho) {
if(verb > 0) {
cat("updating La:")
if(rho.new != rho) cat(" rho=", rho.new, sep="")
if(any(lambda.new != lambda))
cat(" lambda=(", paste(signif(lambda.new,3), collapse=", "),
")", sep="")
cat("\n")
}
new.params <- TRUE
} else new.params <- FALSE
## confirm update of augmented lagrangian
lambda <- lambda.new; rho <- rho.new
if(!slack) ## original AL
ybest <- min(obj + C %*% lambda + Cm^2 %*% rep(1/(2*rho), nc)) ## orig
else ybest <- min(obj + Cm %*% lambda + Cm^2 %*% rep(1/(2*rho), nc)) ## slack
## keep track of lambda and rho
lambdas <- rbind(lambdas, lambda)
rhos <- cbind(rhos, rho)
## rebuild surrogates periodically under new normalized responses
if(t > (start+1) && (t %% urate == 0)) {
## constraint surrogates
Cnorm <- apply(abs(C), 2, max)
for(j in 1:nc) {
if(j %in% cknown) next;
deleteM(Cgpi[j])
d[j,d[j,] < dlim[1]] <- 10*dlim[1]
d[j,d[j,] > dlim[2]] <- dlim[2]/10
Cgpi[j] <- newM(X, C[,j]/Cnorm[j], d[j,], dg.start[2])
d[j,] <- mleM(Cgpi[j], param="d", tmin=dlim[1], tmax=dlim[2],
ab=ab, verb=verb-1)$d
}
ds <- rbind(ds, rowMeans(d, na.rm=TRUE))
## possible objective surrogate
if(fhat) {
deleteM(fgpi)
fnorm <- max(abs(obj))
df[df < dlim[1]] <- 10*dlim[1]
df[df > dlim[2]] <- dlim[2]/10
fgpi <- newM(X, obj/fnorm, df, dg.start[2])
df <- mleM(fgpi, param="d", tmin=dlim[1], tmax=dlim[2], ab=ab,
verb=verb-1)$d
dfs <- rbind(dfs, df)
} else { df <- NULL }
new.params <- FALSE
}
## random candidate grid
ncand <- ncandf(t)
if(!is.finite(m2) || fhat)
XX <- lhs(ncand, B)
else XX <- rbetter(ncand, B, sum(X[which(obj == m2)[1],]))
## NOTE: might be a version of rbetter for fhat
## calculate composite surrogate, and evaluate EI and/or EY
eyei <- alM(XX, fgpi, fnorm, Cgpi, Cnorm, lambda, 1/(2*rho), ybest,
slack, equal, N, fn, Bscale)
eis <- eyei$ei; by <- "ei"
mei <- max(eis)
nzei <- sum(eis > 0)
if(nzei <= ey.tol*ncand) { eis <- -(eyei$ey); by <- "ey"; mei <- Inf }
else if(nzei <= 0.1*ncand) {
if(!is.finite(m2) || fhat)
XX <- rbind(XX, lhs(10*ncand, B))
else XX <- rbind(XX, rbetter(10*ncand, B, sum(X[which(obj == m2)[1],])))
eyei <- alM(XX[-(1:ncand),], fgpi, fnorm, Cgpi, Cnorm, lambda, 1/(2*rho), ybest,
slack, equal, N, fn, Bscale)
eis <- c(eis, eyei$ei)
nzei <- nzei + sum(eis > 0)
if(nzei <= ey.tol*ncand) stop("not sure")
ncand <- ncand + 10*ncand
mei <- max(eis)
}
meis <- c(meis, mei)
## plot progress
if(!is.logical(plotprog) || plotprog) {
par(mfrow=c(1,3+fhat))
plot(prog, type="l", main="progress")
if(is.logical(plotprog)) {
if(length(eis) < 30) { span <- 0.5 } else { span <- 0.1 }
g <- interp.loess(XX[,1], XX[,2], eis, span=span)
} else g <- plotprog(XX[,1], XX[,2], eis)
image(g, xlim=range(X[,1]), ylim=range(X[,2]), main="EI")
valid <- apply(C, 1, function(x) { all(x <= 0) })
if(is.matrix(valid)) valid <- apply(valid, 1, prod)
points(X[1:start,1:2], col=valid[1:start]+3)
points(X[-(1:start),1:2, drop=FALSE], col=valid[-(1:start)]+3, pch=19)
matplot(ds, type="l", lty=1, main="constraint lengthscale")
if(fhat) matplot(dfs, type="l", lty=1, main="objective lengthscale")
}
## calculate next point
m <- which.max(eis)
xstar <- matrix(XX[m,], ncol=ncol(X))
## shift xstar via optim calculation
if(slack > 1) {
out <- optim.alM(xstar, B, alM, by, fgpi, fnorm, Cgpi, Cnorm, lambda,
1/(2*rho), ybest, TRUE, equal, N, fn, Bscale)
## print(cbind(xstar, out$par))
xstar <- out$par
}
## progress meter
if(verb > 0) {
cat("t=", t, " ", sep="")
cat(by, "=", eis[m]/Bscale, " (", nzei, "/", ncand, ")", sep="")
cat("; xbest=[", paste(signif(X[which(obj == m2)[1],],3), collapse=" "), sep="")
cat("]; ybest (v=", m2, ", al=", ybest, ", since=", since, ")\n", sep="")
}
## new run
out <- fn(xstar*Bscale, ...)
ystar <- out$obj; obj <- c(obj, ystar); C <- rbind(C, out$c)
## update GP fits
X <- rbind(X, xstar)
for(j in 1:nc) if(Cgpi[j] >= 0)
updateM(Cgpi[j], xstar, out$c[j]/Cnorm[j], verb=verb-2)
## check if best valid has changed
since <- since + 1
if(all(out$c[!equal] <= 0) && all(abs(out$c[equal]) < ethresh) && ystar < prog[length(prog)]) {
m2 <- ystar; since <- 0
} ## otherwise m2 unchanged; should be the same as prog[length(prog)]
prog <- c(prog, m2)
## check if best auglag has changed
if(!slack) { ## original AL
alstar <- out$obj + lambda %*% drop(out$c)
ce <- as.numeric(equal); ce[ce != 0] <- -Inf
alstar <- drop(alstar + rep(1/(2*rho),nc) %*% pmax(ce, drop(out$c))^2)
} else {
cs <- pmax(-out$c - rho*lambda, 0)
cs[equal] <- 0
cps <- drop(out$c + cs)
alstar <- out$obj + drop(lambda %*% cps + rep(1/(2*rho),nc) %*% cps^2)
}
if(alstar < ybest) { ybest <- alstar; since <- 0 }
}
## delete GP surrogates
for(j in 1:nc) if(Cgpi[j] > 0) deleteM(Cgpi[j])
if(fhat) deleteM(fgpi)
## return output objects
if(!fhat) df <- NULL
return(list(prog=prog, mei=meis, obj=obj, X=X, C=C, d=d, df=df,
lambda=as.matrix(lambdas), rho=as.numeric(rhos)))
}
## optim.efi:
##
## Optimization of known or estimated objective under unknown constraints via
## the EI times the product of validity of the constraint functions, following
## Schonlau et al 1998
optim.efi <- function(fn, B, fhat=FALSE, cknown=NULL,
start=10, end=100, Xstart=NULL, sep=TRUE, ab=c(3/2,8),
urate=10, ncandf=function(t) { t }, dg.start=c(0.1,1e-6),
dlim=sqrt(ncol(B))*c(1/100,10), Bscale=1, plotprog=FALSE,
verb=2, ...)
{
## check start
if(start >= end) stop("must have start < end")
## check sep and determine whether to use GP or GPsep commands
if(sep) { newM <- newGPsep; mleM <- mleGPsep; updateM <- updateGPsep;
efiM <- efiGPsep; deleteM <- deleteGPsep; nd <- nrow(B) }
else { newM <- newGP; mleM <- mleGP; updateM <- updateGP;
efiM <- efiGP; deleteM <- deleteGP; nd <- 1 }
formals(newM)$dK <- TRUE;
## get initial design
X <- dopt.gp(start, Xcand=lhs(10*start, B))$XX
## X <- lhs(start, B)
X <- rbind(Xstart, X)
start <- nrow(X)
## first run to determine dimensionality of the constraint
out <- fn(X[1,]*Bscale, ...)
nc <- length(out$c)
## allocate progress objects, and initialize
prog <- obj <- rep(NA, start)
C <- matrix(NA, nrow=start, ncol=nc)
obj[1] <- out$obj; C[1,] <- out$c
if(all(out$c <= 0)) prog[1] <- out$obj
else prog[1] <- Inf
## remainder of starting run
for(t in 2:start) {
out <- fn(X[t,]*Bscale, ...)
obj[t] <- out$obj; C[t,] <- out$c
## update best so far
if(all(out$c <= 0) && out$obj < prog[t-1]) prog[t] <- out$obj
else prog[t] <- prog[t-1]
}
## best valid so far
m2 <- prog[start]
## initializing constraint surrogates
Cgpi <- rep(NA, nc)
d <- matrix(NA, nrow=nc, ncol=nd)
Cnorm <- rep(NA, nc)
for(j in 1:nc) {
if(j %in% cknown) { Cnorm[j] <- 1; Cgpi[j] <- -1 }
else {
Cnorm[j] <- max(abs(C[,j]))
Cgpi[j] <- newM(X, C[,j]/Cnorm[j], dg.start[1], dg.start[2])
d[j,] <- mleM(Cgpi[j], param="d", tmin=dlim[1], tmax=dlim[2], ab=ab,
verb=verb-1)$d
}
}
ds <- matrix(rowMeans(d, na.rm=TRUE), nrow=1)
## possibly initialize objective surrogate
if(fhat) {
fnorm <- max(abs(obj))
fgpi <- newM(X, obj/fnorm, dg.start[1], dg.start[2])
df <- mleM(fgpi, param="d", tmin=dlim[1], tmax=dlim[2], ab=ab, verb=verb-1)$d
dfs <- matrix(df, nrow=1)
} else { fgpi <- -1; fnorm <- 1 }
## iterating over the black box evaluations
since <- 0
for(t in (start+1):end) {
## rebuild surrogates periodically under new normalized responses
if((t > (start+1)) && (t %% urate == 0)) {
## constraint surrogates
Cnorm <- apply(abs(C), 2, max)
for(j in 1:nc) {
if(j %in% cknown) next;
deleteM(Cgpi[j])
d[j,d[j,] < dlim[1]] <- 10*dlim[1]
d[j,d[j,] > dlim[2]] <- dlim[2]/10
Cgpi[j] <- newM(X, C[,j]/Cnorm[j], d[j,], dg.start[2])
d[j,] <- mleM(Cgpi[j], param="d", tmin=dlim[1], tmax=dlim[2],
ab=ab, verb=verb-1)$d
}
ds <- rbind(ds, rowMeans(d, na.rm=TRUE))
## possible objective surrogate
if(fhat) {
deleteM(fgpi)
fnorm <- max(abs(obj))
fgpi <- newM(X, obj/fnorm, df, dg.start[2])
df <- mleM(fgpi, param="d", tmin=dlim[1], tmax=dlim[2], ab=ab,
verb=verb-1)$d
dfs <- rbind(dfs, df)
} else { df <- NULL }
new.params <- FALSE
}
## random candidate grid
ncand <- ncandf(t)
if(!is.finite(m2) || fhat)
XX <- lhs(ncand, B)
else XX <- rbetter(ncand, B, sum(X[which(obj == m2)[1],]))
## NOTE: might be a version of rbetter for fhat
## calculate composite surrogate, and evaluate EI and/or EY
efi <- efiM(XX, fgpi, fnorm, Cgpi, Cnorm, m2, fn, Bscale)
leis <- efi$lei + rowSums(efi[,-1,drop=FALSE])
## continue plotting
if(!is.logical(plotprog) || plotprog) {
par(mfrow=c(1,3+fhat))
plot(prog, type="l", main="progress")
if(is.logical(plotprog)) {
if(length(leis) < 30) { span <- 0.5 } else { span <- 0.1 }
g <- interp.loess(XX[,1], XX[,2], leis, span=span)
} else g <- plotprog(XX[,1], XX[,2], leis)
image(g, xlim=range(X[,1]), ylim=range(X[,2]), main="EI")
points(X[,1:2], col=valid+3, pch=19)
matplot(ds, type="l", lty=1, main="constraint lengthscale")
if(fhat) matplot(dfs, type="l", lty=1, main="objective lengthscale")
}
## calculate next point
m <- which.max(leis)
xstar <- matrix(XX[m,], ncol=ncol(X))
if(verb > 0) {
cat("t=", t, " ", sep="")
cat("efi=", exp(leis[m])/Bscale, sep="")
cat("; xbest=[", paste(signif(X[which(obj == m2)[1],],3), collapse=" "), sep="")
cat("]; ybest=", m2, ", since=", since, "\n", sep="")
}
## new run
out <- fn(xstar*Bscale, ...)
ystar <- out$obj; obj <- c(obj, ystar); C <- rbind(C, out$c)
## update GP fits
X <- rbind(X, xstar)
for(j in 1:nc) if(Cgpi[j] >= 0)
updateM(Cgpi[j], xstar, out$c[j]/Cnorm[j], verb=verb-2)
## check if best valid has changed
since <- since + 1
valid <- apply(C, 1, function(x) { all(x <= 0) })
if(all(out$c <= 0) && ystar < prog[length(prog)]) {
m2 <- ystar; since <- 0
} ## otherwise m2 unchanged; should be the same as prog[length(prog)]
prog <- c(prog, m2)
}
## delete GP surrogates
for(j in 1:nc) if(Cgpi[j] > 0) deleteM(Cgpi[j])
if(fhat) deleteM(fgpi)
## return output objects
if(!fhat) df <- NULL
return(list(prog=prog, obj=obj, X=X, C=C, d=d, df=df))
}
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.