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

#### Documented in alcGPsepalcoptGPsepalcrayGPsepdalcGPsepdeleteGPsepdeleteGPsepsieciGPsepjmleGPsepjmleGPsep.RllikGPsepmleGPsepmleGPsep.RnewGPseppredGPsepupdateGPsep

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

## newGPsep:
##
## build an initial separable GP representation on the C-side
## using the X-Z data and d/g paramterization.

newGPsep <- function(X, Z, d, g, dK=FALSE)
{
n <- nrow(X)
m <- ncol(X)
if(is.null(n)) stop("X must be a matrix")
if(length(Z) != n) stop("must have nrow(X) = length(Z)")
if(length(d) == 1) d <- rep(d, m)
else if(length(d) != m) stop("must have length(d) = ncol(X)")

out <- .C("newGPsep_R",
m = as.integer(m),
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),
gpsepi = integer(1),
PACKAGE = "laGP")

## return C-side GP index
return(out\$gpsepi)
}

## buildkGPsep:
##
## allocates/calculates the C-side derivative info (only) for particular
## separable GP

buildKGPsep <- function(gpsepi)
{
.C("buildKGPsep_R",
gpisep = as.integer(gpsepi),
PACKAGE = "laGP")
invisible(NULL)
}

## deletedkGPsep:
##
## deletes the C-side derivative info (only) for particular separable GP

deletedkGPsep <- function(gpsepi)
{
.C("deletedKGPsep_R",
gpi = as.integer(gpsepi),
PACKAGE = "laGP")
invisible(NULL)
}

## deleteGPsep:
##
## deletes the C-side of a particular separable GP

deleteGPsep <- function(gpsepi)
{
.C("deleteGPsep_R",
gpsepi = as.integer(gpsepi), PACKAGE="laGP")
invisible(NULL)
}

## deleteGPseps:
##
## deletes all gpseps on the C side

deleteGPseps <- function()
{
.C("deleteGPseps_R", PACKAGE="laGP")
invisible(NULL)
}

## llikGPsep:
##
## calculate the log likelihood of the GP

llikGPsep <- function(gpsepi, dab=c(0,0), gab=c(0,0))
{
r <- .C("llikGPsep_R",
gpsepi = as.integer(gpsepi),
dab = as.double(dab),
gab = as.double(gab),
llik = double(1),
PACKAGE = "laGP")

return(r\$llik)
}

## getmGPsep:
##
## acces the input dimension of a separable GP
##
## totally new to GPsep

getmGPsep <- function(gpsepi)
{
.C("getmGPsep_R", gpsepi = as.integer(gpsepi), m = integer(1), PACKAGE="laGP")\$m
}

## getdGPsep:
##
## acces the separable lengthscale of a separable gp
##
## totally new to GPsep

getdGPsep <- function(gpsepi)
{
m <- getmGPsep(gpsepi)
.C("getdGPsep_R", gpsepi = as.integer(gpsepi), d = double(m), PACKAGE="laGP")\$d
}

## getgGPsep:
##
## acces the input dimension of a separable GP
##
## totally new to GPsep

getgGPsep <- function(gpsepi)
{
.C("getgGPsep_R", gpsepi = as.integer(gpsepi), g = double(1), PACKAGE="laGP")\$g
}

## dllikGPsep:
##
## calculate the first and second derivative of the
## log likelihood of the GP with respect to d, the
## lengthscale parameter
##
## SIMILAR to dllikGP except with vector d and gpsep
## isntead of gp

dllikGPsep <- function(gpsepi, ab=c(0,0), param=c("d", "g"), d2nug=FALSE)
{
param <- match.arg(param)
if(param == "d") {
dim <- getmGPsep(gpsepi)
r <- .C("dllikGPsep_R",
gpsepi = as.integer(gpsepi),
ab = as.double(ab),
d = double(dim),
PACKAGE = "laGP")
return(r\$d)
} else {
if(d2nug) d2 <- 1
else d2 <- 0
r <- .C("dllikGPsep_nug_R",
gpsepi = as.integer(gpsepi),
ab = as.double(ab),
d = double(1),
d2 = as.double(d2),
PACKAGE = "laGP")
if(d2nug) return(list(d=r\$d, d2=r\$r2))
else return(r\$d)
}
}

## newparamsGPsep:
##
## change the separable GP lengthscale and nugget parameerization
## (without destroying the object and creating a new one)

