Nothing
##
## This code is originally from the relations package with minor changes to
## make it compatible with ROI.
##
.expand_types <- function(x, n) {
if ( is.null(x) ) {
## Continuous by default.
rep.int("C", n)
} else {
if ( !is.character(x) || !all(x %in% c("C", "I", "B")) )
stop("Invalid MIP variable types.")
## Be nicer than necessary ...
rep_len(x, n)
}
}
as_roi_control <- function(x) {
x[['order']] <- NULL
x
}
.find_up_to_n_binary_MILP_solutions <-
function(x, nos = 1L, add = FALSE, solver = NULL, control = NULL)
{
## Find up to n binary MILP solutions using a simple branch and cut
## approach (repeatedly splitting the binary variables and cutting
## the non-optimal branches).
objfun <- objective(x)
if(is.na(nos))
nos <- .Machine$integer.max
y <- ROI_solve(x, solver, as_roi_control(control))
if((y$status$code != 0) || (nos == 1L) || !any(x$types == "B"))
return(list(y))
Vopt <- y$objval
tol <- 1e-8
## (Smaller than .Machine$double.eps^0.5 as e.g. used in the
## all.equal() comparisons.)
## We used to have 1e-10, but SYMPHONY was not as precise as this.
v_is_not_optimal <- if(x$maximum) {
function(v) v < Vopt - tol
} else {
function(v) v > Vopt + tol
}
## Improved versions could use relative tolerance, and/or take the
## number of variables into account. Or maybe we can figure out the
## tolerance employed by the solvers when they declare optimality?
if(add) {
## Find up to n solutions by adding binary constraints for each
## split. This is most space efficient, but typically takes
## considerably more time that the default based on successive
## reductions by substitution.
return(.find_up_to_n_binary_MILP_solutions_via_add(x, nos,
solver,
control,
y,
v_is_not_optimal))
}
## Find solutions by recursively splitting binary variable and
## substituting the values into the program.
## Suppose x_j is one of the binary variables. Let x[-j] be the
## vector of variables after dropping x_j. If x_j has value b_j,
## the objective function is
## c[-j]' x[-j] + c_j b_j,
## and the constraints become
## (mat[, -j] %*% x[-j]) dir (rhs - mat[, j] * b_j)
## Note that the new constraint matrix and dir do not depend on the
## value of b_j.
## When recursively splitting, we need to keep track of the b
## values, the sum of the c_j b_j terms to add to the reduced
## objective value to obtain the value of the original problem, and
## the right hand sides.
## Also, binary positions are a nuisance to keep track of, so we
## start by rearranging variables to have the binary variables come
## last in reverse order of splitting.
n_of_variables <- length(x$objective)
types <- .expand_types(x$types, n_of_variables)
binary_positions <- which(types == "B")
n_of_binary_variables <- length(binary_positions)
verbose <- identical(control$verbose, TRUE)
.make_node <- function(b, y, v, r) list(b = b, y = y, v = v, r = r)
.split_single_binary_variable <-
function(node, i, pos_i, obj_i, mat_i, pos_b) {
## Try to avoid unnecessary copying of x via <<-
## manipulations.
## We know that the solution with the current b[i] is
## optimal, so we try the effect of flipping b[i].
b <- node$b
node$y$solution <- node$y$solution[-pos_i]
if(b[i] == 0) {
b1 <- b
b1[i] <- 1
v1 <- node$v + obj_i
r1 <- node$r - mat_i
x$constraints$rhs <<- r1
y1 <- ROI_solve(x, solver, as_roi_control(control))
## Uncomment for debugging ...
## if(verbose) {
## V <- y1$objval + v1
## split <- (y1$status == 0) && !v_is_not_optimal(V)
## message(sprintf("b[i] = 0, flip objval: %f, Delta: %.12f, status: %d, split: %s",
## V, Vopt - V, y1$status, split))
## }
if((y1$status$code != 0) || v_is_not_optimal(y1$objval + v1))
list(node)
else {
## Fill the rest of b with the optimal entries.
b1[-seq_len(i)] <- y1$solution[pos_b]
list(node, .make_node(b1, y1, v1, r1))
}
} else {
b0 <- b
b0[i] <- 0
v0 <- node$v
node$v <- v0 + obj_i
r0 <- node$r
node$r <- r0 - mat_i
x$constraints$rhs <<- r0
y0 <- ROI_solve(x, solver, as_roi_control(control))
## Uncomment for debugging ...
## if(verbose) {
## V <- y0$objval + v0
## split <- (y0$status == 0) && !v_is_not_optimal(V)
## message(sprintf("b[i] = 1, flip objval: %f, Delta: %.12f, status: %d, split: %s",
## V, Vopt - V, y0$status, split))
## }
if((y0$status$code != 0) || v_is_not_optimal(y0$objval + v0))
list(node)
else {
## Fill the rest of b with the optimal entries.
b0[-seq_len(i)] <- y0$solution[pos_b]
list(node, .make_node(b0, y0, v0, r0))
}
}
}
## We allow callers to specify the order in which binary splits
## should be attempted (so the i-th split is for binary variable
## order[i]).
order <- control$order
if(is.null(order) || length(order) != n_of_binary_variables)
order <- seq_len(n_of_binary_variables)
## Rearrange variables to have the binary ones last in reverse split
## order.
ind <- c(which(types != "B"), rev(which(types == "B")[order]))
x$objective <- x$objective[ind]
x$constraints$L <- x$constraints$L[, ind, drop = FALSE]
x$types <- x$types[ind]
y$solution <- y$solution[ind]
pos_i <- n_of_variables
pos_b <- seq.int(from = n_of_variables,
length.out = n_of_binary_variables,
by = -1L)
nodes <- list(.make_node(y$solution[pos_b], y, 0, x$constraints$rhs))
for(i in seq_len(n_of_binary_variables)) {
pos_b <- pos_b[-1L]
obj_i <- x$objective[pos_i]
mat_i <- c(x$constraints$L[, pos_i])
x$objective <- L_objective(x$objective[-pos_i])
x$constraints$L <- x$constraints$L[, -pos_i, drop = FALSE]
x$types <- x$types[-pos_i]
if ( length(x$types) == 0L )
break
nodes <- do.call(c,
lapply(nodes,
.split_single_binary_variable,
i, pos_i, obj_i, mat_i, pos_b))
len <- length(nodes)
if(verbose)
message(gettextf("N_of_binary_variables: %d *** N_of_optimal_branches: %d",
i, len))
if(len >= nos) {
nodes <- nodes[seq_len(nos)]
break
}
pos_i <- pos_i - 1L
}
pos <- order(c(which(types != "B"), which(types == "B")[order]))
finisher <- function(node) {
## Need to reconstruct solutions from the binary and non-binary
## parts and the respective objective values.
y <- node$y
y$solution <- c(y$solution, node$b)[pos]
y$objval <- objfun(y$solution)## y$objval + node$v
y
}
lapply(nodes, finisher)
}
.find_up_to_n_binary_MILP_solutions_via_add <-
function(x, nos = 1L, solver = NULL, control = NULL, y,
v_is_not_optimal)
{
## Recursively add 0/1 constraints for the binary variables.
## Note that one can do this via adding the additional constraints
## to the original ones, or by maintaining them separately (e.g., in
## a vector of values for the binary variables). We do the latter
## which uses as little space as possible, but requires extra time
## for merging both constraints when solving augmented problems.
## Alternatively, one could allow choosing between either approach.
## Note also that we need to keep new constraints and corresponding
## solutions to avoid recomputing solutions when enough were found.
n_of_variables <- length(x$objective)
types <- .expand_types(x$types, n_of_variables)
binary_positions <- which(types == "B")
n_of_binary_variables <- length(binary_positions)
verbose <- identical(control$verbose, TRUE)
mat <- x$constraints$L
.make_additional_constraint_matrix <-
if(is.simple_triplet_matrix(mat)) {
function(bpos) {
len <- length(bpos)
simple_triplet_matrix(seq_len(len), bpos,
rep.int(1, len),
len, n_of_variables)
}
} else {
function(bpos) {
len <- length(bpos)
add <- matrix(0, len, n_of_variables)
add[cbind(seq_len(len), bpos)] <- 1
add
}
}
.solve_MILP_with_additional_binary_constraints <-
function(x, bpos, bval) {
len <- length(bpos)
x$constraints <-
L_constraint(L = rbind(mat,
.make_additional_constraint_matrix(bpos)),
dir = c(x$constraints$dir, rep.int("==", len)),
rhs = c(x$constraints$rhs, bval))
ROI_solve(x, solver, as_roi_control(control))
}
.split_single_binary_variable <- function(y, i) {
oi <- order[i]
ind <- order[seq_len(i)]
pos <- binary_positions[ind]
b <- y$solution[binary_positions]
## Try flipping the i-th binary variable and see whether this
## also delivers optimal solutions.
b[oi] <- 1 - b[oi]
yf <- .solve_MILP_with_additional_binary_constraints(x, pos, b[ind])
if((yf$status$code != 0) || v_is_not_optimal(yf$objval))
list(y)
else
list(y, yf)
}
ylist <- list(y)
## We allow callers to specify the order in which binary splits
## should be attempted (so the i-th split is for binary variable
## order[i], i.e., at binary_positions[order[i]]).
order <- control$order
if(is.null(order))
order <- seq_len(n_of_binary_variables)
for(i in seq_len(n_of_binary_variables)) {
ylist <-
do.call(c, lapply(ylist, .split_single_binary_variable, i))
len <- length(ylist)
if(verbose)
message(gettextf("N_of_binary_variables: %d *** N_of_optimal_branches: %d",
i, len))
if(len >= nos) {
ylist <- ylist[seq_len(nos)]
break
}
}
ylist
}
choose_method <- function(method) {
milp <- OP(1:3, L_constraint(1:3, "<=", 3), types = c("B", "I", "C"))
milp_solvers <- setdiff(ROI_applicable_solvers(milp), "msbinlp")
if ( isTRUE(method %in% milp_solvers) ) {
return( method )
} else {
if ( length(k <- which(ROI_options("solver_selection_table")$MILP %in% milp_solvers)) ) {
new_method <- ROI_options("solver_selection_table")$MILP[head(k, 1)]
} else {
new_method <- head(milp_solvers, 1)
if ( !length(new_method) )
stop("no mixed integer solver found please install at least ",
"one mixed integer solver (e.g. 'ROI.plugin.glpk')")
}
if ( !is.null(method) )
warning("method '", method, "' is not available or applicable ROI uses '", new_method, "' instead!")
}
new_method
}
solve_OP <- function(x, control=list()) {
solver <- ROI_plugin_get_solver_name( getPackageName() )
if ( is.null(control$nsol_max) ) {
nos <- 1L
} else {
if ( is.numeric(control$nsol_max) & isTRUE(control$nsol_max != Inf) ) {
nos <- as.integer(control$nsol_max)
} else {
nos <- NA
}
}
method <- choose_method( control$method )
control[['method']] <- NULL
control[['nsol_max']] <- NULL
out <- .find_up_to_n_binary_MILP_solutions(x, nos = nos,
add = TRUE,
solver = method,
control = control )
class(out) <- c("msbinlp_solution_set", "OP_solution_set")
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.