# Procedures from the new version of the package OptimalDesign (under preparation) # These procedures are enough to run the numerical computations for the AQuA project: # https://arxiv.org/abs/1801.09124 ###################################################################### Fx.cube <- function (formula, lower=NULL, upper=NULL, n.levels=NULL) { # Purpose: # Creates the n times m matrix of all permissible regressors # of a linear regression model on a line, square or cube (up to 9 factors) # # Arguments: # formula ... formula of the model # The rules for creating the formula are standard for R but: # 1) the formula must not contain the dependent variable (it is one-sided) # 2) the d factors (variables) must be labeled x1,x2,x3,... # lower ... d-dimensional vector of the smallest values of factors # if lower==NULL, the program sets lower <- rep(-1, d) # upper ... d-dimensional vector of the largest values of factors # if upper==NULL, the program sets upper <- rep(1, d) # n.levels ... d-dimensional vector of the numbers of levels of each factor # if n.levels==NULL, the program sets n.levels <- rep(2, d) # # Example: # F.cube(~x1+x2+I(x1*x2), n.levels=c(5, 5)) # # Authors: # Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Do some basic verification of the input library(plyr) if(!is.formula(formula)) stop("The first argument does not appear to be a valid R formula.") av <- all.vars(formula); d <- length(av) pot.vars <- c("x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9") if(length(intersect(av, pot.vars[1:d]))!=d) stop("Variables in the formula must be x1,x2,...,xd, d<=9, without skipping any intermediate variable.") if(is.null(lower)) lower <- rep(-1, d) if (!is.numeric(lower) || !is.finite(lower) || !is.vector(lower) || length(lower)!=d) stop("lower must be NULL, a real number, or a real vector of length d, where d is the number of variables x1,x2,...,xd in the formula.") if(is.null(upper)) upper <- rep(1, d) if (!is.numeric(upper) || !is.finite(upper) || !is.vector(upper) || length(upper)!=d) stop("upper must be NULL, a real number, or a real vector of length d, where d is the number of variables x1,x2,...,xd in the formula.") if(is.null(n.levels)) n.levels <- rep(2, d) if (!is.numeric(n.levels) || !is.finite(n.levels) || !is.vector(n.levels) || !all(n.levels==round(n.levels)) || length(n.levels)!=d || min(n.levels)<2) stop("n.levels must be NULL, a natural number >=2, or a vector of length d consisting of natural numbers >=2, where d is the number of variables x1,x2,...,xd in the formula.") if (any(lower >= upper)) stop ("lower[i] must be smaller than upper[i] for all i!") lst <- labs <- c() for (i in 1:d) { lst <- c(lst, list(seq(lower[i], upper[i], length=n.levels[i]))) labs <- c(labs,paste("x", as.character(i), sep="")) } data <- expand.grid(lst) names(data) <- labs mf <- model.frame(formula, data) Fx <- as.data.frame(model.matrix(attr(mf, "terms"), mf)) return(as.matrix(Fx)) } ###################################### Fx.ItoA <- function (Fx, eta = NULL) { # Purpose: # Conversion of non-singular I-optimality to A-optimality # # Arguments: # Fx ... n times m model matrix, where n>=2 is the size of the design space # and m>=2 (m <= n) is the number of model parameters # eta ... a non-negative vector of length n containing the required weights # of the variances of BLUEs of f^T(x)\theta # if eta==NULL the program selects the uniform weights eta <- rep(1, n) # # Output: # The nxm matrix of a transformed model, such that an A-optimal design for # the transformed model is an I-optimal design for the original model # # Authors: # Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Input check if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be greater than or equal to the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be at least 2. The case of a single parameter is trivial.") if (is.null(eta)) eta <- rep(1, n) if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative numeric vector of length nrow(Fx).") L <- infmat(Fx, eta) if(rcond(L) < 1e-9) warning("The problem of I-optimality is badly conditioned and the results may be unreliable.") return(Fx %*% solve(chol(L))) } ########################################### Fx.simp <- function(formula, n.levels=NULL) { # Purpose: # Creates the n times m matrix of all permissible regressors # of a linear regression model on a regular simplex grid (up to 9 factors) # # Arguments: # formula ... formula of the model # The rules for creating the formula are standard for R but: # 1) the formula must not contain the dependent variable (it is one-sided) # 2) the d factors (independent variables) must be labeled x1,x2,x3,... # n.levels ... the number of levels (>=2) of each factor # if n.levels==NULL, the program sets n.levels <- 2*d+1 # # Value: # The n times m matrix where n is the number of combinations of factors levels # and m is the number of parameters of the LRM # # Example: # Fx.simp(~x1+x2+x3+I(x1*x2)+I(x1*x3)+I(x2*x3)-1, 5) # # Authors: Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Do some basic verification of the input library(plyr) if(!is.formula(formula)) stop("The first argument does not appear to be a valid R formula.") av <- all.vars(formula); d <- length(av) if(d < 2) stop("The number of variables must be at least 2.") pot.vars <- c("x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9") if(length(intersect(av, pot.vars[1:d]))!=d) stop("Variables in the formula must be x1,x2,...,xd, d<=9, without skipping any intermediate variable.") if(is.null(n.levels)) n.levels <- 2*d+1 if(!is.numeric(n.levels) || !is.finite(n.levels) || length(n.levels)!=1 || n.levels <= 1) stop("n.levels must be an integer greater than 1.") lst <- labs <- c() cmb <- t(combn(n.levels + d - 2, d - 1)) data <- as.data.frame((cbind(cmb, n.levels + d - 1) - cbind(0, cmb) - 1) / (n.levels - 1)) for (i in 1:d) labs <- c(labs, paste("x", as.character(i), sep="")) names(data) <- labs[1:d] mf <- model.frame(formula, data) Fx <- as.data.frame(model.matrix(attr(mf, "terms"), mf)) return(as.matrix(Fx)) } #################################################### infmat <- function(Fx, w, sparse=TRUE){ # Purpose: # Computes the information matrix of the design w in the model defined by Fx # # Arguments: # Fx ... n times m matrix of real numbers with rows corresponding to model regressors # and columns corresponding to parameters # w ... vector of non-negative real numbers of length n representing the design # Note: The design w need not be normalized (sum(w)==1) # sparse ... boolean variable; TRUE is suggested if the number of non-zero # elements of w tends to be much smaller that the length of w # # Value: # Information matrix of the design w in the model with all # possible regressors given by the rows of Fx # # Authors: Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Check of the input if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") if (!is.vector(w) || !is.numeric(w) || !all(is.finite(w)) || (min(w) < 0)) stop("w must be a vector of non-negative real numbers.") if (length(w) != nrow(Fx)) stop("The length of w must be equal to the number nrow(Fx) of design points.") if (!is.logical(sparse) || !is.vector(sparse) || length(sparse)>1) stop("sparse must be a single logical value TRUE or FALSE.") if(sparse) { supp <- w>0 if(sum(supp) > 1) { return(crossprod(sqrt(w[supp]) * Fx[supp,])) } else { return(Fx[supp,] %*% t(Fx[supp,])) } } else { return(crossprod(sqrt(w) * Fx)) } } ########################################################### od.A.AQuA <- function(Fx, b1, A1, b2, A2, b3, A3, bin, M.anchor, ver, conic, t.max) { # Note: this function should be only run via the wrapper od.AQuA. if(!requireNamespace('gurobi', quietly=TRUE)) stop("AQuA requires the gurobi package.") n <- nrow(Fx); m <- ncol(Fx) k1 <- length(b1); k2 <- length(b2); k3 <- length(b3) b <- c(b1, b2, b3); A <- rbind(A1, A2, A3) sense <- c(rep("<", k1), rep(">", k2), rep("=", k3)) if (is.null(M.anchor)) { if (length(b)==1 && all(A==matrix(1, ncol=n, nrow=1))) { M.anchor <- infmat(Fx, b * od.A.REX(Fx, t.max=t.max/5)$w.best) } else { w.app <- od.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin, t.max=t.max/3, crit="A")$w.best if (is.null(w.app)) { stop("Failure in the computation of M.anchor (and no M.anchor supplied by the user).") } else { M.anchor <- infmat(Fx, w.app) } } } M1 <- solve(M.anchor) M2 <- M1 %*% M1 if(conic) { if(!requireNamespace('Matrix', quietly=TRUE)) stop("Conic version of AQuA requires package Matrix.") if(!requireNamespace('matrixcalc', quietly=TRUE)) stop("Conic version of AQuA requires package matrixcalc.") s <- m*(m+1)/2; Gm <- matrixcalc::duplication.matrix(m); eps <- 1e-12 vec.M2 <- matrixcalc::vec(M2) tilde.h <- t(Gm) %*% vec.M2 if(ver=="+") tilde.Q <- t(Gm) %*% (-vec.M2 %*% t(vec.M2) / sum(diag(M1)) + 0.5 * M1 %x% M2 + 0.5 * M2 %x% M1) %*% Gm if(ver=="-") tilde.Q <- (1/6) * t(Gm) %*% (M1 %x% M2 + M2 %x% M1) %*% Gm if(rcond(tilde.Q) > eps) { tilde.C <- t(chol(tilde.Q)) } else { eig <- eigen(tilde.Q, symmetric=TRUE) mx <- max(eig$values); ts <- sum(eig$values/mx > eps) tilde.C <- eig$vectors[,1:ts] %*% diag(sqrt(eig$values[1:ts])) } diff <- max(abs(tilde.Q - tilde.C%*%t(tilde.C))) print(paste("Decomposition instability:", diff)) H <- matrix(0, nrow=n, ncol=s); k <- 0 for(i1 in 1:m) { for(i2 in i1:m) { k <- k + 1 H[,k] <- Fx[,i1] * Fx[,i2] } } h <- H %*% tilde.h; S <- H %*% tilde.C; u <- ncol(S) model <- list(); k <- nrow(A); np <- n + u + 3 Aeq <- matrix(0, nrow = k + 2 + u, ncol = np); Aeq[1:k, 1:n] <- A Aeq[k + 1, n + 1] <- 1; Aeq[k + 1, n + 3] <- -1/sqrt(2) Aeq[k + 2, n + 2] <- 1; Aeq[k + 2, n + 3] <- 1/sqrt(2) Aeq[(k + 3):(k + 2 + u), 1:n] <- t(S) Aeq[(k + 3):(k + 2 + u), (n + 4):np] <- -diag(u) model$A <- Matrix::Matrix(Aeq, sparse=TRUE) beq <- rep(0, k + 2 + u); beq[1:k] <- b beq[k + 1] <- 1/2/sqrt(2); beq[k + 2] <- 1/2/sqrt(2) model$rhs <- beq model$sense <- c(rep("<", k1), rep(">", k2), rep("=", k3 + u + 2)) model$cones <- list(as.list(c(n+1, n+2, (n+4):np))) lb <- rep(-Inf, np); lb[1:n] <- 0; model$lb <- lb ob <- rep(0, np); ob[1:n] <- -h; ob[n + 3] <- 1 model$obj <- ob; vtypes <- rep("C", np) if(bin) {vtypes[1:n] <- "B"} else {vtypes[1:n] <- "I"} model$vtypes <- vtypes; model$modelsense <- "min" params <- list(TimeLimit=t.max) result <- suppressWarnings(gurobi::gurobi(model, params)); gc() } else { f <- rep(0, n); H <- matrix(0, nrow=n, ncol=n) for (i in 1:n) f[i] <- -t(Fx[i, ]) %*% M2 %*% Fx[i, ] if(ver=="+") { for (i in 1:n) for (j in i:n){ H[j, i] <- H[i, j] <- -(t(Fx[i, ]) %*% M2 %*% Fx[i, ]) * (t(Fx[j, ]) %*% M2 %*% Fx[j, ]) / (sum(diag(M1))) + (t(Fx[i, ]) %*% M2 %*% Fx[j, ]) * (t(Fx[i, ]) %*% M1 %*% Fx[j, ]) } } if(ver=="-") { for (i in 1:n) for (j in i:n){ H[j, i] <- H[i, j] <- (1/3) * (t(Fx[i, ]) %*% M2 %*% Fx[j, ]) * (t(Fx[i, ]) %*% M1 %*% Fx[j, ]) } } model <- list(); model$obj <- f model$Q <- H; model$A <- Matrix::Matrix(A, sparse=TRUE); model$rhs <- c(b) model$sense <- sense; model$lb <- rep(0, n) if(bin) {model$vtypes <- "B"} else {model$vtypes <- "I"} params <- list(TimeLimit=t.max) result <- gurobi::gurobi(model, params); gc() } res <- result$x if (!is.vector(res) || !is.numeric(res) || !all(is.finite(res)) || length(res)1)) w.best <- NULL } return(list(w.best=w.best, status=result$status)) } ################################################################# od.A.KL <- function(Fx, N, Phi.app=NULL, w1=NULL, K=NULL, L=NULL, rest.max=Inf, t.max=60) { # Purpose: # An implementation of the KL exchange algorithm for computing # A-efficient exact designs under the standard "size" constraint # # Arguments: # Fx ... n times m model matrix, where n>=2 is the size of the design space # and m>=2 (m<=n) is the number of model parameters # N ... required size of the design (the number of experimental trials) # Phi.app ... optimal value of the corresponding approximate (relaxed) problem # If NULL, the value is pre-computed by REX # w1 ... initial design used for the first restart of the algorithm # It should be either NULL or a non-singular design of size N # K ... (upper bound on the) number of "least promising" support points # of the current design, for which exchanges are attempted # L ... (upper bound on the) number of "most promising" candidate # design points for which exchanges are attempted # If K is NULL or L is NULL, a reasonable value is assigned automatically # rest.max ... limit on the number of restarts # t.max ... approximate upper limit on the time of computation # If the algorithm stops in a local optimum before the time ellapsed, # the computation is restarted # # Output is a list of: # w.best ... the best design found within the time limit # Phi.best ... the value of the D-criterion of w.best # eff.best ... the efficiency of w.best wrt the approximate optimum # n.ex ... the number of exchanges # n.rest ... the number of restarts executed # t.act ... the actual time taken by the computation # # Notes: # Works reasonably well for most problems with n<=10000 and m<=10 # # Authors: # Radoslav Harman and Lenka Filova eps <- 1e-09; del <- 0.1 n <- nrow(Fx); m <- ncol(Fx); one.m <- rep(1, m) next.sec <- 0; n.ex <- 0; n.rest <- 0 if (is.null(Phi.app)) Phi.app <- N * od.A.REX(Fx, od.PIN(Fx, m, crit="A", rnd=0)$w)$Phi.best start <- as.numeric(proc.time()[3]) info <- paste("Running A.KL for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep="") print(info, quote=FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, " and m=", m, ".", sep="") print(info, quote=FALSE) if (is.null(K)) K <- max(10, min(ceiling(sqrt(c(N, n))))) if (is.null(L)) L <- max(10, ceiling(sqrt(n))) print(paste("Setting K=", K, ", L=", L, sep=""), quote=FALSE) A <- array(0, dim=c(n, m, m)) for (i in 1:n) A[i, , ] <- tcrossprod(Fx[i, ]) finish.all <- FALSE; trMinv.best <- Inf while (!finish.all && n.rest m) { for (k in 1:(N - m)) { i.supp <- sample(1:n, 1) w[i.supp] <- w[i.supp] + 1 M <- M + A[i.supp, , ] } } } else { w <- w1; M <- crossprod(sqrt(w) * Fx) w1 <- NULL # Start the second run from a "pinned" design } trMinv <- sum(diag(solve(M))) if (trMinv.best > trMinv){ w.best <- w; M.best <- M; trMinv.best <- trMinv } finish.all <- finish <- as.numeric(proc.time()[3]) > start + t.max while (!finish){ tm <- as.numeric(proc.time()[3]) - start if (tm > next.sec) { Phi.best <- m / trMinv.best info <- paste("A.KL Time:", round(tm, 1), "Best value:", round(Phi.best, 6), "Efficiency:", round(Phi.best/Phi.app, 6)) print(info, quote=FALSE); next.sec <- ceiling(tm) } a.fun <- ((Fx %*% solve(M))^2) %*% one.m; supp <- (1:n)[w > 0.5] Kact <- min(length(supp), K); Lact <- min(n, L) Kind <- supp[order(a.fun[supp])][1:Kact] Lind <- order(a.fun)[(n-Lact+1):n] imp <- FALSE for (iL in Lact:1){ for (iK in 1:Kact){ M.temp <- M + A[Lind[iL], , ] - A[Kind[iK], , ] if(det(M.temp)>1e-8) { trMinv.temp <- sum(diag(solve(M.temp))) } else { trMinv.temp <- Inf } if(trMinv.temp < trMinv){ w[Lind[iL]] <- w[Lind[iL]] + 1 w[Kind[iK]] <- w[Kind[iK]] - 1 M <- M.temp; trMinv <- trMinv.temp n.ex <- n.ex + 1; imp <- TRUE break } } if (imp) break } if (as.numeric(proc.time()[3]) > start + t.max) finish.all <- TRUE if (finish.all || !imp) finish <- TRUE } M <- crossprod(sqrt(w) * Fx); trMinv <- sum(diag(solve(M))) if (trMinv < trMinv.best){ w.best <- w; M.best <- M; trMinv.best <- trMinv } } t.act <- round(as.numeric(proc.time()[3]) - start, 2) info <- paste("A.KL finished after", t.act, "seconds at", Sys.time()) print(info, quote=FALSE) info <- paste("with", n.rest, "restarts and", n.ex, "exchanges.") print(info, quote=FALSE) Phi.best <- optcrit(Fx, w.best, crit="A") return(list(w.best=w.best, Phi.best=Phi.best, eff.best=Phi.best/Phi.app, n.ex=n.ex, n.rest=n.rest, t.act=t.act)) } ######################################################## od.A.REX <- function (Fx, w1=NULL, ver=1, gamma=1, eff=1-1e-9, it.max=Inf, t.max=30) { # Computes an A-optimal approximate design of experiments on a finite design space # Can also be used for computing I-optimal designs because # the I-optimal design problem can be directly converted to a problem of A-optimality # Based on the paper: https://arxiv.org/abs/1801.05661 # Authors: Radoslav Harman, Lenka Filova, Peter Richtarik # # Arguments: # Fx ... n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters # w1 ... initial design; it must have a non-singular information matrix # and for the sake of efficiency, it should have a small support (e.g., m) # if NULL, a non-singular initial design is generated by PINs. # ver ... version of the algorithm (1 = with and 0 = without the nullity control) # gamma ... positive parameter regulating the size of the exchange batch # eff ... threshold on the minumim design efficiency to stop the computation # it.max ... maximum allowed number of iterations # t.max ... threshold for the maximum computation time # # Output is the list with components: # w.best ... the resulting approximate design # Phi.best ... the A-criterion value of w.best # eff.best ... a lower bound on the eff of w.best wrt an ideal A-optimal design # n.iter ... number of iterations performed # t.act ... the actual time of the computation # Example: An A-optimal design on 50000 random regressors in R^10 # n <- 50000; m <- 10; F.rnd <- matrix(rnorm(n*m), ncol=m) # res <- od.A.REX(F.rnd); supp <- res$w.best>0 # print(cbind((1:n)[supp], res$w.best[supp])) stepsize <- function (M, w, v) { # Computes the A-optimal stepsize alpha for the weight exchange: # w[v[1]] <- w[v[1]] - alpha; w[v[2]] <- w[v[2]] + alpha M.inv <- solve(M) dv <- Fx[v, ] %*% M.inv %*% t(Fx[v, ]) av <- Fx[v, ] %*% M.inv %*% M.inv %*% t(Fx[v, ]) A <- av[2,2] - av[1,1]; C <- dv[2,2] - dv[1,1]; D <- dv[1,1] * dv[2,2] - dv[1,2]^2 B <- 2 * dv[1,2] * av[1,2] - dv[2,2] * av[1,1] - dv[1,1] * av[2,2] G <- A * D + B * C; k <- v[1]; l <- v[2] if ((abs(G) < del.alpha) && (abs(B) > del.alpha)) { r <- -A / (2*B) if ((-w[l] <= r) && (r <= w[k])) return(r) } if (abs(G) > 0){ r <- -(B + sqrt(B^2 - A * G)) / G if ((abs(r) < del.alpha) && (B < 0)) { x <- A * G / B^2 r <- -B * (x^3/16 + x^2/8 + x/2) / G } if ((-w[l] <= r) && (r <= w[k])) return(r) } if (A > del.alpha) { return(w[k]) } else if (A < -del.alpha) { return(-w[l]) } else { return(0) } } start <- as.numeric(proc.time()[3]); del <- 1e-24; del.alpha <- 1e-14 Fx <- as.matrix(Fx); n <- nrow(Fx); m <- ncol(Fx) info <- paste("Running od.A.REX for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, " and m=", m, ".", sep="") print(info, quote=FALSE) eff.inv <- 1/eff; n.iter <- 0; L <- min(n, gamma*m) lx.vec <- rep(0, L); index <- 1:n; one <- rep(1, m) if(is.null(w1)) w1 <- od.PIN(Fx, m, crit="A")$w supp <- (1:n)[w1>0]; K <- length(supp); Fx.supp <- Fx[supp, ] w <- rep(0, n); w[supp] <- 1 / length(supp); w.supp <- w[supp] M <- crossprod(sqrt(w.supp) * Fx.supp) M.inv <- solve(M); a.fun <- ((Fx %*% M.inv)^2) %*% one ord <- order(a.fun, decreasing=TRUE) lx.vec <- sample(ord[1:L]); kx.vec <- sample(supp) while (TRUE) { n.iter <- n.iter + 1; ord1 <- which.min(a.fun[supp]) kb <- supp[ord1]; lb <- ord[1]; v <- c(kb, lb) alpha <- stepsize(M, w, v) w[kb] <- w[kb] - alpha; w[lb] <- w[lb] + alpha M <- M + alpha * (tcrossprod(Fx[lb, ]) - tcrossprod(Fx[kb, ])) if((w[kb] < del) && (ver==1)) { # LBE is nullifying and the version is 1 for(l in 1:L) { lx <- lx.vec[l]; Alx <- tcrossprod(Fx[lx, ]) for (k in 1:K) { kx <- kx.vec[k]; v <- c(kx, lx) alpha <- stepsize(M, w, v) wkx.temp <- w[kx] - alpha; wlx.temp <- w[lx] + alpha if((wkx.temp < del) || (wlx.temp < del)) { w[kx] <- wkx.temp; w[lx] <- wlx.temp M <- M + alpha * (Alx - tcrossprod(Fx[kx, ])) } } } } else { # LBE is non-nullifying or the version is 0 for(l in 1:L) { lx <- lx.vec[l]; Alx <- tcrossprod(Fx[lx, ]) for (k in 1:K) { kx <- kx.vec[k]; v <- c(kx, lx) alpha <- stepsize(M, w, v) w[kx] <- w[kx] - alpha; w[lx] <- w[lx] + alpha M <- M + alpha * (Alx - tcrossprod(Fx[kx, ])) } } } supp <- index[w > del]; K <- length(supp); w.supp <- w[supp] M.inv <- solve(M); a.fun <- ((Fx %*% M.inv)^2) %*% one ord.ind <- (1:n)[a.fun >= -sort(-a.fun, partial=L)[L]] ord <- ord.ind[order(a.fun[ord.ind], decreasing=TRUE)] # The two lines above can be replaced by simpler but usually # somewhat slower ord <- order(a.fun, decreasing=TRUE)[1:L] lx.vec <- sample(ord); kx.vec <- sample(supp) tm <- as.numeric(proc.time()[3]) eff.act <- sum(diag(M.inv)) / a.fun[ord[1]] print(paste("A.REX Time:", round(tm-start, 2), "Efficiency:", round(eff.act, 9)), quote=FALSE) if((1/eff.act < eff.inv) || (n.iter >= it.max) || (tm > start + t.max)) break } t.act <- round(as.numeric(proc.time()[3]) - start, 2) Phi.best <- m/sum(diag(M.inv)); eff.best <- sum(diag(M.inv))/a.fun[ord[1]] print(paste("A.REX finished at", Sys.time()), quote = FALSE) print(paste("A-criterion value:", Phi.best), quote = FALSE) print(paste("Efficiency at least:", eff.best), quote = FALSE) return(list(w.best=w, Phi.best=Phi.best, eff.best=eff.best, iter=n.iter, t.act=t.act)) } ################################################################### od.A.SOCP <- function(Fx, b1, A1, b2, A2, b3, A3, bin, type, gap, t.max){ # Note: this function should be only run via the wrappers od.SOCP or od.MISOCP. if(!requireNamespace('gurobi', quietly=TRUE)) stop("(MI)SOCP requires package gurobi.") if(!requireNamespace('Matrix', quietly=TRUE)) stop("(MI)SOCP requires package Matrix.") n <- nrow(Fx); m <- ncol(Fx) k1 <- length(b1); k2 <- length(b2); k3 <- length(b3) b <- c(b1, b2, b3); A <- rbind(A1, A2, A3) sense <- c(rep("<=", k1), rep(">=", k2), rep("=", k3)) model <- list(); model$modelsense <- "min" da <- dim(A) nc <- da[1] #number of constraints M <- 4 * n + m * n #number of parameters # constraints Aeq <- matrix(0, nrow = m * m + nc + 2 * n, ncol = M) A2 <- matrix(0, nrow = nc, ncol = M) A2[, 1:n] <- A k <- 1 for (i in 1 : (m * m)){ for (j in 1 : n){ u <- i %% m if (u == 0) u <- m Aeq[i, (3 + k) * n + j] = 0.5 * Fx[j, u] } if (i %% m == 0) k <- k + 1 } A3 <- matrix(0, nrow = 2 * n, ncol = M) for (i in 1:n){ A3[i, i] <- 1 A3[i, n + i] <- 1 A3[i, 2 * n + i] <- -1 A3[n + i, i] <- 1 A3[n + i, n + i] <- -1 A3[n + i, 3 * n + i] <- 1 } beq <- rep(0, m * m + nc + 2 * n) beq[1:(m^2)] <- as.vector(diag(m)) for (i in 1:nc){ Aeq[m * m + i, ] <- A2[i, ] beq[m * m + i] <- b[i] } for (i in 1:(2 * n)) Aeq[m * m + nc + i, ] <- A3[i, ] model$A <- Matrix::Matrix(Aeq, sparse=TRUE) model$rhs <- beq sns <- rep("=", m * m + nc + 2 * n) for (i in (m * m + 1):(m * m + nc)) sns[i] <- sense[i-m^2] model$sense <- sns # cone constraints model$quadcon <- vector("list", n) for (i in 1:n){ a <- 4:(3 + m) qc <- list() u <- length(c(2 * n + i, a * n + i, 3 * n + i)) qc$Qc <- spMatrix(M, M, c(2 * n + i, a * n + i, 3 * n + i), c(2 * n + i, a * n + i, 3 * n + i), c(-1, rep(1, u-1))) qc$rhs <- 0 model$quadcon[[i]] <- qc } # model$cones <- vector("list", n) # for (i in 1:n){ # a <- 4:(3 + m) # model$cones[[i]] <- c(2 * n + i, a * n + i, 3 * n + i) # } # lower bound lb <- rep(-Inf, M) lb[1:(3 * n)] <- 0 model$lb <- lb # objective c <- rep(0, M) c[(n + 1):(2 * n)] <- 1 model$obj <- c # variable types vtypes <- rep("C", M) if (type=="exact") { for (i in 1:n) vtypes[i] <- "I" if(bin) vtypes[1:n] <- "B" } model$vtypes <- vtypes if(!is.null(gap)) { params <- list(TimeLimit=t.max, MIPGap=gap) } else { params <- list(TimeLimit=t.max) } result <- gurobi::gurobi(model, params); gc() res <- as.vector(result$x) if (!is.vector(res) || !is.numeric(res) || !all(is.finite(res)) || length(res)1+1e-6)) w.best <- NULL } return(list(w.best=w.best, status=result$status)) } #################################################### od.AA <- function (Fx, crit="D", eta=NULL, w1=NULL, alg="REX", eff=1-1e-06, it.max=Inf, t.max=60) { # Purpose: # Computes a normalized optimal approximate design of experiments # on a finite design space wrt to the D, A, or I criterion # under the standard size constraint # Uses one of computational methods REX, MUL, or VDM # # Arguments: # Fx ... n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters # crit ... optimality criterion to be used # Possible choices are "D", "A", "I" # eta ... the summation weights for the criterion of I-optimality # if crit is not I then eta is ignored # if crit==I and eta=NULL, the program sets eta=(1,...,1) # w1 ... initial design; it must have a non-singular information matrix # it should have a small support (e.g., m) if alg=REX # it should be supported on all design points if alg=MUL # if w1==NULL, a non-singular initial design is generated by od.PINs # alg ... computational method; possible choices are "REX", "MUL", "VDM" # see the procedures od.D.REX, od.A.REX, od.D.MUL, od.A.MUL # od.D.VDM and od.A.VDM for details # eff ... threshold on the minumim design efficiency to stop the computation # it.max ... maximum allowed number of iterations of the method # t.max ... threshold for the maximum computation time # # Output is the list with components: # call ... the call of the function # w.best ... the resulting approximate design # supp ... the indices of the support of w.best # w.supp ... the weights of w.best on the support # M.best ... the information matrix of w.best # Phi.best ... the criterion value of w.best # eff.best ... a lower bound on the eff of w.best wrt an ideal optimal design # n.iter ... number of iterations performed # t.act ... the actual time of the computation # # Note: # Since the result is a normalized approximate design, it only gives # recommended *proportions* of trials in individual design points # To convert it to an exact design with N trials, w.best must be multiplied # by N and the results should be properly rounded to the neighbouring integers # # Authors of the code: # Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Input check if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be greater than or equal to the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be at least 2. The case of a single parameter is trivial.") if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") if (is.null(eta)) eta <- rep(1, n) if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative vector of length nrow(Fx).") if (!is.null(w1)) if(!is.vector(w1) || !is.numeric(w1) || !all(is.finite(w1)) || min(w1)<0 || max(w1)<=0) stop("w1 must be NULL or a non-negative, nonzero finite numeric vector of length nrow(Fx).") if (!is.null(w1)) w1 <- w1 / sum(w1) if (!is.vector(alg) || !is.character(alg) || (length(alg)!=1) || (nchar(alg)<1) || !is.element(alg, c("REX", "MUL", "VDM"))) stop("alg must be 'REX', 'MUL', or 'VDM'.") if (!is.vector(eff) || !is.numeric(eff) || !all(is.finite(eff)) || (length(eff) != 1) || (eff < 0) || (eff > 1)) stop("eff must be a real number in the interval [0,1].") # Use the appropriate engine if (crit=="D") { if (alg=="REX") res <- od.D.REX(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="MUL") res <- od.D.MUL(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="VDM") res <- od.D.VDM(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) } if (crit=="A") { if (alg=="REX") res <- od.A.REX(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="MUL") res <- od.A.MUL(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="VDM") res <- od.A.VDM(Fx=Fx, w1=w1, eff=eff, it.max=it.max, t.max=t.max) } if (crit=="I") { if (alg=="REX") res <- od.A.REX(Fx=Fx.ItoA(Fx, eta), w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="MUL") res <- od.A.MUL(Fx=Fx.ItoA(Fx, eta), w1=w1, eff=eff, it.max=it.max, t.max=t.max) if (alg=="VDM") res <- od.A.VDM(Fx=Fx.ItoA(Fx, eta), w1=w1, eff=eff, it.max=it.max, t.max=t.max) } # Output supp <- (1:n)[res$w.best > 0] res <- list(call=cl, w.best=res$w.best, supp=supp, w.supp=res$w.best[supp], M.best=infmat(Fx, res$w.best), Phi.best=res$Phi.best, eff.best=res$eff.best, n.iter=res$n.iter, t.act=res$t.act) return(res) } ######################################################### od.AQuA <- function(Fx, b1=NULL, A1=NULL, b2=NULL, A2=NULL, b3=NULL, A3=NULL, w0=NULL, bin=FALSE, crit="D", eta=NULL, M.anchor=NULL, ver="+", conic=TRUE, t.max=120) { # Purpose: # Computes an exact experimental design with the aim to maximize crit. # Uses the AQuA principle of local quadratic approximation of crit. # See the papers Harman & Filova 2014 and Filova & Harman 2019 for details # # Authors: # Radoslav Harman, Lenka Filova # # Arguments: # Fx: n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters. # b1, A1: The real vector of length k1 and the k1 times n matrix of real numbers. # If both are not NULL, the constraints A1*w<=b1 are added to the set of all design constraints. # b2, A2: The real vector of length k2 and the k2 times n matrix of real numbers. # If both are not NULL, the constraints A2*w>=b2 are added to the set of all design constraints. # b3, A3: The real vector of length k2 and the k2 times n matrix of real numbers. # If both are not NULL, the constraints A3*w==b3 are added to the set of all design constraints. # Note 1: At least one of b1, b2, b3 must be non-NULL. # Note 2: If bi is non-NULL and Ai is NULL for some i then Ai is set to be the vector of ones. # Note 3: if bi is NULL for some i then Ai is ignored. # w0: The non-negative real vector of length n representing the design to be augmented. # If NULL, w0 is set to the vector of zeros. # If not NULL, the constraints w>=w0 are added to the set of all design constraints. # bin: Should the design be binary, i.e., without repeated use of design points? # crit: The optimality criterion to be used. Possible choices are "D", "A", "I" # eta: The summation weights for the criterion of I-optimality # If crit is not "I" then eta is ignored. # If crit is "I" and eta is NULL, eta is set to be the vector of ones. # M.anchor: a non-singular positive definite "anchor" matrix for the quadratic approximation, # typically the information matrix of an optimal approximate design; see the papers. # If M.anchor is NULL, the optimal approximate matrix M.anchor is pre-computed. # ver: Version of the quadratic approximation; possible values are "+", "0", "-". # conic: If TRUE, the conic form of the problem is used; see the second paper. # Otherwise, the IQP formulation of the problem is used. # Note that "conic" is suitable for relatively larger n and relatively smaller m. # The IQP form is more suitable for relatively smaller n and relatively larger m. # t.max: approximate threshold for the maximum computation time # # Output is the list with components: # call: the call of the function # w.best: the resulting exact design, or NULL if the computation fails # supp: the indices of the support of w.best (or NULL) # w.supp: the weights of w.best on the support (or NULL) # M.best: the information matrix of w.best (or NULL) # Phi.best: the criterion value of w.best (or 0) # eff.best: the ration of Phi(M(w.best)) and Phi(M.anchor) (or 0) # Typically, this is a lower bound on the actual efficiency of w.best # status: the status of the optimization solver # See http://www.gurobi.com/documentation/8.1/refman/optimization_status_codes.html # t.act: the actual time of the computation cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) ## Verify Fx if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be >= the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be >= 2. Use od.m1 for models with a one-dimensional parameter.") ## Verify b1, b2, b3 if (is.null(b1) & is.null(b2) & is.null(b3)) stop("At least one of b1, b2, b3 must be non-NULL.") k1 <- k2 <- k3 <- 0 ## Verify b1 if (!is.null(b1)) { if (!is.vector(b1) || !is.numeric(b1) || !all(is.finite(b1))) stop("b1 must be either NULL, or a real number, or a vector of real numbers.") k1 <- length(b1) if (is.null(A1)) { if (length(b1)!=1) stop("If b1 is non-NULL and A1 is NULL then b1 must be a real number.") A1 <- matrix(1, nrow=k1, ncol=n) } } ## Verify A1 if (!is.null(A1)) { if (!is.matrix(A1) || !is.numeric(A1) || !all(is.finite(A1))) stop("A1 must be NULL, or a matrix of real numbers.") if (nrow(A1)!=length(b1)) stop("The dimension nrow(A1) must be equal to the number length(b1) of constraints.") if (ncol(A1)!=n) stop("The dimension ncol(A1) must be equal to the number nrow(Fx) of design points.") } ## Verify b2 if (!is.null(b2)) { if (!is.vector(b2) || !is.numeric(b2) || !all(is.finite(b2))) stop("b2 must be either NULL, or a real number, or a vector of real numbers.") k2 <- length(b2) if (is.null(A2)) { if (length(b2)!=1) stop("If b2 is non-NULL and A1 is NULL then b1 must be a real number.") A2 <- matrix(1, nrow=k2, ncol=n) } } ## Verify A2 if (!is.null(A2)) { if (!is.matrix(A2) || !is.numeric(A2) || !all(is.finite(A2))) stop("A2 must be NULL, or a matrix of real numbers.") if (nrow(A2)!=length(b2)) stop("The dimension nrow(A2) must be equal to the number length(b2) of constraints.") if (ncol(A2)!=n) stop("The dimension ncol(A2) must be equal to the number nrow(Fx) of design points.") } ## Verify b3 if (!is.null(b3)) { if (!is.vector(b3) || !is.numeric(b3) || !all(is.finite(b3))) stop("b3 must be either NULL, or a real number, or a vector of real numbers.") k3 <- length(b3) if (is.null(A3)) { if (length(b3)!=1) stop("If b3 is non-NULL and A1 is NULL then b3 must be a real number.") A3 <- matrix(1, nrow=k3, ncol=n) } } ## Verify A3 if (!is.null(A3)) { if (!is.matrix(A3) || !is.numeric(A3) || !all(is.finite(A3))) stop("A3 must be NULL, or a matrix of real numbers.") if (nrow(A3)!=length(b3)) stop("The dimension nrow(A3) must be equal to the number length(b3) of constraints.") if (ncol(A3)!=n) stop("The dimension ncol(A3) must be equal to the number nrow(Fx) of design points.") } ## Verify w0 if (is.null(w0)) w0 <- rep(0, n) if (!is.vector(w0) || !is.numeric(w0) || !all(is.finite(w0)) || (min(w0) < 0)) stop("w0 must be NULL, or a vector of non-negative real numbers.") if (length(w0) != n) stop("The length of w0 must be equal to the number nrow(Fx) of design points.") ## Verify bin if (!is.vector(bin) || !is.logical(bin) || length(bin)>1) stop("bin must be a single logical value (TRUE or FALSE).") ## Verify crit if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") ## Verify eta if (!is.null(eta)) { if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative vector of length nrow(Fx).") } ## Verify M.anchor if (!is.null(M.anchor)) { if (!is.matrix(M.anchor) || !is.numeric(M.anchor) || !all(is.finite(M.anchor)) || nrow(M.anchor)!=m || ncol(M.anchor)!=m) stop("M.anchor must be either NULL, or an ncol(Fx) times ncol(Fx) matrix of real numbers.") if (max(abs(M.anchor-t(M.anchor))) > 1e-6 || min(Re(eigen(M.anchor)$values)) < 1e-6) stop("M.anchor must be a symmetric positive definite matrix.") } ## Verify type if (!is.vector(ver) || !is.character(ver) || (length(ver) != 1) || !is.element(ver, c("+", "-", "0"))) stop("ver must be '+', '0', or '-'.") ## Verify conic if (!is.vector(conic) || !is.logical(conic) || length(conic)>1) stop("conic must be a single logical value (TRUE or FALSE).") ## Verify t.max if (!is.vector(t.max) || !is.numeric(t.max) || (length(t.max) != 1) || (t.max <= 0)) stop("t.max must be a positive real number or Inf.") if (!is.null(w0) && sum(w0)>0.5) { A3 <- rbind(A3, diag(n)[w0>0.5,]) b3 <- c(b3, w0[w0>0.5]) } # Print info info <- paste("Running od.AQuA for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, ", m=", m, ",", sep="") info <- paste(info, ", k=", k1+k2+k3, ".", sep="") print(info, quote=FALSE) info <- paste("Other parameters: ", sum(b1), ",", sum(b2), ",", sum(b3), ":", sum(w0), ",", bin, ",", crit, ",", sum(eta), ",", sum(M.anchor), ",", ver, ",", conic, sep="") print(info, quote=FALSE) # Track time start <- as.numeric(proc.time()[3]) # Run the appropriate engine if (crit=="D") res <- od.D.AQuA(Fx, b1, A1, b2, A2, b3, A3, bin, M.anchor, ver, conic, t.max) if (crit=="A") res <- od.A.AQuA(Fx, b1, A1, b2, A2, b3, A3, bin, M.anchor, ver, conic, t.max) if (crit=="I") res <- od.A.AQuA(Fx.ItoA(Fx, eta), b1, A1, b2, A2, b3, A3, bin, M.anchor, ver, conic, t.max) # Calculate the output from the result of the engine w <- res$w.best if (!is.null(w)) { supp <- (1:n)[w > 0.5]; w.supp <- w[supp] M.best <- infmat(Fx, w) Phi.best <- optcrit(Fx, w, crit=crit) } else { supp <- NULL; w.supp <- NULL M.best <- NULL; Phi.best <- 0 } t.act <- round(as.numeric(proc.time()[3]) - start, 2) return(list(call=cl, w.best=w, supp=supp, w.supp=w.supp, M.best=M.best, Phi.best=Phi.best, status=res$status, t.act=t.act)) } ####################################################### od.D.AQuA <- function(Fx, b1, A1, b2, A2, b3, A3, bin, M.anchor, ver, conic, t.max) { # Note: this function should be only run via the wrapper od.AQuA. if(!requireNamespace('gurobi', quietly=TRUE)) stop("AQuA requires the gurobi package.") n <- nrow(Fx); m <- ncol(Fx) k1 <- length(b1); k2 <- length(b2); k3 <- length(b3) b <- c(b1, b2, b3); A <- rbind(A1, A2, A3) sense <- c(rep("<", k1), rep(">", k2), rep("=", k3)) if (is.null(M.anchor)) { if (length(b)==1 && all(A==matrix(1, ncol=n, nrow=1))) { M.anchor <- infmat(Fx, b * od.D.REX(Fx, t.max=t.max/5)$w.best) } else { w.app <- od.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin, t.max=t.max/3)$w.best if (is.null(w.app)) { stop("Failure in the computation of M.anchor (and no M.anchor supplied by the user).") } else { M.anchor <- infmat(Fx, w.app) } } } M1 <- solve(M.anchor) if(conic) { if(!requireNamespace('Matrix', quietly=TRUE)) stop("Conic version of AQuA requires package Matrix.") if(!requireNamespace('matrixcalc', quietly=TRUE)) stop("Conic version of AQuA requires package matrixcalc.") s <- m*(m+1)/2; Gm <- matrixcalc::duplication.matrix(m); eps <- 1e-12 vec.M1 <- matrixcalc::vec(M1) tilde.h <- t(Gm) %*% vec.M1 if(ver=="+") tilde.Q <- 0.5 * t(Gm) %*% (-vec.M1 %*% t(vec.M1) / m + M1 %x% M1) %*% Gm if(ver=="-") tilde.Q <- (1/6) * t(Gm) %*% (vec.M1 %*% t(vec.M1) / m + M1 %x% M1) %*% Gm if(rcond(tilde.Q) > eps) { tilde.C <- t(chol(tilde.Q)) } else { eig <- eigen(tilde.Q, symmetric=TRUE) mx <- max(eig$values); ts <- sum(eig$values/mx > eps) tilde.C <- eig$vectors[,1:ts] %*% diag(sqrt(eig$values[1:ts])) } diff <- max(abs(tilde.Q - tilde.C%*%t(tilde.C))) print(paste("Decomposition instability:", diff)) H <- matrix(0, nrow=n, ncol=s); k <- 0 for(i1 in 1:m) { for(i2 in i1:m) { k <- k + 1; H[, k] <- Fx[, i1] * Fx[, i2] } } h <- H %*% tilde.h; S <- H %*% tilde.C; u <- ncol(S) model <- list(); k <- nrow(A); np <- n + u + 3 Aeq <- matrix(0, nrow = k + 2 + u, ncol = np); Aeq[1:k, 1:n] <- A Aeq[k + 1, n + 1] <- 1; Aeq[k + 1, n + 3] <- -1/sqrt(2) Aeq[k + 2, n + 2] <- 1; Aeq[k + 2, n + 3] <- 1/sqrt(2) Aeq[(k + 3):(k + 2 + u), 1:n] <- t(S) Aeq[(k + 3):(k + 2 + u), (n + 4):np] <- -diag(u) model$A <- Matrix::Matrix(Aeq, sparse=TRUE) beq <- rep(0, k + 2 + u); beq[1:k] <- b beq[k + 1] <- 1/2/sqrt(2); beq[k + 2] <- 1/2/sqrt(2) model$rhs <- beq model$sense <- c(rep("<", k1), rep(">", k2), rep("=", k3 + u + 2)) model$cones <- list(as.list(c(n+1, n+2, (n+4):np))) lb <- rep(-Inf, np); lb[1:n] <- 0; model$lb <- lb ob <- rep(0, np); ob[1:n] <- -h; ob[n + 3] <- 1 model$obj <- ob; vtypes <- rep("C", np) if(bin) {vtypes[1:n] <- "B"} else {vtypes[1:n] <- "I"} model$vtypes <- vtypes; model$modelsense <- "min" params <- list(TimeLimit=t.max) result <- suppressWarnings(gurobi::gurobi(model, params)); gc() } else { f <- rep(0, n); H <- matrix(0, nrow=n, ncol=n) for (i in 1:n) f[i] <- -t(Fx[i, ]) %*% M1 %*% Fx[i, ] if(ver=="+") { for (i in 1:n) for (j in i:n) { H[j, i] <- H[i, j] <- 0.5 * (t(Fx[i, ]) %*% M1 %*% Fx[j, ])^2 - 0.5/m * t(Fx[i, ]) %*% M1 %*% Fx[i, ] * t(Fx[j, ]) %*% M1 %*% Fx[j, ] } } if(ver=="-") { for (i in 1:n) for (j in i:n) { H[j, i] <- H[i, j] <- (1/6) * (t(Fx[i, ]) %*% M1 %*% Fx[j, ])^2 + (1/6/m) * t(Fx[i, ]) %*% M1 %*% Fx[i, ] * t(Fx[j, ]) %*% M1 %*% Fx[j, ] } } model <- list(); model$obj <- f model$Q <- H; model$A <- Matrix::Matrix(A, sparse=TRUE); model$rhs <- c(b) model$sense <- sense; model$lb <- rep(0, n) if(bin) {model$vtypes <- "B"} else {model$vtypes <- "I"} params <- list(TimeLimit=t.max) result <- gurobi::gurobi(model, params) } res <- result$x if (!is.vector(res) || !is.numeric(res) || !all(is.finite(res)) || length(res)1)) w.best <- NULL } return(list(w.best=w.best, status=result$status)) } ############################################################ od.D.KL <- function(Fx, N, Phi.app=NULL, w1=NULL, K=NULL, L=NULL, rest.max=Inf, t.max=60) { # Purpose: # A version of the KL exchange algorithm for computing # D-efficient exact designs under the standard "size" constraint # # Arguments: # Fx ... n times m model matrix, where n>=2 is the size of the design space # and m>=2 (m<=n) is the number of model parameters # N ... required size of the design (the number of experimental trials) # Phi.app ... optimal value of the corresponding approximate (relaxed) problem # If NULL, the value is pre-computed by REX # w1 ... initial design used for the first restart of the algorithm # It should be either NULL or a non-singular design of size N # K ... (upper bound on the) number of "least promising" support points # of the current design, for which exchanges are attempted # L ... (upper bound on the) number of "most promising" candidate # design points for which exchanges are attempted # If K is NULL or L is NULL, a reasonable value is assigned automatically # rest.max ... limit on the number of restarts # t.max ... approximate upper limit on the time of computation # If the algorithm stops in a local optimum before the time ellapsed, # the computation is restarted # # Output is a list of: # w.best ... the best design found within the time limit # Phi.best ... the value of the D-criterion of w.best # eff.best ... the efficiency of w.best wrt the approximate optimum # n.ex ... the number of exchanges # n.rest ... the number of restarts executed # t.act ... the actual time taken by the computation # # Notes: # Works reasonably well at least for problems with n<=10000 and m<=10 # # Authors of this implementation: # Radoslav Harman and Lenka Filova eps <- 1e-09; del <- 0.1 n <- nrow(Fx); m <- ncol(Fx); one.m <- rep(1, m) next.sec <- 0; n.ex <- 0; n.rest <- 0 if (is.null(Phi.app)) { Phi.app <- N * od.D.REX(Fx, od.PIN(Fx, m, rnd=0)$w)$Phi.best } start <- as.numeric(proc.time()[3]) info <- paste("Running od.D.KL for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, " and m=", m, ".", sep="") print(info, quote=FALSE) if (is.null(K)) K <- max(10, min(ceiling(sqrt(c(N, n))))) if (is.null(L)) L <- max(10, ceiling(sqrt(n))) print(paste("Setting K=", K, ", L=", L, sep=""), quote=FALSE) A <- array(0, dim=c(n, m, m)) for (i in 1:n) A[i, , ] <- tcrossprod(Fx[i, ]) finish.all <- FALSE; detM.best <- 0 while (!finish.all && n.rest m) { for (k in 1:(N - m)) { i.supp <- sample(1:n, 1) w[i.supp] <- w[i.supp] + 1 M <- M + A[i.supp, , ] } } } else { w <- w1; M <- crossprod(sqrt(w) * Fx) w1 <- NULL # Start the second run from a "pinned" design } detM <- det(M) if (detM.best < detM){ w.best <- w; M.best <- M; detM.best <- detM } finish.all <- finish <- as.numeric(proc.time()[3]) > start + t.max while (!finish){ tm <- as.numeric(proc.time()[3]) - start if (tm > next.sec){ Phi.best <- detM.best^(1/m) info <- paste("D.KL Time:", round(tm, 1), " Value:", round(Phi.best, 6), "Efficiency:", round(Phi.best/Phi.app, 6)) print(info, quote=FALSE); next.sec <- ceiling(tm) } d.fun <- ((Fx %*% t(chol(solve(M))))^2) %*% one.m; supp <- (1:n)[w > 0.5] Kact <- min(length(supp), K); Lact <- min(n, L) Kind <- supp[order(d.fun[supp])][1:Kact] Lind <- order(d.fun)[(n-Lact+1):n] imp <- FALSE for (iL in Lact:1){ for (iK in 1:Kact){ M.temp <- M + A[Lind[iL], , ] - A[Kind[iK], , ] detM.temp <- det(M.temp) if(detM.temp > detM){ w[Lind[iL]] <- w[Lind[iL]] + 1 w[Kind[iK]] <- w[Kind[iK]] - 1 M <- M.temp; detM <- detM.temp n.ex <- n.ex + 1; imp <- TRUE break } } if (imp) break } if (as.numeric(proc.time()[3]) > start + t.max) finish.all <- TRUE if (finish.all || !imp) finish <- TRUE } M <- crossprod(sqrt(w) * Fx); detM <- det(M) if (detM > detM.best){ w.best <- w; M.best <- M; detM.best <- detM } } t.act <- round(as.numeric(proc.time()[3]) - start, 2) info <- paste("D.KL finished after", t.act, "seconds at", Sys.time()) print(info, quote=FALSE) info <- paste("with", n.rest, "restarts and", n.ex, "exchanges.") print(info, quote=FALSE) Phi.best <- optcrit(Fx, w.best, crit="D") return(list(w.best=w.best, Phi.best=Phi.best, eff.best=Phi.best/Phi.app, n.ex=n.ex, n.rest=n.rest, t.act=t.act)) } ############################################################# od.D.REX <- function (Fx, w1=NULL, ver=1, gamma=4, eff=1-1e-9, it.max=Inf, t.max=30) { # Computes a D-optimal approximate design of experiments on a finite design space # Can also be used to compute the minimum-volume data-enclosing ellipsoid # Based on the paper: https://arxiv.org/abs/1801.05661 # Authors: Radoslav Harman, Lenka Filova, Peter Richtarik # # Arguments: # Fx ... n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters # w1 ... initial design; it must have a non-singular information matrix # and for the sake of efficiency, it should have a small support (e.g., m) # if NULL, a non-singular initial design is generated by PINs. # ver ... version of the algorithm (1 = with and 0 = without the nullity control) # gamma ... positive parameter regulating the size of the exchange batch # eff ... threshold on the minumim design efficiency to stop the computation # it.max ... maximum allowed number of iterations # t.max ... threshold for the maximum computation time # # Output is the list with components: # w.best ... the resulting approximate design # Phi.best ... the D-criterion value of w.best # eff.best ... a lower bound on the eff. of w.best wrt an ideal D-optimal design # n.iter ... number of performed iterations # t.act ... the actual time of the computation # Example: A D-optimal design on 100000 random regressors in R^10 # n <- 100000; m <- 10; F.rnd <- matrix(rnorm(n*m), ncol=m) # res <- od.D.REX(F.rnd); supp <- res$w.best>0 # print(cbind((1:n)[supp], res$w.best[supp])) start <- as.numeric(proc.time()[3]); del <- 1e-14; eps <- 1e-24 Fx <- as.matrix(Fx); n <- nrow(Fx); m <- ncol(Fx) info <- paste("Running od.D.REX for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, " and m=", m, ".", sep="") print(info, quote=FALSE) eff.inv <- 1/eff; n.iter <- 0; L <- min(n, gamma*m) lx.vec <- rep(0, L); index <- 1:n; one <- rep(1, m) if(is.null(w1)) w1 <- od.PIN(Fx, m)$w supp <- (1:n)[w1>0]; K <- length(supp); Fx.supp <- Fx[supp, ] w <- rep(0, n); w[supp] <- 1 / length(supp); w.supp <- w[supp] M <- crossprod(sqrt(w.supp) * Fx.supp) d.fun <- ((Fx %*% t(chol(solve(M))))^2) %*% one / m ord <- order(d.fun, decreasing=TRUE) lx.vec <- sample(ord[1:L]); kx.vec <- sample(supp) while (TRUE) { n.iter <- n.iter + 1; ord1 <- which.min(d.fun[supp]) kb <- supp[ord1]; lb <- ord[1]; v <- c(kb, lb) cv <- Fx[v, ] %*% solve(M, t(Fx[v, ])) alpha <- 0.5 * (cv[2, 2] - cv[1, 1])/(cv[1, 1] * cv[2, 2] - cv[1, 2]^2 + eps) alpha <- min(w[kb], alpha) w[kb] <- w[kb] - alpha; w[lb] <- w[lb] + alpha M <- M + alpha * (tcrossprod(Fx[lb, ]) - tcrossprod(Fx[kb, ])) if((w[kb] < del) && (ver==1)) { # LBE is nullifying and the version is 1 for(l in 1:L) { lx <- lx.vec[l]; Alx <- tcrossprod(Fx[lx, ]) for (k in 1:K) { kx <- kx.vec[k]; v <- c(kx, lx) cv <- Fx[v, ] %*% solve(M, t(Fx[v, ])) alpha <- 0.5 * (cv[2, 2] - cv[1, 1])/(cv[1, 1] * cv[2, 2] - cv[1, 2]^2 + eps) alpha <- min(w[kx], max(-w[lx], alpha)) wkx.temp <- w[kx] - alpha; wlx.temp <- w[lx] + alpha if((wkx.temp < del) || (wlx.temp < del)) { w[kx] <- wkx.temp; w[lx] <- wlx.temp M <- M + alpha * (Alx - tcrossprod(Fx[kx, ])) } } } } else { # LBE is non-nullifying or the version is 0 for(l in 1:L) { lx <- lx.vec[l]; Alx <- tcrossprod(Fx[lx, ]) for (k in 1:K) { kx <- kx.vec[k]; v <- c(kx, lx) cv <- Fx[v, ] %*% solve(M, t(Fx[v, ])) alpha <- 0.5 * (cv[2, 2] - cv[1, 1])/(cv[1, 1] * cv[2, 2] - cv[1, 2]^2 + eps) alpha <- min(w[kx], max(-w[lx], alpha)) w[kx] <- w[kx] - alpha; w[lx] <- w[lx] + alpha M <- M + alpha * (Alx - tcrossprod(Fx[kx, ])) } } } supp <- index[w > del]; K <- length(supp); w.supp <- w[supp] d.fun <- ((Fx %*% t(chol(solve(M))))^2) %*% one / m ord.ind <- (1:n)[d.fun >= -sort(-d.fun, partial=L)[L]] ord <- ord.ind[order(d.fun[ord.ind], decreasing=TRUE)] # The two lines above can be replaced by simpler but usually # somewhat slower ord <- order(d.fun, decreasing=TRUE)[1:L] lx.vec <- sample(ord); kx.vec <- sample(supp) tm <- as.numeric(proc.time()[3]) eff.act <- 1 / d.fun[ord[1]] print(paste("D.REX Time:", round(tm-start, 2), "Efficiency:", round(eff.act, 9)), quote=FALSE) if ((d.fun[ord[1]] < eff.inv) || (n.iter >= it.max) || (tm > start + t.max)) break } t.act <- round(as.numeric(proc.time()[3]) - start, 2) Phi.best <- det(M)^(1/m); eff.best <- 1/d.fun[ord[1]] print(paste("D.REX finished at", Sys.time()), quote = FALSE) print(paste("D-criterion value:", Phi.best), quote = FALSE) print(paste("Efficiency at least:", eff.best), quote = FALSE) return(list(w.best=w, Phi.best=Phi.best, eff.best=eff.best, n.iter=n.iter, t.act=t.act)) } ####################################################### od.D.SOCP <- function(Fx, b1, A1, b2, A2, b3, A3, bin, type, gap, t.max) { # Note: this function should be only run via the wrappers od.SOCP or od.MISOCP. if(!requireNamespace('gurobi', quietly=TRUE)) stop("(MI)SOCP requires package gurobi.") if(!requireNamespace('Matrix', quietly=TRUE)) stop("(MI)SOCP requires package Matrix.") geomeanG <- function(m, n, M) { cones.index <- NULL if(m == 2) { cones.index <- vector("list", 1) cones.index[[1]] <- c(M + 1, M + 2, M + 3) newvars <- 3; ncg <- 2; cn <- 1 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 } if(m == 3) { cones.index <- vector("list", 3) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) ncg <- 6; newvars <- 9; cn <- 3 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, M + 9] <- 1/2; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, M + 9] <- -1/2; Ag[4, M + 5] <- -1 Ag[5, M + 3] <- 1/2; Ag[5, M + 6] <- 1/2; Ag[5, M + 7] <- -1 Ag[6, M + 3] <- 1/2; Ag[6, M + 6] <- -1/2; Ag[6, M + 8] <- -1 } if(m == 4) { cones.index <- vector("list", 3) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) ncg <- 6; newvars <- 9; cn <- 3 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, M + 3] <- 1/2; Ag[5, M + 6] <- 1/2; Ag[5, M + 7] <- -1 Ag[6, M + 3] <- 1/2; Ag[6, M + 6] <- -1/2; Ag[6, M + 8] <- -1 } if(m == 5) { cones.index <- vector("list", 6) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cn <- 6; newvars <- 18; ncg <- 12 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, M + 18] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, M + 18] <- -1; Ag[6, M + 8] <- -1 Ag[7, M + 3] <- 1/2; Ag[7, M + 6] <- 1/2; Ag[7, M + 10] <- -1 Ag[8, M + 3] <- 1/2; Ag[8, M + 6] <- -1/2; Ag[8, M + 11] <- -1 Ag[9, M + 9] <- 1/2; Ag[9, M + 18] <- 1/2; Ag[9, M + 13] <- -1 Ag[10, M + 9] <- 1/2; Ag[10, M + 18] <- -1/2; Ag[10, M + 14] <- -1 Ag[11, M + 12] <- 1/2; Ag[11, M + 15] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 12] <- 1/2; Ag[12, M + 15] <- -1/2; Ag[12, M + 17] <- -1 } if(m == 6) { cones.index <- vector("list", 6) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cn <- 6; newvars <- 18; ncg <- 12 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, n + 5 * m + 6] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, n + 5 * m + 6] <- -1; Ag[6, M + 8] <- -1 Ag[7, M + 3] <- 1/2; Ag[7, M + 6] <- 1/2; Ag[7, M + 10] <- -1 Ag[8, M + 3] <- 1/2; Ag[8, M + 6] <- -1/2; Ag[8, M + 11] <- -1 Ag[9, M + 9] <- 1/2; Ag[9, M + 18] <- 1/2; Ag[9, M + 13] <- -1 Ag[10, M + 9] <- 1/2; Ag[10, M + 18] <- -1/2; Ag[10, M + 14] <- -1 Ag[11, M + 12] <- 1/2; Ag[11, M + 15] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 12] <- 1/2; Ag[12, M + 15] <- -1/2; Ag[12, M + 17] <- -1 } if(m == 7) { cones.index <- vector("list", 7) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cones.index[[7]] <- c(M + 19, M + 20, M + 21) cn <- 7; newvars <- 21; ncg <- 14 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, n + 5 * m + 6] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, n + 5 * m + 6] <- -1; Ag[6, M + 8] <- -1 Ag[7, n + 6 * m + 7] <- 1; Ag[7, M + 21] <- 1; Ag[7, M + 10] <- -1 Ag[8, n + 6 * m + 7] <- 1; Ag[8, M + 21] <- -1; Ag[8, M + 11] <- -1 Ag[9, M + 3] <- 1/2; Ag[9, M + 6] <- 1/2; Ag[9, M + 13] <- -1 Ag[10, M + 3] <- 1/2; Ag[10, M + 6] <- -1/2; Ag[10, M + 14] <- -1 Ag[11, M + 9] <- 1/2; Ag[11, M + 12] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 9] <- 1/2; Ag[12, M + 12] <- -1/2; Ag[12, M + 17] <- -1 Ag[13, M + 15] <- 1/2; Ag[13, M + 18] <- 1/2; Ag[13, M + 19] <- -1 Ag[14, M + 15] <- 1/2; Ag[14, M + 18] <- -1/2; Ag[14, M + 20] <- -1 } if(m == 8) { cones.index <- vector("list", 7) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cones.index[[7]] <- c(M + 19, M + 20, M + 21) cn <- 7; newvars <- 21; ncg <- 14 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, n + 5 * m + 6] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, n + 5 * m + 6] <- -1; Ag[6, M + 8] <- -1 Ag[7, n + 6 * m + 7] <- 1; Ag[7, n + 7 * m + 8] <- 1; Ag[7, M + 10] <- -1 Ag[8, n + 6 * m + 7] <- 1; Ag[8, n + 7 * m + 8] <- -1; Ag[8, M + 11] <- -1 Ag[9, M + 3] <- 1/2; Ag[9, M + 6] <- 1/2; Ag[9, M + 13] <- -1 Ag[10, M + 3] <- 1/2; Ag[10, M + 6] <- -1/2; Ag[10, M + 14] <- -1 Ag[11, M + 9] <- 1/2; Ag[11, M + 12] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 9] <- 1/2; Ag[12, M + 12] <- -1/2; Ag[12, M + 17] <- -1 Ag[13, M + 15] <- 1/2; Ag[13, M + 18] <- 1/2; Ag[13, M + 19] <- -1 Ag[14, M + 15] <- 1/2; Ag[14, M + 18] <- -1/2; Ag[14, M + 20] <- -1 } if(m == 9) { cones.index <- vector("list", 11) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cones.index[[7]] <- c(M + 19, M + 20, M + 21) cones.index[[8]] <- c(M + 22, M + 23, M + 24) cones.index[[9]] <- c(M + 25, M + 26, M + 27) cones.index[[10]] <- c(M + 28, M + 29, M + 30) cones.index[[11]] <- c(M + 31, M + 32, M + 33) cn <- 11; newvars <- 33; ncg <- 22 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, n + 5 * m + 6] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, n + 5 * m + 6] <- -1; Ag[6, M + 8] <- -1 Ag[7, n + 6 * m + 7] <- 1; Ag[7, n + 7 * m + 8] <- 1; Ag[7, M + 10] <- -1 Ag[8, n + 6 * m + 7] <- 1; Ag[8, n + 7 * m + 8] <- -1; Ag[8, M + 11] <- -1 Ag[9, n + 8 * m + 9] <- 1; Ag[9, M + 33] <- 1; Ag[9, M + 13] <- -1 Ag[10, n + 8 * m + 9] <- 1; Ag[10, M + 33] <- -1; Ag[10, M + 14] <- -1 Ag[11, M + 3] <- 1/2; Ag[11, M + 6] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 3] <- 1/2; Ag[12, M + 6] <- -1/2; Ag[12, M + 17] <- -1 Ag[13, M + 9] <- 1/2; Ag[13, M + 12] <- 1/2; Ag[13, M + 19] <- -1 Ag[14, M + 9] <- 1/2; Ag[14, M + 12] <- -1/2; Ag[14, M + 20] <- -1 Ag[15, M + 15] <- 1/2; Ag[15, M + 33] <- 1/2; Ag[15, M + 22] <- -1 Ag[16, M + 15] <- 1/2; Ag[16, M + 33] <- -1/2; Ag[16, M + 23] <- -1 Ag[17, M + 18] <- 1/2; Ag[17, M + 21] <- 1/2; Ag[17, M + 25] <- -1 Ag[18, M + 18] <- 1/2; Ag[18, M + 21] <- -1/2; Ag[18, M + 26] <- -1 Ag[19, M + 24] <- 1/2; Ag[19, M + 33] <- 1/2; Ag[19, M + 28] <- -1 Ag[20, M + 24] <- 1/2; Ag[20, M + 33] <- -1/2; Ag[20, M + 29] <- -1 Ag[21, M + 27] <- 1/2; Ag[21, M + 30] <- 1/2; Ag[21, M + 31] <- -1 Ag[22, M + 27] <- 1/2; Ag[22, M + 30] <- -1/2; Ag[22, M + 32] <- -1 } if(m == 10) { cones.index <- vector("list", 11) cones.index[[1]] <- c(M + 1, M + 2, M + 3) cones.index[[2]] <- c(M + 4, M + 5, M + 6) cones.index[[3]] <- c(M + 7, M + 8, M + 9) cones.index[[4]] <- c(M + 10, M + 11, M + 12) cones.index[[5]] <- c(M + 13, M + 14, M + 15) cones.index[[6]] <- c(M + 16, M + 17, M + 18) cones.index[[7]] <- c(M + 19, M + 20, M + 21) cones.index[[8]] <- c(M + 22, M + 23, M + 24) cones.index[[9]] <- c(M + 25, M + 26, M + 27) cones.index[[10]] <- c(M + 28, M + 29, M + 30) cones.index[[11]] <- c(M + 31, M + 32, M + 33) cn <- 11; newvars <- 33; ncg <- 22 Ag <- matrix(0, nrow = ncg, ncol = M + newvars) Ag[1, n + 1] <- 1; Ag[1, n + m + 2] <- 1; Ag[1, M + 1] <- -1 Ag[2, n + 1] <- 1; Ag[2, n + m + 2] <- -1; Ag[2, M + 2] <- -1 Ag[3, n + 2 * m + 3] <- 1; Ag[3, n + 3 * m + 4] <- 1; Ag[3, M + 4] <- -1 Ag[4, n + 2 * m + 3] <- 1; Ag[4, n + 3 * m + 4] <- -1; Ag[4, M + 5] <- -1 Ag[5, n + 4 * m + 5] <- 1; Ag[5, n + 5 * m + 6] <- 1; Ag[5, M + 7] <- -1 Ag[6, n + 4 * m + 5] <- 1; Ag[6, n + 5 * m + 6] <- -1; Ag[6, M + 8] <- -1 Ag[7, n + 6 * m + 7] <- 1; Ag[7, n + 7 * m + 8] <- 1; Ag[7, M + 10] <- -1 Ag[8, n + 6 * m + 7] <- 1; Ag[8, n + 7 * m + 8] <- -1; Ag[8, M + 11] <- -1 Ag[9, n + 8 * m + 9] <- 1; Ag[9, n + 9 * m + 10] <- 1; Ag[9, M + 13] <- -1 Ag[10, n + 8 * m + 9] <- 1; Ag[10, n + 9 * m + 10] <- -1; Ag[10, M + 14] <- -1 Ag[11, M + 3] <- 1/2; Ag[11, M + 6] <- 1/2; Ag[11, M + 16] <- -1 Ag[12, M + 3] <- 1/2; Ag[12, M + 6] <- -1/2; Ag[12, M + 17] <- -1 Ag[13, M + 9] <- 1/2; Ag[13, M + 12] <- 1/2; Ag[13, M + 19] <- -1 Ag[14, M + 9] <- 1/2; Ag[14, M + 12] <- -1/2; Ag[14, M + 20] <- -1 Ag[15, M + 15] <- 1/2; Ag[15, M + 33] <- 1/2; Ag[15, M + 22] <- -1 Ag[16, M + 15] <- 1/2; Ag[16, M + 33] <- -1/2; Ag[16, M + 23] <- -1 Ag[17, M + 18] <- 1/2; Ag[17, M + 21] <- 1/2; Ag[17, M + 25] <- -1 Ag[18, M + 18] <- 1/2; Ag[18, M + 21] <- -1/2; Ag[18, M + 26] <- -1 Ag[19, M + 24] <- 1/2; Ag[19, M + 33] <- 1/2; Ag[19, M + 28] <- -1 Ag[20, M + 24] <- 1/2; Ag[20, M + 33] <- -1/2; Ag[20, M + 29] <- -1 Ag[21, M + 27] <- 1/2; Ag[21, M + 30] <- 1/2; Ag[21, M + 31] <- -1 Ag[22, M + 27] <- 1/2; Ag[22, M + 30] <- -1/2; Ag[22, M + 32] <- -1 } return(list(cones.index = cones.index, newvars = newvars, ncg = ncg, Ag = Ag, cn = cn)) } n <- nrow(Fx); m <- ncol(Fx) if(m > 10) stop("This procedure is only implemented for m<=10.") k1 <- length(b1); k2 <- length(b2); k3 <- length(b3) b <- c(b1, b2, b3); A <- rbind(A1, A2, A3) sense <- c(rep("<=", k1), rep(">=", k2), rep("=", k3)) model <- list(); model$modelsense <- "max" da <- dim(A); nc <- da[1]; M <- n + m^2 + 4*m*n geo <- geomeanG(m, n, M) varnum <- geo$newvars + M # linear constraints Aeq*x=beq nC <- c(m * m, 1/2 * m * (m - 1), m, nc, m * n, m * n, geo$ncg) sC <- sum(nC) Aeq <- matrix(0, nrow=sum(nC), ncol=varnum) for (i in 1:(m*m)) Aeq[i, i + n] <- -1 for (i in 1:m) { for (j in 1:n) { for (ii in 1:m) { Aeq[(i-1) * m + ii, (n + m^2 + m * n + (j-1) * m + ii)] <- 1/2 * Fx[j, i] } } } ind <- 1 for (i in 1:(m-1)) { for (j in i:(m-1)) { Aeq[m*m + ind, n + (i-1)*m + j + 1] <- 1 ind <- ind + 1 } } for (i in 1:m) { Aeq[m*m + m*(m-1)/2 + i, n + (i-1)*m + i] <- -1 for (j in 1:n) { Aeq[m*m + m*(m-1)/2 + i, n + m*m + (j-1)*m + i] <- 1 } } Aeq[(m ^ 2 + 1/2 * m * (m - 1) + m + 1) : (m ^ 2 + 0.5 * m * (m - 1) + m + nc), 1:n] <- A S <- sum(nC[1:4]) for (i in 1:(n*m)) { Aeq[S + i, n + m^2 + i] <- 1 Aeq[S + i, n + m^2 + 2*n*m + i] <- -1 Aeq[S + i, floor((i-1)/m) + 1] <- 1 } S <- sum(nC[1:5]) for (i in 1:(n*m)) { Aeq[S + i, n + m^2 + i] <- 1 Aeq[S + i, n + m^2 + 3*n*m + i] <- -1 Aeq[S + i, floor((i-1)/m) + 1] <- -1 } S <- sum(nC[1:6]) if(nC[7] > 0) Aeq[(S + 1):sC, ] <- geo$Ag model$A <- Matrix::Matrix(Aeq, sparse=TRUE) rhs <- rep(0, sC); rhs[(sum(nC[1:3]) + 1):(sum(nC[1:4]))] <- b model$rhs <- rhs model$sense <- rep("=", sC) for(i in (m^2 + 1/2 * m * (m-1) + 1) : (m^2 + 1/2 * m * (m-1) + m)) model$sense[i] <- "<=" for(i in (m^2 + 1/2 * m * (m-1) + m + 1) : (m^2 + 1/2 * m * (m-1) + m + nc)) model$sense[i] <- sense[i-(m^2 + 1/2 * m * (m-1) + m)] # cone constraints model$cones <- vector("list", n + geo$cn) model$cones <- geo$cones.index for (i in 1 : (m * n)) model$cones[[i + geo$cn]] <- c(n + m^2 + 2 * m * n + i, n + m ^ 2 + 3 * m * n + i, n + m ^ 2 + m * n + i) # lower and upper bounds blx <- rep(-Inf, varnum); blx[1:n] <- 0 blx[(n + m^2 + 1):(n + m^2 + n*m)] <- 0 blx[(n + m^2 + 2*m*n + 1):(n + m^2 + 3*m*n)] <- 0 # 1st indices of all cones in geomean if (m == 2) blx[M + 1] <- 0 if (m == 3) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 } if (m == 4) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 } if (m == 5) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 } if (m == 6) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 } if (m == 7) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 blx[M + 19] <- 0 } if (m == 8) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0; blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 blx[M + 19] <- 0 } if (m == 9) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 blx[M + 19] <- 0 blx[M + 22] <- 0; blx[M + 25] <- 0; blx[M + 28] <- 0 blx[M + 31] <- 0 } if (m == 10) { blx[M + 1] <- 0; blx[M + 4] <- 0; blx[M + 7] <- 0 blx[M + 10] <- 0; blx[M + 13] <- 0; blx[M + 16] <- 0 blx[M + 19] <- 0 blx[M + 22] <- 0; blx[M + 25] <- 0; blx[M + 28] <- 0 blx[M + 31] <- 0 } bux <- rep(Inf, varnum) model$lb <- blx; model$ub <- bux model$lb[1:n] <- 0 # objective cc <- rep(0, varnum); cc[varnum] <- 1 model$obj <- cc # variable types vtypes <- rep("C", varnum) if(type=="exact") { vtypes[1:n] <- "I" if(bin) vtypes[1:n] <- "B" } if(type=="approximate" && bin) model$ub[1:n] <- pmin(rep(1, n), bux[1:n]) model$vtypes <- vtypes if(!is.null(gap)) { params <- list(TimeLimit=t.max, MIPGap=gap) } else { params <- list(TimeLimit=t.max) } result <- suppressWarnings(gurobi::gurobi(model, params)); gc() res <- result$x if (!is.vector(res) || !is.numeric(res) || !all(is.finite(res)) || length(res)1+1e-6)) { print(min(w.best)) print(max(w.best)) print("Self-check failed") # Note this check should also check viability of (in)equalities with A1,A2,A3 w.best <- NULL } } return(list(w.best=w.best, status=result$status)) } ########################################################## od.KL <- function (Fx, N, Phi.app=NULL, crit="D", eta=NULL, w1=NULL, K=NULL, L=NULL, rest.max=Inf, t.max=120) { # Purpose: # Computes an optimal or near-optimal exact design of experiments on a finite # design space with respect to the criterion of D-, A-, or I-optimality # under the constraint on the size of the experiment (the number of trials) # Uses one specific version of the KL exchange algorithm by Atkinson and Donev # # Arguments: # Fx ... n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters # N ... the required number of trials of the experiment # Phi.app ... the optimal value of the correspondng approximate (relaxed) problem # If Phi.app==NULL, the value is computed using od.AA # crit ... optimality criterion to be used # Possible choices are "D", "A", "I" # eta ... the summation weights for the criterion of I-optimality # if crit is not I then eta is ignored # if crit==I and eta=NULL, the program sets eta=(1,...,1) # w1 ... initial design; it must have a non-singular information matrix # and the size sum(w1) of the design must be N # if w1==NULL, a non-singular initial design is generated by od.PINs # K,L ... integer numbers (or Inf) representing parameters of the method # Various combinations of K and L lead to specific variants # If K=NULL or L=NULL, reasonable values are heuristically chosen # rest.max ... maximum allowed number of restarts of the method # t.max ... threshold for the maximum computation time # # Output is the list with components: # call ... the call of the function # w.best ... the resulting exact design # supp ... the indices of the support of w.best # w.supp ... the weights of w.best on the support # M.best ... the information matrix of w.best # Phi.best ... the criterion value of w.best # eff.best ... a lower bound on the eff of w.best wrt the optimal approximate design # n.rest ... number of restarts performed # n.ex ... number of exchanges performed # t.act ... the actual time of the computation # # Authors of the code: # Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Input check if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be greater than or equal to the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be at least 2. The case of a single parameter is trivial.") if (!is.vector(N) || !is.numeric(N) || !all(is.finite(N)) || length(N)!=1 || N<=0 || N!=round(N)) stop("The required size N of the design must be a natural number.") if (!is.null(Phi.app)) if (!is.vector(Phi.app) || !is.numeric(Phi.app) || length(Phi.app)!=1 || Phi.app<=0) stop("Phi.app must be a positive number.") if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") if (is.null(eta)) eta <- rep(1, n) if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative vector of length nrow(Fx).") if (!is.null(w1)) if(!is.vector(w1) || !is.numeric(w1) || !all(is.finite(w1)) || min(w1)<0 || max(w1)<=0) stop("w1 must be NULL or a non-negative, nonzero finite numeric vector of length nrow(Fx).") if (!is.null(K)) if (!is.vector(K) || !is.numeric(K) || length(K)!=1 || K!=round(K) || K<1) stop("K must be NULL, or a natural number, or Inf.") if (!is.null(L)) if (!is.vector(L) || !is.numeric(L) || length(L)!=1 || L!=round(L) || L<1) stop("L must be NULL, or a natural number, or Inf.") if (!is.vector(rest.max) || !is.numeric(rest.max) || length(rest.max)!=1 || rest.max<=0) stop("rest.max must be a positive number.") if (!is.vector(t.max) || !is.numeric(t.max) || length(t.max)!= 1 || t.max <= 0) stop("t.max must be a positive real number.") # Use the appropriate engine if (crit=="D") res <- od.D.KL(Fx=Fx, N=N, Phi.app=Phi.app, w1=w1, K=K, L=L, rest.max=rest.max, t.max=t.max) if (crit=="A") res <- od.A.KL(Fx=Fx, N=N, Phi.app=Phi.app, w1=w1, K=K, L=L, rest.max=rest.max, t.max=t.max) if (crit=="I") res <- od.A.KL(Fx=Fx.ItoA(Fx, eta), N=N, Phi.app=Phi.app, w1=w1, K=K, L=L, rest.max=rest.max, t.max=t.max) # Output supp <- (1:n)[res$w.best > 0] res <- list(call=cl, w.best=res$w.best, supp=supp, w.supp=res$w.best[supp], M.best=infmat(Fx, res$w.best), Phi.best=res$Phi.best, eff.best=res$eff.best, n.rest=res$n.rest, n.ex=res$n.ex, t.act=res$t.act) return(res) } ############################################################# od.MISOCP <- function(Fx, b1=NULL, A1=NULL, b2=NULL, A2=NULL, b3=NULL, A3=NULL, w0=NULL, bin=FALSE, crit="D", eta=NULL, h=NULL, gap=NULL, t.max=120) { # Purpose: # Computes a (hopefully optimal or at least efficient) exact experimental design. # Uses the MISOCP representation of the problem. # See the paper Sagnol & Harman 2015 for details. # # Authors: # Radoslav Harman & Lenka Filova # # Arguments: # Fx: n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters. # b1, A1: The real vector of length k1 and the k1 times n matrix of real numbers. # If both not NULL, the constraints A1*w<=b1 are added to the set of all design constraints. # b2, A2: The real vector of length k2 and the k2 times n matrix of real numbers. # If both not NULL, the constraints A2*w>=b2 are added to the set of all design constraints. # b3, A3: The real vector of length k3 and the k3 times n matrix of real numbers. # If both not NULL, the constraints A3*w==b3 are added to the set of all design constraints. # Note 1: At least one of b1, b2, b3 must be non-NULL. # Note 2: If bi is non-NULL and Ai is NULL for some i then Ai is set to be the vector of ones. # Note 3: if bi is NULL for some i then Ai is ignored. # w0: The non-negative real vector of length n representing the design to be augmented. # If NULL, w0 is set to the vector of zeros. # If not NULL, the constraints w>=w0 are added to the set of all design constraints. # bin: Should the design be binary, i.e., without repeated use of design points? # crit: The optimality criterion to be used. Possible choices are "D", "A", "I", "c" # eta: The summation weights for the criterion of I-optimality # If crit is not "I" then eta is ignored. # If crit is "I" and eta is NULL, eta is set to be the vector of ones. # h: The vector determining the linear combination of parameters if c-optimality is used # If crit is not "c" then h is ignored. # If crit is "c" and h is NULL, h is set to be the vector c(0,...,0,1). # gap: The gap for the MISOCP solver to stop the computation. # If NULL, the default gap is used. # Note: Setting gap=0 and t.max=Inf will provide the optimal exact design, # but the computation may be extremely long. # t.max: threshold for the maximum computation time (gurobi only) # # Output is the list with components: # call: the call of the function # w.best: the resulting approximate design # supp: the indices of the support of w.best # w.supp: the weights of w.best on the support # M.best: the information matrix of w.best # Phi.best: the criterion value of w.best # status: the status of the optimization solver # See http://www.gurobi.com/documentation/8.1/refman/optimization_status_codes.html # t.act: the actual time of the computation cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) ## Verify Fx if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be >= the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be >= 2. Use od.m1 for models with a one-dimensional parameter.") ## Verify b1, b2, b3 if (is.null(b1) & is.null(b2) & is.null(b3)) stop("At least one of b1, b2, b3 must be non-NULL.") k1 <- k2 <- k3 <- 0 ## Verify b1 if (!is.null(b1)) { if (!is.vector(b1) || !is.numeric(b1) || !all(is.finite(b1))) stop("b1 must be either NULL, or a real number, or a vector of real numbers.") k1 <- length(b1) if (is.null(A1)) { if (length(b1)!=1) stop("If b1 is non-NULL and A1 is NULL then b1 must be a real number.") A1 <- matrix(1, nrow=k1, ncol=n) } } ## Verify A1 if (!is.null(A1)) { if (!is.matrix(A1) || !is.numeric(A1) || !all(is.finite(A1))) stop("A1 must be NULL, or a matrix of real numbers.") if (nrow(A1)!=length(b1)) stop("The dimension nrow(A1) must be equal to the number length(b1) of constraints.") if (ncol(A1)!=n) stop("The dimension ncol(A1) must be equal to the number nrow(Fx) of design points.") } ## Verify b2 if (!is.null(b2)) { if (!is.vector(b2) || !is.numeric(b2) || !all(is.finite(b2))) stop("b2 must be either NULL, or a real number, or a vector of real numbers.") k2 <- length(b2) if (is.null(A2)) { if (length(b2)!=1) stop("If b2 is non-NULL and A1 is NULL then b1 must be a real number.") A2 <- matrix(1, nrow=k2, ncol=n) } } ## Verify A2 if (!is.null(A2)) { if (!is.matrix(A2) || !is.numeric(A2) || !all(is.finite(A2))) stop("A2 must be NULL, or a matrix of real numbers.") if (nrow(A2)!=length(b2)) stop("The dimension nrow(A2) must be equal to the number length(b2) of constraints.") if (ncol(A2)!=n) stop("The dimension ncol(A2) must be equal to the number nrow(Fx) of design points.") } ## Verify b3 if (!is.null(b3)) { if (!is.vector(b3) || !is.numeric(b3) || !all(is.finite(b3))) stop("b3 must be either NULL, or a real number, or a vector of real numbers.") k3 <- length(b3) if (is.null(A3)) { if (length(b3)!=1) stop("If b3 is non-NULL and A1 is NULL then b3 must be a real number.") A3 <- matrix(1, nrow=k3, ncol=n) } } ## Verify A3 if (!is.null(A3)) { if (!is.matrix(A3) || !is.numeric(A3) || !all(is.finite(A3))) stop("A3 must be NULL, or a matrix of real numbers.") if (nrow(A3)!=length(b3)) stop("The dimension nrow(A3) must be equal to the number length(b3) of constraints.") if (ncol(A3)!=n) stop("The dimension ncol(A3) must be equal to the number nrow(Fx) of design points.") } ## Verify w0 if (is.null(w0)) w0 <- rep(0, n) if (!is.vector(w0) || !is.numeric(w0) || !all(is.finite(w0)) || (min(w0) < 0)) stop("w0 must be NULL, or a vector of non-negative real numbers.") if (length(w0) != n) stop("The length of w0 must be equal to the number nrow(Fx) of design points.") ## Verify bin if (!is.vector(bin) || !is.logical(bin) || length(bin)>1) stop("bin must be a single logical value (TRUE or FALSE).") ## Verify crit if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") ## Verify eta if (!is.null(eta)) { if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative vector of length nrow(Fx).") } ## Verify h if (!is.null(h)) { if (!is.vector(h) || !is.numeric(h) || !all(is.finite(h)) || !(length(h)==m)) stop("h must be NULL or a vector of real numbers with the same length as ncol(Fx)") } else { h <- c(rep(0, m-1), 1) } ## Verify gap if (!is.null(gap)) if (!is.vector(gap) || !is.numeric(gap) || !(length(h)==1) || (gap < 0)) stop("gap must be a nonnegative real number or Inf.") ## Verify t.max if (!is.vector(t.max) || !is.numeric(t.max) || (length(t.max) != 1) || (t.max <= 0)) stop("t.max must be a positive real number or Inf.") if (!is.null(w0) && sum(w0)>0.5) { A3 <- rbind(A3, diag(n)[w0>0.5,]) b3 <- c(b3, w0[w0>0.5]) } # Print info info <- paste("Running od.MISOCP for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, ", m=", m, ",", sep="") info <- paste(info, ", k=", k1+k2+k3, ".", sep="") print(info, quote=FALSE) info <- paste("Other parameters: ", sum(b1), ",", sum(b2), ",", sum(b3), ";", sum(w0), ",", bin, ",", crit, ",", sum(eta), ",", sum(h), ",", gap, sep="") print(info, quote=FALSE) # Track time start <- as.numeric(proc.time()[3]) if (crit=="D") res <- od.D.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, type="exact", gap=gap, t.max=t.max) if (crit=="A") res <- od.A.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, type="exact", gap=gap, t.max=t.max) if (crit=="I") res <- od.A.SOCP(Fx.ItoA(Fx, eta), b1, A1, b2, A2, b3, A3, bin=bin, type="exact", gap=gap, t.max=t.max) if (crit=="c") res <- od.c.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, h=h, type="exact", gap=gap, t.max=t.max) w <- res$w.best if (!is.null(w)) { supp <- (1:n)[w > 0.5]; w.supp <- w[supp] M.best <- infmat(Fx, w) Phi.best <- optcrit(Fx, w, crit=crit) } else { supp <- NULL; w.supp <- NULL M.best <- NULL; Phi.best <- 0 } t.act <- round(as.numeric(proc.time()[3]) - start, 2) return(list(call=cl, w.best=w, supp=supp, w.supp=w.supp, M.best=M.best, Phi.best=Phi.best, status=res$status, t.act=t.act)) } ################################################################### od.PIN <- function (Fx, N=NULL, crit="D", eta=NULL, rnd=0) { # Purpose: # Produces a simple and reasonably efficient exact design w that # satisfies the size constraint sum(w)<=N, where N<=m (=ncol(Fx)) # It can be used to initiate various optimal design algorithms, # exact as well as approximate (with uniform weights) # # Details will be published by Radoslav Harman and Samuel Rosa # # Arguments: # Fx ... n times m matrix of model regressors, where n is the number # of design points, m is the number of model parameters # N ... the requires size of the design not exceeding m # if NULL then N is set to be m (the usual choice) # crit ... the optimality criterion ("D", "A" or "I") # Note that I is immediately transformed to A and treated as such # eta ... vector of weights for the criterion of I optimality # eta is ingored if other criterion than I is selected # rnd ... degree of randomization between 0 and 1 # rnd == 0 leads to non-randomized greedy selection of trials # rnd == 1 leads to completely random selection of trials # rnd == 0.05 is usually enough for a good variety on initial designs # # Output is a list of: # w ... the suggested "reasonable" design # supp ... the indices of the support of w # cl <- match.call() # print("Call of the function:", quote=FALSE) # print(cl, quote=FALSE) eps <- 1e-9 # numerical tolerance # Check the input if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); one.n <- rep(1, n) m <- ncol(Fx); one.m <- rep(1, m) if (m < 2) stop("The number m=ncol(Fx) of parameters must be at least 2. The case of a single parameter is trivial.") if (is.null(N)) N <- m if (!is.vector(N) || !is.numeric(N) || !all(is.finite(N)) || (length(N) != 1) || (N <= 0) || (N != round(N))) stop("The required size N of the design must be a natural number.") if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") if (is.null(eta)) eta <- rep(1, n) if (!is.vector(eta) || !is.numeric(eta) || !all(is.finite(eta)) || (length(eta) != n) || (min(eta) < eps) ) stop("eta must be NULL, or a vector of n=nrow(Fx) positive numbers.") if (crit == "I") { Fx <- F.ItoA(Fx, eta) crit <- "A" } w <- rep(0, n); supp <- c() M <- matrix(0, ncol=m, nrow=m) for(i in 1:m) { Mp <- ginv(M) Fp <- Fx %*% (diag(m) - M %*% Mp) v <- (Fp^2) %*% one.m if(crit == "A") v <- v / (one.n + ((Fp %*% Mp) * Fp) %*% one.m) i <- which.max(v); if(v[i] < eps) break if(rnd > eps) { i <- sample(1:n, 1, prob=abs(v / v[i])^(1/rnd - 1)) } w[i] <- w[i] + 1 supp <- c(supp, i) M <- M + tcrossprod(Fx[i,]) } return(list(w=w, supp=supp)) } ############################################################### od.SOCP <- function(Fx, b1=NULL, A1=NULL, b2=NULL, A2=NULL, b3=NULL, A3=NULL, w0=NULL, bin=FALSE, crit="D", eta=NULL, h=NULL, t.max=120) { # Note: This is almost the same as od.MISOCP, but without restricting the integer feasible space. # We decided to keep od.MISOCP and od.SOCP separate in accord with the logic of the package # which computes approximate and exact designs via separate functions. # Compared to od.MISOCP: # od.SOCP ommitts the gap parameter (it always uses the default gurobi value), # the bin parameter means restriction of the design weights to the interval [0,1]. # # Purpose: # Computes an optimal approximate experimental design. # Uses the SOCP representation of the problem. # See the paper Sagnol & Harman 2015 for details. # # Authors: # Radoslav Harman & Lenka Filova # # Arguments: # Fx: n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters. # b1, A1: The real vector of length k1 and the k1 times n matrix of real numbers. # If both not NULL, the constraints A1*w<=b1 are added to the set of all design constraints. # b2, A2: The real vector of length k2 and the k2 times n matrix of real numbers. # If both not NULL, the constraints A2*w>=b2 are added to the set of all design constraints. # b3, A3: The real vector of length k3 and the k3 times n matrix of real numbers. # If both not NULL, the constraints A3*w==b3 are added to the set of all design constraints. # Note 1: At least one of b1, b2, b3 must be non-NULL. # Note 2: If bi is non-NULL and Ai is NULL for some i then Ai is set to be the vector of ones. # Note 3: if bi is NULL for some i then Ai is ignored. # w0: The non-negative real vector of length n representing the design to be augmented. # If NULL, w0 is set to the vector of zeros. # If not NULL, the constraints w>=w0 are added to the set of all design constraints. # bin: Should the design be binary, i.e., without repeated use of design points? # crit: The optimality criterion to be used. Possible choices are "D", "A", "I", "c" # eta: The summation weights for the criterion of I-optimality # If crit is not "I" then eta is ignored. # If crit is "I" and eta is NULL, eta is set to be the vector of ones. # h: The vector determining the linear combination of parameters if c-optimality is used # If crit is not "c" then h is ignored. # If crit is "c" and h is NULL, h is set to be the vector c(0,...,0,1). # t.max: threshold for the maximum computation time (gurobi only) # # Output is the list with components: # call: the call of the function # w.best: the resulting approximate design # supp: the indices of the support of w.best # w.supp: the weights of w.best on the support # M.best: the information matrix of w.best # Phi.best: the criterion value of w.best # status: the status of the optimization solver # See http://www.gurobi.com/documentation/8.1/refman/optimization_status_codes.html # t.act: the actual time of the computation cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) ## Verify Fx if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (n < m) stop("The number nrow(Fx) of design points must be >= the number ncol(Fx) of parameters of the model.") if (m < 2) stop("The number ncol(Fx) of parameters must be >= 2. Use od.m1 for models with a one-dimensional parameter.") ## Verify b1, b2, b3 if (is.null(b1) & is.null(b2) & is.null(b3)) stop("At least one of b1, b2, b3 must be non-NULL.") k1 <- k2 <- k3 <- 0 ## Verify b1 if (!is.null(b1)) { if (!is.vector(b1) || !is.numeric(b1) || !all(is.finite(b1))) stop("b1 must be either NULL, or a real number, or a vector of real numbers.") k1 <- length(b1) if (is.null(A1)) { if (length(b1)!=1) stop("If b1 is non-NULL and A1 is NULL then b1 must be a real number.") A1 <- matrix(1, nrow=k1, ncol=n) } } ## Verify A1 if (!is.null(A1)) { if (!is.matrix(A1) || !is.numeric(A1) || !all(is.finite(A1))) stop("A1 must be NULL, or a matrix of real numbers.") if (nrow(A1)!=length(b1)) stop("The dimension nrow(A1) must be equal to the number length(b1) of constraints.") if (ncol(A1)!=n) stop("The dimension ncol(A1) must be equal to the number nrow(Fx) of design points.") } ## Verify b2 if (!is.null(b2)) { if (!is.vector(b2) || !is.numeric(b2) || !all(is.finite(b2))) stop("b2 must be either NULL, or a real number, or a vector of real numbers.") k2 <- length(b2) if (is.null(A2)) { if (length(b2)!=1) stop("If b2 is non-NULL and A1 is NULL then b1 must be a real number.") A2 <- matrix(1, nrow=k2, ncol=n) } } ## Verify A2 if (!is.null(A2)) { if (!is.matrix(A2) || !is.numeric(A2) || !all(is.finite(A2))) stop("A2 must be NULL, or a matrix of real numbers.") if (nrow(A2)!=length(b2)) stop("The dimension nrow(A2) must be equal to the number length(b2) of constraints.") if (ncol(A2)!=n) stop("The dimension ncol(A2) must be equal to the number nrow(Fx) of design points.") } ## Verify b3 if (!is.null(b3)) { if (!is.vector(b3) || !is.numeric(b3) || !all(is.finite(b3))) stop("b3 must be either NULL, or a real number, or a vector of real numbers.") k3 <- length(b3) if (is.null(A3)) { if (length(b3)!=1) stop("If b3 is non-NULL and A1 is NULL then b3 must be a real number.") A3 <- matrix(1, nrow=k3, ncol=n) } } ## Verify A3 if (!is.null(A3)) { if (!is.matrix(A3) || !is.numeric(A3) || !all(is.finite(A3))) stop("A3 must be NULL, or a matrix of real numbers.") if (nrow(A3)!=length(b3)) stop("The dimension nrow(A3) must be equal to the number length(b3) of constraints.") if (ncol(A3)!=n) stop("The dimension ncol(A3) must be equal to the number nrow(Fx) of design points.") } ## Verify w0 if (is.null(w0)) w0 <- rep(0, n) if (!is.vector(w0) || !is.numeric(w0) || !all(is.finite(w0)) || (min(w0) < 0)) stop("w0 must be NULL, or a vector of non-negative real numbers.") if (length(w0) != n) stop("The length of w0 must be equal to the number nrow(Fx) of design points.") ## Verify bin if (!is.vector(bin) || !is.logical(bin) || length(bin)>1) stop("bin must be a single logical value (TRUE or FALSE).") ## Verify crit if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") ## Verify eta if (!is.null(eta)) { if (!is.vector(eta) || !is.numeric(eta) || length(eta)!=n || min(eta)<0) stop("eta must be NULL, or a non-negative vector of length nrow(Fx).") } ## Verify h if (!is.null(h)) { if (!is.vector(h) || !is.numeric(h) || !all(is.finite(h)) || !(length(h)==m)) stop("h must be NULL or a vector of real numbers with the same length as ncol(Fx)") } else { h <- c(rep(0, m-1), 1) } ## Verify t.max if (!is.vector(t.max) || !is.numeric(t.max) || (length(t.max) != 1) || (t.max <= 0)) stop("t.max must be a positive real number or Inf.") if (!is.null(w0) && sum(w0)>0.5) { A3 <- rbind(A3, diag(n)[w0>0.5,]) b3 <- c(b3, w0[w0>0.5]) } # Print info info <- paste("Running od.SOCP for cca", t.max, "seconds") info <- paste(info, " starting at ", Sys.time(), ".", sep = "") print(info, quote = FALSE) info <- paste("The problem size is n=", n, sep="") info <- paste(info, ", m=", m, ",", sep="") info <- paste(info, ", k=", k1+k2+k3, ".", sep="") print(info, quote=FALSE) info <- paste("Other parameters: ", sum(b1), ",", sum(b2), ",", sum(b3), ";", sum(w0), ",", bin, ",", crit, ",", sum(eta), ",", sum(h), sep="") print(info, quote=FALSE) # Track time start <- as.numeric(proc.time()[3]) if (crit=="D") res <- od.D.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, type="approximate", gap=NULL, t.max=t.max) if (crit=="A") res <- od.A.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, type="approximate", gap=NULL, t.max=t.max) if (crit=="I") res <- od.A.SOCP(Fx.ItoA(Fx, eta), b1, A1, b2, A2, b3, A3, bin=bin, type="approximate", gap=NULL, t.max=t.max) if (crit=="c") res <- od.c.SOCP(Fx, b1, A1, b2, A2, b3, A3, bin=bin, h=h, type="approximate", gap=NULL, t.max=t.max) w <- res$w.best if (!is.null(w)) { supp <- (1:n)[w > 1e-6]; w.supp <- w[supp] M.best <- infmat(Fx, w) Phi.best <- optcrit(Fx, w, crit=crit) } else { supp <- NULL; w.supp <- NULL M.best <- NULL; Phi.best <- 0 } t.act <- round(as.numeric(proc.time()[3]) - start, 2) return(list(call=cl, w.best=w, supp=supp, w.supp=w.supp, M.best=M.best, Phi.best=Phi.best, status=res$status, t.act=t.act)) } ############################################################ optcrit <- function(Fx, w, crit="D", eta=NULL, h=NULL, tol=1e-12) { # Purpose: # Computes the criterion value of design w in model determined by Fx # This procedure is not supposed to be repeatedly run in a cycle # # Arguments: # F ... n times m matrix of all regressors # w ... design (approximate or exact) # crit ... criterion; possible values are "D", "A", "I" and "c" # if crit is "c" then library MASS is needed # eta ... weights for the criterion of I-optimality # if crit is not "I" then eta is ignored # if eta=NULL, uniform measure over the design region is assumed # h ... non-zero vector of coefficients for c-optimality # if crit is not "c" then h is ignored # if h=NULL, then h=(0,...,0,1)^T is assumed # # Value: # Non-negative number corresponding to the criterion value # # Authors: # Radoslav Harman and Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) # Check input if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("F must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx) if (!is.vector(w) || !is.numeric(w) || !all(is.finite(w)) || (min(w) < 0)) stop("w must be a vector of non-negative real numbers.") if (length(w) != n) stop("The length of w must be equal to the number dim(F)[1] of design points.") if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I", "c"))) stop("crit must be 'D', 'A', 'I' or 'c'.") if (crit=="I") { if (is.null(eta)) eta <- rep(1, n) if (!is.vector(eta) || !is.numeric(eta) || (length(eta) != n) || !(min(eta) > 0)) stop("eta must be a positive real vector of length nrow(Fx).") } if (crit=="c") { if (is.null(h)) h <- c(rep(0, m-1), 1) if (!is.vector(h) || !is.numeric(h) || (length(h) != m) || (sum(abs(h)) < tol)) stop("h must be a nonzero real vector of length ncol(Fx).") } if (crit=="I") { M <- infmat(Fx.ItoA(Fx, eta), w) } else { M <- infmat(Fx, w) } if ((min(abs(eigen(M)$values)) <= m*tol)) { if(!(crit=="c")) { return(0) } else { M.plus <- ginv(M) if(max(abs(M %*% M.plus %*% h - h)) > tol) return(0) return(1 / (t(h) %*% M.plus %*% h)[1, 1]) } } else { if (crit=="A" || crit=="I") return(m * sum(diag(solve(M)))^(-1)) if (crit=="D") return((abs(det(M)))^(1/m)) if (crit=="c") return(1 / (t(h) %*% solve(M) %*% h)[1, 1]) } } ######################################################## varfun <- function(Fx, w, crit="D", eta=NULL) { # Purpose: # Computes the 'variance function' of design w in the model determined by Fx # # Arguments: # Fx ... n times m(<=n) matrix containing all possible regressors (as rows), # n is the number of design points, m(>=2) is the number of parameters # w ... non-negative vector of length n representing the design # Note that the information matrix of w must be non-singular # crit ... the criterion (either "D" or "A" or "I") # # Output: # var.fun ... the vector of length nrow(Fx) giving the variance function # # Authors: # Radoslav Harman, Lenka Filova cl <- match.call() print("Call of the function:", quote=FALSE) print(cl, quote=FALSE) eps <- 1e-12 # Numerical tolerance # Check the input if (!is.matrix(Fx) || !is.numeric(Fx) || !all(is.finite(Fx))) stop("Fx must be a matrix of real numbers.") n <- nrow(Fx); m <- ncol(Fx); one <- rep(1, m) if (n < m) stop("The number n=nrow(Fx) of design points must be greater or equal to the number m=ncol(Fx) of parameters.") if (m < 2) stop("The number m=ncol(Fx) of parameters must be at least 2. Use the procedure od.m1 for models with a one-dimensional parameter.") if (!is.vector(w) || !is.numeric(w) || !all(is.finite(w)) || (min(w) < 0)) stop("w must be a vector of non-negative real numbers.") if (length(w) != n) stop("The length of w must be equal to the number nrow(Fx) of design points.") if (!is.vector(crit) || !is.character(crit) || (length(crit) != 1) || !is.element(crit, c("D", "A", "I"))) stop("crit must be 'D', 'A', or 'I'.") if (crit=="I") { if (is.null(eta)) eta <- 1:n if (!is.vector(eta) || !is.numeric(eta) || (length(eta) != n) || !(min(eta) > 0)) stop("eta must be a positive real vector of length nrow(Fx).") } if (crit == "I") Fx <- Fx.ItoA(Fx, eta) supp <- w > 0; w.supp <- w[supp]; Fx.supp <- Fx[supp,] M <- crossprod(sqrt(w.supp) * Fx.supp) if(rcond(M) < eps) stop("Information matrix is singular, i.e., it is impossible to compute the variance function.") if(crit=="D") { var.fun <- ((Fx %*% t(chol(solve(M))))^2) %*% rep(1, m) } else { var.fun <- ((Fx %*% solve(M))^2) %*% rep(1, m) } return(as.vector(var.fun)) }