Sunday, October 27, 2013

Symmetry in Williams designs

Based on some comments I am looking at using symmetry to obtain Williams style designs. Symmetry allows reduction of the number of combinations examined hence faster calculation times. Two avenues are examined. Both work for a low number of treatments. For eight treatments they give different sets, which together define all 12 possible designs.

Horizontal reflection

It is well known that each row in an even number of treatments Williams design is also present reversed. But apn noticed something else, within a row there is a horizontal symmetry; the values of the cell i and n+1-i add up to n+1. That is a bit abstract. Lets look at a six treatment design. Cell [2,2] is 4 and cell [2,5] is 3, they add to 7. And so it continues. If this holds, then the permutations to check is decreased enormously. Especially, since cell [2,2] also by the reversal of rows is equal to cell [5,5] and hence cell [5,2] is 3.
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    4    1    6    3    5
[3,]    3    1    5    2    6    4
[4,]    4    6    2    5    1    3
[5,]    5    3    6    1    4    2
[6,]    6    5    4    3    2    1
After putting this in an algorithm it appeared it worked. Run time for 8 treatments only 100 miliseconds.
microbenchmark(gendesign(4),gendesign(6),gendesign(8))
Unit: milliseconds
         expr        min         lq     median         uq        max neval
 gendesign(4)   1.576118   1.609474   1.638796   1.717969   2.896391   100
 gendesign(6)   8.695041   8.876843   9.054615   9.229821  14.714345   100
 gendesign(8) 100.170001 102.489826 106.144955 108.703397 121.015079   100 
The algorithm gave 8 designs for 8 treatments, 192 for 10 treatments, same as previously.

Symmetrical

I am also looking for a smart way to tackle designs for odd products. There is a nine treatment solution, I want to be able to find it. On a whim, I decided to check if symmetrical designs give is the solution. It is not. However, this approach can also be used for the even number of treatments. Again there are 8 designs for 8 treatments. Four of these have the horizontal semi-reflection, four do not have it. Hence all 12 designs which I found previously by brute force are covered. With 13 seconds run time for 8 treatments, I did not check for the number of solutions for 10 treatments..
Note 1: I do know symmetrical has no meaning in a Williams design. I have only used the designs with rows as subjects, columns as periods/rounds and numbers as treatments. And I do randomize subjects over rows, treatments over numbers. Symmetry is thus just a means.
Note 2: The 9 treatment design? It was found via group theory. Apn gave a link Some New Row-Complete Latin Squares, Archdeacon, Dinitz, Stinson and Tillson

Code

The two sets of code are quite similar. This is only logical since one is derived from the other. Each section should be able to be run separately as displayed here.

Horizontal semi reflection

nextpos <- function(desmat) which(desmat==0,arr.ind=TRUE)

gendesign <- function(n=6) {
  nr <- as.integer(n)
  nc <- nr
  
  desmat <- matrix(0L,nrow=nr,ncol=nc)
  desmat[1,] <- 1L:nc
  desmat[,1] <- 1L:nr
  desmat[n,] <- nc:1L
  desmat[,n] <- nr:1L
  desobject <- list(desmat=desmat,nc=1L:n,n=n,np=n+1L)
  desresult <- list()
  addpoint(desobject,desresult)
}

modify <- function(desobject,row,col,i) {
  desobject$desmat[desobject$np-row,desobject$np-col] <- 
      desobject$desmat[row,col] <- i
  desobject$desmat[desobject$np-row,col] <- 
      desobject$desmat[row,desobject$np-col] <- desobject$np-i
  desobject}

addpoint <- function(desobject,desresult) {
  if (!cocheck(desobject)) return(desresult)
  rc <- nextpos(desobject$desmat)
  if (length(rc)==0) {
    if (var(colSums(desobject$desmat))==0) {
      l <- length(desresult)
      desresult[[l+1]] <- desobject$desmat
  }
    return(desresult)
  }
  row <- rc[1,1]
  col <- rc[1,2]
  avoid <- c(desobject$desmat[row,],
      desobject$desmat[,col])
  nc <- desobject$nc[!is.element(desobject$nc,avoid) ]
  for (i in nc) {
       desresult <- addpoint(modify(desobject,row,col,i)
       ,desresult)
  }
  desresult
}

cocheck <- function(desobject) {
  dm2 <- desobject$desmat*desobject$n
  dm2 <- dm2[,-1]
  dm3 <- desobject$desmat[,-desobject$n]
  dm4 <- dm2+dm3
  dm4 <- c(dm4[dm2!=0 & dm3!=0])  
  all(table(dm4)==1)
}

symmetrical

nextpos <- function(desmat) which(desmat==0,arr.ind=TRUE)

gendesign <- function(n=5) {
  nr <- as.integer(n)
  nc <- nr
  desmat <- matrix(0L,nrow=nr,ncol=nc)
  desmat[1,] <- 1L:nc
  desmat[,1] <- 1L:nr
  desobject <- list(desmat=desmat,nc=1L:n,n=n,np=n+1L)
  desresult <- list()
  addpoint(desobject,desresult)
}

modify <- function(desobject,row,col,i) {
  desobject$desmat[col,row] <- 
      desobject$desmat[row,col] <- i
  desobject}

addpoint <- function(desobject,desresult) {
  if (!cocheck(desobject)) return(desresult)
  rc <- nextpos(desobject$desmat)
  if (length(rc)==0) {
    if (var(colSums(desobject$desmat))==0) {
      l <- length(desresult)
      desresult[[l+1]] <- desobject$desmat
    }
    return(desresult)
  }
  row <- rc[1,1]
  col <- rc[1,2]
  avoid <- c(desobject$desmat[row,],
      desobject$desmat[,col])
  nc <- desobject$nc[!is.element(desobject$nc,avoid) ]
  for (i in nc) {
    desresult <- addpoint(modify(desobject,row,col,i)
        ,desresult)
  }
  desresult
}

cocheck <- function(desobject) {
  dm2 <- desobject$desmat[,-1]*desobject$n
  dm3 <- desobject$desmat[,-desobject$n]
  dm4 <- dm2+dm3
  dm4 <- c(dm4[dm2!=0 & dm3!=0])  
  all(table(dm4)==1)
}

No comments:

Post a Comment