Implemented Projection() operator, .project.net() auxiliary, and Proj…

…1() and Proj2() aliases to evaluate statistics on a network constructed by projecting a bipartite network.

DESCRIPTION

@@ -84,6 +84,7 @@ Collate:
     'InitErgmTerm.indices.R'
     'InitErgmTerm.interaction.R'
     'InitErgmTerm.operator.R'
+    'InitErgmTerm.projection.R'
     'InitErgmTerm.spcache.R'
     'InitErgmTerm.test.R'
     'InitErgmTerm.transitiveties.R'

R/InitErgmTerm.projection.R

@@ -0,0 +1,93 @@
+#  File R/InitErgmTerm.projection.R in package ergm, part of the
+#  Statnet suite of packages for network analysis, https://statnet.org .
+#
+#  This software is distributed under the GPL-3 license.  It is free,
+#  open source, and has the attribution requirements (GPL Section 7) at
+#  https://statnet.org/attribution .
+#
+#  Copyright 2003-2024 Statnet Commons
+################################################################################
+#' @templateVar name Project
+#' @title Evaluation on a projection of a bipartite network
+#'
+#' @description This operator on a bipartite network evaluates the
+#'   formula on the undirected, valued network constructed by
+#'   projecting it onto its specified mode. `Proj1(formula)` and
+#'   `Proj2(formula)` are aliases for `Project(formula, 1)` and
+#'   `Project(formula, 2)`, respectively.
+#'
+#' @usage
+#' # binary: Project(formula, mode)
+#' @template ergmTerm-formula
+#' @param mode the mode onto which to project: 1 or 2
+#'
+#' @template ergmTerm-general
+#'
+#' @concept operator
+#' @concept bipartite
+InitErgmTerm.Project <- function(nw, arglist, ...){
+  a <- check.ErgmTerm(nw, arglist, bipartite = TRUE,
+                      varnames = c("formula", "mode"),
+                      vartypes = c("formula", "numeric"),
+                      defaultvalues = list(NULL, NULL),
+                      required = c(TRUE, TRUE))
+
+  bip <- as.integer(nw %n% "bipartite")
+  n <- as.integer(network.size(nw))
+
+  mode <- a$mode
+  if(! mode %in% 1:2) ergm_Init_stop(sQuote("mode"), " must be 1 or 2.")
+
+  ### Construct an empty network with the correct structure.
+  pnw <- nw
+  pnw %n% "bipartite" <- FALSE
+  pnw[,] <- 0
+  pnw <- switch(mode,
+                get.inducedSubgraph(pnw, seq_len(bip)),
+                get.inducedSubgraph(pnw, bip + seq_len(n-bip)))
+  pnw %ergmlhs% "response" <- structure("w", valued = TRUE)
+
+  m <- ergm_model(a$formula, pnw, ..., offset.decorate=FALSE)
+  ergm_no_ext.encode(m)
+
+  auxiliaries <- trim_env(~.project.net(mode), "mode")
+
+  c(list(name="on_proj_net", iinputs = mode,
+         submodel = m, dependence = TRUE,
+         auxiliaries=auxiliaries),
+    wrap.ergm_model(m, pnw, ergm_mk_std_op_namewrap(paste0('Proj', mode))))
+}
+
+#' @templateVar name Project
+#' @template ergmTerm-rdname
+#' @aliases Proj1-ergmTerm
+#' @usage
+#' # binary: Proj1(formula)
+InitErgmTerm.Proj1 <- function(nw, arglist, ...){
+  arglist[["mode"]] <- 1L
+  f <- InitErgmTerm.Project
+  f(nw, arglist, ...)
+}
+
+#' @templateVar name Project
+#' @template ergmTerm-rdname
+#' @aliases Proj2-ergmTerm
+#' @usage
+#' # binary: Proj2(formula)
+InitErgmTerm.Proj2 <- function(nw, arglist, ...){
+  arglist[["mode"]] <- 2L
+  f <- InitErgmTerm.Project
+  f(nw, arglist, ...)
+}
+
+InitErgmTerm..project.net <- function(nw, arglist, ...){
+  a <- check.ErgmTerm(nw, arglist, bipartite = TRUE,
+                      varnames = c("mode"),
+                      vartypes = c("numeric"),
+                      defaultvalues = list(NULL),
+                      required = c(TRUE))
+
+  mode <- a$mode
+  if(! mode %in% 1:2) ergm_Init_stop(sQuote("mode"), " must be 1 or 2.")
+  list(name=paste0("_proj_net"), iinputs = mode, coef.names = c(), dependence=TRUE)
+}

