---
title: "Vignette for Package LHD"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Vignette for Package LHD}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


This vignette gives a high level overview on how to use the R package LHD to generate different types of efficient Latin hypercube designs (LHDs) with flexible sizes. Popular optimality criteria along with other useful functions will be introduced as well.

\

Load the library first.
```{r}
devtools::load_all()
library(LHD)
```

The overview starts with basic functions. The first function to be introduced is the "rLHD", which generates a random LHD with dimension $n$ by $k$. Suppose we want a random LHD with 6 rows and 3 columns, and denote it as X. We can run the following code.

```{r}
set.seed(1)
X=rLHD(n=6,k=3);X
```

If we want to know the inter-site distance between row $i$ and row $j$ of X ($i \neq j$), we could use function "dij". For example, the inter-site distance of the 2nd and the 4th row of X can be got via the following code.

```{r}
dij(X=X,i=2,j=4) #The default setting is the rectangular distance.

#If Euclidean distance is desired, run the following
dij(X=X,i=2,j=4,q=2)
```

The formula of "dij" is given by the following

\begin{equation}
d(\boldsymbol{x_i}, \boldsymbol{x_j})=  (\sum_{l=1}^{k}|x_{il}-x_{jl}|^q)^{1/q},
\end{equation}

where $\boldsymbol{x_i}$ and $\boldsymbol{x_j}$ denote two design points from LHD matrix X. Having all the inter-site distances calculated, we could further calculate the $\phi_p$ criterion of X according to the formula from Jin, Chen and Sudjianto (2005).

\begin{equation}
\phi_{p}= \bigg\{\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}d(\boldsymbol{x_i}, \boldsymbol{x_j})^{-p} \bigg\} ^{1/p}.
\end{equation}

The function "phi_p" provides calculation for above criterion. We can run the following code.

```{r}
phi_p(X=X,p=15)  

#If Euclidean distance is desired, run the following
phi_p(X=X,p=15,q=2)  
```

$\phi_p$ criterion is popular for maximin distance LHDs. When orthogonal or nearly orthogonal LHDs are desired, average absolute correlation and maximum absolute correlation are widely used (Georgiou (2009)). Their formula are given by the following

\begin{equation}
ave(|q|) = \frac{2 \sum_{i=1}^{k-1} \sum_{j=i+1}^{k}|q_{ij}|}{k(k-1)}, \, \, max(|q|)= max_{ij} |q_{ij}|,
\end{equation}

where $q_{ij}$ is the columnwise correlation. The function "AvgAbsCor" and "MaxAbsCor" provide calculations for above criteria. We can run the following code.

```{r}
AvgAbsCor(X=X)  #The average absolute correlation of X

MaxAbsCor(X=X)  #The maximum absolute correlation of X
```

Another optimality criterion to be introduced is the maximum projection criterion (Joseph, Gul and Ba (2015)), which focuses on the projection property of the design matrix, and its formula is given by the following

\begin{equation}
\psi = \Bigg\{ \frac{1}{{n \choose 2}} \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \frac{1}{\Pi_{l=1}^{k}(x_{il}-x_{jl})^2}  \Bigg\}^{1/k}.
\end{equation}

The function "MaxProCriterion" provides calculation for above criterion. We can run the following code.

```{r}
MaxProCriterion(X=X)  #The maximum projection criterion of X
```

Usually, a random LHD does not guarantee good space-filling property, orthogonality, nor full projection property. Morris and Mitchell (1995) proposed a version of the simulated annealing algorithm (SA) for generating maximin distance LHDs. Our function "SA" implements their algorithm, and we also added OC (optimality criterion) input argument in the "SA" function. Therefore, not only does "SA" generate maximin distance LHDs, but also it generates orthogonal and maximum projection LHDs. This is also true for the other algorithms in our package.

The function "exchange" makes "SA" workable since it is capable of switching two random elements from a given column of X. The following code shows examples of both "exchange" and "SA".

```{r}
#Suppose we want to exchange two random elements from the 1st column of X.
Xnew=exchange(X=X,j=1)

#Look and compare
X;Xnew

#The SA function with default setting
set.seed(1)
trySA=SA(n=6,k=3)
trySA

#The phi_p of trySA is much smaller than the phi_p of X, which was 0.2517586
phi_p(X=trySA,p=15) 

#Now we try SA function with a different optimality criterion
set.seed(1)
trySA2=SA(n=6,k=3,OC="AvgAbsCor")
AvgAbsCor(trySA2)      #The average absolute correlation is about 0.05
```