newparamsGPsep <- function(gpsepi, d, g=-1)
{
if(all(d <= 0) & g < 0) stop("one of d or g must be new")
m <- getmGPsep(gpsepi)
if(length(d) != m) stop("length(d) !=", m)

r <- .C("newparamsGPsep_R",
gpi = as.integer(gpsepi),
d = as.double(d),
g = as.double(g),
PACKAGE = "laGP")

invisible(NULL)
}

## mleGPsep.R:
##
## updates the separable GP to use its MLE lengthscale
## parameterization using the current data;
##
## differs substantially from mleGP in that L-BFGS-B from
## optim is used to optimize over the separable lengthscale;
## an option is also provided to include the nugget in that
## optimization, or do to a mleGP style profile optimization
## for the the nugget instead
##
## in this ".R" version the optim command is used; in mleGPsep
## an internal C-side call to lbfgsb is used

mleGPsep.R <- function(gpsepi, param=c("d", "g"),
tmin=sqrt(.Machine\$double.eps),
tmax=-1, ab=c(0,0), maxit=100, verb=0)
{
param <- match.arg(param)

if(param == "d") { ## lengthscale with L-BFGS-B given nugget

theta <- getdGPsep(gpsepi)
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector")

## objective
f <- function(theta, gpsepi, dab)
{
newparamsGPsep(gpsepi, d=theta)
-llikGPsep(gpsepi, dab=dab)
}
## gradient of objective
g <- function(theta, gpsepi, dab)
{
newparamsGPsep(gpsepi, d=theta)
-dllikGPsep(gpsepi, param="d", ab=dab)
}

## for compatibility with mleGP
tmax[tmax < 0] <- Inf

## possibly print progress meter
if(verb > 0) {
cat("(d=[", paste(signif(theta, 3), collapse=","), "], llik=",
llikGPsep(gpsepi, dab=ab), ") ", sep="")
}

## call R's optim function
out <- optim(theta, fn=f, gr=g, method="L-BFGS-B",
control=list(trace=max(verb-1,0), maxit=maxit), lower=tmin, upper=tmax,
gpsepi=gpsepi, dab=ab)

## sanity check completion of scheme
if(sqrt(mean((out\$par - getdGPsep(gpsepi))^2)) > sqrt(.Machine\$double.eps))
warning("stored d not same as d-hat")

## check that we moved somewhere
if(sqrt(mean((out\$par - theta)^2)) < sqrt(.Machine\$double.eps)) {
out\$convergence <- 0
out\$counts <- c(0,0)
out\$message <- "optim initialized at minima"
}

## possibly print progress meter
if(verb > 0) {
cat("-> ", out\$counts[1], " lbfgs.R its -> (d=[",
paste(signif(theta, 3), collapse=","), "], llik=",
llikGPsep(gpsepi, dab=ab), ")\n", sep="")
}

}
else { ## nugget conditionally on lengthscale

## sanity check
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector");

r <- .C("mleGPsep_nug_R",
gpsepi = as.integer(gpsepi),
verb = as.integer(verb),
tmin = as.double(tmin),
tmax = as.double(tmax),
ab = as.double(ab),
g = double(1),
its = integer(1),
PACKAGE = "laGP")
}

## build object for returning
if(param == "d") return(list(d=out\$par, its=max(out\$counts), msg=out\$message, conv=out\$convergence))
else return(list(g=r\$g, its=r\$its))
}

## mleGPsep:
##
## updates the separable GP to use its MLE lengthscale
## parameterization using the current data;
##
## differs substantially from mleGP in that lbfgsb is used
## to optimize over the separable lengthscale;
## an option is also provided to include the nugget in that
## optimization, or do to a mleGP style profile optimization
## for the the nugget instead
##
## this is a mostly C verision

