# R/solve_game_baseline.R In vpicheny/GPGame: Solving Complex Game Problems using Gaussian Processes

#### Defines functions solve_game_baseline

```#' @noRd
#' @export
solve_game_baseline <- function(
fun, ..., equilibrium="NE", crit="sur", model=NULL, n.init=NULL, n.ite, d, nobj, x.to.obj=NULL, noise.var = NULL,
ncores=1, trace=1, seed=NULL, calibcontrol=NULL, freq.exploit=1e3, baseline_type="1d") {

t1 <- Sys.time()
set.seed(seed)
if(ncores > 1) parallel <- TRUE

####################################################################################################
#### CHECK INPUTS ##################################################################################
####################################################################################################
# Check integcontrol
integ.pts <- integcontrol\$integ.pts
expanded.indices <- integcontrol\$expanded.indices
n.s <- integcontrol\$n.s
init_set <- integcontrol\$init_set
gridtype <- switch(1+is.null(integcontrol\$gridtype), integcontrol\$gridtype, "lhs")
lb <- switch(1+is.null(integcontrol\$lb), integcontrol\$lb, rep(0, d))
ub <- switch(1+is.null(integcontrol\$ub), integcontrol\$ub, rep(1, d))
integcontrol\$renew <- switch(1+is.null(integcontrol\$renew), integcontrol\$renew, equilibrium %in% c("KSE", "CKSE"))

# Check simucontrol
n.ynew <- switch(1+is.null(simucontrol\$n.ynew), simucontrol\$n.ynew, 10)
n.sim <- switch(1+is.null(simucontrol\$n.sim), simucontrol\$n.sim, 10)
IS <- switch(1+is.null(simucontrol\$IS), simucontrol\$IS, TRUE)
cross <- FALSE

# Check kmcontrol
cov.reestim <- switch(1+is.null(kmcontrol\$cov.reestim), kmcontrol\$cov.reestim, TRUE)
model.trend <- switch(1+is.null(kmcontrol\$model.trend), kmcontrol\$model.trend, ~1)
kmlb <- switch(1+is.null(kmcontrol\$lb), kmcontrol\$lb, rep(.1,d))
kmub <- switch(1+is.null(kmcontrol\$ub), kmcontrol\$ub, rep(1,d))
kmnugget <- switch(1+is.null(kmcontrol\$nugget), kmcontrol\$nugget, 1e-8)
control <- switch(1+is.null(kmcontrol\$control), kmcontrol\$control, list(trace=FALSE))
covtype <- switch(1+is.null(kmcontrol\$covtype), kmcontrol\$covtype, "matern5_2")

# Check Noise
if (!is.null(noise.var)) {
if (typeof(noise.var) == "closure") noise.var <- match.fun(noise.var)
else if (typeof(noise.var) == "double" && length(noise.var==1)) noise.var <- rep(noise.var, nobj)
kmnugget <- NULL
}

# Check calibcontrol
if (!is.null(calibcontrol)) {
if (is.null(calibcontrol\$target)){
cat("calibcontrol should contain a target vector \n")
return(NA)
}
if (is.null(calibcontrol\$log)) calibcontrol\$log <- FALSE
if (is.null(calibcontrol\$offset)) calibcontrol\$offset <- 0
}

if (trace>0) cat("--------------------------\n", " Starting", equilibrium, "\n",
" among (", n.s, ") strategies \n", "--------------------------\n", sep = "")

####################################################################################################
#### INITIALIZE VARIABLES AND MODELS ###############################################################
####################################################################################################
if (baseline_type == "RS") {
if (is.null(n.init))     n.init <- 0
n.init <- n.init + n.ite
n.ite <- 0
}

#### Initial design and models ####################
if (is.null(model)){
if (is.null(n.init))     n.init <- 5*d
design <- lhsDesign(n.init, d, seed = seed)\$design
response <- t(apply(design, 1, fun, ... = ...))

if (!is.null(noise.var)) {
if (typeof(noise.var) == "closure") {
newnoise.var <- apply(design, 1, noise.var, ...)
} else if (typeof(noise.var) == "double") {
newnoise.var <- matrix(noise.var, nrow=n.init, ncol=ncol(response), byrow=TRUE)
} else {
tmp <- newnoise.var <- NULL # initialization
for (i in 1:length(response)) {
tmp <- rbind(tmp, response[[i]][[1]])
newnoise.var <- rbind(newnoise.var, response[[i]][[2]])
}
response <- tmp
}
} else {
newnoise.var <- NULL
}
nobj <- ncol(response)

my.km <- function(i) {
km(model.trend, design = design, response = response[, i], covtype=covtype,
control=control, lower=kmlb, upper=kmub, nugget=kmnugget, noise.var=newnoise.var[,i])
}
model <- mclapply(1:nobj, my.km, mc.cores=ncores)
}

#### Integration points ###########################
if (is.null(integcontrol\$integ.pts) || is.null(integcontrol\$expanded.indices)) {
if (equilibrium %in% c("KSE", "CKSE")) include.obs <- TRUE else include.obs <- FALSE
res <- generate_integ_pts(n.s=n.s, d=d, nobj=nobj, x.to.obj=x.to.obj, equilibrium=equilibrium,
gridtype=gridtype, lb = lb, ub = ub, include.obs=include.obs, model=model)
integcontrol\$integ.pts <- integ.pts <- res\$integ.pts
integcontrol\$expanded.indices <- expanded.indices <- res\$expanded.indices
}
if (is.null(expanded.indices)) sorted <- FALSE
else sorted <- !is.unsorted(expanded.indices[, ncol(expanded.indices)])

####################################################################################################
#### DEFINE SEQUENCE OF ACTIONS ####################################################################
####################################################################################################
# J is the individual objective to consider at time ii (0 means all => equilibrium)
J <- rep(c(1, 1:nobj, 1:nobj), ceiling((n.ite+1) / (2*nobj-1)))

if (equilibrium == "KSE") {
task <- rep(c(0, rep(1, nobj), rep(2, nobj)), ceiling((n.ite+1) / (2*nobj-1)))
} else if (equilibrium == "CKSE") {
task <- rep(c(0, rep(3, nobj), rep(3, nobj)), ceiling((n.ite+1) / (2*nobj-1)))
} else {
equilibrium <- "KSE" # Could be improved, e.g., by putting sms in baseline_type
}

####################################################################################################
#### MAIN LOOP STARTS HERE #########################################################################
####################################################################################################

for (ii in 1:(n.ite+1)){

##################################################################
# RENEW INTEGRATION POINTS
##################################################################
if (ii > 1 && integcontrol\$renew) {
include.obs <- equilibrium %in% c("KSE", "CKSE")
res <- generate_integ_pts(n.s=n.s, d=d, nobj=nobj, x.to.obj=x.to.obj, equilibrium=equilibrium,
gridtype=gridtype, lb = lb, ub = ub, include.obs=include.obs, model=model,
init_set=init_set, seed=ii)
integcontrol\$integ.pts <- integ.pts <- res\$integ.pts
integcontrol\$expanded.indices <- expanded.indices <- res\$expanded.indices
}

##################################################################
# MODELS UPDATE
##################################################################
if (ii > 1) {
newmodel <- model
X.new <- matrix(xnew,nrow=1)
my.update <- function(u) {
try(update(object = model[[u]], newX = X.new, newy=ynew[u], newX.alreadyExist=FALSE, newnoise.var = newnoise.var[u],
cov.reestim = cov.reestim, kmcontrol = list(control = list(trace = FALSE))), silent = TRUE)
}
newmodel <- mclapply(1:nobj, my.update, mc.cores=ncores)

for (u in 1:nobj){
if (typeof(newmodel[[u]]) == "character" && cov.reestim) {
cat("Error in hyperparameter estimation - old hyperparameter values used instead for model ", u, "\n")
newmodel[[u]] <- try(update(object = model[[u]], newX = X.new, newy=ynew[u], newnoise.var = newnoise.var[u],
newX.alreadyExist=FALSE, cov.reestim = FALSE), silent = TRUE)
}
if (typeof(newmodel[[u]]) == "character") {
cat("Unable to udpate kriging model ", u, " at iteration", ii-1, "- optimization stopped \n")
cat("lastmodel ", u, " is the model at iteration", ii-1, "\n")
cat("par and values contain the ",ii , "th observation \n \n")
return(list(par = X.new, values = ynew, nsteps = ii, model = model, integcontrol=integcontrol))
} else {
model[[u]] <- newmodel[[u]]
}
}
}
observations <- Reduce(cbind, lapply(model, slot, "y"))

##################################################################
# PREDICTION
##################################################################
if (nrow(integ.pts) < 11e5) {
pred <- mclapply(model, FUN=predict, newdata = integ.pts, checkNames = FALSE, type = "UK", light.return = TRUE, mc.cores=ncores)
predmean <- Reduce(cbind, lapply(pred, function(alist) alist\$mean))
predsd <- Reduce(cbind, lapply(pred, function(alist) alist\$sd))
} else {
start <- seq(1, nrow(integ.pts), 2e5)
stop <- seq(2e5, nrow(integ.pts), 2e5)
if (length(start) > length(stop)) stop <- c(stop, nrow(integ.pts))

predmean <- predsd <- c()
for (b in 1:length(start)) {
integ.pts_red <- integ.pts[start[b]:stop[b], ]
pred <- mclapply(model, FUN=predict, newdata = integ.pts_red, checkNames = FALSE, type = "UK", light.return = TRUE, mc.cores=ncores)
predmean <- rbind(predmean, Reduce(cbind, lapply(pred, function(alist) alist\$mean)))
predsd <- rbind(predsd, Reduce(cbind, lapply(pred, function(alist) alist\$sd)))
}
}

##################################################################
# ACQUISITION
##################################################################
jj <- J[ii]

if (max_predvar) {
################################################################
# Maximise prediction variance
i <- which.max(predsd[,jj])

################################################################
if (trace>1) cat("Looking for individual minimum of objective: ", jj, "\n")
if (is.null(calibcontrol\$target)) {
# EI(min) on each objective (to find potential Utopia points)
xcr <-  (min(model[[jj]]@y) - predmean[,jj])/predsd[,jj]
test <- (min(model[[jj]]@y) - predmean[,jj])*pnorm(xcr) + predsd[,jj] * dnorm(xcr)
} else {
# EI(min) on each objective (to find potential Utopia points)
zmin <- min((model[[jj]]@y - calibcontrol\$target[jj])^2)
mu <- predmean[,jj] - calibcontrol\$target[jj]
sigma   <- predsd[,jj]
a2 <- (sqrt(zmin) - mu)/sigma
a1 <- (-sqrt(zmin) - mu)/sigma
da2 <- dnorm(a2);   da1 <- dnorm(a1)
pa2 <- pnorm(a2);   pa1 <- pnorm(a1)
test <- (zmin - sigma^2 - mu^2) * (pa2 - pa1) + sigma^2 * (a2*da2 - a1*da1) + 2*mu*sigma*(da2 - da1)
}
i <- which.max(test)

################################################################
if (trace>1) cat("Looking for nadir of objective: ", jj, "\n")
PFobs <- nonDom(Reduce(cbind, lapply(model, slot, "y")))
PNDom <- prob.of.non.domination(model=model, integration.points=integ.pts, predictions=pred, target=calibcontrol\$target)

if (is.null(calibcontrol\$target)) {
# EI(max) x Pdom on each objective (to find potential Nadir points) unless Nadir is provided
xcr <-  -(max(PFobs[,jj]) - predmean[,jj])/predsd[,jj]
test <- (-max(PFobs[,jj]) + predmean[,jj])*pnorm(xcr) + predsd[,jj] * dnorm(xcr)
} else {
zmax <- max(PFobs[,jj])
b2 <- (sqrt(zmax) - mu)/sigma
b1 <- (-sqrt(zmax) - mu)/sigma
db2 <- dnorm(b2);   db1 <- dnorm(b1)
pb2 <- pnorm(b2);   pb1 <- pnorm(b1)
test <- (sigma^2 + mu^2 - zmax) * (1 - pb2 + pb1) + sigma^2 * (b2*db2 - b1*db1) + 2*mu*sigma*(db2 - db1)
}
test <- test * PNDom
i <- which.max(test)

################################################################
# SMS-driven
} else if (sms) {
if (trace>1) cat("Looking for max of SMS criterion", "\n")
paretoFront <- nonDom(observations)
PF_range <- apply(paretoFront, 2, range)
refPoint <- matrix(PF_range[2,] + pmax(1, (PF_range[2,] - PF_range[1,]) * 0.2), 1, nobj)

n.pareto <- nrow(paretoFront)
gain <- -qnorm( 0.5*(0.5^(1/nobj)) )
potSol <- predmean - gain*predsd

nondom <- which(nonDomSet(potSol, paretoFront))
Xnd <- potSol[nondom,,drop=FALSE]

get_hypervolume_reduction <- function(newpoint, paretoFront, refPoint) {
subFront <- paretoFront[which(nonDomSet(paretoFront, matrix(newpoint, nrow=1))),]
potFront <- rbind(subFront, newpoint)
return(dominated_hypervolume(points=t(potFront), ref=refPoint))
}

Xnd_l <- lapply(seq_len(nrow(Xnd)), function(cc) Xnd[cc,])
testnd <- unlist(mclapply(Xnd_l, FUN=get_hypervolume_reduction, paretoFront, refPoint, mc.cores=ncores))

# If mclapply blows up
# testnd <- unlist(lapply(Xnd_l, FUN=get_hypervolume_reduction, paretoFront, refPoint))
test <- rep(0, nrow(integ.pts))
test[nondom] <- testnd
i <- which.max(test)

} else {
if (trace>1) cat("Looking for equilibrium \n")
predmean <- Reduce(cbind, lapply(pred, function(alist) alist\$mean))
if (!is.null(calibcontrol\$target)) {
# Calibration mode
predmean <- (predmean - matrix(rep(calibcontrol\$target, nrow(predmean)), byrow=TRUE, nrow=nrow(predmean)))^2
if (calibcontrol\$log) {
predmean <- log(predmean + calibcontrol\$offset)
}
}
currentEq <- try(getEquilibrium(Z = predmean,  equilibrium = equilibrium, nobj = nobj,
expanded.indices=expanded.indices, n.s=n.s, kweights = NULL,
i <- currentEq\$NE
}

##################################################################
# GET NEW OBSERVATION
##################################################################
xnew <- integ.pts[i,]
while (duplicated(rbind(xnew, model[[1]]@X), fromLast = TRUE)[1]) xnew <- integ.pts[sample.int(nrow(integ.pts), 1),]

## Compute new observation
ynew <- try(fun(xnew, ...))

if (typeof(ynew) == "character" ) {
cat("Unable to compute objective function at iteration ", i, "- optimization stopped \n")
cat("Problem occured for the design: ", xnew, "\n")
cat("Last model and problematic design (xnew) returned \n")
return(list(model=model, integcontrol=integcontrol, xnew=xnew))
}

if (!is.null(noise.var)) {
if (typeof(noise.var) == "closure") {
newnoise.var <- noise.var(xnew)
} else if (typeof(noise.var) == "double") {
newnoise.var <- noise.var
} else {#noise.var ="given_by_fn"
newnoise.var <- ynew[[2]]
ynew <- ynew[[1]]
}
} else {
newnoise.var <- NULL
}
}
####################################################################################################
#### MAIN LOOP ENDS HERE ###########################################################################
####################################################################################################
return(list(model = model, integcontrol=integcontrol))
}
```
vpicheny/GPGame documentation built on Sept. 11, 2019, 1:41 p.m.