R/utilcat.R

In FrF2: Fractional Factorial Designs with 2-Level Factors

FrF2.currentlychecked <- function() print(getFrF2(".FrF2.currentlychecked"))

mult.gen <- function(liste){
   ## omitted error check, because slows down very much
   ## must be applied to absolute values of generators!!!
   if (is.list(liste)) hilf <- table(unlist(liste))%%2
   else hilf <- table(liste)%%2
   as.numeric(names(hilf[hilf==1]))
}

mult.gen.a <- function(string.vec,sep=""){
   ## omitted error check, because slows down very much
   ## must be applied to absolute values of generators!!!
   hilf <- table(unlist(strsplit(string.vec,sep)))%%2
   paste(names(hilf[hilf==1]),collapse=sep)
}


words.all <- function(k, gen, design=NULL, max.length=7, select.catlg=catlg){
   ## 2^k is the run number of the design
   ## gen is a list of vectors which generate additional factors from interactions among the first k ones 
   ##   column numbers can be incorporated by using Yates[coln], e.g. Yates[c(7,31)] for gen=list(c(1,2,3), c(1,2,3,4,5))
   ##   but now also directly (gen = vector of column numbers)
   ## max.length is the maximum length of words shown (calculations do become faster by reducing this)
   ##      can be set to NULL in order to use maximum possible word length (may be advisable for use in further calculations)
   if (!(k>0 & floor(k)==k)) stop("k must be a positive integer number.")
   
   if (is.character(design)) gen <- select.catlg[[design]]$gen
   ### check generator entry and transform to list of vector of numbers
   gen <- gen.check(k,gen)
   
   minus <- which(sapply(gen, function(obj) all(obj<0)))
   gen <- lapply(gen, "abs")

   g <- length(gen)
   words <- as.list(rep(0,2^min(g,max.length)-1))
   ## append the letter for factor that is generated by generator
   gen <- mapply(c,gen,as.list((k+1):(k+g)),SIMPLIFY=FALSE)
   if (is.null(max.length)) max.length <- k+g
   words[1:g] <- gen
   words[minus] <- lapply(words[minus],"-")
   if (g>1) {
      cur <- g
      if (max.length >= floor((g+1)/2))
           hilf <- combn(g,floor((g+1)/2))  
      else hilf <- combn(g,max.length)      ## for early breakdown in case of storage problems
      for (i in 2:min(g,max.length)){
            hilf <- combn(g,i)
            for (j in 1:ncol(hilf)){
               cur <- cur+1
               words[[cur]] <- mult.gen(gen[hilf[,j]])
               if (length(which(hilf[,j] %in% minus))%%2 == 1) words[[cur]] <- -words[[cur]]
               }
            }
   }
   if (g<=max.length) words[[2^g-1]] <- mult.gen(gen)
   if (length(minus)%%2==1) words[[2^g-1]] <-  - words[[2^g-1]] 
   WLP <- sapply(words, FUN=length)
      sellang <- WLP<=max.length
   WLP <- WLP[sellang]
   wl <- list(table(WLP),words[sellang][order(WLP)])
   wl[[2]] <- wl[[2]][1:sum(wl[[1]][as.numeric(names(wl[[1]]))<=max.length])]
   wl[[3]] <- max.length

   names(wl) <- c("WLP", paste("words.up.to.length",max.length,sep="."),"max.length")
   class(wl) <- c("list","wordlist")
   wl
}

##wl <- words.all(5,list(c(1,2),c(1,3,4),c(1,3,5),c(1,4,5),c(2,3,4,5)))

words.four <- function(wl, f4, max.length=wl$max.length){
### !!! not adapted to the new possibility of negative sign generators!!!
   if (!inherits(wl,"wordlist")) stop("wl must be a wordlist object.")
   if (!(is.vector(f4) | is.list(f4))) stop("f4 must be a vector of defining factors for 4-level factors or a list of such vectors.")
   if (!is.list(f4) & !length(f4) %in% c(2,3)) stop("The vector f4 must be of length 2 or 3. For more than one 4-level factor, use a list of vectors.")
   if (!is.list(f4)) f4 <- list(f4)
   if (as.numeric(names(wl[[1]])[1])<3) stop("This design confounds main effects with main effects. DO NOT USE IT!") 
   else if (!names(wl[[1]])[1]=="3") stop("Function words.four requires that the third df for the four level factor is also included as a main effect.")

   if (wl[[3]]<max.length) stop("Requested max.length not possible, because object wl has smaller max.length")

   wl[[2]]<-wl[[2]][1:sum(wl[[1]][as.numeric(names(wl[[1]]))<=max.length])]

   drei <- wl[[2]][1:wl[[1]][1]]  ## select three-letter words
   ## find three letter words that correspond to factors in f4
      hilf <- sapply(f4, function(obj) {pos <- which(sapply(drei, function(obj2) all(obj %in% obj2))) 
                                        if (length(pos)==0) pos <- 0 
                                        pos} )
      if (any(hilf==0)) stop(paste("f4 specifications in positions",which(hilf==0),"do not match three-letter words in the design."))
      f4 <- drei[hilf]            ## now all f4 elements have all three entries
   ## reduce word list in order to remove superfluous words
      for (i in 1:length(f4)){
             hilf <- sapply(wl[[2]], function(obj) f4[[i]] %in% obj)
             wl[[2]] <- wl[[2]][colSums(hilf)<2]   ## omit words that are already covered by presence of third factor
          }
     WLP4 <- sapply(wl[[2]],FUN=length)
      WLP4 <- table(WLP4)
      wl[[4]] <- f4
      wl[[5]] <- WLP4
      wl[[3]] <- max.length
      names(wl) <- c("WLP.orig", paste("words.reduced.up.to.length",max.length,sep="."), "max.length", "4.level.generators", "WLP4")
      class(wl) <- c("list","wordlist4")
      wl
}