mleGPsep <- function(gpsepi, param=c("d", "g", "both"),
tmin=rep(sqrt(.Machine\$double.eps), 2),
tmax=c(-1,1), ab=rep(0,4), maxit=100, verb=0)
{
param <- match.arg(param)

if(param == "both") { ## lengthscale and multiple nugget jointly with L-BFGS-B

## sanity checking, and length checking for tmax and tmin
m <- getmGPsep(gpsepi)
if(length(tmax) == 2) tmax <- c(rep(tmax[1], m), tmax[2])
else if(length(tmax) != m+1) stop("length(tmax) should be 2 or m+1")
if(length(tmin) == 2) tmin <- c(rep(tmin[1], m), tmin[2])
else if(length(tmin) != m+1) stop("length(tmin) should be 2 or m+1")
if(length(ab) != 4 || any(ab < 0)) stop("ab should be a positive 4-vector")

## possibly reset params
theta <- c(getdGPsep(gpsepi), getgGPsep(gpsepi))
if(any(theta <= tmin)) {
tmax[tmax < 0] <- sqrt(m)
theta.new <- 0.9*max(tmin, 0) + 0.1*tmax
newparamsGPsep(gpsepi, theta.new[1:m], theta.new[m+1])
return(list(theta=theta.new, its=0, msg="reset due to init on lower boundary", conv=102))
}

out <- .C("mleGPsep_both_R",
gpsepi = as.integer(gpsepi),
maxit = as.integer(maxit),
verb = as.integer(verb),
tmin = as.double(tmin),
tmax = as.double(tmax),
ab = as.double(ab),
par = double(m+1),
counts = integer(2),
msg = paste(rep(0,60), collapse=""),
convergence = integer(1),
PACKAGE = "laGP")

## sanity check completion of scheme
if(sqrt(mean((out\$par - c(getdGPsep(gpsepi), getgGPsep(gpsepi)))^2)) > sqrt(.Machine\$double.eps))
warning("stored theta not same as theta-hat")

} else if(param == "d") { ## lengthscale with L-BFGS-B given nugget

## sanity checking
m <- getmGPsep(gpsepi)
if(length(tmax) == 1 || (length(tmax) == 2 && m != 2)) tmax <- rep(tmax[1], m)
else if(length(tmax) != m) stop("length(tmax) should be 1 or m")
if(length(tmin) == 1 || (length(tmin) == 2 && m != 2)) tmin <- rep(tmin[1], m)
else if(length(tmin) != m) stop("length(tmin) should be 1 or m")
if(length(ab) == 4 && all(ab == 0)) ab <- ab[1:2]
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector")

## possibly reset params
theta <- getdGPsep(gpsepi)
if(any(theta <= tmin)) {
tmax[tmax < 0] <- sqrt(m)
theta.new <- 0.9*tmin + 0.1*tmax
newparamsGPsep(gpsepi, theta.new, -1)
return(list(d=theta.new, its=0, msg="reset due to init on lower boundary", conv=102))
}

out <- .C("mleGPsep_R",
gpsepi = as.integer(gpsepi),
maxit = as.integer(maxit),
verb = as.integer(verb),
dmin = as.double(tmin),
dmax = as.double(tmax),
ab = as.double(ab),
par = double(m),
counts = integer(2),
msg = paste(rep(0,60), collapse=""),
convergence = integer(1),
PACKAGE = "laGP")

## sanity check completion of scheme
if(sqrt(mean((out\$par - getdGPsep(gpsepi))^2)) > sqrt(.Machine\$double.eps))
warning("stored d not same as theta-hat")
}
else { ## nugget conditionally on lengthscale

## sanity check
if(length(ab) == 4 && all(ab == 0)) ab <- ab[1:2]
if(length(ab) != 2 || any(ab < 0)) stop("ab should be a positive 2-vector");

r <- .C("mleGPsep_nug_R",
gpsepi = as.integer(gpsepi),
verb = as.integer(verb),
tmin = as.double(tmin),
tmax = as.double(tmax),
ab = as.double(ab),
g = double(1),
its = integer(1),
PACKAGE = "laGP")
}

## build object for returning
if(param == "both") return(list(theta=out\$par, its=max(out\$counts), msg=out\$msg, conv=out\$convergence))
else if(param == "d") return(list(d=out\$par, its=max(out\$counts), msg=out\$msg, conv=out\$convergence))
else return(list(g=r\$g, its=r\$its))
}

## jmleGPsep.R:
##
## joint MLE for lengthscale (d) and nugget (g) parameters;
## updates the internal GP parameterization (since mleGP does);
## R-only version

