1oa_feasible <- function(nruns, nlevels, strength=2, verbose=TRUE){ 2 aus <- TRUE 3 nlev <- nlevels 4 if (!all(is.numeric(c(nruns,nlev, strength)))) stop("nruns, nlev and strength must be numeric") 5 if (!all(c(nruns,nlev)>0)) stop("nruns, nlev and strength must be positive") 6 if (!all(c(nruns,nlev, strength)%%1==0)) stop("nruns, nlev and strength must be integer") 7 if (!(length(nruns)==1 && length(strength==1))) stop("nruns and strength must be scalar") 8 nfac <- length(nlev) 9 if (!nfac >=2) stop("too few factors") 10 if (nfac < strength) stop("strength larger than number of factors") 11 12 ## tabulation of nlevels would be useful for large nfac 13 ## since this is not the situation for this package, not pursued 14 ## because it sincerely complicates things 15 tablev <- table(nlev) 16 s <- as.numeric(names(tablev)) ## the different levels 17 fs <- tablev ## their frequencies 18 nlevelmix <- length(s) ## their number 19 u <- strength%/%2 ## for calculation of Rao bound 20 21 ## handle simpler pure level situation first 22 if (nlevelmix == 1){ 23 ## pure level designs 24 if (!nruns%% s^strength==0) { 25 if (verbose) cat("nruns is not divisible by", s^strength, fill=TRUE) 26 return(FALSE) 27 } 28 ## Bierbrauer et al.s bound for 2-level (Thms 7.1 and 7.3 with Cor 7.4 (and abstract)) 29 if (s==2){ 30 bound <- round(s^nfac - min(nfac*2^(nfac-1)/(strength+1), (nfac+1)*2^(nfac-2)/(u+1)), 2) 31 if (nruns < bound){ 32 if (verbose) cat("nruns is smaller than Bierbrauer et al.s bound ", bound, " for 2-level designs", fill=TRUE) 33 return(FALSE) 34 } 35 } 36 ## Bierbrauer's bound 37 bound <- round(s^nfac*(1-(s-1)*nfac/(s*(strength+1))), 2) 38 if (nruns < bound) { 39 if (verbose) cat("nruns is smaller than Bierbrauer's bound", bound, " for pure level designs", fill=TRUE) 40 return(FALSE) 41 } 42 ## Rao's bound 43 ii <- 0:(strength%/%2) 44 bound <- sum(choose(nfac, ii)*(s-1)^ii) 45 ## odd strengths 46 if (strength %% 2 == 1) bound <- bound + choose(nfac-1,u)*(s-1)^(u + 1) 47 if (nruns < bound){ 48 if (verbose) cat("nruns is smaller than Rao's bound", bound, "for pure level designs", fill=TRUE) 49 return(FALSE) 50 } 51 ## Bush bound 52 if (nruns==s^strength){ 53 ## HSS Theorem 2.19 54 if (s <= strength) bound <- strength + 1 55 else{ 56 ## s > strength 57 bound <- s + strength - 1 58 if (3 <= strength && s%%2 == 1) bound <- s + strength - 2 59 } 60 if (nfac > bound){ 61 if (verbose) cat("more than ", bound, " factors,", " Bush bound violated", fill=TRUE) 62 return(FALSE) 63 } 64 } 65 ## Bose Bush bounds for strengths 2 and 3, not for 2-level 66 if (nruns > s^strength && !s==2){ 67 ## HSS theorems 2.8 and 2.11 68 ## for s==2, lambda-1 is always a multiple of s-1 69 lambda <- nruns%/%(s^strength) 70 b <- (lambda - 1)%%(s-1) 71 if (b > 0 && strength %in% c(2,3)){ 72 a <- (lambda-1)%/%(s-1) 73 theta <- (sqrt(1+4*s*(s-1-b)) - (2*s-2*b-1))/2 74 if (strength==2) bound <- lambda*(s + 1) + a - floor(theta) - 1 75 if (strength==3) bound <- lambda*(s + 1) + a - floor(theta) 76 if (nfac > bound){ 77 if (verbose) cat("more than ", bound, " factors,", " Bose/Bush bound violated", fill=TRUE) 78 return(FALSE) 79 } 80 } 81 } 82 } 83 else{ 84 ## mixed level designs 85 ## Bierbrauer's mixed level bound, Theorem 2.2 of Diestelkamp 2004 86 ## this bound is stronger than Raos for very large values of strength 87 ## (for which computations of Rao's bound in the current implementation 88 ## will be prohibitive) 89 ## and it is proven only for strength > (SM-1)*nfac/SM - 1 90 sT <- sum(nlev) 91 sM <- max(nlev) 92 sm <- min(nlev) 93 if(strength > (sM-1)*nfac/sM - 1){ 94 bound <- floor(sm^nfac*(1 - max(0, (sT - nfac)/(sT + (strength - nfac + 1)*sM)))) 95 if (nruns < bound){ 96 if (verbose) cat("nruns is smaller than Diestelkamp's mixed level version", 97 " of Bierbrauer's bound", bound, fill=TRUE) 98 return(FALSE) 99 } 100 } 101 102 ## Rao's mixed level bound, Theorem 3.2 of Diestelkamp 2004 103 bound <- sum(nlev - 1) + 1 ## df for up to main effects 104 if (strength >= 2 && nruns < bound){ 105 if (verbose) cat("nruns is smaller than the ", bound, " df needed for main effects", fill=TRUE) 106 return(FALSE) 107 } 108 if (strength == 3){ 109 ## add worst-case df for 2fi 110 bound <- bound + max(sapply(1:nfac, function(obj) (nlev[obj] - 1)*sum(nlev[-obj] - 1))) 111 if (nruns < bound){ 112 if (verbose) cat("nruns is smaller than Rao's bound", bound, "for strength 3", fill=TRUE) 113 return(FALSE) 114 } 115 } 116 if (strength>=4){ 117 for (ii in 2:u){ 118 ## df of strength ii effects 119 sets <- nchoosek(nfac, ii) 120 prods <- sapply(1:ncol(sets), function(obj){ 121 ## obj is a column number of sets 122 ## product of dfs in set 123 prod(nlev[sets[,obj]] - 1) 124 }) 125 dazu <- sum(prods) 126 127 ## increase dazu for the last contribution 128 ## in case of odd strength 129 ## by different summand 130 if (strength == 2 * ii + 1){ 131 hilf <- sapply(1:ncol(sets), function(obj){ 132 (max(nlev[setdiff(1:nfac, sets[,obj])]) - 1)*prods[obj] 133 }) 134 dazu <- dazu + max(hilf) 135 } 136 bound <- bound + dazu 137 } 138 if (nruns < bound){ 139 if (verbose) cat("nruns is smaller than Rao's bound ", bound, " for strength ", strength, fill=TRUE) 140 return(FALSE) 141 } 142 } 143 ## end of Rao's mixed level bound 144 145 ## condition on LCM 146 clcm <- numbers::mLCM(apply(matrix(nlev[nchoosek(nfac,strength)], nrow=strength), 147 2, prod)) 148 if (!nruns %% clcm == 0){ 149 if (verbose) cat("nruns is not divisible by ", clcm, fill=TRUE) 150 return(FALSE) 151 } 152 } 153 ## no violation of criteria for oa with required strength was found 154 if (verbose) cat("no violation of necessary criteria", " for strength ", strength, " was found", fill=TRUE) 155 TRUE 156} 157