Leary, Bhaskar and Keane (2003) modified the work of Morris and Mitchell (1995). They showed that orthogonal-array-based LHD (OALHD) with SA converges faster and generates better designs. Our function "OASA" implements their algorithm. 

The function "OA2LHD" makes "OASA" workable since it is capable of transferring an Orthogonal Array (OA) into an LHD with corresponding size, based on the work of Tang (1993). The following code shows examples of both "OA2LHD" and "OASA".

```{r}
#create an OA(9,2,3,2)
OA=matrix(c(rep(1:3,each=3),rep(1:3,times=3)),ncol=2,nrow=9,byrow = F);OA

#Transfer the OA above into a 9 by 2 LHD
tryOA=OA2LHD(OA=OA)
tryOA

#The OASA function
set.seed(1)
tryOASA=OASA(OA=OA)
tryOASA

#The phi_p of tryOASA is smaller than the phi_p of SA 
phi_p(X=tryOASA,p=15); phi_p(X=SA(n=9,k=2),p=15) 

#Now we try OASA function with a different optimality criterion
tryOASA2=OASA(OA=OA,OC="AvgAbsCor")
AvgAbsCor(tryOASA2)      #The average absolute correlation is 0
```

Joseph and Hung (2008) modified the work of Morris and Mitchell (1995) in another way that is able to reduce both $\phi_{p}$ and column-wise correlations at the same time. Our function "SA2008" implements their algorithm. See the following code and example.

```{r}
set.seed(1)
trySA2008_63=SA2008(n=6,k=3)
trySA2008_63
phi_p(X=trySA2008_63,p=15) 

#Another example with different n and k
trySA2008_92=SA2008(n=9,k=2)
trySA2008_92
phi_p(X=trySA2008_92,p=15) 

#Now we try SA2008 function with a different optimality criterion
set.seed(1)
trySA2008=SA2008(n=6,k=3,OC="AvgAbsCor")
AvgAbsCor(trySA2008)      #The average absolute correlation is about 0.03
```

Ba, Myers and Brenneman (2015) extended the idea of Sliced Latin Hypercube Designs (SLHD) from Qian (2012) and had their own modifications of SA from Morris and Mitchell (1995). They called their algorithm "improved two-stage algorithm". Our function "SLHD" implements their algorithm. See the following code and example.

```{r}
set.seed(1)
trySLHD_63F=SLHD(n=6,k=3,t=2)   #The default setting (stage2 is FALSE)
trySLHD_63F
phi_p(X=trySLHD_63F,p=15) 

#If the second stage is desired, run the following
trySLHD_63T=SLHD(n=6,k=3,t=2,stage2=TRUE)
trySLHD_63T
phi_p(X=trySLHD_63T,p=15)       #Result is improved (phi_p is smaller than above)

#Now we try SLHD function with a different optimality criterion
set.seed(1)
trySLHD=SLHD(n=6,k=3,OC="AvgAbsCor",stage2=TRUE)
AvgAbsCor(trySLHD)      #The average absolute correlation is about 0.07
```

Chen et al. (2013) followed the structure of Particle Swarm Optimization (PSO) framework and proposed a version of it for constructing maximin distance LHDs, and their algorithm is called LaPSO. Our function "LaPSO" implements their algorithm. See the following code and example.

```{r}
set.seed(1)
tryLaPSO_63=LaPSO(n=6,k=3)
tryLaPSO_63
phi_p(X=tryLaPSO_63,p=15) 

#Another example with different n and k
tryLaPSO_92=LaPSO(n=9,k=2)
tryLaPSO_92
phi_p(X=tryLaPSO_92,p=15) 

#Now we try LaPSO function with a different optimality criterion
set.seed(1)
tryLaPSO=LaPSO(n=6,k=3,OC="AvgAbsCor")
AvgAbsCor(tryLaPSO)      #The average absolute correlation is about 0.1
```

Liefvendahl and Stocki (2006) proposed a version of the genetic algorithm (GA) which implemented an elite strategy to focus on the global best directly for constructing maximin distance LHDs. Our function "GA" implements their algorithm. See the following code and example.

