everything is green

DESCRIPTION

@@ -28,7 +28,11 @@ Imports:
     RcppArmadillo,
     RcppEigen,
     osqp,
-    MASS
+    MASS,
+    Matrix,
+    progress,
+    quadprog,
+    rlist
 Remotes:
     google/patrick
 RoxygenNote: 6.1.1
@@ -39,7 +43,13 @@ Suggests:
     rmarkdown,
     R.rsp,
     testthat,
-    patrick
+    patrick,
+    clusterSim,
+    corrplot,
+    igraph,
+    kernlab,
+    pals,
+    viridis
 VignetteBuilder:
     bookdown,
     knitr,

NAMESPACE

@@ -7,7 +7,6 @@ export(cluster_k_component_graph)
 export(learn_bipartite_graph)
 export(learn_bipartite_k_component_graph)
 export(learn_k_component_graph)
-export(prial)
 export(relative_error)
 importFrom(Rcpp,sourceCpp)
 useDynLib(spectralGraphTopology)

R/RcppExports.R

@@ -1,10 +1,6 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-Theta_update <- function(U, lambda, K, M, beta) {
-    .Call('_spectralGraphTopology_Theta_update', PACKAGE = 'spectralGraphTopology', U, lambda, K, M, beta)
-}
-
 eigval_sym <- function(M) {
     .Call('_spectralGraphTopology_eigval_sym', PACKAGE = 'spectralGraphTopology', M)
 }

R/constrLaplacianRank.R

@@ -8,8 +8,8 @@
 #' @param m the maximum number of possible connections for a given node used
 #'        to build an affinity matrix
 #' @param lmd L2-norm regularization hyperparameter
-#' @param eig_tol value below which eigenvalues are considered to be zero
-#' @param edge_tol value below which edge weights are considered to be zero
+#' @param eigtol value below which eigenvalues are considered to be zero
+#' @param edgetol value below which edge weights are considered to be zero
 #' @param maxiter the maximum number of iterations
 #' @return A list containing the following elements:
 #' \item{\code{Laplacian}}{the estimated Laplacian Matrix}
@@ -48,21 +48,21 @@
 #' plot(net, vertex.label = NA, vertex.size = 3)
 #' dev.off()
 #' @export
