Working on new functions to manipulate network objects better

.Rbuildignore

@@ -10,3 +10,4 @@
 ^CODE_OF_CONDUCT\.md$
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^makefile$

DESCRIPTION

@@ -7,8 +7,7 @@ Description: The package implements simulation and estimation algorithms for
   statistics. As a difference from the 'ergm' package, 'ergmito' skips the
   MCMC part of the estimation process and goes straight to MLE.
 Depends: R (>= 3.3.0)
-Authors@R: 
-    c(
+Authors@R:  c(
       person(given = "George", family = "Vega Yon", role = c("cre","aut"),
       email = "g.vegayon@gmail.com", comment = c(ORCID = "0000-0002-3171-0844")),
       person("Army Research Laboratory and the U.S. Army Research Office",

NAMESPACE

@@ -10,6 +10,13 @@ S3method(count_stats,list)
 S3method(ergmito_boot,ergmito)
 S3method(ergmito_boot,formula)
 S3method(logLik,ergmito)
+S3method(matrix_to_network,list)
+S3method(matrix_to_network,matrix)
+S3method(nedges,ergmito)
+S3method(nedges,formula)
+S3method(nedges,list)
+S3method(nedges,matrix)
+S3method(nedges,network)
 S3method(nnets,ergmito)
 S3method(nnets,formula)
 S3method(nnets,list)
@@ -38,6 +45,8 @@ export(ergmito_formulae)
 export(exact_gradient)
 export(exact_loglik)
 export(extract.ergmito)
+export(matrix_to_network)
+export(nedges)
 export(new_rergmito)
 export(nnets)
 export(nvertex)

R/RcppExports.R

@@ -31,6 +31,14 @@ exact_gradient. <- function(x, params, weights, statmat) {
     .Call(`_ergmito_exact_gradient`, x, params, weights, statmat)
 }
 
+init_network <- function(n, directed = TRUE, hyper = FALSE, loops = FALSE, multiple = FALSE, bipartite = FALSE) {
+    .Call(`_ergmito_init_network`, n, directed, hyper, loops, multiple, bipartite)
+}
+
+matrix_to_network. <- function(x, directed, hyper, loops, multiple, bipartite) {
+    .Call(`_ergmito_matrix_to_network`, x, directed, hyper, loops, multiple, bipartite)
+}
+
 make_sets <- function(n) {
     .Call(`_ergmito_make_sets`, n)
 }

R/network.R

@@ -0,0 +1,95 @@
+#' Manipulation of network objects
+#' 
+#' This function implements a vectorized version of [network::network]`.adjmat`.
+#' It allows us to turn regular matrices into network objects quickly.
+#' 
+#' @param x Either a single square matrix (adjacency matrix), or a list of these.
+#' @param directed Logical scalar, if `FALSE` then the function only checks the
+#' upper diagonal of the matrix assuming it is undirected.
+#' @param loops Logical scalar. When `FALSE` (default) it will skip the diagonal
+#' of the adjacency matrix.
+#' @param hyper,multiple,bipartite Currently Ignored. Right now all the network objects
+#' created by this function set these parameters as `FALSE`.
+#' 
+#' @details This version does not support adding the name parameter yet. The 
+#' function in the network package includes the name of the vertices as an
+#' attribute.
+#' 
+#' Just like in the network function, `NA` are checked and added accordingly, i.e.
+#' if there is an `NA` in the matrix, then the value is recorded as a missing edge.
+#' 
+#' @return An object of class `network`. This is a list with the following elements:
+#' - `mel` *Master Edge List*: A named list with length equal to the number of edges in
+#' the network. The list itself has 3 elements: `inl` (tail), `outl` (head), and 
+#' `atl` (attribute). By default `atl`, a list itself, has a single element: `na`.
+#' 
+#' - `gal` *Graph Attributes List*: a named list with the following elements:
+#'   - `n` Number of nodes
+#'   - `mnext` Number of edges + 1
+#'   - `directed`,`hyper`,`loops`,`multiple`,`bipartite` The arguments passed to
+#'   the function.
+#'
+#' - `val` *Vertex Attributes List*
+#' 
+#' - `iel` *In Edgest List*
+#' 
+#' - `oel` *Out Edgest List*
+#' 
+#' @examples 
+#' set.seed(155)
+#' adjmats  <- rbernoulli(rep(5, 20))
+#' networks <- matrix_to_network(adjmats)
+#' 
+#' @export
+matrix_to_network <- function(
+  x,
+  directed  = rep(TRUE, length(x)),
+  hyper     = rep(FALSE, length(x)),
+  loops     = rep(FALSE, length(x)),
+  multiple  = rep(FALSE, length(x)),
+  bipartite = rep(FALSE, length(x))
+  ) UseMethod("matrix_to_network")
+
+#' @export
+#' @rdname matrix_to_network
+matrix_to_network.matrix <- function(
+  x,
+  directed  = rep(TRUE, length(x)),
+  hyper     = rep(FALSE, length(x)),
+  loops     = rep(FALSE, length(x)),
+  multiple  = rep(FALSE, length(x)),
+  bipartite = rep(FALSE, length(x))
+  ) {
+
+  matrix_to_network.(
+    x = list(x),
+    directed = directed,
+    hyper = hyper,
+    loops = loops,
+    multiple = multiple,
+    bipartite = bipartite
+    )[[1L]]
+
+}
+
+#' @export
+#' @rdname matrix_to_network
+matrix_to_network.list <- function(
+  x,
+  directed  = rep(TRUE, length(x)),
+  hyper     = rep(FALSE, length(x)),
+  loops     = rep(FALSE, length(x)),
+  multiple  = rep(FALSE, length(x)),
+  bipartite = rep(FALSE, length(x))
+  ) {
+
+  matrix_to_network.(
+    x = x,
+    directed = directed,
+    hyper = hyper,
+    loops = loops,
+    multiple = multiple,
+    bipartite = bipartite
+  )
+
+}

R/sim.R

@@ -7,8 +7,19 @@
 #' @param x An object of class `ergmito_sampler`.
 #' @param sizes Integer vector. Values between 2 to 5 (6 becomes too intensive).
 #' @param mc.cores Integer. Passed to [parallel::mclapply]
+#' @param force Logical. When `FALSE` (default) will try to use `ergmito`'s stat
+#' count functions (see [count_stats]). This means that if one of the requested
+#' statistics in not avialable in `ergmito`, then we will use `ergm` to compute
+#' them, which is significatnly slower (see details).
 #' @param ... Further arguments passed to [ergm::ergm.allstats].
 #' 
+#' @details 
+#' While the \CRANpkg{ergm} package is very efficient, it was not built to do some
+#' of the computations requiered in the ergmito package. This translates in having
+#' some of the functions of the package (ergm) with poor speed performance. This
+#' led us to "reinvent the wheel" in some cases to speed things up, this includes
+#' calculating observed statistics in a list of networks.
+#' 
 #' @return An environment with the following objects:
 #' 
 #' - `calc_prob`
@@ -25,7 +36,8 @@
 #' 
 #' @export
 #' @importFrom parallel mclapply
-new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
+new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,
+                         force = FALSE, ...) {
 
   # environment(model) <- parent.frame()
 
@@ -58,6 +70,7 @@ new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
   ergm_model_attrs <- which(sapply(ergm_model$attrnames, length) > 0)
 
   # Capturing attributes (if any)
+  sampler_w_attributes <- FALSE
   if (length(ergm_model_attrs) && nnets(net) != 1L) {
 
     stop(
@@ -67,6 +80,8 @@ new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
 
   } else if (length(ergm_model_attrs) && nnets(net) == 1L) {
 
+    sampler_w_attributes <- TRUE
+
     for (a in ergm_model_attrs) {
       ergm_model$attrs[[a]] <-
         if (is.null(ergm_model$attrnames[[a]]))
@@ -134,6 +149,11 @@ new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
   } else {
     # THE ERGM WAY -------------------------------------------------------------
 
+    if (!force)
+      stop("To generate this sampler we need to use statnet's ergm functions since",
+           " not all the requested statistics are available in ergmito. If you",
+           " would like to procede, use the option `force = TRUE`", call. = FALSE)
+
     # Are we addinig attributes?
     if (nnets(net) == 1 && network::is.network(net)) {
       attrs <- list()
@@ -263,7 +283,7 @@ new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
     } 
 
     if (!as_indexes) {
-      ans$networks[[s]][
+      nets <- ans$networks[[s]][
         sample.int(
           n       = length(ans$prob[[s]]),
           size    = n,
@@ -272,6 +292,26 @@ new_rergmito <- function(model, theta = NULL, sizes = NULL, mc.cores = 2L,...) {
           useHash = FALSE
         )
         ]
+
+      # If this is a sampler with attributes, then we need to add the attributes
+      # to the sampled graphs
+      if (sampler_w_attributes) {
+
+        nets <- lapply(nets, function(n) {
+
+          # New baseline, we remove all edges
+          net0    <- network::network.copy(net)
+          net0[,] <- 0
+
+          edgelist <- which(n != 0, arr.ind = TRUE)
+          network::add.edges(net0, edgelist[, 2], edgelist[, 1])
+
+        })
+
+      }
+
+      nets
+
     } else {
       sample.int(
         n       = length(ans$prob[[s]]),

R/utils.R

@@ -1,14 +1,50 @@
 
 #' Utility functions to query network dimensions
-#' @param x Either an object of class [ergmito], [network], or [matrix].
+#' @param x Either an object of class [ergmito], [network], [formula], or [matrix].
+#' @param ... Further arguments passed to the method. Currently only `nedges.network`
+#' receives arguments (see [network::network.edgecount]).
 #' @export
 nvertex <- function(x) UseMethod("nvertex")
 
+#' @export
+#' @rdname nvertex
+nedges <- function(x, ...) UseMethod("nedges")
+
+#' @export
+#' @rdname nvertex
+nedges.network <- function(x, ...) {
+  network::network.edgecount(x, ...)
+}
+
+#' @export
+#' @rdname nvertex
+nedges.list <- function(x, ...) {
+  sapply(x, nedges, ...)
+}
+
+#' @export
+#' @rdname nvertex
+nedges.matrix <- function(x, ...) {
+  sum(x != 0)
+}
+
+#' @export
+#' @rdname nvertex
+nedges.ergmito <- function(x, ...) {
+  nedges(x$network, ...)
+}
+
+#' @export
+#' @rdname nvertex
+nedges.formula <- function(x, ...) {
+  nedges(eval(x[[2]]), envir = environment(x))
+}
+
 #' @export
 #' @rdname nvertex
 nvertex.network <- function(x) {
 
-  x$gal$n
+  network::network.size(x)
 
 }
 

README.Rmd

@@ -68,15 +68,15 @@ gplot(net)
 ```
 
 ```{r comparing-w-ergm}
-model <- net ~ edges + mutual + balance
+model <- net ~ edges + mutual + ctriad
 
 library(ergm)
 ans_ergmito <- ergmito(model)
 ans_ergm  <- ergm(model)
 
-# The ergmito should have a larger value
-ergm.exact(ans_ergmito$coef, model)
-ergm.exact(ans_ergm$coef, model)
+# The ergmito should have a larger value when computing exact loglikelihood
+ergm.exact(ans_ergmito$coef, model) >
+  ergm.exact(ans_ergm$coef, model)
 
 summary(ans_ergmito)
 summary(ans_ergm)

data-raw/fivenets.R

@@ -6,26 +6,31 @@ set.seed(12312)
 
 nets <- replicate(5, {
   network(
-    rbernoulli(4, .5),
-    vertex.attr = list(rpois(4, 20)),
+    rbernoulli(5, .5),
+    vertex.attr = list(floor(rnorm(5, 30, 10))),
     vertex.attrnames = "age")
 }, simplify = FALSE)
 
 # Generating the samplers
 examplers <- lapply(nets, function(net) {
 
-  new_rergmito(net ~ edges + nodeicov("age"), theta = c(-4, .2))
+  new_rergmito(net ~ edges + balance, theta = c(-2, 2))
 
 })
 
 # Generating a random sample of these networks (same size)
-fivenets <- lapply(examplers, function(i) i$sample(1L, nvertex(i$network0))[[1L]])
-ans <- ergmito(fivenets ~ edges + nodeicov("age")) 
+fivenets2 <- lapply(examplers, function(i) {
+  # debug(i$sample)
+  i$sample(1L, nvertex(i$network0))[[1L]]
+})
+ans <- ergmito(fivenets2 ~ edges + balance) 
 confint(ans)
+nedges(ans)
 
 usethis::use_data(fivenets)
 
 fivesamplers <- examplers
 usethis::use_data(fivesamplers)
 
 
+\frac{\exp{\transpose{\theta}\stats{y,x}}}{\transpose{\Mat{W}}\left[\exp{\transpose{\theta}\Mat{W}}\right]}

makefile

@@ -0,0 +1,2 @@
+all:
+	cd ../ && R CMD build ergmito/ && R CMD check --as-cran --use-valgrind ergmito_*.tar.gz