#' 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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.