1ICFTs <- function (design, digits = 3, resk.only = TRUE,
2    kmin = NULL, kmax = ncol(design), detail = FALSE, with.blocks = FALSE,
3    conc = TRUE)
4{
5    if ("design" %in% class(design)) {
6        fn <- names(factor.names(design))
7        if (with.blocks)
8            fn <- c(fn, design.info(design)$block.name)
9        design <- design[, fn]
10        nfac <- length(fn)
11    }
12    else {
13        nfac <- ncol(design)
14        fn <- 1:nfac
15    }
16
17    nlev <- levels.no(design)
18    dfs <- nlev - 1
19
20    if (!is.data.frame(design)) design <- as.data.frame(design)
21    for (i in 1:nfac){
22       design[[i]] <- factor(design[[i]])
23       contrasts(design[[i]]) <- contr.XuWuPoly(nlev[i])
24    }
25
26    ks <- which(round(GWLP(design, kmax = kmax)[-1], 8) > 0)
27    N <- nrow(design)
28    if (length(ks) == 0) {
29        hilf <- list(list(ICFT = cbind(IC = 0, frequency = sum(nlev) -
30            kmax), IC1 = 0))
31        names(hilf) <- kmax
32        return(hilf)
33    }
34    k <- min(ks)
35    if (k < 2)
36        stop("resolution of design must be at least 2")
37    kminset <- FALSE
38    if (is.null(kmin))
39        kmin <- k
40    else {
41        if (!kmin == min(ks)) kminset <- TRUE
42        redu <- ks[ks >= kmin]
43        message(paste("Check sets of sizes ", paste(redu, collapse = ",")))
44        if (length(redu) == 0)
45            return()
46        kmin <- min(redu)
47        k <- kmin
48        ks <- redu
49    }
50    if (k >= kmin) {
51#        if (!"design" %in% class(design)) {
52#            if (!is.data.frame(design))
53#                design <- as.data.frame(design)
54#            faktoren <- sapply(design, is.factor)
55#            keinfaktor <- which(!faktoren)
56#            if (length(keinfaktor) > 0)
57#                for (i in keinfaktor) design[[i]] <- as.factor(design[[i]])
58#        }
59        k <- kmin
60        ns <- choose(nfac, k)
61        auswahl <- 1:ns
62        selproj <- sel <- nchoosek(nfac, k)
63        GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[,
64            obj])[-1]), 4)
65        selproj <- apply(selproj, 2, function(obj) paste(obj,
66            collapse = ":"))
67        names(auswahl) <- selproj
68 #       ergproj <- rep(NA, length(selproj))
69 #       dimnames(erg3) <- list(factor = fn, others = apply(sel,
70 #           2, function(obj) paste(obj, collapse = ":")), 1:(max(nlev) -
71 #           1))
72        if (resk.only) {
73                reskproj <- apply(GWLPs, 2, function(obj) all(obj[-k] ==
74                  0))
75                if (all(!reskproj)) {
76                  message("no projections with resolution ",
77                    k, " or higher")
78                  return()
79          }}
80
81            berechn <- lapply(auswahl, function(obj) {
82                hilf2 <- design[, sel[, obj]]
83                mmX <- model.matrix(formula(substitute(~.^km1,
84                    list(km1 = k))), data = hilf2)
85                mmX <- mmX[,-(1:(ncol(mmX) - prod(dfs[sel[, obj]])))]
86                hilfc <- svd(mmX)
87                hilf2 <- table(round(hilfc$d^2,6))
88                cumcounts <- cumsum(rev(hilf2))
89                from <- c(1, cumcounts[-length(hilf2)]+1)
90                hilf2 <- rep(0, prod(dfs[sel[, obj]])) ## initialize output vector
91                if (conc)
92                hilf2[from] <- sapply(1:length(from),
93                   function(obj) sum((hilfc$d^2*colMeans(hilfc$u)^2)[from[obj]:cumcounts[obj]])
94                   )
95                else {
96                   for (i in 1:length(from)){
97                         bereich <- from[i]:cumcounts[i]
98                         hilf2[bereich] <- mean((hilfc$d^2*colMeans(hilfc$u)^2)[bereich])
99                }}
100                list(hilf2, hilfc$d^2, colMeans(hilfc$u)^2)
101            })
102
103        ICs <- lapply(berechn, function(obj) obj[[1]])
104        sv2s <- lapply(berechn, function(obj) obj[[2]])
105        mean.u2s <- lapply(berechn, function(obj) obj[[3]])
106        rund <- lapply(ICs, function(obj) round(obj,digits))
107
108        ICFT <- table(unlist(rund))
109        ICFT <- cbind(IC = as.numeric(names(ICFT)), frequency = ICFT)
110        rownames(ICFT) <- rep("", nrow(ICFT))
111        aus <- list(ICFT = ICFT)
112
113        if (detail)
114            aus <- c(aus, list(ICs = rund, sv2s = sv2s, mean.u2s = mean.u2s ))
115
116        if (!resk.only || kminset){
117            aus <- list(aus); names(aus) <- k
118            if (!resk.only){
119            ks <- kmin:kmax
120            if (length(ks)>1){
121                 for (k in (kmin+1):kmax){
122                    ns <- choose(nfac, k)
123                    auswahl <- 1:ns
124                    selproj <- sel <- nchoosek(nfac, k)
125                    GWLPs <- round(apply(selproj, 2, function(obj) GWLP(design[,
126                       obj])[-1]), 4)
127                    selproj <- apply(selproj, 2, function(obj) paste(obj,
128                       collapse = ":"))
129                    names(auswahl) <- selproj
130
131            berechn <- lapply(auswahl, function(obj) {
132                hilf2 <- design[, sel[, obj]]
133                mmX <- model.matrix(formula(substitute(~.^km1,
134                    list(km1 = k))), data = hilf2)
135                mmX <- mmX[,-(1:(ncol(mmX) - prod(dfs[sel[, obj]])))]
136                hilfc <- svd(mmX)
137                hilf2 <- table(round(hilfc$d^2,6))
138                cumcounts <- cumsum(rev(hilf2))
139                from <- c(1, cumcounts[-length(hilf2)]+1)
140                hilf2 <- rep(0, prod(dfs[sel[, obj]])) ## initialize output vector
141                if (conc)
142                hilf2[from] <- sapply(1:length(from),
143                   function(obj) sum((hilfc$d^2*colMeans(hilfc$u)^2)[from[obj]:cumcounts[obj]])
144                   )
145                else {
146                   for (i in 1:length(from)){
147                         bereich <- from[i]:cumcounts[i]
148                         hilf2[bereich] <- mean((hilfc$d^2*colMeans(hilfc$u)^2)[bereich])
149                }}
150                list(hilf2, hilfc$d^2, colMeans(hilfc$u)^2)
151            })
152
153           ICs <- lapply(berechn, function(obj) obj[[1]])
154           sv2s <- lapply(berechn, function(obj) obj[[2]])
155           mean.u2s <- lapply(berechn, function(obj) obj[[3]])
156           rund <- lapply(ICs, function(obj) round(obj,digits))
157
158           ICFT <- table(unlist(rund))
159           ICFT <- cbind(IC = as.numeric(names(ICFT)), frequency = ICFT)
160           rownames(ICFT) <- rep("", nrow(ICFT))
161           ausn <- list(ICFT = ICFT)
162           if (detail)
163              ausn <- c(ausn, list(ICs = rund, sv2s = sv2s, mean.u2s = mean.u2s ))
164           ausn <- list(ausn)
165           names(ausn) <- k
166           aus <- c(aus, ausn)
167                 }
168             }
169           }}
170    }
171    else aus <- list(ICFT = NULL)
172    aus
173}
174