ctlab / mwcsr

An R package for solving maximum-weight connected subrgaph problem and it's variants

mwcsr

A package for solving maximum weight connected subgraph (MWCS) problem and its variants.

Supported MWCS variants are:

• classic (simple) MWCS, where only vertices are weighted;
• budget MWCS, where vertices are parametrized by costs and overall budget is limited;
• generalized MWCS (GMWCS), where both vertices and edges are weighted;
• signal generalized MWCS (SGMWCS), where both vertices and edges are marked with weighted “signals”, and a weight of a subgraph is calculated as a sum of weights of its unique signals.

Currently, four solvers are supported:

• heuristic relax-and-cut solver rmwcs_solver for MWCS and Budget MWCS;
• heuristic relax-and-cut solver rnc_solver for MWCS/GMWCS/SGMWCS;
• heuristic simulated annealing solver annealing_solver for MWCS/GMWCS/SGMWCS;
• exact (if CPLEX library is available) or heuristic (without CPLEX) solver virgo_solver for MWCS/GMWCS/SGMWCS.

Installation

The package can be installed from GitHub using devtools:

library(devtools)
install_github("ctlab/mwcsr")

Quick start

Load mwcsr, as well as igraph package, which contains functions for graph manipulations.

library(mwcsr)
library(igraph)

Let’s load an example instance of MWCS problem. The instance is a simple igraph object with weight vertex attribute.

data("mwcs_example")
print(mwcs_example)

## IGRAPH f86c56f UN-- 194 209 -- 
## + attr: name (v/c), label (v/c), weight (v/n), label (e/c)
## + edges from f86c56f (vertex names):
##  [1] C00022_2--C00024_0  C00022_0--C00024_1  C00025_0--C00026_0 
##  [4] C00025_1--C00026_1  C00025_2--C00026_2  C00025_4--C00026_4 
##  [7] C00025_7--C00026_7  C00024_1--C00033_0  C00024_0--C00033_1 
## [10] C00022_0--C00041_0  C00022_1--C00041_1  C00022_2--C00041_2 
## [13] C00036_0--C00049_0  C00036_1--C00049_1  C00036_2--C00049_2 
## [16] C00036_4--C00049_4  C00037_1--C00065_0  C00037_0--C00065_1 
## [19] C00022_0--C00074_5  C00022_1--C00074_6  C00022_2--C00074_7 
## [22] C00024_0--C00083_0  C00024_1--C00083_1  C00026_1--C00091_0 
## + ... omitted several edges

