R/integration_design_robinv.R
In IRSN/RobustInv: "Robust inversion of expensive black-box functions"
Defines functions
integration_design_robinv
Documented in
integration_design_robinv
#' @title Constructing a set of integration and simulation points in robust inversion
#'
#' @description Constructing a set of integration and simulation points in robust inversion
#'
#' @param integcontrol An important list with the following possible fields.
#' \code{distrib}: This field needs to be set to either \code{"sur"} or \code{"surnew"}
#' depending on the sampling criterion that is used.
#'
#' (I) If \code{integcontrol$distrib = "surnew"} (default) the
#' user can set the field \code{n.points} (default: 20) which is the number of integration points. The points
#' live in the space of controlled parameters, called Dinv here.
#' The distribution of these n.points points
#' is controlled through the \code{finaldistrib} field. Possible values are
#' \code{"spec"} to set up manually the integration points through the \code{inv.integration.points}
#' field, \code{"MC"} for
#' uniformly distributed points over Dinv, \code{"sobol"} (default) to use the sobol sequence and \code{"surnew"}
#' to use importance sampling based on a density proportional to pn (1-pn) with pn being the
#' excursion probability of the integration points. In that last case, extra fields are required:
#' \code{n.candidates} is the number of points which are candidate to become integration points (default: n.points*3).
#' The important sampling scheme will choose (with replacement) n.points integration points
#' among the candidates. The distribution of the candidates is set with the field
#' \code{init.distrib} which can be either \code{"MC"} for uniform points,
#' \code{"sobol"} (default) to use the sobol
#' sequence or \code{"spec"} to set up manually the matrix of candidates
#' with the \code{inv.integration.points} field.
#' When scaling is used in the kriging model, the candidate points can be unscaled by setting the
#' \code{unscale.inv.integration.points}
#' to \code{TRUE} (default: \code{FALSE}).
#' Regarding the nuisance parameters, the number of simulation points used, e.g., to compute
#' pn is set with the field \code{n.optpoints} (default: 100). There is the option to choose these points
#' among candidates by setting the field \code{choose_optpoints} to \code{TRUE}
#' (default: \code{FALSE}).
#' In that case, the number of candidates is set with the field \code{n.optpoints.candidates}
#' (default: 10*n.optpoints).
#' The field \code{opt.simulation.points.distrib} refers to the distribution of the
#' simulation points if \code{integcontrol$choose_optpoints=FALSE} and to the
#' candidates otherwise. Possible values are \code{"spec"} for points set up manually with the
#' \code{opt.simulation.points} field and \code{"sobol"} (default) for the sobol sequence.
#' When the sobol sequence is used, there is the option to add the vertices of the input
#' domain (of the nuisance parameters) with the field \code{addvertices} (default: \code{TRUE}).
#' When scaling is used in the kriging model, the simulation points (or the candidates)
#' can be unscaled by setting the \code{unscale.opt.simulation.points} field
#' to \code{TRUE} (default: \code{TRUE}).
#' The number of simulation of Gaussian process sample paths per integration points
#' is controlled with the field
#' \code{nsimu} (default: 1000).
#' Finally, when multiple calls are performed, the \code{seed} argument can be overwritten
#' by setting the \code{seed} field.
#'
#' (II) The option \code{integcontrol$distrib = "sur"} needs to be used only when the
#' sur criterion (and not surnew) is used. In that case, the fields \code{n.points} and
#' \code{finaldistrib} work similarly except that the importance sampling option is set
#' with \code{integcontrol$finaldistrib="sur"} and not \code{"surnew"}. In that case,
#' the field \code{init.distrib} works as in (I). Since there is no Gaussian process simulation here,
#' the excursion probability computation is set with the field \code{pnmethod}.
#' Possible values are \code{"fixed"} and \code{"optimMC"} (recommended). The number
#' of MC trials to get pn (through a multivariate normal cdf call) is set with the
#' field \code{pop.size}. Finally, the discretization in the space of nuisance parameters
#' is also obtained with the field \code{n.optpoints} (maximum value: 10).
#' @param d Dimension (i.e. number of scalar input variables) of the objective function.
#' @param lower Array of size d. Lower bound of the input domain.
#' @param upper Array of size d. Upper bound of the input domain.
#' @param opt.index Array with integers corresponding to the indices of the nuisance parameters.
#' @param inv.index Array with integers corresponding to the indices of the controlled parameters.
#' We must have \code{d = length(opt.index) + length(inv.index)}.
#' @param model The current kriging model. km object.
#' @param T Target threshold.
#' @param min.prob Minimum probability for a candidate integration point to be picked.
#' @param seed The random seed. For repeatability.
#' @return Depeding on the fields in the \code{integcontrol} list, the return is different.
#' See the description on the integcontrol argument and the details.
#' @details (I) If the \code{"surnew"} criterion is used, then the output is a list with
#' the following fields:
#' \cr
#' (i) inv.integration.points: n x dinv matrix containing the integration points.
#' n is the number of integration points and dinv the dimension of the space of controlled parameters.
#' \cr
#' (ii) integration.weights: array of size n containing the weights given to each integration points.
#' \cr
#' (iii) pn: array of size n containing the current excursion probabilities for each integration points. These
#' excursion probabilities are based on nsimu conditional Gaussian process simulations.
#' \cr
#' (iv) allsimu: p x (nsimu*n) matrix where p is the number of simulation points.
#' The first nsimu columns are the nsimu simulations in p points for the first integration point,
#' the next nsimu columns are linked to the next integration point and so on.
#' \cr
#' (v) allsimucentered: p x (nsimu*n) matrix containing the later simulations centered by
#' substracting the kriging mean of the points where the simulations are performed.
#' \cr
#' (vi) allsimupoints: (n*p) x d matrix containing all the points where simulations are performed.
#' The first p rows correspond to the simulation points linked to the first integration points, and so on.
#' \cr
#' (vii) allprecomp: a list with 3 fields obtained from a call to the function
#' \code{precomputeUpdateData} on the points \code{allsimupoints}. Help about \code{precomputeUpdateData}
#' is given in the KrigInv package.
#' \cr
#' (viii) allKn.inv: (n*p) x p matrix containing n different p x p matrices. Matrix number i is the inverse
#' of the p x p non-conditional covariance matrix of the p simulation points associated to the
#' integration point i.
#' \cr
#' (ix) allmn: n*p array with the kriging means of all n*p simulation points stored in \code{allsimupoints}.
#' \cr
#' (II) If the \code{"sur"} criterion is used, then the output is still a list, which however has different fields
#' than in (I):
#' \cr
#' (i) integration.points : (p*n) x d matrix containing all the simulation points. The first p points are
#' associated to the first integration point and so on.
#' \cr
#' (ii) cov: see description of the field allKn.inv (viii) in (I).
#' \cr
#' (iii) mean: see description of the field allmn (ix) in (I).
#' \cr
#' (iv) pn: array of size n containing the current excursion probabilities for each integration points.
#' @export
#' @author Clement Chevalier \email{clement.chevalier@@unine.ch}
#' @examples
#' library(KrigInv)
#' library(randtoolbox)
#' myfun <- branin_robinv
#' d <- 3
#'
#' set.seed(8)
#'
#' n0 <- 30
#' T <- 10
#' opt.index <- c(3)
#' inv.index <- c(1,2)
#' lower <- rep(0,times=d)
#' upper <- rep(1,times=d)
#'
#' design <- matrix(runif(d*n0),nrow=n0)
#' response <- myfun(design)
#' model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2")
#'
#' integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
#' n.candidates=50,nsimu=1000,n.optpoints = 50,
#' choose_optpoints=TRUE,n.optpoints.candidates=500)
#' \dontrun{
#' obj <- integration_design_robinv(integcontrol = integcontrol,d=d,lower=lower,upper=upper,
#' opt.index=opt.index,inv.index=inv.index,model=model,T=T)
#' }
#' #####################################
#' # An example with scaling
#' library(KrigInv)
#' myfun <- function(x){ return(-1*branin_robinv(x) - 50*sin(min(100,1/x[3])) ) }
#' d <- 3
#'
#' set.seed(8)
#'
#' n0 <- 60
#' T <- 30
#' opt.index <- c(3)
#' inv.index <- c(1,2)
#' lower <- rep(0,times=d)
#' upper <- rep(1,times=d)
#'
#' design <- matrix(runif(d*n0),nrow=n0)
#' response <- apply(X = design,FUN = myfun,MARGIN = 1)
#' knots.number <- c(0,0,3)
#' knots <- generate_knots(knots.number = knots.number , d = d)
#'
#' model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2",scaling = TRUE,knots=knots)
#'
#' integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
#' n.candidates=200,nsimu=100,n.optpoints = 50,
#' choose_optpoints=TRUE,n.optpoints.candidates=500,
#' unscale.opt.simulation.points=TRUE)
#' \dontrun{
#' obj <- integration_design_robinv(integcontrol = integcontrol,d=d,lower=lower,upper=upper,
#' opt.index=opt.index,inv.index=inv.index,model=model,T=T)
#' }
integration_design_robinv <- function(integcontrol=NULL,d=NULL,lower,upper,
opt.index,inv.index,
model,T,
min.prob=0.001,seed=NULL){
result<-NULL
if(!is.null(seed)) set.seed(seed)
if(is.null(d)) d <- length(lower)
if (length(lower) != length(upper) ) print("integration_design_robinv: 'lower' and 'upper' must have the same length")
if(is.null(integcontrol$distrib)) integcontrol$distrib <- "surnew"
if(integcontrol$distrib=="sur"){
# only for the "sur" criterion.
n.optpoints <- integcontrol$n.optpoints
if(is.null(n.optpoints)) n.optpoints <- 10
if(n.optpoints > 10) n.optpoints <- 10 # maximum value because we cannot go in dimension higher than 20
finaldistrib <- integcontrol$finaldistrib
if(is.null(finaldistrib)) finaldistrib <- "sobol" # other possible value is "sur"
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]
d.inv <- length(inv.index);d.opt <- length(opt.index)
n.points <- integcontrol$n.points
if(is.null(n.points)) n.points <- 20
if(finaldistrib=="MC") inv.integration.points <- t(lower.inv + t(matrix(runif(d.inv*n.points),ncol=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="sobol") inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="spec") inv.integration.points <- integcontrol$inv.integration.points
# We now deal with the importance sampling case
if(finaldistrib=="sur"){
# construction of the candidates
n.candidates <- integcontrol$n.candidates
if(is.null(n.candidates)) n.candidates <- n.points*5
init.distrib <- integcontrol$init.distrib
if(is.null(init.distrib)) init.distrib <- "sobol"
if(init.distrib=="sobol") initial.inv.integration.points <- t(lower.inv+t(sobol(n=n.candidates,dim=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="MC") initial.inv.integration.points <- t(lower.inv+t(matrix(runif(d.inv*n.candidates),ncol=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="spec") initial.inv.integration.points <- integcontrol$inv.integration.points
inv.integration.points <- initial.inv.integration.points
}
if(d.inv==1) inv.integration.points <- matrix(inv.integration.points,ncol=1)
pnmethod <- integcontrol$pnmethod
if(is.null(pnmethod)) pnmethod <- "fixed" # other possible value is "optimMC"
if(pnmethod=="fixed"){
# in Dopt, some points are generated (and fixed) using the sobol sequence (option 1)
# or uniformly (option 2)
# This option is computationally cheap; but is expected to give excursion probabilities
# that might be widely underestimated.
N <- nrow(inv.integration.points)
integration.points <- matrix(c(0),nrow=N*n.optpoints,ncol=d)
if(is.null(integcontrol$pnmethod.fixed.option)) integcontrol$pnmethod.fixed.option <- 1
if(integcontrol$pnmethod.fixed.option==1){
opt.integration.points <- t(lower.opt + t(sobol(n = n.optpoints, dim = d.opt))*(upper.opt - lower.opt))
if(d.opt==1) opt.integration.points <- matrix(opt.integration.points,ncol=1)
indices <- expand.grid(c(1:n.optpoints),c(1:N))
index1 <- indices[,1]
index2 <- indices[,2]
integration.points[,opt.index] <- as.matrix(opt.integration.points[index1,])
integration.points[,inv.index] <- as.matrix(inv.integration.points[index2,])
}else{
opt.integration.points <- t(lower.opt + t(matrix(runif(n.optpoints*N*d.opt),ncol=d.opt))*(upper.opt - lower.opt))
if(d.opt==1) opt.integration.points <- matrix(opt.integration.points,ncol=1)
indices <- expand.grid(c(1:n.optpoints),c(1:N))
indices[,1] <- c(1:(n.optpoints*N))
index1 <- indices[,1]
index2 <- indices[,2]
integration.points[,opt.index] <- as.matrix(opt.integration.points[index1,])
integration.points[,inv.index] <- as.matrix(inv.integration.points[index2,])
}
#calculation of all cov matrices, means and pn
res <- optim_pn_MC(inv.integration.points=inv.integration.points,model=model,T=T,
opt.index=opt.index,inv.index=inv.index,lower=lower,upper=upper,
n.optpoints=n.optpoints,pop.size=1,
MC=FALSE,deterministicMat=integration.points)
}
if(pnmethod=="optimMC"){
pop.size <- integcontrol$pop.size
if(is.null(pop.size)) pop.size <- 50
res <- optim_pn_MC(inv.integration.points=inv.integration.points,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,
n.optpoints=n.optpoints,pop.size=pop.size,
MC=TRUE)
}
if(finaldistrib!="sur") return(res)
if(finaldistrib=="sur"){
# some importance sampling
pn <- res$pn;Tau.n <- pn*(1-pn);Tau.n.sum <- sum(Tau.n)
if(Tau.n.sum==0) Tau.n.sum <- 1
prob.n <- pmax(Tau.n/Tau.n.sum,min.prob/n.candidates);prob.n <- prob.n/sum(prob.n)
weight.n <- 1/(prob.n*n.candidates*n.points)
prob.n.copy <-c(0,prob.n);prob.n.cum <-cumsum(prob.n.copy)
my.indices <- findInterval(runif(integcontrol$n.points),prob.n.cum,all.inside=TRUE)
pn.my.indices <- pn[my.indices]
integration.weights <- weight.n[my.indices]
new.mean <- rep(0,times=n.optpoints*n.points)
old.mean <- res$mean
new.cov <- matrix(c(0),nrow=n.optpoints*n.points,ncol=n.optpoints)
old.cov <- res$cov
new.intpoints <- matrix(c(0),nrow=n.optpoints*n.points,ncol=d)
old.intpoints <- res$integration.points
for(i in 1:n.points){
index.i <- my.indices[i]
new.mean[(1+(i-1)*n.optpoints):(i*n.optpoints)] <- old.mean[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints)]
new.cov[(1+(i-1)*n.optpoints):(i*n.optpoints),] <- old.cov[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints),]
new.intpoints[(1+(i-1)*n.optpoints):(i*n.optpoints),] <- old.intpoints[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints),]
}
result <- list(integration.points=new.intpoints,
integration.weights=integration.weights,
mean=new.mean,
cov=new.cov,
pn=pn.my.indices)
return(result)
}
}
if(integcontrol$distrib=="surnew"){
## only for the "surnew" criterion.
n.optpoints <- integcontrol$n.optpoints
choose_optpoints <- integcontrol$choose_optpoints
n.optpoints.candidates <- integcontrol$n.optpoints.candidates
nsimu <- integcontrol$nsimu
unscale.inv.integration.points <- integcontrol$unscale.inv.integration.points
if(is.null(nsimu)) nsimu <- 1000
if(is.null(n.optpoints)) n.optpoints <- 100
if(n.optpoints > 1000) n.optpoints <- 1000
if(is.null(choose_optpoints)) choose_optpoints <- FALSE
if(choose_optpoints){if(is.null(n.optpoints.candidates)) n.optpoints.candidates <- min(1000,10*n.optpoints)}
if(is.null(n.optpoints.candidates)) n.optpoints.candidates <- n.optpoints
if(is.null(unscale.inv.integration.points)) unscale.inv.integration.points <- FALSE
finaldistrib <- integcontrol$finaldistrib
if(is.null(finaldistrib)) finaldistrib <- "sobol" #other possible values are "MC" and "surnew"
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]
d.inv <- length(inv.index);d.opt <- length(opt.index)
n.points <- integcontrol$n.points
if(is.null(n.points)) n.points <- 20
if(finaldistrib=="MC") inv.integration.points <- t(lower.inv + t(matrix(runif(d.inv*n.points),ncol=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="sobol") inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="spec") inv.integration.points <- integcontrol$inv.integration.points
# l'option points d'integration choisis a l'avance est traitee
# on traite donc maintenant le cas importance sampling
if(finaldistrib=="surnew"){
# importance sampling will be perfomed later in the function
# Construction of the points which are 'candidates' to become integration points.
n.candidates <- integcontrol$n.candidates
if(is.null(n.candidates)) n.candidates <- n.points*3
init.distrib <- integcontrol$init.distrib
if(is.null(init.distrib)) init.distrib <- "sobol"
if(init.distrib=="sobol") initial.inv.integration.points <- t(lower.inv+t(sobol(n=n.candidates,dim=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="MC") initial.inv.integration.points <- t(lower.inv+t(matrix(runif(d.inv*n.candidates),ncol=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="spec") initial.inv.integration.points <- integcontrol$inv.integration.points
inv.integration.points <- initial.inv.integration.points
}
if ((class(model@covariance)=="covAffineScaling" | class(model@covariance)=="covScaling") & unscale.inv.integration.points) {
# unscaling of integration points (or candidates)
inv.integration.points <- unscalingFun(mat = inv.integration.points,model = model,indices = inv.index,standardize = TRUE,lower = lower,upper = upper)
}
if(d.inv==1) inv.integration.points <- matrix(inv.integration.points,ncol=1)
# opt.simulation.points now
opt.simulation.points.distrib <- integcontrol$opt.simulation.points.distrib
if(is.null(opt.simulation.points.distrib)) opt.simulation.points.distrib <- "sobol"
if(opt.simulation.points.distrib=="spec") opt.simulation.points <- integcontrol$opt.simulation.points
if(opt.simulation.points.distrib=="sobol"){
addvertices <- integcontrol$addvertices
if(is.null(addvertices)) addvertices <- TRUE
if(addvertices) opt.simulation.points <- rbind(allvertices(d = d.opt,lower = lower.opt,upper = upper.opt) , t(lower.opt + t(sobol(n= n.optpoints.candidates-2^d.opt,dim=d.opt))*(upper.opt-lower.opt)) )
if(!addvertices) opt.simulation.points <-t(lower.opt + t(sobol(n= n.optpoints.candidates,dim=d.opt))*(upper.opt-lower.opt))
}
unscale.opt.simulation.points <- integcontrol$unscale.opt.simulation.points
if(is.null(unscale.opt.simulation.points)) unscale.opt.simulation.points <- TRUE
if ((class(model@covariance)=="covAffineScaling" | class(model@covariance)=="covScaling") & unscale.opt.simulation.points) {
# unscaling of candidate simulation points
opt.simulation.points <- unscalingFun(mat = opt.simulation.points,model = model,indices = opt.index,standardize = TRUE,lower = lower,upper = upper)
}
if(d.opt==1) opt.simulation.points <- matrix(opt.simulation.points,ncol=1)
# Getting pn
if(!is.null(integcontrol$seed)) seed <- integcontrol$seed
res <- get_pn_and_simulate(inv.integration.points=inv.integration.points,
opt.simulation.points=opt.simulation.points,
nsimu=nsimu,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,seed=seed,
lowmemory = (finaldistrib=="surnew"),
subsample_simupoints = choose_optpoints,
unscale=unscale.opt.simulation.points,
p = n.optpoints)
if(finaldistrib!="surnew") return(res)
if(finaldistrib=="surnew"){
#some importance sampling
pn <- res$pn;Tau.n <- pn*(1-pn);Tau.n.sum <- sum(Tau.n)
if(Tau.n.sum==0) Tau.n.sum <- 1
prob.n <- pmax(Tau.n/Tau.n.sum,min.prob/n.candidates);
prob.n <- prob.n/sum(prob.n)
weight.n <- 1/(prob.n*n.candidates*n.points)
prob.n.copy <-c(0,prob.n);
prob.n.cum <-cumsum(prob.n.copy)
myrandoms <- runif(integcontrol$n.points)
my.indices <- findInterval(myrandoms,prob.n.cum,all.inside=TRUE)
pn.my.indices <- pn[my.indices]
integration.weights <- weight.n[my.indices]
inv.integration.points.my.indices <- inv.integration.points[my.indices,]
if(d.inv==1) inv.integration.points.my.indices <- matrix(inv.integration.points.my.indices,ncol=1)
testobj <- get_pn_and_simulate(inv.integration.points=inv.integration.points.my.indices,
opt.simulation.points=opt.simulation.points,
nsimu=2,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,
subsample_simupoints = choose_optpoints,
lowmemory = FALSE,
p = n.optpoints)
allprecomp <- testobj$allprecomp
allKn.inv <- testobj$allKn.inv
allsimupoints <- testobj$allsimupoints
allmn <- testobj$allmn
allsimu <- allsimucentered <- matrix( c(0), nrow=n.optpoints , ncol=n.points*nsimu)
for(i in 1:n.points){
index.i <- my.indices[i]
indice.first <- 1+(index.i-1)*nsimu
indice.last <- index.i*nsimu
indice.first2 <- 1+(i-1)*nsimu
indice.last2 <- i*nsimu
allsimu[,indice.first2:indice.last2] <- res$allsimu[,indice.first:indice.last]
allsimucentered[,indice.first2:indice.last2] <- res$allsimucentered[,indice.first:indice.last]
}
result <- list(inv.integration.points=inv.integration.points.my.indices,
integration.weights=integration.weights,
pn=pn.my.indices,
allsimu=allsimu,
allsimucentered=allsimucentered,
allsimupoints=allsimupoints,
allprecomp=allprecomp,
allKn.inv=allKn.inv,
allmn=allmn
)
return(result)
}
}
}
IRSN/RobustInv documentation built on Nov. 20, 2019, 10:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions integration_design_robinv
Documented in integration_design_robinv
#' @title Constructing a set of integration and simulation points in robust inversion
#'
#' @description Constructing a set of integration and simulation points in robust inversion
#'
#' @param integcontrol An important list with the following possible fields.
#' \code{distrib}: This field needs to be set to either \code{"sur"} or \code{"surnew"}
#' depending on the sampling criterion that is used.
#'
#' (I) If \code{integcontrol$distrib = "surnew"} (default) the
#' user can set the field \code{n.points} (default: 20) which is the number of integration points. The points
#' live in the space of controlled parameters, called Dinv here.
#' The distribution of these n.points points
#' is controlled through the \code{finaldistrib} field. Possible values are
#' \code{"spec"} to set up manually the integration points through the \code{inv.integration.points}
#' field, \code{"MC"} for
#' uniformly distributed points over Dinv, \code{"sobol"} (default) to use the sobol sequence and \code{"surnew"}
#' to use importance sampling based on a density proportional to pn (1-pn) with pn being the
#' excursion probability of the integration points. In that last case, extra fields are required:
#' \code{n.candidates} is the number of points which are candidate to become integration points (default: n.points*3).
#' The important sampling scheme will choose (with replacement) n.points integration points
#' among the candidates. The distribution of the candidates is set with the field
#' \code{init.distrib} which can be either \code{"MC"} for uniform points,
#' \code{"sobol"} (default) to use the sobol
#' sequence or \code{"spec"} to set up manually the matrix of candidates
#' with the \code{inv.integration.points} field.
#' When scaling is used in the kriging model, the candidate points can be unscaled by setting the
#' \code{unscale.inv.integration.points}
#' to \code{TRUE} (default: \code{FALSE}).
#' Regarding the nuisance parameters, the number of simulation points used, e.g., to compute
#' pn is set with the field \code{n.optpoints} (default: 100). There is the option to choose these points
#' among candidates by setting the field \code{choose_optpoints} to \code{TRUE}
#' (default: \code{FALSE}).
#' In that case, the number of candidates is set with the field \code{n.optpoints.candidates}
#' (default: 10*n.optpoints).
#' The field \code{opt.simulation.points.distrib} refers to the distribution of the
#' simulation points if \code{integcontrol$choose_optpoints=FALSE} and to the
#' candidates otherwise. Possible values are \code{"spec"} for points set up manually with the
#' \code{opt.simulation.points} field and \code{"sobol"} (default) for the sobol sequence.
#' When the sobol sequence is used, there is the option to add the vertices of the input
#' domain (of the nuisance parameters) with the field \code{addvertices} (default: \code{TRUE}).
#' When scaling is used in the kriging model, the simulation points (or the candidates)
#' can be unscaled by setting the \code{unscale.opt.simulation.points} field
#' to \code{TRUE} (default: \code{TRUE}).
#' The number of simulation of Gaussian process sample paths per integration points
#' is controlled with the field
#' \code{nsimu} (default: 1000).
#' Finally, when multiple calls are performed, the \code{seed} argument can be overwritten
#' by setting the \code{seed} field.
#'
#' (II) The option \code{integcontrol$distrib = "sur"} needs to be used only when the
#' sur criterion (and not surnew) is used. In that case, the fields \code{n.points} and
#' \code{finaldistrib} work similarly except that the importance sampling option is set
#' with \code{integcontrol$finaldistrib="sur"} and not \code{"surnew"}. In that case,
#' the field \code{init.distrib} works as in (I). Since there is no Gaussian process simulation here,
#' the excursion probability computation is set with the field \code{pnmethod}.
#' Possible values are \code{"fixed"} and \code{"optimMC"} (recommended). The number
#' of MC trials to get pn (through a multivariate normal cdf call) is set with the
#' field \code{pop.size}. Finally, the discretization in the space of nuisance parameters
#' is also obtained with the field \code{n.optpoints} (maximum value: 10).
#' @param d Dimension (i.e. number of scalar input variables) of the objective function.
#' @param lower Array of size d. Lower bound of the input domain.
#' @param upper Array of size d. Upper bound of the input domain.
#' @param opt.index Array with integers corresponding to the indices of the nuisance parameters.
#' @param inv.index Array with integers corresponding to the indices of the controlled parameters.
#' We must have \code{d = length(opt.index) + length(inv.index)}.
#' @param model The current kriging model. km object.
#' @param T Target threshold.
#' @param min.prob Minimum probability for a candidate integration point to be picked.
#' @param seed The random seed. For repeatability.
#' @return Depeding on the fields in the \code{integcontrol} list, the return is different.
#' See the description on the integcontrol argument and the details.
#' @details (I) If the \code{"surnew"} criterion is used, then the output is a list with
#' the following fields:
#' \cr
#' (i) inv.integration.points: n x dinv matrix containing the integration points.
#' n is the number of integration points and dinv the dimension of the space of controlled parameters.
#' \cr
#' (ii) integration.weights: array of size n containing the weights given to each integration points.
#' \cr
#' (iii) pn: array of size n containing the current excursion probabilities for each integration points. These
#' excursion probabilities are based on nsimu conditional Gaussian process simulations.
#' \cr
#' (iv) allsimu: p x (nsimu*n) matrix where p is the number of simulation points.
#' The first nsimu columns are the nsimu simulations in p points for the first integration point,
#' the next nsimu columns are linked to the next integration point and so on.
#' \cr
#' (v) allsimucentered: p x (nsimu*n) matrix containing the later simulations centered by
#' substracting the kriging mean of the points where the simulations are performed.
#' \cr
#' (vi) allsimupoints: (n*p) x d matrix containing all the points where simulations are performed.
#' The first p rows correspond to the simulation points linked to the first integration points, and so on.
#' \cr
#' (vii) allprecomp: a list with 3 fields obtained from a call to the function
#' \code{precomputeUpdateData} on the points \code{allsimupoints}. Help about \code{precomputeUpdateData}
#' is given in the KrigInv package.
#' \cr
#' (viii) allKn.inv: (n*p) x p matrix containing n different p x p matrices. Matrix number i is the inverse
#' of the p x p non-conditional covariance matrix of the p simulation points associated to the
#' integration point i.
#' \cr
#' (ix) allmn: n*p array with the kriging means of all n*p simulation points stored in \code{allsimupoints}.
#' \cr
#' (II) If the \code{"sur"} criterion is used, then the output is still a list, which however has different fields
#' than in (I):
#' \cr
#' (i) integration.points : (p*n) x d matrix containing all the simulation points. The first p points are
#' associated to the first integration point and so on.
#' \cr
#' (ii) cov: see description of the field allKn.inv (viii) in (I).
#' \cr
#' (iii) mean: see description of the field allmn (ix) in (I).
#' \cr
#' (iv) pn: array of size n containing the current excursion probabilities for each integration points.
#' @export
#' @author Clement Chevalier \email{clement.chevalier@@unine.ch}
#' @examples
#' library(KrigInv)
#' library(randtoolbox)
#' myfun <- branin_robinv
#' d <- 3
#'
#' set.seed(8)
#'
#' n0 <- 30
#' T <- 10
#' opt.index <- c(3)
#' inv.index <- c(1,2)
#' lower <- rep(0,times=d)
#' upper <- rep(1,times=d)
#'
#' design <- matrix(runif(d*n0),nrow=n0)
#' response <- myfun(design)
#' model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2")
#'
#' integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
#' n.candidates=50,nsimu=1000,n.optpoints = 50,
#' choose_optpoints=TRUE,n.optpoints.candidates=500)
#' \dontrun{
#' obj <- integration_design_robinv(integcontrol = integcontrol,d=d,lower=lower,upper=upper,
#' opt.index=opt.index,inv.index=inv.index,model=model,T=T)
#' }
#' #####################################
#' # An example with scaling
#' library(KrigInv)
#' myfun <- function(x){ return(-1*branin_robinv(x) - 50*sin(min(100,1/x[3])) ) }
#' d <- 3
#'
#' set.seed(8)
#'
#' n0 <- 60
#' T <- 30
#' opt.index <- c(3)
#' inv.index <- c(1,2)
#' lower <- rep(0,times=d)
#' upper <- rep(1,times=d)
#'
#' design <- matrix(runif(d*n0),nrow=n0)
#' response <- apply(X = design,FUN = myfun,MARGIN = 1)
#' knots.number <- c(0,0,3)
#' knots <- generate_knots(knots.number = knots.number , d = d)
#'
#' model <- km(formula = ~1,design = design,response = response,covtype = "matern3_2",scaling = TRUE,knots=knots)
#'
#' integcontrol <- list(distrib = "surnew",n.points = 20,finaldistrib="surnew",
#' n.candidates=200,nsimu=100,n.optpoints = 50,
#' choose_optpoints=TRUE,n.optpoints.candidates=500,
#' unscale.opt.simulation.points=TRUE)
#' \dontrun{
#' obj <- integration_design_robinv(integcontrol = integcontrol,d=d,lower=lower,upper=upper,
#' opt.index=opt.index,inv.index=inv.index,model=model,T=T)
#' }
integration_design_robinv <- function(integcontrol=NULL,d=NULL,lower,upper,
opt.index,inv.index,
model,T,
min.prob=0.001,seed=NULL){
result<-NULL
if(!is.null(seed)) set.seed(seed)
if(is.null(d)) d <- length(lower)
if (length(lower) != length(upper) ) print("integration_design_robinv: 'lower' and 'upper' must have the same length")
if(is.null(integcontrol$distrib)) integcontrol$distrib <- "surnew"
if(integcontrol$distrib=="sur"){
# only for the "sur" criterion.
n.optpoints <- integcontrol$n.optpoints
if(is.null(n.optpoints)) n.optpoints <- 10
if(n.optpoints > 10) n.optpoints <- 10 # maximum value because we cannot go in dimension higher than 20
finaldistrib <- integcontrol$finaldistrib
if(is.null(finaldistrib)) finaldistrib <- "sobol" # other possible value is "sur"
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]
d.inv <- length(inv.index);d.opt <- length(opt.index)
n.points <- integcontrol$n.points
if(is.null(n.points)) n.points <- 20
if(finaldistrib=="MC") inv.integration.points <- t(lower.inv + t(matrix(runif(d.inv*n.points),ncol=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="sobol") inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="spec") inv.integration.points <- integcontrol$inv.integration.points
# We now deal with the importance sampling case
if(finaldistrib=="sur"){
# construction of the candidates
n.candidates <- integcontrol$n.candidates
if(is.null(n.candidates)) n.candidates <- n.points*5
init.distrib <- integcontrol$init.distrib
if(is.null(init.distrib)) init.distrib <- "sobol"
if(init.distrib=="sobol") initial.inv.integration.points <- t(lower.inv+t(sobol(n=n.candidates,dim=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="MC") initial.inv.integration.points <- t(lower.inv+t(matrix(runif(d.inv*n.candidates),ncol=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="spec") initial.inv.integration.points <- integcontrol$inv.integration.points
inv.integration.points <- initial.inv.integration.points
}
if(d.inv==1) inv.integration.points <- matrix(inv.integration.points,ncol=1)
pnmethod <- integcontrol$pnmethod
if(is.null(pnmethod)) pnmethod <- "fixed" # other possible value is "optimMC"
if(pnmethod=="fixed"){
# in Dopt, some points are generated (and fixed) using the sobol sequence (option 1)
# or uniformly (option 2)
# This option is computationally cheap; but is expected to give excursion probabilities
# that might be widely underestimated.
N <- nrow(inv.integration.points)
integration.points <- matrix(c(0),nrow=N*n.optpoints,ncol=d)
if(is.null(integcontrol$pnmethod.fixed.option)) integcontrol$pnmethod.fixed.option <- 1
if(integcontrol$pnmethod.fixed.option==1){
opt.integration.points <- t(lower.opt + t(sobol(n = n.optpoints, dim = d.opt))*(upper.opt - lower.opt))
if(d.opt==1) opt.integration.points <- matrix(opt.integration.points,ncol=1)
indices <- expand.grid(c(1:n.optpoints),c(1:N))
index1 <- indices[,1]
index2 <- indices[,2]
integration.points[,opt.index] <- as.matrix(opt.integration.points[index1,])
integration.points[,inv.index] <- as.matrix(inv.integration.points[index2,])
}else{
opt.integration.points <- t(lower.opt + t(matrix(runif(n.optpoints*N*d.opt),ncol=d.opt))*(upper.opt - lower.opt))
if(d.opt==1) opt.integration.points <- matrix(opt.integration.points,ncol=1)
indices <- expand.grid(c(1:n.optpoints),c(1:N))
indices[,1] <- c(1:(n.optpoints*N))
index1 <- indices[,1]
index2 <- indices[,2]
integration.points[,opt.index] <- as.matrix(opt.integration.points[index1,])
integration.points[,inv.index] <- as.matrix(inv.integration.points[index2,])
}
#calculation of all cov matrices, means and pn
res <- optim_pn_MC(inv.integration.points=inv.integration.points,model=model,T=T,
opt.index=opt.index,inv.index=inv.index,lower=lower,upper=upper,
n.optpoints=n.optpoints,pop.size=1,
MC=FALSE,deterministicMat=integration.points)
}
if(pnmethod=="optimMC"){
pop.size <- integcontrol$pop.size
if(is.null(pop.size)) pop.size <- 50
res <- optim_pn_MC(inv.integration.points=inv.integration.points,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,
n.optpoints=n.optpoints,pop.size=pop.size,
MC=TRUE)
}
if(finaldistrib!="sur") return(res)
if(finaldistrib=="sur"){
# some importance sampling
pn <- res$pn;Tau.n <- pn*(1-pn);Tau.n.sum <- sum(Tau.n)
if(Tau.n.sum==0) Tau.n.sum <- 1
prob.n <- pmax(Tau.n/Tau.n.sum,min.prob/n.candidates);prob.n <- prob.n/sum(prob.n)
weight.n <- 1/(prob.n*n.candidates*n.points)
prob.n.copy <-c(0,prob.n);prob.n.cum <-cumsum(prob.n.copy)
my.indices <- findInterval(runif(integcontrol$n.points),prob.n.cum,all.inside=TRUE)
pn.my.indices <- pn[my.indices]
integration.weights <- weight.n[my.indices]
new.mean <- rep(0,times=n.optpoints*n.points)
old.mean <- res$mean
new.cov <- matrix(c(0),nrow=n.optpoints*n.points,ncol=n.optpoints)
old.cov <- res$cov
new.intpoints <- matrix(c(0),nrow=n.optpoints*n.points,ncol=d)
old.intpoints <- res$integration.points
for(i in 1:n.points){
index.i <- my.indices[i]
new.mean[(1+(i-1)*n.optpoints):(i*n.optpoints)] <- old.mean[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints)]
new.cov[(1+(i-1)*n.optpoints):(i*n.optpoints),] <- old.cov[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints),]
new.intpoints[(1+(i-1)*n.optpoints):(i*n.optpoints),] <- old.intpoints[(1+(index.i-1)*n.optpoints):(index.i*n.optpoints),]
}
result <- list(integration.points=new.intpoints,
integration.weights=integration.weights,
mean=new.mean,
cov=new.cov,
pn=pn.my.indices)
return(result)
}
}
if(integcontrol$distrib=="surnew"){
## only for the "surnew" criterion.
n.optpoints <- integcontrol$n.optpoints
choose_optpoints <- integcontrol$choose_optpoints
n.optpoints.candidates <- integcontrol$n.optpoints.candidates
nsimu <- integcontrol$nsimu
unscale.inv.integration.points <- integcontrol$unscale.inv.integration.points
if(is.null(nsimu)) nsimu <- 1000
if(is.null(n.optpoints)) n.optpoints <- 100
if(n.optpoints > 1000) n.optpoints <- 1000
if(is.null(choose_optpoints)) choose_optpoints <- FALSE
if(choose_optpoints){if(is.null(n.optpoints.candidates)) n.optpoints.candidates <- min(1000,10*n.optpoints)}
if(is.null(n.optpoints.candidates)) n.optpoints.candidates <- n.optpoints
if(is.null(unscale.inv.integration.points)) unscale.inv.integration.points <- FALSE
finaldistrib <- integcontrol$finaldistrib
if(is.null(finaldistrib)) finaldistrib <- "sobol" #other possible values are "MC" and "surnew"
lower.inv <- lower[inv.index];upper.inv <- upper[inv.index]
lower.opt <- lower[opt.index];upper.opt <- upper[opt.index]
d.inv <- length(inv.index);d.opt <- length(opt.index)
n.points <- integcontrol$n.points
if(is.null(n.points)) n.points <- 20
if(finaldistrib=="MC") inv.integration.points <- t(lower.inv + t(matrix(runif(d.inv*n.points),ncol=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="sobol") inv.integration.points <- t(lower.inv + t(sobol(n=n.points,dim=d.inv))*(upper.inv-lower.inv))
if(finaldistrib=="spec") inv.integration.points <- integcontrol$inv.integration.points
# l'option points d'integration choisis a l'avance est traitee
# on traite donc maintenant le cas importance sampling
if(finaldistrib=="surnew"){
# importance sampling will be perfomed later in the function
# Construction of the points which are 'candidates' to become integration points.
n.candidates <- integcontrol$n.candidates
if(is.null(n.candidates)) n.candidates <- n.points*3
init.distrib <- integcontrol$init.distrib
if(is.null(init.distrib)) init.distrib <- "sobol"
if(init.distrib=="sobol") initial.inv.integration.points <- t(lower.inv+t(sobol(n=n.candidates,dim=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="MC") initial.inv.integration.points <- t(lower.inv+t(matrix(runif(d.inv*n.candidates),ncol=d.inv))*(upper.inv-lower.inv))
if(init.distrib=="spec") initial.inv.integration.points <- integcontrol$inv.integration.points
inv.integration.points <- initial.inv.integration.points
}
if ((class(model@covariance)=="covAffineScaling" | class(model@covariance)=="covScaling") & unscale.inv.integration.points) {
# unscaling of integration points (or candidates)
inv.integration.points <- unscalingFun(mat = inv.integration.points,model = model,indices = inv.index,standardize = TRUE,lower = lower,upper = upper)
}
if(d.inv==1) inv.integration.points <- matrix(inv.integration.points,ncol=1)
# opt.simulation.points now
opt.simulation.points.distrib <- integcontrol$opt.simulation.points.distrib
if(is.null(opt.simulation.points.distrib)) opt.simulation.points.distrib <- "sobol"
if(opt.simulation.points.distrib=="spec") opt.simulation.points <- integcontrol$opt.simulation.points
if(opt.simulation.points.distrib=="sobol"){
addvertices <- integcontrol$addvertices
if(is.null(addvertices)) addvertices <- TRUE
if(addvertices) opt.simulation.points <- rbind(allvertices(d = d.opt,lower = lower.opt,upper = upper.opt) , t(lower.opt + t(sobol(n= n.optpoints.candidates-2^d.opt,dim=d.opt))*(upper.opt-lower.opt)) )
if(!addvertices) opt.simulation.points <-t(lower.opt + t(sobol(n= n.optpoints.candidates,dim=d.opt))*(upper.opt-lower.opt))
}
unscale.opt.simulation.points <- integcontrol$unscale.opt.simulation.points
if(is.null(unscale.opt.simulation.points)) unscale.opt.simulation.points <- TRUE
if ((class(model@covariance)=="covAffineScaling" | class(model@covariance)=="covScaling") & unscale.opt.simulation.points) {
# unscaling of candidate simulation points
opt.simulation.points <- unscalingFun(mat = opt.simulation.points,model = model,indices = opt.index,standardize = TRUE,lower = lower,upper = upper)
}
if(d.opt==1) opt.simulation.points <- matrix(opt.simulation.points,ncol=1)
# Getting pn
if(!is.null(integcontrol$seed)) seed <- integcontrol$seed
res <- get_pn_and_simulate(inv.integration.points=inv.integration.points,
opt.simulation.points=opt.simulation.points,
nsimu=nsimu,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,seed=seed,
lowmemory = (finaldistrib=="surnew"),
subsample_simupoints = choose_optpoints,
unscale=unscale.opt.simulation.points,
p = n.optpoints)
if(finaldistrib!="surnew") return(res)
if(finaldistrib=="surnew"){
#some importance sampling
pn <- res$pn;Tau.n <- pn*(1-pn);Tau.n.sum <- sum(Tau.n)
if(Tau.n.sum==0) Tau.n.sum <- 1
prob.n <- pmax(Tau.n/Tau.n.sum,min.prob/n.candidates);
prob.n <- prob.n/sum(prob.n)
weight.n <- 1/(prob.n*n.candidates*n.points)
prob.n.copy <-c(0,prob.n);
prob.n.cum <-cumsum(prob.n.copy)
myrandoms <- runif(integcontrol$n.points)
my.indices <- findInterval(myrandoms,prob.n.cum,all.inside=TRUE)
pn.my.indices <- pn[my.indices]
integration.weights <- weight.n[my.indices]
inv.integration.points.my.indices <- inv.integration.points[my.indices,]
if(d.inv==1) inv.integration.points.my.indices <- matrix(inv.integration.points.my.indices,ncol=1)
testobj <- get_pn_and_simulate(inv.integration.points=inv.integration.points.my.indices,
opt.simulation.points=opt.simulation.points,
nsimu=2,
model=model,T=T,
opt.index=opt.index,inv.index=inv.index,
lower=lower,upper=upper,
subsample_simupoints = choose_optpoints,
lowmemory = FALSE,
p = n.optpoints)
allprecomp <- testobj$allprecomp
allKn.inv <- testobj$allKn.inv
allsimupoints <- testobj$allsimupoints
allmn <- testobj$allmn
allsimu <- allsimucentered <- matrix( c(0), nrow=n.optpoints , ncol=n.points*nsimu)
for(i in 1:n.points){
index.i <- my.indices[i]
indice.first <- 1+(index.i-1)*nsimu
indice.last <- index.i*nsimu
indice.first2 <- 1+(i-1)*nsimu
indice.last2 <- i*nsimu
allsimu[,indice.first2:indice.last2] <- res$allsimu[,indice.first:indice.last]
allsimucentered[,indice.first2:indice.last2] <- res$allsimucentered[,indice.first:indice.last]
}
result <- list(inv.integration.points=inv.integration.points.my.indices,
integration.weights=integration.weights,
pn=pn.my.indices,
allsimu=allsimu,
allsimucentered=allsimucentered,
allsimupoints=allsimupoints,
allprecomp=allprecomp,
allKn.inv=allKn.inv,
allmn=allmn
)
return(result)
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.