# R/minimax.R In minimaxdesign: Minimax and Minimax Projection Designs

#### Documented in minimaxminimax.mapminiMaxPro

```########################################################################
# Main function for mMcPSO
########################################################################
minimax <- function(N,p,q=10,region="hypercube",ini=NA,const=NA,clust_pts=NA,
params_pso=list(w=0.72,c1=1.49,c2=1.49),
npart=5,nclust=1e5,neval=nclust,
itmax_pso=50,itmax_pp=100,itmax_inn=1e4,jit=0.1/sqrt(N)){

#Setting transformation type
regionby=ifelse(p>2,1e-3,-1) #for simplex transformation
bd <- c(0,1) #Default limits are (0,1)

switch(region,
hypercube = {
tf <- function(D,by,num_proc){return(D)}
checkBounds <- function(D){ #function for checking bound
return(D)
}
},
simplex = {
tf <- CtoA;
checkBounds <- function(D){ #function for checking bound
return(D)
}
},
ball = {
tf <- CtoB;
bd <- c(-1,1);
checkBounds <- function(D){ #function for checking bound
good.idx <- which(rowSums(D^2)<=1)
return(D[good.idx,])
}
},
ineq = {
tf <- function(D,by,num_proc){ #randomly sample until enough points
num_pts <- nrow(D)
samp <- matrix(NA,nrow=num_pts,ncol=p)
cur_pts <- 0
while (cur_pts < num_pts){
# print(paste0("cur_pts: ",cur_pts))
xx <- stats::runif(p)
if (const(xx)){
samp[cur_pts+1,] <- xx
cur_pts <- cur_pts + 1
}
}
return(samp)
}
checkBounds <- function(D){ #function for checking bound
good.idx <- apply(D,1,const)
return(D[good.idx,])
}
},
custom = { #custom clustering points
tf <- function(D,by,num_proc){return(D)}
checkBounds <- function(D){ #function for checking bound
return(D)
}
bd <- c(-1e8,1e8)
}
)
by <- regionby

# Setting clustering points
if (any(is.na(clust_pts))){
print("Generating clustering points ...")
# point = tf(randtoolbox::sobol(nclust,p),by,parallel::detectCores()) #clustering points
clust_pts = tf(as.matrix(Lattice(max(conf.design::primes(nclust)),p))) #clustering points
clust_pts <- checkBounds(clust_pts)

#Generate eval_pts
# eval_pts = tf(randtoolbox::sobol(neval,p),by,parallel::detectCores()) #minimax approximation points for post-processing
eval_pts = tf(as.matrix(Lattice(max(conf.design::primes(neval)),p,shift=TRUE))) #minimax approximation points for post-processing
# eval_pts = tf(as.matrix(Lattice(max(conf.design::primes(neval)),p))) #minimax approximation points for post-processing
if (region=="hypercube"){
eval_pts = rbind(eval_pts,tf(gtools::permutations(3,p,c(0.0,0.5,1.0),repeats.allowed=TRUE))) #add 3^p design
}
eval_pts <- checkBounds(eval_pts)
}else{
print("Using custom clustering points ...")
bd <- c(-1e8,1e8)
eval_pts <- clust_pts
}

#   offset = seq(0,1-1/part_num,1/part_num)

# Generate initial center particles

#automatically set initialization flag
by.clust <- 1e-4
if (any(is.na(ini))){
if (region=="hypercube"){
if ( (is.whole(log(N)/log(2))) && (N <= 2^p) ){ # ff initialization
ini <- "ff"
}else if( (N>2^p) && (N<=1.25*2^p) ){
ini <- "ff"
}else{
ini <- "sobol"
}
}else{
ini <- "sobol"
}
}
# Initialize
if (region=="custom"){
cluster_center = clust_pts[sample(1:nrow(clust_pts),N),] #initial cluster centers
for (i in 1:(npart-1)){
cluster_center = cbind(cluster_center,clust_pts[sample(1:nrow(clust_pts),N),])
}
}else if (any(is.matrix(ini))){
cluster_center = tf(ini,by.clust,parallel::detectCores()) #initial cluster centers
for (i in 1:(npart-1)){ #add scrambled randtoolbox::sobol sets
cluster_center = cbind(cluster_center,tf(jitter(ini),by.clust,parallel::detectCores()))
}
}else if(any(is.list(ini))){
npart <- length(ini)
cluster_center <- c()
for (i in 1:npart){
cluster_center = cbind(cluster_center,tf(ini[[i]],by.clust,parallel::detectCores()))
}
}else if (ini=="sobol"){ # initialize via randomized sobol'
cluster_center = tf(randtoolbox::sobol(N,p),by.clust,parallel::detectCores()) #initial cluster centers
for (i in 1:(npart-1)){ #add scrambled randtoolbox::sobol sets
cluster_center = cbind(cluster_center,tf(randtoolbox::sobol(N,p,init=FALSE),by.clust,parallel::detectCores()))
}
}else if (ini=="ff"){ # initialize via 2^(k-p) fractional factorial
if ( N >= 2*(2^p) ){
stop("Design size too large for FF initialization!")
}else{
pwf = floor(log(N)/log(2));
cur_center = tf(0.25*DoE.base::desnum(FrF2::FrF2(2^pwf,p))+0.5,by.clust,parallel::detectCores()) #initial cluster centers
cur_center = rbind(cur_center,tf(jitter(randtoolbox::sobol(N-2^pwf,p)),by.clust,parallel::detectCores()))
cluster_center = cur_center
for (i in 1:(npart-1)){
cur_center = tf(0.25*DoE.base::desnum(FrF2::FrF2(2^pwf,p))+0.5,by.clust,parallel::detectCores()) #initial cluster centers
cur_center = rbind(cur_center,tf(jitter(randtoolbox::sobol(N-2^pwf,p,init=FALSE)),by.clust,parallel::detectCores()))
cluster_center = cbind(cluster_center,cur_center)
}
}
}

# print("4")
# mMc-PSO: function coded in C++
t1 = Sys.time()
# print(dim(clust_pts))
# print(dim(eval_pts))
# print(dim(cluster_center))
# Sys.sleep(0.5)
# print(bd)
D = kmeanspso(clust_pts, eval_pts, cluster_center, q, 2.0,
params_pso\$w,params_pso\$c1,params_pso\$c2,
2*npart, itmax_pso, itmax_pp, 1000, 1000, 1e-4, 1e-4,
1e-4, itmax_inn,
parallel::detectCores(), jit, bd[1], bd[2], rep(0,N))
fitpre <- D\$gbes_centers
fitpost <- fitpre;

t2 = Sys.time()
tm <- difftime(t2,t1)
# print(5)

# D\$gbes_centers <- fitpost
# D\$time <- tm
return(D\$gbes_centers)
# return(eval_pts)
}

########################################################################
# Functions for minimax projection designs
########################################################################

#Main function for generating minimax projection designs
miniMaxPro <- function(N,p,mMdes=NA, mMtol=1e-3*p,
neval=1e5, itmax_refine=100, ...){
#Generate minimax design from mMc-PSO
if (is.na(mMdes)){
mMdes <- minimax(N,p,...)
}
# if (is.na(refine_pts)){
# refine_pts <- randtoolbox::sobol(refine_num,p)
# }
eval_pts = as.matrix(Lattice(max(conf.design::primes(neval)),p,shift=TRUE)) #minimax approximation points for post-processing
# eval_pts = as.matrix(Lattice(max(conf.design::primes(neval)),p)) #minimax approximation points for post-processing
eval_pts = rbind(eval_pts,gtools::permutations(3,p,c(0.0,0.5,1.0),repeats.allowed=TRUE)) #add 3^p design

#Refinement
ret <- refine(mMdes,eval_pts=eval_pts,pp_itmax=itmax_refine,mM_tol=mMtol)

return( list(minimax=mMdes,miniMaxPro=ret) )

}

#Refinement step for miniMaxPro designs
refine <- function(mMdes,eval_pts=NA,
pp_itmax=50,mM_tol=1e-3*(ncol(mMdes)),jitter.scl=5){

mp=function(xx)
{
#Blockwise MaxPro objective
dif <- DD - rep(1,dim(DD)[1])%*%t(xx)
d=apply(dif, 1, prod)
val=sum(1/(d^2+10^(-10))) #For numerical stability
lfn <- log(val)

B <- (t(1/dif)%*%matrix(d,ncol=1))
G <- 2*B/val

}

con=function(xx)
{
#Constraint function
return(list("constraints"= sum((xx-cur_pt)^2) - toler^2,
"jacobian"= 2*(xx-cur_pt) ))
}

#Compute minimax criterion
N <- nrow(mMdes)
p <- ncol(mMdes)
orig_vec <- mMcrit_allpts( mMdes, eval_pts )
orig_mM <- max(mMcrit_allpts( mMdes, eval_pts )) #Current minimax distance

numCores <- parallel::detectCores()
cluster = parallel::makeCluster(numCores,type="SOCK")
doParallel::registerDoParallel(cluster)

cur_des <- mMdes
itcur <- 0

print(paste0("MaxPro refinement ..."))
while (itcur < pp_itmax){
if(itcur>1) printBar(itcur/pp_itmax)

#Blockwise MaxPro for each design point
dst_vec <- mMcrit_allpts( cur_des, eval_pts ) #First compute the minimum distance for each design point
dst_vec <- pmax(orig_mM + mM_tol - dst_vec,0)
mM_ind <- which(dst_vec==0)
prev_des <- cur_des

indices <- gtools::permute(setdiff(1:N,mM_ind))
flg_cont <- F

while (!flg_cont){
biglist <- foreach::'%dopar%'(foreach::foreach (i = indices, .packages="nloptr"),
# biglist <- foreach::'%do%'(foreach::foreach (i = indices, .packages="nloptr"),
{
toler <- dst_vec[i]/jitter.scl # Ensures the jittered initial points do not exceed minimax
# print(toler)
cur_pt <- cur_des[i,] #Current point
DD <- cur_des[-i,] #Design without current point

#         ini_pt <- pmin(pmax(cur_des[i,]+toler/sqrt(p)*stats::runif(p,min=-1),0),1)
ini_pt <- pmin(pmax(cur_pt+toler*stats::runif(p,min=-1),0),1)
# ini_pt <- cur_pt
tryCatch({a=nloptr::nloptr( x0=ini_pt,
mp,lb=rep(-1,p),ub=rep(1,p),eval_g_ineq = con,
opts = list("algorithm"="NLOPT_LD_MMA","maxeval"=100))
list("obj"=a\$objective,"cur_sol"=a\$sol,"prev_obj"=mp(cur_des[i,])\$objective)}
,error = function(err){
list("obj"=mp(cur_des[i,])\$objective,"cur_sol"=cur_des[i,],"prev_obj"=mp(cur_des[i,])\$objective)
})

})

prev_obj <- rep(0,length(indices))
cur_obj <- rep(0,length(indices))
solvec <- vector("list",length(indices))
for (i in 1:length(indices)){
#If there is an improvement, then change. Otherwise, stick with original point
if (biglist[[i]]\$obj < biglist[[i]]\$prev_obj){
cur_des[indices[i],] <- biglist[[i]]\$cur_sol
}
}
flg_cont=T

#       #Check if minimax distance is violated
#       cur_dst_vec <- mMcrit_allpts( cur_des, eval_pts )
#       if (max(cur_dst_vec)-orig_mM <= mM_tol){
#         #... terminate loop
#         flg_cont = T
#         scl <- jitter.scl
#       }else{
#         cur_des <- prev_des
#         scl <- scl + 5
# #         print(paste0("Decreasing jitter: ", scl))
#
#         #If scl hits maxscl, terminate algorithm
#         if (scl >= maxscl){
#           cur_des <- prev_des
#           flg_cont = T
#           itcur <- pp_itmax
#         }
#
#       }

}

#     plot(mMdes,xlim=c(0,1),ylim=c(0,1))
#     points(cur_des,col="red",pch=4)

#increment
itcur <- itcur + 1
}

parallel::stopCluster(cluster)

return(cur_des)
}

########################################################################
# Functions for generating minimax designs from images
########################################################################

minimax.map = function(N,img,p=2,q=10,
params_pso=list(w=0.72,c1=1.49,c2=1.49),
npart=5,nclust=1e5,neval=nclust,
itmax_pso=50,itmax_pp=100,itmax_inn=1e4,jit=0.1/sqrt(N)){

mapfile <- round(img)

#Get evaluation points for minimax PSO
eval_pts <- which(mapfile==0,arr.ind=T)
eval_pts[,1] <- eval_pts[,1]/nrow(mapfile)
eval_pts[,2] <- eval_pts[,2]/ncol(mapfile)

# Get clustering data
clust_pts <- eval_pts[sample.int(nrow(eval_pts),nclust),]
# inv_pt <- cbind(point[,2],point[,1])
# inv_pt[,2] <- 1 - inv_pt[,2]
# plot(inv_pt)

#Generate initial centers
cluster_center <- eval_pts[sample.int(nrow(eval_pts),N),]
for (i in 1:(npart[1]-1)){
cluster_center = cbind(cluster_center, eval_pts[sample.int(nrow(eval_pts),N),])
}

# Do mMc-PSO
bd <- c(0,1)
D <- kmeanspso(clust_pts, eval_pts, cluster_center, q, 2.0,
params_pso\$w,params_pso\$c1,params_pso\$c2,
2*npart, itmax_pso, itmax_pp, 1000, 1000, 1e-4, 1e-4,
1e-4, itmax_inn,
parallel::detectCores(), jit, bd[1], bd[2], rep(0,N))
# D\$gbes_centers <- cbind(D\$gbes_centers[,2],1-D\$gbes_centers[,1])
return(D\$gbes_centers)
}
```

## Try the minimaxdesign package in your browser

Any scripts or data that you put into this service are public.

minimaxdesign documentation built on July 13, 2021, 1:06 a.m.