jmleGPsep.R <- function(gpsepi, N=100, drange=c(sqrt(.Machine\$double.eps), 10),
grange=c(sqrt(.Machine\$double.eps), 1), dab=c(0,0), gab=c(0,0), maxit=100,
mleGPsep=mleGPsep.R, verb=0)
{
## sanity check N
if(length(N) != 1 && N > 0)
stop("N should be a positive scalar integer")
m <- getmGPsep(gpsepi)
dmle <- matrix(NA, nrow=N, ncol=m)
gmle <- dits <- dconv <- 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 <- mleGPsep(gpsepi, param="d", tmin=drange[1], tmax=drange[2],
ab=dab, maxit=maxit, verb=verb)
dmle[i,] <- d\$d; dits[i] <- d\$its; dconv[i] <- d\$conv
g <- mleGPsep(gpsepi, param="g", tmin=grange[1], tmax=grange[2],
ab=gab, verb=verb)
gmle[i] <- g\$g; gits[i] <- g\$its
if((gits[i] <= 2 && (dits[i] <= m+1 && dconv[i] == 0)) || dconv[i] > 1) break;
}

## check if not converged
if(i == N) warning("max outer its (N=", N, ") reached", sep="")
else {
dmle <- dmle[1:i,]; dits <- dits[1:i]; dconv <- dconv[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,,drop=FALSE], g=gmle[i], tot.its=totits,
conv=dconv[i]), prog=data.frame(dmle=dmle, dits=dits, dconv=dconv, gmle=gmle,
gits=gits)))
}

## jmleGPsep
##
## interface to C-version for jmleGPsep;
## right now doesn't take an N argument -- the C-side hard-codes
## N=100

jmleGPsep <- function(gpsepi, drange=c(sqrt(.Machine\$double.eps), 10),
grange=c(sqrt(.Machine\$double.eps), 1), dab=c(0,0), gab=c(0,0), maxit=100,
verb=0)
{
## sanity check tmin and tmax
m <- getmGPsep(gpsepi)
if(length(drange) != 2) stop("drange should be a two vector for c(dmin, dmax)")
dmin <- rep(drange[1], m)
dmax <- rep(drange[2], m)
if(length(grange) != 2) stop("grange should be a 2-vector for c(gmin, gmax)")

## 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("jmleGPsep_R",
gpsepi = as.integer(gpsepi),
maxit = as.integer(maxit),
verb = as.integer(verb),
dmin = as.double(dmin),
dmax = as.double(dmax),
grange = as.double(grange),
dab = as.double(dab),
gab = as.double(gab),
d = double(m),
g = double(1),
dits = integer(1),
gits = integer(1),
dconv = integer(1),
PACKAGE = "laGP")

return(data.frame(d=t(r\$d), g=r\$g, tot.its=r\$dits+r\$gits,
dits=r\$dits, gits=r\$gits, dconv=r\$dconv))
}

## predGPsep
##
## obtain the parameters to a multivariate-t
## distribution describing the predictive surface
## of the fitted GP model

predGPsep <- function(gpsepi, XX, lite=FALSE, nonug=FALSE)
{
nn <- nrow(XX)
if(is.null(nn) || nn == 0) stop("XX bad dims")

if(lite) {  ## lite means does not compute full Sigma, only diag
out <- .C("predGPsep_R",
gpsepi = as.integer(gpsepi),
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 { ## compute full predictive covariance matrix

out <- .C("predGPsep_R",
gpsepi = as.integer(gpsepi),
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))
}
}

## updateGPsep:
##
## add X-Z pairs to the C-side GPsep represnetation
## using only O(n^2) for each pair

updateGPsep <- function(gpsepi, X, Z, verb=0)
{
if(length(Z) != nrow(X))

out <- .C("updateGPsep_R",
gpsepi = as.integer(gpsepi),
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)
}

## alGPsep:
##
## calculate the E(Y) and EI(Y) for an augmented Lagrangian
## composite objective function with linear objective (in X), or
## estimate objective (fhat) and constraint separable GP (gpsepi)
## predictive surfaces

alGPsep <- function(XX, fgpsepi, fnorm, Cgpsepis, Cnorms, lambda, alpha, ymin,
slack=FALSE, equal=rep(FALSE, length(Cgpsepis)), N=100, fn=NULL, Bscale=1)
{
## dims
m <- ncol(XX)
nn <- nrow(XX)
nCgpseps <- length(Cgpsepis)

## checking lengths for number of constraint gps
if(length(Cnorms) != nCgpseps) stop("length(Cgpsepis) != length(Cnorms)")
if(length(lambda) != nCgpseps) stop("length(Cgpsepis) != length(lambda)")
if(length(alpha) != 1) stop("length(alpha) != 1")

## checking scalars
if(length(equal) != length(Cgpsepis)) stop("equal should be a vector of length(Cgpsepis)")
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(fgpsepi < 0 || any(Cgpsepis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpsepi or Cgpsepis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpsepi < 0) {
if(is.null(out\$obj)) stop("fgpsepi < 0 but out\$obj from fn() is NULL")
obj <- out\$obj
} else obj <- NULL
if(any(Cgpsepis < 0)) {
C <- out\$c
if(ncol(C) != sum(Cgpsepis < 0)) stop("ncol(C) != sum(Cgpsepis < 0)")
} else C <- NULL
} else { obj <- C <- NULL }

## call the C-side
out <- .C("alGPsep_R",
m = as.integer(m),
XX = as.double(t(XX)),
nn = as.integer(nn),
fgpsepi = as.integer(fgpsepi),
ff = as.double(obj),
fnorm = as.double(fnorm),
nCgpseps = as.integer(nCgpseps),
Cgpsepis = as.integer(Cgpsepis),
CC = as.double(C),
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))
}

