[Pruning] Refactor phasing

R/balance.R

@@ -0,0 +1,63 @@
+#' @name balance2
+#' @title balance2
+#'
+#' @description
+#'
+#' beta reimplementation of balance
+#'
+#' @param gg (gGraph) with node fields
+#' - $cn (numeric) CN guess
+#' - $weight (numeric)
+#' @param lambda
+#' @param marginal
+#' @param phased (logical) default FALSE
+#' @param ism (logical) default TRUE
+#' @param M (numeric) default 1000
+#' @param epgap (numeric) default 0.1
+#' @param trelim (numeric) max size of uncompressed tree in GB
+#' @param nodefileind (numeric) one of 0 (no node file) 1 (in memory compressed) 2 (on disk uncompressed) 3 (on disk compressed) default 1
+#' @return junction-balanced gGraph
+balance2 = function(gg = NULL,
+                    lambda = 100,
+                    marginal = NULL,
+                    phased = FALSE,
+                    ism = TRUE,
+                    M = 1000,
+                    epgap = 0.1,
+                    trelim = 16,
+                    nodefileind = 1) {
+
+}
+
+#' @name check_balance_inputs
+#' @title check_balance_inputs
+#'
+#' @description
+#'
+#' Validate input gGraph for balance
+#'
+#' @param gg (gGraph)
+check_balance_inputs = function(gg) {
+}
+
+#' @name create_node_variables
+#' @title create_node_variables
+#'
+#' @description
+#' create node CN, node residual, and helper variables for nodes
+#' 
+#' @param gg
+#' @return data.table
+create_node_variables = function(gg) {
+}
+
+#' @name create_edge_variables
+#' @title create_edge_variables
+#'
+#' @description
+#' edge CN, edge residual, and helper variables for edges
+#'
+#' @param gg
+#' @return data.table
+create_edge_variables = function(gg) {
+}

R/converters.R