```{r}
set.seed(1)
tryGA_63=GA(n=6,k=3)
tryGA_63
phi_p(X=tryGA_63,p=15) 

#Another example with different n and k.
tryGA_92=GA(n=9,k=2)
tryGA_92
phi_p(X=tryGA_92,p=15) 

#Now we try GA function with a different optimality criterion
set.seed(1)
tryGA=GA(n=6,k=3,OC="AvgAbsCor")
AvgAbsCor(tryGA)      #The average absolute correlation is about 0.07
```

In addition to above algorithms, Wang, L., Xiao, Q., and Xu, H. (2018) proposed a construction method for maximin $L_1$ distance LHDs based on good lattice point designs. One big advantage is that their method gives instant solutions, and this is very helpful for practitioners who seek immediate results. See the following code and example.

```{r}
#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".
```

Besides maximin distance LHDs, orthogonal Latin hypercube designs (OLHDs) are popular as well. Ye (1998) proposed a construction method for generating OHLDs. Our function "OLHD.Y1998" implements this method. See the following code and example.

```{r}
#create an orthogonal LHD with m=4. So n=2^m+1=17 and k=2*m-2=6
tryOLHD=OLHD.Y1998(m=4)
tryOLHD
MaxAbsCor(tryOLHD)  #zero columnwise correlations
```

Cioppa and Lucas (2007) extended Ye's idea, and their method is able to accommodate more factors with the same run size. Our function "OLHD.C2007" implements this method. See the following code and example. 

```{r}
#create an orthogonal LHD with m=4. So n=2^m+1=17 and k=4+3=7
tryOLHD=OLHD.C2007(m=4)
tryOLHD
MaxAbsCor(tryOLHD)  #zero columnwise correlations
```

Lin et al. (2009) proposed to couple OLHDs with OAs to generate OLHDs with large factor size. Our function "OLHD.L2009" implements this method. See the following code and example. 

```{r}
#create a 5 by 2 OLHD
OLHD=OLHD.C2007(m=2)

#create an OA(25,6,5,2)
OA=matrix(c(2,2,2,2,2,1,2,1,5,4,3,5,3,2,1,5,4,5,1,5,4,3,2,5,
4,1,3,5,2,3,1,2,3,4,5,2,1,3,5,2,4,3,1,1,1,1,1,1,4,3,2,1,5,5,
5,5,5,5,5,1,4,4,4,4,4,1,3,1,4,2,5,4,3,3,3,3,3,1,3,5,2,4,1,3,
3,4,5,1,2,2,5,4,3,2,1,5,2,3,4,5,1,2,2,5,3,1,4,4,1,4,2,5,3,4,
4,2,5,3,1,4,2,4,1,3,5,3,5,3,1,4,2,4,5,2,4,1,3,3,5,1,2,3,4,2,
4,5,1,2,3,2),ncol=6,nrow=25,byrow=T)

#Construct a 25 by 12 OLHD
tryOLHD=OLHD.L2009(OLHD,OA)
tryOLHD
MaxAbsCor(tryOLHD)  #zero columnwise correlations
```

Sun et al. (2010) proposed to construct OLHDs with flexible run sizes via matrices manipulations. Our function "OLHD.S2010" implements this method. See the following code and example. 

```{r}
#create an orthogonal LHD with C=3, r=3, type="odd".
#So n=3*2^4+1=49 and k=2^3=8
tryOLHD1=OLHD.S2010(C=3,r=3,type="odd")
tryOLHD1

#create an orthogonal LHD with C=3, r=3, type="even".
#So n=3*2^4=48 and k=2^3=8
tryOLHD2=OLHD.S2010(C=3,r=3,type="even")
tryOLHD2

MaxAbsCor(tryOLHD1)  #zero columnwise correlations
MaxAbsCor(tryOLHD2)  #zero columnwise correlations
```

Butler (2001) proposed to construct OLHDs using Williams transformation. Our function "OLHD.B2001" implements this method. See the following code and example. 

```{r}
#create an orthogonal LHD with n=11 and k=5
OLHD.B2001(n=11,k=5)

#create an orthogonal LHD with n=7 and k=6
OLHD.B2001(n=7,k=6)
```
