Wiekvoet: Carry-over balanced designs for 8 treatments

Those are Williams designs you might say, but it has become clear to me that Williams designs are just a subset of all carry-over balanced designs. Not through hard work of mine, comments by Apn on my previous post Creating Williams designs with even number of products lead to this. Among other things, Apn claimed there were designs which did not confirm to the well known symmetry that each row is the reflection of another row. Related to that was a claim that there is a nine treatment square carry-over balanced design. I made my algorithm faster and in this post I can confirm the former, but the latter is too long a calculation for the algorithm I used.

Designs

There are 12 designs, as given below. Numbers 4 and 5 can be obtained from each other by certain recoding of the treatments. The last four are the non-symmetrical ones.
[[1]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 4 1 6 3 8 5 7
[3,] 3 1 5 2 7 4 8 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 8 4 7 2 5 1 3
[7,] 7 5 8 3 6 1 4 2
[8,] 8 7 6 5 4 3 2 1

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

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

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

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

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

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

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

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

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

[[11]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 8 5 1 7 3 6 4
[3,] 3 5 4 6 1 8 2 7
[4,] 4 1 6 8 3 7 5 2
[5,] 5 7 1 3 2 4 8 6
[6,] 6 3 8 7 4 2 1 5
[7,] 7 6 2 5 8 1 4 3
[8,] 8 4 7 2 6 5 3 1

[[12]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 8 5 7 4 1 3 6
[3,] 3 5 2 6 8 7 1 4
[4,] 4 7 6 2 1 5 8 3
[5,] 5 4 8 1 6 3 2 7
[6,] 6 1 7 5 3 8 4 2
[7,] 7 3 1 8 2 4 6 5
[8,] 8 6 4 3 7 2 5 1

Calculation speed

Up to 5 treatments 100 repeats are done in the blink of an eye. For 6 it is doable (0.025 sec per run). For 7 treatments I did 10 repeats (1.6 sec per run), while for 8 one run took almost 800 sec. Clearly 9 treatments is too much.
> microbenchmark(gendesign(3),gendesign(4),gendesign(5),gendesign(6))
Note: no visible global function definition for 'error'
Unit: microseconds
expr min lq median uq max neval
gendesign(3) 107.029 113.2605 119.492 125.3565 240.449 100
gendesign(4) 343.081 355.5430 369.472 405.7585 544.677 100
gendesign(5) 1141.403 1163.0290 1202.615 1258.6955 1472.755 100
gendesign(6) 24753.851 25352.7755 25595.058 26917.8980 35721.433 100
> system.time(gendesign(7))
user system elapsed
1.64 0.00 1.67
> microbenchmark(gendesign(7),times=10L)
Unit: seconds
expr min lq median uq max neval
gendesign(7) 1.624484 1.647051 1.657468 1.669608 1.768613 10
> system.time(gendesign(8))
user system elapsed
793.72 0.17 805.82

Code

Code is below. Running the JIT saves half of the time. It is basically a simplified version of my old algorithm, but with a lot of small modifications, basically avoiding loops, ifs and a load of trial and error for speed.

library(microbenchmark)

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

carover <- matrix(TRUE,nrow=nr,ncol=nc)

for (i in 1L:(nc-1L)) carover[i,i+1] <- FALSE

todo <- nextpos(desmat)

desobject <- list(desmat=desmat,carover = carover,nc=1L:n,n=n,

index =1L,npos=nrow(todo),

row=todo[,1L],

col=todo[,2L])

desresult <- list()

addpoint(desobject,desresult)

}

modify <- function(desobject,row,col,i,previous) {

desobject$desmat[row,col] <- i

desobject$carover[previous,i] <- FALSE

desobject$index <- desobject$index + 1L

desobject}

addpoint <- function(desobject,desresult) {

if (desobject$index>desobject$npos) {

l <- length(desresult)

desresult[[l+1]] <- desobject$desmat

# cat('#')

return(desresult)

}

row <- desobject$row[desobject$index]

col <- desobject$col[desobject$index]

previous <- desobject$desmat[row,col-1L]

avoid <- c(desobject$desmat[row,],

desobject$desmat[,col])

nc <- desobject$nc[!is.element(desobject$nc,avoid) ]

nc <- nc[desobject$carover[previous,nc]]

for (i in nc) {

desresult <- addpoint(modify(desobject,row,col,i,previous)

,desresult)

}

desresult

}

library(compiler)

enableJIT(3)

Wiekvoet

Saturday, October 19, 2013

Carry-over balanced designs for 8 treatments

Designs

Calculation speed

Code

1 comment: