From 0830084b38df17fa598aeb1370f5291a9915a8d1 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Thu, 28 Nov 2024 21:58:02 +1100 Subject: [PATCH] Implemented Projection() operator, .project.net() auxiliary, and Proj1() and Proj2() aliases to evaluate statistics on a network constructed by projecting a bipartite network. --- DESCRIPTION | 1 + R/InitErgmTerm.projection.R | 92 +++++++++++++++++ man/Project-ergmTerm-36962a06.Rd | 34 +++++++ src/changestats_projection.c | 153 +++++++++++++++++++++++++++++ tests/testthat/test-term-project.R | 29 ++++++ 5 files changed, 309 insertions(+) create mode 100644 R/InitErgmTerm.projection.R create mode 100644 man/Project-ergmTerm-36962a06.Rd create mode 100644 src/changestats_projection.c create mode 100644 tests/testthat/test-term-project.R diff --git a/DESCRIPTION b/DESCRIPTION index a89fa7bfb..fede22b0d 100644 --- a/DESCRIPTION +++ b/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' diff --git a/R/InitErgmTerm.projection.R b/R/InitErgmTerm.projection.R new file mode 100644 index 000000000..2c94af5cb --- /dev/null +++ b/R/InitErgmTerm.projection.R @@ -0,0 +1,92 @@ +# 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 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_error(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 Proj2-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_error(sQuote("mode"), " must be 1 or 2.") + list(name=paste0("_proj_net"), iinputs = mode, coef.names = c(), dependence=TRUE) +} diff --git a/man/Project-ergmTerm-36962a06.Rd b/man/Project-ergmTerm-36962a06.Rd new file mode 100644 index 000000000..8679dd9a5 --- /dev/null +++ b/man/Project-ergmTerm-36962a06.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/InitErgmTerm.projection.R +\name{Project-ergmTerm} +\alias{Project-ergmTerm} +\alias{InitErgmTerm.Project} +\alias{InitErgmTerm.Proj1} +\alias{Proj2-ergmTerm} +\alias{InitErgmTerm.Proj2} +\title{Evaluation on a projection of a bipartite network} +\usage{ +# binary: Project(formula, mode) + +# binary: Proj1(formula) + +# binary: Proj2(formula) +} +\arguments{ +\item{formula}{a one-sided \code{\link[=ergm]{ergm()}}-style formula with the terms to be evaluated} + +\item{mode}{the mode onto which to project: 1 or 2} +} +\description{ +This operator on a bipartite network evaluates the +formula on the undirected, valued network constructed by +projecting it onto its specified mode. \code{Proj1(formula)} and +\code{Proj2(formula)} are aliases for \code{Project(formula, 1)} and +\code{Project(formula, 2)}, respectively. +} +\seealso{ +\code{\link{ergmTerm}} for index of model terms currently visible to the package. + +\Sexpr[results=rd,stage=render]{ergm:::.formatTermKeywords("ergmTerm", "Project", "subsection")} +} +\concept{operator bipartite} diff --git a/src/changestats_projection.c b/src/changestats_projection.c new file mode 100644 index 000000000..6d7d32fee --- /dev/null +++ b/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; +} diff --git a/tests/testthat/test-term-project.R b/tests/testthat/test-term-project.R new file mode 100644 index 000000000..54c27ec3a --- /dev/null +++ b/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)) +})