src/changestats_projection.c

@@ -0,0 +1,153 @@
+/*  File src/changestats_projection.c in package ergm, part of the
+ *  Statnet suite of packages for network analysis, https://statnet.org .
+ *
+ *  This software is distributed under the GPL-3 license.  It is free,
+ *  open source, and has the attribution requirements (GPL Section 7) at
+ *  https://statnet.org/attribution .
+ *
+ *  Copyright 2003-2024 Statnet Commons
+ */
+#define STRICT_Wt_HEADERS
+#include "ergm_wtchangestat_operator.h"
+#include "ergm_changestat_operator.h"
+#include "ergm_storage.h"
+
+typedef struct {
+  WtModel *m;
+  Vertex *t, *h;
+  double *w;
+} StoreWtModelAndWtChanges;
+
+
+I_CHANGESTAT_FN(i_on_proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  ALLOC_STORAGE(1, StoreWtModelAndWtChanges, storage);
+  // No need to allocate it: we are only storing a pointer to a model.
+  storage->m = WtModelInitialize(getListElement(mtp->R, "submodel"), mtp->ext_state, pnwp, FALSE);
+
+  /* WtSELECT_C_OR_D_BASED_ON_SUBMODEL(storage->m); */
+  WtDELETE_IF_UNUSED_IN_SUBMODEL(x_func, storage->m);
+  WtDELETE_IF_UNUSED_IN_SUBMODEL(z_func, storage->m);
+
+  unsigned int mode = IINPUT_PARAM[0],
+    maxchanges = mode == 1 ? BIPARTITE-1 : N_NODES-BIPARTITE-1;
+
+  storage->t = R_Calloc(maxchanges, Vertex);
+  storage->h = R_Calloc(maxchanges, Vertex);
+  storage->w = R_Calloc(maxchanges, double);
+}
+
+/* TODO: A d_function? */
+
+C_CHANGESTAT_FN(c_on_proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  GET_STORAGE(StoreWtModelAndWtChanges, storage);
+
+  int echange = edgestate ? -1 : 1;
+
+  unsigned int nt = 0, mode = IINPUT_PARAM[0];
+
+  switch(mode){
+  case 1:
+    EXEC_THROUGH_FINEDGES(head, e, tail2, {
+        if(tail!=tail2){
+          storage->t[nt] = MIN(tail, tail2);
+          storage->h[nt] = MAX(tail, tail2);
+          storage->w[nt] = WtGETWT(tail, tail2, pnwp) + echange;
+          nt++;
+        }
+      }); break;
+  case 2:
+    EXEC_THROUGH_FOUTEDGES(tail, e, head2, {
+        if(head!=head2){
+          storage->t[nt] = MIN(head-BIPARTITE, head2-BIPARTITE);
+          storage->h[nt] = MAX(head-BIPARTITE, head2-BIPARTITE);
+          storage->w[nt] = WtGETWT(head-BIPARTITE, head2-BIPARTITE, pnwp) + echange;
+          nt++;
+        }
+      }); break;
+  default: error("We should never be here.");
+  }
+
+  WtChangeStats(nt, storage->t, storage->h, storage->w, pnwp, storage->m);
+
+  memcpy(CHANGE_STAT, storage->m->workspace, N_CHANGE_STATS*sizeof(double));
+}
+
+Z_CHANGESTAT_FN(z_on_proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  GET_STORAGE(StoreWtModelAndWtChanges, storage)
+
+  WtZStats(pnwp, storage->m, FALSE);
+
+  memcpy(CHANGE_STAT, storage->m->workspace, N_CHANGE_STATS*sizeof(double));
+}
+
+X_CHANGESTAT_FN(x_on_proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  GET_STORAGE(StoreWtModelAndWtChanges, storage);
+  ModelTerm *_mymtp = mtp;
+  WtSEND_X_SIGNAL_INTO(pnwp, storage->m, NULL, _mymtp->dstats, type, data);
+}
+
+F_CHANGESTAT_FN(f_on_proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  GET_STORAGE(StoreWtModelAndWtChanges, storage);
+
+  R_Free(storage->t);
+  R_Free(storage->h);
+  R_Free(storage->w);
+  WtModelDestroy(pnwp, storage->m);
+}
+
+
+
+I_CHANGESTAT_FN(i__proj_net){
+  unsigned int mode = IINPUT_PARAM[0];
+
+  WtNetwork *pnwp = AUX_STORAGE = WtNetworkInitialize(NULL, NULL, NULL, 0, mode == 1 ? BIPARTITE : N_NODES - BIPARTITE, DIRECTED, FALSE, FALSE, 0, NULL);
+
+  switch(mode){
+  case 1:
+    EXEC_THROUGH_NET_EDGES_PRE(tail, head, e1, {
+        EXEC_THROUGH_FINEDGES(head, e2, tail2, {
+            if(tail < tail2)
+              WtSETWT(tail, tail2, WtGETWT(tail, tail2, pnwp) + 1, pnwp);
+          });
+      }); break;
+  case 2:
+    EXEC_THROUGH_NET_EDGES_PRE(tail, head, e1, {
+        EXEC_THROUGH_FOUTEDGES(tail, e2, head2, {
+            if(head < head2)
+              WtSETWT(head-BIPARTITE, head2-BIPARTITE, WtGETWT(head-BIPARTITE, head2-BIPARTITE, pnwp) + 1, pnwp);
+          });
+      }); break;
+  default: error("We should never be here.");
+  }
+}
+
+U_CHANGESTAT_FN(u__proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  int echange = edgestate ? -1 : 1;
+  unsigned int mode = IINPUT_PARAM[0];
+
+  switch(mode){
+  case 1:
+    EXEC_THROUGH_FINEDGES(head, e, tail2, {
+        if(tail!=tail2)
+          WtSETWT(tail, tail2, WtGETWT(tail, tail2, pnwp) + echange, pnwp);
+      }); break;
+  case 2:
+    EXEC_THROUGH_FOUTEDGES(tail, e, head2, {
+        if(head!=head2)
+          WtSETWT(head-BIPARTITE, head2-BIPARTITE, WtGETWT(head-BIPARTITE, head2-BIPARTITE, pnwp) + echange, pnwp);
+      }); break;
+  default: error("We should never be here.");
+  }
+}
+
+F_CHANGESTAT_FN(f__proj_net){
+  GET_AUX_STORAGE(WtNetwork, pnwp);
+  WtNetworkDestroy(pnwp);
+  AUX_STORAGE = NULL;
+}