#words.four(wl, c(1,2,6))

print.wordlist <- function(x, ...){
   if (!inherits(x,"wordlist")) stop("Function print.wordlist prints wordlist objects only.")
   print(x[[1]])
   cat("\n\nWord list (up to length ",x$max.length,"):\n",sep="")
   if (max(unlist(x[[2]]))>50) print(sapply(x[[2]],function(obj) 
                          paste(sign(obj[1]),"(",paste(abs(obj),collapse=","),")",sep=""),quote=FALSE))
   else print(sapply(x[[2]],function(obj) paste(if (obj[1]<0) "-" else "",paste(Letters[abs(obj)],collapse=""),sep="")),quote=FALSE)
}

print.wordlist4 <- function(x, ...){
### !!! not adapted to the new possibility of negative sign generators!!!
   if (!inherits(x,"wordlist4")) stop("Function print.wordlist4 prints wordlist4 objects only.")
   print(x[[1]])
   cat("\n\nWord list with 4-level factor internal words removed:\n")
   if (max(unlist(x[[2]]))>50) print(sapply(x[[2]],function(obj) paste("(",paste(obj,collapse=","),")",sep=""),quote=FALSE))
   else print(sapply(x[[2]],function(obj) paste(Letters[obj],collapse="")),quote=FALSE)
   cat("\n\nWord length pattern for reduced word list:\n")
   print(x[[5]])
   cat("\n\nContrast combinations for 4-level factors:\n")
   if (max(unlist(x[[4]]))>50) print(sapply(x[[4]],function(obj) paste("(",paste(obj,collapse=","),")",sep=""),quote=FALSE))
   else print(sapply(x[[4]],FUN=function(obj) paste(Letters[obj],collapse=""),simplify=FALSE),quote=FALSE)
}

