Stuttgart_April_13_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
vrepChecks  <- c(3,3)

ptm         <- proc.time()
pDesign <- genBlkDesign(nTreatments, nRows, nCols, nChecks, vrepChecks, flag=0)
cputime <- proc.time() - ptm
print(cputime)


# copyDesign <- pDesign
# 
# t00  <- AOpt(pDesign, nTreatments, nRows, nCols, repli=-1, 0)
# t0   <- t00
# # Swapping two treatments
# # The first treatment
# niter <- 100
# Temp  <- 1.0
# T_min <- 0.00001
# alpha <- 0.9
# 
# while (Temp>=T_min){
#   i <- 1
#   while (i <= niter){
#     #new_solution = neighbor(solution)
#     #new_cost     = cost(new_solution)
#     #ap = acceptance_probability(old_cost, new_cost, T)
#     #if ap > random():
#     #  solution = new_solution
#     #old_cost = new_cost
#     #i += 1
#     #T = T*alpha
#   }
# }
# for (i in 1:(nCols-1)){
#   print(i)
#   for(j in (i+1):nCols){
#   iter <- 0
#   while (iter < niter){
#     iPDesign <- pDesign 
#     irow1 <- sample(1:nRows, 1, replace = F)
#     iTrt1 <- iPDesign[irow1,i]
#     # The second treatment
#     irow2 <- sample(1:nRows, 1, replace = F)
#     iTrt2 <- iPDesign[irow2,j]
#     iPDesign[irow1,i] <- iTrt2
#     iPDesign[irow2,j] <- iTrt1
#     t    <- AOpt(iPDesign, nTreatments, nRows, nCols, repli=-1, 0)
#     if (t>t0){
#       t0      <- t
#       pDesign <- iPDesign
#    } 
#     
#     iter <- iter + 1
#   }
#   }
# }
# t0

cat("\014")

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

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

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

Cn1      <- infoMat1$A
InvCO1   <- infoMat1$InvA

Cn2      <- infoMat2$A
InvCO2   <- infoMat2$InvA

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



XX <- matrix( 
     c(0, 0, 0, -1, 1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0), # the data elements 
     nrow=5,              # number of rows 
     ncol=3,              # number of columns 
     byrow = TRUE)        # fill matrix by rows 

Nb <- infoMat1$N
Na <- infoMat2$N
kdelta <- infoMat1$kdelta

LLL <- Nb %*% t(XX)+XX %*% t(Nb)+ XX %*% t(XX)