Relax and cut solver for sgmwcs problem, documentation, licenses, vignettes.

Author: alexloboda

.Rbuildignore

@@ -3,5 +3,8 @@
 ^\.travis\.yml$
 ^data-raw$
 ^.idea$
+^README.Rmd$
+^src/.idea$
+^src/.+\.o$
 ^lib$
 ^vignettes/tutorial_cache$

DESCRIPTION

@@ -1,14 +1,18 @@
 Package: mwcsr
-Title: Solvers for Various Maximum Weight Connected Subgraph Problems
+Title: Solvers for Maximum Weight Connected Subgraph Problem and its variants
 Version: 0.1.0
 Authors@R: c(person("Alexander", "Loboda", email = "[email protected]", role = "aut"),
     person("Nikolay", "Poperechnyi", email = "[email protected]", role = "aut"),
     person("Eduardo", "Alvarez-Miranda", email = "[email protected]", role = "aut"), 
     person("Markus", "Sinnl", email = "[email protected]", role = "aut"),
-    person("Alexey", "Sergushichev", email = "[email protected]", role = c("aut", "cre")))
-Description: The package implements algorithms for solving various Maximum
+    person("Alexey", "Sergushichev", email = "[email protected]", role = c("aut", "cre")),
+    person("Paul", "Hosler Jr.", role = "cph", comment = "Bundled JOpt Simple package"),
+    person("www.hamcrest.org", role = "cph", comment = "Bundled hamcrest package"),
+    person("Barak Naveh and Contributors", role = "cph", comment = "Bundled JGraphT package"),
+    person("The Apache Software Foundation", role = "cph", comment = "Bundled Apache Commons Math package"))
+Description: Algorithms for solving various Maximum
     Weight Connected Subgraph Problems, including variants with 
-    budget constraints, cardinallity constraints, weighted edges.
+    budget constraints, cardinallity constraints and weighted edges.
 Depends: R (>= 3.5)
 Imports:
     methods,
@@ -17,13 +21,16 @@ Imports:
 Suggests:
     knitr,
     rmarkdown,
+    mathjaxr,
     testthat,
     BioNet, 
+    roxygen2, 
     DLBCL   
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
+Roxygen: list(markdown = TRUE)
 VignetteBuilder: knitr
 URL: https://github.com/ctlab/mwcsr
 BugReports: https://github.com/ctlab/mwcsr/issues

LICENSE

@@ -1,2 +1,2 @@
 YEAR: 2017
-COPYRIGHT HOLDER: Alexey Sergushichev
+COPYRIGHT HOLDER: Alexey Sergushichev, Alexander Loboda, Nikolay Poperechnyi, Markus Sinnl, Eduardo Alvarez-Miranda

LICENSE.note

@@ -4,11 +4,13 @@ S3method("[",mwcs_solver)
 S3method("[<-",mwcs_solver)
 S3method(parameters,default)
 S3method(parameters,rmwcs_solver)
+S3method(parameters,rnc_solver)
 S3method(parameters,simulated_annealing_solver)
 S3method(parameters,virgo_solver)
 S3method(print,mwcs_solver)
 S3method(print,mwcs_solver_params)
 S3method(solve_mwcsp,rmwcs_solver)
+S3method(solve_mwcsp,rnc_solver)
 S3method(solve_mwcsp,simulated_annealing_solver)
 S3method(solve_mwcsp,virgo_solver)
 export("timelimit<-")
@@ -18,6 +20,7 @@ export(get_weight)
 export(normalize_sgmwcs_instance)
 export(parameters)
 export(rmwcs_solver)
+export(rnc_solver)
 export(set_parameters)
 export(solve_mwcsp)
 export(virgo_solver)

NAMESPACE

@@ -5,6 +5,10 @@ sa_solve <- function(instance, solver) {
     .Call('_mwcsr_sa_solve', PACKAGE = 'mwcsr', instance, solver)
 }
 
+sgmwcs_solve <- function(instance, solver_params) {
+    .Call('_mwcsr_sgmwcs_solve', PACKAGE = 'mwcsr', instance, solver_params)
+}
+
 rmwcs_solve <- function(network, params) {
     .Call('_mwcsr_rmwcs_solve', PACKAGE = 'mwcsr', network, params)
 }

R/RcppExports.R

