Stuttgart_May_03_2017.R

In nhavt/Block3AugDesigns:

#' Testing update formulae
#' 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)
library(Matrix)
devtools::load_all()

nTreatments <- 39
nChecks     <- 3
nRows       <- 6
nCols       <- 9
repli       <- -1
flag        <- 0
rowgroup    <- 1
colgroup    <- 1
subrgroup   <- c(1) # beta
subcgroup   <- c(1) # alpha
vrepChecks  <- c(6,6,6)
totalchecks <- sum(vrepChecks)
fullchecks  <- rep(c(1:nChecks),vrepChecks)

loops        <- 1
optDesignarr <- rep(0,loops)
Designarr    <- array(list(), loops)
randomcheck  <- 1
#mlist       <- startDesignV2(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks)

while(randomcheck<=loops){
  mlist       <- startDesignV2(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks)
  nDesign     <- mlist$nDesign
  optDesign   <- mlist$optDesign
  copyDesign <- nDesign
  inform1    <- info2Matrix(copyDesign, nTreatments, nRows, nCols, 1)
  niter <-loops
  for (i in 1:(nCols-1)){
     for(j in (i+1):nCols){
     iter <- 0
     while (iter < niter){
       irow1 <- sample(1:nRows, 1, replace = F)
       irow2 <- sample(1:nRows, 1, replace = F)
       #temp    <- aRCDesign(iPDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom = 1)
       if(copyDesign[irow1, i]!=copyDesign[irow2, j]){
           check    <- updateFormulaeT(copyDesign, inform1, irow1, i, irow2,j)
           if (check$Aopt>optDesign){
             optDesign  <- check$Aopt
             # We need to swap treatments
    		 iPDesign <- copyDesign
             xx <- iPDesign[irow1, i]
             yy <- iPDesign[irow2, j]
             iPDesign[irow1, i] <- yy
             iPDesign[irow2, j] <- xx
             copyDesign <- iPDesign
             inform1    <- check
          }
           iter <- iter + 1
        }
      }
     }
  }
  #
    for (i in 1:(nRows-1)){
     for(j in (i+1):nRows){
     iter <- 0
     while (iter < niter){
       icol1 <- sample(1:nCols, 1, replace = F)
       icol2 <- sample(1:nCols, 1, replace = F)
       #temp    <- aRCDesign(iPDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom = 1)
       if(copyDesign[i,icol1]!=copyDesign[j, icol1]){
           check    <- updateFormulaeT(copyDesign, inform1, i, icol1, j, icol2)
           if (check$Aopt>optDesign){
             optDesign  <- check$Aopt
             # We need to swap treatments
    		     iPDesign <- copyDesign
             xx <- iPDesign[i, icol1]
             yy <- iPDesign[j, icol2]
             iPDesign[i, icol1] <- yy
             iPDesign[j, icol2] <- xx
             copyDesign <- iPDesign
             inform1    <- check
          }
           iter <- iter + 1
        }
      }
     }
  }
  optDesignarr[randomcheck] <- optDesign
  Designarr[[randomcheck]]  <- copyDesign
  randomcheck <- randomcheck+1
  print(randomcheck)
}
cat("\f")

#t    <- aRCDesign(copyDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom = 1)

#$inform1 <- info2Matrix(copyDesign, nTreatments, nRows, nCols, 1)
#t2 <- updateFormulaeT(copyDesign, inform1, irow1, i, irow2,j)
#t       <- aRCDesign(iPDesign, nTreatments, nRows, nCols, repli, ncheck=nChecks, flag=0, cRandom = 1)
cat("\f")
max(optDesignarr)

ptm <- proc.time()

inform1 <- info2Matrix(copyDesign, nTreatments, nRows, nCols, 1)
print(proc.time() - ptm)
iPDesign <- copyDesign
iPDesign[irow1, i] <- copyDesign[irow2, j]
iPDesign[irow2, j] <- copyDesign[irow1, i]
# construct an alpha matrix
# for row-component designs
xx <- copyDesign[irow1, i]
yy <- copyDesign[irow2, j]
alpha1 <- matrix(0,nTreatments,nRows)
alpha1[xx,irow1] <- alpha1[xx,irow1] - 1
alpha1[xx,irow2] <- alpha1[xx,irow2] + 1
alpha1[yy,irow1] <- alpha1[yy,irow1] + 1
alpha1[yy,irow2] <- alpha1[yy,irow2] - 1

# for column-component designs
alpha2 <- matrix(0,nTreatments,nCols)
alpha2[xx,i] <- alpha2[xx,i] - 1
alpha2[xx,j] <- alpha2[xx,j] + 1
alpha2[yy,i] <- alpha2[yy,i] + 1
alpha2[yy,j] <- alpha2[yy,j] - 1

inform2 <- info2Matrix(iPDesign, nTreatments, nRows, nCols,1)
test    <- (1/nCols)*(inform2$YY %*% t(inform2$YY)-inform1$YY%*%t(inform1$YY)) +
  (1/nRows)*(inform2$ZZ%*% t(inform2$ZZ)-inform1$ZZ%*% t(inform1$ZZ))

rtt   <- inform1$YY %*% t(alpha1) + alpha1 %*%t(inform1$YY) + alpha1 %*% t(alpha1)
ctt   <- inform1$ZZ %*% t(alpha2) + alpha2 %*%t(inform1$ZZ) + alpha2 %*% t(alpha2)
final <- (1/nCols)*rtt+(1/nRows)*ctt
rdeltapowerinv <- inform1$rdeltapowerinv
#final <- rdeltapowerinv %*% final %*%rdeltapowerinv
selColumn <- sort(setdiff(c(1:nTreatments), c(xx,yy)))
newdesign <- final[,c(xx,yy,selColumn)]
newdesign <- newdesign[c(xx,yy,selColumn),]
ordervec  <- c(xx,yy,selColumn)
newordervec <- c(which(ordervec==1),which(ordervec==2),which(ordervec==3),which(ordervec==4),which(ordervec==5))

tz <- newdesign[3:nTreatments,1]
a1 <- -newdesign[1,2]
a2 <- newdesign[1,1] + newdesign[1,2]

lamda1 <- a1 + sqrt(a1*a1+a2*a2+2*t(tz)%*%tz)
lamda2 <- a1 - sqrt(a1*a1+a2*a2+2*t(tz)%*%tz)

q11 <- (a1*lamda1+ t(tz)%*%tz)/(lamda1-a2)
q21 <- q11 - lamda1


q12 <- (a1*lamda2+ t(tz)%*%tz)/(lamda2-a2)
q22 <- q12 - lamda2

p1     <- c(q11,q21,t(tz))
normp1 <- sqrt(1/(p1%*%p1))*p1

p2  <- c(q12,q22,t(tz))
normp2 <- sqrt(1/(p2%*%p2))*p2

#------------------------------------------------------------------#
InvCO1  <- inform1$InvCn
zz      <- c(lamda1,lamda2)
omega   <- diag(zz)
X       <- cbind(matrix(normp1),matrix(normp2))
X       <- X[newordervec,]
siM     <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM  <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
ginv(rdeltapowerinv)%*%checkInvCO2%*%ginv(rdeltapowerinv)

#Echeck      <- (nTrts-1)/sum(diag(checkInvCO2))