Merge pull request #31 from dpmcsuss/dev

Merge from dev

DESCRIPTION

@@ -40,4 +40,4 @@ Remotes:
     dpmcsuss/rlapjv
 Encoding: UTF-8
 LazyData: TRUE
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.1

NAMESPACE

@@ -4,11 +4,11 @@ export(bari_start)
 export(best_matches)
 export(center_graph)
 export(check_seeds)
+export(do_lap)
 export(edge_match_info)
 export(get_perm)
 export(graph_match_ExpandWhenStuck)
 export(graph_match_FW)
-export(graph_match_FW_multi)
 export(graph_match_IsoRank)
 export(graph_match_PATH)
 export(graph_match_Umeyama)
@@ -21,6 +21,7 @@ export(largest_common_cc)
 export(match_plot_igraph)
 export(match_plot_matrix)
 export(match_report)
+export(matrix_list)
 export(pad)
 export(rds_from_sim_start)
 export(rds_perm_bari_start)
@@ -35,6 +36,7 @@ export(sample_correlated_ieg_pair)
 export(sample_correlated_rdpg)
 export(sample_correlated_sbm_pair)
 export(sample_correlated_sbm_pair_w_junk)
+export(split_igraph)
 export(splr)
 export(splr_sparse_plus_constant)
 export(splr_to_sparse)

R/center_graph.R

@@ -71,7 +71,7 @@ center_graph <- function(A, scheme = c(-1, 1), use_splr = TRUE){
     } 
   } else {
     stop("scheme must be either 'center', 'naive', ",
-      "a positive integer, or a pair of scales.")
+      "a positive integer, or a pair of scalars.")
   }
   g
 }
@@ -109,6 +109,10 @@ pad <- function(m, nr, nc = nr){
       a = Matrix(rbind2(m@a, Matrix(0, nr, da)), sparse = TRUE),
       b = Matrix(rbind2(m@b, Matrix(0, nc, da)), sparse = TRUE))
   }
+  else if(is(m, "matrix_list")) {
+    m <- lapply(m, function(ml) ml[])
+    matrix_list(lapply(m, pad, nr = nr, nc = nc))
+  }
   else{
     Matrix::bdiag(m, Matrix(0, nr, nc))
   }

R/check_graphs.R