@@ -9,8 +9,7 @@ check_sa_solver <- function(solver) {
 
 #' @export
 parameters.simulated_annealing_solver <- function(solver) {
-    params(parameter("normalization", type = "logical"),
-        parameter("schedule", type = "mc", mc = schedules),
+    params(parameter("schedule", type = "mc", mc = schedules),
         parameter("initial_temperature", type = "float", positive = TRUE),
         parameter("final_temperature", type = "float", positive = TRUE),
         parameter("verbose", type = "logical"))
@@ -22,72 +21,70 @@ parameters.simulated_annealing_solver <- function(solver) {
 #' Typically, it allows to find some good solution in a short time. This
 #' implementation doesn't compute any upper bound on solution, so there is no
 #' guarantee of optimality of solution provided.
-#' @param normalization adjust all weights so that all possible changes have
-#' mean square 1.0
+#'
+#' Algorithm maintains connected subgraph staring with empty subgraph.
+#' Each iteration one random action is considered. It may be a removal of a
+#' vertex or an edge which does not separate graph or addition of an extra vertex or
+#' an edge neighboring existing graph. If the subgraph is empty all vertices
+#' are considered as candidates to form a subgaph. After candidate is chosen two
+#' subgraph scores are considered: for a new subgraph and for an old one. Simulated
+#' annealing operates with a notion of temperature. The candidate action is
+#' accepted with probability: p(S'|S) = exp(-E / T), where E is weight difference
+#' between subgraphs and T is current temperature.
+#'
+#' Temperature is calculated in each iteration. in `mwcsr` there are two
+#' temperature schedules supported. So called Boltmann annealing uses the formula:
+#' T(k) = T0 / (ln(1 + k)), while in case of fast annealing this one is used:
+#' T(k) = T0 / k, where k is iteration number.
+#'
+#' To tune the algorithm it is useful to realize how typical changes in the goal
+#' function for single actions are distributed. Calculating the acceptance probabilities
+#' at initial temperature and final temperatures may help to choose schedule and
+#' temperatures.
+#'
+#' @seealso [rnc_solver] will probably be a better choice with minimal tuning necessary
 #' @param schedule boltzmann annealing or fast annealing
 #' @param initial_temperature initial value for the temperature
 #' @param final_temperature final value for the temperature
 #' @param verbose whether be verbose or not
 #' @export
-annealing_solver <- function(normalization = TRUE,
-                             schedule = c("fast", "boltzmann"),
+annealing_solver <- function(schedule = c("fast", "boltzmann"),
                              initial_temperature = 1.0,
                              final_temperature = 1e-6,
                              verbose = FALSE) {
     solver_ctor(c(sa_class, mwcs_solver_class))
 }
 
-normalize_weights <- function(instance) {
-    weights <- c(V(instance)$weight, E(instance)$weight)
-    support <- sum(weights[weights > 0]) / 2
-    if (any(weights < -EPS)) {
-        support <- min(support, -min(weights))
-    }
-    instance
-}
-
-#' Solve generalized maximum weight subgraph problem using simulated annealing
-#'
-#' @inheritParams solve_mwcsp.virgo_solver
-#' @param solver An annealing solver object.
-#' @param instance An igraph instance to work with.
-#' @return An object of class mwcsp_solution.
+#' @rdname solve_mwcsp
+#' @order 5
+#' @param warm_start warm start solution, an object of the class mwcsp_solution.
 #' @export
-solve_mwcsp.simulated_annealing_solver <- function(solver, instance, ...) {
+solve_mwcsp.simulated_annealing_solver <- function(solver, instance, warm_start, ...) {
     if (!inherits(solver, sa_class)) {
         stop("Not a simulated annealing solver")
     }
 
-    inst_type <- get_instance_type(instance)
-
-
-    if (inst_type$type == "SGMWCS" && inst_type$valid) {
-        signal_instance <- instance
-    } else if (inst_type$type == "GMWCS" && inst_type$valid) {
-        signal_instance <- normalize_sgmwcs_instance(instance,
-                                                     nodes.group.by = NULL,
-                                                     edges.group.by = NULL)
-    } else if (inst_type$type == "MWCS" && inst_type$valid) {
-        inst2 <- instance
-        E(inst2)$weight <- 0
-        signal_instance <- normalize_sgmwcs_instance(inst2,
-                                                     nodes.group.by = NULL,
-                                                     edges.group.by = NULL)
-    } else {
-        stop("Not a valid MWCS, GMWCS or SGMWCS instance")
-    }
+    signal_instance <- to_signal_instance(instance)
 
     inst_rep <- instance_from_graph(signal_instance)
     inst_rep[["vertex_signals"]] <- match(igraph::V(signal_instance)$signal, names(signal_instance$signals)) - 1
     inst_rep[["edge_signals"]] <- match(igraph::E(signal_instance)$signal, names(signal_instance$signals)) - 1
     inst_rep[["signal_weights"]] <- signal_instance$signals
 
+    if (!missing(warm_start)) {
+        ws_sol <- warm_start$warm_start_solution
+        inst_rep[["warm_start_vertices"]] <- ws_sol$vertices
+        inst_rep[["warm_start_edges"]] <- ws_sol$edges
+        inst_rep[["warm_start_weight"]] <- ws_sol$weight
+    }
+
     res <- sa_solve(inst_rep, solver)
     if (length(res$edges) == 0) {
         g <- igraph::induced_subgraph(instance, vids = res$vertices)
     } else {
         g <- igraph::subgraph.edges(instance, eids = res$edges)
     }
     weight <- get_weight(g)
-    solution(g, weight, solved_to_optimality = FALSE)
+    res$weight <- weight
+    solution(g, weight, solved_to_optimality = FALSE, warm_start_solution = res)
 }

R/annealing.R

@@ -28,42 +28,80 @@ parameters.rmwcs_solver <- function(solver) {
          parameter("verbose", type = "logical"))
 }
 
-#' Generates a rmwcs solver with corresponding parameters
+#' Generate a rmwcs solver
+#'
+#' The method is based on relax-and-cut approach and allows to solve
+#' Maximum Weight Subgraph Probleam and its budget and cardinality variants.
+#' By constructing lagrangian
+#' relaxation of MWCS problem necessary graph connectivity constraints are introduced
+#' in the objective function giving upper bound on the weight of the optimal
+#' solution. On the other side, primal heuristic uses individul contribution
+#' of the variables to lagrangian relaxation to find possible solution of the initial
+#' problem. The relaxation is then optimized by using iterative subgradient method.
+#'
+#' One iteration of algorithm includes solving lagrangian relaxation and updating
+#' lagrange multipliers. It may also contain cuts (or connectivity constraints) separation process, run of
+#' heuristic method, variable fixing routine. The initial step size for
+#' subgradient method can be passed as `beta` argument. If there is no improvement in
+#' upper bound in consequtive `beta_iterations` iterations the step size is
+#' halved. There are three possible strategies for updating multipliers. See the references
+#' section for the article where differences are discussed.
+#'
+#' Some initial cuts are added at the start of the algorithm if `start_constraints`
+#' is set to \code{TRUE}. Other constraints are separated on the fly and are
+#' unaffected in the next `sep_iter_freeze` iterations of the subgradient mehod.
+#' Then the corresponding lagrange mutipliers are updated from iteration to iteration.
+#' Aging procedure for cuts is incorporated in the algorithm meaning constraint multipliers
+#' are updated for non-violated cuts for up to `max_age` iterations from
+#' the point where a cut was violated last time. There are two separation methods
+#' implemented: fast and strong, where tha latter is supposed to minimize number of
+#' variables used in generated constraint while in the former case there is no need to explore
+#' whole graph to construct a constraint.
+#'
+#' A variant of MST approximation of PCSTP is used as Primal Heuristic.
+#' See references for more details.
+#'
+#' The frequences
+#' of running separation process and heuristic are specified in
+#' `sep_iterations` and `heur_iterations`.
+#'
 #' @param timelimit Timelimit in seconds
-#' @param max_iterations Maximum number of subgradient iterations
+#' @param max_iterations Maximum number of iterations
+#' @param start_constraints Whether to add flow-conservation/degree constraints at start
+#' @param separation Separation: "strong" or "fast"
+#' @param sep_iterations Frequency of separating cuts (in iterations)
 #' @param beta_iterations Number of nonimproving iterations until beta is halved
-#' @param separation Separation: "strong" of "fast"
-#' @param sep_iterations Extending the life of non-violated inequalities
-#' @param start_constraints Whether to add flow-conservation/degree cons at start
-#' @param pegging Pegging (variable fixing)
-#' @param max_age extending the life of non-violated inequalities
-#' @param sep_iter_freeze After how many iterations we are checking added ineqs
-#' @param heur_iterations After how many iterations we are doing heuristics
+#' @param pegging variable fixing
+#' @param max_age number of iterations in aging procedure for non-violated  cuts
+#' @param sep_iter_freeze Number of iterations when a newly separated cut is anaffected by subgradient algorithm.
+#' @param heur_iterations Frequency of calling heuristic method (in iterations)
 #' @param subgradient Subgradient: "classic", "average", "cft"
-#' @param beta Beta for subgradient
-#' @param verbose Whether to print solving progress
+#' @param beta Initial step size of subgradient algorithm
+#' @param verbose Should the solving progress and stats be printed?
+#' @references Álvarez-Miranda, E., & Sinnl, M. (2017). A Relax-and-Cut framework
+#' for large-scale maximum weight connected subgraph problems. Computers & Operations Research, 87, 63-82.
 #' @export
 #' @import igraph
 rmwcs_solver <- function(timelimit = 1800L,
                   max_iterations = 1000L,
-                  beta_iterations = 50L,
+                  beta_iterations = 5L,
                   separation = "strong",
                   start_constraints = TRUE,
                   pegging = TRUE,
-                  max_age = 3,
-                  sep_iterations= 1L,
-                  sep_iter_freeze = 1L,
-                  heur_iterations = 1L,
+                  max_age = 10,
+                  sep_iterations= 10L,
+                  sep_iter_freeze = 50L,
+                  heur_iterations = 10L,
                   subgradient = "classic",
                   beta = 2.0,
                   verbose = FALSE) {
     solver_ctor(c(rmwcs_class, mwcs_solver_class))
 }
 
-#' Solves MWCS problem using relax-and-cut approach
-#' @inheritParams solve_mwcsp
+#' @rdname solve_mwcsp
 #' @param max_cardinality integer maximum number of vertices in solution.
 #' @param budget numeric maximum budget of solution.
+#' @order 3
 #' @export
 solve_mwcsp.rmwcs_solver <- function(solver, instance, max_cardinality = NULL,
                                      budget = NULL, ...) {
@@ -103,5 +141,6 @@ solve_mwcsp.rmwcs_solver <- function(solver, instance, max_cardinality = NULL,
     subgraph <- igraph::induced_subgraph(instance, vs$graph)
     weight <- sum(instance_rep$vertex_weights[vs$graph])
     opt <- isTRUE(all.equal(weight, vs$ub))
-    solution(subgraph, weight, opt, upper_bound = vs$ub)
+    solution(subgraph, weight, opt, upper_bound = vs$ub, incumbent = vs$lb,
+             solved_to_optimality = abs(vs$lb - vs$ub) < EPS)
 }

R/rmwcs.R

@@ -0,0 +1,64 @@
+rnc_class <- "rnc_solver"
+
+check_rnc_solver <- function(solver) {
+    if (!inherits(solver, rnc_class)) {
+        stop("Invalid relax&cut sgmwcs solver provided.")
+    }
+}
+
+#' @export
+parameters.rnc_solver <- function(solver) {
+    params(parameter("max_iterations", type = "integer", positive = TRUE),
+           parameter("beta_iterations", type = "integer", positive = TRUE),
+           parameter("sep_iterations", type = "integer", positive = TRUE),
+           parameter("heur_iterations", type = "integer", positive = TRUE),
+           parameter("max_age", type = "integer", positive = TRUE),
+           parameter("verbose", type = "logical"))
+}
+
+#' Construct relax-and-cut SGMWCS solver
+#'
+#' The solver is based on the same approach as rmwcs solver. Modifications to
+#' the original scheme are introduced to tackle problems arising with introduction
+#' of edge weights and signals. It is recommended to use rmwcs solver to solve MWCS
+#' problems, while due to differences in primal heuristic it may be a good pratice
+#' to run both solvers on the same problem.
+#' @inheritParams rmwcs_solver
+#' @seealso [rmwcs_solver]
+#' @export
+rnc_solver <- function(max_iterations = 1000L,
+                       beta_iterations = 50L,
+                       heur_iterations = 10L,
+                       sep_iterations = 10L,
+                       verbose = FALSE) {
+    solver_ctor(c(rnc_class, mwcs_solver_class))
+}
+
+#' @rdname solve_mwcsp
+#' @order 4
+#' @export
+solve_mwcsp.rnc_solver <- function(solver, instance, ...) {
+    check_rnc_solver(solver)
+
+    signal_instance <- to_signal_instance(instance)
+
+    inst_rep <- instance_from_graph(signal_instance)
+    inst_rep[["vertex_signals"]] <- match(igraph::V(signal_instance)$signal, names(signal_instance$signals)) - 1
+    inst_rep[["edge_signals"]] <- match(igraph::E(signal_instance)$signal, names(signal_instance$signals)) - 1
+    inst_rep[["signal_weights"]] <- signal_instance$signals
+
+    res <- sgmwcs_solve(inst_rep, solver)
+
+    if (length(res$edges) == 0) {
+        g <- igraph::induced_subgraph(instance, vids = res$vertices)
+    } else {
+        g <- igraph::subgraph.edges(instance, eids = res$edges)
+    }
+
+    weight <- get_weight(g)
+
+    stopifnot(abs(weight - res$lb) < EPS)
+
+    solution(g, weight, solved_to_optimality = abs(res$lb - res$ub) < EPS,
+             upper_bound = res$ub)
+}