alias3fi <- function(k, gen, order=3){
   ## k number of factors spanning basis full factorial
   ## gen list of vectors giving generators (in terms of position of full factorial factors, 
   ##       e.g. list element c(1,2) is AB
   if (!order %in% c(2,3)) stop("order must be 2 or 3")
   gen <- gen.check(k, gen)
   if (!(k>0 & floor(k)==k)) stop("k must be a positive integer number.")
   if (!is.list(gen)) stop("gen must be a list of generator vectors.")
   g <- length(gen)  ## number of generators
   struc <- NULL
   sep <- ""

   if (k + g > 50) {
      Letters <- paste("F",1:(k+g),sep="")
      sep <- ":"
      }

   minus <- c(rep("",k),sapply(gen, function(obj) if (all(obj<0)) "-" else ""))
   gen <- lapply(gen, "abs")
   
   ## vector of main effects in terms of generators in character form
   absall1 <- sapply(c(as.list(1:k),gen),function(obj) paste(Letters[obj],collapse=sep))
   all1 <- paste(minus,absall1,sep="")
   names(all1) <- Letters[1:(g+k)]

   ## vector of 2fis in terms of generators in character form
   all2 <- rep(0,choose(k+g,2))
   sel <- combn(k+g,2)
   n2 <- ncol(sel)
   for (i in 1:n2) all2[i] <- mult.gen.a(all1[sel[,i]],sep=sep)
   names(all2) <- apply(sel,2,function(obj) paste(Letters[obj],collapse=sep))

   all3 <- NULL
   if (order==3){
   ## vector of 3fis in terms of generators in character form
   all3 <- rep(0,choose(k+g,3))
   sel <- combn(k+g,3)
   n3 <- ncol(sel)
   for (i in 1:n3) all3[i] <- mult.gen.a(all1[sel[,i]],sep=sep)
       names(all3) <- apply(sel,2,function(obj) paste(Letters[obj],collapse=sep))}

   ## pattern contains number of occurrences of each effect among main to 3fis
   ## or 2fis only in case of order=2
   all3 <- c(all1,all2,all3)
   pattern <- table(absall3 <- sub("-","",all3))
       ## now all3 contains all lengths
   if (max(pattern)==1){
          if (order==2)
          aus <- "no aliasing among main effects and 2fis"
          else
          aus <- "no aliasing of main effects or 2fis with effects up to order 3"
   }
   else {
          struc <- as.list(rep(0,length(pattern)))
          for (i in 1:length(pattern)){ 
               ## make sure that alias information apears in systematic order
                  struc[[i]] <- c(names(all3)[all3==names(pattern)[i]], 
                        paste("-",names(all3)[which(all3==paste("-",names(pattern)[i],sep=""))],sep=""))
                  struc[[i]] <- struc[[i]][sort(gsub("-","",struc[[i]]),index.return=TRUE)$ix]
                  hilf <- gsub("-","",struc[[i]])
                  struc[[i]] <- struc[[i]][which(nchar(hilf)>0)]
                  hilf <- hilf[which(nchar(hilf)>0)]
                  hilf <- sort(nchar(hilf),index.return=TRUE)$ix
                  struc[[i]] <- struc[[i]][hilf]
                  }
          struc <- sapply(struc[which(sapply(struc,length)>1)], "paste",collapse="=")
          ## make every equation start positive
          struc <- sapply(struc, function(obj){
               if (length(grep("^-",obj))==1){
                   obj <- gsub("-","~",obj)
                   obj <- gsub("=","=-",obj)
                   obj <- gsub("=-~","=",obj)
                   obj <- gsub("~","",obj)
               } 
               obj
               })
          names(struc) <- NULL
   ## sort in ascending order
   ## with main effect
   if (k+g<=50){
   wme <- grep("^[[:alpha:]]=[[:alpha:][:punct:]]*",struc)
   ## with 2fi
   wme2 <- grep("^[[:alpha:]]{2}=[[:alpha:][:punct:]]*",struc)
   ## with 3fi
   if (order==3) wme3 <- grep("^[[:alpha:]]{3}=[[:alpha:][:punct:]]*",struc)
   ## sort order not yet good for more than 25 factors (a before A etc., and also unfortunately locale dependent)
   ## !!! improve by making sort depend on factor number rather than letter
   }
   else{
   wme <- grep("^F[[:digit:]]+=F[[:digit:][:punct:]]*",struc)
   ## with 2fi
   wme2 <- grep("^F[[:digit:]]+:F[[:digit:]]+=F[[:digit:][:punct:]]*",struc)
   ## with 3fi
   if (order==3) wme3 <- grep("^F[[:digit:]]+:F[[:digit:]]+:F[[:digit:]]+=F[[:digit:][:punct:]]*",struc)
   }
   if (order==2) aus <- list(main=sort(struc[wme]),fi2=sort(struc[wme2]))
   else aus <- list(main=sort(struc[wme]),fi2=sort(struc[wme2]),fi3=sort(struc[wme3]))
   }
   aus
   }

   gen.check <- function(k,gen){   if (!is.list(gen)) {
                 if (!(is.numeric(gen) | is.character(gen))) 
                      stop("gen must be a list of generator vectors, a vector of column numbers or a character vector of generators.")
             ## August 2019: handle full factorial
                 if (length(gen)==0) return(vector(mode="list"))
                 if (is.character(gen)){
                     absgen <- sub("-","",gen)
                     minus <- grep("^-",gen)
                     gen <- lapply(strsplit(absgen,""), function(obj) which(Letters %in% obj))
                 }
                 else {
                 ## numeric gen
                 absgen <- abs(gen)
                 minus <- which(gen<0)
                 if (any(!gen==floor(gen))) stop("All entries in gen must be integer numbers.")
                 if (any(2^(0:(k-1)) %in% absgen)) stop("This design is of resolution II and is not permitted in package FrF2.")
                 if (min(absgen)<1 | max(absgen)>2^k-1) stop("Column numbers in gen must be in the range of 3 to 2^k-1.")
                 ## gen <- Yates[absgen]}
                 ## usage of Yates replaced by indexcalc in order to work for larger designs
                 gen <- indexcalc(absgen)}
                 gen[minus] <- lapply(gen[minus],"-")
              }
              if (any(sapply(gen,function(obj) any(abs(obj)<1 | obj>k | !floor(obj)==obj) )))
                   stop(paste("All generators must contain integer numbers from 1 to", k, 
                     "\n or letters from",Letters[1],"to", Letters[k], "only."))
              gen
              }
 
   slowfast <- function(k){
   ## yields a permutation vector 
   ## for switching back and forth between slow and fast ordering
   ## of the base factors of a design
   ind <- 0:(2^k-1)
   ord(as.data.frame(lapply(1:k, function(obj) ind%%(2^obj))))
}