## efiGPsep:
##
## calculate EI(f) and p(Y(c) <= 0) for known linear or esitmated
## objective f and vectorized constraints C via separable GP (gpsepi)
## predictive surfaces; returns log probabilities (lplex) and
## and log EIs

efiGPsep <- function(XX, fgpsepi, fnorm, Cgpsepis, Cnorms, fmin, fn=NULL, Bscale=1)
{
## doms
m <- ncol(XX)
nn <- nrow(XX)
nCgpseps <- length(Cgpsepis)

## checking lengths for number of constraint gps
if(length(Cnorms) != nCgpseps) stop("length(Cgpsepis) != 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(fgpsepi < 0 || any(Cgpsepis < 0)) {
if(is.null(fn)) stop("fn must be provided when fgpsepi or Cgpsepis < -1")
out <- fn(XX*Bscale, known.only=TRUE)
if(fgpsepi < 0) {
if(is.null(out\$obj)) stop("fgpsepi < 0 but out\$obj from fn() is NULL")
obj <- out\$obj
} else obj <- NULL
if(any(Cgpsepis < 0)) {
C <- out\$c
if(ncol(C) != sum(Cgpsepis < 0)) stop("ncol(C) != sum(Cgpsepis < 0)")
} else C <- NULL
}

## calculate expected improvement part
if(fgpsepi < 0) {
obj <- rowSums(XX) * fnorm
if(!is.finite(fmin)) fmin <- quantile(obj, p=0.9)
I <- fmin - obj
ei <- pmax(I, 0)
} else {
p <- predGPsep(fgpsepi, 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), nCgpseps)
ik <- 1
for(j in 1:nCgpseps) {
if(Cgpsepis[j] < 0) {
lplez[,j] <- log(C[,ik] <= 0)
ik <- ik + 1
} else {
pc <- predGPsep(Cgpsepis[j], XX=XX, lite=TRUE)
lplez[,j] <- pnorm(0, pc\$mean, sqrt(pc\$s2), log.p=TRUE)
}
}

## done
return(data.frame(lei=log(ei), lplez=lplez))
}

## alcGPsep:
##
## wrapper used to calculate the ALCs in C using
## the pre-stored separable GP representation.
## Note that this only calculates the s2' component
## of ds2 = s2 - s2'

alcGPsep <- function(gpsepi, 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("alcGPsep_omp_R")) stop("OMP not supported in this build; please re-compile")

out <- .C("alcGPsep_omp_R",
gpsepi = as.integer(gpsepi),
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") { stop("alcGPsep_gpu_R not implemented")
} else {

out <- .C("alcGPsep_R",
gpsepi = as.integer(gpsepi),
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)
}

## alcoptGPsep:
##
## interface to C version of alcoptGPsep.R which continuously optimizes
## ALC based on derivatives, using the starting locations and bounding
## boxes and (stored) gpsep provided; ... has arguments to optim including
## trace/verb level

alcoptGPsep <- function(gpsepi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGPsep(gpsepi)
if(ncol(Xref) != m) stop("gpsepi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpsepi 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("alcoptGPsep_R",
gpsepi = as.integer(gpsepi),
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))
}

## alcoptGPsep.R:
##
## continuously optimizes ALC based on derivatives, using the
## starting locations and (stored) gpsep provided; ... has arguments
## to optim including trace/verb level