-cluster_k_component_graph <- function(Y, k = 1, m = 5, lmd = 1, eig_tol = 1e-9,
-                                      edge_tol = 1e-6, maxiter = 1000) {
+cluster_k_component_graph <- function(Y, k = 1, m = 5, lmd = 1, eigtol = 1e-9,
+                                      edgetol = 1e-6, maxiter = 1000) {
   time_seq <- c(0)
   start_time <- proc.time()[3]
   A <- build_initial_graph(Y, m)
   n <- ncol(A)
-  if (is.null(S0))
-    S <- matrix(1/n, n, n)
-  else
-    S <- S0
+  S <- matrix(1/n, n, n)
   DS <- diag(.5 * colSums(S + t(S)))
   LS <-  DS - .5 * (S + t(S))
   DA <- diag(.5 * colSums(A + t(A)))
   LA <- DA - .5 * (A + t(A))
-  F <- matrix(eigvec_sym(LA)[, 1:k])
+  if (k == 1)
+    F <- matrix(eigvec_sym(LA)[, 1:k])
+  else
+    F <- eigvec_sym(LA)[, 1:k]
   # bounds for variables in the QP solver
   bvec <- c(1, rep(0, n))
   Amat <- cbind(rep(1, n), diag(n))
@@ -80,7 +80,7 @@ cluster_k_component_graph <- function(Y, k = 1, m = 5, lmd = 1, eig_tol = 1e-9,
     LS <- DS - .5 * (S + t(S))
     F <- eigvec_sym(LS)[, 1:k]
     eig_vals <- eigval_sym(LS)
-    n_zero_eigenvalues <- sum(abs(eig_vals) < eig_tol)
+    n_zero_eigenvalues <- sum(abs(eig_vals) < eigtol)
     time_seq <- c(time_seq, proc.time()[3] - start_time)
     pb$tick(token = list(lmd = lmd, null_eigvals = n_zero_eigenvalues))
     if (k < n_zero_eigenvalues)
@@ -91,13 +91,12 @@ cluster_k_component_graph <- function(Y, k = 1, m = 5, lmd = 1, eig_tol = 1e-9,
       break
     lmd_seq <- c(lmd_seq, lmd)
   }
-  LS[abs(LS) < edge_tol] <- 0
+  LS[abs(LS) < edgetol] <- 0
   AS <- diag(diag(LS)) - LS
-  return(list(Laplacian = LS, Adjacency = AS, eigenvalues = eigval_sym(LS),
+  return(list(Laplacian = LS, Adjacency = AS, eigenvalues = eig_vals,
               lmd_seq = lmd_seq, elapsed_time = time_seq))
 }
 
-
 build_initial_graph <- function(Y, m) {
   n <- nrow(Y)
   A <- matrix(0, n, n)

R/learnGraphTopology.R

@@ -79,7 +79,7 @@
 #' @export
 learn_k_component_graph <- function(S, is_data_matrix = FALSE, k = 1, w0 = "naive", lb = 0, ub = 1e4, alpha = 0,
                                     beta = 1e4, beta_max = 1e6, fix_beta = TRUE, rho = 1e-2, m = 7,
-                                    maxiter = 1e4, abstol = 1e-6, reltol = 1e-4, eig_tol = 1e-9,
+                                    maxiter = 1e4, abstol = 1e-6, reltol = 1e-4, eigtol = 1e-9,
                                     record_objective = FALSE, record_weights = FALSE, verbose = TRUE) {
   if (is_data_matrix || ncol(S) != nrow(S)) {
     A <- build_initial_graph(S, m = m)
@@ -146,7 +146,7 @@ learn_k_component_graph <- function(S, is_data_matrix = FALSE, k = 1, w0 = "naiv
                            relerr = 2 * max(werr / (w + w0), na.rm = 'ignore')))
     }
     if (!fix_beta) {
-      n_zero_eigenvalues <- sum(abs(eigvals) < eig_tol)
+      n_zero_eigenvalues <- sum(abs(eigvals) < eigtol)
       if (k <= n_zero_eigenvalues)
         beta <- (1 + rho) * beta
       else if (k > n_zero_eigenvalues)
@@ -200,8 +200,6 @@ learn_k_component_graph <- function(S, is_data_matrix = FALSE, k = 1, w0 = "naiv
 #' @param maxiter the maximum number of iterations
 #' @param abstol absolute tolerance on the weight vector w
 #' @param reltol relative tolerance on the weight vector w
-#' @param record_objective whether to record the objective function values at
-#'        each iteration
 #' @param record_weights whether to record the edge values at each iteration
 #' @param verbose whether to output a progress bar showing the evolution of the
 #'        iterations
@@ -441,7 +439,7 @@ learn_bipartite_graph <- function(S, is_data_matrix = FALSE, z = 0, nu = 1e4, al
 #' bipartite <- graph_from_adjacency_matrix(Atrue, mode = "undirected", weighted = TRUE)
 #' n <- ncol(Laplacian)
 #' Y <- MASS::mvrnorm(40 * n, rep(0, n), MASS::ginv(Laplacian))
-#' graph <- learn_k_component_bipartite(cov(Y), k = 2, beta = 1e2, nu = 1e2, verbose = FALSE)
+#' graph <- learn_bipartite_k_component_graph(cov(Y), k = 2, beta = 1e2, nu = 1e2, verbose = FALSE)
 #' graph$Adjacency[graph$Adjacency < 1e-2] <- 0
 #' # Plot Adjacency matrices: true, noisy, and estimated
 #' corrplot(Atrue / max(Atrue), is.corr = FALSE, method = "square", addgrid.col = NA, tl.pos = "n", cl.cex = 1.25)
@@ -472,7 +470,7 @@ learn_bipartite_k_component_graph <- function(S, is_data_matrix = FALSE, z = 0,
                                               w0 = "naive", m = 7, alpha = 0., beta = 1e4,
                                               rho = 1e-2, fix_beta = TRUE, beta_max = 1e6, nu = 1e4,
                                               lb = 0, ub = 1e4, maxiter = 1e4, abstol = 1e-6,
-                                              reltol = 1e-4, eig_tol = 1e-9,
+                                              reltol = 1e-4, eigtol = 1e-9,
                                               record_weights = FALSE, record_objective = FALSE, verbose = TRUE) {
   if (is_data_matrix || ncol(S) != nrow(S)) {
     A <- build_initial_graph(S, m = m)
@@ -538,7 +536,7 @@ learn_bipartite_k_component_graph <- function(S, is_data_matrix = FALSE, z = 0,
     if (verbose)
       pb$tick(token = list(beta = beta, kth_eigval = eigvals[k], relerr = 2*max(werr/(w + w0), na.rm = 'ignore')))
     if (!fix_beta) {
-      n_zero_eigenvalues <- sum(abs(eigvals) < eig_tol)
+      n_zero_eigenvalues <- sum(abs(eigvals) < eigtol)
       if (k < n_zero_eigenvalues)
         beta <- (1 + rho) * beta
       else if (k > n_zero_eigenvalues)

R/utils.R

@@ -22,7 +22,6 @@ relative_error <- function(Ltrue, Lest) {
 #' @param Lest estimated Laplacian matrix
 #' @param Lscm estimated Laplacian matrix via the generalized inverse of the
 #'        of the sample covariance matrix
-#' @export
 prial <- function(Ltrue, Lest, Lscm) {
   return(100 * (1 - (norm(Lest - Ltrue, type = "F") /
                      norm(Lnaive - Ltrue, type = "F"))^2))

src/RcppExports.cpp

@@ -7,21 +7,6 @@
 
 using namespace Rcpp;
 
-// Theta_update
-Eigen::MatrixXd Theta_update(const Eigen::MatrixXd& U, const Eigen::VectorXd& lambda, const Eigen::MatrixXd& K, const Eigen::MatrixXd& M, const double beta);
-RcppExport SEXP _spectralGraphTopology_Theta_update(SEXP USEXP, SEXP lambdaSEXP, SEXP KSEXP, SEXP MSEXP, SEXP betaSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type U(USEXP);
-    Rcpp::traits::input_parameter< const Eigen::VectorXd& >::type lambda(lambdaSEXP);
-    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type K(KSEXP);
-    Rcpp::traits::input_parameter< const Eigen::MatrixXd& >::type M(MSEXP);
-    Rcpp::traits::input_parameter< const double >::type beta(betaSEXP);
-    rcpp_result_gen = Rcpp::wrap(Theta_update(U, lambda, K, M, beta));
-    return rcpp_result_gen;
-END_RCPP
-}
 // eigval_sym
 arma::vec eigval_sym(arma::mat M);
 RcppExport SEXP _spectralGraphTopology_eigval_sym(SEXP MSEXP) {
@@ -235,7 +220,6 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_spectralGraphTopology_Theta_update", (DL_FUNC) &_spectralGraphTopology_Theta_update, 5},
     {"_spectralGraphTopology_eigval_sym", (DL_FUNC) &_spectralGraphTopology_eigval_sym, 1},
     {"_spectralGraphTopology_eigvec_sym", (DL_FUNC) &_spectralGraphTopology_eigvec_sym, 1},
     {"_spectralGraphTopology_inv_sympd", (DL_FUNC) &_spectralGraphTopology_inv_sympd, 1},

tests/testthat/test-learnGraphTopology.R

@@ -51,14 +51,14 @@ test_that("learn_adjancecy_and_laplacian can learn k-component bipartite graph",
 })
 
 
-test_that("learn_bipartite_k_component_graph_graph converges with simple bipartite graph", {
+test_that("learn_bipartite_k_component_graph converges with simple bipartite graph", {
   w <- c(1, 0, 0, 1, 0, 1)
   Adjacency <- A(w)
   n <- ncol(Adjacency)
   Y <- MASS::mvrnorm(n * 500, rep(0, n), MASS::ginv(L(w)))
   res <- learn_bipartite_k_component_graph(cov(Y))
   expect_that(res$convergence, is_true())
-  expect_that(relative_error(Adjacency, res$Adjacency) < 1e-1, is_true())
+  expect_that(relative_error(Adjacency, res$Adjacency) < 1.5e-1, is_true())
   expect_that(metrics(Adjacency, res$Adjacency, 1e-1)[1] > .9, is_true())
 })