Merge pull request #36 from statnet/i21-gw-cutoff

Refactor *spL() terms to work like those of 'ergm' and remove *degreeL() terms.

Author: krivit

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ergm.multi
-Version: 0.3.0-4176
-Date: 2025-04-08
+Version: 0.3.0-4188
+Date: 2025-06-13
 Title: Fit, Simulate and Diagnose Exponential-Family Models for Multiple or Multilayer Networks
 Authors@R: c(person(c("Pavel", "N."), "Krivitsky", role=c("aut","cre"), email="pavel@statnet.org", comment=c(ORCID="0000-0002-9101-3362")),
              person(c("Mark", "S."), "Handcock", role=c("ctb"), email="handcock@stat.ucla.edu"),

R/InitErgmConstraint.R

@@ -58,8 +58,8 @@
 ##          n <- network.size(lhs.nw)
 ##          storage.mode(n) <- "integer"
 ##          a <- c(ergm_get_vattr(attr, lhs.nw)) # Strip attributes, which confuse rle().
-##          if(NVL(lhs.nw%n%"bipartite",0)){
-##            bip <- lhs.nw %n% "bipartite"
+##          if(NVL(lhs.b1.size(nw),0)){
+##            bip <- lhs.b1.size(nw)
 ##            ea <- a[seq_len(bip)]
 ##            aa <- a[bip+seq_len(n-bip)]
 ##            if(length(rle(ea)$lengths)!=length(unique(rle(ea)$values)) || length(rle(aa)$lengths)!=length(unique(rle(aa)$values))) stop("Current implementation of block-diagonal sampling requires that the blocks of the egos and the alters be contiguous. See ", sQuote("ergmConstraint?blockdiag"), " for more information.")

R/InitErgmProposal.blockdiag.R

@@ -109,7 +109,7 @@
 #' @keywords internal
 #' @export
 ergm_block_diag_samp_info <- function(nw, a){
-  bip <- nw %n% "bipartite"
+  bip <- b1.size(nw)
   if(bip){
     ea <- a[seq_len(bip)]
     aa <- a[bip+seq_len(network.size(nw)-bip)]
@@ -135,7 +135,9 @@ ergm_block_diag_samp_info <- function(nw, a){
     w <- cumsum(a$lengths*(a$lengths-1)) # cumulative block weights ~ # dyads in the block
     w <- w/max(w)
     # Note that this automagically takes care of singleton blocks by giving them weight 0.
-    out <- list(nblk=as.integer(length(w)),  pos=as.integer(b), cumwt=as.double(w), ndyads=sum(as.double(a$lengths*(a$lengths-1)/(if(is.directed(nw)) 1 else 2))))
+    out <- list(nblk = as.integer(length(w)),  pos = as.integer(b),
+                cumwt = as.double(w),
+                ndyads = sum(choose(a$lengths, 2L) * (1L + is.directed(nw))))
   }
   structure(out, class="ergm_block_diag_samp_info")
 }

R/InitErgmTerm.R

@@ -217,7 +217,9 @@ direct.network <- function(x, rule=c("both", "upper", "lower")){
                upper = cbind(pmin(el[,1],el[,2]),pmax(el[,1],el[,2])),
                lower = cbind(pmax(el[,1],el[,2]),pmin(el[,1],el[,2])))
 
-  o <- network.initialize(network.size(x), directed=TRUE, bipartite=x%n%"bipartite", loops=has.loops(x), hyper=is.hyper(x), multiple=is.multiplex(x))
+  o <- network.initialize(network.size(x), directed = TRUE,
+                          bipartite = b1.size(x), loops = has.loops(x),
+                          hyper = is.hyper(x), multiple = is.multiplex(x))
   o <- network.edgelist(el, o)
   nvattr.copy.network(o, x)
 }
@@ -482,7 +484,7 @@ Layer <- function(..., .symmetric=NULL, .bipartite=NULL, .active=NULL){
   }
 
   # nwl may now be a list with networks of heterogeneous bipartitedness.
-  bip <- sapply(nwl, `%n%`, "bipartite") %>% sapply(NVL, 0L)
+  bip <- sapply(nwl, b1.size) %>% sapply(NVL, 0L)
   blockout <- if(all_identical(bip)) rep(FALSE, length(nwl)) else bip
 
   nwl <- Map(function(nw, b){
@@ -1033,3 +1035,120 @@ InitErgmTerm.twostarL<-function(nw, arglist,  ...) {
   iinputs <- c(typeID, a$distinct)
   list(name="twostarL", coef.names=coef.names, iinputs=iinputs, auxiliaries=auxiliaries, minval=0, dependence=TRUE)
 }
+
+################################################################################
+
+#' @templateVar name mutualL
+#' @title Mutuality
+#' @description In binary ERGMs, equal to the number of
+#'   pairs of actors \eqn{i} and \eqn{j} for which \eqn{(i{\rightarrow}j)}{(i,j)}
+#'   and \eqn{(j{\rightarrow}i)}{(j,i)} both exist.
+#'
+#' @details This term can only be used with directed networks.
+#'
+#' @usage
+#' # binary: mutualL(same=NULL, diff=FALSE, by=NULL, keep=NULL, Ls=NULL)
+#' @param same optional argument. If passed the name of a vertex attribute,
+#'   only mutual pairs that match on the attribute are counted. Only one of `same`
+#'   or `by` may be used. If both parameters are passed, `same` takes precedent. This
+#'   parameter is affected by `diff`.
+#' @param diff separate counts for each unique matching value can be obtained by using
+#'   `diff=TRUE` with `same`.
+#' @param by each node is counted separately for each mutual pair in which it
+#'   occurs and the counts are tabulated by unique values of the attribute if
+#'   passed the name of a vertex attribute. This means that the sum of the mutual statistics when `by` is used
+#'   will equal twice the standard mutual statistic. Only one of `same`
+#'   or `by` may be used. If both parameters are passed, `same` takes precedent. This
+#'   parameter is not affected by `diff`.
+#' @param keep a numerical vector to specify which statistics should be kept whenever the `mutual` term would
+#'   ordinarily result in multiple statistics.
+#' @templateVar Ls.howmany one or two
+#' @templateVar Ls.interp . If given, the statistic will count the number of dyads where a tie in `Ls[[1]]` reciprocates a tie in `Ls[[2]]` and vice versa. (Note that dyad that has mutual ties in `Ls[[1]]` and in `Ls[[2]]` will add 2 to this statistic.) If a formula is given, it is replicated.
+#' @template ergmTerm-Ls
+#'
+#' @template ergmTerm-general
+#'
+#' @concept directed
+#' @concept frequently-used
+#' @concept layer-aware
+InitErgmTerm.mutualL<-function (nw, arglist, ...) {
+  ## Check the network and arguments to make sure they are appropriate.
+  a <- check.ErgmTerm(nw, arglist, directed=TRUE, bipartite=NULL,
+                      varnames = c("same", "by", "diff", "keep", "Ls"),
+                      vartypes = c("character", "character", "logical", "numeric", "formula,list"),
+                      defaultvalues = list(NULL, NULL, FALSE, NULL, NULL),
+                      required = c(FALSE, FALSE, FALSE, FALSE, FALSE))
+
+  assert_LHS_Layer(nw)
+
+  ## Process the arguments
+  if (!is.null(a$same) || !is.null(a$by)) {
+    if (!is.null(a$same)) {
+     attrname <- a$same
+     if (!is.null(a$by)) 
+       warning("Ignoring 'by' argument to mutual because 'same' exists", call.=FALSE)
+    }else{
+     attrname <- a$by
+    }
+    nodecov <- get.node.attr(nw, attrname)
+    u <- NVL(a$levels, sort(unique(nodecov)))
+    if (!is.null(a$keep)) {
+      u <- u[a$keep]
+    }
+    #   Recode to numeric
+    nodecov <- match(nodecov,u,nomatch=length(u)+1)
+    # All of the "nomatch" should be given unique IDs so they never match:
+    dontmatch <- nodecov==(length(u)+1)
+    nodecov[dontmatch] <- length(u) + (1:sum(dontmatch))
+    ui <- seq(along=u)
+  }
+
+  ### Construct the list to return
+  if (!is.null(a$same) || !is.null(a$by)) {
+    if (is.null(a$same)) {
+      coef.names <- paste("mutual.by", attrname, u, sep=".")
+      inputs <- c(ui, nodecov)
+    }else{
+     if (a$diff) {
+      coef.names <- paste("mutual.same", attrname, u, sep=".")
+      inputs <- c(ui, nodecov)
+     }else{ 
+      coef.names <- paste("mutual", attrname, sep=".")
+      inputs <- nodecov
+     }
+    }
+    if (is.null(a$same) && !is.null(a$by)) {
+     name <- "mutual_by_attr"
+    }else{
+     name <- "mutual"
+    }
+  }else{
+     name <- "mutual"
+     coef.names <- "mutual"
+     inputs <- NULL
+  }
+
+  maxval <- network.dyadcount(nw,FALSE)/2
+
+  Ls <- .set_layer_namemap(a$Ls, nw)
+  if(is(Ls,"formula")) Ls <- list(Ls)
+  L1 <- Ls[[1]]
+  L2 <- Ls[[2]]
+  if(!is.null(L1) || !is.null(L2)){
+    NVL(L1) <- L2
+    NVL(L2) <- L1
+    auxiliaries <- .mk_.layer.net_auxform(L1)
+    aux2 <- .mk_.layer.net_auxform(L2)
+    auxiliaries[[2]] <- call("+", auxiliaries[[2]], aux2[[2]])
+    name <- paste(name, "ML", sep="_")
+    coef.names <- .lspec_coef.namewrap(if(L1==L2) list(L1) else list(L1,L2))(coef.names)
+    maxval <- maxval*2
+  }else auxiliaries <- NULL
+  
+  list(name=name,                      #name: required
+       coef.names = coef.names,        #coef.names: required
+       inputs=inputs,
+       auxiliaries = auxiliaries,
+       minval = 0,
+       maxval = maxval) 
+}

R/InitErgmTerm.dgw_sp.R

@@ -19,10 +19,11 @@ InitErgmTerm..spcache.netL<-function(nw, arglist, ...){
 
   type <- match.arg(tolower(a$type), c("otp","osp","isp","utp")) # ITP not included, because it's just OTP with direction reversed.
 
-  if(is.directed(nw)==(type=="utp")) stop("Type UTP may only be used with undirected networks, the others only with directed.")
+  if (is.directed(nw) == (type == "utp")
+      && !(is.bipartite(nw) && type %in% c("osp", "isp")))
+    stop("Type UTP may only be used with undirected networks, OSP and ISP with bipartite or directed, and the rest only with directed.")
 
   dname <- paste0("_",type,"_wtnet")
-  
   linfo <- .sp.handle_layers(nw, a, type, FALSE)
 
 

R/InitErgmTerm.multilayer.R

@@ -234,7 +234,7 @@ combine_networks <- function(nwl, ignore.nattr=c("mnext"), ignore.vattr=c(), ign
   attrset <- if(keep.unshared.attr) union else intersect
 
   ns <- sapply(nwl, network.size)
-  es <- sapply(nwl, "%n%", "bipartite")
+  es <- sapply(nwl, b1.size)
   eblks <- c(0, cumsum(es))
   bip <- eblks[length(eblks)]
   ablks <- cumsum(c(bip, ns-es))
@@ -537,7 +537,7 @@ subnetwork_templates <- function(nw, split.vattr=nw%n%".blockID.vattr", names.va
 .block_vertexmap <- function(nw, by=nw %n% ".blockID.vattr", same_dim=FALSE){
   a <- .peek_vattrv(nw, by)
   n <- length(a)
-  bip <- nw %n% "bipartite"
+  bip <- b1.size(nw)
   if(NVL(bip,0)){
     ea <- a[seq_len(bip)]
     aa <- a[bip+seq_len(n-bip)]

R/InitErgmTerm.spcache.multilayer.R

@@ -0,0 +1,11 @@
+## These will eventually be moved to 'statnet.common'.
+
+dbl_along <- function(x) numeric(length(x))
+int_along <- function(x) integer(length(x))
+chr_along <- function(x) character(length(x))
+lgl_along <- function(x) logical(length(x))
+
+## This should be in 'network'.
+
+b1.size <- function(x) if(is.bipartite(x)) x %n% "bipartite" else FALSE
+b2.size <- function(x) if(is.bipartite(x)) network.size(x) - b1.size(x) else FALSE

R/combine.networks.R

@@ -0,0 +1,2 @@
+#' @note This term does not support multilayer networks with
+#'   heterogeneous bipartedness. This may change in the future.