# R/gp.R In laGP: Local Approximate Gaussian Process Regression

#### Documented in alcGPalcoptGPalcrayGPdalcGPdeleteGPdeleteGPsfishGPieciGPjmleGPjmleGP.RllikGPmleGPmspeGPnewGPpredGPupdateGP

```#*******************************************************************************
#
# 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
# 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))

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") {

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)

## 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")

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")

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.