R/gp.R

#*******************************************************************************
#
# 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 ([email protected])
#
#*******************************************************************************


## 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(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),
              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, 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))
    
    ## 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 May 2, 2019, 9:20 a.m.