## ## Two main functions for generating designs ## find.SSP(c(n1, n2), c(w, s)) ## ## find.BSP(c(n1, n2), c(b, w, s)) ## ## ## Four functions for computing various criterion for a given design ## SSP(nfat, nUnit, p.col, b.col) p.col is the defining columns and b.col is the splitting columns ## ## SSP.wd(nfat, nUnit, wd, bkwe) wd is the defining words and bkwd is the splitting words ## ## BSP(nfat, nUnit, p.col, b.col) p.col is the defining columns and b.col is the splitting columns ## ## BSP.wd(nfat, nUnit, wd1, wd2, bkwe) wd1 is the WP defining words, wd2 is the subplot defining words ## and bkwd is the splitting words ## ## fun1=function(vec, A) sum(colSums(A[vec,])%%2) fun2=function(a) { tmp=rowSums(a) for (i in 2:nrow(a)) { out=t(combn(nrow(a), i)) b=sort(apply(out, 1, fun1, a)) tmp=c(tmp, b) } return(tmp) } fun3=function(a, Mat) { flag=1 n=nrow(Mat) k=1 while( flag==1 & k<=n ) { if (all(a==Mat[k, ])) flag=0 k=k+1 } return(flag) } ret.matrix=function(nMatrix) { #generate all defining words for designs of size 2^(nMatrix) maxCount=1 ceiling = (maxCount + 1) ^ nMatrix cv = rep(0, nMatrix) ret = matrix(NA, nMatrix, (ceiling-1)) for(counter in 1:(ceiling-1)) { val = counter for(i in 1:nMatrix) { cv[i] = val %% (maxCount + 1) val = val %/% (maxCount + 1) } ret[,counter]= cv } ret } fin.alias=function(vec1, vec) { #calculate the element-by-element sum (mod 2) of two vectors (vec1+vec)%%2 } ele.sum=function(ind, col.select) { #calculate the sum (mod 2) of col.select row-wise if (length(ind)==1) { tmp=col.select[ind,] } else { vec=as.matrix(col.select[ind,]) tmp=apply(vec, 2, sum)%%2 } tmp } fin.sum=function(ind, col.select) { #calculate the sum (mod 2) of col.select column-wise vec=as.matrix(col.select[,ind]) tmp=apply(vec, 1, sum)%%2 } cal.df=function(ind, nUnit) { #calculate the degree of freedom in each stratum #for a design with nUnit block structure #ind=length(nUnit) number of strata tmp=1 for (k in 1: ind) { if (k==ind) a=nUnit[k]-1 else a=nUnit[k] tmp=tmp*a } tmp } gen2ind=function(vec,nMatrix=5) { ret=ret.matrix(nMatrix) nbase=c(1, cumprod(rep(2, nMatrix-1))) tmp=sum(vec*nbase) tmp } #wd2index=function(wd,n1,n2,p1,p2){ # a1=n1-p1 # a2=n2-p2 # nMatrix=a1+a2 # ba=c(LETTERS[1:a1],LETTERS[15+(1:a2)]) # ind=pmatch(unlist(strsplit(wd,"")),ba) # tmp=rep(0,nMatrix) # tmp[ind]=1 # index=gen2ind(tmp, nMatrix) # return(index) #} wd2index=function(wd,n1,n2,p1,p2){ #find the corresponding index for a given independent word(s) # n1, n2: number of whole-plot and subplot factors # n1=0 is a block design where p1 is the number of independent blocking words # n2 is the number of treatment factors and p2 is the number of independent treatment words #n1>0 is a split-plot design where n1 is no of whole-plot factors and n2 is no of subplot factors # p1, p2 are the degree of fractionations in whole-plot and subplot strata # n1-p1>1, n2-p2>1 for split-plot design #wd: for block desins- a set of defining words using the Captial letters for treatment factors # for block defining words ignoring the blocking factors #wd: for split-plot designs- a set of defining words using # the Captial letters starting with A for whole-plot independing words, and # the samll letter starting with p for sub-plot independing words #use pmatch to find the corresponding 0-1 vector(s) #call gen2ind to find the corresponding index for the 0-1 vector(s) if (n1==0) { nMatrix=n2-p2 ba=c(LETTERS[1:n2]) } else { a1=n1-p1 a2=n2-p2 nMatrix=a1+a2 ba=c(LETTERS[1:a1],LETTERS[15+(1:a2)]) } ind=pmatch(unlist(strsplit(wd,"")),ba) tmp=rep(0,nMatrix) tmp[ind]=1 index=gen2ind(tmp, nMatrix) return(index) } wd2index.ssp=function(wd,nfat,nUnit){ #index for a given wd for split-plot design with few WP factors m=log2(nUnit) nMatrix=sum(m) n1=nfat[1] a2=nMatrix-n1 ba=c(LETTERS[1:n1],LETTERS[15+(1:a2)]) ind=pmatch(unlist(strsplit(wd,"")),ba) tmp=rep(0,nMatrix) tmp[ind]=1 index=gen2ind(tmp, nMatrix) return(index) } keytomi=function(key, nUnit) { N=prod(nUnit) m=log2(nUnit) n=nrow(key) nStrata=length(nUnit) mi=rep(0, (N-1)) base=c( 2^((sum(m)-1):1),1) L=rep(list(c()), nStrata+1) if (nStrata==2){ L=rep(list(c()), 3) L[[1]]=c(1:(2^m[1]-1)) L[[2]]=c(2^m[1]:(N-1)) L[[3]]=0 } else { a=0 for (depth in 1:(nStrata)) { tmp=m[depth] ind=c(cumprod(2^a):(cumprod(2^(a+tmp))-1)) L[[depth]]=ind a=a+tmp } L[[depth+1]]=0 } out=t(combn(nrow(key), 2)) col.2fi=matrix(0, ncol=sum(m), nrow=nrow(out)) for (i in 1:nrow(out)) col.2fi[i,]=colSums(key[out[i,], ])%%2 kk=col.2fi%*%base for (i in 1:(N-1)) mi[i]=sum(kk==i) mi[N]=length(kk[kk==0]) mi.vec=rep(list(c()), nStrata) L.main=key%*%base for (depth in 1:nStrata) { ind=L[[depth]] ind=ind[!(ind %in% L.main)] mi.vec[[depth]]=mi[ind] } mst=matrix(0, ncol=nStrata, nrow=2) for (i in 1: nStrata) { mst[,i]=c(sum(unlist(mi.vec[c(i:nStrata)])), sum(unlist(mi.vec[c(i:nStrata)])^2)) } list(key=key, mi.vec=mi.vec, mst=mst) } sort.order.sum=function(A) { #A[,1], A[,2],- A[,3],-A[,4] #sum s, sum s+w, sum s2 sum s2+w2 sz=dim(A)[2] cmd = paste( "order(" ) for (i in (1:(sz/2))) { if (i>1) cmd=paste( cmd, ", ", sep="") cmd = paste( cmd, "A[,", i, "]", sep="") } for (i in (1:(sz/2)+(sz/2))) { cmd=paste( cmd, ", ", sep="") cmd = paste( cmd, "-1*A[,", i, "]", sep="") } cmd = paste( cmd, ", decreasing=TRUE)") cmd_expr = parse(text=cmd) ind = eval(cmd_expr) } comp=function(vec1, vec2) { tmp=1 if (vec1[1]vec2[1]) tmp=2 if (vec1[1]==vec2[1]) { if (vec1[2]>vec2[2]) tmp=0 else if (vec1[2]1){ for (j in 2:sz) { #compare with other old ones mij=matrix(M.mi[j, c(1:(2*nStrata)) ], byrow=T, nrow=2) flag=rep(1, (k+1)) for (step in 1:(k+1)) flag[step]=comp(mst[,step], mij[,step]) if (any(flag==0)) m.tmp=c(m.tmp, j) #old is at least not worse than the new one, } M.mi<<-rbind(mcrt, M.mi[m.tmp,]) M.op<<-rbind(key.col, M.op[m.tmp,]) } } else { if (any(flag==2) & any(flag==0) ) { #not admissible m.tmp=c(1) if (sz==1) { M.mi<<- rbind(mcrt, M.mi[m.tmp, ]) M.op<<- rbind(key.col, M.op[m.tmp, ]) } else { for (j in 2:sz) { #compare with other old ones mij=matrix(M.mi[j , c(1:(2*nStrata)) ], byrow=T, nrow=2) flag=rep(1, (k+1)) for (step in 1:(k+1)) flag[step]=comp(mst[,step], mij[,step]) if (all(flag!=2)==TRUE) { M.mi<<- M.mi M.op<<- M.op break #new one is worse than one of the old ones } else { if (any(flag==0)) m.tmp=c(m.tmp, j) } if (j==sz) { M.mi<<- rbind(mcrt, M.mi[m.tmp, ]) M.op<<- rbind(key.col, M.op[m.tmp, ]) } } } } } } } } keytodefwd<- function(key, nfat,nUnit){ m=log2(nUnit) n1=nfat[1];n2=nfat[2] nMatrix=sum(m) nStrata=length(nUnit) ind=1:nMatrix nStrata=length(nUnit) BF=key[ind,] sp.ind=ind[rowSums(BF)!=1] sbkwd=list(c()) #index for splitting words for (i in 1:length(sp.ind)) { tmp=ind[!(ind %in% sp.ind)] if (nStrata==2){ flag=1 for (j in 1:length(tmp)) { cn=t(combn(length(tmp), j)) for (k in 1:nrow(cn)) { A=BF[c(cn[k,], sp.ind[i]),] sa=colSums(A)%%2 if (sum(sa[1:(nMatrix-m[1])])==0 ) { flag=0 sbkwd[[i]]=row.names(A) } if (flag==0) break } if (flag==0) break } } if (nStrata==3) { if (n1m[2] flag=1 for (j in 1:length(tmp)) { cn=t(combn(length(tmp), j)) for (k in 1:nrow(cn)) { A=BF[c(cn[k,], sp.ind[i]),] sa=colSums(A)%%2 if (sum(sa[1:(nMatrix-m[1])])==0 ) { flag=0 sbkwd[[i]]=row.names(A) } if (flag==0) break } if (flag==0) break } } } } wd=list(c()) sp.ind=(nMatrix+1):nrow(key) for (i in 1:length(sp.ind)) { flag=1 for (j in 2:nMatrix) { cn=t(combn(nMatrix, j)) for (k in 1:nrow(cn)) { A=key[c(cn[k,], sp.ind[i]),] sa=colSums(A)%%2 if (sum(sa)==0) { flag=0 wd[[i]]=row.names(A) } if (flag==0) break } if (flag==0) break } } list(sbkwd, wd) } keytowlp<-function(key){ m=ncol(key) n=nrow(key) A3=3 ret=rep(0, (n-A3+1)) for (i in A3:(n-1)){ vec=t(combn(n, i)) ind=apply(vec, 1, fun1, key)==0 ret[i-A3+1]=sum(ind) } if (sum(colSums(key)%%2)==0) ret[n-A3+1]=1 (ret) } bkwd2mi=function(nfat,nUnit, p.col, b.col){ N=prod(nUnit) nStrata=length(nUnit) nf=sum(nfat) n1=nfat[1] n2=nfat[2] nMatrix=log2(N) ret=ret.matrix(nMatrix) p=length(p.col) #no of independent treatment words b=length(b.col) #no of independent block words npcol=c(p.col, b.col) np=p+b index=c(nf,p,b) nb=2^(b) nsb=2^(nf-np) nblk=log2(nUnit) nDF=rep(NA, nStrata) for (depth in 1:nStrata) nDF[depth]=cal.df(depth, nUnit) L=2^(np)-1 if (b==0) stop("no blocking") gen= cbind( matrix(ret[, npcol], byrow=TRUE, nrow=np),diag(1,np) ) def.sub=c() for (i in 1:np){ out=combn(np,i) def.sub=rbind(def.sub, t(apply(out, 2, ele.sum, gen))) } dm_suf=matrix(0, ncol=2, nrow=L) for (i in 1:L) dm_suf[i,]=c(sum(def.sub[i,1:(nf)]), sum(def.sub[i, nf+(1:b)])) Aij=matrix(0, nrow=(nf+1), ncol=(b+1)) for (i in 1:L) { Aij[(dm_suf[i,1]+1),(dm_suf[i,2]+1)]=Aij[(dm_suf[i,1]+1),(dm_suf[i,2]+1)]+1 } Ai=Aij[,1] Bi=apply(as.matrix(Aij[,-1]),1,sum) flag=1 if (sum(Bi[c(1:2)]!=0)) { flag=0 } ma.scf=c(rbind(Ai[-1],Bi[-(nf+1)]))[-c(1:4)] #A_3,B_2,A_4,B_3... ma.w2=c(ma.scf[1:3]) ma.wcc=c(3*ma.scf[1]+ma.scf[2],ma.scf[3]) ma.w1=c(Ai[3+1:2],Bi[3+1],Ai[6]) TW=as.matrix(def.sub[,c(1:nf)]) if ( any(rowSums(TW[-c(1:np),])==1) ) stop("aliased main effect") B=as.matrix(def.sub[,-c(1:nf)]) nbase=c(1, cumprod(rep(2, (b-1)))) fbase=c(1, cumprod(rep(2, (nf-1)))) b.base=nbase%*%t(B) m=diag(nf) TWsub=matrix(TW[b.base==0,], ncol=nf) #defining words alias.set=rep(list(c()), N-1) mi.vec=rep(list(c()),nStrata) eff.used=fbase%*%t(TWsub) for (i in 1:(nb-1)) { TWB=matrix(TW[b.base==i,], ncol=nf) alias.set[[i]]=TWB len.word=apply(TWB,1, sum) if (sum( len.word[-1]==1) >=1) stop("main effect aliased") if (any(len.word==1)) { } else { mi.vec[[1]]=c(mi.vec[[1]], length(len.word[len.word==2])) } eff.used=c(eff.used, fbase%*%t(TWB) ) } k=(nb-1) for (i in 1:nf) { eff=m[i,] eff.ind=sum(eff*fbase) if (!(eff.ind %in% eff.used)){ k=k+1 alias.set[[k]]=rbind(m[i,], t(apply(TWsub, 1, fin.alias, m[i,]) )) eff.used=c(eff.used, fbase%*%t(alias.set[[k]])) } } ind.sta1=k for (i in 2:(nf)){ out=combn(nf,i) a=fbase%*%apply(out,2, fin.sum, m) #if ind in eff.used the interaction is aliased with some effects if (k<(N-1) & !(all(a%in%eff.used))) { ind=t(out[,!(a%in%eff.used)]) for (j in 1:dim(ind)[1]) { eff=apply(m[ind[j,],],2,sum)%%2 if (sum(eff)==1) stop("main effect alias") eff.ind=sum(eff*fbase) if (!(eff.ind %in% eff.used) ) { k=k+1 alias.set[[k]]=rbind(eff, t(apply(TWsub, 1, fin.alias, eff ))) len.word=apply(alias.set[[k]],1, sum) mi.vec[[2]][(k-ind.sta1)]=length(len.word[len.word==2]) eff.used=c(eff.used, fbase%*%t(alias.set[[k]])) if(k>(N-1)) break } } } if(k>(N-1)) break } mst=matrix(0, ncol=nStrata, nrow=2) for (i in 1: nStrata) { mst[,i]=c(sum(unlist(mi.vec[c(i:nStrata)])), sum(unlist(mi.vec[c(i:nStrata)])^2)) } list(index=index,gen=gen,# Aij=Aij,def.sub=def.sub,dm_suf=dm_suf, # alias.set=alias.set, mi.vec=mi.vec, ma=Ai[-c(1:3)], ma.scf=ma.scf , ma.w2=ma.w2, ma.w1=ma.w1, ma.wcc=ma.wcc, mst=mst,flag) } SSP=function(nfat,nUnit, def.ind, sp.ind){ ## nfat=c(n1, n2) nUnit=c(w,s) ## def.ind: subplot columns ## sp.ind: splitting columns N=prod(nUnit) nStrata=length(nUnit) n=sum(nfat) n1=nfat[1] n2=nfat[2] nMatrix=log2(N) ret=ret.matrix(nMatrix) base=c(1, cumprod(rep(2, (nMatrix-1)))) p.col=def.ind b.col=c(base[1:n1], sp.ind) q=nMatrix-n1 p=length(c(def.ind, sp.ind)) cat("(n1, n2; w s) = (", c(n1, n2), ";", nUnit, ")\n") cat("subplot columns:", def.ind, "\n") cat("splitting column:", sp.ind) cat(" \n") a=bkwd2mi(nfat, nUnit, p.col, b.col) mst=a$mst mi.vec=a$mi.vec wlp=a$ma cat("mi at whole-plot stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi at subplot stratum: ", sort(mi.vec[[2]], decreasing=T),"\n") mst=c(mst) cat("(b) and (c) : (", mst[c(1,2)],"); (", mst[c(3,4)], ")\n") cat("word length pattern: ", wlp, "\n") cat("\n") } SSP.wd=function(nfat,nUnit, wd, bkwd){ ## nfat=c(n1, n2) nUnit=c(w,s) ## def.ind: subplot columns ## sp.ind: splitting columns nStrata=length(nUnit) N=prod(nUnit) n1=nfat[1];n2=nfat[2] nMatrix=log2(N) ret=ret.matrix(nMatrix) base=c(1,2,4,8, 16, 32, 64)# cumprod(rep(2, (nMatrix-1)))) p1.col=c() b.col=c() if (!is.null(wd)) { for (i in 1:length(wd)){ tmp=wd2index.ssp(wd[i], nfat, nUnit) p1.col=c(p1.col,tmp) } } if (!is.null(bkwd)) { for (i in 1:length(bkwd)){ tmp=wd2index.ssp(bkwd[i], nfat, nUnit) b.col=c(b.col,tmp) } } a=bkwd2mi(nfat, nUnit, p1.col, c(base[1:n1],b.col) ) mst=a$mst mi.vec=a$mi.vec wlp=a$ma cat("(n1, n2; w s) = (", nfat, ";", nUnit, ")\n") cat("subplot words: ", wd, "(", p1.col, ") \n") cat("splitting words: ", bkwd, "(", b.col, ")\n") cat("mi in whole-plot stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi in subplot stratum: ", sort(mi.vec[[2]], decreasing=T),"\n") mst=c(mst) cat("(b) and (c) : (", mst[c(1,2)],"); (", mst[c(3,4)], ")\n") cat("word length pattern: ", wlp, "\n") cat("\n") } find.SSP=function(nfat, nUnit){ ## nfat=c(n1, n2) nUnit=c(w,s) t0 <- proc.time() N=prod(nUnit) m=log2(nUnit) nStrata=length(nUnit) nMatrix=sum(m) M.ind=c() K=diag(m[1]+m[2]) nf=sum(nfat) key=matrix(NA, nrow=nf, ncol=nMatrix) if (!(nfat[1]1 & rowSums(as.matrix(BD[, c(1:m[2])]))!=0 ), ]) #can be set for SP sz=nrow(BD) Indice=rep(1, sz) for (j in 1:sz) { for (k in (m[2]+nfat[1]+1):(m[1]+m[2])){ if (all(BD[j,]==K[k,])) Indice[j]=0 } } BD=BD[Indice==1,] Vec=t(combn(nrow(BD), no.add)) if(nrow(Vec)>20000) sam.ind=sample(nrow(Vec), 20000) else sam.ind=c(1:nrow(Vec)) for (j in sam.ind ) { b.col=matrix(BD[Vec[j,], ], nrow=no.add) if (sum(rowSums(t(b.col))%%2)!=0) { # key=rbind(K, b.col ) key[1:nMatrix,]=K key[(nMatrix+1):nf,]=b.col mst=keytomi(key, nUnit)$mst if (is.null(M.mst)) { M.mst <<- matrix(c(t(mst)), nrow=1) M.vec <<- matrix(c(key), nrow=1) } else { flag=fun3(c(t(mst)), M.mst) if (flag==1) { M.mst <<- rbind(M.mst, c(t(mst))) M.vec <<- rbind(M.vec, c(key)) } } } } } ind=1 if (nrow(M.mst)>1) { ind=sort.order.sum(M.mst) M.mst <<- M.mst[ind,] M.vec <<- M.vec[ind,] M.mi <<- matrix(M.mst[1,], nrow=1) M.op <<- matrix(M.vec[1,], nrow=1) for (i in 2:nrow(M.mst)) { key=M.vec[i,] mmi=M.mst[i,] comp.mst(mmi, key) } } else { M.mi <<- M.mst M.op <<- M.vec } cat("nfat:", nfat, " nUnit:", nUnit, "\n") cat("number of admissible design(s): ", nrow(M.mi), "\n") for (k in 1:nrow(M.mi)) { key<-matrix(M.op[k,], ncol=nMatrix) colnames(key)=c(paste("S", 1:m[2], sep=""), paste("W", 1:m[1], sep="")) rownames(key)=c(LETTERS[15+1:min(m[2], nfat[2])], LETTERS[1:nfat[1]], c(LETTERS,letters)[15+(min(m[2],nfat[2])+1):nfat[2]] ) a=keytomi(key, nUnit) mst=a$mst mi.vec=a$mi.vec tmp=keytodefwd(key,nfat, nUnit) wd=list(c(), c()) p.col=list(c(),c()) for (i in 1:nStrata) { vec=tmp[[i]] for (j in 1:length(vec)){ b.wd<- paste(sort(vec[[j]]),sep="", collapse="") wd[[i]][j] <- b.wd p.col[[i]][j]=wd2index.ssp(b.wd, nfat,nUnit) } } cat("defining words: ", wd[[2]], "(", p.col[[2]], ") \n") cat("splitting words: ", wd[[1]], "(", p.col[[1]], ")\n") cat("mi in whole-plot stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi in subplot stratum: ", sort(mi.vec[[2]], decreasing=T), " \n") cat("(b) and (c) : (", mst[,1],"); (", mst[,2], ")\n") wlp=keytowlp(key) cat("word length pattern: ",wlp, "\n") cat("\n") } proc.time() - t0 } fun.add=function(nfat, nUnit, K){ m=log2(nUnit) nMatrix=nrow(K) no.add=nfat[1]+nfat[2]-nMatrix BD=expand.grid(rep(list(c(0,1)),(nMatrix))) BD=as.matrix(BD[ (rowSums(BD)>1 & rowSums(as.matrix(BD[, c(1:m[3])]))!=0) , ]) #can be set for SP sz=nrow(BD) Indice=rep(1, sz) for (j in 1:sz) { for (k in (m[2]+m[3]+1):(nMatrix)){ if (all(BD[j,]==K[k,])) Indice[j]=0 } } BD=BD[Indice==1,] Vec=t(combn(nrow(BD), no.add)) if(nrow(Vec)>20000) sam.ind=sample(nrow(Vec), 20000) else sam.ind=c(1:nrow(Vec)) for (j in sam.ind ) { b.col=matrix(BD[Vec[j,], ], nrow=no.add) if (sum(rowSums(t(b.col))%%2)!=0) { #check if the added columns are ind key=rbind(K, b.col ) mst=keytomi(key, nUnit)$mst if (is.null(M.mst)) { M.mst <<- matrix(c(t(mst)), nrow=1) M.vec <<- matrix(c(key), nrow=1) } else { flag=fun3(c(t(mst)), M.mst) if (flag==1) { M.mst <<- rbind(M.mst, c(t(mst))) M.vec <<- rbind(M.vec, c(key)) } } } } } bkbsp2mi=function(nfat, nUnit, p1, p2, def.ind, sp.ind){ N=prod(nUnit) nblk=log2(nUnit) n1=nfat[1];n2=nfat[2] bz=nUnit[1];nwp=nUnit[2];nsp=nUnit[3] nStrata=length(nUnit) nf=n1+n2 np=length(def.ind) npcol=def.ind b.col=sp.ind bw=length(sp.ind) nr=np+bw L=2^(nr)-1 nMatrix=log2(N) ret=ret.matrix(nMatrix) nb=2^bw nrcol=c(npcol, b.col) index=c(n1,n2,p1,p2,bw) if (np==0) { gen=cbind(matrix(ret[, nrcol], byrow=T, nrow=nr), diag(nr)) } else if (p1==0 & p2>0) { m1=rbind(diag(p2),matrix(0, nrow=bw,ncol=p2)) m2=rbind(matrix(0, nrow=p2, ncol=bw), diag(bw)) gen=cbind(matrix(ret[1:(nf-np), nrcol],byrow=T, nrow=nr), m1, m2) } else if (p1>0 & p2==0) { m1=rbind(diag(p1),matrix(0, nrow=bw,ncol=p1)) m2=rbind(matrix(0, nrow=p1, ncol=bw), diag(bw)) gen=cbind(matrix(ret[1:(n1-p1), nrcol],byrow=T, nrow=nr), m1, matrix(ret[-c(1:(n1-p1)), nrcol],byrow=T, nrow=nr), m2) } else { m1=rbind(diag(p1),matrix(0, nrow=(p2+bw),ncol=p1)) m2_1=cbind(diag(p2),matrix(0, nrow=p2, ncol=bw)) m2_2=cbind(matrix(0, nrow=bw, ncol=(p2)), diag(bw)) m2=rbind(matrix(0, nrow=p1, ncol=(p2+bw)),m2_1,m2_2) gen=cbind(matrix(ret[1:(n1-p1), nrcol],byrow=T, nrow=nr), m1, matrix(ret[-c(1:(n1-p1)), nrcol],byrow=T, nrow=nr), m2) } def.sub=c() for (i in 1:nr){ out=combn(nr,i) def.sub=rbind(def.sub, t(apply(out, 2, ele.sum, gen))) } dm_suf=matrix(0, nrow=L,ncol=3) for (i in 1:L) dm_suf[i,]=c(sum(def.sub[i,1:n1]), sum(def.sub[i, n1+(1:n2)]), sum(def.sub[i,nf+(1:bw)])) A_suf=cbind(apply(dm_suf[,c(1,2)],1,sum),dm_suf[,3]) Aij=matrix(0, nrow=(nf+1), ncol=(bw+1)) for (i in 1:L) { Aij[(A_suf[i,1]+1),(A_suf[i,2]+1)]=Aij[(A_suf[i,1]+1),(A_suf[i,2]+1)]+1 } Ai=Aij[,1] Bi=apply(as.matrix(Aij[,-1]),1,sum) if (sum(Bi[c(1:2)]>n1)) { flag=0 return(list(flag=0, gen=gen, mst=matrix(0,2,2))) } ma.ma=c(c(rbind(Ai[-1],Bi[-(nf+1)]))[-c(1:4)],Bi[(nf+1)]) #A_3,B_2,A_4,B_3... wlp=Ai ma.scf=c(rbind(Ai[-1],Bi[-(nf+1)]))[-c(1:4)] #A_3,B_2,A_4,B_3.. TW=as.matrix(def.sub[,c(1:nf)]) B=as.matrix(def.sub[,-c(1:nf)]) nbase=c(1, cumprod(rep(2, (bw-1)))) fbase=c(1, cumprod(rep(2, (nf-1)))) b.base=nbase%*%t(B) Twsp=apply(dm_suf[b.base!=0,c(1:2)],1,sum) if (any(Twsp==1)) stop("some main effect is confounded with blocks") TWsub=matrix(TW[b.base==0,], ncol=nf) #defining words BWsub=matrix(TW[b.base!=0,], ncol=nf) eff.used=fbase%*%t(TWsub) alias.set=rep(list(c()), N-1) for (i in 1:(nb-1)) { TWB=matrix(TW[b.base==i,], ncol=nf) alias.set[[i]]=TWB if (sum(rowSums(alias.set[[i]])==1)!=0) stop("a main effect in block stratum\n") eff.used=c(eff.used, fbase%*%t(TWB) ) } m=diag(nf) m.wp=matrix(m[1:n1,],nrow=n1) for (i in 1:n1) { alias.set[[i+nb-1]]=rbind(m.wp[i,], t(apply(TWsub, 1, fin.alias, m.wp[i,]) )) eff.used=c(eff.used, fbase%*%t(alias.set[[i+nb-1]])) } k=n1+nb-1 df.wp.ind=nb*(nwp-1) wp.ind=0 #index for no of whole-plot interactions in wholeplot stratum which is not in the block if (n1>=2){ for (i in 2:(n1)){ out=combn(n1,i) a=fbase%*%apply(out,2, fin.sum, t(m.wp)) #if ind in eff.used the interaction is aliased with some effects if (wp.ind=df.wp.ind) break } } } if(wp.ind>=df.wp.ind) break } } sz.mwp=wp.ind+n1 #numbers of initial whole plot effects if(wp.ind<=df.wp.ind){ for (i in 1:sz.mwp) { w.eff=alias.set[[i+nb-1]][1,] for (j in 1:(nb-1)) { b.eff=alias.set[[j]][1,] eff=(w.eff+b.eff)%%2 eff.ind=sum(eff*fbase) if (!(eff.ind %in% eff.used) ) { k=k+1 wp.ind=wp.ind+1 alias.set[[k]]=rbind(eff, t(apply(TWsub, 1, fin.alias, eff ))) eff.used=c(eff.used, fbase%*%t(alias.set[[k]])) if (sum(rowSums(alias.set[[k]])==1)!=0) { cat("a subplot factor in whole-plot stratum\n") } if(wp.ind>=df.wp.ind) break #} else { #stop("not a good setting") } } if(wp.ind>=df.wp.ind) break } } if ((wp.ind+n1) != nUnit[1]*(nUnit[2]-1)) cat("odd setting\n") m.sp=matrix(m[-c(1:n1),], nrow=n2) for (i in 1:n2) { eff=m.sp[i,] eff.ind=sum(eff*fbase) if ( !(eff.ind%in%eff.used)) { k=k+1 alias.set[[k]]=rbind(eff, t(apply(TWsub, 1, fin.alias, eff ))) eff.used=c(eff.used, fbase%*%t(alias.set[[k]])) } } sp.ind=0 if (n2>=2) { for (i in 2:nf) { out=combn(nf,i) a=fbase%*%apply(out,2, fin.sum, m) #if ind in eff.used the interaction is aliased with some effects if (k<(N-1) & !(all(a%in%eff.used))) { ind=t(out[,!(a%in%eff.used)]) for (j in 1:dim(ind)[1]) { eff=apply(m[ind[j,],],2,sum)%%2 eff.ind=sum(eff*fbase) if (!(eff.ind %in% eff.used) ) { k=k+1 sp.ind=sp.ind+1 alias.set[[k]]=rbind(eff, t(apply(TWsub, 1, fin.alias, eff ))) eff.used=c(eff.used, fbase%*%t(alias.set[[k]])) if (sum(rowSums(alias.set[[k]])==1)!=0) { cat("another main effect\n") } if(k>=(N-1)) break } } } if(k>=(N-1)) break } } df=c((nUnit[1]-1), nUnit[1]*(nUnit[2]-1), nUnit[1]*nUnit[2]*(nUnit[3]-1)) mi.vec=list(rep(NA,df[1]), rep(NA, df[2]-n1), rep(NA,df[3]-n2)) for (i in 1:df[1]) { tmp=rowSums(alias.set[[i]]) mi.vec[[1]][i]=sum(tmp==2) } if (length(mi.vec[[2]])!=0) { for (i in 1:(df[2]-n1)) { ii=i+df[1]+n1 tmp=rowSums(alias.set[[ii]]) mi.vec[[2]][i]=sum(tmp==2) } } for (i in 1:(df[3]-n2)) { ii=i+df[1]+df[2]+n2 tmp=rowSums(alias.set[[ii]]) mi.vec[[3]][i]=sum(tmp==2) } mst=matrix(0, ncol=nStrata, nrow=2) for (i in 1: nStrata) { mst[,i]=c(sum( unlist(mi.vec[c(i:nStrata)]), na.rm=T ), sum(unlist(mi.vec[c(i:nStrata)])^2, na.rm=T ) ) } mi=unlist(mi.vec) mi.a=length(which(mi==1)) list(nUnit=nUnit, mi.vec=mi.vec, mst=mst,wlp=wlp, wlp.scf=ma.scf) #Defwd=TWsub, Bkwd=BWsub, #alias.set=alias.set,# dm_suf=dm_suf, Ai=Ai,Bi=Bi, } BSP.wd=function(nfat, nUnit, wd1, wd2, wdbk) { n1=nfat[1];n2=nfat[2] N=prod(nUnit) n=n1+n2 nMatrix=log2(N) ret=ret.matrix(nMatrix) p1=length(wd1) p2=length(wd2) p1.col=c() p2.col=c() sp.ind=c() if (p1!=0) { for (i in 1:p1){ tmp=wd2index(wd1[i], n1, n2, p1, p2) p1.col=c(p1.col,tmp) } } if (p2!=0) { for (i in 1:p2){ tmp=wd2index(wd2[i], n1, n2, p1, p2) p2.col=c(p2.col,tmp) } } def.ind=c(p1.col, p2.col) for (i in 1:length(wdbk)){ tmp=wd2index(wdbk[i], n1, n2, p1, p2) sp.ind=c(sp.ind,tmp) } a=bkbsp2mi(nfat, nUnit, p1,p2, def.ind, sp.ind) mst=a$mst mi.vec=a$mi.vec wlp=a$wlp wlp.scf=a$wlp.scf cat("(n1, n2; w s) = (", c(n1, n2), ";", nUnit, ")\n" ) cat("defining words:", c(wd1, wd2), " (", def.ind, ") \n") cat("splitting/blocking words:",wdbk, " (", sp.ind, ") \n") cat("mi at block stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi at whole-plot stratum: ", sort(mi.vec[[2]], decreasing=T), "\n") cat("mi at subplot stratum: ", sort(mi.vec[[3]], decreasing=T),"\n") cat("(a), (b) and (c) : (", mst[,1],"); (", mst[,2],"); (", mst[,3], ")\n") cat("word length pattern: ", wlp[-c(1:3)], "\n") cat("MB word length pattern: ", wlp.scf, "\n") cat("\n") } BSP=function(nfat, nUnit, p.col, b.col) { n1=nfat[1];n2=nfat[2] N=prod(nUnit) n=n1+n2 m=log2(nUnit) nMatrix=log2(N) nn1=min(n1,nMatrix-m[3]) ret=ret.matrix(nMatrix) base=c(1,2,4,8,16,32) p1=sum(p.col<=base[nn1+1]) p2=length(p.col)-p1 a=bkbsp2mi(nfat, nUnit, p1,p2, p.col, b.col) mst=a$mst mi.vec=a$mi.vec wlp=a$wlp wlp.scf=a$wlp.scf cat("(n1, n2; b w s) = (", c(n1, n2), ";", nUnit, ")\n" ) cat("defining words:", p.col ) cat(" splitting/blocking words:", b.col,"\n") cat("mi at block stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi at whole-plot stratum: ", sort(mi.vec[[2]], decreasing=T), "\n") cat("mi at subplot stratum: ", sort(mi.vec[[3]], decreasing=T),"\n") cat("(a), (b) and (c) : (", mst[,1],"); (", mst[,2],"); (", mst[,3], ")\n") cat("word length pattern: ", wlp[-c(1:3)], "\n") cat("MB word length pattern: ", wlp.scf, "\n") cat("\n") } find.BSP=function(nfat, nUnit){ #t0 <- proc.time() n1=nfat[1];n2=nfat[2] N=prod(nUnit) m=log2(nUnit) nMatrix=sum(m) M.ind=c() K=diag(nMatrix) nStrata=length(nUnit) colnames(K)=c(paste("S", 1:m[3], sep=""), paste("W", 1:m[2], sep=""), paste("B", 1:m[1], sep="")) M.mst <<- c() #column for the key S1 S2... W1 W2 M.vec <<- c() B=as.matrix(expand.grid(rep(list(c(0,1)),m[3]))) BC=as.matrix(expand.grid(rep(list(c(0,1)),m[2]))) if (n1m[2]){ if ((n1+m[3])1 & rowSums(as.matrix(BD[, c(1:m[3])]))!=0) , ]) #can be set for SP if (nrow(BC)<=m[1]) { out2=as.matrix(expand.grid(rep(list(c(1:nrow(BC))),m[1]))) } else { out2=t(combn(nrow(BC), (m[1]))) } for (j in 1:nrow(out2)){ K[(m[2]+m[3]+1):nMatrix, (m[3]+1):(m[2]+m[3])]=BC[out2[j,],] wpbase=nMatrix-m[3] wpadd=(n1-wpbase) M.ind=c() for (ii in 2:wpbase) { out3=t(combn(wpbase, ii)) for (kk in 1:nrow(out3)) { eff=apply(K[m[3]+out3[kk,],] ,2,sum)%%2 if (sum(eff[1:(m[2]+m[3])])!=0) M.ind=rbind(M.ind, eff) } } key=matrix(NA, nrow=sum(nfat), ncol=nMatrix) vec=t(combn(nrow(M.ind), wpadd)) for (ii in 1:nrow(vec)) { key[1:nMatrix,]=K key[nMatrix+1:wpadd,]= M.ind[vec[ii,], ] if (m[3]20000) sam.ind=sample(nrow(Vec), 20000) else sam.ind=c(1:nrow(Vec)) for (j in sam.ind ) { b.col=matrix(BD[Vec[j,], ], nrow=(n2-m[3])) if (sum(rowSums(t(b.col))%%2)!=0) { #check if the added columns are ind key[(m[3]+n1+1):sum(nfat),]=b.col mst=keytomi(key, nUnit)$mst if (is.null(M.mst)) { M.mst <<- matrix(c(t(mst)), nrow=1) M.vec <<- matrix(c(key), nrow=1) } else { flag=fun3(c(t(mst)), M.mst) if (flag==1) { M.mst <<- rbind(M.mst, c(t(mst))) M.vec <<- rbind(M.vec, c(key)) } } } } } else { mst=keytomi(key, nUnit)$mst if (is.null(M.mst)) { M.mst <<- matrix(c(t(mst)), nrow=1) M.vec <<- matrix(c(key), nrow=1) } else { flag=fun3(c(t(mst)), M.mst) if (flag==1) { M.mst <<- rbind(M.mst, c(t(mst))) M.vec <<- rbind(M.vec, c(key)) } } } } } } } ind=1 if (nrow(M.mst)>1) { ind=sort.order.sum(M.mst) M.mst <<- M.mst[ind,] M.vec <<- M.vec[ind,] M.mi <<- matrix(M.mst[1,], nrow=1) M.op <<- matrix(M.vec[1,], nrow=1) for (i in 2:nrow(M.mst)) { key=M.vec[i,] mmi=M.mst[i,] comp.mst(mmi, key, nStrata=3) } } else { M.mi <<- M.mst M.op <<- M.vec } # cat("nfat:", nfat, " nUnit:", nUnit, "\n") # cat("number of admissible design(s): ", nrow(M.mi), "\n") for (k in 1:nrow(M.mi)) { key<-matrix(M.op[k,], ncol=nMatrix) colnames(key)=c(paste("S", 1:m[3], sep=""), paste("W", 1:m[2], sep=""), paste("B", 1:m[1], sep="")) if (m[3]==n2) { rownames(key)=c(LETTERS[15+1:n2], LETTERS[1:n1]) }else { rownames(key)=c(LETTERS[15+1:min(m[3], nfat[2])], LETTERS[1:nfat[1]], LETTERS[15+(min(m[3],nfat[2])+1):nfat[2]]) } a = keytomi(key, nUnit) mst=a$mst mi.vec=a$mi.vec tmp=keytodefwd(key,nfat, nUnit) wd =c() for (ii in 1:length(tmp[[2]])) wd=c(wd, paste(sort(tmp[[2]][[ii]]), sep="", collapse="")) p1=0 for (ii in 1:length(wd)) { vec=wd[ii] ind=pmatch(unlist(strsplit(vec,"")),LETTERS[1:n1]) if (sum(is.na(ind))==0) p1=p1+1 } p2=length(wd)-p1 p.col=c() for (ii in 1:length(wd)) { p.col=c(p.col, wd2index(wd[ii],n1,n2,p1,p2)) } bkwd =c() b.col=c() for (ii in 1:length(tmp[[1]])) { vec=paste(sort(tmp[[1]][[ii]]), sep="", collapse="") bkwd=c(bkwd, vec) b.col=c(b.col, wd2index(vec,n1,n2,p1,p2)) } cat("\n") cat("defining words: ", wd, "(", p.col, ") ") cat(" splitting/blocking words: ", bkwd, "(", b.col, ") \n") cat("mi in block stratum: ", sort(mi.vec[[1]], decreasing=T), "\n") cat("mi in whole-plot stratum: ", sort(mi.vec[[2]], decreasing=T), "\n") cat("mi in subplot stratum: ", sort(mi.vec[[3]], decreasing=T), " \n") cat("(a), (b) and (c) : (", mst[,1],"); (", mst[,2],"); (", mst[,3], ")\n") wlp=keytowlp(key) cat("word length pattern: ",wlp, "\n") cat("\n") } cat("\n") # proc.time() - t0 }