diff --git a/DESCRIPTION b/DESCRIPTION
index 6119c4e..a50906b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: wdnet
 Title: Weighted and Directed Networks
 Version: 1.1.0
-Date: 2023-05-14
+Date: 2023-05-15
 Authors@R: c(
     person("Yelie", "Yuan", email = "yelie.yuan@uconn.edu",
            role = c("aut", "cre")),
@@ -11,6 +11,7 @@ Authors@R: c(
            role = c("aut")),
     person("Panpan", "Zhang", email = "panpan.zhang@vumc.org",
            role = "aut"))
+Maintainer: Yelie Yuan <yelie.yuan@uconn.edu>
 Description: Implementations of network analysis including
     (1) assortativity coefficient of weighted and directed networks, 
     Yuan, Yan and Zhang (2021) <doi:10.1093/comnet/cnab017>,
diff --git a/R/cls_rpacontrol.R b/R/cls_rpacontrol.R
index 91b8c24..11aaf80 100644
--- a/R/cls_rpacontrol.R
+++ b/R/cls_rpacontrol.R
@@ -34,41 +34,31 @@ is_rpacontrol <- function(control) {
 #' 
 print_control_preference <- function(control, directed = NULL) {
     cat("Preference functions:\n")
-    spref <- " - spref: "
-    tpref <- " - tpref: "
-    pref <- " - pref: "
+    cat(" - ftype: ", control$preference$ftype, "\n", sep = "")
     if (control$preference$ftype == "default") {
-        tmp <- control$preference$sparams
         spref <- paste0(
-            spref,
-            tmp[1], " * outs^", tmp[2], " + ", tmp[3], " * ins^", tmp[4],
-            " + ", tmp[5],
-            sep = ""
+            " - sparams: ",
+            paste(control$preference$sparams, collapse  = " ")
         )
-        tmp <- control$preference$tparams
         tpref <- paste0(
-            tpref,
-            tmp[1], " * outs^", tmp[2], " + ", tmp[3], " * ins^", tmp[4],
-            " + ", tmp[5],
-            sep = ""
+            " - sparams: ",
+            paste(control$preference$tparams, collapse  = " ")
         )
-        tmp <- control$preference$params
         pref <- paste0(
-            pref,
-            "s^", tmp[1], " + ", tmp[2],
-            sep = ""
+            " - params: ",
+            paste(control$preference$params, collapse  = " ")
         )
     } else if (control$preference$ftype == "customized") {
         spref <- paste0(
-            spref,
+            " - spref: ",
             control$preference$spref
         )
         tpref <- paste0(
-            tpref,
+            " - tpref: ",
             control$preference$tpref
         )
         pref <- paste0(
-            pref,
+            " - pref: ",
             control$preference$pref
         )
     } else {
@@ -81,7 +71,6 @@ print_control_preference <- function(control, directed = NULL) {
     } else {
         cat(pref, "\n")
     }
-    cat("\n")
 
     invisible(NULL)
 }
@@ -97,6 +86,7 @@ print_control_preference <- function(control, directed = NULL) {
 print_control_details <- function(x, control_name, control_description) {
     if (control_name == "preference") {
         print_control_preference(control = x, directed = NULL)
+        cat("\n")
         return(invisible(NULL))
     }
     cat(control_description[[control_name]], ":\n", sep = "")
diff --git a/R/cls_wdnet.R b/R/cls_wdnet.R
index e9ed081..16dc72d 100644
--- a/R/cls_wdnet.R
+++ b/R/cls_wdnet.R
@@ -29,20 +29,27 @@ NULL
 #'   weight 1.
 #' @param directed Logical, whether the network is directed (TRUE) or undirected
 #'   (FALSE).
+#' @param nodegroup A numeric vector of node groups.
 #' @param ... Additional components to be added to the \code{wdnet} object.
 #' @return A \code{wdnet} object with the specified \code{edgelist},
 #'   \code{edgeweight} and \code{directed}.
 #' @examples
 #' edgelist <- matrix(c(1, 2, 2, 3, 3, 1), ncol = 2, byrow = TRUE)
 #' edgeweight <- c(1, 2, 3)
-#' netwk <- edgelist_to_wdnet(edgelist, edgeweight)
-#' 
+#' nodegroup <- c(1, 1, 2)
+#' netwk <- edgelist_to_wdnet(
+#'   edgelist = edgelist,
+#'   edgeweight = edgeweight,
+#'   directed = TRUE,
+#'   nodegroup = nodegroup)
+#'
 #' @export
 #' 
 edgelist_to_wdnet <- function(
     edgelist,
     edgeweight,
     directed,
+    nodegroup,
     ...) {
   if (missing(directed) || is.null(directed)) {
     # cat("Assume the network is directed.\n\n")
@@ -64,11 +71,21 @@ edgelist_to_wdnet <- function(
   if (ncol(edgelist) != 2) {
     stop('"edgelist" must have exactly 2 columns.')
   }
+
+  if (missing(nodegroup)) {
+    nodegroup <- NULL
+  }
   
   max_node <- max(edgelist)
   min_node <- min(edgelist)
   num_unique_edges <- length(unique(c(edgelist)))
   
+  if (!is.null(nodegroup)) {
+    if (length(nodegroup) != max_node) {
+      stop('Length of "nodegroup" must match the number of nodes in "edgelist".')
+    }
+  }
+
   if (max_node != num_unique_edges || min_node != 1) {
     stop("Node index must be consecutive integers starting from 1.")
   }
@@ -103,16 +120,13 @@ edgelist_to_wdnet <- function(
       "s" = tmp$s
     )
   }
+  netwk$node.attr$group <- nodegroup
   
   additional_components <- list(...)
   if (length(additional_components) > 0) {
     netwk <- c(netwk, additional_components)
     class(netwk) <- "wdnet"
   }
-  if (!is.null(netwk$nodegroup)) {
-    netwk$node.attr$group <- netwk$nodegroup
-    netwk$nodegroup <- NULL
-  }
   
   if (!is_wdnet(netwk)) {
     stop('Failed to create a valid "wdnet" object.')
@@ -129,6 +143,7 @@ edgelist_to_wdnet <- function(
 #'   (FALSE). If \code{adj} is asymmetric, the network is directed.
 #' @param weighted Logical, whether the network is weighted (TRUE) or unweighted
 #'   (FALSE).
+#' @param nodegroup A numeric vector of node groups.
 #' @param ... Additional components to be added to the \code{wdnet} object.
 #' @return A \code{wdnet} object with the specified \code{adj}.
 #' @export
@@ -140,6 +155,7 @@ adj_to_wdnet <- function(
     adj,
     directed = TRUE,
     weighted = TRUE,
+    nodegroup,
     ...) {
   if (missing(adj) || is.null(adj)) {
     stop('Please provide "adj".')
@@ -165,6 +181,16 @@ adj_to_wdnet <- function(
   
   stopifnot(is.logical(directed))
   stopifnot(is.logical(weighted))
+
+  if (missing(nodegroup)) {
+    nodegroup <- NULL
+  }
+
+  if(!is.null(nodegroup)) {
+    if (length(nodegroup) != nrow(adj)) {
+      stop('Length of "nodegroup" must match the number of nodes in "adj".')
+    }
+  }
   
   tmp <- adj_to_edgelist(
     adj = adj,
@@ -180,6 +206,7 @@ adj_to_wdnet <- function(
     edgelist = edgelist,
     edgeweight = edgeweight,
     directed = directed,
+    nodegroup = nodegroup,
     ...
   )
 }
@@ -196,6 +223,7 @@ adj_to_wdnet <- function(
 #'   create a new \code{wdnet} object.
 #' @param edgelist A two-column matrix representing edges.
 #' @param edgeweight A vector representing the weights of the edges.
+#' @param nodegroup A numeric vector of node groups.
 #' @param directed A logical value indicating whether the network is directed.
 #'   Required if \code{netwk} is \code{NULL}.
 #' @param adj An adjacency matrix.
@@ -210,6 +238,7 @@ create_wdnet <- function(
     netwk,
     edgelist,
     edgeweight,
+    nodegroup,
     directed,
     adj,
     weighted,
@@ -223,6 +252,7 @@ create_wdnet <- function(
           adj = adj,
           directed = directed,
           weighted = weighted,
+          nodegroup = nodegroup,
           ...
         )
       }
@@ -232,6 +262,7 @@ create_wdnet <- function(
         edgelist = edgelist,
         edgeweight = edgeweight,
         directed = directed,
+        nodegroup = nodegroup,
         ...
       )
     }
diff --git a/src/rewire.cpp b/src/rewire.cpp
index ced2c1e..85c3b6b 100644
--- a/src/rewire.cpp
+++ b/src/rewire.cpp
@@ -33,7 +33,7 @@ Rcpp::List dprewire_directed_cpp(
     arma::uvec index_t,
     arma::mat eta, 
     bool rewire_history) {
-  GetRNGstate();
+  // GetRNGstate();
   arma::vec outout(iteration, arma::fill::zeros);
   arma::vec outin(iteration, arma::fill::zeros);
   arma::vec inout(iteration, arma::fill::zeros);
@@ -111,7 +111,7 @@ Rcpp::List dprewire_directed_cpp(
     // r_in_in[n] = (arma::cor(r_sourceIn, r_targetIn)).eval()(0, 0);
   }
   
-  PutRNGstate();
+  // PutRNGstate();
   Rcpp::List ret;
   ret["tnode"] = tnode;
   if (rewire_history) {
@@ -161,7 +161,7 @@ Rcpp::List dprewire_undirected_cpp(
     arma::vec index2,
     arma::mat e, 
     bool rewire_history) {
-  GetRNGstate();
+  // GetRNGstate();
   arma::vec rho(iteration, arma::fill::zeros);
   int nedge = index1.size();
   int e1, e2, temp, count = 0;
@@ -253,7 +253,7 @@ Rcpp::List dprewire_undirected_cpp(
     }
     rho[n] = (arma::cor(degree1, degree2)).eval()(0, 0);
   }
-  PutRNGstate();
+  // PutRNGstate();
   Rcpp::List ret;
   if (rewire_history) {
     ret["history"] = history;
diff --git a/src/rpanet_bag.cpp b/src/rpanet_bag.cpp
index d234e72..c39a9f7 100644
--- a/src/rpanet_bag.cpp
+++ b/src/rpanet_bag.cpp
@@ -26,7 +26,7 @@ Rcpp::List rpanet_bag_cpp(arma::vec snode,
                           double delta_in,
                           bool directed)
 {
-  GetRNGstate();
+  // GetRNGstate();
   int n = scenario.size();
   double u, v;
   int j;
@@ -165,7 +165,7 @@ Rcpp::List rpanet_bag_cpp(arma::vec snode,
     }
     nedge++;
   }
-  PutRNGstate();
+  // PutRNGstate();
 
   Rcpp::List ret;
   ret["snode"] = snode;
diff --git a/src/rpanet_binary_directed.cpp b/src/rpanet_binary_directed.cpp
index c971ea2..956a48e 100644
--- a/src/rpanet_binary_directed.cpp
+++ b/src/rpanet_binary_directed.cpp
@@ -395,7 +395,7 @@ Rcpp::List rpanet_binary_directed(
     updatePrefD(node1, func_type, sparams, tparams, custmSourcePref, custmTargetPref);
   }
   // sample edges
-  GetRNGstate();
+  // GetRNGstate();
   for (i = 0; i < nstep; i++)
   {
     n_reciprocal = 0;
@@ -618,7 +618,7 @@ Rcpp::List rpanet_binary_directed(
       q1.pop();
     }
   }
-  PutRNGstate();
+  // PutRNGstate();
   // free memory (queue)
   queue<node_d *>().swap(q);
   queue<node_d *>().swap(q1);
diff --git a/src/rpanet_binary_undirected.cpp b/src/rpanet_binary_undirected.cpp
index 374f341..244bb6c 100644
--- a/src/rpanet_binary_undirected.cpp
+++ b/src/rpanet_binary_undirected.cpp
@@ -274,7 +274,7 @@ Rcpp::List rpanet_binary_undirected_cpp(
     updatePrefUnd(node1, func_type, params, custmPref);
   }
   // sample edges
-  GetRNGstate();
+  // GetRNGstate();
   for (i = 0; i < nstep; i++)
   {
     m_error = false;
@@ -403,7 +403,7 @@ Rcpp::List rpanet_binary_undirected_cpp(
       q1.pop();
     }
   }
-  PutRNGstate();
+  // PutRNGstate();
   // free memory (queue)
   queue<node_und *>().swap(q);
   queue<node_und *>().swap(q1);
diff --git a/src/rpanet_linear_directed.cpp b/src/rpanet_linear_directed.cpp
index 91357bd..bffb9b2 100644
--- a/src/rpanet_linear_directed.cpp
+++ b/src/rpanet_linear_directed.cpp
@@ -182,7 +182,7 @@ Rcpp::List rpanet_linear_directed_cpp(
 
   // sample edges
   queue<int> q1;
-  GetRNGstate();
+  // GetRNGstate();
   for (i = 0; i < nstep; i++)
   {
     n_reciprocal = 0;
@@ -444,7 +444,7 @@ Rcpp::List rpanet_linear_directed_cpp(
     // checkDiffD(spref, total_spref);
     // checkDiffD(tpref, total_tpref);
   }
-  PutRNGstate();
+  // PutRNGstate();
 
   Rcpp::List ret;
   ret["m"] = m;
diff --git a/src/rpanet_linear_undirected.cpp b/src/rpanet_linear_undirected.cpp
index 1cc4e17..02986f9 100644
--- a/src/rpanet_linear_undirected.cpp
+++ b/src/rpanet_linear_undirected.cpp
@@ -151,7 +151,7 @@ Rcpp::List rpanet_linear_undirected_cpp(
 
   // sample edges
   queue<int> q1;
-  GetRNGstate();
+  // GetRNGstate();
   for (i = 0; i < nstep; i++)
   {
     m_error = false;
@@ -302,7 +302,7 @@ Rcpp::List rpanet_linear_undirected_cpp(
     }
     // checkDiffUnd(pref, total_pref);
   }
-  PutRNGstate();
+  // PutRNGstate();
 
   Rcpp::List ret;
   ret["m"] = m;
diff --git a/src/utils.cpp b/src/utils.cpp
index e3db495..586a13b 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -37,7 +37,7 @@ Rcpp::List find_node_undirected_cpp(arma::vec node1,
                        arma::vec node2, 
                        arma::vec start_edge, 
                        arma::vec end_edge) {
-  GetRNGstate();
+  // GetRNGstate();
   int n = node1.size(), n1 = 0, n2 = 0;
   double u;
   for (int j = 0; j < n; j++) {
@@ -60,7 +60,7 @@ Rcpp::List find_node_undirected_cpp(arma::vec node1,
       n2++;
     }
   }
-  PutRNGstate();
+  // PutRNGstate();
   
   Rcpp::List ret;
   ret["node1"] = node1;
@@ -117,13 +117,13 @@ Rcpp::List node_strength_cpp(arma::vec snode,
 //'
 // [[Rcpp::export]]
 arma::vec sample_node_cpp(arma::vec total_node) {
-  GetRNGstate();
+  // GetRNGstate();
   int n = total_node.size();
   arma::vec nodes(n, arma::fill::zeros);
   for (int i = 0; i < n; i++) {
     nodes[i] = Rcpp::sample(total_node[i], 1)[0];
   }
-  PutRNGstate();
+  // PutRNGstate();
   return nodes;
 }
 
diff --git a/tests/testthat/test-rpanet.R b/tests/testthat/test-rpanet.R
index 10d2228..0636d8a 100644
--- a/tests/testthat/test-rpanet.R
+++ b/tests/testthat/test-rpanet.R
@@ -1,7 +1,7 @@
 test_that("rpanet with default preference functions", {
   # sample PA networks
   set.seed(1234)
-  nstep <- 1e4
+  nstep <- 1e3
   for (method in c("linear", "binary", "bag", "bagx")) {
     if (method == "linear" | method == "binary") {
       control <- rpa_control_preference(
@@ -170,6 +170,9 @@ test_that("rpanet initial network", {
     distribution = rgamma, dparams = list(shape = 5, scale = 0.2)
   ) + rpa_control_newedge(
     distribution = rpois, dparams = list(lambda = 1), shift = 1
+  ) + rpa_control_reciprocal(
+    group.prob = c(0.2, 0.4, 0.4),
+    recip.prob = matrix(rep(0.5, 9), nrow = 3)
   )
 
   netwk1 <- rpanet(1e3, directed = TRUE, control = ctr1)
@@ -185,15 +188,17 @@ test_that("rpanet initial network", {
 
   # check initial netwk
   check_initial_network <- function(netwk1, netwk2) {
-    n <- nrow(netwk1$edgelist)
+    nedge <- nrow(netwk1$edgelist)
+    nnode <- nrow(netwk1$node.attr)
     netwk1$edge.attr$scenario <- 0
     expect_true(all(
-      identical(netwk1$edgelist, netwk2$edgelist[1:n, ]),
-      netwk2$edge.attr$scenario[1:n] == 0,
-      identical(netwk1$edge.attr$weight, netwk2$edge.attr$weight[1:n]),
+      identical(netwk1$edgelist, netwk2$edgelist[1:nedge, ]),
+      netwk2$edge.attr$scenario[1:nedge] == 0,
+      identical(netwk1$edge.attr$weight, netwk2$edge.attr$weight[1:nedge]),
       identical(netwk1$directed, netwk2$directed),
       identical(netwk1$weighted, netwk2$weighted),
-      identical(netwk1$control, netwk2$control)
+      identical(netwk1$control, netwk2$control),
+      identical(netwk1$node.attr$group, netwk2$node.attr$group[1:nnode])
     ))
   }
   check_initial_network(netwk1, netwk2)