@@ -25,14 +25,14 @@ check_graph <- function(A, B,
 
   # this will make the graphs be matrices if they are igraph objects
   if (is.list(A) && !igraph::is.igraph(A)) {
-    A <- lapply(A, function(Al) Al[])
+    A <- matrix_list(lapply(A, function(Al) Al[]))
   } else {
-    A <- list(A[])
+    A <- matrix_list(list(A[]))
   }
   if ( is.list(B) && !igraph::is.igraph(B)) {
-    B <- lapply(B, function(Bl) Bl[])
+    B <- matrix_list(lapply(B, function(Bl) Bl[]))
   } else {
-    B <- list(B[])
+    B <- matrix_list(list(B[]))
   }
 
   totv1 <- ncol(A[[1]])
@@ -51,7 +51,7 @@ check_graph <- function(A, B,
         "square matrices for matching.")
     }
     if (any(sapply(B, function(Bl) nrow(Bl) != totv2))) {
-      stop("B is not square. graph_match_FW only supports ",
+      stop("B is not square. This method only supports ",
         "square matrices for matching.")
     }
   } else {
@@ -60,28 +60,33 @@ check_graph <- function(A, B,
 
 
   try({
-    A <- lapply(A, function(Al) as(Al, "dgCMatrix"))
-    B <- lapply(B, function(Bl) as(Bl, "dgCMatrix"))
+    A <- matrix_list(lapply(A, function(Al) as(Al, "dgCMatrix")))
+    B <- matrix_list(lapply(B, function(Bl) as(Bl, "dgCMatrix")))
   }, silent = TRUE)
-  try({B <- as(B, "dgCMatrix")}, silent = TRUE)
+  # try({
+  #   A <- as(A, "dgCMatrix")
+  #   B <- as(B, "dgCMatrix")
+  # }, silent = TRUE)
 
   if (same_order) {
     if (totv1 > totv2) {
       diff <- totv1 - totv2
-      B <- lapply(B, function(Bl)
-        pad(Bl[], diff))
+      B <- pad(B, diff)
+      # B <- lapply(B, function(Bl)
+      #   pad(Bl[], diff))
     }else if (totv1 < totv2) {
       diff <- totv2 - totv1
-      A <- lapply(A, function(Al)
-        pad(Al[], diff))
+      A <- pad(A, diff)
+      # A <- lapply(A, function(Al)
+      #   pad(Al[], diff))
     }
   }
 
-  if (! as_list) {
+  if (!as_list) {
     if (length(A) > 1) {
       stop("A is multi-layer and must be converted to single layer.\
        (check_graph: is_list = FALSE)")
-    } else if (length(A) > 1) {
+    } else if (length(B) > 1) {
       stop("B is multi-layer and must be converted to single layer.\
        (check_graph: is_list = FALSE)")
     } else {

R/check_sim.R

@@ -43,8 +43,7 @@ check_sim <- function(sim, seeds, nonseeds, totv1, totv2){
   } 
 
   # otherwise, things seem wrong
-  stop(paste0("Similarity matrix must be either NULL or ", 
-      "a square matrix of dimension equal to the number of nonseeds, ", 
+  stop(paste0("Square similarity matrices must have dimension equal to the number of nonseeds, ", 
       nn, ", or the total number of vertices, ", nv, "."))
 
 }

R/gradient_seed_to_nonseed.R

@@ -1,6 +1,6 @@
 # Function to compute the seed to non-seed portion of the gradient
 get_s_to_ns <- function(Alist, Blist, seeds, nonseeds,
-    perm = seq(sum(seeds))){
+    perm = seq(sum(seeds))) {
 
   # NEED TO CHANGE ???
   nns <- nrow(nonseeds)
@@ -16,11 +16,7 @@ get_s_to_ns <- function(Alist, Blist, seeds, nonseeds,
     Bsn <- B[seeds$B,nonseeds$B] %*% t(pmat)
     Bns <- pmat %*% B[nonseeds$B,seeds$B]
 
-    if (ns == 1){
-      outer(Ans, Bns) + outer(Asn, Bsn)
-    } else {
-      (Ans %*% t(Bns)) + (t(Asn) %*% Bsn)
-    }
+    tcrossprod(Ans, Bns) + crossprod(Asn, Bsn)
   }
 
   if (!is(Alist, "list")){

R/graph_match_methods.R

@@ -58,7 +58,7 @@
 #'
 graph_match_FW <- function(A, B, seeds = NULL,
   start = "convex", max_iter = 20,
-  similarity = NULL, lap_method = NULL){
+  similarity = NULL, lap_method = NULL) {
 
 
   graph_pair <- check_graph(A, B)
@@ -88,7 +88,6 @@ graph_match_FW <- function(A, B, seeds = NULL,
 
   # seed to non-seed info
   s_to_ns <- get_s_to_ns(A, B, seeds, nonseeds, rp)
-
   P <- P[, rp]
 
   zero_mat <- Matrix::Matrix(0, nn, nn)
@@ -97,16 +96,15 @@ graph_match_FW <- function(A, B, seeds = NULL,
   similarity <- similarity %*% Matrix::t(rpmat)
 
   # keep only nonseeds
-  A <- lapply(A, function(Al) Al[nonseeds$A, nonseeds$A])
-  B <- lapply(B, function(Bl) Bl[nonseeds$B, nonseeds$B][rp, rp])
+  A <- A[nonseeds$A, nonseeds$A]
+  B <- B[nonseeds$B, nonseeds$B][rp, rp]
   nc <- length(A)
 
   lap_method <- set_lap_method(lap_method, totv1, totv2)
 
 
 
   while(toggle && iter < max_iter){
-
     iter <- iter + 1
     # non-seed to non-seed info
     tAnn_P_Bnn <- zero_mat
@@ -169,7 +167,7 @@ graph_match_FW <- function(A, B, seeds = NULL,
   corr[seeds$A] <- seeds$B
   P <- Matrix::Diagonal(nv)[corr, ]
   D <- P
-  D[nonseeds$A, nonseeds$B] <- as.matrix(D_ns %*% rpmat)
+  D[nonseeds$A, nonseeds$B] <- D_ns %*% rpmat
 
   cl <- match.call()
   list(
@@ -225,7 +223,7 @@ graph_match_convex <- function(A, B, seeds = NULL,
 
   # make a random permutation
   rp <- sample(nn)
-  rpmat <- Matrix::Diagonal(nn)[rp, ]
+  rpmat <- Matrix::Diagonal(nn)[rp, ] 
 
   Asn <- A[seeds$A,nonseeds$A]
   Ann <- A[nonseeds$A,nonseeds$A]
@@ -261,11 +259,15 @@ graph_match_convex <- function(A, B, seeds = NULL,
 
 
   lap_method <- set_lap_method(lap_method, totv1, totv2)
-
+  alpha_seq <- NULL
+  Pseq <- list()
   while(toggle && iter < max_iter){
     f_old <- f
     iter <- iter + 1
-    Grad <- ml_sum(AtA%*%P + P%*%BBt - ABns_sn - t(Ann)%*%P%*%Bnn - Ann%*%P%*%t(Bnn) + similarity)
+    Grad <- -ml_sum(
+      AtA %*% P + P %*% BBt - ABns_sn +
+      -t(Ann) %*% P %*% Bnn - Ann %*% P %*% t(Bnn) +
+      similarity)
 
 
     corr <- do_lap(Grad, lap_method)
@@ -276,26 +278,28 @@ graph_match_convex <- function(A, B, seeds = NULL,
     Cnn <- Ann %*% (P - Pdir) - (P - Pdir) %*% Bnn
     Dnn <- Ann %*% Pdir - Pdir %*% Bnn
 
-    if(ns > 0){
-      Cns <- -(P-Pdir) %*% Bns
-      Csn <- Asn %*% (P-Pdir)
+    if(ns > 0) {
+      Cns <- - (P - Pdir) %*% Bns
+      Csn <- Asn %*% (P - Pdir)
 
       Dns <- Ans - Pdir %*% Bns
       Dsn <- Asn %*% Pdir - Bsn
-    }else{
-      Dns <- Dsn <-Cns <- Csn <- 0
+    } else {
+      Dns <- Dsn <- Cns <- Csn <- 0
     }
-    aq <- innerproduct(Cnn, Cnn) +
+    aq <- innerproduct(Cnn, Cnn) +    
       innerproduct(Cns, Cns) +
       innerproduct(Csn, Csn)
     bq <- innerproduct(Cnn, Dnn) +
       innerproduct(Cns, Dns) +
       innerproduct(Csn, Dsn)
-    aopt <- ifelse(aq == 0 && bq == 0, 0, -bq/aq)
-
+    aopt <- ifelse(aq == 0 && bq == 0, 0,
+      ifelse(-bq / aq > 1, 1, -bq/aq))
+    alpha_seq <- c(alpha_seq, aopt)
+    Pseq <- c(Pseq, Pdir)
     P_new <- aopt * P + (1 - aopt) * Pdir
     f <- innerproduct(Ann %*% P_new - P_new %*% Bnn,
-      Ann %*% P_new - P_new %*% Bnn)
+      Ann %*% P_new - P_new %*% Bnn) + innerproduct(ABns_sn, P_new)
 
     f_diff <- abs(f - f_old)
     P_diff <- norm(P - P_new, "f")
@@ -316,8 +320,28 @@ graph_match_convex <- function(A, B, seeds = NULL,
   corr[nonseeds$A] <- nonseeds$B[corr_ns]
   corr[seeds$A] <- seeds$B
   P <- Matrix::Diagonal(nv)[corr, ]
-  D <- P
-  D[nonseeds$A, nonseeds$B] <- as.matrix(D_ns %*% rpmat)
+  # D <- P
+  # D[nonseeds$A, nonseeds$B] <- D_ns %*% rpmat
+  reorderA <- order(c(nonseeds$A, seeds$A))
+  reorderB <- order(c(nonseeds$B, seeds$B))
+
+  D <- pad(D_ns %*% rpmat, ns)[reorderA, reorderB]
+  if (is(D, "splrMatrix")) {
+    D@x[seeds$A, seeds$B] <- P[seeds$A, seeds$B]  
+  } else {
+    D[seeds$A, seeds$B] <- P[seeds$A, seeds$B]
+  }
+
+
+  # get_f <- function(a){
+  #   PP <- a * P + (1 - a) * Pdir
+  #   innerproduct(Ann %*% PP - PP %*% Bnn,
+  #     Ann %*% PP - PP %*% Bnn)
+  # }
+  # tibble(a = seq(0, 1, 0.01)) %>%
+  #   mutate(f = map_dbl(a, get_f)) %>%
+  #   ggplot(aes(a,f)) + geom_point() +
+  #   geom_vline(xintercept = aopt)
 
   cl <- match.call()
   z <- list(
@@ -326,7 +350,10 @@ graph_match_convex <- function(A, B, seeds = NULL,
     ns = ns, 
     P = P,
     D = D,
-    num_iter = iter)
+    num_iter = iter
+    # ,
+    # seq = list(alpha_seq = alpha_seq, Pseq = Pseq)
+  )  
   z
 }