##################################################################
#
## updated revision on 2016/05/20
## functions for computer information in Tsai and Gilmour (2010)
##
##  cal.Q3L(d) -- Q-value for a three-level design
##  cal.gwc(d) -- generalized word count for a  three-level design
##  cal.Qb3L(d, p) -- Qb-value for three-level designs coded with a given prior p
##  cal.Db3L(d, p, tau) -- Q-value for three-level designs coded  with a given
##  prior p and tau
##
## details can be found in gwc_supp.pdf
## data files used in the supplimentary
##
## Col_F6R19.out for all six three-level designs in 18 runs
## pb16.out: 16-run two-level Placket-Burman design
## Qb18T6.out: the four six-factor designs studied
##
##################################################################


dim2=function(mat) {
  dim(mat)[2]
}


myFun = function(arg){
  prod(arg)
}



proj.calc = function(ind.col, Design) {
  sub.matrix=as.matrix(Design[, ind.col] )
  prodV = apply(sub.matrix, 1, myFun)
  tmp=sum(prodV)^2
  return(tmp)
}


matrix.calc.fin = function(data, noRun) {
  prodV = apply(data, 1, myFun)
  (sum(prodV))^2/(noRun*noRun)
}

num.mod2 = function(a.min, noCol, noRun) {
#number of model for first order maximal model 
   tmp=0
   a.bound.min=min((noRun-1), noCol)

   if (a.min<=a.bound.min)  {
       for (a in a.min:a.bound.min) {
        tmp = tmp + choose(noCol-a.min, a-a.min)
        }   
   }
}


F.mod2.int = function(c.val, a, a.min, c.min, c.max, noCol, noRun) { 
   if (c.val+a < noRun ) {
     ret = choose(noCol-a.min, a-a.min) * choose(c.max-c.min, c.val-c.min)
      } else {
       ret=0    
     }  
}


 
delta2.int = function(ind, noCol, noRun) {
#for two-level second-order maximal model
#apply function F.mod to find delta for each combination of a,c 
  a.bound.min=min(noRun-1, noCol)
  tmp=0
  if (ind[1] <=noCol) {
      for(a in ind[1]:a.bound.min) {
        c.max=choose(a,2)
        c.ind=ind[2]:c.max
        tmp=tmp+sum(unlist(lapply(c.ind,  F.mod2.int, a, ind[1], ind[2], c.max, noCol, noRun)))
     }
   }
  (tmp)
}

matrix.calc.real = function(nMatrix, combList, vec, myArray, noRun) {
  matrix.selected = c()
  for(i in 1:nMatrix) {
    combSlice = combList[[i]]
    matrix.selected = cbind(matrix.selected, 
                            myArray[, combSlice[vec[i],], i])
  }
  myRes = matrix.calc.fin(matrix.selected, noRun)
  myRes
}


matrix.calc.rec = function(nMatrix, depth, combList, vec, colUnion, myArray, noRun) {
  ret = c()
  if(depth==nMatrix+1) {
    ret = matrix.calc.real(nMatrix, combList, vec, myArray, noRun)
  } else {
    combSlice = combList[[depth]]
    if(dim(combSlice)[1]==0) {
      ret = c(ret, matrix.calc.rec(nMatrix, depth+1, combList, c(vec,0), colUnion, myArray, noRun))
    } else {
      for(i in 1:dim(combSlice)[1]) {
         if( length(colUnion)==0 || 
            length(intersect(colUnion, combSlice[i,]))==0)
            ret = c(ret, matrix.calc.rec(nMatrix, depth+1, combList, c(vec, i), 
                   union(colUnion, combSlice[i,]), myArray, noRun))
      }
    }
  }
  ret
}




matrix.choose.vector = function(nMatrix, maxCount, cv, myArray, noRun) {
  for(i in 1:nMatrix) {
    if(i==1) {
      if(cv[i]==0) combList = list(matrix(0,0,0))
      else combList = list(t(combn(maxCount, cv[i])))
    } else {
      if(cv[i]==0) combList = append(combList, list(matrix(0,0,0)))
      else combList = append(combList, list(t(combn(maxCount, cv[i]))))
    }
  }
  ret = matrix.calc.rec(nMatrix, 1, combList, c(), c(), myArray, noRun)
  ind = unlist(lapply(combList, dim2))
  c(sum(ind), ind, sum(ret))
}