tests/testthat/test-term-project.R

@@ -0,0 +1,29 @@
+n <- 20
+b <- 5
+
+nw0 <- network.initialize(n, bipartite = b, directed = FALSE)
+nw0 %v% "b1" <- c(rnorm(b), rep(NA, n-b))
+nw0 %v% "b2" <- c(rep(NA, b), rnorm(n-b))
+nw1 <- simulate(nw0 ~ edges, coef = 0)
+nw2 <- simulate(nw0 ~ edges, coef = 0)
+m1 <- as.matrix(nw1)
+m2 <- as.matrix(nw2)
+
+nw1.absdiff.b1 <- sum(abs(outer(na.omit(nw1 %v% "b1"), na.omit(nw1 %v% "b1"), FUN = "-")) * tcrossprod(m1))/2
+nw1.absdiff.b2 <- sum(abs(outer(na.omit(nw1 %v% "b2"), na.omit(nw1 %v% "b2"), FUN = "-")) * crossprod(m1))/2
+nw2.absdiff.b1 <- sum(abs(outer(na.omit(nw2 %v% "b1"), na.omit(nw2 %v% "b1"), FUN = "-")) * tcrossprod(m2))/2
+nw2.absdiff.b2 <- sum(abs(outer(na.omit(nw2 %v% "b2"), na.omit(nw2 %v% "b2"), FUN = "-")) * crossprod(m2))/2
+
+f <- nw1 ~ Project(~ absdiff("b1"), 1) + Proj1(~absdiff("b1")) + Project(~absdiff("b2"), 2) + Proj2(~absdiff("b2"))
+
+test_that("Projection(), Proj1(), and Proj2() summary", {
+  expect_equal(summary(f),
+               c(`Proj1~absdiff.sum.b1` = nw1.absdiff.b1, `Proj1~absdiff.sum.b1` = nw1.absdiff.b1,
+                 `Proj2~absdiff.sum.b2` = nw1.absdiff.b2, `Proj2~absdiff.sum.b2` = nw1.absdiff.b2))
+})
+
+test_that("Projection(), Proj1(), and Proj2() change", {
+  changes <- as.matrix(xor(nw1, nw2), matrix.type = "edgelist")
+  expect_equal(c(ergm.godfather(f, changes = changes)),
+               rep(c(nw2.absdiff.b1, nw2.absdiff.b2), each = 2))
+})