alcoptGPsep.R <- function(gpsepi, Xref, start, lower, upper, maxit=100, verb=0)
{
m <- getmGPsep(gpsepi)
if(ncol(Xref) != m) stop("gpsepi stored X and Xref have mismatched cols")
if(length(start) != m) stop("gpsepi stored X and start have mismatched cols")

## objective (and derivative saved)
deriv <- NULL
f <- function(x, gpsepi, Xref) {
out <- dalcGPsep(gpsepi, 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, gpsepi, 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, gpsepi=gpsepi, Xref=Xref, lower=lower, upper=upper,
method="L-BFGS-B", control=control)
## version without derivatives
## opt <- optim(start, f, gpsepi=gpsepi, Xref=Xref, lower=lower, upper=upper,
## method="L-BFGS-B", control=control)

## keep track of progress in derivatives
## f(opt\$par, gpsepi, Xref)
## grads <<- rbind(grads, df(opt\$par, gpsepi, Xref))

## for now return the whole optim output
return(opt)
}

## lalcoptGPsep.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
## lalcoptGPsep, since the starts are random from 1:offset

lalcoptGPsep.R <- function(gpsepi, 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 <- alcoptGPsep(gpsepi, Xref, Xstart[i,], rect[1,], rect[2,], verb=verb)
## opt <- alcoptGPsep.R(gpsepi, 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 <- alcGPsep(gpsepi, 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)
}

## lalcoptGPsep:
##
## 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 lalcoptGPsep.R, since the starts are
## determined by a deterministic round robin similar to lalcrayGPsep

lalcoptGPsep <- function(gpsepi, 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("lalcoptGPsep_R",
gpsepi = as.integer(gpsepi),
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)
}

## getmGPsep:
##
## access the input dimension of a GPsep

getmGPsep <- function(gpsepi)
{
.C("getmGPsep_R", gpsepi = as.integer(gpsepi), m = integer(1), PACKAGE="laGP")\$m
}

## dalcGPsep:
##
## wrapper used to calculate the derivative of ALCs in C using
## the pre-stored separable GP representation.  Note that this only
## calculates the s2' component of ds2 = s2 - s2'

dalcGPsep <- function(gpsepi, Xcand, Xref=Xcand, verb=0)
{
m <- ncol(Xcand)
if(ncol(Xref) != m) stop("Xcand and Xref have mismatched cols")
ncand <- nrow(Xcand)

out <- .C("dalcGPsep_R",
gpsepi = as.integer(gpsepi),
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)))
}

## ieciGPsep:
##
## wrapper used to calculate the IECIs in C using
## the pre-stored separable GP representation.

ieciGPsep <- function(gpsepi, 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("ieciGPsep_R",
gpsepi = as.integer(gpsepi),
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)
}

## mleGPsep.switch:
##
## switch function for mle calculaitons by laGPsep.R

mleGPsep.switch <- function(gpsepi, 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 != "efi") buildKGPsep(gpsepi)

if(d\$mle && g\$mle) { ## joint lengthscale and nugget
return(jmleGPsep(gpsepi, 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 <- mleGPsep(gpsepi, param="d", d\$min, d\$max, d\$ab, verb=verb)
return(data.frame(d=matrix(dmle\$d, nrow=1), dits=dmle\$its))
}
if(g\$mle) { ## nugget only
gmle <- mleGPsep(gpsepi, param="g", g\$min, g\$max, g\$ab, verb=verb)
return(data.frame(g=gmle\$g, gits=gmle\$its))
}
}
}

## alcrayGPsep:
##
## wrapper used to optimize AIC via a ray search using
## the pre-stored separable GP representation.  Return
## the convex combination s in (0,1) between Xstart and Xend;

alcrayGPsep <- function(gpsepi, 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("alcrayGPsep_R",
gpsepi = as.integer(gpsepi),
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))

## return the convex combination
return(list(r=out\$r, s=out\$s))
}

## lalcrayGPsep.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
## lalcrayGPsep, since the starts of the rays are random from
## 1:offset -- otherwise this is nearly identical to lalcrayGP.R

lalcrayGPsep.R <- function(gpsepi, 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 <- alcrayGPsep(gpsepi, 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)
}

## lalcrayGPsep:
##
## 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 alcrayGPsep (on the C-side) is called to
## optimize over the ray.  The candidate in Xcand which is closest
## to the solution is returned -- nearly identical to lalcrayGP

lalcrayGPsep <- function(gpsepi, 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("lalcrayGPsep_R",
gpsepi = as.integer(gpsepi),
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)
}
```

## 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.