summary(V(mwcs_example)$weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.7379 -0.7379  1.9294  5.9667  7.2931 38.1546

Now let us initialize a heuristic relax-and-cut MWCS solver (Alvarez-Miranda and Sinnl, 2017):

rcsolver <- rmwcs_solver()

Now we can use this solver to solve the example instance:

m <- solve_mwcsp(rcsolver, mwcs_example)
print(m$graph)

## IGRAPH 4fe8482 UN-- 162 164 -- 
## + attr: name (v/c), label (v/c), weight (v/n)
## + edges from 4fe8482 (vertex names):
##  [1] C00022_0--C00024_1  C00022_0--C00041_0  C00022_0--C00074_5 
##  [4] C00022_0--C00149_0  C00022_1--C00041_1  C00022_1--C00074_6 
##  [7] C00022_1--C00149_2  C00022_2--C00024_0  C00022_2--C00041_2 
## [10] C00022_2--C00074_7  C00022_2--C00149_1  C00024_0--C00033_1 
## [13] C00024_0--C00222_0  C00024_1--C00033_0  C00024_1--C00222_2 
## [16] C00025_0--C00026_0  C00025_1--C00026_1  C00025_2--C00026_2 
## [19] C00025_4--C00026_4  C00025_7--C00026_7  C00026_0--C00091_1 
## [22] C00026_0--C00311_1  C00026_1--C00091_0  C00026_1--C00311_0 
## + ... omitted several edges

print(m$weight)

## [1] 1178.432

Using exact CPLEX-based Virgo solver

The mwcsr package also provide and interface to exact CPLEX-based Virgo solver
(https://github.com/ctlab/virgo-solver) which can be used to solve MWCS, GMWCS and SGMWCS instances to provable optimality.

To setup this solver CPLEX libraries has to be available. CPLEX can be downloaded from the official web-site: https://www.ibm.com/products/ilog-cplex-optimization-studio. Free licence can be obtained for academic purposes.

First, initialize cplex_dir variable to contain path to CPLEX libraries (for example, /opt/ibm/ILOG/CPLEX_Studio129/).

cplex_dir <- '<replace with path to cplex folder>'

Then initialize the solver:

vsolver <- virgo_solver(cplex_dir=cplex_dir)

And run it on the same MWCS instance:

m <- solve_mwcsp(vsolver, mwcs_example)
print(m$graph)

## IGRAPH 9e2522d UN-- 164 166 -- 
## + attr: name (v/c), label (v/c), weight (v/n), label (e/c)
## + edges from 9e2522d (vertex names):
##  [1] C00022_2--C00024_0  C00022_0--C00024_1  C00025_0--C00026_0 
##  [4] C00025_1--C00026_1  C00025_2--C00026_2  C00025_4--C00026_4 
##  [7] C00025_7--C00026_7  C00024_1--C00033_0  C00024_0--C00033_1 
## [10] C00022_0--C00041_0  C00022_1--C00041_1  C00022_2--C00041_2 
## [13] C00036_0--C00049_0  C00036_1--C00049_1  C00036_2--C00049_2 
## [16] C00036_4--C00049_4  C00037_1--C00065_0  C00037_0--C00065_1 
## [19] C00022_0--C00074_5  C00022_1--C00074_6  C00022_2--C00074_7 
## [22] C00026_1--C00091_0  C00042_2--C00091_0  C00026_0--C00091_1 
## + ... omitted several edges

While the solution is a bit different its weight is the same as found before. The solutions differs only in zero-weight vertices.

get_weight(m$graph)

## [1] 1178.432

However, Virgo guarantees that the result is optimal, unless the solver was interrupted on time limit.

m$solved_to_optimality

## [1] TRUE

Next, consider a GMWCS instance which additionally has edge weights:

data("gmwcs_example")
gmwcs_example

## IGRAPH f86c56f UNW- 194 209 -- 
## + attr: name (v/c), label (v/c), weight (v/n), label (e/c), weight
## | (e/n)
## + edges from f86c56f (vertex names):
##  [1] C00022_2--C00024_0 C00022_0--C00024_1 C00025_0--C00026_0 C00025_1--C00026_1
##  [5] C00025_2--C00026_2 C00025_4--C00026_4 C00025_7--C00026_7 C00024_1--C00033_0
##  [9] C00024_0--C00033_1 C00022_0--C00041_0 C00022_1--C00041_1 C00022_2--C00041_2
## [13] C00036_0--C00049_0 C00036_1--C00049_1 C00036_2--C00049_2 C00036_4--C00049_4
## [17] C00037_1--C00065_0 C00037_0--C00065_1 C00022_0--C00074_5 C00022_1--C00074_6
## [21] C00022_2--C00074_7 C00024_0--C00083_0 C00024_1--C00083_1 C00026_1--C00091_0
## [25] C00042_2--C00091_0 C00026_0--C00091_1 C00042_1--C00091_1
## + ... omitted several edges

summary(E(gmwcs_example)$weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.2715 -0.5762 -0.1060  0.5452  1.0458  8.5829

The same solver can be used to solve this instance:

m <- solve_mwcsp(vsolver, gmwcs_example)
print(m$graph)

## IGRAPH c12fb8d UNW- 182 181 -- 
## + attr: name (v/c), label (v/c), weight (v/n), label (e/c), weight
## | (e/n)
## + edges from c12fb8d (vertex names):
##  [1] C00022_2--C00024_0  C00022_0--C00024_1  C00025_0--C00026_0 
##  [4] C00025_1--C00026_1  C00025_2--C00026_2  C00025_4--C00026_4 
##  [7] C00025_7--C00026_7  C00024_1--C00033_0  C00024_0--C00033_1 
## [10] C00022_0--C00041_0  C00022_1--C00041_1  C00022_2--C00041_2 
## [13] C00036_0--C00049_0  C00036_1--C00049_1  C00036_2--C00049_2 
## [16] C00036_4--C00049_4  C00037_1--C00065_0  C00037_0--C00065_1 
## [19] C00022_0--C00074_5  C00022_1--C00074_6  C00022_2--C00074_7 
## + ... omitted several edges

get_weight(m$graph)

## [1] 1296.41

Finally, let consider an SGMWCS instance. The weights of nodes and edges are defined not directly, but through the signals attribute:

data("sgmwcs_example")
sgmwcs_example

## IGRAPH f86c56f UN-- 194 209 -- 
## + attr: signals (g/n), name (v/c), label (v/c), signal (v/c), label
## | (e/c), signal (e/c)
## + edges from f86c56f (vertex names):
##  [1] C00022_2--C00024_0 C00022_0--C00024_1 C00025_0--C00026_0 C00025_1--C00026_1
##  [5] C00025_2--C00026_2 C00025_4--C00026_4 C00025_7--C00026_7 C00024_1--C00033_0
##  [9] C00024_0--C00033_1 C00022_0--C00041_0 C00022_1--C00041_1 C00022_2--C00041_2
## [13] C00036_0--C00049_0 C00036_1--C00049_1 C00036_2--C00049_2 C00036_4--C00049_4
## [17] C00037_1--C00065_0 C00037_0--C00065_1 C00022_0--C00074_5 C00022_1--C00074_6
## [21] C00022_2--C00074_7 C00024_0--C00083_0 C00024_1--C00083_1 C00026_1--C00091_0
## [25] C00042_2--C00091_0 C00026_0--C00091_1 C00042_1--C00091_1
## + ... omitted several edges

str(V(sgmwcs_example)$signal)

##  chr [1:194] "s1" "s1" "s1" "s2" "s3" "s4" "s4" "s4" "s4" "s4" "s5" "s5" ...

str(E(sgmwcs_example)$signal)

##  chr [1:209] "s103" "s104" "s105" "s106" "s107" "s108" "s109" "s110" "s110" ...

head(sgmwcs_example$signals)

##        s1        s2        s3        s4        s5        s6 
##  5.008879 -0.737898 -0.737898 20.112627 19.890279  2.069292

Again, we can solve this instance with Virgo solver:

m <- solve_mwcsp(vsolver, sgmwcs_example)
print(m$graph)

## IGRAPH f81ae9a UN-- 40 39 -- 
## + attr: signals (g/n), name (v/c), label (v/c), signal (v/c), index
## | (v/n), label (e/c), signal (e/c), index (e/n)
## + edges from f81ae9a (vertex names):
##  [1] C00022_0--C00024_1  C00025_1--C00026_1  C00024_1--C00033_0 
##  [4] C00022_0--C00041_0  C00036_1--C00049_1  C00037_1--C00065_0 
##  [7] C00022_0--C00074_5  C00026_1--C00091_0  C00042_2--C00091_0 
## [10] C00042_2--C00091_51 C00111_6--C00118_2  C00117_6--C00119_4 
## [13] C00042_2--C00122_1  C00022_0--C00149_0  C00122_1--C00149_0 
## [16] C00036_1--C00149_1  C00122_1--C00149_1  C00111_6--C00184_0 
## [19] C00117_6--C00199_1  C00118_2--C00231_1  C00199_1--C00231_1 
## + ... omitted several edges

get_weight(m$graph)

## [1] 270.7676

Running Virgo heuristics without CPLEX

In case CPLEX is not available, Virgo solver can be run in the heuristic mode. Just set cplex_dir parameter to NULL:

vhsolver <- virgo_solver(cplex_dir=NULL)

While the results are not optimal, sometimes they can be enough for practical applications:

m <- solve_mwcsp(vhsolver, mwcs_example)
get_weight(m$graph)

## [1] 1174.737

m$solved_to_optimality

## [1] FALSE

m <- solve_mwcsp(vhsolver, gmwcs_example)
get_weight(m$graph)

## [1] 1276.772

m <- solve_mwcsp(vhsolver, sgmwcs_example)
get_weight(m$graph)

## [1] 227.3432

mwcsr's Issues

virgo_solver crypticly fails when CPLEX Community Edition is used

hello, I am trying to replicate some results produced from the GATOM shiny you developed.
I have followed the recommended installation of cplex and I managed to run the tutorial in the repo

cplex_dir <- "/opt/ibm/ILOG/CPLEX_Studio_Community221/"
vsolver <- virgo_solver(cplex_dir=cplex_dir)
vsolver

MWCS Solver:  virgo_solver 

Parameters:
 name      value                                                         
 cplex_bin /opt/ibm/ILOG/CPLEX_Studio_Community221/cplex/bin/x86-64_linux
 cplex_jar /opt/ibm/ILOG/CPLEX_Studio_Community221/cplex/lib/cplex.jar   
 threads   20                                                            
 timelimit NULL                                                          
 penalty   0                                                             
 memory    2G                                                            
 mst       FALSE                                                         
 log       0  

m <- solve_mwcsp(vsolver, mwcs_example)
print(m$graph)

IGRAPH d543f48 UN-- 164 166 -- 
+ attr: name (v/c), label (v/c), weight (v/n), label (e/c)
+ edges from d543f48 (vertex names):
 [1] C00022_2--C00024_0  C00022_0--C00024_1  C00025_0--C00026_0  C00025_1--C00026_1  C00025_2--C00026_2 
 [6] C00025_4--C00026_4  C00025_7--C00026_7  C00024_1--C00033_0  C00024_0--C00033_1  C00022_0--C00041_0 
[11] C00022_1--C00041_1  C00022_2--C00041_2  C00036_0--C00049_0  C00036_1--C00049_1  C00036_2--C00049_2 
[16] C00036_4--C00049_4  C00037_1--C00065_0  C00037_0--C00065_1  C00022_0--C00074_5  C00022_1--C00074_6 
[21] C00022_2--C00074_7  C00026_1--C00091_0  C00042_2--C00091_0  C00026_0--C00091_1  C00042_1--C00091_1 
[26] C00026_4--C00091_51 C00042_2--C00091_51 C00026_7--C00091_52 C00042_1--C00091_52 C00111_7--C00118_1 
[31] C00111_6--C00118_2  C00111_5--C00118_3  C00042_1--C00122_0  C00042_2--C00122_1  C00022_0--C00149_0 
[36] C00036_0--C00149_0  C00122_1--C00149_0  C00022_2--C00149_1  C00036_1--C00149_1  C00122_1--C00149_1 
+ ... omitted several edges

get_weight(m$graph)

[1] 1178.432

when I try to run the same solver on the graph I am trying to analyze, I get the following error

vsolver <- virgo_solver(cplex_dir=cplex_dir)
res <- solve_mwcsp(solver = vsolver, instance = gs)

Error occurred while solving:Worker 0failed
Error in solver$run_main(args, solver$log) : 
  Failed to run java-based solver.

I can run the approximate solution on the same gs object using the following

vsolver <- virgo_solver(cplex_dir=NULL)
res <- solve_mwcsp(solver = vsolver, instance = gs)
m <- res$graph
print(m)

IGRAPH ada5e05 UN-- 199 211 -- 
+ attr: signals (g/n), name (v/c), metabolite (v/c), element (v/c), label (v/c), url (v/c), pval
| (v/n), origin (v/n), score (v/n), signal (v/c), index (v/n), label (e/c), pval (e/n), origin
| (e/n), Symbol (e/c), gene (e/c), enzyme (e/c), reaction_name (e/c), reaction_equation (e/c), url
| (e/c), reaction (e/c), log2FC (e/n), baseMean (e/n), logPval (e/n), signal (e/c), signalRank
| (e/n), score (e/n), index (e/n)
+ edges from ada5e05 (vertex names):
[1] C00019_5.9438_9.1752 --C00021_5.9438_9.1752  C00025_-0.3248_2.8125--C00026_-0.3248_2.8125
[3] C00025_-2.9228_2.8125--C00026_-0.3248_2.8125 C00025_-0.3248_2.8125--C00026_-2.9228_2.8125
[5] C00025_-2.9228_2.8125--C00026_-2.9228_2.8125 C00025_-4.2219_3.5625--C00026_-4.2219_3.5625
[7] C00025_0.9743_3.5625 --C00026_-4.2219_3.5625 C00025_-4.2219_3.5625--C00026_0.9743_3.5625 
+ ... omitted several edges

this prevents me from reproducing the results from the shiny app, can you help troubleshooting?

here my sessionInfo

R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.9       purrr_0.3.4       readr_2.1.2       tidyr_1.2.0      
 [7] tibble_3.1.8      ggplot2_3.3.6     tidyverse_1.3.1   msigdbr_7.5.1     fgsea_1.22.0      mwcsr_0.1.5      
[13] igraph_1.2.11     data.table_1.14.2 gatom_0.0.0.9000 

loaded via a namespace (and not attached):
 [1] BioNet_1.56.0          bitops_1.0-7           fs_1.5.2               lubridate_1.8.0       
 [5] bit64_4.0.5            RColorBrewer_1.1-3     httr_1.4.3             GenomeInfoDb_1.32.2   
 [9] tools_4.2.1            backports_1.4.1        utf8_1.2.2             R6_2.5.1              
[13] DBI_1.1.3              BiocGenerics_0.42.0    colorspace_2.0-3       withr_2.5.0           
[17] tidyselect_1.1.2       gridExtra_2.3          GGally_2.1.2           bit_4.0.4             
[21] compiler_4.2.1         cli_3.3.0              rvest_1.0.2            Biobase_2.56.0        
[25] network_1.17.2         intergraph_2.0-2       xml2_1.3.3             scales_1.2.0          
[29] digest_0.6.29          XVector_0.36.0         pkgconfig_2.0.3        dbplyr_2.2.1          
[33] fastmap_1.1.0          rlang_1.0.4            readxl_1.4.0           pryr_0.1.5            
[37] rstudioapi_0.13        RSQLite_2.2.15         farver_2.1.1           generics_0.1.3        
[41] jsonlite_1.8.0         statnet.common_4.6.0   vroom_1.5.7            BiocParallel_1.30.3   
[45] RCurl_1.98-1.8         magrittr_2.0.3         GenomeInfoDbData_1.2.8 Matrix_1.4-1          
[49] Rcpp_1.0.9             munsell_0.5.0          S4Vectors_0.34.0       fansi_1.0.3           
[53] babelgene_22.3         lifecycle_1.0.1        stringi_1.7.8          zlibbioc_1.42.0       
[57] plyr_1.8.7             grid_4.2.1             blob_1.2.3             parallel_4.2.1        
[61] ggrepel_0.9.1          crayon_1.5.1           lattice_0.20-45        Biostrings_2.64.0     
[65] haven_2.5.0            hms_1.1.1              KEGGREST_1.36.3        sna_2.7               
[69] pillar_1.8.0           codetools_0.2-18       stats4_4.2.1           fastmatch_1.1-3       
[73] reprex_2.0.1           XML_3.99-0.10          glue_1.6.2             renv_0.15.5           
[77] BiocManager_1.30.18    modelr_0.1.8           png_0.1-7              vctrs_0.4.1           
[81] tzdb_0.3.0             cellranger_1.1.0       gtable_0.3.0           reshape_0.8.9         
[85] assertthat_0.2.1       cachem_1.0.6           broom_1.0.0            coda_0.19-4           
[89] AnnotationDbi_1.58.0   memoise_2.0.1          IRanges_2.30.0         ellipsis_0.3.2 

virgo_solver: libcplex files

Hi! I tried setting the cplex path for the virgo_solver, but I obtain the following error message:

> cplex_dir <-  "C:\\cplex"
> vhsolver <- virgo_solver(cplex_dir)
Error in virgo_solver(cplex_dir) :
  Could not find libcplex files in C:\cplex

This is on Windows 10 and I am using the free community edition of cplex. I do not have files that match pattern = "libcplex\\d+.(so|jnilib|dll)" - I can find those files on UNIX, though. Using Windows, my callable library files seem to match the pattern cplex2010.dll (http://www-eio.upc.edu/lceio/manuals/cplex-11/html/gscplex/gscplexpreface3.html).

Recipe for 'r_helper.o' failed

Hi! I am trying to install the package mwcsr. I tried installing from both GitHub as well as from CRAN. In both cases, I receive the following error message (e.g. if I run install.packages("mwcsr") or install_github("ctlab/mwcsr")).

recipe for target 'r_helper.o' failed
make: *** [r_helper.o] Error 1
ERROR: compilation failed for package ‘mwcsr’

I also tried removing all user-installed packages first and then have mwcsr as the first installed package. However, I still receive the same error.