#' This function constructs a symphony solver
#' @param verbose would you like to see verbose messages?
#' @import ompr
#' @export
configure_symphony_solver <- function(verbose = FALSE) {
function(model) {
# define the order of all variables
vars <- model@variables
var_list <- unlist(lapply(names(vars), function(key) {
var <- vars[[key]]
if (length(var@instances) == 1 && nchar(var@instances) == 0) {
key
} else {
paste0(key, "_", var@instances)
}
}))
column_types <- unlist(lapply(names(vars), function(key) {
var <- vars[[key]]
rep.int(x = var@type, times = length(var@instances))
}))
ncols <- length(var_list)
nrows <- length(model@constraints)
# build column bounds
column_bounds_l <- unlist(lapply(names(vars), function(key) {
var <- vars[[key]]
lb <- if (length(var@lb) == 0) -Inf else var@lb
if (var@type == "binary") lb <- 0
rep.int(x = lb, times = length(var@instances))
}))
column_bounds_u <- unlist(lapply(names(vars), function(key) {
var <- vars[[key]]
ub <- if (length(var@ub) == 0) Inf else var@ub
if (var@type == "binary") ub <- 1
rep.int(x = ub, times = length(var@instances))
}))
# build objective coeffcient vector
objective <- model@objective
build_coefficent_vector <- function(extracted_coefficients) {
coef_vector <- rep.int(0, ncols)
coefficients <- extracted_coefficients
names(coefficients) <- NULL
obj_constant <- extracted_coefficients$constant
bound_coefs <- unlist(Map(function(var_coef) {
var_ast <- var_coef$ast
if (is.call(var_ast) && length(var_ast) > 1) {
var_name <- as.character(var_ast[[2]])
search_key <- paste0(c(var_name, as.character(var_ast[3:length(var_ast)])), collapse = "_")
} else {
var_name <- as.character(var_ast)
search_key <- var_name
}
setNames(var_coef$coef, search_key)
}, coefficients))
coef_positions <- which(var_list %in% names(bound_coefs))
if (length(coef_positions) > 0) {
coef_vector[coef_positions] <- as.numeric(bound_coefs)
}
coef_vector
}
if (!is.null(objective)) {
coefficients <- extract_coefficients(model@objective@expression[[1]])
obj_constant <- coefficients$constant
if (!is.numeric(obj_constant)) obj_constant <- 0
coefficients <- coefficients$coefficients
names(coefficients) <- NULL
obj_vector <- build_coefficent_vector(coefficients)
} else {
obj_constant <- 0
obj_vector <- rep.int(0, ncols)
}
# build constraint matrix
matrices <- lapply(model@constraints, function(constraint) {
coefficients_lhs <- extract_coefficients(constraint@lhs[[1]])
coefficients_rhs <- extract_coefficients(constraint@rhs[[1]])
direction <- constraint@direction
list(
lhs = build_coefficent_vector(coefficients_lhs$coefficients),
rhs = build_coefficent_vector(coefficients_rhs$coefficients),
direction = direction,
lhs_constant = coefficients_lhs$constant,
rhs_constant = coefficients_rhs$constant
)
})
constraint_matrix <- t(rbind(sapply(matrices, function(constraint) {
constraint$lhs - constraint$rhs
})))
# build row upper bound (aka b)
constraint_rhs <- sapply(matrices, function(constraint) {
constraint$rhs_constant - constraint$lhs_constant
})
constraint_dir <- sapply(matrices, function(constraint) {
constraint$direction
})
binary_vars <-
row_ub <- sapply(matrices, function(constraint) {
if (constraint$direction != ">=") {
constraint$rhs_constant - constraint$lhs_constant
} else {
Inf
}
})
# build row lower bound
row_lb <- sapply(matrices, function(constraint) {
if (constraint$direction != "<=") {
constraint$rhs_constant - constraint$lhs_constant
} else {
-Inf
}
})
symph_column_types <- column_types
symph_column_types[symph_column_types == "binary"] <- "B"
symph_column_types[symph_column_types == "integer"] <- "I"
symph_column_types[symph_column_types == "continuous"] <- "C"
result <- Rsymphony::Rsymphony_solve_LP(obj_vector,
constraint_matrix,
constraint_dir,
constraint_rhs,
bounds = list(
lower = list(
ind = 1:length(column_bounds_l),
val = column_bounds_l
),
upper = list(
ind = 1:length(column_bounds_u),
val = column_bounds_u
)
),
types = symph_column_types,
verbosity = if (verbose) 1 else -2,
max = model@objective@direction == "max")
status <- if (!is.na(result$status) && result$status == 0) "optimal" else "infeasible"
solution <- result$solution
# the solution should be named
names(solution) <- sapply(var_list, function(var) {
if (grepl(x = var, pattern = "_", fixed = TRUE)) {
splited_els <- strsplit(var, "_", fixed = TRUE)[[1]]
paste0(splited_els[1], "[",
paste0(splited_els[2:length(splited_els)], collapse = ","),
"]")
} else {
var
}
})
solution <- new("Solution",
status = status,
model = model,
objective_value = result$objval + obj_constant,
solution = solution)
solution
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.