Stuttgart_April_12_2017.R

In nhavt/Block3AugDesigns:

#' Stuttgart 05 April 2017
#' @author Nha Vo-Thanh
#' reading data from text files
#'
rm(list = ls())
setwd(
  "D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/UpdateFormular/Block3AugDesigns/"
)
#Sys.setenv(JAVA_HOME='C:\\Program Files\\Java\\jre1.8.0_121')
#install.packages("rJava")
library("rJava")
library(xlsx)

library(devtools)
library(roxygen2)
library(MASS)
devtools::load_all()

#-------------------------------------#
source(
  "D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/A_Eff_RowColumnDesign.r"
)
source("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/A_Eff_BlockDesign.r")
nTreatments <- 5
nChecks     <- 2
nRows       <- 3
nCols       <- 3
repli       <- -1
flag        <- 0
repliz    <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- c(3 , 3)
pDesign  <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)

check <- AOpt(pDesign, nTreatments, nRows, nCols , repli, flag)
t     <- 0
ptm   <- proc.time()
iter  <- 0
niter <- 100
olditer <- 0

rowgroup    <- 1
colgroup    <- 1
#nRows       <- 9
#nCols       <- 18
repli       <- -1
ncheck      <- 2
flag        <- 0
subrgroup   <- c(3) # beta
subcgroup   <- c(3) # alpha

eff <- 0
while (t == 0) {
  iter <- iter +1
  for (ii in 1:nRows) {
    indexx <- sample(1:nCols, nCols, replace = F)
    pDesign[ii, ] <-  pDesign[ii, indexx]
  }
  
  t   <- aRCDesign(pDesign, nTreatments, nRows, nCols, repli, nChecks, 0, cRandom = 1)
  #eff <- A3BlkDesign(pDesign, nTreatments, nRows=22, nCols=10, rowgroup, colgroup, subrgroup, subcgroup, repli, ncheck, flag, cRandom=1, reps=1)
  if((iter==niter) && (t==0)){
     repliz    <- rep(1, nTreatments)
     randIndex <- sample(1:nTreatments, nChecks, replace = F)
     repliz[randIndex] <- c(3 , 3)
     pDesign  <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
     iter     <- 0
     olditer  <- niter
     
     print("The first iterations")
  }
}
cputime <- proc.time() - ptm
print(cputime)
#cat("\014")
#check <- AOpt(t(pDesign), nTreatments, nCols , nRows, repli, flag)


#adding observation, and deleting observation

nTreatments <- nTreatments
nRows       <- nRows
nCols       <- nCols
irow        <- 1
icol        <- 9
itrt        <- 1

# information matrix of an orginal design

infoMat <- infoMatrix(pDesign,nTreatments, nRows, nCols)
Cn      <- infoMat$Cn
AnBn    <- infoMat$AnBn
InvCO   <- infoMat$InvCn

#deleting a treatment#
#dlist  <- delInfoMatrix(pDesign,nTreatments, nRows, nCols, irow, icol)
#alist  <- addInfoMatrix(pDesign,nTreatments, nRows, nCols, irow, icol, itrt)

#write.xlsx(InvCO, file="./mydata.xlsx", sheetName="Orginal")
#write.xlsx(alist$InvCn, file="./mydata.xlsx", sheetName="Update", append=TRUE)

# Swapping two treatments
# The first treatment
irow1 <- 2
icol1 <- 1
iTrt1 <- pDesign[irow1,icol1]
# The second treatment
irow2 <- 2
icol2 <- 3
iTrt2 <- pDesign[irow2,icol2]

nDesign <- pDesign
nDesign[irow1,icol1] <- iTrt2
nDesign[irow2,icol2] <- iTrt1

infoMat1 <- infoMatrix(pDesign,nTreatments, nRows, nCols)
infoMat2 <- infoMatrix(nDesign,nTreatments, nRows, nCols)

Cn1      <- infoMat1$Cn
AnBn1    <- infoMat1$AnBn
InvCO1   <- infoMat1$InvCn

Cn2      <- infoMat2$Cn
AnBn2    <- infoMat2$AnBn
InvCO2   <- infoMat2$InvCn

deltaM1 <- Cn1-Cn2
zz     <- round(eigen(deltaM1)$values*10^10)/10^10
zz     <- zz[which(zz!=0)]
omega  <- diag(zz)

X      <- eigen(deltaM1)$vectors[,c(1,5)]

deltaM2<- X %*%omega %*%t(X)

siM    <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)

checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1

# k plots
# b blocks
k <-3
b <-3  
  
Kb <- matrix(1/b,b,b)
Kk <- matrix(1/k,k,k)
Ib <- diag(b)
Ik <- diag(k)


C0 <- kronecker(Kb,Kk)
C1 <- kronecker(Ib-Kb,Kk)
C2 <- kronecker(Ib,Ik-Kk)