matrix.choose.count = function(nMatrix, maxCount, totalCount, myArray, noRun) {
  ceiling = (maxCount + 1) ^ nMatrix
  cv = rep(0, nMatrix)
  ret = c()
  for(counter in 1:(ceiling-1)) {
    val = counter
    for(i in nMatrix:1) {
      cv[i] = val %% (maxCount + 1)
      val = val %/% (maxCount + 1)
    }

    if(sum(cv)==totalCount) {
      ret = rbind(ret, matrix.choose.vector(nMatrix, maxCount, cv, myArray, noRun))
    }
  }
  ret
}


mod.num.calc=function(c.val, a.val, c.max, a.min, c.min, noF){
    mod.num.calc=choose(noF-a.min, a.val-a.min)*choose(c.max-c.min, c.val-c.min)
}

mod.num=function(ac, a.min, c.min, noF){
    #two-level design
    a.val=ac[1]
    c.max=ac[2]
    unlist(lapply(c(0:c.max), mod.num.calc, a.val, c.max, a.min, c.min, noF))
}


F.mod = function(bc, a, a.min, b.min, c.min,  c.max, noCol, noRun) {
# number of models with a.min, b.min, c.min effects included and 
# whose number of parameters are less than run size
   if (a + sum(bc) < noRun) {
     ret = c(a, bc,  choose(noCol-a.min, a-a.min) * choose(a-b.min, bc[1]-b.min) * choose(c.max-c.min, bc[2]-c.min))
  } else {
	ret = numeric(4)
  }
  ret
}


mod.delta=function(abc, noCol, noRun){
  #for second-order maximal model
#apply function F.mod to find delta for each combination of a,b,c 
  a.min=abc[1]
  b.min=abc[2]
  c.min=abc[3]
  a.bound.min=min(noRun-1,noCol)
  tmp=0
  if (a.min<=noCol) {
  
   for(a in a.min:a.bound.min) {
      b.max=a
     c.max=choose(a,2)
      bc.max = noRun - a
      a1 = matrix( rep(b.min:b.max, c.max - c.min + 1), 
		c.max - c.min + 1, b.max - b.min + 1, byrow=T)
      a2 = matrix( rep(c.min:c.max, b.max - b.min + 1), 
		c.max - c.min + 1, b.max - b.min + 1)
      A = array(c(a1,a2), c(c.max - c.min + 1, b.max - b.min + 1, 2) )
      R= apply(A, c(1,2), F.mod, a, a.min, b.min, c.min, c.max, noCol, noRun)
  
      tmp= tmp + sum(R[4,,])
     }
  } else {tmp=0}
  tmp
}




delta.3.qual=function(a.min ,noCol, noRun) {
   a=1:noCol
   a.bound=length(a[(2*a+choose(a,2)*4)<noRun])
   tmp=0
   if (a.min<=a.bound) {
    for (i in a.min:a.bound) {
       tmp=tmp+choose(noCol-a.min,i-a.min)
    }
   }
  return(tmp)
}

total.effect = function(arg, len.arg){
   sum(arg*seq(len.arg)) }


word.sum = function(arg, ind, worList) {
  tmp=dim2(worList)
  sum(worList[ind==arg, tmp]) }


b.word =function(arg,worList){
    true.index=c()
    tmp=c() 
    ind=dim2(worList)
    for (i in 1:dim(worList)[1]) {
      if (all(worList[i,seq(length(arg))]==arg)=="TRUE")
           true.index=i
      }
    tmp=worList[true.index, ind]
    if (is.null(true.index)) tmp=0
     return(tmp)  
} 


 

pr.mod=function(bc,a, p1, p2, p3, noF) {
    b=bc[1]
    c=bc[2] 
    c.max=a*(a-1)/2
    tmp=p1^a*(1-p1)^(noF-a)*p2^b*(1-p2)^(a-b)*p3^c*(1-p3)^(c.max-c)
    nu.mod=choose(noF, a) * choose(a, b)* choose(c.max, c)
    ret=c(a, bc, nu.mod, tmp, tmp*nu.mod)
#    print(ret)
}