@@ -937,6 +937,7 @@ read.juncs = function(rafile,
             while (grepl("^((#)|(chrom)|(chr))", thisline)) {
                 headers = c(headers, thisline)
                 thisline = readLines(f, 1)
+                if(length(thisline)==0) break
             }
             ln = sum(length(headers), length(thisline))
             while (length(thisline) > 0) {

R/eventCallers.R

@@ -2523,8 +2523,10 @@ amp = function(gg, jcn.thresh = 8, cn.thresh = 2, fbi.cn.thresh = 0.5,  n.jun.hi
     ploidy = gg$nodes$dt[!is.na(cn), sum(cn*as.numeric(width))/sum(as.numeric(width))]
     keep = (gg$nodes$dt$cn/ploidy) > cn.thresh
     gg$clusters(keep)
-    if (!any(!is.na(gg$nodes$dt$cluster)))
+
+    if (!any(!is.na(gg$nodes$dt$cluster))) {
         return(gg)
+    }
 
   tiny = gg$edges$mark(tiny = gg$edges$dt$class %in% c('DEL-like', 'DUP-like') & gg$edges$span <1e4)
   ucl = gg$nodes$dt[!is.na(cluster), .(wid = sum(width)), by = cluster][wid > width.thresh, cluster] %>% sort
@@ -2536,6 +2538,7 @@ amp = function(gg, jcn.thresh = 8, cn.thresh = 2, fbi.cn.thresh = 0.5,  n.jun.hi
       return(NULL)
     if (!length(cl.edges)) 
       return(NULL)
+
     data.table(cluster = cl,
                nodes = paste(cl.nodes$dt$node.id,
                              collapse = ","),
@@ -2545,6 +2548,14 @@ amp = function(gg, jcn.thresh = 8, cn.thresh = 2, fbi.cn.thresh = 0.5,  n.jun.hi
                n.jun = length(cl.edges),
                n.jun.high = sum(cl.edges$dt[, sum(cn > 3)]), 
                max.jcn = max(c(0, cl.edges$dt$cn)),
+               n.nodes = cl.nodes$dt[, .N],
+               n.inv = cl.edges$dt[class == "INV-like", .N],
+               n.del = cl.edges$dt[class == "DEL-like", .N],
+               n.dup = cl.edges$dt[class == "DUP-like", .N],
+               n.tra = cl.edges$dt[class == "TRA-like", .N],
+               width = sum(cl.nodes$dt[, width], na.rm = TRUE),
+               n.chr = length(unique(cl.nodes$dt[, seqnames])),
+               max.loose.cn = max(c(0, cl.nodes$dt[, loose.cn.left], cl.nodes$dt[, loose.cn.right]), na.rm = TRUE),
                max.cn = max(cl.nodes$dt$cn),
                footprint = paste(gr.string(cl.nodes$footprint),
                                  collapse = ","))
@@ -2553,17 +2564,20 @@ amp = function(gg, jcn.thresh = 8, cn.thresh = 2, fbi.cn.thresh = 0.5,  n.jun.hi
   if (nrow(amps))
   {
       if (!mark.nos) {
-          amps = amps[max.jcn >= jcn.thresh,]
+          amps = amps[max.jcn >= jcn.thresh,
+                      .(nodes, edges, fbi.cn, n.jun, n.jun.high, max.jcn)]
           ##amps[max.jcn >= jcn.thresh, ]
       } else {
           ## keep only clusters with a sufficient number of nodes but don't filter by jcn
-          amps = amps[max.jcn >= jcn.thresh | 
-                      (n.jun >= min.jun &
-                       (strsplit(nodes, ",") %>% lapply(length) %>% unlist) >= min.nodes),]
+          amps = amps[max.jcn >= jcn.thresh | (n.jun >= min.jun & n.nodes >= min.nodes)]
+          amps[, fbi.frac := fbi.cn / max.cn]
+          amps[, inv.frac := n.inv / n.jun]
+          amps[, dup.frac := n.dup / n.jun]
+          amps[, del.frac := n.del / n.jun]
+          amps[, tra.frac := n.tra / n.jun]
       }
   }
 
-
   ## implementing decision tree in https://tinyurl.com/srlbkh2
   if (nrow(amps))
   {

R/gGnome.R

@@ -5054,16 +5054,17 @@ gGraph = R6::R6Class("gGraph",
                              ## cap individual signed edge.ids to their capacity
                              ub[1:nrow(ed)] = pmin(ub[1:nrow(ed)], ed[[cfield]])
                            }
-                           ## sol = Rcplex::Rcplex(Amat = Amat,
-                           sol = Rcplex2(Amat = Amat,
-                                          lb = rep(0, ncol(Amat)),
-                                          ub = ub,
-                                          bvec = b$bvec,
-                                          sense = b$sense,
-                                          vtype = 'I', 
-                                          control = list(trace = verbose > 1),
-                                          objsense = ifelse(max, 'max', 'min'),
-                                          cvec = cvec)
+
+                           sol = Rcplex2(cvec = cvec,
+                                         Amat = Amat,
+                                         bvec = b$bvec,
+                                         lb = rep(0, ncol(Amat)),
+                                         ub = ub,
+                                         sense = b$sense,
+                                         vtype = 'I', 
+                                         control = list(trace = verbose > 1),
+                                         objsense = ifelse(max, 'max', 'min'))
+
 
                            ## find edges used in the solution
                            opted = ed[round(sol$xopt[1:.N])!=0, ]
@@ -7023,21 +7024,12 @@ lengths = function(x, use.names = T) {
   UseMethod("lengths")
 }
 
-#' @name lengths
-#' @title lengths
-#' @description
-#' A vector of walk lengths associated with this walk
-#'
-#' @param gWalk a \code{gWalk} object
-#'
-#' @return the number of nodes per walk in the gWalk
-#' @export
+
 `lengths.gWalk` = function(x, use.names = FALSE){
   return(x$lengths)
 }
 
 
-
 #' @name length
 #' @title length
 #' @description
@@ -8352,7 +8344,7 @@ gWalk = R6::R6Class("gWalk", ## GWALKS
                               sense = ll$sense
                               vtype = ll$vtype
                           }
-                          ## sol = Rcplex::Rcplex(cvec = c,
+
                           sol = Rcplex2(cvec = c,
                                         Amat = A,
                                         bvec = b,

R/utils.R

@@ -1,6 +1,73 @@
 ## appease R CMD check vs data.table 
 sid=side1=side2=side_1=side_2=silent=snid=splice.variant=splicevar=str1=str2=strand_1=strand_2=subject.id=suffix=tag=threep.cc=threep.coord=threep.exon=threep.frame=threep.pc=threep.sc=threep.sc.frame=to=transcript.id.x=transcript.id.y=transcript_associated=transcript_id=twidth=tx.cc=tx.ec=tx_strand=tx_strand.x=tx_strand.y=txid=type=uid=uids=val=walk.id=walk.iid=walk.iid.x=walk.iid.y=wkid=NULL
 
+#' @name run.gurobi
+#' @title run.gurobi
+#'
+#' @description
+#' wrapper to run gurobi with CPLEX-like function call for easy switching bw optimizers
+#'
+#' @param cvec (numeric)
+#' @param Amat (sparse matrix)
+#' @param bvec (numeric)
+#' @param Qmat (sparse matrix)
+#' @param lb (numeric)
+#' @param ub (numeric)
+#' @param sense (character)
+#' @param vtype (variable type)
+#' @param objsense (character) default min
+#' @param control (list) should have epgap, tilim, trace, ideally
+#' @param threads (numeric) 
+#' 
+#' @return sol - list with names $x, $epgap, $status
+run.gurobi = function(cvec = NULL,
+                      Amat = NULL,
+                      bvec = NULL,
+                      Qmat = NULL,
+                      lb = NULL,
+                      ub = NULL,
+                      sense = NULL,
+                      vtype = NULL,
+                      objsense = 'min',
+                      control = list(epgap = 1e-2, tilim = 360, trace = 2),
+                      threads = 32) {
+
+    ## build model
+    model = list(
+        obj = cvec,
+        A = Amat,
+        rhs = bvec,
+        Q = Qmat,
+        lb = lb,
+        ub = ub,
+        vtype = vtype,
+        sense = c("E"="=", "G"=">", "L"="<")[sense], ## inequalities are leq, geq (e.g. not strict)
+        modelsense = objsense)
+
+    ## params
+    params = list()
+    if (!is.null(control$epgap)) {
+        params$MIPGap = control$epgap
+    }
+    if (!is.null(control$tilim)) {
+        params$TimeLimit = control$tilim
+    }
+    if (!is.null(control$trace)) {
+        params$LogToConsole = ifelse(control$trace > 0, 1, 0)
+    }
+
+    ## TODO: set up env list for running on compute cluster
+
+    ## run gurobi
+    sol = gurobi::gurobi(model = model, params = params)
+
+    ## make solution consistent with Rcplex output
+    ## but return all the things
+    sol$epgap = sol$mipgap
+
+    return(sol)
+}
+
 
 #' @name duplicated.matrix
 #' @title R-3.5.1 version of duplicated.matrix