Nothing
#' Fast Maximin Distance LHD
#'
#' \code{FastMmLHD} returns a \code{n} by \code{k} maximin distance LHD matrix generated by the construction method of Wang, L., Xiao, Q., and Xu, H. (2018)
#'
#' @param n A positive integer, which stands for the number of rows (or run size).
#' @param k A positive integer, which stands for the number of columns (or factor size).
#' @param method A distance measure method. The default setting is "manhattan", and it could be one of the following: "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski". Any unambiguous substring can be given.
#' @param t1 A tunning parameter, which determines how many repeats will be implemented to search for the optimal design. The default is set to be 10.
#'
#' @return If all inputs are logical, then the output will be a \code{n} by \code{k} maximin distance LHD under under the maximin L_1 distance criterion..
#'
#' @references Wang, L., Xiao, Q., and Xu, H. (2018) Optimal maximin $L_{1}$-distance Latin hypercube designs based on good lattice point designs. \emph{The Annals of Statistics}, \strong{46}(6B), 3741-3766.
#'
#' @examples
#' #n by n design when 2n+1 is prime
#' try=FastMmLHD(8,8)
#' try
#' phi_p(try) #calculate the phi_p of "try".
#'
#' #n by n design when n+1 is prime
#' try2=FastMmLHD(12,12)
#' try2
#' phi_p(try2) #calculate the phi_p of "try2".
#'
#' #n by n-1 design when n is prime
#' try3=FastMmLHD(7,6)
#' try3
#' phi_p(try3) #calculate the phi_p of "try3".
#'
#' #General cases
#' try4=FastMmLHD(24,8)
#' try4
#' phi_p(try4) #calculate the phi_p of "try4".
#'
#' @export
FastMmLHD=function(n,k,method="manhattan",t1=10){
#n: number of runs
#k: number of columns
#method: the distance method, with default setting is "manhattan".
#t1: tunning parameter
if (is.numeric(n) != TRUE || !(n%%1 == 0) || (n < 3))
{
stop("'n' must be a postive integer no less than 3")
}
if (is.numeric(k) != TRUE || !(k%%1 == 0) || (k < 2))
{
stop("'k' must be a postive integer no less than 2")
}
if (n < k) {stop("n must be no less than k")}
if (!is.na(pmatch(method, "euclidian")) || (method == "l2") || (method == "L2"))
{method = "euclidean"}
if (!is.na(pmatch(method, "manhattan"))|| (method == "l1") || (method == "L1"))
{method = "manhattan"}
METHODS = c("euclidean", "manhattan")
count = pmatch(method, METHODS)
if (is.na(count)==TRUE)
stop("invalid distance method")
if (count == -1)
stop("ambiguous distance method")
###
#the first special case of n \times n design when m = 2n+1 is prime. (2.3)
if ((n == k) && numbers::isPrime(2*n + 1)==TRUE)
{
m = 2*n + 1
GD=GLP(n=m,k=m-1,h=1:(m-1)) #GLP design
RD=GD #result design
for (i in 1:m){
for (j in 1:(m-1)){
if(GD[i,j]<(m/2)){RD[i,j]=2*GD[i,j]}
if(GD[i,j]>=(m/2)){RD[i,j]=2*(m-GD[i,j])}
}
}
resultdesign = RD/2
return(resultdesign[1:n, 1:n]-1)
}
#the second special case of n \times (k<=n) design when n+1 is prime. (2.2)
else if ((n >= k) && numbers::isPrime(n + 1)==TRUE)
{
finaldesign = NULL
finalvalue = -Inf
for(i in 1:t1)
{
GD=GLP(n=n+1,k=k,h=sample(seq(from=1,to=n),k)) #starting design: D_b with b=0. GLP design
TD=WT(GD,baseline=0)[-(n+1),] #E_b with b=0. Transformed design
minvalue=min(stats::dist(TD, method = method)) #current minimum distance
for (c in 1:n) #for b=1, ... , (n+1)-1
{
newdesign = (GD + c) %% (n+1) #new D_b
newresult = WT(newdesign,baseline=0) #new E_b
cutv = newresult[(n+1),1] #cut value: this is value of the constant row
newresult=newresult[-(n+1),] #get rid of the last row since it is constant
for(i in 1:n){ #bring one number down for remaining elements > cutv
for(j in 1:k){
if(newresult[i,j] > cutv){
newresult[i,j] = newresult[i,j] - 1
}
}
}
newminvalue = min(stats::dist(newresult, method = method))
if (newminvalue > minvalue) #updating minimum distance
{
resultdesign = newresult
minvalue = newminvalue
}
}
if (minvalue > finalvalue)
{
finaldesign = resultdesign
finalvalue = minvalue
}
}
return(finaldesign)
}
#the third specical case of n \times (k<n) design when n is prime. (2.1)
else if ((n > k) && numbers::isPrime(n)==TRUE)
{
finaldesign = NULL
finalvalue = -Inf
for(i in 1:t1)
{
GD=GLP(n=n,k=k,h=sample(seq(from=1,to=(n-1)),k)) #starting design: D_b with b=0. GLP design
TD=WT(GD,baseline=0) #E_b with b=0. Transformed design
minvalue=min(stats::dist(TD, method = method)) #current minimum distance
for(c in 1:(n-1)) #for b=1, ... , n-1
{
newdesign = (GD + c) %% n #new D_b
newresult = WT(newdesign,baseline=0) #new E_b
newminvalue = min(stats::dist(newresult, method = method))
if (newminvalue > minvalue) #updating minimum distance
{
resultdesign = newresult
minvalue = newminvalue
}
}
if (minvalue > finalvalue)
{
finaldesign = resultdesign
finalvalue = minvalue
}
}
return(finaldesign)
}
else {
#General case, when above special cases are not true.
#coprime vector to n
copn=NULL
for(i in 1:(n-1)){
if (numbers::coprime(i,n)==TRUE)
{
copn=c(copn,i)
}
}
#Euler function value of n
eul=length(copn)
if (k>eul){
stop(paste0("the number of column for this run size must be no greater than ", eul))
}
else {
GD = GLP(n=n,k=eul,h=sample(copn,eul))
TD = WT(GD,baseline=0)
minvalue = min(stats::dist(TD, method = method))
for (c in 1:(n-1))
{
newdesign = (GD + c) %% n
newresult = WT(newdesign,baseline=0)
newminvalue = min(stats::dist(newresult, method = method))
if (newminvalue > minvalue)
{
resultdesign = newresult
minvalue = newminvalue
}
}
#This is the searching part: search for the best subset if k!=eul.
if(k==eul){
return(resultdesign)
}
else {
Temp=NULL
Temp.value=Inf
for (i in 1:t1) {
rcol=sample(1:eul,k) #random columns
X=resultdesign[,rcol]
X.value=phi_p(X)
if(X.value<=Temp.value){
Temp=X
Temp.value=X.value
}
}
print("Design has been generated, but it may not be the best one for this size.")
return(Temp)
}
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.