#' Select SMED points from a set of options
#'
#' SMED_select: R implementation
#'
#' @useDynLib SMED, .registration=TRUE
#' @importFrom Rcpp evalCpp
#' @param f Function
#' @param n Number of points to select
#' @param X0 Design matrix of already selected points (points in rows)
#' @param Xopt Matrix where each row is a candidate point
#' @param theta Correlation parameters in each dimension for distance scaling
#' @param useC Should SMED_select call SMED_selectC?
#'
#' @return Vector with indices of selected points
#' @export
#'
#' @examples
#' X0 <- matrix(c(.45,.1,.55,.9),ncol=2,byrow=TRUE)
#' Xopt <- matrix(c(.45,.5,.5,.2),ncol=2,byrow=TRUE)
#' plot(NULL, xlim=c(0,1), ylim=c(0,1))
#' points(X0)
#' text(Xopt, col=2)
#' SMED_select(function(a) {1}, 1, X0, Xopt)
#'
#' # Checking if using theta works
#' X0 <- matrix(c(.1,.1,.5,.5,.9,.9), ncol=2, byrow=TRUE)
#' Xopt <- matrix(c(.1,.9,.3,.5,.5,.3), ncol=2, byrow=TRUE)
#' plot(X0,xlim=0:1,ylim=0:1, pch=19)
#' text(Xopt, col=2, pch=19)
#' f <- function(x){1} # all have same charge, so only theta matters
#' # Each dimension equal, pick 1 since it is most spread out
#' SMED_select(f=f, n=3, X0=X0, Xopt=Xopt)
#' # y dimension has more variability, pick 3 since it is furthest from X0 in y direction
#' SMED_select(f=f, n=3, X0=X0, Xopt=Xopt, theta=c(1,1e3))
#' # x dimension has more variability, pick 2 since it is furthest from X0 in x direction
#' SMED_select(f=f, n=3, X0=X0, Xopt=Xopt, theta=c(1,1e-3))
SMED_select <- function(f,n=1, X0=NULL, Xopt=NULL, theta=rep(1,ncol(X0)), useC=F) {
if (useC) {
return(SMED_selectC(f=f, n=n, X0=X0, Xopt=Xopt, theta=theta))
}
#p <- d # dimension
p <- ncol(X0)
k <- 4*p # MED distance thing
X <- X0
n0 <- if(is.null(X0)) 0 else nrow(X0)
# Charge function qq
qq <- function(xx){f(xx)^-(1/(2*p))}
# Function we will optimize
#f_min <- function(xnew,xall,kk) {
# qq(xnew)^kk*sum(apply(xall,1,function(xx){(qq(xx)/(sqrt(sum((xx-xnew)^2))))^kk}))
#}
f_min <- function(xnew,xall, qqall,kk) {
#qq(xnew)^kk*sum(apply(xall,1,function(xx){(qq(xx)/(sqrt(sum((xx-xnew)^2))))^kk}))
#qq(xnew)^kk*sum(sapply(1:nrow(xall),function(ii){(qqall[ii]/(sqrt(sum((xall[ii,]-xnew)^2))))^kk}))
qq(xnew)^kk*sum(sapply(1:nrow(xall),function(ii){(qqall[ii]/(sqrt(sum((xall[ii,]-xnew)^2 * theta))))^kk}))
}
if (p==2) {
# Get contour plot
#my.filled.contour.func(f,nlevels=5)
#contourfilled.func(f,nlevels=5)
}
#points(X, pch=19)
Y <- apply(X,1,f)
qqX <- apply(X,1,qq)
Delta <- .01 * max(Y)
keep.Delta <- (Y > Delta)
#browser()
Xopt.inds <- 1:nrow(Xopt)
Xopt.selected <- c()
#Xopt.consider <- rep(T, nrow(X.opt)) # to only check good points when taking more than 1
# Get rest of points
for(i in 1:n) {#print(keep.Delta)
# Use log scale for optimization, had trouble before when numbers were 1e88
opt.func <- function(xx)log(f_min(xx,X[keep.Delta, , drop=F],
qqX[keep.Delta],kk=k))
func.vals <- apply(Xopt, 1, opt.func)
best <- which.min(func.vals)
xnew <- Xopt[best,]
Xopt <- Xopt[-best, , drop=F]
Xopt.selected <- c(Xopt.selected, Xopt.inds[best])
Xopt.inds <- Xopt.inds[-best]
#Xopt.consider <- Xopt.consider[-best]
X <- rbind(X,unname(xnew))
#Y <- apply(X,1,f) # could just append
Y <- c(Y, f(xnew))
#qqX <- apply(X,1,qq)
qqX <- c(qqX, qq(xnew))
# Delta is used to 'hide' points from the potential function that are not useful
# Without this the points selected were terrible
Delta <- .01 * max(Y) #.01 * (n0 / (n + i)) ^ (1 / p) * max(Y)
#keep.Delta <- (Y > Delta)
#browser()
keep.Delta <- ifelse(1:nrow(X) <= n0, Y > Delta, T)
if (p==2) {
#text(x=xnew[1],y=xnew[2],labels=i,col=1)
}
}
if (p>2) {
#pairs(X)
}
# Return design matrix
#return(X)
return(Xopt.selected)
}
if (F) {
#setwd("C:/Users/cbe117/School/DOE/SMED/SMED-Code")
# source('TestFunctions.R')
# source('C:/Users/cbe117/School/DOE/Codes/contour/contourfilled/R/contourfilled.R')
SMED(TestFunctions::banana,p=2,n=10,max.time=.2)
SMED(function(xx){xx[1]+xx[2]^2-sin(2*pi*xx[3])},p=3,n=10,max.time=.2)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.