pr.mat = function(pp, noF, noRun=18) {
  p1=pp[1]
  p2=pp[2]
  p3=pp[3]
  tmp=c()
  a.bound.min=min(noRun-1,noF)
  for(a in 0:a.bound.min) {
      b.max=a
      c.max=choose(a,2)
      bc.max = noRun - a
      b.max=a
      c.max=choose(a,2)
      bc.max = noRun - a

      a1 = matrix( rep(0:b.max, c.max + 1), 
		c.max + 1, b.max + 1, byrow=T)
      a2 = matrix( rep(0:c.max, b.max + 1), 
		c.max  + 1, b.max  + 1)
      A = array(c(a1,a2), c(c.max + 1, b.max + 1, 2) )
      R=matrix(unlist(apply(A, c(1,2), pr.mod, a, p1, p2, p3, noF)),byrow=T, ncol=6)
     tmp=rbind(tmp, R)
  }
  tmp
}


pr.adj=function(prmat, noRun){
   ind=apply(prmat[,c(1,2,3)],1, sum) >= noRun
       pr.gamma=sum(prmat[ind,6])
       mod.gamma=sum(prmat[ind,4])
   if (any(ind)==TRUE)   adj=pr.gamma/mod.gamma   else adj=0

   pr=prmat[,5]
   ind2=apply(prmat[,c(1,2,3)],1, sum)==(noRun-1)
   pr[ind2]=pr[ind2]+adj
   pr[ind]=0
   pradj=cbind(prmat[,c(1:4)], pr) 
   list(pradj=pradj, adj=adj)
}


plot.hist=function(pp, noF=6, noRun=18) {

  prmat=pr.mat(pp, noF, noRun)
  pr.hist=prmat[,4]*prmat[,5]
  nopara=2*noF+noF*(noF-1)/2+1
  pr=numeric(nopara-1)
  for (i in 1:(nopara-1)) {
     pr[i]=sum(pr.hist[apply(prmat[,c(1,2,3)],1, sum)==i])
   }
  pr=c(0,pr) 
  lab1=sum(pr[1:(noRun)])
  barplot(pr,col="wheat", ylim=c(0,(max(pr)+.05)), xlim=c(0,nopara+1),
           xlab="No of parameters", ylab="Probability", 
           main="Probability for different models being the true model, P(Size<18)")
  abline(v=(noRun+1), col=2,lwd=2)
  text(noRun-1, max(pr)+0.03, paste(round(100*lab1,2), "%",sep=""))
}


F.pr.mod=function(bc, a, a.min, b.min, c.min, c.max, adj, pp, noF, noRun) {
   p1=pp[1]
   p2=pp[2]
   p3=pp[3]
  b=bc[1]
    c=bc[2] 
   if (a + sum(bc) < (noRun-1)) {
         tmp=p1^a*(1-p1)^(noF-a)*p2^b*(1-p2)^(a-b)*p3^c*(1-p3)^(c.max-c)
         nu.mod = choose(noF-a.min, a-a.min) * choose(a-b.min, bc[1]-b.min) * choose(c.max-c.min, bc[2]-c.min)
         ret=c(a, bc, nu.mod, tmp, tmp*nu.mod)
   } else  if  (a + sum(bc) == (noRun-1)) {
              tmp=p1^a*(1-p1)^(noF-a)*p2^b*(1-p2)^(a-b)*p3^c*(1-p3)^(c.max-c)+adj
              nu.mod = choose(noF-a.min, a-a.min) * choose(a-b.min, bc[1]-b.min) * choose(c.max-c.min, bc[2]-c.min)
              ret=c(a, bc, nu.mod, tmp, tmp*nu.mod)
        }
      else
          ret=numeric(6)
   	
  ret
}


