Tuesday, July 24, 2012

Williams designs with 5 products

In a previous post I created small Williams designs for an even number of products. This worked very well, also because the number of permutations could be restricted significantly due to symmetry. Unfortunately this does not work so well with an odd number of products. The symmetry is not so obvious. In fact, using brute force I can find non-symmetric solutions for five products. However, I am lacking force to handle more than five products in a reasonable time.

design setup

The design fixed parameters

To get the design I used a few assumptions. The number of rows is twice the number of columns. This is because I remember reading there are no solutions with the same number of rows as columns. The first row reads 1:n. The first column is all elements of 1:n twice.

Filling in the design

The approach is actually very simple.
  • Look at the first not yet know point. 
  • Examine which value it could have. 
    • not yet used in current row
    • at most once used in current column
    • at maximum once used after the product in current row, one column to the left (carry over)
    • if this is the second column and the first column of this row is the same as the first column of the previous row, then the current value must be higher than the value of one row up.
  • Try the possible values in a loop
  • Within that loop go back to step one
Implemented recursively. In all actually, it is quite brute force. I should be doing this in C or Julia rather than in R which is relatively slow. 

Tricks

Compared to the version with an even number of products I implemented some things a bit smarter. Removed some if within a for loop by selecting the objects over which the loop sequences. I also noticed some things which surprised me in making things faster. For instance, creating temporary variables (e.g. row <- desconst$todo[desobject$count,2L]) is faster than using desconst$todo[desobject$count,2L] to index matrices later on. Merging two simple lines of code into one complex line of code is slower. Apparently it is cheaper to store a small result than have the interpreter figure out a complex line of code. Hence I used quite a lot of local variables, certainly more than I find 'nice'. 

Designs

The result if 90 Williams designs for five products. Without doubt there is quite some duplication from permutations of the objects. The designs can be separated in three classes. Designs symmetric in columns (72 of those), designs with four rows symmetric in columns (five of those) and designs where are not symmetric in columns at all (13 of those).

design symmetric in columns

The first design has five unique rows each copied reversed to make the ten rows. For ease I colored the matching rows.
      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    2    3    4    5
 [2,]    1    3    2    5    4
 [3,]    2    1    4    3    5
 [4,]    2    4    1    5    3
 [5,]    3    1    5    2    4
 [6,]    3    5    1    4    2
 [7,]    4    2    5    1    3
 [8,]    4    5    2    3    1
 [9,]    5    3    4    1    2
[10,]    5    4    3    2    1

designs four rows symmetric in columns

The colored rows have a reversed row in the same color. 

      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    2    3    4    5
 [2,]    1    3    4    2    5
 [3,]    2    3    5    1    4
 [4,]    2    4    1    5    3
 [5,]    3    1    5    4    2
 [6,]    3    5    2    1    4
 [7,]    4    1    2    5    3
 [8,]    4    5    1    3    2
 [9,]    5    2    4    3    1
[10,]    5    4    3    2    1

design not symmetric in columns


      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    2    3    4    5
 [2,]    1    3    5    2    4
 [3,]    2    4    1    5    3
 [4,]    2    5    3    4    1
 [5,]    3    1    4    2    5
 [6,]    3    2    1    5    4
 [7,]    4    3    5    1    2
 [8,]    4    5    2    1    3
 [9,]    5    1    4    3    2
[10,]    5    4    2    3    1

R code

generation

gendesign <- function(n=3) {
nc <- as.integer(n)
nr <- nc+nc
desmat <- matrix(NA,nrow=nr,ncol=nc)
desmat[1,] <- 1L:nc
desmat[,1] <- rep(1L:nc,each=2)
carover <- matrix(0L,nrow=nr,ncol=nc)
for (i in 1L:(nc-1L)) carover[i,i+1]  <- 1L
todo <- which(is.na(t(desmat)),arr.ind=TRUE)
desobject <- list(desmat=desmat,carover = carover,count=0)
desconst <- list(todo=todo,ntodo=nrow(todo),nc=nc,totest = 1L:nc)
desresult <- list()
addpoint(desobject,desresult,desconst)
}

addpoint <- function(desobject,desresult,desconst) {
desobject$count <- desobject$count+1L
if (desobject$count > desconst$ntodo) {
l <- length(desresult)
desresult[[l+1]] <- desobject$desmat
return(desresult)
}
row <- desconst$todo[desobject$count,2L]
col <- desconst$todo[desobject$count,1L]
dob <- desobject
currow <- desobject$desmat[row,]
currow <- currow[!is.na(currow)]
totest <- desconst$totest[-currow]
prev <- desobject$desmat[row,col-1L]
totest <- totest[desobject$carover[prev,totest]<2L]
counts <- rowSums(outer(totest,desobject$desmat[,col],function(x,y) as.integer(x==y)),na.rm=TRUE)
totest <- totest[counts<2L]
if (col==2 & desobject$desmat[row,1]==desobject$desmat[row-1L,1L]) {
totest <- totest[totest > (desobject$desmat[row-1L,2])]
}
for (i in totest) {
desobject$carover[prev,i] <- desobject$carover[prev,i] + 1L
desobject$desmat[row,col] <- i
desresult <- addpoint(desobject,desresult,desconst)
desobject <- dob
}
desresult
}

