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)
}