adj.xi=function(abc, adj, pp, noF, noRun) {
  a.min=abc[1]
  b.min=abc[2]
  c.min=abc[3]
  a.bound.min=min(noRun-1,noF)
  tmp=0
  if (a.min<=noF) {

    for(a in a.min:a.bound.min) {
      b.max=a
      c.max=choose(a,2)
      bc.max = noRun - a
      b.max=a
      c.max=choose(a,2)
      bc.max = noRun - a

      a1 = matrix( rep(b.min:b.max, c.max - c.min + 1), 
		c.max - c.min + 1, b.max - b.min + 1, byrow=T)
      a2 = matrix( rep(c.min:c.max, b.max - b.min + 1), 
		c.max - c.min + 1, b.max - b.min + 1)
      A = array(c(a1,a2), c(c.max - c.min + 1, b.max - b.min + 1, 2) )
      R=apply(A, c(1,2), F.pr.mod, a, a.min, b.min, c.min, c.max, adj, pp, noF, noRun)
      tmp= tmp + sum(R[6,,])
     }
  } else {tmp=0}
  tmp
# print(tmp)
}


cal.Wmat=function(xi.s, noF, noRun) {
   s=rep(0,14)
   s[1]=xi.s[1]
   s[2]=xi.s[2]
   s[3]=xi.s[3]
   s[4]=xi.s[3]
   s[5]=xi.s[4]
   s[6]=xi.s[6]
   s[7]=xi.s[5]
   s[8]=xi.s[6]
   s[9]=xi.s[7]
   s[10]=xi.s[8]
   s[11]=xi.s[9]
   s[12]=xi.s[10]
   s[13]=xi.s[11]
   s[14]=xi.s[12]
   n=noF 
   sz2=(n+1)*(n+2)/2
   W.mat=matrix(0, nrow=sz2, ncol=sz2)

   out=combn(noF, 2)
   ind_a=out[1,]
   ind_b=out[2,]
#  i1=1
#  ind_a=c() 
#  ind_b=c()
#  while (i1<n) {
#    kk=c(1:i1)
#    b=rep(1, i1)*(i1+1)
#    ind_a=c(ind_a, kk) 
#    ind_b=c(ind_b, b)
#    i1=i1+1
#  }


    ni=length(ind_a)
    pl=c(1:n);
    pq=pl+n;
    pi=c((2*n+1):(2*n+ni));


   W.mat[(pl+1),(pi+1)]=matrix(1, nrow=n,ncol=ni)*s[11];   # 301
   W.mat[(pq+1),(pi+1)]=matrix(1, nrow=n,ncol=ni)*s[12];  # 311
   for (i in pl) {
      j=i+1;
      W.mat[j,j]=s[2];     # 100
      W.mat[(j+n), (j+n)]=s[3];   # 110
      W.mat[j, (pq+1)]=matrix(1, nrow=1, ncol=n)*s[7];  #210
      W.mat[j, (j+n)]=s[4];  #110
      while (j<n+1) {
        W.mat[(i+1), (j+1)]=s[5];  #200
        W.mat[(i+1+n), (j+1+n)]=s[10]; #220
       j=j+1;
      }
    ind=c(which(ind_a==i), which(ind_b==i));
    ind_i=pi[ind];
    no_ind=length(ind_i)
    W.mat[(i+1),(ind_i+1)]=rep(1, no_ind)*s[8]; #201
    W.mat[(i+1+n),(ind_i+1)]=rep(1, no_ind)*s[9]; #211
    }

  for (ii in c(1:ni)) {
    jj=ii+2*n+1;
    ij=jj;
    W.mat[jj,jj]=s[6];  #201
    while (ij<sz2){ 
      W.mat[jj,(ij+1)]=s[13]; #302
      ij=ij+1;
    }
 
       a1v=rep(1,ni);
       a2v=rep(1,ni);
       a1=ind_a[ii];
       b1=ind_b[ii];
       ind1=c(which(ind_a==a1), which(ind_a==b1));
       ind2=c(which(ind_b==a1), which(ind_b==b1));
       a1v[ind1]=0*ind1;
       a2v[ind2]=0*ind2;
       av=a1v*a2v;
       ind=which(av!=0);
       if (length(ind)!=0){
           W.mat[jj,(ind+2*n+1)]=((0*ind)+1)*s[14];  #402
       }
  }
  W.mat[1,1]=s[1];
  for (i in 2:sz2) {
    W.mat[1,i]=W.mat[i,i]
    for (j in 1:(i-1)) {
        W.mat[i,j]=W.mat[j,i]
    }
  }
  W.mat;
}