g5 <- gendesign(5)

checking

# check that each column has each object twice
sapply(g5,function(sub) all(sapply(1:5,function(x)   
      all(table(sub[,x])==rep(2,5)))))

# check that each product is followed by each other product twice
sub <- g5[[1]]
target <- Reduce('+',lapply(1:4,function(x) table(sub[,x],sub[,x+1])))
target  # shows a symmetric matrix with 0 on diagonal, 2 elsewhere.
sapply(g5,function(sub) identical(target,
      Reduce('+',lapply(1:4,function(x) table(sub[,x],sub[,x+1])))))

#check the number of reversed rows
sa <- sapply(g5,function(sub1) {
10-sum(1:10 %in% sapply(1:10,function(i) {
r1 <- rev(sub1[i,])
mat <- matrix(rep(r1,each=10),nrow=10,ncol=5)
which(rowSums(sub1==mat)==5)
}))
})
table(sa)
# number of reversed rows
# 0  4 10 
#13  5 72 

prints

g5[[1]]
g5[[4]]
g5[[8]]

Thursday, July 12, 2012

Creating Williams designs with even number of products

A Williams design is a special Latin square with the additional property of first order carry over (each product is followed equally often by each other product). In R the package crossdes can be used to create them.

> williams(4)
     [,1] [,2] [,3] [,4]
[1,]    1    2    4    3
[2,]    2    3    1    4
[3,]    3    4    2    1
[4,]    4    1    3    2
As a consequence of the carry over restriction, the design has the property that a row in this design is also reversed present in the design. Example, row 3 is row 1 reversed. For small designs this property can be used to generate the designs by brute force. Example; a four by four design has without loss of generality the first row and column designated 1 to 4 (using . as unknown).
1 2 3 4
2 . . .
3 . . .
4 . . .
Adding the reversal of the first row gives:
1 2 3 4
2 . . .
3 . . .
4 3 2 1
In practice, for an even number of products the last column can be created as reversed first column.
1 2 3 4
2 . . 3
3 . . 2
4 3 2 1
It is rather obvious that only one solution remains for the 2,2 location and the rest is simple filling in the blanks. The resulting design is the same as the solution of williams(4) with '3' and '4' permuted.
1 2 3 4
2 4 . 3
3 . . 2
4 3 2 1
It is possible to use the same approach for small designs with an even number of products. For an odd number of products the number of rows has to be doubled so this will follow in a later post. To create the design with 6 products a program it is more convenient than manual work. It appears, unknown to many, that using this approach there are two solutions for 6 products. These two solutions are not permutations of each other!
[[1]]
     [,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

[[2]]
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    4    6    1    3    5
[3,]    3    6    2    5    1    4
[4,]    4    1    5    2    6    3
[5,]    5    3    1    6    4    2
[6,]    6    5    4    3    2    1

There are 8 solutions for 8 products and 192 solutions for 10 products. Above that, calculation time is prohibitive. Even to get to 10, the program had to be streamlined considerably.

R code

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

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

checkdes <- function(desobject,row,col) {
# test for only once carry over 
all(desobject$carover<=1) & 
# each product only once in each row and each column
!any(desobject$desmat[row,-col]== desobject$desmat[row,col],na.rm=TRUE )  &
!any(desobject$desmat[-row,col]== desobject$desmat[row,col],na.rm=TRUE )  
}

addpoint <- function(desobject,desresult) {
todo <- nextpos(desobject$desmat)
if (length(todo)==0) {
l <- length(desresult)
desresult[[l+1]] <- desobject$desmat
return(desresult)
row <- todo[1L,1L]
col <- todo[1L,2L]
nc <- ncol(desobject$desmat)
dob <- desobject
for (i in 1L:nc) {
desobject$desmat[row,col] <- i
desobject$desmat[nc-row+1,nc-col+1] <- i
desobject$carover[desobject$desmat[row,col-1L],i] <- desobject$carover[desobject$desmat[row,col-1],i] + 1L
desobject$carover[i,desobject$desmat[row,col-1L]] <- desobject$carover[i,desobject$desmat[row,col-1L]] + 1L
other <- desobject$desmat[row,col+1L]
if (!is.na(other)) {
desobject$carover[other,i] <- desobject$carover[other,i] + 1
desobject$carover[i,other] <- desobject$carover[i,other] + 1
}
if (checkdes(desobject,row,col)) desresult <- addpoint(desobject,desresult)
desobject <- dob
}
desresult
}
gendesign(6)