Nothing
R/estimateSUR.R
Defines functions t.list estimateSUR
Documented in estimateSUR
#' EstimateSUR
#'
#' A function for estimating a SUR criterion with a realisation of a PPP
# @importFrom Rcpp sourceCpp evalCpp
# @useDynLib BMP
estimateSUR <- function(PPP,
#' @param PPP the Poisson point process generated to get alpha.
xi_PPP_X,
#' @param xi_PPP_X the output of xi(cbind(PPP$X, PPP$final_X)).
integrated = TRUE,
#' @param integrated boolean to specify of SUR criterion is standard or
#' integrated.
N_ppp,
#' @param N_ppp the number of Poisson processes used for the SUR criterion estimation.
method = "discrete",
#' @param method eiter "genoud" for an optimisation using the package \code{rgenoud} or
#' "discrete" for a discrete search over \code{SUR_pop}.
SUR_pop,
#' @param SUR_pop if \code{optimcontrol$method=="discrete"}, SUR_pop is the population
#' onto which minimizer is sought. Should
#' be a matrix d x n.
r = N.batch,
#' @param r number of points to be added to the DoE.
optimcontrol = list(pop.size = 50*d,
max.generations = 10*d),
#' @param optimcontrol a list of control parameters for the optimisation
#' of the SUR criterion using the \code{rgenoud} package.
approx.pnorm,
#' @param approx.pnorm (optional) an approximation of base pnorm function running faster.
J = 0,
#' @param J the center of an interval of size 8 for pnorm approximation.
N.batch = foreach::getDoParWorkers(),
#' @param N.batch Number of batchs for parallel computation.
verbose = 0,
#' @param verbose to control the print level of the algorithm
...
#' @param ... further arguments to be passed to fSUR.
){
# fix foreach NOTE
X <- NULL
N <- length(PPP$final_U)
d <- dim(PPP$final_X)[1]
n <- dim(xi_PPP_X$Kinv_kn)[1]
if(method=='discrete' && missing(SUR_pop)){
SUR_pop = cbind(PPP$X, PPP$final_X)
u_SUR_pop <- c(PPP$U, PPP$final_U)
N_pop = dim(SUR_pop)[2]
xi_SUR_POP <- xi_PPP_X
SUR_aug <- rbind(SUR_pop, matrix(xi_SUR_POP$mean, nrow = 1), matrix(xi_SUR_POP$sd, nrow = 1), xi_SUR_POP$kn, xi_SUR_POP$Kinv_kn)
}
if(!missing(N_ppp) && N_ppp<N) {
ind <- sample(x = 1:N, size = N_ppp, replace = FALSE)
sel1 <- names(PPP$U)%in%ind
PPP$X <- PPP$X[,sel1]
PPP$U <- PPP$U[sel1]
sel2 <- names(PPP$final_U)%in%ind
PPP$final_X <- PPP$final_X[,sel2]
PPP$final_U <- PPP$final_U[sel2]
xi_PPP_X$mean <- xi_PPP_X$mean[c(sel1, sel2)]
xi_PPP_X$sd <- xi_PPP_X$sd[c(sel1, sel2)]
xi_PPP_X$kn <- as.matrix(xi_PPP_X$kn[,c(sel1, sel2)])
xi_PPP_X$Kinv_kn <- as.matrix(xi_PPP_X$Kinv_kn[,c(sel1, sel2)])
N <- N_ppp
}
x = PPP$final_X
u = PPP$final_U
u_lpa <- NULL
xi_x <- xi_PPP_X
if(integrated == TRUE){
x = cbind(PPP$X, x)
u = c(PPP$U, u)
u_lpa <- rep(PPP$U, each = N)
} else {
xi_x$mean = tail(xi_PPP_X$mean, N)
xi_x$sd = tail(xi_PPP_X$sd, N)
xi_x$kn = t(tail(t(xi_PPP_X$kn), N))
xi_x$Kinv_kn = t(tail(t(xi_PPP_X$Kinv_kn), N))
}
args = list(...)
args$x = x
args$u = u
args$u_lpa = u_lpa
args$xi_x = xi_x
args$integrated = integrated
args$N = N
if(approx.pnorm) {
x_pnorm <- seq(from = -4, to = 4, by = 0.2)
args$pnorm = stats::approxfun(x_pnorm, pnorm(x_pnorm), yleft = 0, yright = 1)
}
args$Xnew <- Xnew <- NULL
args$xi_Xnew <- xi_Xnew <- NULL
cat(" * Evaluation of SUR criterion: integrated = ",integrated, ", r = ", ifelse(is.null(args$r), 1, args$r), ", approx = ",ifelse(is.null(args$approx), FALSE, args$approx),", approx.pnorm = ",approx.pnorm,", optim = ", method,", N_ppp = ", ifelse(missing(N_ppp),N, N_ppp)," \n", sep = "")
sur = list(x = NULL, u = NULL, t = NULL)
for(k in 1:r){
if(method=="discrete"){
SUR <- foreach::foreach(X = iterators::iter(SUR_aug, by ='col', chunksize = ceiling(N_pop/N.batch)),
.combine = 'c',
.export = 'fSUR',
.options.multicore = list(set.seed = TRUE),
.errorhandling = "pass") %dopar% {
tmp <- c(apply(X, 2, function(xaug){
args$xnew <- xaug[1:d]
args$xi_xnew <- list(mean = xaug[d+1], sd = xaug[d+2], kn = xaug[(d+3):(d+2+n)], Kinv_kn = tail(xaug, n))
do.call(fSUR, args)
}))
return(tmp)
}
# if(class(SUR)!="numeric"){
if(!(inherits(SUR,"numeric"))){
message(' ! memory issue with parallel computing, approximated SUR used insteead !')
args$approx = TRUE
x_pnorm <- seq(from = -4, to = 4, by = 0.2)
args$pnorm = stats::approxfun(x_pnorm, pnorm(x_pnorm), yleft = 0, yright = 1)
SUR <- foreach::foreach(X = iterators::iter(SUR_aug, by ='col', chunksize = ceiling(N_pop/N.batch)),
.combine = 'c',
.export = 'fSUR',
.options.multicore = list(set.seed = TRUE),
.errorhandling = "pass") %dopar% {
tmp <- c(apply(X, 2, function(xaug){
args$xnew <- xaug[1:d]
args$xi_xnew <- list(mean = xaug[d+1], sd = xaug[d+2], kn = xaug[(d+3):(d+2+n)], Kinv_kn = tail(xaug, n))
do.call(fSUR, args)
}))
return(tmp)
}
}
# if(class(SUR)!="numeric"){
if(!(inherits(SUR,"numeric"))){
message(' ! memory issue with parallel computing, standard SUR used insteead !')
args$integrated = FALSE
SUR <- foreach::foreach(X = iterators::iter(SUR_aug, by ='col', chunksize = ceiling(N_pop/N.batch)),
.combine = 'c',
.export = 'fSUR',
.options.multicore = list(set.seed = TRUE),
.errorhandling = "pass") %dopar% {
tmp <- c(apply(X, 2, function(xaug){
args$xnew <- xaug[1:d]
args$xi_xnew <- list(mean = xaug[d+1], sd = xaug[d+2], kn = xaug[(d+3):(d+2+n)], Kinv_kn = tail(xaug, n))
do.call(fSUR, args)
}))
return(tmp)
}
}
# if(class(SUR)!="numeric"){
if(!(inherits(SUR,"numeric"))){
print(SUR)
}
ind_min = which.min(SUR)
SUR_point <- SUR_pop[,ind_min]
SUR_aug <- SUR_aug[,-ind_min]
sur$u = c(sur$u, u_SUR_pop[ind_min])
sur$t = c(sur$t, ind_min/length(SUR))
} else {
f_gen <- function(xnew) {
args$xnew = xnew
do.call(fSUR, args)
}
SUR_point <- do.call(rgenoud::genoud, c(optimcontrol, list(fn = f_gen, nvars = d, print.level = verbose)))$par
}
Xnew = cbind(Xnew, as.matrix(SUR_point))
sur$x = cbind(sur$x, SUR_point)
}
return(sur)
#' @return a list containing the points minimising the criterion
}
t.list <- function(l){
lapply(split(do.call("c", l), names(l[[1]])), unname)
}
Try the mistral package in your browser
Any scripts or data that you put into this service are public.