cal.gwc<-function(d,noLev=3) {

#d: a noRun*noF three-level design matrix with coding -1, 0, 1 
#calaculate the generalized word count
 noRun=nrow(d)
 noF=ncol(d)
 assign("Global.res", noRun, envir = .GlobalEnv)
 assign("Global.res", noF, envir = .GlobalEnv)
 assign("Global.res", noLev, envir = .GlobalEnv)
 
 colnames(d)=LETTERS[1:noF]    #give A,B, .. as the factor names
 rownames(d)=c(1:noRun)       #numbering the runs
 NM=noLev-1

 if (noLev==3) {
   myArray=array(NA, c(noRun,noF,NM))
   myArray[,,1] = d
   myArray[,,2] = sqrt(1/2)*(3 * myArray[,,1] ^2 - 2)
   myArray[,,1] = sqrt(3/2)*myArray[,,1]
  # print(myArray)

   
   ind=matrix(c(0,0,0, 1,0,0, 1,1,0,   #index for b_k(i,j)
	2,0,0,2,1,0, 2,0,1,2,1,1,2,2,0,
	3,0,1,3,1,1,3,0,2,4,0,2), byrow=T, ncol=3)
   txt1=c("Factors", "Lin", "Quad", "b_k(i,j)")
 
   worList=c()
   resList=c()

   for(TC in 1:noF) {
     word=matrix.choose.count(NM, noF, TC, myArray, noRun)
     resList = append(resList, list(word))
     worList = rbind(worList, word)
  }
  colnames(worList)=txt1

  ind.A=unlist(apply(worList[,-1][,c(1:NM)], 1, sum))
  ind.B=unlist(apply(worList[,-1][,c(1:NM)], 1, total.effect, NM))
  bList=worList[sort.list(ind.B),]
  rownames(bList)=ind.B[sort.list(ind.B)]
  word.A=as.matrix(unlist(lapply(c(1:noF), word.sum, ind.A, worList[,-1])))
  word.B=as.matrix(unlist(lapply(c(1:(NM*noF)), word.sum, ind.B, worList[,-1])))
 
  gwc=list(worList=worList, bList=bList, word.A=word.A, word.B=word.B)
 }
 if (noLev==2) { 
   bij=numeric(noF)
   for( TC in 1:noF)  {
     out=t(combn(noF,TC))
     bij[TC]=sum(apply(out, 1, proj.calc, d))/(noRun*noRun)
   }
   worList=rbind(1:noF, bij)
   txt1=c("Factors", "k", "b_k(k)")
   gwc=list(worList=worList)
 }
 
 gwc
}



