1.listMat <- function(M, f) { 2 l <- vector("list", ncol(M)) 3 names(l) <- dimnames(M)[[2]] 4 for(i in 1:ncol(M)) 5 l[[i]] <- f(M[, i]) 6 l 7} 8 9.paste0 <- function(...) paste(..., sep = "") 10 11.zf <- function(x) { 12 x <- as.character(x) 13 m <- max(n <- nchar(x)) 14 z <- paste(rep(0, m), collapse="") 15 .paste0(substring(z, 0, m-n), x) 16} 17 18conf.design <- function(G, p, block.name = "Blocks", 19 treatment.names = NULL) { 20 if(!is.matrix(G)) 21 G <- rbind(G) 22 if(is.null(treatment.names)) 23 treatment.names <- if(is.null(nam <- dimnames(G)[[2]])) 24 .paste0("T", .zf(1:ncol(G))) else nam 25 26 stopifnot(is.character(treatment.names), length(treatment.names) == ncol(G), 27 is.numeric(G), all(G >= 0), all(G %% 1 == 0), 28 is.numeric(p), length(p) == 1, p > 0, p %% 1 ==0, 29 p %in% primes(p), 30 is.character(block.name), length(block.name) == 1) 31 32 D <- as.matrix(expand.grid(rep(list(0:(p-1)), ncol(G)))) 33 B <- .listMat((D %*% t(G)) %% p, format) 34 o <- do.call(order, B) 35 B <- do.call(".paste0", B) 36 D <- cbind(B, format(D))[o, ] 37 dimnames(D) <- list(B[o], c(block.name, treatment.names)) 38 data.frame(.listMat(D, factor)) 39} 40 41.space <- function(G, p) { 42### 43### Generate all distinct linear combinations of the rows of G over GF(p) 44### 45 x <- 0:(p - 1) 46 M <- as.matrix(x) 47 k <- nrow(G) 48 if(k > 1) { 49 for(i in 2:k) { 50 N <- NULL 51 for(j in x) 52 N <- rbind(N, cbind(M, j)) 53 M <- N 54 } 55 } 56 M <- (M %*% G) %% p 57### 58### if the rows of G can be assumed linearly independent the rest can 59### be omitted. 60### 61 m <- 0 62 for(j in 1:ncol(M)) 63 m <- p * m + M[, j] 64 M[!duplicated(m), , drop = FALSE] 65} 66 67conf.set <- function(G, p) { 68 S <- .space(G, p) 69 dimnames(S)[[2]] <- dimnames(G)[[2]] 70 S[apply(S, 1, function(x) 71 any(t <- x > 0) && x[t][1] == 1), ] 72} 73 74direct.sum <- function(D1, ..., tiebreak = letters) { 75 l <- list(...) 76 if(length(l) == 0) return(D1) 77 D2 <- if(length(l) == 1) l[[1]] else { 78 Recall(..., tiebreak = c(tiebreak[-1], tiebreak[1])) 79 } 80 stopifnot(is.data.frame(D1), is.data.frame(D2)) 81 n1 <- nrow(D1) 82 n2 <- nrow(D2) 83 D1 <- D1[rep(seq_len(n1), each = n2), ] 84 D2 <- D2[rep(seq_len(n2), times = n1), ] 85 if(any(dups <- names(D2) %in% names(D1))) 86 names(D2)[dups] <- .paste0(names(D2)[dups], tiebreak[1]) 87 D <- cbind(D1, D2) 88 row.names(D) <- format(seq_len(nrow(D))) 89 D 90} 91 92factorize <- function(x, ...) 93 UseMethod("factorize") 94 95factorize.default <- function(x, divisors = primes(max(x)), ...) { 96 stopifnot(is.numeric(x)) 97 if (length(x) > 1) { 98 l <- vector("list", length(x)) 99 names(l) <- as.character(x) 100 for (i in seq_along(x)) 101 l[[i]] <- Recall(x[i], divisors = divisors, ...) 102 return(l) 103 } 104 if (x %% 1 > 0 || x < 2) 105 return(x) 106 tab <- divisors 107 fac <- numeric(0) 108 while(length(tab <- tab[x %% tab == 0]) > 0) { 109 x <- x/prod(tab) 110 fac <- c(fac, tab) 111 } 112 sort(fac) 113} 114 115factorize.factor <- function(x, name = deparse(substitute(x)), 116 extension = letters, ...) { 117 llev <- factorize.default(length(levels(x))) 118 if(length(llev) == 1) 119 return(structure(data.frame(x), names = name)) 120 121 D <- expand.grid(lapply(llev, seq_len))[x, ] - 1 122 123 D <- lapply(D, factor) 124 names(D) <- .paste0(name, rep(extension, length.out = length(D))) 125 D <- data.frame(D) 126 row.names(D) <- format(seq_len(nrow(D))) 127 D 128} 129 130.makeList <- function(x) if(is.list(x)) x else list(x) 131 132join <- function(...) { 133 l <- list(...) 134 l <- do.call(c, lapply(l, .makeList)) 135 l <- lapply(l, function(f) 136 format(as.character(f))) 137 factor(do.call(paste, c(l, list(sep = ":")))) 138} 139 140rjoin <- function(..., part.name = "Part") { 141 l <- lapply(list(...), as.data.frame) ### for some safety... 142 checkNames <- all(sapply(l, function(bit) 143 all(names(bit) == names(l[[1]])))) 144 if(!checkNames) stop("The column names differ between components!") 145 bf <- factor(.paste0(part.name, .zf(rep(1:length(l), sapply(l, nrow))))) 146 D <- data.frame(Part = bf, do.call("rbind", l)) 147 names(D)[1] <- part.name 148 D 149} 150 151primes <- function(n) { 152### Find all primes less than n (or max(n) if length(n) > 1). 153### Uses an obvious sieve method. Nothing flash. 154### 155### 2013: This now uses a slightly improved coding of the version of 156### the algorithm used in the pracma package primes() function 157### 158 stopifnot(is.numeric(n), all(n %% 1 == 0)) 159 if ((M2 <- max(n)) <= 1) 160 return(numeric(0)) 161 x <- seq(1, M2, 2) 162 np <- length(x) 163 x[1] <- 2 164 if(M2 > 8) { 165 top <- floor(sqrt(M2)) 166 p <- 1 167 while((p <- p+2) <= top) 168 if(x[(p + 1)/2] > 0) 169 x[seq((p*p + 1)/2, np, p)] <- 0 170 } 171 x[x > 0] 172} 173 174 175 176