Nothing
R/gp.R
In laGP: Local Approximate Gaussian Process Regression
Defines functions
getmGP
fishGP
dmus2GP
mspeGP
efiGP
alGP
alcrayGP
lalcrayGP
lalcrayGP.R
ray.end
dalcGP
lalcoptGP
lalcoptGP.R
alcoptGP.R
alcoptGP
alcGP
ieciGP
predGP
updateGP
mleGP.switch
dllikGP
mleGP
jmleGP
jmleGP.R
llikGP
newparamsGP
copyGP
deleteGPs
deleteGP
deletedkGP
buildKGP
newGP
Documented in
alcGP
alcoptGP
alcrayGP
dalcGP
deleteGP
deleteGPs
fishGP
ieciGP
jmleGP
jmleGP.R
llikGP
mleGP
mspeGP
newGP
predGP
updateGP
#*******************************************************************************
#
# 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)
#
#*******************************************************************************
## newGP:
##
## build an initial GP representation on the C-side
## using the X-Z data and d/g paramterization. Calling
## this function writes over the previous GP representation
newGP <- function(X, Z, d, g, dK=FALSE)
{
n <- nrow(X)
if(is.null(n)) stop("X must be a matrix")
if(length(Z) != n) stop("must have nrow(X) = length(Z)")
out <- .C("newGP_R",
m = as.integer(ncol(X)),
n = as.integer(n),
X = as.double(t(X)),
Z = as.double(Z),
d = as.double(d),
g = as.double(g),
dK = as.integer(dK),
gpi = integer(1),
PACKAGE = "laGP")
## return C-side GP index
return(out$gpi)
}
## buildkGP:
##
## allocates/calculates the C-side derivative info (only) for particular GP
buildKGP <- function(gpi)
{
.C("buildKGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deletedkGP:
##
## deletes the C-side derivative info (only) for particular GP
deletedkGP <- function(gpi)
{
.C("deletedKGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deleteGP:
##
## deletes the C-side of a particular GP
deleteGP <- function(gpi)
{
.C("deleteGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deleteGPs:
##
## deletes all gps on the C side
deleteGPs <- function()
{
.C("deleteGPs_R", PACKAGE="laGP")
invisible(NULL)
}
## copyGP:
##
## allocate a new GP with a copy of the contents of an
## old one
copyGP <- function(gpi)
{
r <- .C("copyGP_R",
gpi = as.integer(gpi),
newgpi = integer(1),
PACKAGE = "laGP")
return(r$newgpi)
}
## newparamsGP:
##
## change the GP lengthscale and nugget parameerization
## (without destroying the object and creating a new one)
newparamsGP <- function(gpi, d=-1, g=-1)
{
if(d <= 0 & g < 0) stop("one of d or g must be new")
r <- .C("newparamsGP_R",
gpi = as.integer(gpi),
d = as.double(d),
g = as.double(g),
PACKAGE = "laGP")
invisible(NULL)
}
## llikGP:
##
## calculate the log likelihood of the GP
llikGP <- function(gpi, dab=c(0,0), gab=c(0,0))
{
r <- .C("llikGP_R",
gpi = as.integer(gpi),
dab = as.double(dab),
gab = as.double(gab),
llik = double(1),
PACKAGE = "laGP")
return(r$llik)
}
## jmleGP.R:
##
## joint MLE for lengthscale (d) and nugget (g) parameters;
## updates the internal GP parameterization (since mleGP does);
## R-only version
jmleGP.R <- function(gpi, N=100, drange=c(sqrt(.Machine$double.eps), 10),
grange=c(sqrt(.Machine$double.eps), 1), dab=c(0,0), gab=c(0,0), verb=0)
{
## sanity check N
if(length(N) != 1 && N > 0)
stop("N should be a positive scalar integer")
dmle <- gmle <- rep(NA, N)
dits <- gits <- rep(NA, N)
## sanity check tmin and tmax
if(length(drange) != 2) stop("drange should be a 2-vector for c(min,max)")
if(length(grange) != 2) stop("grange should be a 2-vector for c(min,max)")
## loop over outer interations
for(i in 1:N) {
d <- mleGP(gpi, param="d", tmin=drange[1], tmax=drange[2],
ab=dab, verb=verb)
dmle[i] <- d$d; dits[i] <- d$its
g <- mleGP(gpi, param="g", tmin=grange[1], tmax=grange[2],
ab=gab, verb=verb)
gmle[i] <- g$g; gits[i] <- g$its
if(gits[i] <= 1 && dits[i] <= 1) break;
}
## check if not convergedf
if(i == N) warning("max outer its (N=", N, ") reached", sep="")
else {
dmle <- dmle[1:i]; dits <- dits[1:i]
gmle <- gmle[1:i]; gits <- gits[1:i]
}
## total iteration count
totits <- sum(c(dits, gits), na.rm=TRUE)
## assemble return objects
return(list(mle=data.frame(d=dmle[i], g=gmle[i], tot.its=totits),
prog=data.frame(dmle=dmle, dits=dits, gmle=gmle, gits=gits)))
}
## jmleGP
##
## interface to C-version for jmleGP;
## right now doesn't take an N argument -- the C-side hard-codes
## N=100
jmleGP <- function(gpi, drange=c(sqrt(.Machine$double.eps), 10),
grange=c(sqrt(.Machine$double.eps), 1), dab=c(0,0), gab=c(0,0), verb=0)
{
## sanity check tmin and tmax
if(length(drange) != 2) stop("drange should be a 2-vector for c(d,g)")
if(length(grange) != 2) stop("grange should be a 2-vector for c(d,g)")
## sanity check ab
if(length(dab) != 2 || any(dab < 0)) stop("dab should be a positive 2-vector")
if(length(gab) != 2 || any(gab < 0)) stop("gab should be a positive 2-vector")
## call the C-side function
r <- .C("jmleGP_R",
gpi = as.integer(gpi),
verb = as.integer(verb),
drange = as.double(drange),
grange = as.double(grange),
dab = as.double(dab),
gab = as.double(gab),
d = double(1),
g = double(1),
dits = integer(1),
gits = integer(1),
PACKAGE = "laGP")
return(data.frame(d=r$d, g=r$g, tot.its=r$dits+r$gits,
dits=r$dits, gits=r$gits))
}
## mleGP:
##
## updates the GP to use its MLE lengthscale
## parameterization using the current data
mleGP <- function(gpi, param=c("d", "g"),
tmin=sqrt(.Machine$double.eps),
tmax=-1, ab=c(0,0), verb=0)
{
param <- match.arg(param)
if(param == "d") param <- 1
else param <- 2
## sanity check
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector")
r <- .C("mleGP_R",
gpi = as.integer(gpi),
param = as.integer(param),
verb = as.integer(verb),
tmin = as.double(tmin),
tmax = as.double(tmax),
ab = as.double(ab),
theta = double(1),
its = integer(1),
PACKAGE = "laGP")
if(param == 1) return(list(d=r$theta, its=r$its))
else return(list(g=r$theta, its=r$its))
}
## dllikGP:
##
## calculate the first and second derivative of the
## log likelihood of the GP with respect to d, the
## lengthscale parameter
dllikGP <- function(gpi, ab=c(0,0), param=c("d", "g"))
{
param <- match.arg(param)
if(param == "d") {
r <- .C("dllikGP_R",
gpi = as.integer(gpi),
ab = as.double(ab),
d = double(1),
d2 = double(1),
PACKAGE = "laGP")
} else {
r <- .C("dllikGP_nug_R",
gpi = as.integer(gpi),
ab = as.double(ab),
d = double(1),
d2 = double(1),
PACKAGE = "laGP")
}
return(data.frame(d=r$d, d2=r$d2))
}
## mleGP.switch:
##
## switch function for mle calculaitons by laGP.R
mleGP.switch <- function(gpi, method, d, g, verb)
{
## do nothing if no MLE required
if(!(d$mle || g$mle)) return(NULL)
## calculate derivatives
if(d$mle && method != "mspe" && method != "fish") buildKGP(gpi)
## switch
if(d$mle && g$mle) {
## joint lengthscale and nugget
return(jmleGP(gpi, drange=c(d$min,d$max), grange=c(g$min, g$max),
dab=d$ab, gab=g$ab))
} else { ## maybe one or the other
if(d$mle) { ## lengthscale only
dmle <- mleGP(gpi, param="d", d$min, d$max, d$ab, verb=verb)
return(data.frame(d=dmle$d, dits=dmle$its))
}
if(g$mle) { ## nugget only
gmle <- mleGP(gpi, param="g", g$min, g$max, g$ab, verb=verb)
return(data.frame(g=gmle$g, gits=gmle$its))
}
}
}
## updateGP:
##
## add X-Z pairs to the C-side GP represnetation
## using only O(n^2) for each pair
updateGP <- function(gpi, X, Z, verb=0)
{
if(length(Z) != nrow(X))
stop("bad dims")
out <- .C("updateGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(X)),
n = as.integer(nrow(X)),
X = as.double(t(X)),
Z = as.double(Z),
verb = as.integer(verb),
PACKAGE = "laGP")
invisible(NULL)
}
## predGP
##
## obtain the parameters to a multivariate-t
## distribution describing the predictive surface
## of the fitted GP model
predGP <- function(gpi, XX, lite=FALSE, nonug=FALSE)
{
nn <- nrow(XX)
if(is.null(nn) || nn == 0) stop("XX bad dims")
if(lite) {
out <- .C("predGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
lite = as.integer(TRUE),
nonug = as.integer(nonug),
mean = double(nn),
s2 = double(nn),
df = double(1),
llik = double(1),
PACKAGE = "laGP")
## coerce matrix output
return(list(mean=out$mean, s2=out$s2, df=out$df, llik=out$llik))
} else {
out <- .C("predGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
lite = as.integer(FALSE),
nonug = as.integer(nonug),
mean = double(nn),
Sigma = double(nn*nn),
df = double(1),
llik = double(1),
PACKAGE = "laGP")
## coerce matrix output
Sigma <- matrix(out$Sigma, ncol=nn)
## return parameterization
return(list(mean=out$mean, Sigma=Sigma, df=out$df, llik=out$llik))
}
}
## ieciGP:
##
## wrapper used to calculate the IECIs in C using
## the pre-stored isotropic GP representation.
ieciGP <- function(gpi, Xcand, fmin, Xref=Xcand, w=NULL, nonug=FALSE, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
nref <- nrow(Xref)
if(is.null(w)) wb <- 0
else {
wb <- 1
if(length(w) != nref || any(w < 0))
stop("w must be a non-negative vector of length nrow(Xref)")
}
out <- .C("ieciGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
fmin = as.double(fmin),
Xref = as.double(t(Xref)),
nref = as.integer(nref),
w = as.double(w),
wb = as.integer(wb),
nonug = as.integer(nonug),
verb = as.integer(verb),
iecis = double(ncand),
PACKAGE = "laGP")
return(out$iecis)
}
## alcGP:
##
## wrapper used to calculate the ALCs in C using
## the pre-stored GP representation. Note that this only
## calculates the s2' component of ds2 = s2 - s2'
alcGP <- function(gpi, Xcand, Xref=Xcand, parallel=c("none", "omp", "gpu"),
verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
parallel <- match.arg(parallel)
if(parallel == "omp") {
if(!is.loaded("alcGP_omp_R")) stop("OMP not supported in this build; please re-compile")
out <- .C("alcGP_omp_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
} else if(parallel == "gpu") {
out <- .C("alcGP_gpu_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
} else {
out <- .C("alcGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
}
return(out$alcs)
}
## alcoptGP:
##
## interface to C version of alcoptGP.R which continuously optimizes
## ALC based on derivatives, using the starting locations and bounding
## boxes and (stored) gp provided; ... has arguments to optim including
## trace/verb level
alcoptGP <- function(gpi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGP(gpi)
if(ncol(Xref) != m) stop("gpi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpi stored X and start have mismatched cols")
## check lower and upper arguments
if(length(lower) == 1) lower <- rep(lower, m)
else if(length(lower) != m) stop("lower should be a vector of length ncol(Xref)")
if(length(upper) == 1) upper <- rep(upper, m)
else if(length(upper) != m) stop("upper should be a vector of length ncol(Xref)")
if(any(lower >= upper)) stop("some lower >= upper")
out <- .C("alcoptGP_R",
gpi = as.integer(gpi),
maxit = as.integer(maxit),
verb = as.integer(verb),
start = as.double(start),
lower = as.double(lower),
upper = as.double(upper),
m = as.integer(m),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
par = double(m),
value = double(1),
counts = integer(2),
msg = paste(rep(0,60), collapse=""),
convergence = integer(1),
PACKAGE = "laGP")
## for now return the whole optim output
return(list(par=out$par, value=out$value, its=out$counts, msg=out$msg, convergence=out$convergence))
}
## alcoptGP.R:
##
## continuously optimizes ALC based on derivatives, using the
## starting locations and (stored) gp provided; ... has arguments
## to optim including trace/verb level
alcoptGP.R <- function(gpi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGP(gpi)
if(ncol(Xref) != m) stop("gpi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpi stored X and start have mismatched cols")
## objective (and derivative saved)
deriv <- NULL
f <- function(x, gpi, Xref) {
out <- dalcGP(gpi, matrix(x, nrow=1), Xref, verb=0)
deriv <<- list(x=x, df=-out$dalcs/out$alcs)
return(- log(out$alcs))
}
## derivative read from global variable
df <- function(x, gpi, Xref) {
if(any(x != deriv$x)) stop("xs don't match for successive f and df calls")
return(deriv$df)
}
## set up control
control <- list(maxit=maxit, trace=verb, pgtol=1e-1)
## call optim with derivative and global variable
opt <- optim(start, f, df, gpi=gpi, Xref=Xref, lower=lower, upper=upper,
method="L-BFGS-B", control=control)
## version without derivatives
## opt <- optim(start, f, gpi=gpi, Xref=Xref, lower=lower, upper=upper,
## method="L-BFGS-B", control=control)
## keep track of derivative progress
## f(opt$par, gpi, Xref)
## grads <<- rbind(grads, df(opt$par, gpi, Xref))
## for now return the whole optim output
return(opt)
}
## lalcoptGP.R:
##
## optimizes ALC continuously from an initial random nearest (of start)
## neighbor(s) to Xref in Xcand. The candidate in Xcand which is closest
## to the solution is returned. This works differently than
## lalcoptGP, since the starts are random from 1:offset
lalcoptGP.R <- function(gpi, Xref, Xcand, rect=NULL, offset=1, numstart=1, verb=0)
{
## sanity checks
m <- ncol(Xref)
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(length(offset) != 1 || offset < 1 || offset > nrow(Xcand))
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numstart) != 1 || numstart < 1)
stop("numstart should be an integer scalar >= 1")
## adjust numstart
if(numstart > nrow(Xcand)) numstart <- nrow(Xcand)
## calculate bounding rectangle from candidates
if(is.null(rect)) rect <- apply(Xcand, 2, range)
else if(nrow(rect) != 2 || ncol(rect) != ncol(Xref))
stop("bad rect dimensions, must be 2 x ncol(Xref)")
## get starting and ending point of ray
Xstart <- Xcand[offset:(offset + numstart - 1),,drop=FALSE]
## multi-start scheme for searching via derivatives
best.obj <- -Inf; best.w <- NA
for(i in 1:nrow(Xstart)) {
opt <- alcoptGP(gpi, Xref, Xstart[i,], rect[1,], rect[2,], verb=verb)
## opt <- alcoptGP.R(gpi, Xref, Xstart[i,], rect[1,], rect[2,], verb=verb)
## calculate the index of the closest Xcand to opt$par and evaluate
## the ALC criteria there
w <- which.min(distance(matrix(opt$par, nrow=1), Xcand)[1,])
obj <- alcGP(gpi, Xcand[w,,drop=FALSE], Xref)
## determine if that location has the best ALC so far
if(obj > best.obj) { best.obj <- obj; best.w <- w }
}
return(best.w)
}
## lalcoptGP:
##
## wrapper to a C-side function used to optimize ALC continuously
## from an initial neighbor(s) to Xref in Xcand.
## The candidate in Xcand which is closest to the solution is returned.
## This works differently than lalcoptGP.R, since the starts are
## determined by a deterministic round robin similar to lalcrayGP
lalcoptGP <- function(gpi, Xref, Xcand, rect=NULL, offset=1, numstart=1, maxit=100,
verb=0)
{
## sanity checks
m <- ncol(Xref)
ncand <- nrow(Xcand)
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(length(offset) != 1 || offset < 1 || offset > ncand)
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numstart) != 1 || numstart < 1)
stop("numstart should be an integer scalar >= 1")
## calculate bounding rectangle from candidates
if(is.null(rect)) rect <- apply(Xcand, 2, range)
else if(nrow(rect) != 2 || ncol(rect) != ncol(Xref))
stop("bad rect dimensions, must be 2 x ncol(Xref)")
out <- .C("lalcoptGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
offset = as.integer(offset-1),
numstart = as.integer(numstart),
rect = as.double(t(rect)),
maxit = as.integer(maxit),
verb = as.integer(verb),
w = integer(1),
PACKAGE = "laGP")
return(out$w+1)
}
## dalcGP:
##
## wrapper used to calculate the derivative of ALCs in C using
## the pre-stored GP representation. Note that this only
## calculates the s2' component of ds2 = s2 - s2'
dalcGP <- function(gpi, Xcand, Xref=Xcand, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
out <- .C("dalcGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
dalcs = double(ncand*m),
PACKAGE = "laGP")
return(list(alcs=out$alcs, dalcs=matrix(out$dalcs, ncol=m, byrow=TRUE)))
}
## ray.bounds:
##
## get the end of the ray emanating from Xstart that goes
## away (in the direction) from Xref which is 10 times the
## distance but not beyond the bounding rectangle
ray.end <- function(numrays, Xref, Xstart, rect)
{
for(r in 1:numrays) {
Xend[r,] <- as.numeric(10*(Xstart[r,] - Xref) + Xstart[r,])
while(any(Xend[r,] < rect[1,]) | any(Xend[r,] > rect[2,])) {
w <- which(Xend[r,] < rect[1,])
if(length(w) > 0) { col <- 1
} else { w <- which(Xend[r,] > rect[2,]); col <- 2 }
w <- w[1]
sc <- (rect[col,w] - Xstart[r,w])/(Xend[r,w] - Xstart[r,w])
Xend[r,] <- (Xend[r,] - Xstart[r,])*sc + Xstart[r,]
}
}
return(Xend)
}
## lalcrayGP.R:
##
## calculates a ray emanating from a random nearest (of start)
## neighbor(s) to Xref in Xcand. The ending point of the ray
## is 10 times the (opposite) distance from Xstart to Xref,
## then alcrayGP (either C or R version) is called to optimize
## over the ray. The candidate in Xcand which is closest
## to the solution is returned. This works differently than
## lalcrayGP, since the starts of the rays are random from
## 1:offset
lalcrayGP.R <- function(gpi, Xref, Xcand, rect, offset=1, numrays=ncol(Xref),
verb=0)
{
## sanity checks
m <- ncol(Xref)
if(nrow(Xref) != 1) stop("alcray only applies for one Xref")
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(ncol(rect) != m) stop("ncol(rect) != ncol(Xref)")
if(length(offset) != 1 || offset < 1 || offset > nrow(Xcand))
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numrays) != 1 || numrays < 1)
stop("numrays should be an integer scalar >= 1")
## adjust numrays
if(numrays > nrow(Xcand)) numrays <- nrow(Xcand)
## get starting and ending point of ray
Xstart <- Xcand[sample(1:offset, numrays),,drop=FALSE]
Xend <- ray.end(numrays, Xref, Xstart, rect)
## solve for the best convex combination of Xstart and Xend
so <- alcrayGP(gpi, Xref, Xstart, Xend, verb)
Xstar <- matrix((1-so$s)*Xstart[so$r,] + so$s*Xend[so$r,], nrow=1)
## return the index of the closest Xcand to Xstar
w <- which.min(distance(Xstar, Xcand)[1,])
return(w)
}
## lalcrayGP:
##
## wrapper to a C-side function used to calculate a ray emanating
## from a random nearest (of start) neighbor(s) to Xref in Xcand.
## The ending point of the ray is 10 times the (opposite) distance
## from Xstart to Xref, then alcrayGP (on the C-side) is called to
## optimize over the ray. The candidate in Xcand which is closest
## to the solution is returned
lalcrayGP <- function(gpi, Xref, Xcand, rect, offset=1, numrays=ncol(Xref), verb=0)
{
## sanity checks
m <- ncol(Xref)
ncand <- nrow(Xcand)
if(nrow(Xref) != 1) stop("alcray only applies for one Xref")
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(ncol(rect) != m) stop("ncol(rect) != ncol(Xref)")
if(length(offset) != 1 || offset < 1 || offset > ncand)
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numrays) != 1 || numrays < 1)
stop("numrays should be an integer scalar >= 1")
out <- .C("lalcrayGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
offset = as.integer(offset-1),
numrays = as.integer(numrays),
rect = as.double(t(rect)),
verb = as.integer(verb),
w = integer(1),
PACKAGE = "laGP")
return(out$w+1)
}
## alcrayGP:
##
## wrapper used to optimize AIC via a ray search using
## the pre-stored GP representation. Return the convex
## combination s in (0,1) between Xstart and Xend
alcrayGP <- function(gpi, Xref, Xstart, Xend, verb=0)
{
## coerse to matrices
if(is.null(ncol(Xref))) Xref <- matrix(Xref, nrow=1)
if(is.null(ncol(Xstart))) Xstart <- matrix(Xstart, nrow=1)
if(is.null(ncol(Xend))) Xend <- matrix(Xend, nrow=1)
## check dimensions of matrices
m <- ncol(Xstart)
if(ncol(Xref) != m) stop("Xstart and Xref have mismatched cols")
if(ncol(Xend) != m) stop("Xend and Xref have mismatched cols")
if(nrow(Xref) != 1) stop("only one reference location allowed for ray search")
numrays <- nrow(Xstart)
if(nrow(Xend) != numrays) stop("must have same number of starting and ending locations")
## call the C routine
out <- .C("alcrayGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xref = as.double(t(Xref)),
numrays = as.integer(numrays),
Xstart = as.double(t(Xstart)),
Xend = as.double(t(Xend)),
verb = as.integer(verb),
s = double(1),
r = integer(1),
PACKAGE = "laGP")
## return the convex combination
return(list(r=out$r, s=out$s))
}
## alGP:
##
## calculate the E(Y) and EI(Y) for an augmented Lagrangian
## composite objective function with linear objective (in X)
## and constraint GP (gpi) predictive surfaces
alGP <- function(XX, fgpi, fnorm, Cgpis, Cnorms, lambda, alpha, ymin,
slack=FALSE, equal=rep(FALSE,length(Cgpis)), N=100, fn=NULL, Bscale=1)
{
## dimensions
m <- ncol(XX)
nn <- nrow(XX)
nCgps <- length(Cgpis)
## checking lengths for number of gps
if(length(Cnorms) != nCgps) stop("length(Cgpis) != length(Cnorms)")
if(length(lambda) != nCgps) stop("length(Cgpis) != length(lambda)")
if(length(alpha) != 1) stop("length(alpha) != 1")
## checking scalars
if(length(equal) != length(Cgpis)) stop("equal should be a vector of length(Cgpis)")
if(length(N) != 1 || N <= 0) stop("N should be a positive integer scalar")
if(length(ymin) != 1) stop("ymin should be a scalar")
if(length(fnorm) != 1) stop("fnorm should be a scalar")
## run fn to get cheap objectives and constraints
if(fgpi < 0 || any(Cgpis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpi or Cgpis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpi < 0) {
if(is.null(out$obj)) stop("fgpi < 0 but out$obj from fn() is NULL")
obj <- out$obj
} else obj <- NULL
if(any(Cgpis < 0)) {
C <- out$c
if(ncol(C) != sum(Cgpis < 0)) stop("ncol(C) != sum(Cgpis < 0)")
} else C <- NULL
} else { obj <- C <- NULL }
## call the C-side
out <- .C("alGP_R",
m = as.integer(m),
XX = as.double(t(XX)),
nn = as.integer(nn),
fgpi = as.integer(fgpi),
fnorm = as.double(fnorm),
nCgps = as.integer(nCgps),
Cgpis = as.integer(Cgpis),
Cnorms = as.double(Cnorms),
lambda = as.double(lambda),
alpha = as.double(alpha),
ymin = as.double(ymin),
slack = as.integer(slack),
equal = as.double(equal),
N = as.integer(N),
eys = double(nn),
eis = double(nn),
PACKAGE = "laGP")
## done
return(data.frame(ey=out$eys, ei=out$eis))
}
## efiGP:
##
## calculate EI(f) and p(Y(c) <= 0) for known linear or esitmated
## objective f and vectorized constraints C via isotropic GP (gpsi)
## predictive surfaces; returns log probabilities (lplex) and
## EIs on the original scale
efiGP <- function(XX, fgpi, fnorm, Cgpis, Cnorms, fmin, fn=NULL, Bscale=1)
{
## doms
m <- ncol(XX)
nn <- nrow(XX)
nCgps <- length(Cgpis)
## checking lengths for number of constraint gps
if(length(Cnorms) != nCgps) stop("length(Cgpis) != length(Cnorms)")
## checking scalars
if(length(fmin) != 1) stop("ymin should be a scalar")
if(length(fnorm) != 1) stop("fnorm should be a scalar")
## run fn to get cheap objectives and constraints
if(fgpi < 0 || any(Cgpis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpi or Cgpis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpi < 0) {
if(is.null(out$obj)) stop("fgpi < 0 but out$obj from fn() is NULL")
obj <- out$obj
} else obj <- NULL
if(any(Cgpis < 0)) {
C <- out$c
if(ncol(C) != sum(Cgpis < 0)) stop("ncol(C) != sum(Cgpis < 0)")
} else C <- NULL
}
## calculate expected improvement part
if(fgpi <= 0) {
obj <- rowSums(XX) * fnorm
if(!is.finite(fmin)) fmin <- quantile(obj, p=0.9)
I <- fmin - obj
ei <- pmax(I, 0)
} else {
p <- predGP(fgpi, XX=XX, lite=TRUE)
pm <- p$mean * fnorm
ps <- sqrt(p$s2) * fnorm
if(!is.finite(fmin)) fmin <- quantile(pm, p=0.9)
u <- (fmin - pm)/ps
ei <- ps*dnorm(u) + (fmin-pm)*pnorm(u)
}
## calculate constraint part
lplez <- matrix(NA, nrow=nrow(XX), nCgps)
for(j in 1:nCgps) {
pc <- predGP(Cgpis[j], XX=XX, lite=TRUE)
lplez[,j] <- pnorm(0, pc$mean, sqrt(pc$s2), log.p=TRUE)
}
## done
return(data.frame(ei=ei, lplez=lplez))
}
## mspeGP:
##
## wrapper used to calculate the MSPEs in C using
## the pre-stored GP representation.
mspeGP <- function(gpi, Xcand, Xref=Xcand, fi=TRUE, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
out <- .C("mspeGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
fi = as.integer(fi),
verb = as.integer(verb),
mspes = double(ncand),
PACKAGE = "laGP")
return(out$mspes)
}
## dmus2GP:
##
## obtain the derivative of the predictive scale
## of the fitted GP model
dmus2GP <- function(gpi, XX)
{
nn <- nrow(XX)
out <- .C("dmus2GP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
mu = double(nn),
dmu = double(nn),
d2mu = double(nn),
s2 = double(nn),
ds2 = double(nn),
d2s2 = double(nn),
PACKAGE = "laGP")
return(data.frame(mu=out$mu, dmu=out$dmu, d2mu=out$d2mu,
s2=out$s2, ds2=out$ds2, d2s2=out$d2s2))
}
## fishGP:
##
## obtain the expected (approx) Fisher information for
## the fitted GP model; returns the absolute value (i.e.,
## determinant)
fishGP <- function(gpi, Xcand)
{
nn <- nrow(Xcand)
out <- .C("efiGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(Xcand)),
nn = as.integer(nn),
Xcand = as.double(t(Xcand)),
efi = double(nn),
PACKAGE = "laGP")
## remove silly values
return(out$efi)
}
## getmGP:
##
## access the input dimension of a GP
getmGP <- function(gpi)
{
.C("getmGP_R", gpi = as.integer(gpi), m = integer(1), PACKAGE="laGP")$m
}
Try the laGP package in your browser
Any scripts or data that you put into this service are public.
laGP documentation built on March 31, 2023, 9:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions getmGP fishGP dmus2GP mspeGP efiGP alGP alcrayGP lalcrayGP lalcrayGP.R ray.end dalcGP lalcoptGP lalcoptGP.R alcoptGP.R alcoptGP alcGP ieciGP predGP updateGP mleGP.switch dllikGP mleGP jmleGP jmleGP.R llikGP newparamsGP copyGP deleteGPs deleteGP deletedkGP buildKGP newGP
Documented in alcGP alcoptGP alcrayGP dalcGP deleteGP deleteGPs fishGP ieciGP jmleGP jmleGP.R llikGP mleGP mspeGP newGP predGP updateGP
#*******************************************************************************
#
# 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)
#
#*******************************************************************************
## newGP:
##
## build an initial GP representation on the C-side
## using the X-Z data and d/g paramterization. Calling
## this function writes over the previous GP representation
newGP <- function(X, Z, d, g, dK=FALSE)
{
n <- nrow(X)
if(is.null(n)) stop("X must be a matrix")
if(length(Z) != n) stop("must have nrow(X) = length(Z)")
out <- .C("newGP_R",
m = as.integer(ncol(X)),
n = as.integer(n),
X = as.double(t(X)),
Z = as.double(Z),
d = as.double(d),
g = as.double(g),
dK = as.integer(dK),
gpi = integer(1),
PACKAGE = "laGP")
## return C-side GP index
return(out$gpi)
}
## buildkGP:
##
## allocates/calculates the C-side derivative info (only) for particular GP
buildKGP <- function(gpi)
{
.C("buildKGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deletedkGP:
##
## deletes the C-side derivative info (only) for particular GP
deletedkGP <- function(gpi)
{
.C("deletedKGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deleteGP:
##
## deletes the C-side of a particular GP
deleteGP <- function(gpi)
{
.C("deleteGP_R",
gpi = as.integer(gpi),
PACKAGE = "laGP")
invisible(NULL)
}
## deleteGPs:
##
## deletes all gps on the C side
deleteGPs <- function()
{
.C("deleteGPs_R", PACKAGE="laGP")
invisible(NULL)
}
## copyGP:
##
## allocate a new GP with a copy of the contents of an
## old one
copyGP <- function(gpi)
{
r <- .C("copyGP_R",
gpi = as.integer(gpi),
newgpi = integer(1),
PACKAGE = "laGP")
return(r$newgpi)
}
## newparamsGP:
##
## change the GP lengthscale and nugget parameerization
## (without destroying the object and creating a new one)
newparamsGP <- function(gpi, d=-1, g=-1)
{
if(d <= 0 & g < 0) stop("one of d or g must be new")
r <- .C("newparamsGP_R",
gpi = as.integer(gpi),
d = as.double(d),
g = as.double(g),
PACKAGE = "laGP")
invisible(NULL)
}
## llikGP:
##
## calculate the log likelihood of the GP
llikGP <- function(gpi, dab=c(0,0), gab=c(0,0))
{
r <- .C("llikGP_R",
gpi = as.integer(gpi),
dab = as.double(dab),
gab = as.double(gab),
llik = double(1),
PACKAGE = "laGP")
return(r$llik)
}
## jmleGP.R:
##
## joint MLE for lengthscale (d) and nugget (g) parameters;
## updates the internal GP parameterization (since mleGP does);
## R-only version
jmleGP.R <- function(gpi, N=100, drange=c(sqrt(.Machine$double.eps), 10),
grange=c(sqrt(.Machine$double.eps), 1), dab=c(0,0), gab=c(0,0), verb=0)
{
## sanity check N
if(length(N) != 1 && N > 0)
stop("N should be a positive scalar integer")
dmle <- gmle <- rep(NA, N)
dits <- gits <- rep(NA, N)
## sanity check tmin and tmax
if(length(drange) != 2) stop("drange should be a 2-vector for c(min,max)")
if(length(grange) != 2) stop("grange should be a 2-vector for c(min,max)")
## loop over outer interations
for(i in 1:N) {
d <- mleGP(gpi, param="d", tmin=drange[1], tmax=drange[2],
ab=dab, verb=verb)
dmle[i] <- d$d; dits[i] <- d$its
g <- mleGP(gpi, param="g", tmin=grange[1], tmax=grange[2],
ab=gab, verb=verb)
gmle[i] <- g$g; gits[i] <- g$its
if(gits[i] <= 1 && dits[i] <= 1) break;
}
## check if not convergedf
if(i == N) warning("max outer its (N=", N, ") reached", sep="")
else {
dmle <- dmle[1:i]; dits <- dits[1:i]
gmle <- gmle[1:i]; gits <- gits[1:i]
}
## total iteration count
totits <- sum(c(dits, gits), na.rm=TRUE)
## assemble return objects
return(list(mle=data.frame(d=dmle[i], g=gmle[i], tot.its=totits),
prog=data.frame(dmle=dmle, dits=dits, gmle=gmle, gits=gits)))
}
## jmleGP
##
## interface to C-version for jmleGP;
## right now doesn't take an N argument -- the C-side hard-codes
## N=100
jmleGP <- function(gpi, drange=c(sqrt(.Machine$double.eps), 10),
grange=c(sqrt(.Machine$double.eps), 1), dab=c(0,0), gab=c(0,0), verb=0)
{
## sanity check tmin and tmax
if(length(drange) != 2) stop("drange should be a 2-vector for c(d,g)")
if(length(grange) != 2) stop("grange should be a 2-vector for c(d,g)")
## sanity check ab
if(length(dab) != 2 || any(dab < 0)) stop("dab should be a positive 2-vector")
if(length(gab) != 2 || any(gab < 0)) stop("gab should be a positive 2-vector")
## call the C-side function
r <- .C("jmleGP_R",
gpi = as.integer(gpi),
verb = as.integer(verb),
drange = as.double(drange),
grange = as.double(grange),
dab = as.double(dab),
gab = as.double(gab),
d = double(1),
g = double(1),
dits = integer(1),
gits = integer(1),
PACKAGE = "laGP")
return(data.frame(d=r$d, g=r$g, tot.its=r$dits+r$gits,
dits=r$dits, gits=r$gits))
}
## mleGP:
##
## updates the GP to use its MLE lengthscale
## parameterization using the current data
mleGP <- function(gpi, param=c("d", "g"),
tmin=sqrt(.Machine$double.eps),
tmax=-1, ab=c(0,0), verb=0)
{
param <- match.arg(param)
if(param == "d") param <- 1
else param <- 2
## sanity check
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector")
r <- .C("mleGP_R",
gpi = as.integer(gpi),
param = as.integer(param),
verb = as.integer(verb),
tmin = as.double(tmin),
tmax = as.double(tmax),
ab = as.double(ab),
theta = double(1),
its = integer(1),
PACKAGE = "laGP")
if(param == 1) return(list(d=r$theta, its=r$its))
else return(list(g=r$theta, its=r$its))
}
## dllikGP:
##
## calculate the first and second derivative of the
## log likelihood of the GP with respect to d, the
## lengthscale parameter
dllikGP <- function(gpi, ab=c(0,0), param=c("d", "g"))
{
param <- match.arg(param)
if(param == "d") {
r <- .C("dllikGP_R",
gpi = as.integer(gpi),
ab = as.double(ab),
d = double(1),
d2 = double(1),
PACKAGE = "laGP")
} else {
r <- .C("dllikGP_nug_R",
gpi = as.integer(gpi),
ab = as.double(ab),
d = double(1),
d2 = double(1),
PACKAGE = "laGP")
}
return(data.frame(d=r$d, d2=r$d2))
}
## mleGP.switch:
##
## switch function for mle calculaitons by laGP.R
mleGP.switch <- function(gpi, method, d, g, verb)
{
## do nothing if no MLE required
if(!(d$mle || g$mle)) return(NULL)
## calculate derivatives
if(d$mle && method != "mspe" && method != "fish") buildKGP(gpi)
## switch
if(d$mle && g$mle) {
## joint lengthscale and nugget
return(jmleGP(gpi, drange=c(d$min,d$max), grange=c(g$min, g$max),
dab=d$ab, gab=g$ab))
} else { ## maybe one or the other
if(d$mle) { ## lengthscale only
dmle <- mleGP(gpi, param="d", d$min, d$max, d$ab, verb=verb)
return(data.frame(d=dmle$d, dits=dmle$its))
}
if(g$mle) { ## nugget only
gmle <- mleGP(gpi, param="g", g$min, g$max, g$ab, verb=verb)
return(data.frame(g=gmle$g, gits=gmle$its))
}
}
}
## updateGP:
##
## add X-Z pairs to the C-side GP represnetation
## using only O(n^2) for each pair
updateGP <- function(gpi, X, Z, verb=0)
{
if(length(Z) != nrow(X))
stop("bad dims")
out <- .C("updateGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(X)),
n = as.integer(nrow(X)),
X = as.double(t(X)),
Z = as.double(Z),
verb = as.integer(verb),
PACKAGE = "laGP")
invisible(NULL)
}
## predGP
##
## obtain the parameters to a multivariate-t
## distribution describing the predictive surface
## of the fitted GP model
predGP <- function(gpi, XX, lite=FALSE, nonug=FALSE)
{
nn <- nrow(XX)
if(is.null(nn) || nn == 0) stop("XX bad dims")
if(lite) {
out <- .C("predGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
lite = as.integer(TRUE),
nonug = as.integer(nonug),
mean = double(nn),
s2 = double(nn),
df = double(1),
llik = double(1),
PACKAGE = "laGP")
## coerce matrix output
return(list(mean=out$mean, s2=out$s2, df=out$df, llik=out$llik))
} else {
out <- .C("predGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
lite = as.integer(FALSE),
nonug = as.integer(nonug),
mean = double(nn),
Sigma = double(nn*nn),
df = double(1),
llik = double(1),
PACKAGE = "laGP")
## coerce matrix output
Sigma <- matrix(out$Sigma, ncol=nn)
## return parameterization
return(list(mean=out$mean, Sigma=Sigma, df=out$df, llik=out$llik))
}
}
## ieciGP:
##
## wrapper used to calculate the IECIs in C using
## the pre-stored isotropic GP representation.
ieciGP <- function(gpi, Xcand, fmin, Xref=Xcand, w=NULL, nonug=FALSE, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
nref <- nrow(Xref)
if(is.null(w)) wb <- 0
else {
wb <- 1
if(length(w) != nref || any(w < 0))
stop("w must be a non-negative vector of length nrow(Xref)")
}
out <- .C("ieciGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
fmin = as.double(fmin),
Xref = as.double(t(Xref)),
nref = as.integer(nref),
w = as.double(w),
wb = as.integer(wb),
nonug = as.integer(nonug),
verb = as.integer(verb),
iecis = double(ncand),
PACKAGE = "laGP")
return(out$iecis)
}
## alcGP:
##
## wrapper used to calculate the ALCs in C using
## the pre-stored GP representation. Note that this only
## calculates the s2' component of ds2 = s2 - s2'
alcGP <- function(gpi, Xcand, Xref=Xcand, parallel=c("none", "omp", "gpu"),
verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
parallel <- match.arg(parallel)
if(parallel == "omp") {
if(!is.loaded("alcGP_omp_R")) stop("OMP not supported in this build; please re-compile")
out <- .C("alcGP_omp_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
} else if(parallel == "gpu") {
out <- .C("alcGP_gpu_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
} else {
out <- .C("alcGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
PACKAGE = "laGP")
}
return(out$alcs)
}
## alcoptGP:
##
## interface to C version of alcoptGP.R which continuously optimizes
## ALC based on derivatives, using the starting locations and bounding
## boxes and (stored) gp provided; ... has arguments to optim including
## trace/verb level
alcoptGP <- function(gpi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGP(gpi)
if(ncol(Xref) != m) stop("gpi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpi stored X and start have mismatched cols")
## check lower and upper arguments
if(length(lower) == 1) lower <- rep(lower, m)
else if(length(lower) != m) stop("lower should be a vector of length ncol(Xref)")
if(length(upper) == 1) upper <- rep(upper, m)
else if(length(upper) != m) stop("upper should be a vector of length ncol(Xref)")
if(any(lower >= upper)) stop("some lower >= upper")
out <- .C("alcoptGP_R",
gpi = as.integer(gpi),
maxit = as.integer(maxit),
verb = as.integer(verb),
start = as.double(start),
lower = as.double(lower),
upper = as.double(upper),
m = as.integer(m),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
par = double(m),
value = double(1),
counts = integer(2),
msg = paste(rep(0,60), collapse=""),
convergence = integer(1),
PACKAGE = "laGP")
## for now return the whole optim output
return(list(par=out$par, value=out$value, its=out$counts, msg=out$msg, convergence=out$convergence))
}
## alcoptGP.R:
##
## continuously optimizes ALC based on derivatives, using the
## starting locations and (stored) gp provided; ... has arguments
## to optim including trace/verb level
alcoptGP.R <- function(gpi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGP(gpi)
if(ncol(Xref) != m) stop("gpi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpi stored X and start have mismatched cols")
## objective (and derivative saved)
deriv <- NULL
f <- function(x, gpi, Xref) {
out <- dalcGP(gpi, matrix(x, nrow=1), Xref, verb=0)
deriv <<- list(x=x, df=-out$dalcs/out$alcs)
return(- log(out$alcs))
}
## derivative read from global variable
df <- function(x, gpi, Xref) {
if(any(x != deriv$x)) stop("xs don't match for successive f and df calls")
return(deriv$df)
}
## set up control
control <- list(maxit=maxit, trace=verb, pgtol=1e-1)
## call optim with derivative and global variable
opt <- optim(start, f, df, gpi=gpi, Xref=Xref, lower=lower, upper=upper,
method="L-BFGS-B", control=control)
## version without derivatives
## opt <- optim(start, f, gpi=gpi, Xref=Xref, lower=lower, upper=upper,
## method="L-BFGS-B", control=control)
## keep track of derivative progress
## f(opt$par, gpi, Xref)
## grads <<- rbind(grads, df(opt$par, gpi, Xref))
## for now return the whole optim output
return(opt)
}
## lalcoptGP.R:
##
## optimizes ALC continuously from an initial random nearest (of start)
## neighbor(s) to Xref in Xcand. The candidate in Xcand which is closest
## to the solution is returned. This works differently than
## lalcoptGP, since the starts are random from 1:offset
lalcoptGP.R <- function(gpi, Xref, Xcand, rect=NULL, offset=1, numstart=1, verb=0)
{
## sanity checks
m <- ncol(Xref)
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(length(offset) != 1 || offset < 1 || offset > nrow(Xcand))
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numstart) != 1 || numstart < 1)
stop("numstart should be an integer scalar >= 1")
## adjust numstart
if(numstart > nrow(Xcand)) numstart <- nrow(Xcand)
## calculate bounding rectangle from candidates
if(is.null(rect)) rect <- apply(Xcand, 2, range)
else if(nrow(rect) != 2 || ncol(rect) != ncol(Xref))
stop("bad rect dimensions, must be 2 x ncol(Xref)")
## get starting and ending point of ray
Xstart <- Xcand[offset:(offset + numstart - 1),,drop=FALSE]
## multi-start scheme for searching via derivatives
best.obj <- -Inf; best.w <- NA
for(i in 1:nrow(Xstart)) {
opt <- alcoptGP(gpi, Xref, Xstart[i,], rect[1,], rect[2,], verb=verb)
## opt <- alcoptGP.R(gpi, Xref, Xstart[i,], rect[1,], rect[2,], verb=verb)
## calculate the index of the closest Xcand to opt$par and evaluate
## the ALC criteria there
w <- which.min(distance(matrix(opt$par, nrow=1), Xcand)[1,])
obj <- alcGP(gpi, Xcand[w,,drop=FALSE], Xref)
## determine if that location has the best ALC so far
if(obj > best.obj) { best.obj <- obj; best.w <- w }
}
return(best.w)
}
## lalcoptGP:
##
## wrapper to a C-side function used to optimize ALC continuously
## from an initial neighbor(s) to Xref in Xcand.
## The candidate in Xcand which is closest to the solution is returned.
## This works differently than lalcoptGP.R, since the starts are
## determined by a deterministic round robin similar to lalcrayGP
lalcoptGP <- function(gpi, Xref, Xcand, rect=NULL, offset=1, numstart=1, maxit=100,
verb=0)
{
## sanity checks
m <- ncol(Xref)
ncand <- nrow(Xcand)
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(length(offset) != 1 || offset < 1 || offset > ncand)
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numstart) != 1 || numstart < 1)
stop("numstart should be an integer scalar >= 1")
## calculate bounding rectangle from candidates
if(is.null(rect)) rect <- apply(Xcand, 2, range)
else if(nrow(rect) != 2 || ncol(rect) != ncol(Xref))
stop("bad rect dimensions, must be 2 x ncol(Xref)")
out <- .C("lalcoptGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
offset = as.integer(offset-1),
numstart = as.integer(numstart),
rect = as.double(t(rect)),
maxit = as.integer(maxit),
verb = as.integer(verb),
w = integer(1),
PACKAGE = "laGP")
return(out$w+1)
}
## dalcGP:
##
## wrapper used to calculate the derivative of ALCs in C using
## the pre-stored GP representation. Note that this only
## calculates the s2' component of ds2 = s2 - s2'
dalcGP <- function(gpi, Xcand, Xref=Xcand, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
out <- .C("dalcGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
verb = as.integer(verb),
alcs = double(ncand),
dalcs = double(ncand*m),
PACKAGE = "laGP")
return(list(alcs=out$alcs, dalcs=matrix(out$dalcs, ncol=m, byrow=TRUE)))
}
## ray.bounds:
##
## get the end of the ray emanating from Xstart that goes
## away (in the direction) from Xref which is 10 times the
## distance but not beyond the bounding rectangle
ray.end <- function(numrays, Xref, Xstart, rect)
{
for(r in 1:numrays) {
Xend[r,] <- as.numeric(10*(Xstart[r,] - Xref) + Xstart[r,])
while(any(Xend[r,] < rect[1,]) | any(Xend[r,] > rect[2,])) {
w <- which(Xend[r,] < rect[1,])
if(length(w) > 0) { col <- 1
} else { w <- which(Xend[r,] > rect[2,]); col <- 2 }
w <- w[1]
sc <- (rect[col,w] - Xstart[r,w])/(Xend[r,w] - Xstart[r,w])
Xend[r,] <- (Xend[r,] - Xstart[r,])*sc + Xstart[r,]
}
}
return(Xend)
}
## lalcrayGP.R:
##
## calculates a ray emanating from a random nearest (of start)
## neighbor(s) to Xref in Xcand. The ending point of the ray
## is 10 times the (opposite) distance from Xstart to Xref,
## then alcrayGP (either C or R version) is called to optimize
## over the ray. The candidate in Xcand which is closest
## to the solution is returned. This works differently than
## lalcrayGP, since the starts of the rays are random from
## 1:offset
lalcrayGP.R <- function(gpi, Xref, Xcand, rect, offset=1, numrays=ncol(Xref),
verb=0)
{
## sanity checks
m <- ncol(Xref)
if(nrow(Xref) != 1) stop("alcray only applies for one Xref")
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(ncol(rect) != m) stop("ncol(rect) != ncol(Xref)")
if(length(offset) != 1 || offset < 1 || offset > nrow(Xcand))
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numrays) != 1 || numrays < 1)
stop("numrays should be an integer scalar >= 1")
## adjust numrays
if(numrays > nrow(Xcand)) numrays <- nrow(Xcand)
## get starting and ending point of ray
Xstart <- Xcand[sample(1:offset, numrays),,drop=FALSE]
Xend <- ray.end(numrays, Xref, Xstart, rect)
## solve for the best convex combination of Xstart and Xend
so <- alcrayGP(gpi, Xref, Xstart, Xend, verb)
Xstar <- matrix((1-so$s)*Xstart[so$r,] + so$s*Xend[so$r,], nrow=1)
## return the index of the closest Xcand to Xstar
w <- which.min(distance(Xstar, Xcand)[1,])
return(w)
}
## lalcrayGP:
##
## wrapper to a C-side function used to calculate a ray emanating
## from a random nearest (of start) neighbor(s) to Xref in Xcand.
## The ending point of the ray is 10 times the (opposite) distance
## from Xstart to Xref, then alcrayGP (on the C-side) is called to
## optimize over the ray. The candidate in Xcand which is closest
## to the solution is returned
lalcrayGP <- function(gpi, Xref, Xcand, rect, offset=1, numrays=ncol(Xref), verb=0)
{
## sanity checks
m <- ncol(Xref)
ncand <- nrow(Xcand)
if(nrow(Xref) != 1) stop("alcray only applies for one Xref")
if(m != ncol(Xcand)) stop("ncol(Xref) != ncol(Xcand)")
if(ncol(rect) != m) stop("ncol(rect) != ncol(Xref)")
if(length(offset) != 1 || offset < 1 || offset > ncand)
stop("offset should be a scalar integer >= 1 and <= nrow(Xcand)")
if(length(numrays) != 1 || numrays < 1)
stop("numrays should be an integer scalar >= 1")
out <- .C("lalcrayGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
offset = as.integer(offset-1),
numrays = as.integer(numrays),
rect = as.double(t(rect)),
verb = as.integer(verb),
w = integer(1),
PACKAGE = "laGP")
return(out$w+1)
}
## alcrayGP:
##
## wrapper used to optimize AIC via a ray search using
## the pre-stored GP representation. Return the convex
## combination s in (0,1) between Xstart and Xend
alcrayGP <- function(gpi, Xref, Xstart, Xend, verb=0)
{
## coerse to matrices
if(is.null(ncol(Xref))) Xref <- matrix(Xref, nrow=1)
if(is.null(ncol(Xstart))) Xstart <- matrix(Xstart, nrow=1)
if(is.null(ncol(Xend))) Xend <- matrix(Xend, nrow=1)
## check dimensions of matrices
m <- ncol(Xstart)
if(ncol(Xref) != m) stop("Xstart and Xref have mismatched cols")
if(ncol(Xend) != m) stop("Xend and Xref have mismatched cols")
if(nrow(Xref) != 1) stop("only one reference location allowed for ray search")
numrays <- nrow(Xstart)
if(nrow(Xend) != numrays) stop("must have same number of starting and ending locations")
## call the C routine
out <- .C("alcrayGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xref = as.double(t(Xref)),
numrays = as.integer(numrays),
Xstart = as.double(t(Xstart)),
Xend = as.double(t(Xend)),
verb = as.integer(verb),
s = double(1),
r = integer(1),
PACKAGE = "laGP")
## return the convex combination
return(list(r=out$r, s=out$s))
}
## alGP:
##
## calculate the E(Y) and EI(Y) for an augmented Lagrangian
## composite objective function with linear objective (in X)
## and constraint GP (gpi) predictive surfaces
alGP <- function(XX, fgpi, fnorm, Cgpis, Cnorms, lambda, alpha, ymin,
slack=FALSE, equal=rep(FALSE,length(Cgpis)), N=100, fn=NULL, Bscale=1)
{
## dimensions
m <- ncol(XX)
nn <- nrow(XX)
nCgps <- length(Cgpis)
## checking lengths for number of gps
if(length(Cnorms) != nCgps) stop("length(Cgpis) != length(Cnorms)")
if(length(lambda) != nCgps) stop("length(Cgpis) != length(lambda)")
if(length(alpha) != 1) stop("length(alpha) != 1")
## checking scalars
if(length(equal) != length(Cgpis)) stop("equal should be a vector of length(Cgpis)")
if(length(N) != 1 || N <= 0) stop("N should be a positive integer scalar")
if(length(ymin) != 1) stop("ymin should be a scalar")
if(length(fnorm) != 1) stop("fnorm should be a scalar")
## run fn to get cheap objectives and constraints
if(fgpi < 0 || any(Cgpis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpi or Cgpis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpi < 0) {
if(is.null(out$obj)) stop("fgpi < 0 but out$obj from fn() is NULL")
obj <- out$obj
} else obj <- NULL
if(any(Cgpis < 0)) {
C <- out$c
if(ncol(C) != sum(Cgpis < 0)) stop("ncol(C) != sum(Cgpis < 0)")
} else C <- NULL
} else { obj <- C <- NULL }
## call the C-side
out <- .C("alGP_R",
m = as.integer(m),
XX = as.double(t(XX)),
nn = as.integer(nn),
fgpi = as.integer(fgpi),
fnorm = as.double(fnorm),
nCgps = as.integer(nCgps),
Cgpis = as.integer(Cgpis),
Cnorms = as.double(Cnorms),
lambda = as.double(lambda),
alpha = as.double(alpha),
ymin = as.double(ymin),
slack = as.integer(slack),
equal = as.double(equal),
N = as.integer(N),
eys = double(nn),
eis = double(nn),
PACKAGE = "laGP")
## done
return(data.frame(ey=out$eys, ei=out$eis))
}
## efiGP:
##
## calculate EI(f) and p(Y(c) <= 0) for known linear or esitmated
## objective f and vectorized constraints C via isotropic GP (gpsi)
## predictive surfaces; returns log probabilities (lplex) and
## EIs on the original scale
efiGP <- function(XX, fgpi, fnorm, Cgpis, Cnorms, fmin, fn=NULL, Bscale=1)
{
## doms
m <- ncol(XX)
nn <- nrow(XX)
nCgps <- length(Cgpis)
## checking lengths for number of constraint gps
if(length(Cnorms) != nCgps) stop("length(Cgpis) != length(Cnorms)")
## checking scalars
if(length(fmin) != 1) stop("ymin should be a scalar")
if(length(fnorm) != 1) stop("fnorm should be a scalar")
## run fn to get cheap objectives and constraints
if(fgpi < 0 || any(Cgpis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpi or Cgpis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpi < 0) {
if(is.null(out$obj)) stop("fgpi < 0 but out$obj from fn() is NULL")
obj <- out$obj
} else obj <- NULL
if(any(Cgpis < 0)) {
C <- out$c
if(ncol(C) != sum(Cgpis < 0)) stop("ncol(C) != sum(Cgpis < 0)")
} else C <- NULL
}
## calculate expected improvement part
if(fgpi <= 0) {
obj <- rowSums(XX) * fnorm
if(!is.finite(fmin)) fmin <- quantile(obj, p=0.9)
I <- fmin - obj
ei <- pmax(I, 0)
} else {
p <- predGP(fgpi, XX=XX, lite=TRUE)
pm <- p$mean * fnorm
ps <- sqrt(p$s2) * fnorm
if(!is.finite(fmin)) fmin <- quantile(pm, p=0.9)
u <- (fmin - pm)/ps
ei <- ps*dnorm(u) + (fmin-pm)*pnorm(u)
}
## calculate constraint part
lplez <- matrix(NA, nrow=nrow(XX), nCgps)
for(j in 1:nCgps) {
pc <- predGP(Cgpis[j], XX=XX, lite=TRUE)
lplez[,j] <- pnorm(0, pc$mean, sqrt(pc$s2), log.p=TRUE)
}
## done
return(data.frame(ei=ei, lplez=lplez))
}
## mspeGP:
##
## wrapper used to calculate the MSPEs in C using
## the pre-stored GP representation.
mspeGP <- function(gpi, Xcand, Xref=Xcand, fi=TRUE, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)
out <- .C("mspeGP_R",
gpi = as.integer(gpi),
m = as.integer(m),
Xcand = as.double(t(Xcand)),
ncand = as.integer(ncand),
Xref = as.double(t(Xref)),
nref = as.integer(nrow(Xref)),
fi = as.integer(fi),
verb = as.integer(verb),
mspes = double(ncand),
PACKAGE = "laGP")
return(out$mspes)
}
## dmus2GP:
##
## obtain the derivative of the predictive scale
## of the fitted GP model
dmus2GP <- function(gpi, XX)
{
nn <- nrow(XX)
out <- .C("dmus2GP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(XX)),
nn = as.integer(nn),
XX = as.double(t(XX)),
mu = double(nn),
dmu = double(nn),
d2mu = double(nn),
s2 = double(nn),
ds2 = double(nn),
d2s2 = double(nn),
PACKAGE = "laGP")
return(data.frame(mu=out$mu, dmu=out$dmu, d2mu=out$d2mu,
s2=out$s2, ds2=out$ds2, d2s2=out$d2s2))
}
## fishGP:
##
## obtain the expected (approx) Fisher information for
## the fitted GP model; returns the absolute value (i.e.,
## determinant)
fishGP <- function(gpi, Xcand)
{
nn <- nrow(Xcand)
out <- .C("efiGP_R",
gpi = as.integer(gpi),
m = as.integer(ncol(Xcand)),
nn = as.integer(nn),
Xcand = as.double(t(Xcand)),
efi = double(nn),
PACKAGE = "laGP")
## remove silly values
return(out$efi)
}
## getmGP:
##
## access the input dimension of a GP
getmGP <- function(gpi)
{
.C("getmGP_R", gpi = as.integer(gpi), m = integer(1), PACKAGE="laGP")$m
}
Try the laGP package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.