cal.del<-function(noF,noRun, noLev) { 
 assign("Global.res", noRun, envir = .GlobalEnv)
 if (noF<=3) stop("number of factors is less than 3")
 del.mat=c()
 if (noLev==3) {
    ind=matrix(c(0,0,0, 1,0,0, 1,1,0,   #index for b_k(i,j)
	2,0,0,2,1,0, 2,0,1,2,1,1,2,2,0,
	3,0,1,3,1,1,3,0,2,4,0,2), byrow=T, ncol=3)
 
    for (i in 1:dim(ind)[1]){
       abc=ind[i,]
       ret2=mod.delta(abc, noF, noRun)
       del=c(abc, ret2)
       del.mat=rbind(del.mat, del)
    }
 }
 if (noLev==2) {
     ind=matrix(c(0,0,1,0, 2,0, 2,1,
                 3,1,3,2,4,2), byrow=T, ncol=2)

     for (i in 1:dim(ind)[1]){
       abc=ind[i,]
       ret2=delta2.int(abc, noF, noRun)
       del=c(abc, ret2)
       del.mat=rbind(del.mat, del)
    }

 }
 delta=del.mat
 delta
}



     
cal.Q3L<- function(d, type=T) {
#Type=F: calculate the Q criterion used the scaling in Tsai et al (2000)
  
   noRun=nrow(d)
   noF=ncol(d)
   noLev=3
   e1=rep(1, noRun) 
   NM=noLev-1
   myArray=array(NA, c(noRun,noF,NM))

  
   if (type==T) { 
     myArray[,,1] = d
     myArray[,,2] = sqrt(1/2)*(3 * myArray[,,1] ^2 - 2)
     myArray[,,1] = sqrt(3/2)*myArray[,,1]

     out1=t(combn(noF, 2))
      yi=c()
      for (i in 1:dim(out1)[1]) {
         yi=cbind(yi, apply(myArray[,,1][,out1[i,]],1,prod))  
       }

      U=cbind(e1, myArray[,,1], myArray[,,2], yi)
   } else {
     lin = d
     quad=3/2*lin^2-1
     out1=t(combn(noF, 2))
     yi=c()
     for (i in 1:dim(out1)[1]) {
        yi=cbind(yi, apply(lin[,out1[i,]],1,prod))       
     }
     U=cbind(e1, lin, quad, yi)
   }    
 
  
  
   ind=matrix(c(0,0,0, 1,0,0, 1,1,0, 
	2,0,0,2,1,0, 2,0,1,2,1,1,2,2,0,
	3,0,1,3,1,1,3,0,2,4,0,2), byrow=T, ncol=3)
 
   del.mat=c()
   for (i in 1:dim(ind)[1]){
      abc=ind[i,]
      ret2=mod.delta(abc, noF, noRun)
      del=c(abc, ret2)
      del.mat=rbind(del.mat, del)
   }
   W.mat=cal.Wmat(del.mat[,4], noF,noRun)
  

  xtx=t(U)%*%U
  sz2=nrow(xtx)
  a=matrix(rep(diag(xtx), sz2),sz2,sz2)
  rxtx=xtx*xtx/(a*a*t(a))

  qq=sum((rxtx*W.mat)[,-1])/(del.mat[1,4]-1)    
  names(qq)="Q"
  list(d=d,  Q=qq)
}


cal.Qb3L <- function(d, pp){
#d: a n*noF three-level design matrix with coding -1, 0, 1.
#pp: the prior values for linear, quadratic and 2-factor interaction  effects
#call pr.mat() and pr.adj for calculating the adjustment prior probability for 
#models being the best model 
# return xi.mat for \xi_{abc}

   noRun=nrow(d)
   noF=ncol(d)
   colnames(d)=LETTERS[1:noF]    #give A,B, .. as the factor names
   rownames(d)=c(1:noRun)       #numbering the runs
   noLev=3
   NM=noLev-1
   myArray=array(NA, c(noRun,noF,NM))

   myArray[,,1] = d
   myArray[,,2] = sqrt(1/2)*(3 * myArray[,,1] ^2 - 2)
   myArray[,,1] = sqrt(3/2)*myArray[,,1]

    ind=matrix(c(0,0,0, 1,0,0, 1,1,0, 
	2,0,0,2,1,0, 2,0,1,2,1,1,2,2,0,
	3,0,1,3,1,1,3,0,2,4,0,2), byrow=T, ncol=3)

   
 

  prmat=pr.mat(pp, noF, noRun) 
  pr.adj.c=pr.adj(prmat, noRun)
  pradj=pr.adj.c$pradj
  adj=pr.adj.c$adj
  del.mat=c()
  xi.mat=c()
  for (i in 1:dim(ind)[1]){
     abc=ind[i,]
      ret2=mod.delta(abc, noF, noRun)
      del=c(abc, ret2)
      del.mat=rbind(del.mat, del)
      ret=adj.xi(abc, adj, pp, noF, noRun)
      xi=c(abc, ret)
      xi.mat=rbind(xi.mat , xi)
  }
 
   P.mat=cal.Wmat(xi.mat[,4], noF, noRun)
   W.mat=cal.Wmat(del.mat[,4], noF, noRun)

    e1=rep(1, noRun)
    out1=t(combn(noF, 2))
    yi=c()
    for (i in 1:dim(out1)[1]) {
      yi=cbind(yi, apply(myArray[,,1][,out1[i,]],1,prod))    
     
    }

    U=cbind(e1, myArray[,,1], myArray[,,2], yi)
    xtx=t(U)%*%U
    sz2=nrow(xtx)
    a=matrix(rep(diag(xtx), sz2),sz2,sz2)
    rxtx=xtx*xtx/(a*a*t(a))
  
    qq=sum((rxtx*W.mat)[,-1])/(del.mat[1,4]-1)    

    qqb=sum((rxtx*P.mat)[,-1])

   (list(d=d,  Q=qq, Q.b=qqb, xi.mat=xi.mat ))
}



