#' Develop by Nha Vo-Thanh
#' Stuttgart on 22 March 2017
#' Augmented row-column designs
#' Version 1.0
#'
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2Matrix <- function(nTreatments, nRows, nCols, A)
{
trtModel <- matrix(0, nRows * nCols, nTreatments)
for (i in 1:nRows) {
for (j in 1:nCols) {
# set level for treatments in a model matrix
for (trt in 1:nTreatments) {
if (A[i, j] == trt) {
trtModel[(i - 1) * nCols + j, trt] = 1
}
}
}
}
## a model for rows
init <- rep(1,nCols)
mlist <- rep(list(init),nRows)
rModel <- as.matrix(bdiag(mlist))
## a model for columns
initCol <- rep(1,nCols)
initRow <- rep(1,nRows)
initblock <- diag(initCol)
cModel <- kronecker(initRow, initblock)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
#' create a model matrix
#' @param nTreatments number of treatments
#' @param nRows number of rows
#' @param nCols number of columns
#' @return a model matrix
#' @export
model2MatrixV2 <- function(nTreatments, nRows, nCols, A)
{
trtModel <- matrix(0, nRows * nCols, nTreatments)
vecA <- c(t(A))
mylist <- rep(list(vecA),nTreatments)
fxn <- function(var1,var2){
as.numeric(var1==var2)
}
var2 <- as.list(c(1:nTreatments))
trtModel <- mapply(fxn,mylist,var2)
## a model for rows
init <- rep(1,nCols)
mlist <- rep(list(init),nRows)
rModel <- as.matrix(bdiag(mlist))
## a model for columns
initCol <- rep(1,nCols)
initRow <- rep(1,nRows)
initblock <- diag(initCol)
cModel <- kronecker(initRow, initblock)
model <- list(trtModel,rModel,cModel)
names(model) <- c("X","R","C")
return(model)
}
##################################
# return an information matrix
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param flag a flag to decide whether to include inverse matrix
#' @return information matrix
#' @export
info2Matrix <- function(Design,nTreatments, nRows, nCols, flag=0){
# connectedFlag
connectedFlag <- 0
#modelA <- model2Matrix(nTreatments, nRows, nCols, Design)
modelA <- model2MatrixV2(nTreatments, nRows, nCols, Design)
# information for treatments
X <- modelA$X
# information for rows and columns
R <- modelA$R
C <- modelA$C
# tranpose matrix
tX <- t(X)
tR <- t(R)
tC <- t(C)
# incidence matrix
NR <- tX %*% R
NC <- tX %*% C
# tranpose matrix
tNR <- t(NR)
tNC <- t(NC)
# creates a vector of ones
onevec <- matrix(1, nrow = dim(tX)[2], ncol = 1)
r <- tX %*% onevec
tr <- t(r)
rdelta <- tX %*% X
rdeltapower <- rdelta^(1/2)
rdeltapowerinv <- solve(rdeltapower)
# information matrix
A <- rdelta - (1/nCols) * NR %*% tNR -(1/nRows)*NC %*% tNC + (1/(nRows*nCols)) * r %*% tr
A_star <- rdeltapowerinv %*% A %*% rdeltapowerinv
AnBn <- cbind(NR,NC,r)
YY <- NR # for the row blocking factor
ZZ <- NC # for the column blocking factor
if(flag){
InvA <- ginv(A)
InvAStar <- ginv(A_star)
}else{
InvA <- 0
InvAStar <- 0
}
if(rankMatrix(A)[1]==nTreatments-1){
connectedFlag <- 1
}
infoMat <- list(A,A_star,AnBn,InvA,InvAStar,YY,ZZ,rdeltapowerinv,connectedFlag,r)
names(infoMat) <- c("Cn","CnStar","AnBn","InvCn","InvAStar","YY","ZZ","rdeltapowerinv","connectedness","repli")
return(infoMat)
}
##################################
# return an information matrix after deleting a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @return informatrix
#' @export
delInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol){
# Cn1
# get infomration for the orginal matrix
infoMat <- info2Matrix(Design,nTreatments, nRows, nCols)
Cn <- infoMat$Cn
AnBn <- infoMat$AnBn
InvCn <- infoMat$InvCn
# get new vectors
a1 <- rep(0,nTreatments)
a1[Design[irow,icol]] <- 1
b1 <- rep(0,nRows+nCols+1)
b1[irow] <- 1
b1[nRows+icol] <- 1
b1[nRows+nCols+1]<- 1
Ik <- (1/nCols) * diag(x = 1, nRows, nRows)
Ib <- (1/nRows) * diag(x = 1, nCols, nCols)
InvB <- bdiag(Ik,Ib,-1/(nCols*nRows))
#
xxx1 <- a1-AnBn %*% InvB %*% b1
xxx2 <- sqrt(1-t(b1) %*% InvB %*% b1)
c1 <- as.matrix(xxx1/xxx2[1])
Cn1 <- Cn - c1 %*% t(c1)
InvCn1 <- InvCn + (InvCn %*% c1 %*% t(c1) %*% InvCn)/(1-t(c1) %*% InvCn %*% c1)[1]
mlist <- list(Cn,Cn1,a1,b1,InvB,InvCn1)
names(mlist) <- c("Cn","Cn1","a1","b1","InvB","InvCn1")
return(mlist)
}
##################################
# return an information matrix after adding a treatment in irow and icol
# input :
# output:
# Version 1.1
##################################
#' @param Design a given design
#' @param nTreatments the number of treatments
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @param irow the element in a row will be deleted
#' @param icol the element in a column will be deleted
#' @param itrt the treatment will be deleted
#' @return informatrix
#' @export
addInfoMatrix <- function(Design,nTreatments, nRows, nCols, irow, icol, itrt){
mlist <- delInfoMatrix(Design,nTreatments, nRows, nCols, irow, icol)
Cn1 <- mlist$Cn1
a1 <- mlist$a1
b1 <- mlist$b1
InvB <- as.matrix(mlist$InvB)
InvCn1 <- mlist$InvCn1
a2 <- rep(0,nTreatments)
a2[itrt]<- 1
b2 <- rep(0,nRows+nCols+1)
b2[irow] <- 1
b2[nRows+icol] <- 1
b2[nRows+nCols+1]<- 1
An_1Bn_1 <- AnBn-a1%*%t(b1)
InvB_1 <- InvB+ (InvB %*% b1 %*% t(b1) %*% InvB)/(1-t(b1) %*% InvB %*% b1)[1]
xxx1 <- a2-An_1Bn_1 %*% InvB_1 %*% b2
xxx2 <- sqrt(1+t(b2) %*% InvB_1 %*% b2)
c2 <- xxx1/xxx2[1]
Cn <- Cn1 + c2 %*% t(c2)
#InvCn1 <- ginv(Cn1)
InvCn <- InvCn1 - (InvCn1 %*% c2 %*% t(c2) %*% InvCn1)/(1+t(c2) %*% InvCn1 %*% c2)[1]
mlist <- list(Cn,Cn1,c2,a2,b2,InvB_1,InvCn)
names(mlist) <- c("Cn","Cn1","c2","a2","b2","InvB_1","InvCn")
return(mlist)
}
##################################
# Searching optimal row-column designs using simulated annealing algorithm
# input :
# output:
# Version 1.1
# Stuttgart 11 April 2017
##################################
#' @param nTreatments the number of treatments
#' @param nCheks the number of checks
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @return return an optimal design
#' @export
startDesignV1 <- function(nTreatments, nChecks, nRows, nCols, repli){
# Construct an incidence matrix N = X' %*% R
# Starting with a design
startD <- matrix(0, nRows, nCols)
# Allocate the different treatments to one block at a time
nObs <- nRows*nCols
newObs <- nObs
N <- matrix(0,nTreatments,nRows)
newrepli <-repli
listtrt <- 1:nTreatments
for(ib in 1:nRows){
kn <- nCols/newObs
for (i in 1:nTreatments){
N[i,ib] <- floor(newrepli[i]*kn)
}
rfill <- nCols-sum(N[,ib])
apha <- newrepli-N[,ib]
ixx <- listtrt[apha>0]
selTrts <- sample(ixx, rfill, replace=F)
N[selTrts,ib] <- N[selTrts,ib]+1
newObs <- newObs-nCols
newrepli <- newrepli-N[,ib]
}
##--- Constructing a design from a given incidence matrix
# N = t(X) %*% Z
startD <- matrix(0, nRows, nCols)
for (ir in 1:nRows){
startD[ir,] <- rep(listtrt,N[,ir])
}
return(startD)
}
#' The function helps to generate a starting connected design
#' @param nTreatments the number of treatments
#' @param nCheks the number of checks
#' @param nRows the number of rows
#' @param nCols the number of columns
#' @return return an optimal design
#' @export
startDesignV2<- function(nTreatments, nChecks, nRows, nCols, repli, nchecks, vrepChecks, weighted=0, flag){
# Construct an incidence matrix N = X' %*% R
# Starting with a design
#--------------------------------------------------------#
nTrts <- nTreatments-nChecks+1
repliz <- rep(1, nTrts)
repliz[1] <- sum(vrepChecks)
totalchecks <- sum(vrepChecks)
fullchecks <- rep(c(1:nChecks),vrepChecks)
optDesign <- 0
checksM <- 0
while((optDesign<=0)&&(checksM<=101)){
checksM <- 0
while(optDesign<=0){
pDesign <- startDesignV1(nTrts, nChecks, nRows, nCols, repliz)
#write.xlsx(pDesign, "pDesign1.xlsx",sheetName = "Step1")
for(i in 1:nRows){
indexx <- sample(1:nCols, nCols, replace = F)
pDesign[i,] <-pDesign[i,indexx]
}
optDesign <- aRCDesign(pDesign, nTrts, nRows, nCols, repliz, ncheck=1, flag=0, cRandom=0, reps=1,weighted)
}
#write.xlsx(pDesign, "pDesign2.xlsx",sheetName = "Step2")
nDesign <- pDesign
nDesign[which(nDesign==1)] <- rep(-3,sum(vrepChecks))
nDesign <- nDesign + (nChecks-1)
xxxx <- (nChecks-1) + -3
repchecks <- fullchecks
indexx <- sample(1:totalchecks, totalchecks, replace = F)
repchecks <- repchecks[indexx]
nDesign[which(nDesign==xxxx)] <- repchecks
optDesign <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
# randomly assign checks to qseudo check
while((optDesign<=0)&&(checksM<=100)){
checksM <- checksM +1
repchecks <- fullchecks
indexx <- sample(1:totalchecks, totalchecks, replace = F)
repchecks <- repchecks[indexx]
nDesign[which(nDesign==-1)] <- repchecks
optDesign <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=0, cRandom=0, reps=1,weighted)
}
#write.xlsx(nDesign, "pDesign3.xlsx",sheetName = "Step3")
}
optDesign <- aRCDesign(nDesign, nTreatments, nRows, nCols, repliz, ncheck=nChecks, flag=flag, cRandom=0, reps=1,weighted)
mlist <- list(nDesign,optDesign)
names(mlist) <- c("nDesign","optDesign")
return(mlist)
}
#' @param nTreatments
#' @param nRows
#' @param nCols
#' @param repliz
#' @param flag type of design, CRD or row-column design, or ...
#' @return a connected design
#' @export
genBlkDesign <- function(nTreatments, nRows, nCols, nChecks, vrepChecks, flag=0){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
## fixed flag for Aopt function-- can be changed
t <- AOpt(pDesign, nTreatments, nRows, nCols , repli, 0)
t <- 0
iter <- 0
niter <- 100
olditer <- 0
# Generating a connected design
t <- 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 <- AOpt(pDesign, nTreatments, nRows, nCols, repli=-1, 0)
if((iter==niter) && (t==0)){
repliz <- rep(1, nTreatments)
randIndex <- sample(1:nTreatments, nChecks, replace = F)
repliz[randIndex] <- vrepChecks
pDesign <- StartDesign(nTreatments, nChecks, nRows, nCols, repliz)
iter <- 0
olditer <- niter
}
}
return(pDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
updateFormulaeT<- function(pDesign, infoMat, iRow1, iCol1, iRow2,iCol2, flag=0){
# get information of a given matrix
Cn1 <- infoMat$CnStar
InvCO1 <- infoMat$InvAStar
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
nTrts <- length(unique(c(pDesign)))
# a design after exchanging treatments
nDesign <- pDesign
nDesign[iRow1,iCol1] <- pDesign[iRow2,iCol2]
nDesign[iRow2,iCol2] <- pDesign[iRow1,iCol1]
infoMat2 <- info2Matrix(nDesign,nTrts, nRows, nCols)
zz <- 0
if(infoMat2$connectedness==1){
Cn2 <- infoMat2$CnStar
deltaM0 <- infoMat$Cn - infoMat2$Cn
deltaM1 <- Cn1-Cn2
if(sum(abs(deltaM1))!=0){
zz0 <- round(eigen(deltaM1)$values*10^10)/10^10
zz <- zz0[which(zz0!=0)]
if(length(zz)==2){
omega <- diag(zz)
}else{
omega <- zz
}
X <- eigen(deltaM1)$vectors[,which(zz0!=0)]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
Echeck <- (nTrts-1)/sum(diag(checkInvCO2))
}else{
Echeck <- 0
checkInvCO2 <- 0
}
}else{
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
zz <- 0
deltaM0 <- 0
}
mlist <- list(Echeck,checkInvCO2,Cn2,zz,deltaM0)
names(mlist) <- c("Aopt","InvAStar","CnStar","Eigen","deltaM0")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
Swapping<- function(pDesign,irow1,icol1,irow2,icol2){
copyDesign <- pDesign
xx <- pDesign[irow1, icol1]
yy <- pDesign[irow2, icol2]
copyDesign[irow1, icol1] <- yy
copyDesign[irow2, icol2] <- xx
return(copyDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param infoMat
#' @param nTreatments
#' @param twocombn
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @param weighted = 0
#' @param flag = 0, 'u', 'c'
#' @export
updateFormulaeT3 <- function(pDesign, infoMat, nTreatments, nChecks, twoCombn, check2Combn , iRow1, iCol1, iRow2,iCol2, weighted = 0, flag=0){
#initialization
Echeck <- 0
checkInvCO2 <- 0
Cn2 <- 0
r <- 0
YY <- 0
ZZ <- 0
x0 <- 0
final <- 0
#initialization
# get information of a given matrix
r <- infoMat$repli
r_a <- sum(r)/nTreatments
r_square <- r_a*r_a
Cn1 <- infoMat$Cn
InvCO1 <- infoMat$InvCn
nRows <- dim(pDesign)[1]
nCols <- dim(pDesign)[2]
# get two different treatments will be exchanged
xx <- pDesign[iRow1, iCol1]
yy <- pDesign[iRow2, iCol2]
# for row-component design
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 design
alpha2 <- matrix(0,nTreatments,nCols)
alpha2[xx,iCol1] <- alpha2[xx,iCol1] - 1
alpha2[xx,iCol2] <- alpha2[xx,iCol2] + 1
alpha2[yy,iCol1] <- alpha2[yy,iCol1] + 1
alpha2[yy,iCol2] <- alpha2[yy,iCol2] - 1
# calculate differences
#I need to check here
#print(sum(abs(alpha1)==0))
#print(sum(abs(alpha2)==0))
#if(sum(abs(alpha1)==0)){
#rtt <- matrix(0,nTreatments,nTreatments)
#}else{
rtt <- infoMat$YY %*% t(alpha1) + alpha1 %*%t(infoMat$YY) + alpha1 %*% t(alpha1)
#}
#if(sum(abs(alpha2)==0)){
#ctt <- matrix(0,nTreatments,nTreatments)
#}else{
ctt <- infoMat$ZZ %*% t(alpha2) + alpha2 %*%t(infoMat$ZZ) + alpha2 %*% t(alpha2)
#}
final <- (1/nCols)*rtt+(1/nRows)*ctt
YY <- infoMat$YY + alpha1 # for the row blocking factor
ZZ <- infoMat$ZZ + alpha2 # for the column blocking factor
Cn2 <- Cn1 - final
# check connectedness of a new design
diffdesign <- sum(abs(final))!=0
connected <- matrix.rank(Cn2)==(nTreatments-1)
if(diffdesign&&connected){
selColumn <- setdiff(c(1:nTreatments), c(xx,yy))
ordervec <- c(xx,yy,selColumn)
newdesign <- final[,ordervec]
newdesign <- newdesign[ordervec,]
newordervec <- matrix(0,1,nTreatments)
for(iii in 1:nTreatments){
newordervec[iii] <- c(which(ordervec==iii))
}
#get elements in a differences matrix
tz <- newdesign[3:nTreatments,1]
a1 <- -newdesign[1,2]
a2 <- newdesign[1,1] + newdesign[1,2]
tvalue <- t(tz)%*%tz
# calculate eigen values
lambda1 <- as.numeric(a1 + sqrt(a1*a1+a2*a2+2*tvalue))
lambda2 <- as.numeric(a1 - sqrt(a1*a1+a2*a2+2*tvalue))
# calculate eigen vectors
a <- newdesign[1,1]
b <- newdesign[1,2]
# There will be two different eigen vectors, which are denoted by p1 and p2
if(round(lambda1*10^10)/10^10!=0){
if((a+b+lambda1)==0){
x11 <- 0
x12 <- as.numeric((tvalue)/b)
x13 <- -t(tz)*x12/lambda1
}else{
rate1 <- (a+b-lambda1)/(a+b+lambda1)
if((lambda1-a-b*rate1)==0){
x11 <- as.numeric(sqrt((1-(tvalue))/(1+rate1*rate1)))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}else{
x11 <- as.numeric((tvalue)/(lambda1-a-b*rate1))
x12 <- as.numeric(rate1*x11)
x13 <- (t(tz)*x11-t(tz)*x12)/lambda1
}
}
p1 <- c(x11,x12,x13)
}else{
p1 <- rep(0,nTreatments)
p1[1] <- 1
}
# vector p2
if(round(lambda2*10^10)/10^10!=0){
if((a+b+lambda2)==0){
x21 <- 0
x22 <- as.numeric((tvalue)/b)
x23 <- -t(tz)*x22/lambda2
}else{
rate2 <- (a+b-lambda2)/(a+b+lambda2)
if((lambda2-a-b*rate2)==0){
x21 <- as.numeric(sqrt((1-(tvalue))/(1+rate2*rate2)))
x22 <- as.numeric(rate2*x21)
x23 <- (t(tz)*x21-t(tz)*x22)/lambda2
}else{
x21 <- (tvalue)/(lambda2-a-b*rate2)
x22 <- rate2*x21[1]
x23 <- (t(tz)*x21[1]-t(tz)*x22[1])/lambda2
}
}
p2 <- c(x21,x22,x23)
}else{
p2 <- rep(0,nTreatments)
p2[1] <- 1
}
lambda <- round(c(lambda1,lambda2)*10^10)/10^10
ieigen <- which(lambda!=0)
reallambda <- lambda[ieigen]
ifelse(length(reallambda)==2,omega <- diag(reallambda),omega <- reallambda)
# normalized vectors
normp1 <- sqrt(1/(p1%*%p1))*p1
normp2 <- sqrt(1/(p2%*%p2))*p2
#get update for an inverse matrix
InvCO1 <- infoMat$InvCn
x0 <- cbind(matrix(normp1),matrix(normp2))
X <- x0[,ieigen]
ifelse(is.null(dim(X)),X <- t(t(X)), X <- X)
#get update for an inverse matrix
X <- X[newordervec,]
siM <- solve(omega)- t(X) %*%InvCO1 %*% X
invSiM <- solve(siM)
checkInvCO2 <- InvCO1 + InvCO1 %*% X %*% invSiM %*% t(X) %*% InvCO1
# add on 23 May 2017
# calculate A-efficiency factor
temp <- 0
#for (ii in c(1:(nTreatments-1))) {
# for(jj in c((ii+1):nTreatments)){
# #v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj])
# v <- (r[ii]*r[jj]/r_square)*(checkInvCO2[ii,ii]+checkInvCO2[jj,jj]-2*checkInvCO2[ii,jj])
# temp2 <- temp2 + v
# }
#}
## incorporate with flag option
if ( weighted == 1 ) {
if (flag=="0") {
rarb <- r[twocombn[,1]]*r[twocombn[,2]]
diagElements <- diag(checkInvCO2)
AAA <- rarb*(diagElements[twocombn[,1]]+diagElements[twocombn[,2]])
BBB <- diag(as.vector(r))%*%checkInvCO2%*%diag(as.vector(r))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(flag=="u") {
nEntries <- nTreatments - nChecks
InvCO2Unrep <- checkInvCO2[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCO2Unrep))-sum(InvCO2Unrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else {
InvCO2Checks <- checkInvCO2[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
if(Echeck<0){
save(checkInvCO2,file = "INV.r")
print(InvCO2Checks)
print(solve(siM))
print(connected)
print(Echeck)
}else{
#print(InvCO2Checks)
}
}
} else {
temp <- nTreatments*sum(diag(checkInvCO2))-sum(checkInvCO2)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
} # diffdesign
mlist <- list(Echeck,checkInvCO2,Cn2,r,YY,ZZ,x0,final)
names(mlist) <- c("Aopt","InvCn","Cn","repli","YY","ZZ","Eigen","final")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
layoutAugDesign <- function(pDesign,nTreatments,nChecks){
layoutDesign <- pDesign
for(i in (nChecks+1):nTreatments){
layoutDesign[pDesign==i] <- 0
}
return(layoutDesign)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
NeighboorUpdate <- function(optDeg,optA, inforDesign, iRow1, iCol1, iRow2, iCol2){
check <- updateFormulaeT3(optDeg, inforDesign, iRow1, iCol1, iRow2, iCol2)
optNew <- check$Aopt
ap <- exp((optNew-optA)/temper)
if ((ap > runif(1, 0, 1))&&(optNew>0)) {
optDeg <- Swapping(optDeg, iRow1, iCol1, iRow2, iCol2)
optA <- optNew
inforDesign <- check
}
mlist <- list(optDeg,optA,inforDesign)
names(mlist) <- c("optDeg","optA","inforDesign")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CoordinateExchange<- function(pDesign, nTreatments, nChecks , twocombn, check2Combn, nRows, nCols, infoDesign, cOpt, pos, weighted=0, flag){
# obtain row and column positions
irow <- pos[1]
icol <- pos[2]
tempCOpt <- cOpt
tempDesign <- pDesign
tempInfoDesign <- infoDesign
colPos <- icol+1
## ----> move to right middle
while ( colPos <= nCols ) {
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos + 1
}
## ----> move to left middle
colPos <- icol-1
while(colPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, irow, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
colPos <- colPos - 1
}
## ----> move to forward
rowPos <- irow-1
while(rowPos>=1){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
}
## -----> move to backward
rowPos <- irow + 1
while(rowPos<=nRows){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, icol, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
}
# backward right
rowPos <- irow+1
colPos <- icol+1
while(((rowPos<=nRows))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos + 1
}
# forward left
rowPos <- irow-1
colPos <- icol-1
while(((rowPos>=1))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos - 1
}
# forward right
rowPos <- irow-1
colPos <- icol+1
while(((rowPos>=1))&&(colPos<=nCols)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos - 1
colPos <- colPos + 1
}
# backward left
rowPos <- irow+1
colPos <- icol-1
while(((rowPos<=nRows))&&(colPos>=1)){
result <- CompareSwappingDesign(pDesign, infoDesign, tempDesign, tempCOpt, tempInfoDesign, nTreatments, nChecks , twocombn, check2Combn, irow, icol, rowPos, colPos, weighted, flag)
tempCOpt <- result$Aopt
tempDesign <- result$pDesign
tempInfoDesign <- result$infoDesign
rowPos <- rowPos + 1
colPos <- colPos - 1
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' The function helps to update an inversve of an information matrix
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
CompareSwappingDesign <- function(pDesign, infoDesign, cDesign, cOpt, cInfoDesign , nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- cDesign
tempCOpt <- cOpt
tempInfoDesign <- cInfoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag) #print(uptDesign$Aopt)
if (uptDesign$Aopt> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param pDesign
#' @param inforMatrix
#' @param iRow1
#' @param iCol1
#' @param iRow2
#' @param iCol2
#' @export
SwappingTwoPlots <- function(pDesign, infoDesign, AOpt, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted, flag){
iTrt1 <- pDesign[iRow1, iCol1]
iTrt2 <- pDesign[iRow2, iCol2]
# set variables
tempDesign <- pDesign
tempCOpt <- round(AOpt*10^10)/10^10
tempInfoDesign <- infoDesign
if(iTrt1!=iTrt2){
uptDesign <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks , twocombn, check2Combn, iRow1, iCol1, iRow2, iCol2, weighted,flag)
#print(uptDesign$Aopt)
if (round(uptDesign$Aopt*10^10)/10^10> tempCOpt) {
tempCOpt <- uptDesign$Aopt
tempDesign <- Swapping(pDesign, iRow1, iCol1, iRow2, iCol2)
tempInfoDesign <- uptDesign
}
}
mlist <- list(tempDesign,tempInfoDesign,tempCOpt)
names(mlist) <- c("pDesign","infoDesign","Aopt")
return(mlist)
}
RandomSearch <- function(pDesign,nTreatments, nRows, nCols, optDesign, infoDesign, twocombn, check2Combn, weighted, flag, niter){
copyDesign <- pDesign
iter <- 0
recordImp <- rep(0,niter)
while (iter < niter) {
cTreatments <- as.matrix(rbind(which(pDesign==1,arr.ind=T),which(pDesign==2,arr.ind=T),which(pDesign==3,arr.ind=T),which(pDesign==4,arr.ind=T)))
icheck <- sample(1:lcTreatment, 1, replace = T)
siter <- 0
while (siter < 1000) {
irow <- c(as.double(cTreatments[icheck,1]),sample(1:nRows, 1, replace = T))
icol <- c(as.double(cTreatments[icheck,2]),sample(1:nCols, 1, replace = T))
if (pDesign[irow[1], icol[1]] != pDesign[irow[2], icol[2]]) {
check <- updateFormulaeT3(pDesign, infoDesign, nTreatments, nChecks, twocombn, check2Combn, irow[1], icol[1], irow[2], icol[2], weighted,flag)
if (round(check$Aopt*10^10)/10^10 > round(optDesign*10^10)/10^10) {
#print('I am here')
optDesign <- check$Aopt
# We need to swap treatments
iPDesign <- Swapping(pDesign, irow[1], icol[1], irow[2], icol[2])
pDesign <- iPDesign
infoDesign <- check
}
}
siter <- siter + 1
}
recordImp[iter] <- optDesign
iter <- iter + 1
}
mlist <- list(optDesign,pDesign,infoDesign,recordImp)
names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(mlist)
}
#' Swapping two different plots in a given design
#' @param InvInfoMatrix
#' @param TwoCombn
#' @param Repli
#' @param Weight
#' @param Flag
#' @export
AeffRowColumnDesign <- function(InvInfoMatrix, nTreatments, nChecks, r_a, ...){
# get variables
arguments <- list(...)
TwoCombn <- arguments$TwoCombn
ifelse(!is.null(arguments$TwoCombn),TwoCombn <- arguments$TwoCombn,TwoCombn <- -1)
ifelse(!is.null(arguments$Repli),Repli <- arguments$Repli, Repli <- -1)
ifelse(!is.null(arguments$Weight),Weight <- arguments$Weight, Weight <- -1)
ifelse(!is.null(arguments$Flag),Flag <- arguments$Flag,Flag <- -1)
# set variables
v_bar <- 0
Aopt <- 0
Echeck <- 0
r_square <- r_a*r_a
if ( Weight == 1 ) {
if (Flag=="0") {
rarb <- Repli[TwoCombn[,1]]*Repli[TwoCombn[,2]]
diagElements <- diag(InvInfoMatrix)
AAA <- rarb*(diagElements[TwoCombn[,1]]+diagElements[TwoCombn[,2]])
BBB <- diag(as.vector(Repli))%*%InvInfoMatrix%*%diag(as.vector(Repli))
CCC <- sum(BBB)-sum(diag(BBB))
temp <- (sum(AAA) - CCC)/r_square
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
} else if(Flag=="u") {
nEntries <- nTreatments - nChecks
InvCOUnrep <- InvInfoMatrix[c((nChecks+1):nTreatments),c((nChecks+1):nTreatments)]
temp <- nEntries*sum(diag(InvCOUnrep))-sum(InvCOUnrep)
v_bar <- temp*2/(nEntries*(nEntries-1))
Aopt <- 2/v_bar
Echeck <- Aopt
} else if(Flag=="c") {
InvCO2Checks <- InvInfoMatrix[c(1:nChecks),c(1:nChecks)]
rCheck <- sum(r[1:nChecks])/nChecks
rSquareCheck <- rCheck*rCheck
temp <- 0
for (ii in c(1:(nChecks-1))) {
for(jj in c((ii+1):nChecks)){
v <- (r[ii]*r[jj]/rSquareCheck)*(InvCO2Checks[ii,ii]+InvCO2Checks[jj,jj]-2*InvCO2Checks[ii,jj])
temp <- temp + v
}
}
v_bar <- temp*2/(nChecks*(nChecks-1))
Aopt <- 2/(rCheck*v_bar)
Echeck <- Aopt
}
} else {
temp <- nTreatments*sum(diag(InvInfoMatrix))-sum(InvInfoMatrix)
v_bar <- temp*2/(nTreatments*(nTreatments-1))
Aopt <- 2/(r_a*v_bar)
Echeck <- Aopt
}
#mlist <- list(optDesign,pDesign,infoDesign,recordImp)
#names(mlist) <- c("optDesign","pDesign","infoDesign","recordImp")
return(Echeck)
}
##----------------------------------##
## A place to test your codes
##----------------------------------##
#setwd("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/UpdateFormular")
#source("D:/user/vthanh/Google Drive/Project Stuttgart/R_programing/A_Eff_RowColumnDesign.r")
#data <- split(read.table("./Designs/Design.txt"), gl(1, 4))
#repli <- 3
#pDesign <- data$`1`
#modelD <- modelMatrix(10, 4, 4, pDesign)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.