https://github.com/cran/GPGame
Tip revision: becef560c88451a1d5de0ef4209f74e7d9114b50 authored by Victor Picheny on 10 June 2017, 05:17:19 UTC
version 1.0.0
version 1.0.0
Tip revision: becef56
solve_game.R
##' Main function to solve games.
##' @title Main solver
##' @details If \code{noise.var="given_by_fn"}, \code{fn} returns a list of two vectors, the first being the objective functions and the second
##' the corresponding noise variances.
##'
##' \code{integcontrol} controls the way the design space is discretized. One can directly provide a set of points \code{integ.pts} with
##' corresponding indices \code{expanded.indices} (for \code{NE}). Otherwise, the points are generated according to the number of strategies \code{n.s}.
##' If \code{n.s} is a scalar, it corresponds to the total number of strategies (to be divided equally among players),
##' otherwise it corresponds to the nb of strategies per player. In addition, one may choose the type of discretization with \code{gridtype}.
##' Options are '\code{lhs}' or '\code{cartesian}'. Finally, \code{lb} and \code{ub} are vectors specifying the bounds for the design variables.
##' By default the design space is \code{[0,1]^d}.
##'
##' \code{simucontrol} controls options on conditional GP simulations. Options are \code{IS}: if \code{TRUE}, importance sampling is used for \code{ynew};
##' \code{n.ynew} number of samples of \code{Y(x_{n+1})} and \code{n.sim} number of sample path generated.
##'
##' \code{plotcontrol} can be used to generate plots during the search. Options are \code{plots} (Boolean, \code{FALSE} by default), \code{compute.actual}
##' (Boolean, \code{FALSE} by default, to draw the actual problem, only for inexpensive \code{fun}), and \code{pbname} (string, for figure title and pdf export).
##'
##' \code{filtercontrol} controls filtering options. \code{filter} sets how to select a subset of simulation and candidate points,
##' either either a single value or a vector of two to use different filters for simulation and candidate points.
##' Possible values are '\code{window}', '\code{Pnash}' (for \code{NE}), '\code{PND}' (probability of non domination), '\code{none}'.
##' \code{nsimPoints} and \code{ncandPoints} set the maximum number of simulation/candidate points wanted
##' (use with filter '\code{Pnash}' for now). Default values are \code{800} and \code{200}, resp.
##' \code{randomFilter} (\code{TRUE} by default) sets whereas the filter acts randomly or deterministically.
##'
##' \code{kmcontrol} Options for handling \code{nobj} \code{\link[DiceKriging]{km}} models.
##' \code{cov.reestim} (Boolean, \code{TRUE} by default) specifies if the kriging hyperparameters
##' should be re-estimated at each iteration,
##'
##' \code{returncontrol} sets options for the last iterations and what is returned by the algorithm.
##' \code{return.Eq} (Boolean, \code{TRUE} by default) specifies
##' if a final search for the equilibrium is performed at the end. \code{finalcrit} sets a different criterion for the last iteration.
##' \code{track.Eq} allows to estimate the equilibrium at each iteration; options are '\code{none}' to do nothing,
##' "\code{mean}" (default) to compute the equilibrium of the prediction mean (all candidates),
##' "\code{empirical}" (for \code{KSE}) and "\code{pex}"/"\code{psim}" (\code{NE} only)
##' for using \code{Pnash} estimate (along with mean estimate, on integ.pts only, NOT reestimated if \code{filter.simu} or \code{crit} is \code{Pnash}).
##'
##' @param fun fonction with vectorial output
##' @param equilibrium either '\code{NE}', '\code{KSE}' or '\code{NKSE}' for Nash/Kalai-Smoridinsky/Nash-Kalai-Smoridinsky equilibria
##' @param crit '\code{sur}' (default) is available for all equilibria, '\code{psim}' and '\code{pex}' are available for Nash
##' @param model list of \code{\link[DiceKriging]{km}} models
##' @param n.init number of points of the initial design of experiments if no model is given
##' @param n.ite number of iterations of sequential optimization
##' @param d variable dimension
##' @param nobj number of objectives (players)
##' @param x.to.obj for \code{NE} and \code{NKSE}, which variables for which objective
##' @param noise.var noise variance. Either a scalar (same noise for all objectives), a vector (constant noise, different for each objective),
##' a function (type closure) with vectorial output (variable noise, different for each objective) or \code{"given_by_fn"}, see Details.
##' If not provided, \code{noise.var} is taken as the average of \code{model@noise.var}.
##' @param integcontrol optional list for handling integration points. See Details.
##' @param simucontrol optional list for handling conditional simulations. See Details.
##' @param plotcontrol optional list for handling during-optimization plots. See Details.
##' @param filtercontrol optional list for handling filters. See Details.
##' @param kmcontrol optional list for handling \code{\link[DiceKriging]{km}} models. See Details.
##' @param returncontrol optional list for choosing return options. See Details.
##' @param ncores number of CPU available (> 1 makes mean parallel \code{TRUE})
##' @param trace controls the level of printing: \code{0} (no printing), \code{1} (minimal printing), \code{3} (detailed printing)
##' @param seed to fix the random variable generator
##' @param ... additional parameter to be passed to \code{fun}
##' @return
##' A list with components:
##' \itemize{
##' \item{\code{model}}{: a list of objects of class \code{\link[DiceKriging]{km}} corresponding to the last kriging models fitted.}
##' \item{\code{Jplus}}{: recorded values of the acquisition function maximizer}
##' \item{\code{integ.pts} and \code{expanded.indices}}{: the discrete space used,}
##' \item{\code{predEq}}{: a list containing the recorded values of the estimated best solution,}
##' \item{\code{Eq.design, Eq.poff}}{: estimated equilibrium and corresponding pay-off (if \code{return.Eq==TRUE})}
##' }
##'
##' @export
##' @importFrom grDevices dev.off pdf rainbow
##' @importFrom graphics axis pairs par points title
##' @importFrom stats pnorm qnorm rnorm dnorm
##' @importFrom methods slot
##' @importFrom graphics filled.contour
##' @import DiceKriging DiceDesign parallel
##' @importFrom MASS mvrnorm
##' @importFrom grDevices terrain.colors
##' @importFrom graphics legend
##' @useDynLib GPGame, .registration = TRUE
##' @references
##' V. Picheny, M. Binois, A. Habbal (2016+), A Bayesian optimization approach to find Nash equilibria,
##' \emph{https://arxiv.org/abs/1611.02440}.
##' @examples
##' \dontrun{
##'
##' ##############################################
##' # Example 1: 2 variables, 2 players, no filter
##' ##############################################
##' # Define objective function (R^2 -> R^2)
##' fun <- function (x)
##' {
##' if (is.null(dim(x))) x <- matrix(x, nrow = 1)
##' b1 <- 15 * x[, 1] - 5
##' b2 <- 15 * x[, 2]
##' return(cbind((b2 - 5.1*(b1/(2*pi))^2 + 5/pi*b1 - 6)^2 + 10*((1 - 1/(8*pi)) * cos(b1) + 1),
##' -sqrt((10.5 - b1)*(b1 + 5.5)*(b2 + 0.5)) - 1/30*(b2 - 5.1*(b1/(2*pi))^2 - 6)^2-
##' 1/3 * ((1 - 1/(8 * pi)) * cos(b1) + 1)))
##' }
##'
##' # To use parallel computation (turn off on Windows)
##' library(parallel)
##' parallel <- FALSE #TRUE #
##' if(parallel) ncores <- detectCores() else ncores <- 1
##'
##' # Simple configuration: no filter, discretization is a 21x21 grid
##'
##' # Grid definition
##' n.s <- rep(21, 2)
##' x.to.obj <- c(1,2)
##' gridtype <- 'cartesian'
##'
##' # Run solver with 6 initial points, 4 iterations
##' # Increase n.ite to at least 10 for better results
##' res <- solve_game(fun, equilibrium = "NE", crit = "sur", n.init=6, n.ite=4,
##' d = 2, nobj=2, x.to.obj = x.to.obj,
##' integcontrol=list(n.s=n.s, gridtype=gridtype),
##' ncores = ncores, trace=1, seed=1)
##'
##' # Get estimated equilibrium and corresponding pay-off
##' NE <- res$Eq.design
##' Poff <- res$Eq.poff
##'
##' # Draw results
##' plotGame(res)
##'
##' ##############################################
##' # Example 2: 4 variables, 2 players, filtering
##' ##############################################
##' fun <- function(x, nobj = 2){
##' if (is.null(dim(x))) x <- matrix(x, 1)
##' y <- matrix(x[, 1:(nobj - 1)], nrow(x))
##' z <- matrix(x[, nobj:ncol(x)], nrow(x))
##' g <- rowSums((z - 0.5)^2)
##' tmp <- t(apply(cos(y * pi/2), 1, cumprod))
##' tmp <- cbind(t(apply(tmp, 1, rev)), 1)
##' tmp2 <- cbind(1, t(apply(sin(y * pi/2), 1, rev)))
##' return(tmp * tmp2 * (1 + g))
##' }
##'
##' # Grid definition: player 1 plays x1 and x2, player 2 x3 and x4
##' # The grid is a lattice made of two LHS designs of different sizes
##' n.s <- c(44, 43)
##' x.to.obj <- c(1,1,2,2)
##' gridtype <- 'lhs'
##'
##' # Set filtercontrol: window filter applied for integration and candidate points
##' # 500 simulation and 200 candidate points are retained.
##' filtercontrol <- list(nsimPoints=500, ncandPoints=200,
##' filter=c("window", "window"))
##'
##' # Set km control: lower bound is specified for the covariance range
##' # Covariance type and model trend are specified
##' kmcontrol <- list(lb=rep(.2,4), model.trend=~1, covtype="matern3_2")
##'
##' # Run solver with 20 initial points, 4 iterations
##' # Increase n.ite to at least 20 for better results
##' res <- solve_game(fun, equilibrium = "NE", crit = "psim", n.init=20, n.ite=2,
##' d = 4, nobj=2, x.to.obj = x.to.obj,
##' integcontrol=list(n.s=n.s, gridtype=gridtype),
##' filtercontrol=filtercontrol,
##' kmcontrol=kmcontrol,
##' ncores = 1, trace=1, seed=1)
##'
##' # Get estimated equilibrium and corresponding pay-off
##' NE <- res$Eq.design
##' Poff <- res$Eq.poff
##'
##' # Draw results
##' plotGame(res)
##'
##' }
solve_game <- function(
fun, ..., equilibrium="NE", crit="sur", model=NULL, n.init=NULL, n.ite, d, nobj, x.to.obj=NULL, noise.var = NULL,
integcontrol=NULL, simucontrol=NULL, plotcontrol=NULL, filtercontrol=NULL, kmcontrol=NULL, returncontrol=NULL,
ncores=1, trace=1, seed=NULL) {
t1 <- Sys.time()
### Initialize variables
set.seed(seed)
if(ncores > 1) parallel <- TRUE
# Check critcontrols
if(crit == 'pex') PnashMethod <- 'exact'
else PnashMethod <- 'simu'
# Check integcontrol
if (!is.null(integcontrol$integ.pts)) integ.pts <- integcontrol$integ.pts else integ.pts <- NULL
if (!is.null(integcontrol$expanded.indices)) expanded.indices <- integcontrol$expanded.indices else expanded.indices <- NULL
if (!is.null(integcontrol$n.s)) n.s <- integcontrol$n.s else n.s <- NULL
if (!is.null(integcontrol$gridtype)) gridtype <- integcontrol$gridtype else gridtype <- "lhs"
if (!is.null(integcontrol$lb)) lb <- integcontrol$lb else lb <- rep(0, d)
if (!is.null(integcontrol$ub)) ub <- integcontrol$ub else ub <- rep(1, d)
# Check simucontrol
if (!is.null(simucontrol$n.ynew)) n.ynew <- simucontrol$n.ynew else n.ynew <- 10
if (!is.null(simucontrol$n.sim)) n.sim <- simucontrol$n.sim else n.sim <- 10
if (!is.null(simucontrol$IS)) IS <- simucontrol$IS else IS <- TRUE
cross <- FALSE
# Check plotcontrol
if (!is.null(plotcontrol$compute.actual)) compute.actual <- plotcontrol$compute.actual else compute.actual <- FALSE
if (!is.null(plotcontrol$plots)) plots <- plotcontrol$plots else plots <- FALSE
if (!is.null(plotcontrol$pbname)) pbname <- plotcontrol$pbname else pbname <- NULL
if (!is.null(pbname)) exportPDF <- TRUE else exportPDF <- FALSE
# Check filtercontrol
if (!is.null(filtercontrol$filter)) filter <- filtercontrol$filter else filter <- 'none'
if (!is.null(filtercontrol$nsimPoints)) nsimPoints <- filtercontrol$nsimPoints else nsimPoints <- 800
if (!is.null(filtercontrol$ncandPoints)) ncandPoints <- filtercontrol$ncandPoints else ncandPoints <- 200
if (!is.null(filtercontrol$randomFilter)) randomFilter <- filtercontrol$randomFilter else randomFilter <- TRUE
if (length(filter)==1) {
simu.filter <- cand.filter <- filter
} else {
simu.filter <- filter[1]
cand.filter <- filter[2]
}
if (equilibrium=="KSE") {
if ("Pnash" %in% c(simu.filter, cand.filter)) cat("Pnash filter only available for NE; switching to PND \n")
if (simu.filter == "Pnash") simu.filter <- 'PND'
if (cand.filter == "Pnash") cand.filter <- 'PND'
}
if (length(randomFilter)==1) randomFilter <- rep(randomFilter, 2)
# Check kmcontrol
if (!is.null(kmcontrol$cov.reestim)) plots <- kmcontrol$cov.reestim else cov.reestim <- TRUE
if (!is.null(kmcontrol$model.trend)) model.trend <- kmcontrol$model.trend else model.trend <- ~1
if (!is.null(kmcontrol$lb)) kmlb <- kmcontrol$lb else kmlb <- rep(.1,d)
if (!is.null(kmcontrol$ub)) kmub <- kmcontrol$ub else kmub <- rep(1,d)
if (!is.null(kmcontrol$nugget)) kmnugget <- kmcontrol$nugget else kmnugget <- 1e-8
if (!is.null(kmcontrol$control)) control <- kmcontrol$control else control <- list(trace=FALSE)
if (!is.null(kmcontrol$covtype)) covtype <- kmcontrol$covtype else covtype <- "matern5_2"
# Check returncontrol
if (!is.null(returncontrol$return.Eq)) return.Eq <- returncontrol$return.Eq else return.Eq <- TRUE
if (!is.null(returncontrol$finalcrit)) finalcrit <- returncontrol$finalcrit else finalcrit <- crit
if (!is.null(returncontrol$track.Eq)) track.Eq <- returncontrol$track.Eq else track.Eq <- 'mean'
if (equilibrium=="KSE") {
if (!is.null(track.Eq)) {
if (track.Eq %in% c("psim", "pex")) {
cat("track.Eq must be set to none or mean for KSE - switching to mean \n")
track.Eq <- "mean"
}
}
finalcrit <- "sur"
} else {
if (track.Eq == "empirical") {
cat("track.Eq= empirical option only available for KSE; switching to psim \n")
track.Eq <- "psim"
}
}
# 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
}
### Initialize variables
window <- NULL
all.Jplus <- c()
include.obs <- FALSE
Eq.design.estimate <- Eq.poff.estimate <- trackPnash <- NULL
predEq <- list()
if (crit!="sur" && equilibrium!="NE"){
cat("pex and psim available only for Nash equilibria; crit switched to sur \n")
crit <- "sur"
}
if (trace>0) cat("--------------------------\n Starting", equilibrium, "search with:", "\n",
crit, "strategy,","\n",
simu.filter, "simulation point filter,", "\n",
cand.filter, "candidate point filter,", "\n",
"among (", n.s, ") strategies \n --------------------------\n " )
#### Initial design and models ####################
if (is.null(model)){
if (is.null(n.init)) n.init <- 5*d
design <- lhsDesign(n.init, d, seed = 42)$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 {#noise.var ="given_by_fn"
# newnoise.var <- response[[2]]
# ynew <- response[[1]]
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)) {
res <- generate_integ_pts(n.s=n.s, d=d, nobj=nobj, x.to.obj=x.to.obj, gridtype=gridtype, lb = lb, ub = ub)
integcontrol$integ.pts <- integ.pts <- res$integ.pts
integcontrol$expanded.indices <- expanded.indices <- res$expanded.indices
}
n.integ.pts <- nrow(integ.pts)
cand.pts <- integ.pts
if (is.null(expanded.indices)) {
sorted <- FALSE
} else {
sorted <- !is.unsorted(expanded.indices[,nobj])
}
if (compute.actual) {
if (plots) {
if (nobj==2) {
if(exportPDF) pdf(paste0(pbname, "_obj.pdf"), width=5, height=5)
par(mfrow=c(1,1), mar=c(3,3,2,1), mgp=c(2,1,0))
groundTruth <- plotGameGrid(fun = fun, domain = matrix(rep(c(0,1), each = d), d), graphs = "objective",
n.grid = 31, x.to.obj = x.to.obj, integcontrol=integcontrol, equilibrium = equilibrium, ... = ...)
if(exportPDF) dev.off()
}
if(exportPDF) pdf(paste0(pbname, "_design.pdf"), width=5, height=5)
par(mfrow=c(1,1), mar=c(3,3,2,1), mgp=c(2,1,0))
plotGameGrid(fun = fun, domain = matrix(rep(c(0,1), each = d), d), graphs = "design", n.grid = 31, x.to.obj = x.to.obj,
integcontrol=my.integcontrol, type = equilibrium, ... = ...)
if(exportPDF) dev.off()
}
}
if (trace>2) cat("Time for initialization: \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
#############################################################
#### MAIN LOOP STARTS HERE ##################################
#############################################################
if (plots) par(mfrow=c(min(n.ite+1, 3),2), mar=c(2,2,2,2))
Jnres <- rep(NA, n.ite + 1)
for (ii in 1:(n.ite+1)){
if (ii < (n.ite+1) && trace>0){cat("--- Iteration #", ii, " ---\n")}
t0 <- t1 <- Sys.time()
##################################################################
# MODELS UPDATE
##################################################################
if (ii > 1) {
if (ii == n.ite){
crit <- finalcrit
cand.filter <- 'none'
}
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")
# if (ii > 1) allX.new <- rbind(model[[u]]@X[(model[[u]]@n + 1):(model[[u]]@n+ii-1),, drop=FALSE], X.new)
return(list(
par = X.new,
values = ynew,
nsteps = ii,
model = model,
Jplus = Jplus, integcontrol=integcontrol))
} else {
model[[u]] <- newmodel[[u]]
}
}
}
if (trace>2) cat("Time for models updates: \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
observations <- Reduce(cbind, lapply(model, slot, "y"))
#--------------------------------------------------------#
# Regular loop : find infill point
#--------------------------------------------------------#
if (ii < (n.ite+1) || return.Eq){
if (ii==(n.ite+1)) include.obs <- TRUE
pred <- mclapply(model, FUN=predict, newdata = integ.pts, checkNames = FALSE, type = "UK", light.return = TRUE, mc.cores=ncores)
if (trace>2) cat("Time for predictions: \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
##################################################################
# FILTERING INTEGRATION POINTS
##################################################################
if (simu.filter == "window") {
if (is.null(window) || crit != "sur") {
# Special case for initialization - otherwise do not change the old "window" value
predmean <- Reduce(cbind, lapply(pred, function(alist) alist$mean))
Eq_simu <- getEquilibrium(Z = predmean, equilibrium = equilibrium, nobj = nobj, expanded.indices=expanded.indices,
n.s=n.s, sorted = sorted, cross = cross)
if(nrow(Eq_simu) < 1 || all(is.na(Eq_simu))){
## if no window has been defined in previous iterations, use the mean, otherwise keep the old one
if(is.null(window)){
window <- colMeans(observations)
}
} else if(nrow(Eq_simu) == 1) {
window <- as.vector(Eq_simu)
} else {
window <- colMeans(Eq_simu, na.rm = TRUE)
}
}
options <- list(window = window)
} else if (simu.filter == "Pnash") {
options <- list(method=PnashMethod, nsim=100)
}
if (trace > 1 && simu.filter == "window") cat("window for integ.pts", window, "\n")
# Reduce the number of integration points first
if (simu.filter != 'none') {
filt <- filter_for_Game(n.s.target = pmin(n.s, round(nsimPoints^(1/nobj))), integcontrol=integcontrol,
model = model, predictions=pred, type=simu.filter, options = options, random=randomFilter[1],
include.obs=include.obs, ncores = ncores, equilibrium=equilibrium)
I <- filt$I
} else {
I <- 1:nrow(integ.pts)
}
my.integ.pts <- integ.pts[I,]
my.expanded.indices <- expanded.indices[I,]
my.pred <- lapply(pred, function(alist) list(mean=alist$mean[I], sd=alist$sd[I]))
my.n.s <- apply(my.expanded.indices, 2, function(x) length(unique(x)))
if (trace>1) cat("Number of strategies after simu filtering: (", my.n.s, ")\n")
my.expanded.indices2 <- my.expanded.indices
for (i in 1:ncol(my.expanded.indices)) {
unik_i <- unique(my.expanded.indices[,i])
for (j in 1:length(unik_i)) {
my.expanded.indices2[which(my.expanded.indices[,i]==unik_i[j]),i] <- j
}
}
my.expanded.indices <- my.expanded.indices2
my.integcontrol <- list(expanded.indices=my.expanded.indices, integ.pts=my.integ.pts, n.s=my.n.s)
if (trace>2) cat("Time for filter 1 (simu): \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
##################################################################
# Precalculations and simulations (if not performed already)
##################################################################
precalc.data <- mclapply(model, FUN=precomputeUpdateData, integration.points=my.integ.pts, mc.cores=ncores)
if(crit == 'sur' || cand.filter=="window"){
my.Simu <- t(Reduce(rbind, mclapply(model, simulate, nsim=n.sim, newdata=my.integ.pts, cond=TRUE,
checkNames=FALSE, nugget.sim = 10^-8, mc.cores=ncores)))
}
if (trace>2) cat("Time for precalc + simu: \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
##################################################################
# FILTERING CANDIDATE POINTS
##################################################################
cand.pts <- my.integ.pts
if ((simu.filter == "window" && crit == 'sur') || cand.filter == "window" || plots) {
Eq_simu <- getEquilibrium(my.Simu, equilibrium = equilibrium, nobj = nobj, expanded.indices=my.expanded.indices,
n.s=my.n.s, sorted = sorted, cross = cross)
Eq_simu <- Eq_simu[which(!is.na(Eq_simu[,1])),, drop = FALSE]
Jnres[ii] <- det(cov(Eq_simu))
temp <- apply(Eq_simu, 2, range, na.rm = TRUE)
if (!any(is.na(temp))) window <- temp
options <- list(window = window)
} else if (cand.filter == "Pnash") {
options <- list(method=PnashMethod, nsim = 100)
}
# Plots
if (nobj==2 && plots) {
if(exportPDF) pdf(paste0(pbname, "_objspace_ite_", ii, ".pdf"), width=5, height=5)
par(mfrow=c(1,1), mar=c(3,3,2,1), mgp=c(2,1,0))
if(compute.actual){
xrange <- range(groundTruth$response.grid[,1])
yrange <- range(groundTruth$response.grid[,2])
xrange[1] <- xrange[1] - .1*(xrange[2] - xrange[1])
yrange[1] <- yrange[1] - .1*(yrange[2] - yrange[1])
plot(NA, xlim=xrange, ylim=yrange,
main=paste0("Iteration ", ii), xlab=expression(f[1]), ylab=expression(f[2]))
points(groundTruth$trueEqPoff[1], groundTruth$trueEqPoff[2], pch = 24, bg = "green", cex=2)
# plotParetoEmp(trueParetoFront , col="green")
points(observations[,1], observations[,2], col="blue", pch=19, cex=1.2)
# plotParetoEmp(t(nondominated_points(t(observations))), col="blue")
points(Eq_simu[,1], Eq_simu[,2], col="red", pch=3, lwd=2)
legend("topright", legend = c("Actual Eq", "Simulated Eq", "Current obs", "New obs"), lwd = rep(NA,4),
col = c("black", "red", "blue", "red"), pch=c(24, 3, 19, 19), pt.cex= c(2, 1.2, 1, 2),
pt.lwd=c(1,2,1,1), pt.bg=c("green", NA, NA, NA))
}else{
UQ_eq <- TRUE
if(crit == "sur") UQ_eq <- FALSE
plotGame(res = list(model = model, Jplus = all.Jplus, integcontrol=integcontrol,
predEq = predEq, trackPnash = trackPnash),
equilibrium = equilibrium, add = FALSE, UQ_eq = UQ_eq, simucontrol = simucontrol)
title(paste0("Iteration ", ii))
if(crit == "sur") points(Eq_simu[,1], Eq_simu[,2], col="red", pch=3, lwd=2)
}
if (ii > n.ite) if(exportPDF) dev.off()
}
if (trace>1 && cand.filter == "window") cat("window for cand.pts", window, "\n")
if (cand.filter != 'none') {
if (cand.filter == "Pnash" && simu.filter == "Pnash") {
# Special case if Pnash used twice: do not redo calculations
if (randomFilter[2]) {
J <- sample.int(length(I), size=min(length(I), ncandPoints), replace=FALSE, prob=filt$crit[I])
} else {
J <- I[1:min(length(I), ncandPoints)]
}
} else {
# Regular case
filt2 <- filter_for_Game(n.s.target = pmin(my.n.s, round(ncandPoints^(1/nobj))), integcontrol=my.integcontrol,
model = model, predictions=my.pred, type=cand.filter, options = options,
random=randomFilter[2], include.obs=include.obs, ncores = ncores, equilibrium=equilibrium)
J <- filt2$I
}
} else {
# No filter case
J <- 1:nrow(cand.pts)
}
## Remove already evaluated designs from candidates
if((ii > 1) && (ii<(n.ite+1))){
my.pred.s <- my.pred[[1]]$sd[J]
J <- J[which(my.pred.s >= sqrt(model[[1]]@covariance@sd2)/1e6)]
}
my.cand.pts <- cand.pts[J,, drop=FALSE]
cand.s <- my.expanded.indices[J,, drop=FALSE]
# my.pred.s <- my.pred[[1]]$sd[J]
if (trace>2) cat("Time for filter 2 (cand): \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
if (trace>0) cat("Nb of integration / candidate pts:", length(I), "/", length(J), "\n")
##################################################################
# CRITERION MINIMIZATION
##################################################################
my.integ.pts <- data.frame(my.integ.pts)
if (crit=="sur") {
Jplus <- unlist(mclapply(X=J, FUN=crit_SUR_Eq, model=model, integcontrol=my.integcontrol,
Simu=my.Simu, equilibrium = equilibrium, precalc.data=precalc.data,
n.ynew=n.ynew, IS=IS, cross=cross, mc.cores = ncores))
} else if (crit == 'pex') {
Jplus <- -crit_PNash(idx=J, integcontrol=my.integcontrol,
type = 'exact', model = model, ncores = ncores)
} else if (crit=="psim") {
Jplus <- -crit_PNash(idx=J, integcontrol=my.integcontrol,
type = 'simu', model = model, ncores = ncores, control = list(nsim = 100))
} else {
cat("wrong crit \n")
break;
}
if(all(is.na(Jplus))){
cat("No equilibrium for this problem \n")
return(list(model = model, Jplus = Jplus, integcontrol=integcontrol))
}
Jplus[is.na(Jplus)] <- max(Jplus, na.rm = TRUE)
if (trace>2) cat("Time for optim: \n")
if (trace>2) print(Sys.time() - t1)
if (trace > 0 && crit=="sur") cat("Jplus:", min(Jplus), '\n')
if (trace > 0 && crit!="sur") cat("Pmax:", -min(Jplus), '\n')
if (ii < n.ite+1) {
# Get best solution
i <- which.min(Jplus)
xnew <- my.cand.pts[i,]
ynew <- fun(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
}
all.Jplus <- c(Jplus, min(Jplus))
if (trace>1) cat("New design evaluated: \n")
if (trace>1) cat(xnew, "\n")
if (trace>1) cat("Corresponding new observation: \n")
if (trace>1) cat(ynew, "\n")
if (trace>0) cat("Total iteration time: \n")
if (trace>0) print(Sys.time() - t0)
##################################################################
# PLOTS
##################################################################
if (nobj==2 && plots) {
# Representation du nouveau point
points(ynew[1], ynew[2], col="red", pch=19, cex=2)
if(exportPDF) dev.off()
}
if (d==2 && plots) {
# Representation du critere
if (gridtype =="cartesian") {
Jgrid <- rep(max(Jplus), prod(n.s))
Jgrid[J] <- Jplus
if(exportPDF) pdf(paste0(pbname, "_crit_ite_", ii, ".pdf"), width=5, height=5)
par(mfrow=c(1,1), mar=c(3,3,2,1), mgp=c(2,1,0))
# TO CHECK: that integ.pts corresponds to the structure needed in filled.contour
filled.contour(seq(0, 1, length.out = n.s[1]), seq(0, 1, length.out = n.s[2]), nlevels = 20,
matrix(Jgrid, n.s[1], n.s[2]), main = paste0("Iteration ", ii),
xlab = expression(x[1]), ylab = expression(x[2]),
plot.axes = {axis(1); axis(2);
points(xnew[1], xnew[2], pch = 21, bg = "blue", cex=2);
points(model[[1]]@X[,1], model[[1]]@X[,2], pch = 21, bg = "white");
}
)
if(compute.actual) points(groundTruth$trueEqdesign[1], groundTruth$trueEqdesign[2], pch = 24, bg = "green", cex=1.5)
if(exportPDF) dev.off()
}
}
}
# else {
# #--------------------------------------------------------#
# # If optimization is over: find best Eq estimate
# #--------------------------------------------------------#
# Eq.design.estimate <- cand.pts
# Eq.poff.estimate <- -Jplus
# }
##################################################################
# Tracking of estimated best solution
##################################################################
if (is.null(track.Eq)) track.Eq <- "psim"
if (track.Eq != 'none' || (ii == n.ite+1 && return.Eq)){
if (track.Eq == "empirical") {
observations <- Reduce(cbind, lapply(model, slot, "y"))
currentEq <- getEquilibrium(Z=observations, equilibrium = equilibrium, nobj=nobj, return.design=TRUE, cross=cross)
# if (ii == 1) {
# predEq <- list(Eq_emp)
# } else {
# predEq <- c(predEq, list(Eq_emp))
# }
} else if (track.Eq == "psim" || track.Eq == "pex") {
if (ii == 1) trackPnash <- list()
if (simu.filter == "Pnash") {
ImaxPnash <- which.max(filt$crit)
maxPnash <- max(filt$crit)
# trackPnash <- c(trackPnash, list(list(c(predEqPoff = unlist(lapply(my.pred, function(x) x$mean[ImaxPnash])),
# Eq = my.integ.pts[ImaxPnash,, drop = F]), Pnash = maxPnash)))
currentEq <- list(c(predEqPoff = unlist(lapply(my.pred, function(x) x$mean[ImaxPnash])),
Eq = my.integ.pts[ImaxPnash,, drop = F]), Pnash = maxPnash)
} else {
if (track.Eq == 'pex') {
estPnash <- crit_PNash(idx = 1:nrow(my.integcontrol$expanded.indices), integcontrol=my.integcontrol,
type = 'exact', model = model, ncores = ncores)
} else {
estPnash <- crit_PNash(idx = 1:nrow(my.integcontrol$expanded.indices), integcontrol=my.integcontrol,
type = 'simu', model = model, ncores = ncores,
control = list(nsim = 100))
}
ImaxPnash <- which.max(estPnash)
maxPnash <- max(estPnash)
currentEq <- list(predEqPoff = unlist(lapply(my.pred, function(x) x$mean[ImaxPnash])),
Eq = my.integ.pts[ImaxPnash,, drop = F], Pnash = maxPnash)
# trackPnash <- c(trackPnash, list(list(predEqPoff = unlist(lapply(my.pred, function(x) x$mean[ImaxPnash])),
# Eq = my.integ.pts[ImaxPnash,, drop = F], Pnash = maxPnash)))
}
if (trace>2) cat("Maximum estimated Pnash:", maxPnash," \n")
if (trace>2) print(Sys.time() - t1)
} else {
## Eq of the GP predictive means
predmean <- Reduce(cbind, lapply(pred, function(alist) alist$mean))
currentEq <- getEquilibrium(Z = predmean, equilibrium = equilibrium, nobj = nobj,
expanded.indices=expanded.indices, n.s=n.s,
sorted = sorted, cross = cross, return.design = T)
}
predEq <- c(predEq, list(currentEq))
}
if (trace>2) cat("Time for estimating best candidate: \n")
if (trace>2) print(Sys.time() - t1)
t1 <- Sys.time()
}
}
if (return.Eq) {
#--------------------------------------------------------#
# When optimization is over: find best Eq estimate
#--------------------------------------------------------#
Eq.design.estimate <- integ.pts[predEq[[length(predEq)]]$NE,]
Eq.poff.estimate <- predEq[[length(predEq)]]$NEPoff
return(list(model = model, Eq.design=Eq.design.estimate, Eq.poff=Eq.poff.estimate,
Jplus = all.Jplus, integcontrol=integcontrol,
predEq = predEq))
} else {
return(list(model = model, Jplus = all.Jplus, integcontrol=integcontrol,
predEq = predEq))
}
}