cal.Db3L <- function(d, pp, tau){

#d: a n*noF three-level design matrix with coding 1 2 3
#tau: the turning constant for the primary and potential effects in Bayeain D-optimal design

  if (sum(duplicated(pp))==0 | sum(duplicated(pp))==3  ){ 
     D.b=0
  } else {  #with some same prior

    noRun=nrow(d)
    noF=ncol(d)
    e1=rep(1, noRun)
    x.lin=d
    x.quad=x.lin^2-2/3
    out1=t(combn(noF, 2))
    x.int=c()
    for (i in 1:dim(out1)[1]) {
     x.int=cbind(x.int, apply(x.lin[,out1[i,]],1,prod))         
    }
    x.int=x.int/2

    if (pp[1]==pp[2]) {
     x.prime=cbind(e1, x.lin, x.quad)
     x.pot=cbind(x.int) 
    } 

    if (pp[2]==pp[3]) {
     x.prime=cbind(e1, x.lin)
     x.pot=cbind(x.quad, x.int)
    }

    if (pp[1]==pp[3]) {
     x.prime=cbind(e1, x.lin, x.int)
     x.pot=cbind(x.quad)
    }

      no.prime=ncol(x.prime)
      no.pot=ncol(x.pot)


      U=cbind(x.prime, x.pot)       
      xtx=t(U)%*%U
      K=rbind(cbind(diag(0,no.prime), matrix(0, nrow=no.prime, ncol=no.pot)), 
            cbind(matrix(0,  nrow=no.pot, ncol=no.prime), diag(1, no.pot)))
      Z.star=xtx+K/(tau^2)
      D.b=log(det(Z.star))

   }
 D.b
}


cal.Q2L = function(x) {

  noCol=ncol(x)
  noRun=nrow(x) 
  b.word=numeric(noCol)
  for( TC in 1:noCol)  {
    out=t(combn(noCol,TC))
    b.word[TC]=sum(apply(out, 1, proj.calc, x)/(noRun*noRun))
  }

  
  Q.1st=( b.word[1] * num.mod2(2, noCol, noRun)+
	  b.word[2] * 2 * num.mod2(3, noCol, noRun))/num.mod2(1, noCol, noRun)

  if (noCol==3) { 
      Q.2nd=( b.word[1] * ( delta2.int(c(1,0), noCol, noRun) + 
		2 * (noCol-1) * delta2.int(c(2,1), noCol, noRun) )  + 
	b.word[2] * ( 2 * delta2.int(c(2,0), noCol, noRun) + 
		delta2.int(c(2,1), noCol, noRun) + 
		2 * (noCol-2)*delta2.int(c(3,2), noCol, noRun)) + 
	b.word[3] * 6 * delta2.int(c(3,1), noCol, noRun) )/delta2.int(c(0,0), noCol, noRun)
   } else {
   Q.2nd=(b.word[1] * (delta2.int(c(1,0), noCol, noRun) + 
		2*(noCol-1) * delta2.int(c(2,1), noCol, noRun)  )+ 
	b.word[2] * ( 2 * delta2.int(c(2,0), noCol, noRun) + 
		delta2.int(c(2,1), noCol, noRun) + 
		2 * (noCol-2)*delta2.int(c(3,2), noCol, noRun) ) +
	b.word[3] * 6 * delta2.int(c(3,1), noCol, noRun) + 
	b.word[4] * 6 * delta2.int(c(4,2), noCol, noRun) )/delta2.int(c(0,0), noCol, noRun)
	}

   c( Q.2nd)
}

#probability of a model with a particular set of a.val main effects and c.val being the fitting model
prb.calc=function(c.val, a.val, c.max, noF, p1, p2){
  if (c.val<=c.max) {
  prb.calc=p1^a.val*(1-p1)^(noF-a.val)*p2^c.val*(1-p2)^(c.max-c.val) }
 else {
   prb.calc=0}
}


