################################################### ### Preliminaries ################################################### suppressWarnings(RNGversion("3.5")) set.seed(1782) inst_pkgs <- installed.packages() general_pkgs <- c("xtable", "jsonlite", "slam") roi_pkgs <- c("ROI", "ROI.plugin.alabama", "ROI.plugin.deoptim", "ROI.plugin.ecos", "ROI.plugin.glpk", "ROI.plugin.lpsolve", "ROI.plugin.msbinlp", "ROI.plugin.neos", "ROI.plugin.qpoases", "ROI.plugin.scs", "ROI.plugin.symphony", "ROI.models.netlib", "ROI.models.miplib") missing_pkgs <- setdiff(c(general_pkgs, roi_pkgs), inst_pkgs[, "Package"]) msg <- "" if (length(missing_pkgs)) { msg <- paste("Some packages are missing! Use the code below to install", "the missing packages:\n\n", collapse = " ") msg <- sprintf("%sinstall.packages(%s)\n", msg, paste(deparse(missing_pkgs), collapse = "\n")) } if (!("ROI.tests" %in% inst_pkgs[, "Package"])) { if (!nchar(msg)) { msg <- paste("Some packages are missing! Use the code below to install", "the missing packages:\n\n", collapse = " ") } msg <- sprintf("%s%s\n", msg, "install.packages(\"ROI.tests\", repos = \"http://R-Forge.R-project.org\")") } if (nchar(msg)) stop(msg) library("xtable") ################################################### ### 4. A general optimization infrastructure for R ################################################### library("ROI") ################################################### ### 4.4. Bounds ################################################### V_bound(li = 1:4, ui = 1:4, lb = c(-Inf, 0, 2, 0), ub = c(4, 100, Inf, Inf)) V_bound(lb = c(-Inf, 0, 2, 0), ub = c(4, 100, Inf, Inf)) V_bound(li = c(1L, 3L), ui = c(1L, 2L), lb = c(-Inf, 2), ub = c(4, 100), nobj = 4L) V_bound(li = 3, lb = 0, ld = -20, ud = 20, nobj = 4L) ################################################### ### 4.5. Optimization problem ################################################### A <- rbind(c(5, 7, 2), c(3, 2, -9), c(1, 3, 1)) dir <- c("<=", "<=", "<=") rhs <- c(61, 35, 31) lp <- OP(objective = L_objective(c(3, 7, -12)), constraints = L_constraint(A, dir = dir, rhs = rhs), bounds = V_bound(li = 3, ui = 3, lb = -10, ub = 10, nobj = 3), maximum = TRUE) lp <- OP() objective(lp) <- L_objective(c(3, 7, -12)) constraints(lp) <- L_constraint(A, dir = c("<=", "<=", "<="), rhs = rhs) bounds(lp) <- V_bound(li = 3, ui = 3, lb = -10, ub = 10, nobj = 3) maximum(lp) <- TRUE lp param <- rep.int(1, length(objective(lp))) objective(lp)(param) terms(objective(lp)) ################################################### ### 5. Package ROI ################################################### ### 5.1. Solving optimization problems ################################################### (lp_sol <- ROI_solve(lp, solver = "glpk")) ################################################### ### 5.2. Solution and status code ################################################### solution(lp_sol) blp <- OP(objective = L_objective(c(-1, -1, -1, -1, -99)), constraints = L_constraint(L = rbind(c(1, 1, 0, 0, 0), c(0, 0, 1, 1, 0), c(0, 0, 0, 1, 1)), dir = c("<=", "<=", "<="), rhs = rep.int(1, 3)), types = rep("B", 5L)) (blp_sol <- ROI_solve(blp, solver = "msbinlp", method = "glpk", nsol_max = 32)) solution(blp_sol) lp_inaccurate_sol <- ROI_solve(lp, solver = "scs", tol = 1e-32) solution(lp_inaccurate_sol) solution(lp_inaccurate_sol, force = TRUE) solution(lp_sol, type = "dual") solution(lp_sol, type = "aux") solution(lp_sol, type = "msg") solution(lp_sol, type = "objval") solution(lp_sol, type = "status") solution(lp_sol, type = "status_code") ################################################### ### 5.3. Reformulations ################################################### Q <- rbind(c(0, 3, 0), c(0, 0, 1), c(0, 0, 0)) bqp <- OP(Q_objective(Q = Q + t(Q), L = c(-1, -4, -1)), types = rep("B", 3)) glpk_signature <- ROI_solver_signature("glpk") head(glpk_signature, 3) milp <- ROI_reformulate(x = bqp, to = glpk_signature) ROI_solve(milp, solver = "glpk") ################################################### ### 5.4. ROI solvers ################################################### ROI_available_solvers(bqp)[, c("Package", "Repository")] head(ROI_registered_solvers(), 3) ROI_applicable_solvers(lp) head(ROI_installed_solvers(), 3) ################################################### ### 5.5. ROI read/write ################################################### lp_file <- tempfile() ROI_write(lp, lp_file, "lp_lpsolve") writeLines(readLines(lp_file)) ROI_read(lp_file, "lp_lpsolve") ################################################### ### 5.6. ROI models ################################################### library("ROI.models.miplib") if (length(miplib()) == 0L) { miplib_download_benchmark(quiet = TRUE) miplib_download_metainfo() } ops <- miplib("ns1766074") ops library("ROI.models.netlib") agg <- netlib("agg") ROI_solve(agg, "glpk") library("jsonlite") nested_unclass <- function(x) { x <- unclass(x) if (is.list(x)) x <- lapply(x, nested_unclass) x } agg_json <- toJSON(nested_unclass(agg)) tmp_file <- tempfile() writeLines(agg_json, tmp_file) ################################################### ### 5.7. ROI settings ################################################### simple_gradient <- function(func, x, ...) { numDeriv::grad(func, x, method = "simple", ...) } ROI_options("gradient", simple_gradient) simple_jacobian <- function(func, x, ...) { numDeriv::jacobian(func, x, method = "simple", ...) } ROI_options("jacobian", simple_jacobian) ROI_options("default_solver") ROI_options("default_solver", "glpk") ROI_options("default_solver", "auto") Sys.setenv(ROI_LOAD_PLUGINS = FALSE) ################################################### ### 6. Examples ################################################### ### 6.1. Linear optimization problems ################################################### lp <- OP(c(1, 2), maximum = TRUE, L_constraint(L = rbind(c(1, 8), c(5, 1)), dir = c("==", "<="), rhs = c(9, 6)), bounds = V_bound(lb = c(-9, -7), ub = c(9, 7))) (lp_sol <- ROI_solve(lp, "glpk")) solution(lp_sol) ################################################### ### 6.2. Quadratic optimization problems ################################################### qp <- OP(Q_objective(Q = diag(2), L = c(-1, 0)), L_constraint(c(4, 6), ">=", 10)) (qp_sol <- ROI_solve(qp, "qpoases")) solution(qp_sol) qcqp <- qp constraints(qcqp) <- rbind(constraints(qp), Q_constraint(-diag(c(2, 2)), L = c(7, 11), ">=", 40)) ROI_applicable_solvers(qcqp) (qcqp_sol <- ROI_solve(qcqp, "alabama", start = c(5, 5))) solution(qcqp_sol) ## This part is commented out to avoid to much traffic ## on the neos server. # (qcqp_sol <- ROI_solve(qcqp, "neos", method = "mosek")) # solution(qcqp_sol) ################################################### ### 6.3. Conic optimization problems ################################################### cpeq <- C_constraint(c(1, 8), K_zero(1), 9) cpleq <- C_constraint(c(5, 1), K_lin(1), 6) zlcp <- lp constraints(zlcp) <- c(cpeq, cpleq) (zlcp_sol <- ROI_solve(zlcp, solver = "ecos")) solution(zlcp_sol) soc1 <- OP(c(1, 1, 0), C_constraint(L = rbind(c(0, 0, -7), c(-3, 0, 0), c(0, -5, 0)), cone = K_soc(3), rhs = c(6, 2, 4)), maximum = TRUE, bounds = V_bound(ld = -Inf, ui = 3, ub = 9, nobj = 3)) (soc1_sol <- ROI_solve(soc1, solver = "ecos")) solution(soc1_sol) A <- rbind(c(0, 0, -1), c(-1, 0, 0), c(0, -1, 0)) soc2 <- OP(objective = L_objective(c(0, 0, 1)), constraints = c(C_constraint(A, K_soc(3), c(0, 0, 0)), L_constraint(c(1, 1, 0), "==", 2))) (soc2_sol <- ROI_solve(soc2, solver = "ecos")) solution(soc2_sol) cexpp <- C_constraint(L = rbind(c(-3, -5, 0, 0), c(0, 0, 0, 0), c(0, 0, -11, -12)), cone = K_expp(1), rhs = c(7, 1, 9)) expp1 <- OP(c(1, 2, 0, 0), cexpp, bounds = V_bound(ld = -Inf, ub = c(20, 20, 50, 50)), maximum = TRUE) (expp1_sol <- ROI_solve(expp1, solver = "ecos")) solution(expp1_sol) A <- rbind(c(0, -1), c(0, 0), c(-7, 0)) log1 <- OP(L_objective(c(0, 1), c("x", "t")), C_constraint(A, K_expp(1), rhs = c(0, 1, 9)), bounds = V_bound(lb = c(0, -Inf), ub = c(1, Inf)), maximum = TRUE) (log1_sol <- ROI_solve(log1, solver = "ecos")) solution(log1_sol) A <- rbind(c(-1, 0), c(0, 0), c(0, -1)) cpowp <- C_constraint(A, K_powp(1/4), rhs = c(5, 1, 2)) powp1 <- OP(c(3, 5), cpowp, bounds = V_bound(lb = c(0, 2))) (powp1_sol <- ROI_solve(powp1, solver = "scs", max_iter = 1e6)) solution(powp1_sol) (A <- matrix(c(1, 2, 3, 2, 4, 5, 3, 5, 6), nrow = 3)) vech(A) F1 <- rbind(c(10, 3), c(3, 10)) F2 <- rbind(c(6, -4), c(-4, 10)) F3 <- rbind(c(8, 1), c(1, 6)) F0 <- rbind(c(16, -13), c(-13, 60)) psd <- OP(objective = L_objective(c(1, 1, -1)), constraints = C_constraint(L = vech(F1, F2, F3), cone = K_psd(3), rhs = vech(F0))) (psd_sol <- ROI_solve(psd, solver = "scs")) solution(powp1_sol) ################################################### ### 6.4. General nonlinear optimization problems ################################################### nlp_1 <- OP(maximum = TRUE, bounds = V_bound(ud = 42, nobj = 3L)) objective(nlp_1) <- F_objective(F = function(x) prod(x), n = 3, G = function(x) c(prod(x[-1]), prod(x[-2]), prod(x[-3]))) constraint <- function(x) x[1] + 2 * x[2] + 2 * x[3] constraints(nlp_1) <- F_constraint(F = constraint, dir = "<=", rhs = 72, J = function(x) c(1, 2, 2)) nlp_1 nlp_2 <- nlp_1 constraints(nlp_2) <- L_constraint(L = c(1, 2, 3), "<=", 72) nlp_2 (nlp_1_sol_1 <- ROI_solve(nlp_1, "alabama", start = c(10, 10, 10))) solution(nlp_1_sol_1) (nlp_1_sol_2 <- ROI_solve(nlp_1, "alabama", start = c(20, 20, 20))) solution(nlp_1_sol_2) (nlp_1_sol_3 <- ROI_solve(nlp_1, "alabama", start = c(10, 10, 10), tol = 1E-24)) solution(nlp_1_sol_3, force = TRUE) (nlp_1_sol_4 <- ROI_solve(nlp_1, "deoptimr", start = c(20, 20, 20), max_iter = 400, tol = 1E-6)) solution(nlp_1_sol_4) ################################################### ### 6.5. Mixed integer problems ################################################### A <- rbind(c(5, 3), c(7, 1)) milp <- OP(c(5, 7), constraints = L_constraint(A, c(">=", ">="), c(7, 9)), types = c("B", "I")) (milp_sol <- ROI_solve(milp, solver = "glpk")) solution(milp_sol) ################################################### ### 7. Extending ROI ################################################### ### 7.1. Signatures ################################################### OP_signature(lp) glpk_signature <- ROI_plugin_make_signature(objective = "L", constraints = "L", types = c("C", "I", "B", "CI", "CB", "IB", "CIB"), bounds = c("X", "V"), maximum = c(TRUE, FALSE)) ################################################### ### 7.2. Writing a new solver method ################################################### glpk_solve_OP <- function(x, control = list()) { control$canonicalize_status <- FALSE glpk <- list(Rglpk_solve_LP, obj = terms(objective(x))[["L"]], mat = constraints(x)$L, dir = constraints(x)$dir, rhs = constraints(x)$rhs, bounds = bounds(x), types = types(x), max = maximum(x), control = control) mode(glpk) <- "call" if (isTRUE(control$dry_run)) return(glpk) out <- eval(glpk) ROI_plugin_canonicalize_solution(solution = out$solution, optimum = out$optimum, status = out$status, solver = "glpk", message = out) } ################################################### ### 7.3. Registering solver methods ################################################### # ROI_plugin_register_solver_method(glpk_signature, "glpk", glpk_solve_OP) ################################################### ### 7.4. Adding additional information ################################################### # ROI_plugin_add_status_code_to_db(solver = "glpk", code = 5L, # symbol = "GLP_OPT", message = "Solution is optimal.", roi_code = 0L) # ROI_plugin_register_solver_control(solver = "glpk", args = "tm_limit", # roi_control = "max_time") ################################################### ### 7.5. Registering reformulations ###################################################b qp_signature <- ROI_plugin_make_signature(objective = "Q", constraints = c("X", "L"), types = c("B"), bounds = c("X", "V"), cones = c("X"), maximum = c(TRUE, FALSE)) milp_signature <- ROI_plugin_make_signature(objective = "L", constraints = c("X", "L"), types = c("C", "I", "B", "CI", "CB", "IB", "CIB"), bounds = c("X", "V"), maximum = c(TRUE, FALSE), cones = c("X")) # ROI_plugin_register_reformulation( # from = bqp_signature, to = milp_signature, method_name = "bqp_to_lp", # method = bqp_to_lp, description = "", cite = "", author = "") ################################################### ### 7.6. Registering reader/writer ################################################### library("slam") json_reader_lp <- function(file, ...) { stopifnot(is.character(file)) y <- read_json(file, simplifyVector = TRUE) to_slam <- function(x) do.call(simple_triplet_matrix, x) x <- OP() objective(x) <- L_objective(to_slam(y[["objective"]][["L"]]), y[["objective"]][["names"]]) constraints(x) <- L_constraint(to_slam(y[["constraints"]][["L"]]), y[["constraints"]][["dir"]], y[["constraints"]][["rhs"]], y[["constraints"]][["names"]]) types(x) <- y[["types"]] bounds(x) <- structure(y[["bounds"]], class = c("V_bound", "bound")) maximum(x) <- as.logical(y[["maximum"]]) x } json_writer_lp <- function(x, file, ...) { writeLines(toJSON(nested_unclass(x), null = "null"), con = file) } plugin_name <- "io" ROI_plugin_register_writer("json", plugin_name, milp_signature, json_writer_lp) ROI_plugin_register_reader("json", plugin_name, json_reader_lp) fname <- tempfile() file <- ROI_write(lp, file = fname, type = "json") (lp_json <- ROI_read(file = fname, type = "json")) ################################################### ### 7.7. ROI tests ################################################### library("ROI.tests") test_solver("glpk") ################################################### ### 8. Applications ################################################### ### 8.1. L_1 regression ################################################### create_L1_problem <- function(x, y) { p <- ncol(x) + 1L m <- 2 * nrow(x) L <- cbind(1, x, diag(nrow(x)), -diag(nrow(x))) bnds <- V_bound(li = seq_len(p), lb = rep(-Inf, p), nobj = p + m) OP(objective = L_objective(c(rep(0, p), rep(1, m))), constraints = L_constraint(L, dir = eq(nrow(x)), rhs = y), bounds = bnds) } data("stackloss", package = "datasets") l1p <- create_L1_problem(x = as.matrix(stackloss)[,-4], y = stackloss[,4]) L1_res <- ROI_solve(l1p, solver = "glpk") solution(L1_res)[1:ncol(stackloss)] ################################################### ### 8.3. Relative risk regression ################################################### logbin_gps <- function(y, X) { loglikelihood <- function(beta) { xb <- drop(X %*% beta) if (any(xb > 0)) NaN else sum(y * xb + (1 - y) * log(1 - exp(xb))) } gradient <- function(beta) { exb <- exp(drop(X %*% beta)) drop(crossprod(X, (y - exb) / (1 - exb))) } OP(F_objective(loglikelihood, n = ncol(X), G = gradient), L_constraint(L = X, dir = leq(nrow(X)), rhs = double(nrow(X))), bounds = V_bound(ld = -Inf, nobj = ncol(X)), maximum = TRUE) } logbin_cp <- function(y, X, rhs_eps = 1e-7) { y_is_0 <- y == 0L n_y_is_0 <- sum(y_is_0) o <- OP(c(y %*% X, double(n_y_is_0), rep(1, n_y_is_0)), maximum = TRUE) L1 <- cbind(X, matrix(0, nrow(X), 2 * n_y_is_0)) log1exp <- function(xi, j, n_y_is_0) { M <- matrix(0, nrow = 6, ncol = length(xi) + 2 * n_y_is_0) M[1, seq_along(xi)] <- -xi M[3, length(xi) + j] <- -1 M[4, length(xi) + n_y_is_0 + j] <- -1 M[6, length(xi) + j] <- 1 M } L2 <- mapply(log1exp, split(X[y_is_0,], seq_len(n_y_is_0)), seq_len(n_y_is_0), MoreArgs = list(n_y_is_0 = n_y_is_0), SIMPLIFY = FALSE) rhs <- c(c(0, 1, 0), c(0, 1, 1)) rhs <- c(rep(-rhs_eps, nrow(X)), rep(rhs, n_y_is_0)) cones <- c(K_lin(nrow(X)), K_expp(2 * n_y_is_0)) L <- do.call(rbind, c(list(L1), L2)) constraints(o) <- C_constraint(L, cones, rhs) bounds(o) <- V_bound(ld = -Inf, nobj = length(objective(o))) o } generate_data <- function(n) { treat <- factor(sample(c("a", "b", "c"), n, TRUE)) num.diseases <- sample(0:4, n, TRUE) age <- rnorm(n, 50L, 10L) cholesterol <- rnorm(n, 200L, 25L) weight <- rnorm(n, 150L, 20L) sex <- factor(sample(c("female", "male"), n, TRUE)) # Specify population model for log probability that Y = 1 L <- (-1 + 0.1 * (num.diseases - 2) + 0.045 * (age - 70) + (log(cholesterol - 10) - 5.2) - 2 * (treat == "a") + 0.5 * (treat == "b") - 0.5 * (treat == "c")) # Simulate binary y to have Prob(y = 1) = exp(L) y <- as.double(runif(n) < exp(L)) A <- cbind(intercept = 1, age, sex, weight, logchol = log(cholesterol - 10), num.diseases, treatb = (treat == "b"), treatc = (treat == "c")) return(list(y = y, A = A)) } suppressWarnings(RNGversion("3.5")) set.seed(1234) dat <- generate_data(1500L) start <- c(log(0.2), double(ncol(dat$A) - 1)) prob_login_bin_gps <- logbin_gps(dat$y, dat$A) s1 <- ROI_solve(prob_login_bin_gps, "alabama", start = start) solution(s1) prob_login_bin_cp <- logbin_cp(dat$y, dat$A) s2 <- ROI_solve(prob_login_bin_cp, solver = "ecos") head(solution(s2), ncol(dat$A)) obj_fun <- objective(prob_login_bin_gps) obj_fun(head(solution(s2), ncol(dat$A))) - obj_fun(solution(s1))