prb=function(ac, noF, p1, p2){
#probability of  a model with a.val main effects and  each of 0 to c.max  interactions being the fitting model
#calling prb function of real calculation
  a.val=ac[1]
  c.max=ac[2]
  unlist(lapply(c(0:c.max), prb.calc, a.val, c.max, noF, p1, p2))
}
  


prob.lev2=function(p1, p2, noF, noRun){
#a list of probability for model with i-1 main effects and 0 to c.max interaction being the fitting model (without adjustment)

  xi.size=7
  xi=numeric(xi.size)
  xi.ind=matrix(c(0,1,2,2,3,3,4,
	0,0,0,1,1,2,2), 2, 7, byrow=T)
  pij=c()
  a.ind=c(0:noF)
  c.ind=choose(a.ind,2)
  prob=apply(rbind(a.ind,c.ind), 2, prb, noF, p1, p2)
  mod=apply(rbind(a.ind, c.ind), 2, mod.num, 0, 0, noF)
  delta=apply(xi.ind, 2, delta.ij, noF, noRun)

  #mod is different from delta2.int in Q2.r since here non-eligible models still count
  #mod[[i]][c.ind] the number of model with (i-1) main effects and (c.ind-1) interaction

  A=matrix(0,0,0)
  for (i in a.ind) {
     prb.m=mod[[i+1]]*prob[[i+1]]
     A=append(A, list(prb.m))
  }
  
  tmp.prob=0
  tmp.num=0
  for (i in 1:noF) {
    for (j in 0:c.ind[i+1]) {
      if ((i+j)>=noRun) {
         tmp.prob=tmp.prob+A[[i+1]][j+1]
      }
     if ((i+j)==(noRun-1))  {
      tmp.num=tmp.num+mod[[i+1]][j+1]}
     }
  }

  prob.adj=tmp.prob/tmp.num

  for (i in 1:noF) {
    for (j in 0:c.ind[i+1]) {
      if ((i+j)>=noRun) {
          prob[[i+1]][j+1]=0
      }
     if ((i+j)==(noRun-1))  {
             prob[[i+1]][j+1]= prob[[i+1]][j+1]+prob.adj
	}
     }
  } #adjust gamma to model of a size smaller.

  A=matrix(0,0,0)
  for (i in a.ind) {
    prb.m=mod[[i+1]]*prob[[i+1]]
    A=append(A, list(prb.m))
  }
  
  for (i.ind in 1:xi.size){
    a.min=xi.ind[1,i.ind]
    c.min=xi.ind[2,i.ind]
    mod.ad=apply(rbind(a.ind, c.ind), 2, mod.num, a.min, c.min, noF)
    A.ad=matrix(0,0,0)
    for (i in a.ind) {
      prb.m=mod.ad[[i+1]]*prob[[i+1]]
      A.ad=append(A.ad, list(prb.m))
     }  
    pij=rbind(pij, sum(unlist(A.ad),na.rm=TRUE))
  }
  return(pij)
}



cal.Qb2L=function(p1,p2,d){
  noCol=ncol(d)
  noRun=nrow(d)
  b.word=numeric(noCol)
  for( TC in 1:noCol)  {
    out=t(combn(noCol,TC))
    b.word[TC]=sum(apply(out, 1, proj.calc, d))/(noRun*noRun)
  }

  pij=prob.lev2(p1, p2, noCol, noRun)

  if (noCol==3) {
     Q.B=( b.word[1]*( pij[2] + 2*(noF-1)*pij[4] )+ 
	b.word[2] * ( 2*pij[3] + pij[4]+ 2*(noF-2)*pij[6] ) +
	b.word[3] * 6 * pij[5] ) / pij[1]
   } else {
   Q.B =( b.word[1]*( pij[2] + 2*(noF-1)*pij[4] )+ 
	b.word[2] * ( 2*pij[3] + pij[4]+ 2*(noF-2)*pij[6] ) +
	b.word[3] * 6 * pij[5] + b.word[4] * 6 * pij[7] ) / pij[1]
   } 

 return(list(Q.b=Q.B, b.word=b.word, pij=pij))
}