From 2c322bcf53c40621d39357fbe72d79bbb142d6e2 Mon Sep 17 00:00:00 2001 From: "Carter T. Butts" Date: Mon, 27 Jul 2026 02:16:21 -0700 Subject: [PATCH 01/18] These changes add the distance term to ergm; the package passes R CMD check, the term passes its checks, and the examples run correctly. A NEWS item has been added, and the term is also documented (and a test has been added). --- R/InitErgmTerm.R | 239 ++++++++++++++++++++++++++++++ inst/NEWS.Rd | 3 + man/distance-ergmTerm-9bc71012.Rd | 197 ++++++++++++++++++++++++ src/changestats.c | 85 +++++++++++ tests/testthat/test-distance.R | 108 ++++++++++++++ 5 files changed, 632 insertions(+) create mode 100644 man/distance-ergmTerm-9bc71012.Rd create mode 100644 tests/testthat/test-distance.R diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index 5d927d602..ad65f6aa4 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -25,6 +25,7 @@ # C: = # D: # +# # E: # G: # @@ -2412,6 +2413,244 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) ) } +################################################################################ + +#' @templateVar name distance +#' @title Inter-point distances +#' @description This term adds a single statistic to the model whose value is +#' the sum over all edge variables of the edge variable value times either +#' the distance between the respective vertices, or its log (if +#' \code{log==TRUE}). The \code{coord} agument must contain a vector, +#' matrix, or \code{data.frame} of coordinates (vertex x dimension), or +#' else the name of a network attribute with such coordinates. +#' +#' @usage +#' # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +#' # log=TRUE, mindist=1e-5, distoff=0) +#' +#' @param coord node by dimension coordinate matrix, or name of a network +#' attribute containing the coordinate matrix; for spherical +#' coordinates, must contain two columns of angular coordinates +#' in lat/lon order +#' @param metric power to use for the Minkowski metric +#' @param sphere logical; should great circle distances be used (rather than +#' Minkowski distances in free space)? +#' @param radius for spherical distances, the radius of the sphere to use; +#' defaults to the IUGG mean Earth radius, in km +#' @param log logical; use log rather than raw distances? +#' @param mindist for log distances, a lower threshold; for dyads such that +#' the observed distance plus any offset is less than the +#' minimum, the log of the minimum is used instead. Ignored for +#' raw distances +#' @param distoff for log distances, an offset to be added to observed +#' distances before taking the logarithm +#' +#' @details Either spherical (\code{sphere=TRUE}) or Minkowski metrics +#' (\code{sphere=FALSE}) may be selected. For the latter, any number of +#' dimensions may be supplied, but if \code{coord} is given as a vector the +#' space in question is assumed to be one-dimensional. The choice of +#' Minkowski metric is determined by \code{metric}, with the distance being +#' given by +#' \deqn{ +#' D(i,j) = (sum_d |x_{id}-x_{jd}|^p)^(1/p) +#' } +#' where the sum is over the dimensions of the space, \eqn{x} is the +#' coordinate matrix, and \eqn{p} is the metric parameter (i.e. +#' \code{metric}). The default value of \code{metric} is 2, yielding the +#' Euclidean distance. Note that \code{metric==1}, and \code{log==FALSE} +#' in the one-dimensional case is identical to the \code{absdiff} term. +#' +#' When \code{sphere=TRUE}, \code{coord} must be a two-dimensional matrix of +#' angular coordinates in lat/lon form (i.e., the first column must contain +#' units of decimal degrees between -90 and 90, and the second must contain +#' units of decimal degrees betweein -180 and 180). Distance is then +#' computed on between the specified angular coordinates on the surface of a +#' sphere of radius \code{radius}; by default, this is the IUGG mean Earth +#' radius in km, and hence supplying lat/lon coordinates yields geospherical +#' distances in kilometers. As \code{metric} has no meaning here, it is +#' ignored. +#' +#' When \code{log==TRUE}, the logarithm of the distance rather than the +#' distance itself is used to form the statistic. This produces interaction +#' functions that roughly approximate the power law spatial interaction +#' functions (SIFs) from Butts and Acton (2011). To prevent divergence of +#' the log distance for points at identical positions, \code{distoff} can be +#' used to provide an offset to the raw distance (with the computed value +#' becoming \code{log(d+distoff)}. (This can have an effect that is similar +#' to an \dQuote{attenuated power law} SIF, per Butts and Acton.) As an +#' additional failsafe, the (offset) distances are thresholded from below by +#' \code{mindist}. Setting both \code{mindist} and \code{distoff} to zero is +#' allowed, but may produce exciting results if points precisely overlap. +#' Note that both \code{distoff} and \code{mindist} are ignored when +#' \code{log==FALSE}. +#' +#' This term can be used for directed or undirected networks. +#' +#' @examples +#' \dontrun{ #We turn these off by default, b/c they are a bit slow +#' #Create an example network, in a two-dimensional space +#' # Effective SIF is 1/(1 + exp(-3 + 2 log(d))); apx inverse square +#' n <- 300 +#' d <- 2 +#' x <- matrix(runif(d*n,0,50),ncol=d) +#' net <- simulate(network.initialize(n,directed=FALSE)~edges +#' + distance(x), coef=c(3,-2), +#' control=control.simulate.formula(MCMC.burnin=n^3)) +#' +#' #Examine the network +#' plot(net,coord=x) #At first, seems not so spatial.... +#' dis <- as.matrix(dist(x)) +#' plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]), +#' xlab="Distance", ylab="Edge State") #But actually quite spatial! +#' lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF... +#' y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2) +#' lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))), +#' lwd=2, lty=3, col=3) #Theoretical SIF +#' +#' #Recover the parameters +#' summary(ergm(net~edges+distance(x))) +#' +#' #What happens if we use raw instead of log distances? +#' # Effective SIF is 1/(1 + exp(-3 + d)); apx exponential decay +#' net <- simulate(network.initialize(n,directed=FALSE)~edges +#' + distance(x, log=FALSE), coef=c(3,-1), +#' control=control.simulate.formula(MCMC.burnin=n^3)) +#' plot(net,coord=x) #Ties are much more local - no long edges +#' +#' plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]), +#' xlab="Distance", ylab="Edge State") +#' lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF... +#' y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2) +#' lines((0:100)/100*max(dis), 1/(1+exp(-3+1*((0:100)/100*max(dis)))), +#' lwd=2, lty=3, col=3) #New SIF (green) +#' lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))), +#' lwd=2, lty=3, col=4) #Previous SIF (blue) +#' +#' summary(ergm(net~edges+distance(x, log=FALSE))) #Recover parameters +#' +#' #A small lat/lon example; begin with some US cities +#' co<-rbind( +#' AustinTX=c(30.26694,-97.74278), +#' BaltimoreMD=c(39.29028,-76.6125), +#' BatonRougeLA=c(30.45056,-91.15444), +#' CharlotteNC=c(35.22694,-80.84333), +#' ChicagoIL=c(41.85,-87.65), +#' ColumbusOH=c(39.96111,-82.99889), +#' DetroitMI=c(42.33139,-83.04583), +#' DurhamNC=c(35.99389,-78.89889), +#' HoustonTX=c(29.76306,-95.36306), +#' IndianapolisIN=c(39.76833,-86.15806), +#' IrvineCA=c(33.66944,-117.8222), +#' IthacaNY=c(42.44056,-76.49694), +#' JacksonvilleFL=c(30.33194,-81.65583), +#' KonaHI=c(19.64083,-155.9856), +#' LasVegasNV=c(36.175,-115.1364), +#' LosAngelesCA=c(34.05222,-118.2428), +#' MemphisTN=c(35.14944,-90.04889), +#' MiamiFL=c(25.77389,-80.19389), +#' MilwaukeeWI=c(43.03889,-87.90639), +#' MobileAL=c(30.69417,-88.04306), +#' NewOrleansLA=c(29.95444,-90.075), +#' NewYorkNY=c(40.71417,-74.00639), +#' PhiladelphiaPA=c(39.95222,-75.16417), +#' PhoenixAZ=c(33.44833,-112.0733), +#' PittsburghPA=c(40.39556,-79.83889), +#' SanAntonioTX=c(29.42389,-98.49333), +#' SanDiegoCA=c(32.71528,-117.1564), +#' SanFranciscoCA=c(37.775,-122.4183), +#' SanJoseCA=c(37.33944,-121.8939), +#' SeattleWA=c(47.60639,-122.3308) +#' ) +#' +#' #By default, spherical distances give us great circle distances +#' #on the geosphere, in kilometers +#' net <- network.initialize(30, directed = FALSE) +#' net[1,8] <- 1 #Create a tie from Austin, TX to Durham, NC +#' summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 1863 km +#' net[1,8] <- 0 +#' net[14,22] <- 1 #Now try Kona, HI to NYC, NY +#' summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 7940 km +#' +#' #Model a network among these fine cities +#' net <- simulate(net ~ edges + distance(co, sphere = TRUE), +#' coef=c(5,-1), control=control.simulate.formula(MCMC.burnin=1e5)) +#' +#' #Plot the cities in lat/lon space +#' plot(net, coord=co[,2:1], suppress.axes=FALSE, xlab="Longitude", +#' ylab="Latitude") +#' +#' #Recover the parameters +#' summary(ergm(net~edges+distance(co, sphere = TRUE))) +#'} +#' +#' @references Butts, Carter T. and Acton, Ryan M. (2011). \dQuote{Spatial +#' Modeling of Social Networks.} In Timothy Nyerges, Helen Couclelis, and +#' Robert McMaster (Eds.), \emph{The Sage Handbook of GIS and Society +#' Research}, 222--250. SAGE Publications. +#' +#' @template ergmTerm-general +#' +#' @concept dyad-independent +#' @concept directed +#' @concept undirected +#' @concept quantitative nodal attribute +InitErgmTerm.distance <- function(nw, arglist, ...) { + a <- check.ErgmTerm(nw, arglist, + varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff"), + vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric"), + required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), + defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0)) + #Process arguments + if(is.character(a$coord)){ + a$coord<-get.network.attribute(nw,a$coord) + if(is.null(a$coord)) + stop("Distance term requires either a coordinate matrix, or the name of a network attribute containing one.") + } + if(length(dim(a$coord))==2){ + coord<-a$coord + if(is.data.frame(coord)) + coord<-as.matrix(coord) + }else if(length(dim(a$coord))==0){ + coord<-matrix(rep(a$coord,length=network.size(nw)),ncol=1) + }else{ + coord<-as.matrix(coord) + } + if(NROW(coord)!=network.size(nw)) + stop("Distance term requires that coordinates be provided for all vertices.") + if(any(!apply(coord,1:2,is.numeric))) + stop("Distance term requires numeric coordinates.") + if(a$sphere){ + if(NCOL(coord)!=2) + stop("Spherical coordinates chosen for distance term, but more or fewer than two dimensions provided; coordinates must be lat,lon or equivalent (in that order, and in angular units).") + if(any(coord[,1]< -90)||any(coord[,1]>90)) + stop("Illegal latitude passed to distance term.") + if(any(coord[,2]< -180)||any(coord[,2]>180)) + stop("Illegal longitude passed to distance term.") + if(a$radius<=0) + stop("Radius of sphere to be used for distance calculations must be positive.") + } + if(a$log&&(a$mindist<0)) + stop("mindist cannot be negative in distance term.") + if(a$log&&(a$distoff<0)) + stop("distoff cannot be negative in distance term.") + #Set up the call + if(a$log) + logstr<-".log" + else + logstr<-"" + if(a$sphere) + basestr<-"dist" + else + basestr<-paste0("dist.L",a$metric) + list(name = "distance", + coef.names = paste0(basestr,logstr), + inputs = c(is.directed(nw), network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, coord), + dependence = FALSE, + emptynwstats = 0 + ) +} + + ################################################################################ #' @templateVar name dyadcov diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 7210c5bba..c3f711e37 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -79,6 +79,9 @@ \subsection{NEW FEATURES}{ \itemize{ + \item{ + New term \ergmTerm{ergm}{distance}{()} for Minkowski and great circle distances between vertices now available. + } \item{ \ergmTerm{ergm}{nodemix}{()} and \ergmTerm{ergm}{mm}{()} are now more user-friendly in how \code{levels2=} arguments are processed. } diff --git a/man/distance-ergmTerm-9bc71012.Rd b/man/distance-ergmTerm-9bc71012.Rd new file mode 100644 index 000000000..b3346b3ad --- /dev/null +++ b/man/distance-ergmTerm-9bc71012.Rd @@ -0,0 +1,197 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/InitErgmTerm.R +\name{distance-ergmTerm} +\alias{distance-ergmTerm} +\alias{InitErgmTerm.distance} +\title{Inter-point distances} +\usage{ +# binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +# log=TRUE, mindist=1e-5, distoff=0) +} +\arguments{ +\item{coord}{node by dimension coordinate matrix, or name of a network +attribute containing the coordinate matrix; for spherical +coordinates, must contain two columns of angular coordinates +in lat/lon order} + +\item{metric}{power to use for the Minkowski metric} + +\item{sphere}{logical; should great circle distances be used (rather than +Minkowski distances in free space)?} + +\item{radius}{for spherical distances, the radius of the sphere to use; +defaults to the IUGG mean Earth radius, in km} + +\item{log}{logical; use log rather than raw distances?} + +\item{mindist}{for log distances, a lower threshold; for dyads such that +the observed distance plus any offset is less than the +minimum, the log of the minimum is used instead. Ignored for +raw distances} + +\item{distoff}{for log distances, an offset to be added to observed +distances before taking the logarithm} +} +\description{ +This term adds a single statistic to the model whose value is +the sum over all edge variables of the edge variable value times either +the distance between the respective vertices, or its log (if +\code{log==TRUE}). The \code{coord} agument must contain a vector, +matrix, or \code{data.frame} of coordinates (vertex x dimension), or +else the name of a network attribute with such coordinates. +} +\details{ +Either spherical (\code{sphere=TRUE}) or Minkowski metrics +(\code{sphere=FALSE}) may be selected. For the latter, any number of +dimensions may be supplied, but if \code{coord} is given as a vector the +space in question is assumed to be one-dimensional. The choice of +Minkowski metric is determined by \code{metric}, with the distance being +given by +\deqn{ + D(i,j) = (sum_d |x_{id}-x_{jd}|^p)^(1/p) + } +where the sum is over the dimensions of the space, \eqn{x} is the +coordinate matrix, and \eqn{p} is the metric parameter (i.e. +\code{metric}). The default value of \code{metric} is 2, yielding the +Euclidean distance. Note that \code{metric==1}, and \code{log==FALSE} +in the one-dimensional case is identical to the \code{absdiff} term. + +When \code{sphere=TRUE}, \code{coord} must be a two-dimensional matrix of +angular coordinates in lat/lon form (i.e., the first column must contain +units of decimal degrees between -90 and 90, and the second must contain +units of decimal degrees betweein -180 and 180). Distance is then +computed on between the specified angular coordinates on the surface of a +sphere of radius \code{radius}; by default, this is the IUGG mean Earth +radius in km, and hence supplying lat/lon coordinates yields geospherical +distances in kilometers. As \code{metric} has no meaning here, it is +ignored. + +When \code{log==TRUE}, the logarithm of the distance rather than the +distance itself is used to form the statistic. This produces interaction +functions that roughly approximate the power law spatial interaction +functions (SIFs) from Butts and Acton (2011). To prevent divergence of +the log distance for points at identical positions, \code{distoff} can be +used to provide an offset to the raw distance (with the computed value +becoming \code{log(d+distoff)}. (This can have an effect that is similar +to an \dQuote{attenuated power law} SIF, per Butts and Acton.) As an +additional failsafe, the (offset) distances are thresholded from below by +\code{mindist}. Setting both \code{mindist} and \code{distoff} to zero is +allowed, but may produce exciting results if points precisely overlap. +Note that both \code{distoff} and \code{mindist} are ignored when +\code{log==FALSE}. + +This term can be used for directed or undirected networks. +} +\examples{ +\dontrun{ #We turn these off by default, b/c they are a bit slow +#Create an example network, in a two-dimensional space +# Effective SIF is 1/(1 + exp(-3 + 2 log(d))); apx inverse square +n <- 300 +d <- 2 +x <- matrix(runif(d*n,0,50),ncol=d) +net <- simulate(network.initialize(n,directed=FALSE)~edges + + distance(x), coef=c(3,-2), + control=control.simulate.formula(MCMC.burnin=n^3)) + +#Examine the network +plot(net,coord=x) #At first, seems not so spatial.... +dis <- as.matrix(dist(x)) +plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]), + xlab="Distance", ylab="Edge State") #But actually quite spatial! +lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF... + y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2) +lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))), + lwd=2, lty=3, col=3) #Theoretical SIF + +#Recover the parameters +summary(ergm(net~edges+distance(x))) + +#What happens if we use raw instead of log distances? +# Effective SIF is 1/(1 + exp(-3 + d)); apx exponential decay +net <- simulate(network.initialize(n,directed=FALSE)~edges + + distance(x, log=FALSE), coef=c(3,-1), + control=control.simulate.formula(MCMC.burnin=n^3)) +plot(net,coord=x) #Ties are much more local - no long edges + +plot(dis[upper.tri(dis)], jitter(as.matrix(net)[upper.tri(dis)]), + xlab="Distance", ylab="Edge State") +lines(smooth.spline(x=dis[upper.tri(dis)], #Empirical SIF... + y=as.matrix(net)[upper.tri(dis)]), col=2, lwd=2) +lines((0:100)/100*max(dis), 1/(1+exp(-3+1*((0:100)/100*max(dis)))), + lwd=2, lty=3, col=3) #New SIF (green) +lines((0:100)/100*max(dis), 1/(1+exp(-3+2*log((0:100)/100*max(dis)))), + lwd=2, lty=3, col=4) #Previous SIF (blue) + +summary(ergm(net~edges+distance(x, log=FALSE))) #Recover parameters + +#A small lat/lon example; begin with some US cities +co<-rbind( + AustinTX=c(30.26694,-97.74278), + BaltimoreMD=c(39.29028,-76.6125), + BatonRougeLA=c(30.45056,-91.15444), + CharlotteNC=c(35.22694,-80.84333), + ChicagoIL=c(41.85,-87.65), + ColumbusOH=c(39.96111,-82.99889), + DetroitMI=c(42.33139,-83.04583), + DurhamNC=c(35.99389,-78.89889), + HoustonTX=c(29.76306,-95.36306), + IndianapolisIN=c(39.76833,-86.15806), + IrvineCA=c(33.66944,-117.8222), + IthacaNY=c(42.44056,-76.49694), + JacksonvilleFL=c(30.33194,-81.65583), + KonaHI=c(19.64083,-155.9856), + LasVegasNV=c(36.175,-115.1364), + LosAngelesCA=c(34.05222,-118.2428), + MemphisTN=c(35.14944,-90.04889), + MiamiFL=c(25.77389,-80.19389), + MilwaukeeWI=c(43.03889,-87.90639), + MobileAL=c(30.69417,-88.04306), + NewOrleansLA=c(29.95444,-90.075), + NewYorkNY=c(40.71417,-74.00639), + PhiladelphiaPA=c(39.95222,-75.16417), + PhoenixAZ=c(33.44833,-112.0733), + PittsburghPA=c(40.39556,-79.83889), + SanAntonioTX=c(29.42389,-98.49333), + SanDiegoCA=c(32.71528,-117.1564), + SanFranciscoCA=c(37.775,-122.4183), + SanJoseCA=c(37.33944,-121.8939), + SeattleWA=c(47.60639,-122.3308) +) + +#By default, spherical distances give us great circle distances +#on the geosphere, in kilometers +net <- network.initialize(30, directed = FALSE) +net[1,8] <- 1 #Create a tie from Austin, TX to Durham, NC +summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 1863 km +net[1,8] <- 0 +net[14,22] <- 1 #Now try Kona, HI to NYC, NY +summary(net ~ distance(co, sphere = TRUE, log = FALSE)) #About 7940 km + +#Model a network among these fine cities +net <- simulate(net ~ edges + distance(co, sphere = TRUE), + coef=c(5,-1), control=control.simulate.formula(MCMC.burnin=1e5)) + +#Plot the cities in lat/lon space +plot(net, coord=co[,2:1], suppress.axes=FALSE, xlab="Longitude", + ylab="Latitude") + +#Recover the parameters +summary(ergm(net~edges+distance(co, sphere = TRUE))) +} + +} +\references{ +Butts, Carter T. and Acton, Ryan M. (2011). \dQuote{Spatial +Modeling of Social Networks.} In Timothy Nyerges, Helen Couclelis, and +Robert McMaster (Eds.), \emph{The Sage Handbook of GIS and Society +Research}, 222--250. SAGE Publications. +} +\seealso{ +\code{\link{ergmTerm}} for index of model terms currently visible to the package. + +\Sexpr[results=rd,stage=render]{ergm:::.formatTermKeywords("ergmTerm", "distance", "subsection")} +} +\concept{directed} +\concept{dyad-independent} +\concept{quantitative nodal attribute} +\concept{undirected} diff --git a/src/changestats.c b/src/changestats.c index ee0ba70d8..c3a1a9c73 100644 --- a/src/changestats.c +++ b/src/changestats.c @@ -1463,6 +1463,91 @@ C_CHANGESTAT_FN(c_degree_w_homophily) { } } + +/***************** + changestat: d_distance +*****************/ +/* +Utility function to calculate distances on a sphere of radius r, between two points +in spherical coordinates. Our notation betrays our assumption that these are +lat/lon points on the geosphere, but technically you could use any sphere you +wanted, by passing an alternative radius (so long as your coordinates were given +in angular units). +*/ +double spheredist(double lat0, double lon0, double lat1, double lon1, double r) { + double rlat0,rlat1,rlon0,rlon1,cosrl0,cosrl1,cosdl,sinrl0,sinrl1,sindl; + + rlat0=lat0/180.0*M_PI; + rlat1=lat1/180.0*M_PI; + rlon0=lon0/180.0*M_PI; + rlon1=lon1/180.0*M_PI; + cosrl0=cos(rlat0); + cosrl1=cos(rlat1); + cosdl=cos(rlon0-rlon1); + sinrl0=sin(rlat0); + sinrl1=sin(rlat1); + sindl=sin(rlon0-rlon1); + return r*atan2(sqrt((cosrl0*sindl) * (cosrl0*sindl) + (cosrl1*sinrl0-sinrl1*cosrl0*cosdl) * (cosrl1*sinrl0-sinrl1*cosrl0*cosdl)), sinrl0*sinrl1+cosrl0*cosrl1*cosdl); +} + +/*Actual distance changestat function itself (calls spheredist)*/ +CHANGESTAT_FN(d_distance) { + Vertex t, h; + int i,j; + int isdir,nv,dim,logd,sphd; + double dis,dexp,*coord,sphr,logmind,logdoff; + + /*Set things up*/ + isdir = (int)INPUT_PARAM[0]; /*Is the graph directed?*/ + nv = (int)INPUT_PARAM[1]; /*Number of vertices*/ + dim = (int)INPUT_PARAM[2]; /*Dimension of coordinate space*/ + dexp = INPUT_PARAM[3]; /*Minkowski exponent*/ + logd = (int)INPUT_PARAM[4]; /*Should we use the log distance?*/ + sphd = (int)INPUT_PARAM[5]; /*Should we use lat/lon distances on the geosphere?*/ + sphr = INPUT_PARAM[6]; /*Radius to use for spherical coordinates*/ + logmind = INPUT_PARAM[7]; /*Distance minimum for the log case*/ + logdoff = INPUT_PARAM[8]; /*Distance offset for the log case*/ + coord = INPUT_PARAM+9; /*Pointer to the coordinate matrix*/ + + /*Compute the changescores*/ + ZERO_ALL_CHANGESTATS(i); + FOR_EACH_TOGGLE(i) { + t = TAIL(i); h = HEAD(i); + if(sphd){ /*Compute distances on the sphere*/ + dis=spheredist(coord[t-1],coord[t-1+nv],coord[h-1],coord[h-1+nv],sphr); + }else{ /*Compute Minkowski distances in free space*/ + dis=0.0; + for(j=0;j0]), ignore_attr=TRUE) +}) +test_that("L1 works", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE)), sum(d1[g>0]), ignore_attr=TRUE) +}) +test_that("L2 works", { + expect_equal(summary(net~distance(co)), sum(log(d2[g>0])), ignore_attr=TRUE) +}) +test_that("log L1 works", { + expect_equal(summary(net~distance(co,metric=1)), sum(log(d1[g>0])), ignore_attr=TRUE) +}) + +#Verify that alternative specifications work +test_that("character matches matrix", { + expect_equal(summary(net~distance(co)), summary(net~distance("co")), ignore_attr=TRUE) +}) +test_that("character matches data.frame", { + expect_equal(summary(net~distance(as.data.frame(co))), summary(net~distance("co")), ignore_attr=TRUE) +}) + +#Verify that values are correct +summary(net~distance(co)) + +#Test for lat/lon calculations +co<-rbind( + AustinTX=c(30.26694,-97.74278), + BaltimoreMD=c(39.29028,-76.6125), + BatonRougeLA=c(30.45056,-91.15444), + CharlotteNC=c(35.22694,-80.84333), + ChicagoIL=c(41.85,-87.65), + ColumbusOH=c(39.96111,-82.99889), + DetroitMI=c(42.33139,-83.04583), + DurhamNC=c(35.99389,-78.89889), + HoustonTX=c(29.76306,-95.36306), + IndianapolisIN=c(39.76833,-86.15806), + IrvineCA=c(33.66944,-117.8222), + IthacaNY=c(42.44056,-76.49694), + JacksonvilleFL=c(30.33194,-81.65583), + KonaHI=c(19.64083,-155.9856), + LasVegasNV=c(36.175,-115.1364), + LosAngelesCA=c(34.05222,-118.2428), + MemphisTN=c(35.14944,-90.04889), + MiamiFL=c(25.77389,-80.19389), + MilwaukeeWI=c(43.03889,-87.90639), + MobileAL=c(30.69417,-88.04306), + NewOrleansLA=c(29.95444,-90.075), + NewYorkNY=c(40.71417,-74.00639), + PhiladelphiaPA=c(39.95222,-75.16417), + PhoenixAZ=c(33.44833,-112.0733), + PittsburghPA=c(40.39556,-79.83889), + SanAntonioTX=c(29.42389,-98.49333), + SanDiegoCA=c(32.71528,-117.1564), + SanFranciscoCA=c(37.775,-122.4183), + SanJoseCA=c(37.33944,-121.8939), + SeattleWA=c(47.60639,-122.3308) +) + +#By default, spherical distances give us great circle distances +#on the geosphere, in kilometers +net <- network.initialize(30, directed = FALSE) +net[1,8] <- 1 #Create a tie from Austin, TX to Durham, NC +test_that("first great circle distance is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 1862.879, ignore_attr=TRUE) +}) +net[1,8] <- 0 +net[14,22] <- 1 #Now try Kona, HI to NYC, NY +test_that("second great circle distance is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 7940.019, ignore_attr=TRUE) +}) +net[1,8] <- 1 +test_that("great circle distance sum is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = FALSE)),3), 9802.898, ignore_attr=TRUE) +}) +test_that("log great circle distance sum is correct", { + expect_equal(round(summary(net ~ distance(co, sphere = TRUE, log = TRUE)),3), 16.51, ignore_attr=TRUE) +}) From f90cf9ed8eb72d1d01f6a33331df32f392d363da Mon Sep 17 00:00:00 2001 From: "Carter T. Butts" Date: Tue, 28 Jul 2026 02:27:49 -0700 Subject: [PATCH 02/18] Changed the changestat function for distance from d_ to c_, added new functionality (including vertex based attributes, scaling, and the ability to combine spherical and Minkowski metrics), updated docs and tests, and tested everything. Seems to be working fine. --- R/InitErgmTerm.R | 91 ++++++++++++++++++++------- man/distance-ergmTerm-9bc71012.Rd | 57 +++++++++++++---- src/changestats.c | 100 +++++++++++++++--------------- tests/testthat/test-distance.R | 20 +++++- 4 files changed, 181 insertions(+), 87 deletions(-) diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index ad65f6aa4..932d6119e 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -2421,20 +2421,24 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' the sum over all edge variables of the edge variable value times either #' the distance between the respective vertices, or its log (if #' \code{log==TRUE}). The \code{coord} agument must contain a vector, -#' matrix, or \code{data.frame} of coordinates (vertex x dimension), or -#' else the name of a network attribute with such coordinates. +#' matrix, or \code{data.frame} of coordinates (vertex x dimension), the +#' name of a network attribute with such coordinates, or a vector of vertex +#' attributes to be used as coordinates. #' #' @usage #' # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, -#' # log=TRUE, mindist=1e-5, distoff=0) +#' # log=TRUE, mindist=1e-5, distoff=0, scale=1) #' -#' @param coord node by dimension coordinate matrix, or name of a network -#' attribute containing the coordinate matrix; for spherical -#' coordinates, must contain two columns of angular coordinates -#' in lat/lon order +#' @param coord node by dimension coordinate matrix, name of a network +#' attribute containing the coordinate matrix, or vector of +#' vertex attribute names containing coordinates (to be used +#' jointly); for \code{sphere==TRUE}, first two coordinates +#' must be angular units, in lat/lon order #' @param metric power to use for the Minkowski metric #' @param sphere logical; should great circle distances be used (rather than -#' Minkowski distances in free space)? +#' Minkowski distances in free space)? If \code{TRUE}, the first +#' two coordinates are respectively interpreted as latitude and +#' longitude coordinates on a sphere of specified radius #' @param radius for spherical distances, the radius of the sphere to use; #' defaults to the IUGG mean Earth radius, in km #' @param log logical; use log rather than raw distances? @@ -2444,6 +2448,8 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' raw distances #' @param distoff for log distances, an offset to be added to observed #' distances before taking the logarithm +#' @param scale factor by which raw distances should be rescaled prior to other +#' operations (notably, offsetting, logging, and thresholding) #' #' @details Either spherical (\code{sphere=TRUE}) or Minkowski metrics #' (\code{sphere=FALSE}) may be selected. For the latter, any number of @@ -2453,6 +2459,8 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' given by #' \deqn{ #' D(i,j) = (sum_d |x_{id}-x_{jd}|^p)^(1/p) +#' }{ +#' D(i,j) = \left(\sum_d |x_{id}-x_{jd}|^p\right)^{1/p} #' } #' where the sum is over the dimensions of the space, \eqn{x} is the #' coordinate matrix, and \eqn{p} is the metric parameter (i.e. @@ -2460,7 +2468,7 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' Euclidean distance. Note that \code{metric==1}, and \code{log==FALSE} #' in the one-dimensional case is identical to the \code{absdiff} term. #' -#' When \code{sphere=TRUE}, \code{coord} must be a two-dimensional matrix of +#' When \code{sphere=TRUE}, the first two dimensions of \code{coord} must be #' angular coordinates in lat/lon form (i.e., the first column must contain #' units of decimal degrees between -90 and 90, and the second must contain #' units of decimal degrees betweein -180 and 180). Distance is then @@ -2470,6 +2478,29 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' distances in kilometers. As \code{metric} has no meaning here, it is #' ignored. #' +#' If \code{sphere=TRUE} and more than two coordinates are provided per +#' vertex, then the first two are used to compute great circle distances, +#' with the remainder being used to compute the selected Minkowski metric. +#' These two distances are then added to form the final distance. While +#' this has exotic uses, the most obvious is to add elevation. It is +#' important to ensure that the coordinates are specified such that the +#' two distances have the same units! (By default, spherical distances are +#' in kilometers, so unless \code{radius} is changed, the user should ensure +#' that additional coordinates are also in kilometers.) Note that elevation +#' effects specified in this way do not take into account factors such as the +#' need to go down to the street and then up to a target floor within another +#' building (as in the case of artificial elevation in cities), nor the need +#' to e.g. get up before you get down in hilly terrain. Life is full of such +#' disappointments. +#' +#' The overall scale of the computed distances can be modified by setting the +#' \code{scale} argument. This is most often useful for numerical or +#' interpretational reasons (e.g., to avoid overflow/underflow issues, or to +#' change to familiar units). However, it can also be used in conjunction +#' with offsets and log scaling to change the form of the distance effect, +#' as described below. In particular, note that rescaling is performed prior +#' to application of offsets or thresholding, where applicable. +#' #' When \code{log==TRUE}, the logarithm of the distance rather than the #' distance itself is used to form the statistic. This produces interaction #' functions that roughly approximate the power law spatial interaction @@ -2479,10 +2510,12 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' becoming \code{log(d+distoff)}. (This can have an effect that is similar #' to an \dQuote{attenuated power law} SIF, per Butts and Acton.) As an #' additional failsafe, the (offset) distances are thresholded from below by -#' \code{mindist}. Setting both \code{mindist} and \code{distoff} to zero is +#' \code{mindist}. With all elements considered, the net effect is thus +#' \code{log(max(mindist, distoff + scale*d))}, where \eqn{d} is the raw +#' distance. Setting both \code{mindist} and \code{distoff} to zero is #' allowed, but may produce exciting results if points precisely overlap. #' Note that both \code{distoff} and \code{mindist} are ignored when -#' \code{log==FALSE}. +#' \code{log==FALSE} (but \code{scale} is not). #' #' This term can be used for directed or undirected networks. #' @@ -2593,18 +2626,26 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' @concept dyad-independent #' @concept directed #' @concept undirected +#' @concept bipartite #' @concept quantitative nodal attribute InitErgmTerm.distance <- function(nw, arglist, ...) { a <- check.ErgmTerm(nw, arglist, - varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff"), - vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric"), - required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), - defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0)) + varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale"), + vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric"), + required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), + defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1)) #Process arguments - if(is.character(a$coord)){ - a$coord<-get.network.attribute(nw,a$coord) + if(is.character(a$coord)){ #Some sort of attribute + if(length(a$coord)>1){ #Vertex attributes + a$coord<-ergm_get_vattr(a$coord,nw,accept="numeric",multiple="matrix") + }else{ #Either vertex or network attribute + if(a$coord%in%list.vertex.attributes(nw)) + a$coord<-ergm_get_vattr(a$coord,nw,accept="numeric") + else + a$coord<-get.network.attribute(nw,a$coord) + } if(is.null(a$coord)) - stop("Distance term requires either a coordinate matrix, or the name of a network attribute containing one.") + stop("Distance term requires either a coordinate matrix, the name of a network attribute containing one, or one or more names of vertex attributes containing coordinates.") } if(length(dim(a$coord))==2){ coord<-a$coord @@ -2619,9 +2660,11 @@ InitErgmTerm.distance <- function(nw, arglist, ...) { stop("Distance term requires that coordinates be provided for all vertices.") if(any(!apply(coord,1:2,is.numeric))) stop("Distance term requires numeric coordinates.") + if(any(is.na(coord))) + stop("Missing coordinate values not allowed in distance term.") if(a$sphere){ - if(NCOL(coord)!=2) - stop("Spherical coordinates chosen for distance term, but more or fewer than two dimensions provided; coordinates must be lat,lon or equivalent (in that order, and in angular units).") + if(NCOL(coord)<2) + stop("Spherical coordinates chosen for distance term, but fewer than two dimensions provided; first two coordinates must be lat,lon or equivalent (in that order, and in angular units).") if(any(coord[,1]< -90)||any(coord[,1]>90)) stop("Illegal latitude passed to distance term.") if(any(coord[,2]< -180)||any(coord[,2]>180)) @@ -2639,12 +2682,14 @@ InitErgmTerm.distance <- function(nw, arglist, ...) { else logstr<-"" if(a$sphere) - basestr<-"dist" + basestr<-"dist.S" else - basestr<-paste0("dist.L",a$metric) + basestr<-"dist" + if((!a$sphere)||(NCOL(coord)>2)) + basestr<-paste0(basestr,".L",a$metric) list(name = "distance", coef.names = paste0(basestr,logstr), - inputs = c(is.directed(nw), network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, coord), + inputs = c(network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, a$scale, coord), dependence = FALSE, emptynwstats = 0 ) diff --git a/man/distance-ergmTerm-9bc71012.Rd b/man/distance-ergmTerm-9bc71012.Rd index b3346b3ad..ed989b9ee 100644 --- a/man/distance-ergmTerm-9bc71012.Rd +++ b/man/distance-ergmTerm-9bc71012.Rd @@ -6,18 +6,21 @@ \title{Inter-point distances} \usage{ # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, -# log=TRUE, mindist=1e-5, distoff=0) +# log=TRUE, mindist=1e-5, distoff=0, scale=1) } \arguments{ -\item{coord}{node by dimension coordinate matrix, or name of a network -attribute containing the coordinate matrix; for spherical -coordinates, must contain two columns of angular coordinates -in lat/lon order} +\item{coord}{node by dimension coordinate matrix, name of a network +attribute containing the coordinate matrix, or vector of +vertex attribute names containing coordinates (to be used +jointly); for \code{sphere==TRUE}, first two coordinates +must be angular units, in lat/lon order} \item{metric}{power to use for the Minkowski metric} \item{sphere}{logical; should great circle distances be used (rather than -Minkowski distances in free space)?} +Minkowski distances in free space)? If \code{TRUE}, the first +two coordinates are respectively interpreted as latitude and +longitude coordinates on a sphere of specified radius} \item{radius}{for spherical distances, the radius of the sphere to use; defaults to the IUGG mean Earth radius, in km} @@ -31,14 +34,18 @@ raw distances} \item{distoff}{for log distances, an offset to be added to observed distances before taking the logarithm} + +\item{scale}{factor by which raw distances should be rescaled prior to other +operations (notably, offsetting, logging, and thresholding)} } \description{ This term adds a single statistic to the model whose value is the sum over all edge variables of the edge variable value times either the distance between the respective vertices, or its log (if \code{log==TRUE}). The \code{coord} agument must contain a vector, -matrix, or \code{data.frame} of coordinates (vertex x dimension), or -else the name of a network attribute with such coordinates. +matrix, or \code{data.frame} of coordinates (vertex x dimension), the +name of a network attribute with such coordinates, or a vector of vertex +attributes to be used as coordinates. } \details{ Either spherical (\code{sphere=TRUE}) or Minkowski metrics @@ -49,6 +56,8 @@ Minkowski metric is determined by \code{metric}, with the distance being given by \deqn{ D(i,j) = (sum_d |x_{id}-x_{jd}|^p)^(1/p) + }{ + D(i,j) = \left(\sum_d |x_{id}-x_{jd}|^p\right)^{1/p} } where the sum is over the dimensions of the space, \eqn{x} is the coordinate matrix, and \eqn{p} is the metric parameter (i.e. @@ -56,7 +65,7 @@ coordinate matrix, and \eqn{p} is the metric parameter (i.e. Euclidean distance. Note that \code{metric==1}, and \code{log==FALSE} in the one-dimensional case is identical to the \code{absdiff} term. -When \code{sphere=TRUE}, \code{coord} must be a two-dimensional matrix of +When \code{sphere=TRUE}, the first two dimensions of \code{coord} must be angular coordinates in lat/lon form (i.e., the first column must contain units of decimal degrees between -90 and 90, and the second must contain units of decimal degrees betweein -180 and 180). Distance is then @@ -66,6 +75,29 @@ radius in km, and hence supplying lat/lon coordinates yields geospherical distances in kilometers. As \code{metric} has no meaning here, it is ignored. +If \code{sphere=TRUE} and more than two coordinates are provided per +vertex, then the first two are used to compute great circle distances, +with the remainder being used to compute the selected Minkowski metric. +These two distances are then added to form the final distance. While +this has exotic uses, the most obvious is to add elevation. It is +important to ensure that the coordinates are specified such that the +two distances have the same units! (By default, spherical distances are +in kilometers, so unless \code{radius} is changed, the user should ensure +that additional coordinates are also in kilometers.) Note that elevation +effects specified in this way do not take into account factors such as the +need to go down to the street and then up to a target floor within another +building (as in the case of artificial elevation in cities), nor the need +to e.g. get up before you get down in hilly terrain. Life is full of such +disappointments. + +The overall scale of the computed distances can be modified by setting the +\code{scale} argument. This is most often useful for numerical or +interpretational reasons (e.g., to avoid overflow/underflow issues, or to +change to familiar units). However, it can also be used in conjunction +with offsets and log scaling to change the form of the distance effect, +as described below. In particular, note that rescaling is performed prior +to application of offsets or thresholding, where applicable. + When \code{log==TRUE}, the logarithm of the distance rather than the distance itself is used to form the statistic. This produces interaction functions that roughly approximate the power law spatial interaction @@ -75,10 +107,12 @@ used to provide an offset to the raw distance (with the computed value becoming \code{log(d+distoff)}. (This can have an effect that is similar to an \dQuote{attenuated power law} SIF, per Butts and Acton.) As an additional failsafe, the (offset) distances are thresholded from below by -\code{mindist}. Setting both \code{mindist} and \code{distoff} to zero is +\code{mindist}. With all elements considered, the net effect is thus +\code{log(max(mindist, distoff + scale*d))}, where \eqn{d} is the raw +distance. Setting both \code{mindist} and \code{distoff} to zero is allowed, but may produce exciting results if points precisely overlap. Note that both \code{distoff} and \code{mindist} are ignored when -\code{log==FALSE}. +\code{log==FALSE} (but \code{scale} is not). This term can be used for directed or undirected networks. } @@ -191,6 +225,7 @@ Research}, 222--250. SAGE Publications. \Sexpr[results=rd,stage=render]{ergm:::.formatTermKeywords("ergmTerm", "distance", "subsection")} } +\concept{bipartite} \concept{directed} \concept{dyad-independent} \concept{quantitative nodal attribute} diff --git a/src/changestats.c b/src/changestats.c index c3a1a9c73..b364d9ff5 100644 --- a/src/changestats.c +++ b/src/changestats.c @@ -1465,7 +1465,7 @@ C_CHANGESTAT_FN(c_degree_w_homophily) { /***************** - changestat: d_distance + changestat: c_distance *****************/ /* Utility function to calculate distances on a sphere of radius r, between two points @@ -1491,60 +1491,60 @@ double spheredist(double lat0, double lon0, double lat1, double lon1, double r) } /*Actual distance changestat function itself (calls spheredist)*/ -CHANGESTAT_FN(d_distance) { - Vertex t, h; +C_CHANGESTAT_FN(c_distance) { + Vertex t,h; int i,j; int isdir,nv,dim,logd,sphd; - double dis,dexp,*coord,sphr,logmind,logdoff; + double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis; /*Set things up*/ - isdir = (int)INPUT_PARAM[0]; /*Is the graph directed?*/ - nv = (int)INPUT_PARAM[1]; /*Number of vertices*/ - dim = (int)INPUT_PARAM[2]; /*Dimension of coordinate space*/ - dexp = INPUT_PARAM[3]; /*Minkowski exponent*/ - logd = (int)INPUT_PARAM[4]; /*Should we use the log distance?*/ - sphd = (int)INPUT_PARAM[5]; /*Should we use lat/lon distances on the geosphere?*/ - sphr = INPUT_PARAM[6]; /*Radius to use for spherical coordinates*/ - logmind = INPUT_PARAM[7]; /*Distance minimum for the log case*/ - logdoff = INPUT_PARAM[8]; /*Distance offset for the log case*/ - coord = INPUT_PARAM+9; /*Pointer to the coordinate matrix*/ - - /*Compute the changescores*/ - ZERO_ALL_CHANGESTATS(i); - FOR_EACH_TOGGLE(i) { - t = TAIL(i); h = HEAD(i); - if(sphd){ /*Compute distances on the sphere*/ - dis=spheredist(coord[t-1],coord[t-1+nv],coord[h-1],coord[h-1+nv],sphr); - }else{ /*Compute Minkowski distances in free space*/ - dis=0.0; - for(j=0;j2)*/ + mdis=0.0; + for(j=0+2*sphd;j0]), ignore_attr=TRUE) }) @@ -41,6 +45,12 @@ test_that("L2 works", { test_that("log L1 works", { expect_equal(summary(net~distance(co,metric=1)), sum(log(d1[g>0])), ignore_attr=TRUE) }) +test_that("scaling works", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE,scale=3)), 3*sum(d1[g>0]), ignore_attr=TRUE) +}) +test_that("thresholding works", { + expect_equal(summary(net~distance(co,metric=1,log=TRUE,scale=1e-10,mindist=5)), log(5)*sum(g>0), ignore_attr=TRUE) +}) #Verify that alternative specifications work test_that("character matches matrix", { @@ -49,9 +59,13 @@ test_that("character matches matrix", { test_that("character matches data.frame", { expect_equal(summary(net~distance(as.data.frame(co))), summary(net~distance("co")), ignore_attr=TRUE) }) +test_that("vertex attributes matches matrix", { + expect_equal(summary(net~distance(co)), summary(net~distance(c("x","y","z"))), ignore_attr=TRUE) +}) +test_that("single attribute defaults to vertex", { + expect_equal(summary(net~distance(co[,1])), summary(net~distance("x")), ignore_attr=TRUE) +}) -#Verify that values are correct -summary(net~distance(co)) #Test for lat/lon calculations co<-rbind( From 602c05072d4e200403a0d2ce84ac17c4c3bdf4af Mon Sep 17 00:00:00 2001 From: "Carter T. Butts" Date: Tue, 28 Jul 2026 02:42:55 -0700 Subject: [PATCH 03/18] The tests weren't happy about a concept label. --- R/InitErgmTerm.R | 1 - 1 file changed, 1 deletion(-) diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index 932d6119e..c560c6d88 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -2626,7 +2626,6 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' @concept dyad-independent #' @concept directed #' @concept undirected -#' @concept bipartite #' @concept quantitative nodal attribute InitErgmTerm.distance <- function(nw, arglist, ...) { a <- check.ErgmTerm(nw, arglist, From 934b8274daf50b1fa1d4be27059388937b9a80ab Mon Sep 17 00:00:00 2001 From: "Carter T. Butts" Date: Tue, 28 Jul 2026 02:50:45 -0700 Subject: [PATCH 04/18] Forgot to roxygenize the last change. This completes it. --- man/distance-ergmTerm-9bc71012.Rd | 1 - 1 file changed, 1 deletion(-) diff --git a/man/distance-ergmTerm-9bc71012.Rd b/man/distance-ergmTerm-9bc71012.Rd index ed989b9ee..5cd5ceefe 100644 --- a/man/distance-ergmTerm-9bc71012.Rd +++ b/man/distance-ergmTerm-9bc71012.Rd @@ -225,7 +225,6 @@ Research}, 222--250. SAGE Publications. \Sexpr[results=rd,stage=render]{ergm:::.formatTermKeywords("ergmTerm", "distance", "subsection")} } -\concept{bipartite} \concept{directed} \concept{dyad-independent} \concept{quantitative nodal attribute} From c1daceb7af4fb6bcd1cc779b126b65f8594f4baa Mon Sep 17 00:00:00 2001 From: "Carter T. Butts" Date: Fri, 31 Jul 2026 14:59:05 -0700 Subject: [PATCH 05/18] Added the pow option to the distance term, and documented it. Added a test, as well. Seems to work fine, and passes tests. --- R/InitErgmTerm.R | 22 +++++++++++++++------- man/distance-ergmTerm-9bc71012.Rd | 11 +++++++++-- src/changestats.c | 15 +++++++++++---- tests/testthat/test-distance.R | 3 +++ 4 files changed, 38 insertions(+), 13 deletions(-) diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index c560c6d88..a02c43af5 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -2427,7 +2427,7 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' #' @usage #' # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, -#' # log=TRUE, mindist=1e-5, distoff=0, scale=1) +#' # log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1) #' #' @param coord node by dimension coordinate matrix, name of a network #' attribute containing the coordinate matrix, or vector of @@ -2450,6 +2450,8 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' distances before taking the logarithm #' @param scale factor by which raw distances should be rescaled prior to other #' operations (notably, offsetting, logging, and thresholding) +#' @param pow power to which scaled distances should be raised prior to other +#' other operations (notably, offsetting, logging, and thresholding) #' #' @details Either spherical (\code{sphere=TRUE}) or Minkowski metrics #' (\code{sphere=FALSE}) may be selected. For the latter, any number of @@ -2499,7 +2501,11 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' change to familiar units). However, it can also be used in conjunction #' with offsets and log scaling to change the form of the distance effect, #' as described below. In particular, note that rescaling is performed prior -#' to application of offsets or thresholding, where applicable. +#' to application of offsets or thresholding, where applicable. Similarly, +#' the \code{pow} argument can be used to raise the scaled distance to an +#' arbitrary power prior to further calculation. (Note, however, that power +#' law spatial dependence is not realized via this mechanism, but by working +#' with log distances, as discussed below.) #' #' When \code{log==TRUE}, the logarithm of the distance rather than the #' distance itself is used to form the statistic. This produces interaction @@ -2629,10 +2635,10 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' @concept quantitative nodal attribute InitErgmTerm.distance <- function(nw, arglist, ...) { a <- check.ErgmTerm(nw, arglist, - varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale"), - vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric"), - required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), - defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1)) + varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale", "pow"), + vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric", "numeric"), + required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), + defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1,1)) #Process arguments if(is.character(a$coord)){ #Some sort of attribute if(length(a$coord)>1){ #Vertex attributes @@ -2686,9 +2692,11 @@ InitErgmTerm.distance <- function(nw, arglist, ...) { basestr<-"dist" if((!a$sphere)||(NCOL(coord)>2)) basestr<-paste0(basestr,".L",a$metric) + if(a$pow!=1) + basestr<-paste0(basestr,".pow",a$pow) list(name = "distance", coef.names = paste0(basestr,logstr), - inputs = c(network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, a$scale, coord), + inputs = c(network.size(nw), NCOL(coord), a$metric, a$log, a$sphere, a$radius, a$mindist, a$distoff, a$scale, a$pow, coord), dependence = FALSE, emptynwstats = 0 ) diff --git a/man/distance-ergmTerm-9bc71012.Rd b/man/distance-ergmTerm-9bc71012.Rd index 5cd5ceefe..f48420796 100644 --- a/man/distance-ergmTerm-9bc71012.Rd +++ b/man/distance-ergmTerm-9bc71012.Rd @@ -6,7 +6,7 @@ \title{Inter-point distances} \usage{ # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, -# log=TRUE, mindist=1e-5, distoff=0, scale=1) +# log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1) } \arguments{ \item{coord}{node by dimension coordinate matrix, name of a network @@ -37,6 +37,9 @@ distances before taking the logarithm} \item{scale}{factor by which raw distances should be rescaled prior to other operations (notably, offsetting, logging, and thresholding)} + +\item{pow}{power to which scaled distances should be raised prior to other +other operations (notably, offsetting, logging, and thresholding)} } \description{ This term adds a single statistic to the model whose value is @@ -96,7 +99,11 @@ interpretational reasons (e.g., to avoid overflow/underflow issues, or to change to familiar units). However, it can also be used in conjunction with offsets and log scaling to change the form of the distance effect, as described below. In particular, note that rescaling is performed prior -to application of offsets or thresholding, where applicable. +to application of offsets or thresholding, where applicable. Similarly, +the \code{pow} argument can be used to raise the scaled distance to an +arbitrary power prior to further calculation. (Note, however, that power +law spatial dependence is not realized via this mechanism, but by working +with log distances, as discussed below.) When \code{log==TRUE}, the logarithm of the distance rather than the distance itself is used to form the statistic. This produces interaction diff --git a/src/changestats.c b/src/changestats.c index b364d9ff5..9b9bc8d6e 100644 --- a/src/changestats.c +++ b/src/changestats.c @@ -1493,9 +1493,9 @@ double spheredist(double lat0, double lon0, double lat1, double lon1, double r) /*Actual distance changestat function itself (calls spheredist)*/ C_CHANGESTAT_FN(c_distance) { Vertex t,h; - int i,j; - int isdir,nv,dim,logd,sphd; - double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis; + int j; + int nv,dim,logd,sphd; + double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis,pw; /*Set things up*/ nv = (int)INPUT_PARAM[0]; /*Number of vertices*/ @@ -1507,7 +1507,8 @@ C_CHANGESTAT_FN(c_distance) { logmind = INPUT_PARAM[6]; /*Distance minimum for the log case*/ logdoff = INPUT_PARAM[7]; /*Distance offset for the log case*/ scale = INPUT_PARAM[8]; /*Scaling adjustment for raw distance*/ - coord = INPUT_PARAM+9; /*Pointer to the coordinate matrix*/ + pw = INPUT_PARAM[9]; /*Power to which distance should be raised*/ + coord = INPUT_PARAM+10; /*Pointer to the coordinate matrix*/ /*Compute the changescore*/ t=tail-1; /*0-indexed values, for C convenience*/ @@ -1536,6 +1537,12 @@ C_CHANGESTAT_FN(c_distance) { } dis+=mdis; /*Add the two distance types*/ dis*=scale; /*Apply the scaling coefficient*/ + if(pw==0.5) /*Apply the power transformation*/ + dis=sqrt(dis); + else if(pw==2.0) + dis*=dis; + else if(pw!=1.0) + dis=pow(dis,pw); /*If using log distances, log 'em. But to prevent divergence, threshold from below by logmind (after adding the offset).*/ if(logd) diff --git a/tests/testthat/test-distance.R b/tests/testthat/test-distance.R index b3df3bcba..912214cd5 100644 --- a/tests/testthat/test-distance.R +++ b/tests/testthat/test-distance.R @@ -48,6 +48,9 @@ test_that("log L1 works", { test_that("scaling works", { expect_equal(summary(net~distance(co,metric=1,log=FALSE,scale=3)), 3*sum(d1[g>0]), ignore_attr=TRUE) }) +test_that("powers work", { + expect_equal(summary(net~distance(co,metric=1,log=FALSE,pow=3)), sum(d1[g>0]^3), ignore_attr=TRUE) +}) test_that("thresholding works", { expect_equal(summary(net~distance(co,metric=1,log=TRUE,scale=1e-10,mindist=5)), log(5)*sum(g>0), ignore_attr=TRUE) }) From b9e6dcda31b0e2c1c6628c6188b9dd83642bf893 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:26:54 -0400 Subject: [PATCH 06/18] --Duplicate changestats.c history into changestat_spheredist.c --- src/{changestats.c => changestat_spheredist.c} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestats.c => changestat_spheredist.c} (100%) diff --git a/src/changestats.c b/src/changestat_spheredist.c similarity index 100% rename from src/changestats.c rename to src/changestat_spheredist.c From 0dacc82479ce04eb385b2927d8e01112faa86691 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:26:54 -0400 Subject: [PATCH 07/18] --Restore changestats.c --- src/changestats.c | 3515 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 3515 insertions(+) create mode 100644 src/changestats.c diff --git a/src/changestats.c b/src/changestats.c new file mode 100644 index 000000000..9b9bc8d6e --- /dev/null +++ b/src/changestats.c @@ -0,0 +1,3515 @@ +/* File src/changestats.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-2026 Statnet Commons + */ +#include "changestats.h" +#include "ergm_storage.h" +#include "ergm_dyad_hashmap.h" +#include "ergm_edgelist.h" + +/******************** changestats: A ***********/ + +/***************** + changestat: d_adegcor +*****************/ +D_CHANGESTAT_FN(d_adegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i),HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ +// CHANGE_STAT[0] = mtp->dstats[0] - current; +// Rprintf("c %f p %f",current,mtp->dstats[0]); + mtp->dstats[0] -= current; +// Rprintf(" p-c %f\n",mtp->dstats[0]); +FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_adegcor) { + Vertex tail, head, taildeg, headdeg; + Edge e; + double mu, mu2, sigma2, cross; + + mu = 0.0; + mu2 = 0.0; + cross = 0.0; + for(tail=1; tail <= N_NODES; tail++) { + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + mu += (double)(taildeg + headdeg); + mu2 += (double)(taildeg*taildeg + headdeg*headdeg); + cross += 2.0*taildeg*headdeg; + } + } + mu = mu / (2.0*N_EDGES); + sigma2 = mu2/(2.0*N_EDGES) - mu*mu; + CHANGE_STAT[0] = (cross / (2.0*N_EDGES) - mu*mu) / sigma2; +} + +/***************** + changestat: d_altkstar +*****************/ +C_CHANGESTAT_FN(c_altkstar) { + double lambda, oneexpl, change; + Vertex taild, headd=0; + + change = 0.0; + lambda = INPUT_PARAM[0]; + oneexpl = 1.0-1.0/lambda; + + /* *** don't forget tail -> head */ + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + if(taild!=0){ + change += (edgestate?-1:+1)*(1.0-pow(oneexpl,(double)taild)); + } + if(headd!=0){ + change += (edgestate?-1:+1)*(1.0-pow(oneexpl,(double)headd)); + } + CHANGE_STAT[0] = change*lambda; +} + +/***************** + changestat: d_asymmetric +*****************/ +C_CHANGESTAT_FN(c_asymmetric) { + double matchval, change; + int j, ninputs, noattr; + + ninputs = N_INPUT_PARAMS - N_NODES; + noattr = (N_INPUT_PARAMS == 0); + + /* *** don't forget tail -> head */ + change = (edgestate==IS_OUTEDGE(head, tail) ? 1.0 : -1.0) ; + if (noattr) { /* "plain vanilla" asymmetric, without node attributes */ + CHANGE_STAT[0] += change; + } else { /* Only consider asymmetrics where node attributes match */ + matchval = INPUT_PARAM[tail+ninputs-1]; + if (matchval == INPUT_PARAM[head+ninputs-1]) { /* We have a match! */ + if (ninputs==0) {/* diff=F in network statistic specification */ + CHANGE_STAT[0] += change; + } else { /* diff=T */ + for (j=0; j head */ + b1 = tail; + echange = IS_OUTEDGE(b1,head) ? -1 : 1; + b1deg = OUT_DEG[b1]; + CHANGE_STAT[0] += (b1deg + echange > 1) - (b1deg > 1); +} + +/***************** + changestat: d_b1concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_b1concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. */ + int j, echange, b1attr; + Vertex b1, b1deg; + + /* *** don't forget tail -> head */ + b1 = tail; + echange = IS_OUTEDGE(b1,head) ? -1 : 1; + b1deg = OUT_DEG[b1]; + b1attr = INPUT_PARAM[N_CHANGE_STATS + b1 - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (b1attr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (b1deg + echange > 1) - (b1deg > 1); + } + } +} + +/***************** + changestat: d_b1nodematch +*****************/ +C_CHANGESTAT_FN(c_b1nodematch) { + + Vertex node3, node4, ninputs; + int count, exponenttype, matchval, b2attrsize, attrval1, attrval2, diffstatus; + /* int j, numofstats; */ + Edge e, e2; + double beta, alpha, change=0.0, exponent; + const int BetaType=1, AlphaType=2; + + b2attrsize = INPUT_PARAM[0]; + + if(b2attrsize > 0){ + ninputs = N_INPUT_PARAMS - N_NODES - b2attrsize;/*have 2 sets of node attributes and b2attrvals */ + } + else{ + ninputs = N_INPUT_PARAMS - BIPARTITE; + } + + diffstatus = !(ninputs == 3); /* 1 if Diff = T and 0 if Diff = F */ + /* numofstats = diffstatus ? (b2attrsize == 0 ? (ninputs - 3): (ninputs - 3) * b2attrsize) : (b2attrsize == 0 ? 1 : b2attrsize); */ + + exponent = beta = INPUT_PARAM[1]; /* exponent on nodematch count */ + exponenttype = BetaType; + alpha = INPUT_PARAM[2]; + + if (beta >= 1.0 && alpha < 1.0) { + exponent = alpha; + exponenttype = AlphaType; + } + // Rprintf("N_INPUT_PARAMS = %d, N_NODES=%d\n", N_INPUT_PARAMS, N_NODES); + // Rprintf("ninputs = %d, beta=%f, alpha=%f, exponenttype=%d, exponent=%f\n", + // ninputs, beta, alpha, exponenttype, exponent); + + matchval = INPUT_PARAM[tail + ninputs - 1]; + + /* Now count the neighbors of head whose attribute value equals matchval */ + /* All neighbors of head are inedges because this is a bipartite network */ + count = 0; + change = 0.0; + + if(b2attrsize == 0){ + + STEP_THROUGH_INEDGES(head, e, node3) { + if (INPUT_PARAM[node3 + ninputs - 1] == matchval && tail != node3) { /* match! */ + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + count = 0; + + STEP_THROUGH_OUTEDGES(tail, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != head) { /* RPB */ + count += IS_OUTEDGE(node3, node4); /* add 1 if node4 connects node3 with tail */ + } + } + + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + change += (count==0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + + if (exponenttype == BetaType && count>0) { + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, exponent); + change -= 0.5*count*(exponent==0.0? (count==1? 0.0 : 1.0) : pow((count-1), exponent)); + } + + if (diffstatus) { /* diff=T */ + if(ninputs==4) /* keep != NULL*/ + CHANGE_STAT[0] += edgestate ? -change : change; + else + CHANGE_STAT[matchval-1] += edgestate ? -change : change; + + } else { /* diff=F */ + CHANGE_STAT[0] += edgestate ? -change : change; + } + + } else { + + attrval1 = INPUT_PARAM[head + ninputs + b2attrsize - 1]; + + STEP_THROUGH_INEDGES(head, e, node3) { + + if (INPUT_PARAM[node3 + ninputs - 1] == matchval && tail != node3) { /* match! */ + + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + + count = 0; + + STEP_THROUGH_OUTEDGES(tail, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != head) { + attrval2 = INPUT_PARAM[node4 + ninputs + b2attrsize - 1]; + if(attrval2 == attrval1) count += IS_OUTEDGE(node3, node4); + } + } + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + /* setting the change stat for each parameter */ + change += (count== 0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + if (exponenttype == BetaType && count > 0) { + + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + if(diffstatus){ + CHANGE_STAT[b2attrsize*(matchval-1) + attrval1 - 1] += edgestate ? -change : change; + } else{ + CHANGE_STAT[attrval1 - 1] += edgestate ? -change : change; + } + + } +} + +/***************** + changestat: d_b1starmix +*****************/ +C_CHANGESTAT_FN(c_b1starmix) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, taild; + int nstats; + double tailattr, headattr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + taild = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + if(headattr == INPUT_ATTRIB[node3-1]){++taild;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr && + INPUT_ATTRIB[nnodes+nstats+j] == headattr) { + change = CHOOSE(taild, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b1starmixhomophily +*****************/ +C_CHANGESTAT_FN(c_b1starmixhomophily) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, taild; + double tailattr, headattr; + + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + taild = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + if(headattr == INPUT_ATTRIB[node3-1]){++taild;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr) { + change = CHOOSE(taild, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b1twostar +*****************/ +C_CHANGESTAT_FN(c_b1twostar) { + double change; + int j; + Edge e; + Vertex node3, nnodes; + int nstats; + double tailattr, headattr, n3attr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + + /* *** don't forget tail -> head */ + change = IS_OUTEDGE(tail = tail, head)? -1.0 : 1.0 ; + tailattr = INPUT_PARAM[tail-1]; + headattr = INPUT_PARAM[head-1]; + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + n3attr = INPUT_PARAM[node3-1]; + for(j=0; j < N_CHANGE_STATS; j++) { + if (node3 != head && INPUT_PARAM[nnodes + j] == tailattr && + INPUT_PARAM[nnodes + nstats + j] == MIN(headattr, n3attr) && + INPUT_PARAM[nnodes + 2*nstats + j] == MAX(headattr, n3attr)) { + CHANGE_STAT[j] += change; + } + } + } +} + +/***************** + changestat: d_b2concurrent +*****************/ +C_CHANGESTAT_FN(c_b2concurrent) { + int echange; + Vertex b2, b2deg; + + /* *** don't forget tail -> head */ + b2 = head; + echange = IS_OUTEDGE(tail, b2) ? -1 : 1; + b2deg = IN_DEG[b2]; + CHANGE_STAT[0] += (b2deg + echange > 1) - (b2deg > 1); +} + +/***************** + changestat: d_b2concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_b2concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes.*/ + int j, echange, b2attr; + Vertex b2, b2deg; + + + /* *** don't forget tail -> head */ + b2 = head; + echange = IS_OUTEDGE(tail, b2) ? -1 : 1; + b2deg = IN_DEG[b2]; + b2attr = INPUT_PARAM[N_CHANGE_STATS + b2 - 1 - BIPARTITE]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (b2attr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (b2deg + echange > 1) - (b2deg > 1); + } + } +} + +/***************** + changestat: d_b2nodematch +*****************/ +C_CHANGESTAT_FN(c_b2nodematch) { + + Vertex node3, node4, ninputs; + int count, exponenttype, matchval, b1attrsize, attrval1, attrval2, diffstatus; + /* int j, ind, numofstats; */ + Edge e, e2; + double beta, alpha, change=0.0, exponent; + const int BetaType=1, AlphaType=2; + + b1attrsize = INPUT_PARAM[0]; + + if(b1attrsize > 0){ + ninputs = N_INPUT_PARAMS - N_NODES - b1attrsize;/*have 2 sets of node attributes and b2attrvals */ + } + else{ + ninputs = N_INPUT_PARAMS - N_NODES + BIPARTITE; + } + + diffstatus = !(ninputs == 3); /* 1 if Diff = T and o if Diff = F - RPB */ + /* numofstats = diffstatus ? (b1attrsize == 0 ? (ninputs - 3): (ninputs - 3) * b1attrsize) : (b1attrsize == 0 ? 1 : b1attrsize); */ + + exponent = beta = INPUT_PARAM[1]; /* exponent on nodematch count */ + exponenttype = BetaType; + alpha = INPUT_PARAM[2]; + if (beta >= 1.0 && alpha < 1.0) { + exponent = alpha; + exponenttype = AlphaType; + } + // Rprintf("N_INPUT_PARAMS = %d, N_NODES=%d\n", N_INPUT_PARAMS, N_NODES); + // Rprintf("ninputs = %d, beta=%f, alpha=%f, exponenttype=%d, exponent=%f\n", + // ninputs, beta, alpha, exponenttype, exponent); + + matchval = INPUT_PARAM[head + ninputs - BIPARTITE - 1]; + /* Now count the neighbors of tail whose attribute value equals matchval */ + /* All neighbors of tail are outedges because this is a bipartite network */ + count=0; + change = 0.0; + + /* RPB */ + /* double CHANGE[b1attrsize]; */ + + + if(b1attrsize == 0){ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { + if (INPUT_PARAM[node3 + ninputs - BIPARTITE - 1] == matchval && head != node3) { /* match! */ + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting head and node3 */ + count = 0; + + STEP_THROUGH_INEDGES(head, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != tail) { + count += IS_OUTEDGE(node4, node3); /* add 1 if node4 connects node3 with head */ + } + } + /* if count==0, then the statistic is always 1 */ + // Rprintf("count is %d\n", count); + change += (count==0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + + if (exponenttype == BetaType && count>0) { + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + + if (diffstatus) { /* diff=T */ + if(ninputs==4) /* keep != NULL*/ + CHANGE_STAT[0] += edgestate ? -change : change; + else + CHANGE_STAT[matchval-1] += edgestate ? -change : change; + + } else { /* diff=F */ + CHANGE_STAT[0] += edgestate ? -change : change; + } + } else { + + attrval1 = INPUT_PARAM[tail + ninputs + N_NODES + b1attrsize - BIPARTITE - 1]; + + STEP_THROUGH_OUTEDGES(tail, e, node3) { + + if (INPUT_PARAM[node3 + ninputs - BIPARTITE - 1] == matchval && head != node3) { /* match! */ + + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + + count = 0; + + STEP_THROUGH_INEDGES(head, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != tail) { + attrval2 = INPUT_PARAM[node4 + ninputs + N_NODES + b1attrsize - BIPARTITE - 1]; + if(attrval2 == attrval1) count += IS_OUTEDGE(node4, node3); + } + } + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + /* setting the change stat for each parameter */ + change += (count== 0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + if (exponenttype == BetaType && count > 0) { + + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + if(diffstatus){ + CHANGE_STAT[b1attrsize*(matchval-1) + attrval1 - 1] += edgestate ? -change : change; + } else{ + CHANGE_STAT[attrval1 - 1] += edgestate ? -change : change; + } + + } +} + +/***************** + changestat: d_b2starmix +*****************/ +C_CHANGESTAT_FN(c_b2starmix) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, headd; + int nstats; + double tailattr, headattr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + headd = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr && + INPUT_ATTRIB[nnodes+nstats+j] == headattr) { + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b2starmixhomophily +*****************/ +C_CHANGESTAT_FN(c_b2starmixhomophily) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, headd; + double tailattr, headattr; + + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + headd = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == headattr) { + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b2twostar +*****************/ +C_CHANGESTAT_FN(c_b2twostar) { + double change; + int j; + Edge e; + Vertex node3, nnodes; + int nstats; + double tailattr, headattr, n3attr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + + + /* *** don't forget tail -> head */ + change = edgestate? -1.0 : 1.0 ; + tailattr = INPUT_PARAM[tail-1]; + headattr = INPUT_PARAM[head-1]; + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + n3attr = INPUT_PARAM[node3-1]; + for(j=0; j < N_CHANGE_STATS; j++) { + if (node3 != tail && INPUT_PARAM[nnodes + j] == headattr && + INPUT_PARAM[nnodes + nstats + j] == MIN(tailattr, n3attr) && + INPUT_PARAM[nnodes + 2*nstats + j] == MAX(tailattr, n3attr)) { + CHANGE_STAT[j] += change; + } + } + } +} + +/***************** + changestat: d_balance +*****************/ +C_CHANGESTAT_FN(c_balance) { + int a, b, c, d, e, edgecount, t300, + t210, t120C, t120U, t120D, t201, t030C, t030T, t111U, + t111D, t021C, t021U, t021D, t102, t012; /* , t003; */ + Vertex node3; + + + /* *** don't forget tail -> head */ + if (DIRECTED) { /* directed version */ + t300 = 0; + t210 = 0; + t120C = 0; t120U = 0; t120D = 0; t201 = 0; + t030C = 0; t030T = 0; t111U = 0; t111D = 0; + t021C = 0; t021U = 0; t021D = 0; t102 = 0; + t012 = 0; + + if (MIN_OUTEDGE(head)!=0 || MIN_INEDGE(head)!=0 || + MIN_OUTEDGE(tail)!=0 || MIN_INEDGE(tail)!=0) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = IS_OUTEDGE(head, tail); + b = IS_OUTEDGE(head, node3); + c = IS_OUTEDGE(node3, head); + d = IS_OUTEDGE(node3, tail); + e = IS_OUTEDGE(tail, node3); + edgecount = (a + b + c + d + e); + + switch(edgecount) { + case 0: /* 012 */ + ++t012; + + case 1: { /* 021C, 021U, 021D, 102 */ + if ((b == 1) || (d == 1)) + ++t021C; + if (c == 1) + ++t021U; + if (e == 1) + ++t021D; + if (a == 1) + ++t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if ((b + d) == 2) + ++t030C; + if (((b + e) == 2) || ((c + d) == 2) || ((c + e) == 2)) + ++t030T; + if (((a + b) == 2) || ((a + e) == 2) || ((d + e) == 2)) + ++t111U; + if (((a + c) == 2) || ((a + d) == 2) || ((b + c) == 2)) + ++t111D; + } + break; + + case 3: { /* 120C, 120U, 120D, 201 */ + if (a == 1) { + if (((b + d) == 2) || ((c + e) == 2)) + ++t120C; + if ((b + e) == 2) + ++t120U; + if ((c + d) == 2) + ++t120D; + if (((b + c) == 2) || ((d + e) == 2)) + ++t201; + } else { + if (b == 1) { + if (((c + d) == 2) || ((d + e) == 2)) + ++t120C; + if ((c + e) == 2) + ++t120D; + } else { + ++t120U; + } + } + } + break; + + case 4: /* 210 */ + ++t210; + break; + + case 5: /* 300 */ + ++t300; + break; + } + + switch(edgecount) { + case 1: /* 102, 021D, 021U, 021C */ + --t012; + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if (((a + c) == 2) || ((a + e) == 2) || ((b + d) == 2) || + ((c + e) == 2)) + --t021C; + if (((a + d) == 2) || ((b + e) == 2)) + --t021U; + if (((a + b) == 2) || ((c + d) == 2)) + --t021D; + if (((b + c) == 2) || ((d + e) == 2)) + --t102; + } + break; + + case 3: { /* 201, 120D, 120U, 120C */ + if (a == 1) { + if ((c + e) == 2) + --t030C; + if (((c + d) == 2) || ((b + e) == 2) || ((b + d) == 2)) + --t030T; + if ((b + c) == 2) + --t111U; + if ((d + e) == 2) + --t111D; + } else { + if (b == 1) { + if ((c + d) == 2) + --t111U; + if (((c + e) == 2) || ((d + e) == 2)) + --t111D; + } + else + --t111U; + } + } + break; + + case 4: { /* 210 */ + if (a == 1) { + if (((b + c + e) == 3) || ((c + d + e) == 3)) + --t120C; + if ((b + c + d) == 3) + --t120U; + if ((b + d + e) == 3) + --t120D; + } else { + if ((b + c + d + e) == 4) + --t201; + } + } + break; + + case 5: /* 300 */ + --t210; + break; + } + } + } /* ****** move to next node3 ******** */ + } + else + t012 = t012 + (N_NODES - 2); + + /* t003 = (t300+t210+t120C+t120U+t120D+t201+t030C+t030T); + t003 = t003+(t111U+t111D+t021C+t021U+t021D+t102+t012); */ + b = t102 + t300; + CHANGE_STAT[0] += edgestate ? -(double)b : (double)b; + + /* *** don't forget tail -> head */ + }else{ /* undirected */ + t300 = 0; t201 = 0; t102 = 0; t012 = 0; + + if (MIN_OUTEDGE(head)!=0 || MIN_INEDGE(head)!=0 || + MIN_OUTEDGE(tail)!=0 || MIN_INEDGE(tail)!=0) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = IS_UNDIRECTED_EDGE(node3, head); + b = IS_UNDIRECTED_EDGE(node3, tail); + edgecount = (a + b); + + switch(edgecount){ + case 0: { /* 012 */ + ++t102; + --t012; + } + break; + + case 1: { /* 021C, 021U, 021D, 102 */ + ++t201; + --t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + ++t300; + --t201; + } + break; + } + } + + } /* ****** move to next node3 ******** */ + } else + t102 = t102 + (N_NODES - 2); + + /* t003 = (t102+t201+t300); */ + b = t102 + t300; + CHANGE_STAT[0] += edgestate ? -(double)b : (double)b; + } +} + +/***************** + changestat: d_boundeddegree +*****************/ +C_CHANGESTAT_FN(c_boundeddegree) { + int j, echange; + Vertex taild, headd=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = DEG(tail); + headd = DEG(head); + for(j = 0; j+1 < nstats; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + CHANGE_STAT[j] += (headd + echange == deg) - (headd == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); + CHANGE_STAT[nstats-1] += (headd + echange >= bound) - (headd >= bound); +} + +/***************** + changestat: d_boundedidegree +*****************/ +C_CHANGESTAT_FN(c_boundedidegree) { + int j, echange; + Vertex taild=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = IN_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); +} + +/***************** + changestat: d_boundedistar +*****************/ +C_CHANGESTAT_FN(c_boundedistar) { + double change, headod; + double newheadod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + headod = IN_DEG[head]; + newheadod = headod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = MIN(bound,CHOOSE(newheadod, k))-MIN(bound,CHOOSE(headod, k)); + CHANGE_STAT[j] += change; + } +} + +/***************** + changestat: d_boundedkstar +*****************/ +C_CHANGESTAT_FN(c_boundedkstar) { + double change, tailod, headod; + double newtailod, newheadod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + tailod = DEG(tail); + newtailod = tailod + (edgestate ? -1 : 1); + headod = DEG(head); + newheadod = headod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = (MIN(bound,CHOOSE(newtailod, k))-MIN(bound,CHOOSE(tailod, k))) + + (MIN(bound,CHOOSE(newheadod, k))-MIN(bound,CHOOSE(headod, k))); + + CHANGE_STAT[j] += change; /* (edgestate ? - change : change); */ + } +} + +/***************** + changestat: d_boundedodegree +*****************/ +C_CHANGESTAT_FN(c_boundedodegree) { + int j, echange; + Vertex taild=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = OUT_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); +} + +/***************** + changestat: d_boundedostar +*****************/ +C_CHANGESTAT_FN(c_boundedostar) { + double change, tailod; + double newtailod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + tailod = OUT_DEG[tail]; + newtailod = tailod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = MIN(bound,CHOOSE(newtailod, k))-MIN(bound,CHOOSE(tailod, k)); + CHANGE_STAT[j] += change; + } + } + +/***************** + changestat: d_boundedtriangle +*****************/ +Vertex CountTriangles (Vertex tail, Vertex head, int outcount, + int incount, Network *nwp); +C_CHANGESTAT_FN(c_boundedtriangle) { + Edge e; + Vertex node3; + double boundedchange, htcount; + Vertex tailtri, headtri; + int bound = (int)INPUT_PARAM[0]; + + /* *** don't forget tail -> head */ + tailtri=0; + headtri=0; + STEP_THROUGH_OUTEDGES(tail, e, node3) { + tailtri += CountTriangles(tail, node3, 1, 1, nwp); + } + STEP_THROUGH_INEDGES(tail, e, node3) { + tailtri += CountTriangles(tail, node3, 1, 1, nwp); + } + STEP_THROUGH_OUTEDGES(head, e, node3) { + headtri += CountTriangles(head, node3, 1, 1, nwp); + } + STEP_THROUGH_INEDGES(head, e, node3) { + headtri += CountTriangles(head, node3, 1, 1, nwp); + } + tailtri = tailtri/2; + headtri = headtri/2; + htcount = CountTriangles(tail, head, 1, 1, nwp); + boundedchange = (MIN(headtri+(edgestate ? -1:1)*htcount,bound)-MIN(headtri,bound)+ + MIN(tailtri+(edgestate ? -1:1)*htcount,bound)-MIN(tailtri,bound)); + CHANGE_STAT[0] += boundedchange; +} + +/***************** + CountTriangles: called by d_boundedtriangle +*****************/ +Vertex CountTriangles (Vertex tail, Vertex head, int outcount, int incount, + Network *nwp) { + Edge e; + Vertex change; + Vertex k; + + /* *** don't forget tail -> head */ + change=0; + if(outcount){ + STEP_THROUGH_OUTEDGES(head, e, k) /* step through outedges of head */ + { + if (IS_UNDIRECTED_EDGE(k,tail)) + ++change; + } + } + + if(incount){ + STEP_THROUGH_INEDGES(head, e, k) /* step through inedges of head */ + { + if (IS_UNDIRECTED_EDGE(k,tail)) + ++change; + } + } + return(change); +} + + + +/******************** changestats: C ***********/ +/***************** + changestat: d_concurrent +*****************/ +C_CHANGESTAT_FN(c_concurrent) { + int echange; + Vertex taildeg, headdeg; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + if(!DIRECTED){ + taildeg += IN_DEG[tail]; + headdeg += OUT_DEG[head]; + } + CHANGE_STAT[0] += (taildeg + echange > 1) - (taildeg > 1); + CHANGE_STAT[0] += (headdeg + echange > 1) - (headdeg > 1); +} + +/***************** + changestat: d_concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex taildeg, headdeg; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + if(!DIRECTED){ + taildeg += IN_DEG[tail]; + headdeg += OUT_DEG[head]; + } + tailattr = INPUT_PARAM[N_CHANGE_STATS + tail - 1]; + headattr = INPUT_PARAM[N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (tailattr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (taildeg + echange > 1) - (taildeg > 1); + } + if (headattr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (headdeg + echange > 1) - (headdeg > 1); + } + } +} + +/***************** + changestat: d_ctriple +*****************/ +C_CHANGESTAT_FN(c_ctriple) { + Edge e; + Vertex change, node3; + int j; + double tailattr, edgemult; + + /* *** don't forget tail -> head */ + edgemult = edgestate ? -1.0 : 1.0; + change = 0; + if(N_INPUT_PARAMS > 0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]) { + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_OUTEDGE(node3, tail); + } + if(N_CHANGE_STATS > 1) { /* diff = TRUE; matches must be tabled */ + for (j=0; j=(from) && (x)<(to)) + +/***************** + changestat: d_degrange +*****************/ +C_CHANGESTAT_FN(c_degrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex taildeg = DEG(tail), headdeg = DEG(head); + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } +} + +/***************** + changestat: d_degrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_degrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 3*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex taildeg = DEG(tail), headdeg = DEG(head); + int tailattr = INPUT_PARAM[3*N_CHANGE_STATS + tail - 1], + headattr = INPUT_PARAM[3*N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } +} + +/***************** + changestat: d_degrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_degrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of degrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + Vertex taildeg, headdeg, v; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + taildeg=headdeg=-1; /* since tailattr==headattr, subtract the automatic match */ + taildeg=headdeg=0; + STEP_THROUGH_OUTEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_INEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_OUTEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + STEP_THROUGH_INEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_degree +*****************/ +C_CHANGESTAT_FN(c_degree) { + int j, echange; + Vertex taildeg, headdeg, deg; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + taildeg = DEG(tail); + headdeg = DEG(head); + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } +} + +/***************** + changestat: c_degdist +*****************/ +C_CHANGESTAT_FN(c_degdist) { + int echange = edgestate ? -1:+1; + + Vertex otd = DEG(tail), ohd = DEG(head), + ntd = otd + echange, nhd = ohd + echange; + + if(ntd > N_CHANGE_STATS || nhd > N_CHANGE_STATS) cutoff_error(mtp); + + if(otd) CHANGE_STAT[otd-1]--; + if(ohd) CHANGE_STAT[ohd-1]--; + if(ntd) CHANGE_STAT[ntd-1]++; + if(nhd) CHANGE_STAT[nhd-1]++; +} + +/***************** + changestat: d_degreepopularity +*****************/ +C_CHANGESTAT_FN(c_degreepopularity) { + double change; + + /* *** don't forget tail -> head */ + change = 0.0; + Vertex tdeg = DEG(tail); + Vertex hdeg = DEG(head); + if(edgestate){ + change -= sqrt(tdeg); + change += (tdeg-1.0)*(sqrt(tdeg-1.0)-sqrt(tdeg)); + change -= sqrt(hdeg); + change += (hdeg-1.0)*(sqrt(hdeg-1.0)-sqrt(hdeg)); + }else{ + change += sqrt(tdeg+1.0); + change += tdeg*(sqrt(tdeg+1.0)-sqrt(tdeg)); + change += sqrt(hdeg+1.0); + change += hdeg*(sqrt(hdeg+1.0)-sqrt(hdeg)); + } + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_degree_by_attr +*****************/ +C_CHANGESTAT_FN(c_degree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr, testattr; + Vertex taildeg, headdeg, d; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:1; + taildeg = DEG(tail); + headdeg = DEG(head); + tailattr = INPUT_PARAM[2*N_CHANGE_STATS + tail - 1]; + headattr = INPUT_PARAM[2*N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += (taildeg + echange == d) - (taildeg == d); + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += (headdeg + echange == d) - (headdeg == d); + } +} + +/***************** + changestat: d_degree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_degree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex taildeg, headdeg, deg, v; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + tailattr = (int)nodeattr[tail]; + headattr = (int)nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + echange = edgestate ? -1:1; + taildeg=headdeg=-1; /* since tailattr==headattr, subtract the automatic match */ + taildeg=headdeg=0; + STEP_THROUGH_OUTEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_INEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_OUTEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + STEP_THROUGH_INEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } + } +} + + +/***************** + changestat: c_distance +*****************/ +/* +Utility function to calculate distances on a sphere of radius r, between two points +in spherical coordinates. Our notation betrays our assumption that these are +lat/lon points on the geosphere, but technically you could use any sphere you +wanted, by passing an alternative radius (so long as your coordinates were given +in angular units). +*/ +double spheredist(double lat0, double lon0, double lat1, double lon1, double r) { + double rlat0,rlat1,rlon0,rlon1,cosrl0,cosrl1,cosdl,sinrl0,sinrl1,sindl; + + rlat0=lat0/180.0*M_PI; + rlat1=lat1/180.0*M_PI; + rlon0=lon0/180.0*M_PI; + rlon1=lon1/180.0*M_PI; + cosrl0=cos(rlat0); + cosrl1=cos(rlat1); + cosdl=cos(rlon0-rlon1); + sinrl0=sin(rlat0); + sinrl1=sin(rlat1); + sindl=sin(rlon0-rlon1); + return r*atan2(sqrt((cosrl0*sindl) * (cosrl0*sindl) + (cosrl1*sinrl0-sinrl1*cosrl0*cosdl) * (cosrl1*sinrl0-sinrl1*cosrl0*cosdl)), sinrl0*sinrl1+cosrl0*cosrl1*cosdl); +} + +/*Actual distance changestat function itself (calls spheredist)*/ +C_CHANGESTAT_FN(c_distance) { + Vertex t,h; + int j; + int nv,dim,logd,sphd; + double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis,pw; + + /*Set things up*/ + nv = (int)INPUT_PARAM[0]; /*Number of vertices*/ + dim = (int)INPUT_PARAM[1]; /*Dimension of coordinate space*/ + dexp = INPUT_PARAM[2]; /*Minkowski exponent*/ + logd = (int)INPUT_PARAM[3]; /*Should we use the log distance?*/ + sphd = (int)INPUT_PARAM[4]; /*Should we use lat/lon distances on the geosphere?*/ + sphr = INPUT_PARAM[5]; /*Radius to use for spherical coordinates*/ + logmind = INPUT_PARAM[6]; /*Distance minimum for the log case*/ + logdoff = INPUT_PARAM[7]; /*Distance offset for the log case*/ + scale = INPUT_PARAM[8]; /*Scaling adjustment for raw distance*/ + pw = INPUT_PARAM[9]; /*Power to which distance should be raised*/ + coord = INPUT_PARAM+10; /*Pointer to the coordinate matrix*/ + + /*Compute the changescore*/ + t=tail-1; /*0-indexed values, for C convenience*/ + h=head-1; + if(sphd){ /*Compute great circle distances on the sphere*/ + dis=spheredist(coord[t],coord[t+nv],coord[h],coord[h+nv],sphr); + }else{ + dis=0.0; + } + /*Compute Minkowski distances in free space (adding extra dimensions to + great circle distance if sphd and we were given dim>2)*/ + mdis=0.0; + for(j=0+2*sphd;j 0){ + nrow = (N_NODES)-(long int)(INPUT_PARAM[0]); + }else{ + nrow = (long int)(INPUT_PARAM[0]); + } + +/* Rprintf("nrow %d noffset %d\n",nrow, noffset); + Rprintf("attrib: "); + for(i=0;i<1000;i++) + Rprintf("%1.0f",INPUT_ATTRIB[i]); + + Rprintf("\n;"); */ + + if(DIRECTED){ + /* directed version */ + /*Get the initial state of the edge and its reflection*/ + refedgestate = (IS_OUTEDGE(head, tail)); + + /* *** don't forget tail -> head */ + + /*Get the dyadic covariate*/ + /* val = INPUT_ATTRIB[(head-1-nrow)+(tail-1)*ncols]; */ + index = (head-1-noffset)*nrow+(tail-1); + if(index >= 0 && index <= nrow*nrow){ + val = INPUT_ATTRIB[(head-1-noffset)*nrow+(tail-1)]; + /* Rprintf("tail %d head %d nrow %d ncols %d val %f\n",tail, head, nrow, ncols, val); */ + + /*Update the change statistics, as appropriate*/ + if(refedgestate){ /* Reflected edge is present */ + if(edgestate){ /* Toggled edge _was_ present */ + if(head>tail){ /* Mut to low->high */ + CHANGE_STAT[0] -= val; + CHANGE_STAT[1] += val; + }else{ /* Mut to high->low */ + CHANGE_STAT[0] -= val; + CHANGE_STAT[2] += val; + } + }else{ /* Toggled edge _was not_ present */ + if(head>tail){ /* Low->high to mut */ + CHANGE_STAT[1] -= val; + CHANGE_STAT[0] += val; + }else{ /* High->low to mut */ + CHANGE_STAT[2] -= val; + CHANGE_STAT[0] += val; + } + } + }else{ /* Reflected edge is absent */ + if(edgestate){ /* Toggled edge _was_ present */ + if(head>tail){ /* High->low to null */ + CHANGE_STAT[2] -= val; + }else{ /* Low->high to null */ + CHANGE_STAT[1] -= val; + } + }else{ /* Toggled edge _was not_ present */ + if(head>tail){ /* Null to high->low */ + CHANGE_STAT[2] += val; + }else{ /* Null to low->high */ + CHANGE_STAT[1] += val; + } + } + } + } +}else{ + /* undirected case (including bipartite) */ + + /* *** don't forget tail -> head */ + /*Get the initial edge state*/ + /*Get the covariate value*/ + /* val = INPUT_ATTRIB[(head-1-nrow)+(tail-1)*ncols]; */ + index = (head-1-noffset)*nrow+(tail-1); + if(index >= 0 && index <= nrow*((long int)(INPUT_PARAM[0]))){ + val = INPUT_ATTRIB[(head-1-noffset)*nrow+(tail-1)]; + /*Update the change statistic, based on the toggle type*/ + /* Rprintf("tail %d head %d nrow %d noffset %d val %f\n",tail, head, nrow, noffset, val); */ + /*Update the change statistic, based on the toggle type*/ + CHANGE_STAT[0] += edgestate ? -val : val; + } + } +} + + +/******************** changestats: G ***********/ + +/***************** + changestat: d_gwdegree +*****************/ + +#define GWD0(d0) ((decay) ? exp(loneexpd*(d0)) : (d0)==0) + +C_CHANGESTAT_FN(c_gwdegree) { + int echange=0; + double decay, loneexpd, change; + Vertex taild, headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + change = 0.0; + echange = edgestate ? -1:+1; + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + change += echange*(GWD0(taild) + GWD0(headd)); + + CHANGE_STAT[0] = change; + +} + +/***************** + changestat: d_gwdegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwdegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 200?) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int tailattr, headattr, echange=0; + double decay, loneexpd; + Vertex taild, headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:+1; + taild = DEG(tail) - edgestate; + tailattr = INPUT_PARAM[tail]; + CHANGE_STAT[tailattr-1] += echange*GWD0(taild); + + headd = DEG(head) - edgestate; + headattr = INPUT_PARAM[head]; + CHANGE_STAT[headattr-1] += echange*GWD0(headd); + +} + +/***************** + changestat: d_gwidegree +*****************/ +C_CHANGESTAT_FN(c_gwidegree) { + double decay, loneexpd, change; + Vertex headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + change = 0.0; + + /* *** don't forget tail -> head */ + headd = IN_DEG[head] - edgestate; + change += (edgestate? -1.0 : 1.0) * GWD0(headd); + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_gwidegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwidegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 2008) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int headattr, echange; + double decay, loneexpd; + Vertex headd; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + headd = IN_DEG[head] - edgestate; + headattr = INPUT_PARAM[head - BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + CHANGE_STAT[headattr-1] += echange*GWD0(headd); +} + +/***************** + changestat: d_gwodegree +*****************/ +C_CHANGESTAT_FN(c_gwodegree) { + double decay, loneexpd, change; + Vertex taild; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + change = 0.0; + + /* *** don't forget tail -> head */ + taild = OUT_DEG[tail] - edgestate; + change += (edgestate? -1 : 1) * GWD0(taild); + CHANGE_STAT[0] = change; +} + +/***************** + changestat: d_gwodegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwodegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 2008) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int tailattr, echange; + double decay, loneexpd; + Vertex taild; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = OUT_DEG[tail] - edgestate; + tailattr = INPUT_PARAM[tail]; + CHANGE_STAT[tailattr-1] += echange*GWD0(taild); +} + +/******************** changestats: H ***********/ +/***************** + changestat: d_hamming +*****************/ +/* This function must be passed two networks in the forms of edgelists: + One is the network from which hamming distances are calculated and the + other is the network of weights on the dyads, which is an edgelist with a + third column of weights. Note that all non-edges in the second network + will have the value given by defaultval; thus, the "unweighted" hamming + distance is obtained when the default is 1.0 and the second network is + empty. */ +C_CHANGESTAT_FN(c_hamming) { + int discord; + + Edge wt_net_start= INPUT_PARAM[0]*2+2; + double defaultval = INPUT_PARAM[wt_net_start-1]; /* Hamming wt for non-edges in cov nwp */ + double *wt_net = INPUT_PARAM+wt_net_start; + + /* *** don't forget tail -> head */ + + discord = XOR(dEdgeListSearch(tail, head, INPUT_PARAM), edgestate); + + /* Second, search second network to see if the weight is different from + defaultval. In unweighted case, this network is empty. */ + + Edge wt_pos = dEdgeListSearch(tail, head, wt_net); + double val = wt_pos ? wt_net[wt_pos+2*(unsigned int)wt_net[0]] : defaultval; + + CHANGE_STAT[0] += (discord ? -val : val); + +} + +/******************** changestats: I ***********/ + +// A macro indicating whether x is in [from,to) +#define FROM_TO(x, from, to) ((x)>=(from) && (x)<(to)) + +/***************** + changestat: d_idegrange +*****************/ +C_CHANGESTAT_FN(c_idegrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex headideg = IN_DEG[head]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } +} + +/***************** + changestat: d_idegrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_idegrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex headideg = IN_DEG[head]; + int headattr = INPUT_PARAM[3*N_CHANGE_STATS + head - 1 - BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + for(j = 0; j < N_CHANGE_STATS; j++){ + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (headattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } +} + +/***************** + changestat: d_idegrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_idegrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of idegrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (headattr == tailattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + Vertex headideg=0, v; + STEP_THROUGH_INEDGES(head, e, v) { headideg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_idegree +*****************/ +C_CHANGESTAT_FN(c_idegree) { + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1 : +1; + Vertex headd = IN_DEG[head]; + + for(j=0; j < N_CHANGE_STATS; j++){ + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] += (headd + echange == deg) - (headd == deg); + } +} + + +/***************** + changestat: d_idegdist +*****************/ +C_CHANGESTAT_FN(c_idegdist) { + int echange = edgestate ? -1 : +1; + Vertex ohd = IN_DEG[head], + nhd = ohd + echange; + + if(nhd > N_CHANGE_STATS) cutoff_error(mtp); + + if(ohd) CHANGE_STAT[ohd-1]--; + if(nhd) CHANGE_STAT[nhd-1]++; +} + + +/***************** + changestat: d_idegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_idegree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, headattr, testattr; + Vertex headdeg, d; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1 : +1; + headdeg = IN_DEG[head]; + headattr = INPUT_PARAM[2*N_CHANGE_STATS + head - 1- BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += (headdeg + echange == d) - (headdeg == d); + } +} + +/***************** + changestat: d_idegree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_idegree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex headdeg, deg, tmp; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + tailattr = (int)nodeattr[tail]; + headattr = (int)nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + echange=edgestate ? -1 : +1; + headdeg=0; + STEP_THROUGH_INEDGES(head, e, tmp){ + headdeg += (nodeattr[tmp]==headattr); + } + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } + } +} + +/***************** + changestat: d_idegreepopularity +*****************/ +C_CHANGESTAT_FN(c_idegreepopularity) { + double change; + Vertex deg=0; + + /* *** don't forget tail -> head */ + change = 0.0; + deg = (double)(IN_DEG[head]); + if(edgestate){ + change -= sqrt(deg); + change += (deg-1.0)*(sqrt(deg-1.0)-sqrt(deg)); + }else{ + change += sqrt(deg+1.0); + change += deg*(sqrt(deg+1.0)-sqrt(deg)); + } + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_intransitive +*****************/ +C_CHANGESTAT_FN(c_intransitive) { + Edge e; + Vertex node2; + double change; + + /* *** don't forget tail -> head */ + change = 0.0; + STEP_THROUGH_OUTEDGES(head, e, node2) { + if (node2 != tail){ + if (!IS_OUTEDGE(tail,node2)){ + change = change + 1.0; + } + } + } + STEP_THROUGH_INEDGES(head, e, node2) { + if (node2 != tail){ + if (IS_OUTEDGE(tail, node2)){ + change = change - 1.0; + } + } + } + STEP_THROUGH_INEDGES(tail, e, node2) { + if (node2 != head){ + if (!IS_OUTEDGE(node2,head)){ + change = change + 1.0; + } + } + } + CHANGE_STAT[0] += edgestate ? -change : change; +/* Rprintf("tail %d head %d edgestate %d change %f\n",tail,head, change); */ +} + +/***************** +changestat: d_isolatededges +*****************/ +D_CHANGESTAT_FN(d_isolatededges) { + int i, edgestate; + Vertex tail, head, neighbor, taild, headd; + Edge e; + + /* *** don't forget tail -> head */ + ZERO_ALL_CHANGESTATS(i); + FOR_EACH_TOGGLE(i) { + // is there an edge tail -> head? + edgestate = IS_OUTEDGE(tail=TAIL(i), head=HEAD(i)); + + taild = DEG(tail); + headd = DEG(head); + + if(edgestate) { // we are removing an edge + + // if head and tail both have degree one, then + // we are removing an isolated edge + if(taild == 1 && headd == 1) + CHANGE_STAT[0] -= 1; + + // if tail has degree 2 and has a degree one node other than head as a neighbor, + // then we are making a non-isolated edge into an isolated edge by removing + // the edge tail -> head + if(taild == 2) { + STEP_THROUGH_OUTEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != head) + CHANGE_STAT[0] += 1; + } + STEP_THROUGH_INEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != head) + CHANGE_STAT[0] += 1; + } + } + + // ditto head + if(headd == 2) { + STEP_THROUGH_OUTEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != tail) + CHANGE_STAT[0] += 1; + } + STEP_THROUGH_INEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != tail) + CHANGE_STAT[0] += 1; + } + } + } else { // we are adding an edge + + // if head and tail both have degree zero, then + // we are adding an isolated edge + if(taild == 0 && headd == 0) + CHANGE_STAT[0] += 1; + + // if tail has degree 1 and so does its neighbor, then we + // are making an isolated edge into a non-isolated edge; + // note that for undirected graphs, this neighbor cannot + // be head, as the current toggle is to turn on the edge + // tail -> head + if(taild == 1) { + STEP_THROUGH_OUTEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + STEP_THROUGH_INEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + } + + // ditto head + if(headd == 1) { + STEP_THROUGH_OUTEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + STEP_THROUGH_INEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + } + } + + TOGGLE_IF_MORE_TO_COME(i); + } + + UNDO_PREVIOUS_TOGGLES(i); +} + +/***************** + changestat: d_isolates +*****************/ +C_CHANGESTAT_FN(c_isolates) { + int echange; + Vertex taild, headd=0; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:+1; + taild = DEG(tail); + headd = DEG(head); + CHANGE_STAT[0] += (taild + echange == 0) - (taild == 0); + CHANGE_STAT[0] += (headd + echange == 0) - (headd == 0); + +} + +S_CHANGESTAT_FN(s_isolates) { + /* *** don't forget tail -> head */ + CHANGE_STAT[0] = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++){ + if(DEG(tail) == 0) + CHANGE_STAT[0] ++; + } +} + +/***************** + changestat: d_istar +*****************/ +C_CHANGESTAT_FN(c_istar) { + double change, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double tailattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + headd = -(int)edgestate; + STEP_THROUGH_INEDGES(head, e, node3) {/* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headd = IN_DEG[head] - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/******************** changestats: K ***********/ +/***************** + changestat: d_kstar +*****************/ +C_CHANGESTAT_FN(c_kstar) { + double change, taild, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double tailattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + taild = -(int)edgestate; + STEP_THROUGH_OUTEDGES(tail, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++taild;} + } + STEP_THROUGH_INEDGES(tail, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++taild;} + } + headd = -(int)edgestate; + STEP_THROUGH_OUTEDGES(head, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + STEP_THROUGH_INEDGES(head, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; +/* if (kmo==0) { + change=1; + } else { */ + change = CHOOSE(taild, kmo) + CHOOSE(headd, kmo); +/* } uncomment these few lines to define 1-stars as equivalent to + edges (currently, each edge is counted as two 1-stars) */ + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) + { + kmo = ((int)INPUT_PARAM[j]) - 1; +/* if (kmo==0) { + change=1; + } else { */ + change = CHOOSE(taild, kmo) + CHOOSE(headd, kmo); +/* } uncomment these few lines to define 1-stars as equivalent to + edges (currently, each edge is counted as two 1-stars) */ + CHANGE_STAT[j] += (edgestate ? - change : change); + } + + } +} + + +/******************** changestats: L ***********/ +/***************** + changestat: d_localtriangle +*****************/ +C_CHANGESTAT_FN(c_localtriangle) { + Edge e; + Vertex node3, nmat; + double change; + + nmat = (Vertex)(INPUT_PARAM[0]); + + /* *** don't forget tail -> head */ + change = 0.0; + + if(INPUT_PARAM[1+(head-1)+(tail-1)*nmat] == 1.0){ + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(INPUT_PARAM[1+(node3-1)+(tail-1)*nmat] == 1.0 && + INPUT_PARAM[1+(node3-1)+(head-1)*nmat] == 1.0 ){ + if (DIRECTED){ + if (IS_INEDGE(node3,tail) ) ++change; + if (IS_OUTEDGE(node3,tail)) ++change; + }else{ + if (IS_UNDIRECTED_EDGE(node3,tail)) ++change; + } + } + } + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(INPUT_PARAM[1+(node3-1)+(tail-1)*nmat] == 1.0 && + INPUT_PARAM[1+(node3-1)+(head-1)*nmat] == 1.0 ){ + if (DIRECTED) + { + if (IS_INEDGE(node3,tail) ) ++change; + if (IS_OUTEDGE(node3,tail)) ++change; + } + else + { + if (IS_UNDIRECTED_EDGE(node3,tail)) ++change; + } + } + } + + CHANGE_STAT[0] += edgestate ? - change : change; + + } +} + +/******************** changestats: M ***********/ +/***************** + changestat: d_m2star +*****************/ +C_CHANGESTAT_FN(c_m2star) { + int tailid, headod, change; + int backedgestate; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if the edge from tail to head */ + /* exists and will disappear */ + /* edgestate is 0 if the edge does not exist */ + backedgestate = (IS_OUTEDGE(head, tail)); + + tailid = IN_DEG[tail]; + headod = OUT_DEG[head]; + change = tailid + headod - 2*backedgestate; + CHANGE_STAT[0] += (edgestate ? -change : change); + +} + +/***************** + changestat: d_mutual + + (1,1) -> anything = -1 + anything -> (1,1) = +1 +*****************/ +C_CHANGESTAT_FN(c_mutual) { + double matchval, change; + int j, ninputs, noattr; + + ninputs = N_INPUT_PARAMS - N_NODES; + noattr = (N_INPUT_PARAMS == 0); + + /* *** don't forget tail -> head */ + if (IS_OUTEDGE(head,tail)) { /* otherwise, no change occurs */ + change = edgestate ? -1.0 : 1.0 ; + if (noattr) { /* "plain vanilla" mutual, without node attributes */ + CHANGE_STAT[0] += change; + } else { /* Only consider mutuals where node attributes match */ + matchval = INPUT_PARAM[tail+ninputs-1]; + if (matchval == INPUT_PARAM[head+ninputs-1]) { /* We have a match! */ + if (ninputs==0) {/* diff=F in network statistic specification */ + CHANGE_STAT[0] += change; + } else { /* diff=T */ + for (j=0; j head */ + if (IS_OUTEDGE(head,tail)) { /* otherwise, no change occurs */ + change = edgestate ? -1.0 : 1.0 ; + for (j=0; j head */ + edgestateth = (!IS_OUTEDGE(head,tail)); + + for(node3=1;node3<=N_NODES;node3++){ + if((node3!=tail)&&(node3!=head)){ + sc = edgestateth + (!IS_OUTEDGE(node3,tail)); + if(sc < 2){ + sc += (!IS_OUTEDGE(tail,node3)); + if(sc < 2){ + sc += (!IS_OUTEDGE(node3,head)); + if(sc < 2){ + sc += (!IS_OUTEDGE(head,node3)); + if(sc < 2){ + change=0.0; + if (sc == 0 && edgestate == 0 ){--change;} + if (sc == 0 && edgestate == 1 ){++change;} + if (sc == 1 && edgestate == 0 ){++change;} + if (sc == 1 && edgestate == 1 ){--change;} + CHANGE_STAT[0] += change; + } + } + } + } + } + } + +} + +/******************** changestats: O ***********/ + +// A macro indicating whether x is in [from,to) +#define FROM_TO(x, from, to) ((x)>=(from) && (x)<(to)) + +/***************** + changestat: d_odegrange +*****************/ +C_CHANGESTAT_FN(c_odegrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex tailodeg = OUT_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } +} + +/***************** + changestat: d_odegrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_odegrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex tailodeg = OUT_DEG[tail]; + int tailattr = INPUT_PARAM[3*N_CHANGE_STATS + tail - 1]; + for(j = 0; j < N_CHANGE_STATS; j++){ + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } +} + +/***************** + changestat: d_odegrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_odegrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of odegrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + Vertex tailodeg=0, v; + STEP_THROUGH_OUTEDGES(tail, e, v) { tailodeg += (nodeattr[v]==tailattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_odegree +*****************/ +C_CHANGESTAT_FN(c_odegree) { + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1 : 1; + Vertex taild = OUT_DEG[tail]; + + for(j=0; j < N_CHANGE_STATS; j++) { + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] = (taild + echange == deg) - (taild == deg); + } +} + + +/***************** + changestat: d_odegdist +*****************/ +C_CHANGESTAT_FN(c_odegdist) { + int echange = edgestate ? -1 : +1; + Vertex otd = OUT_DEG[tail], + ntd = otd + echange; + + if(ntd > N_CHANGE_STATS) cutoff_error(mtp); + + if(otd) CHANGE_STAT[otd-1]--; + if(ntd) CHANGE_STAT[ntd-1]++; +} + + +/***************** + changestat: d_odegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_odegree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, testattr; + Vertex taildeg, d; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1 : +1; + taildeg = OUT_DEG[tail]; + tailattr = INPUT_PARAM[2*N_CHANGE_STATS + tail - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (tailattr == testattr) { /* we have tail attr match */ + CHANGE_STAT[j] += (taildeg + echange == d) - (taildeg == d); + } + } +} + +/***************** + changestat: d_odegree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_odegree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange=edgestate ? -1 : +1; + Vertex taildeg=0, tmp; + STEP_THROUGH_OUTEDGES(tail, e, tmp){ + taildeg += (nodeattr[tmp]==tailattr); + } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + } + } +} + +/***************** + changestat: d_opentriad +*****************/ +C_CHANGESTAT_FN(c_opentriad) { + + /* *** don't forget tail -> head */ + Vertex node3; + Edge change = 0, e; + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + + // -3 * triangles + + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + change += IS_UNDIRECTED_EDGE(node3,tail); + } + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + change += IS_UNDIRECTED_EDGE(node3,tail); + } + CHANGE_STAT[0] += change * (edgestate ? 3.0 : -3.0); + + + // +1 * 2-stars + + Vertex taild = DEG(tail) - edgestate; + Vertex headd = DEG(head) - edgestate; + change = taild + headd; + CHANGE_STAT[0] += (edgestate ? -change : change); + +} + +/***************** + changestat: d_ostar +*****************/ +C_CHANGESTAT_FN(c_ostar) { + double change, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double headattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headattr = INPUT_ATTRIB[head-1]; + if(headattr == INPUT_ATTRIB[tail-1]){ + headd = -(int)edgestate; + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of head */ + if(headattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headd = OUT_DEG[tail] - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_odegreepopularity +*****************/ +C_CHANGESTAT_FN(c_odegreepopularity) { + double change; + Vertex deg=0; + + /* *** don't forget tail -> head */ + change = 0.0; + deg = (double)(OUT_DEG[tail]); + if(edgestate){ + change -= sqrt(deg); + change += (deg-1.0)*(sqrt(deg-1.0)-sqrt(deg)); + }else{ + change += sqrt(deg+1.0); + change += deg*(sqrt(deg+1.0)-sqrt(deg)); + } +CHANGE_STAT[0]=change; +} + +/******************** changestats: P ***********/ +/***************** + changestat: d_pdegcor +*****************/ +D_CHANGESTAT_FN(d_pdegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + mtp->dstats[0] -= current; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_pdegcor) { + Vertex taildeg, headdeg; + Edge e; + double mu, mu2, mutail, mutail2, sigma2, sigmatail2, cross; + + mu = 0.0; + mu2 = 0.0; + mutail = 0.0; + mutail2 = 0.0; + cross = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + mu += (double)(headdeg); + mutail += (double)(taildeg); + mu2 += (double)(headdeg*headdeg); + mutail2 += (double)(taildeg*taildeg); + cross += taildeg*headdeg; + } + } + mu = mu / (N_EDGES); + mutail = mutail / (N_EDGES); + sigma2 = mu2/(N_EDGES) - mu*mu; + sigmatail2 = mutail2/(N_EDGES) - mutail*mutail; + CHANGE_STAT[0] = (cross / (N_EDGES) - mutail*mu) / sqrt(sigma2*sigmatail2); +} + +/******************** changestats: R ***********/ +/***************** + changestat: d_rdegcor +*****************/ +D_CHANGESTAT_FN(d_rdegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ +// CHANGE_STAT[0] = mtp->dstats[0] - current; +// Rprintf("c %f p %f",current,mtp->dstats[0]); + mtp->dstats[0] -= current; +// Rprintf(" p-c %f\n",mtp->dstats[0]); + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_rdegcor) { + Vertex taildeg, headdeg; + Edge e; + double mu, mu2, sigma2, cross; + Vertex tailrank, headrank; + Vertex *ndeg=R_Calloc(N_NODES+1, Vertex); + + for(Vertex tail=0; tail <= N_NODES; tail++) { ndeg[tail]=0; } + for(Vertex tail=0; tail < N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + ndeg[taildeg+1]++; + ndeg[headdeg+1]++; + } + } +for(Vertex tail=1; tail <= N_NODES; tail++) { + ndeg[tail] += ndeg[tail-1]; +} +// Rprintf("tail %d taildeg[tail] %d \n",tail,ndeg[tail]);} + + mu = 0.0; + mu2 = 0.0; + cross = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + tailrank = (ndeg[taildeg+1]+ndeg[taildeg+2]+1)*0.5; + headrank = (ndeg[headdeg+1]+ndeg[headdeg+2]+1)*0.5; + mu += (double)(tailrank + headrank); + mu2 += (double)(tailrank*tailrank + headrank*headrank); + cross += 2.0*tailrank*headrank; + } + } + mu = mu / (2.0*N_EDGES); + sigma2 = mu2/(2.0*N_EDGES) - mu*mu; + CHANGE_STAT[0] = (cross / (2.0*N_EDGES) - mu*mu) / sigma2; + R_Free(ndeg); +} + +/***************** + changestat: d_simmelian +*****************/ +C_CHANGESTAT_FN(c_simmelian) { + Edge e; + Vertex change, node3; + + /* *** don't forget tail -> head */ + + if(IS_OUTEDGE(head, tail)){ + change = 0; + + STEP_THROUGH_OUTEDGES(head, e, node3) /* step through outedges of head */ + { + if (node3 != tail + && IS_OUTEDGE(node3, tail) + && IS_OUTEDGE(tail, node3) + && IS_OUTEDGE(node3, head) + ){++change;} + } + + CHANGE_STAT[0] += edgestate ? -(double)change : (double)change; + } + +} + +/***************** + changestat: d_simmelianties +*****************/ +C_CHANGESTAT_FN(c_simmelianties) { + Edge e, e2; + Vertex change, node3, node4, first, htflag; + + /* *** don't forget tail -> head */ + + if(IS_OUTEDGE(head, tail)){ + change = htflag = 0; + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if (node3 != tail + && IS_OUTEDGE(node3, tail) && IS_OUTEDGE(tail, node3) && IS_OUTEDGE(node3, head)){ + htflag=1; /* tail, head is itself in a simmelian triple (along with head, tail)*/ + first = 1; + /* Find out whether (tail, node3) is in any other simmelian triple */ + STEP_THROUGH_OUTEDGES (tail, e2, node4) { /* step through outedges of tail */ + if (node4 != head && node4 != node3 && IS_OUTEDGE(node4, tail) + && IS_OUTEDGE(node4, node3) && IS_OUTEDGE(node3, node4)){ + first = 0; + } + } + if(first){++change;} + first = 1; + /* Find out whether (head, node3) is in any other simmelian triple */ + STEP_THROUGH_OUTEDGES (head, e2, node4) { /* step through outedges of head */ + if (node4 != tail && node4 != node3 && IS_OUTEDGE(node4, head) + && IS_OUTEDGE(node4, node3) && IS_OUTEDGE(node3, node4)) { + first = 0; + } + } + if(first){++change;} + } + } + change += htflag; + change = 2*change; /* All changes must happen in pairs here; no tie can + be counted without its opposite */ + CHANGE_STAT[0] += edgestate ? -(double)change : (double)change; + } +} + +/******************** changestats: T ***********/ + +/***************** + changestat: d_threetrail +*****************/ +C_CHANGESTAT_FN(c_threetrail) { + int j, k, change, dchange[4]; + Edge e; + Vertex node3; + /* The four values of dchange represent the four different types of + directed threetrails oriented so that the middle step is always + "right" (R). In order: RRR, RRL, LRR, LRL + i.e., >>> >>< <>> <>< */ + + + /* *** don't forget tail -> head */ + /* Step A: Count threetrails in which tail->head is the middle edge */ + dchange[0] = IN_DEG[tail] * OUT_DEG[head]; /* R then R; may count head->tail->head->tail */ + dchange[1] = IN_DEG[tail] * (IN_DEG[head]-edgestate); /* R then L */ + dchange[2] = (OUT_DEG[tail]-edgestate) * OUT_DEG[head]; /* L then R */ + dchange[3] = (OUT_DEG[tail]-edgestate) * (IN_DEG[head]-edgestate); /* L then L */ + /* Step B: Count threetrails where tail is one endpoint */ + STEP_THROUGH_OUTEDGES(head, e, node3) { /* tail->head->node3-x which means -RL */ + dchange[1] += IN_DEG[node3]-1; /* RRL; subtract 1 for head itself */ + dchange[0] += OUT_DEG[node3]; /* RRR; possibly counted tail->head->tail->head */ + } + STEP_THROUGH_INEDGES(head, e, node3) { /* x-node3->head<-tail which means -RL*/ + if (node3 != tail) { + dchange[3] += OUT_DEG[node3]-1; /* LRL; subtract 1 for head itself */ + dchange[1] += IN_DEG[node3]; /* RRL */ + } + } + /* Step C: Count threetrails where head is one endpoint */ + STEP_THROUGH_INEDGES(tail, e, node3) { /* x-node3->tail->head which means -RR */ + dchange[2] += OUT_DEG[node3]-1; /* LRR; subtract 1 for tail itself */ + dchange[0] += IN_DEG[node3]; /* RRR; possibly counted tail->head->tail->head */ + } + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* head<-tail->node3-x which means LR- */ + if (node3 != head) { + dchange[3] += IN_DEG[node3]-1; /* LRL; subtract 1 for tail itself */ + dchange[2] += OUT_DEG[node3]; /* LRR */ + } + } + /* Finally, correct for overcounted head->tail->head->tail and tail->head->tail->head */ + if (DIRECTED) { + dchange[0] -= IS_INEDGE(tail, head) * (1 + 2 * edgestate); + /* head->tail->head->tail is counted in A whenever IS_INEDGE(tail,head) but + TT->head->tail->head is only counted in B and C when also edgestate */ + for (j = 0; j < N_INPUT_PARAMS; j++) { + k = (int) INPUT_PARAM[j]; + CHANGE_STAT[j] += (edgestate ? -dchange[k-1] : dchange[k-1]); + } + } + else { /* Undirected case; don't need head->tail->head->tail correction */ + change = dchange[0] + dchange[1] + dchange[2] + dchange[3]; + CHANGE_STAT[0] += (edgestate ? -change : change); + } +} + +/***************** + changestat: d_transitive +*****************/ +C_CHANGESTAT_FN(c_transitive) { + Edge e; + Vertex node2; + double change; + + /* *** don't forget tail -> head */ + change = 0.0; /* change should become the number of transitive triples + a->b, b->c, a->c in which tail->head is found */ + + STEP_THROUGH_OUTEDGES(head, e, node2) { /* step through outedges of head */ + if (tail != node2 && IS_OUTEDGE(tail, node2)){ + change = change + 1.0; /* Here we have tail->head, head->node2, tail->node2 */ + } + } + STEP_THROUGH_INEDGES(head, e, node2) { /* step through inedges of head */ + if (tail != node2) { + change = change + IS_OUTEDGE(tail, node2) + IS_OUTEDGE(node2, tail); + /* Here we have tail->head and node2->head, with either node2->tail or tail->node2 */ + } + } +// STEP_THROUGH_INEDGES(tail, e, node2) { /* step through inedges of tail */ +// if (node2 != head){ +// if (!IS_OUTEDGE(node2, head)){ +// change = change - 1.0; +// } +// } +// } + CHANGE_STAT[0] += edgestate ? -change : change; +// Rprintf("tail %d head %d edgestate %d change %f C_S[0]=%f\n", tail, head, change,CHANGE_STAT[0]); +} + +C_CHANGESTAT_FN(c_transitiveties) { + int echange, ochange; + int L2th, L2tu, L2uh; + double cumchange; + double tailattr; + + + /* *** don't forget tail -> head */ + cumchange=0.0; + L2th=0; + ochange = GETWT(tail, head) ? -1 : 0; + echange = 2*ochange + 1; + if(N_INPUT_PARAMS>0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(tail, u) && (tailattr == INPUT_ATTRIB[u-1])){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v) && (tailattr == INPUT_ATTRIB[v-1])){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_INEDGES(head, e, u, { + if (GETWT(tail, u) && (tailattr == INPUT_ATTRIB[u-1])){ + L2th++; + } + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head) && (tailattr == INPUT_ATTRIB[v-1])){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + }else{ /* no attributes */ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(tail, u)){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v)){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_INEDGES(head, e, u, { + if (GETWT(tail, u)){ + L2th++; + } + if (GETWT(u, tail)){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head)){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + + cumchange += (L2th>0) ; +// Rprintf("L2th %d echange %d cumchange %f tail %d head %d\n", L2th, echange, cumchange,tail,head); + cumchange = echange*cumchange; + (CHANGE_STAT[0]) += cumchange; +} + +C_CHANGESTAT_FN(c_cyclicalties) { + int echange, ochange; + int L2th, L2tu, L2uh; + double cumchange; + double tailattr; + + + /* *** don't forget tail -> head */ + cumchange=0.0; + L2th=0; + ochange = GETWT(tail, head) ? -1 : 0; + echange = 2*ochange + 1; + if(N_INPUT_PARAMS>0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v) && (tailattr == INPUT_ATTRIB[v-1])){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2th++; + } + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head) && (tailattr == INPUT_ATTRIB[v-1])){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + }else{ /* no attributes */ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail)){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v)){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through outedges of head */ + + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail)){ + L2th++; + } + if (GETWT(u, tail)){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head)){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + + cumchange += (L2th>0) ; +// Rprintf("L2th %d echange %d cumchange %f tail %d head %d\n", L2th, echange, cumchange,tail,head); + cumchange = echange*cumchange; + (CHANGE_STAT[0]) += cumchange; +} + +/***************** + changestat: d_triadcensus +*****************/ +C_CHANGESTAT_FN(c_triadcensus) { + int j, a, b, c, d, e, edgecount, t300, + t210, t120C, t120U, t120D, t201, t030C, t030T, t111U, + t111D, t021C, t021U, t021D, t102, t012, t003; + Vertex triadtype, node3; + + /* *** don't forget tail -> head */ + if (DIRECTED) { + /* directed version */ + t300 = 0; + t210 = 0; + t120C = 0; t120U = 0; t120D = 0; t201 = 0; + t030C = 0; t030T = 0; t111U = 0; t111D = 0; + t021C = 0; t021U = 0; t021D = 0; t102 = 0; + t012 = 0; + + if ( (MIN_OUTEDGE(head) != 0) || + (MIN_INEDGE(head) != 0) || + (MIN_OUTEDGE(tail) != 0) || + (MIN_INEDGE(tail) != 0) ) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = (IS_OUTEDGE(head, tail)); + b = (IS_OUTEDGE(head, node3)); + c = (IS_OUTEDGE(node3, head)); + d = (IS_OUTEDGE(node3, tail)); + e = (IS_OUTEDGE(tail, node3)); + edgecount = (a + b + c + d + e); + + switch(edgecount) { + case 0: /* 012 */ + ++t012; + + case 1: { /* 021C, 021U, 021D, 102 */ + if ((b == 1) || (d == 1)) + ++t021C; + if (c == 1) + ++t021U; + if (e == 1) + ++t021D; + if (a == 1) + ++t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if ((b + d) == 2) + ++t030C; + if (((b + e) == 2) || ((c + d) == 2) || ((c + e) == 2)) + ++t030T; + if (((a + b) == 2) || ((a + e) == 2) || ((d + e) == 2)) + ++t111U; + if (((a + c) == 2) || ((a + d) == 2) || ((b + c) == 2)) + ++t111D; + } + break; + + case 3: { /* 120C, 120U, 120D, 201 */ + if (a == 1) { + if (((b + d) == 2) || ((c + e) == 2)) + ++t120C; + if ((b + e) == 2) + ++t120U; + if ((c + d) == 2) + ++t120D; + if (((b + c) == 2) || ((d + e) == 2)) + ++t201; + } else { + if (b == 1) { + if (((c + d) == 2) || ((d + e) == 2)) + ++t120C; + if ((c + e) == 2) + ++t120D; + } else { + ++t120U; + } + } + } + break; + + case 4: /* 210 */ + ++t210; + break; + + case 5: /* 300 */ + ++t300; + break; + } + + switch(edgecount) { + case 1: /* 102, 021D, 021U, 021C */ + --t012; + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if (((a + c) == 2) || ((a + e) == 2) || ((b + d) == 2) || + ((c + e) == 2)) + --t021C; + if (((a + d) == 2) || ((b + e) == 2)) + --t021U; + if (((a + b) == 2) || ((c + d) == 2)) + --t021D; + if (((b + c) == 2) || ((d + e) == 2)) + --t102; + } + break; + + case 3: { /* 201, 120D, 120U, 120C */ + if (a == 1) { + if ((c + e) == 2) + --t030C; + if (((c + d) == 2) || ((b + e) == 2) || ((b + d) == 2)) + --t030T; + if ((b + c) == 2) + --t111U; + if ((d + e) == 2) + --t111D; + } else { + if (b == 1) { + if ((c + d) == 2) + --t111U; + if (((c + e) == 2) || ((d + e) == 2)) + --t111D; + } else { + --t111U; + } + } + } + break; + + case 4: { /* 210 */ + if (a == 1) + { + if (((b + c + e) == 3) || ((c + d + e) == 3)) + --t120C; + if ((b + c + d) == 3) + --t120U; + if ((b + d + e) == 3) + --t120D; + } else { + if ((b + c + d + e) == 4) + --t201; + } + } + break; + + case 5: /* 300 */ + --t210; + break; + } + } + } /* ****** move to next node3 ******** */ + }else{ + t012 = t012 + (N_NODES - 2); + } + + for(j = 0; j < N_CHANGE_STATS; j++) { + triadtype = (Vertex)INPUT_PARAM[j]; + + switch(triadtype) { /* SEARCH_ON_THIS_TO_TRACK_DOWN_TRIADCENSUS_CHANGE + to undo triadcensus change, change - to plus in + next two lines: */ + case 1: t003 = -(t300+t210+t120C+t120U+t120D+t201+t030C+t030T); + t003 = t003-(t111U+t111D+t021C+t021U+t021D+t102+t012); + CHANGE_STAT[j] += edgestate ? -(double)t003 : (double)t003; + break; + case 2: CHANGE_STAT[j] += edgestate ? -(double)t012 : (double)t012; + break; + case 3: CHANGE_STAT[j] += edgestate ? -(double)t102 : (double)t102; + break; + case 4: CHANGE_STAT[j] += edgestate ? -(double)t021D : (double)t021D; + break; + case 5: CHANGE_STAT[j] += edgestate ? -(double)t021U : (double)t021U; + break; + case 6: CHANGE_STAT[j] += edgestate ? -(double)t021C : (double)t021C; + break; + case 7: CHANGE_STAT[j] += edgestate ? -(double)t111D : (double)t111D; + break; + case 8: CHANGE_STAT[j] += edgestate ? -(double)t111U : (double)t111U; + break; + case 9: CHANGE_STAT[j] += edgestate ? -(double)t030T : (double)t030T; + break; + case 10: CHANGE_STAT[j] += edgestate ? -(double)t030C : (double)t030C; + break; + case 11: CHANGE_STAT[j] += edgestate ? -(double)t201 : (double)t201; + break; + case 12: CHANGE_STAT[j] += edgestate ? -(double)t120D : (double)t120D; + break; + case 13: CHANGE_STAT[j] += edgestate ? -(double)t120U : (double)t120U; + break; + case 14: CHANGE_STAT[j] += edgestate ? -(double)t120C : (double)t120C; + break; + case 15: CHANGE_STAT[j] += edgestate ? -(double)t210 : (double)t210; + break; + case 16: CHANGE_STAT[j] += edgestate ? -(double)t300 : (double)t300; + break; + } + } + } else { + /* undirected */ + + /* *** don't forget tail -> head */ + t300 = 0; t201 = 0; t102 = 0; t012 = 0; + + if ( (MIN_OUTEDGE(head) != 0) || + (MIN_INEDGE(head) != 0) || + (MIN_OUTEDGE(tail) != 0) || + (MIN_INEDGE(tail) != 0) ) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = (IS_UNDIRECTED_EDGE(node3, head)); + b = (IS_UNDIRECTED_EDGE(node3, tail)); + edgecount = (a + b); + + switch(edgecount) { + case 0: { /* 012 */ + ++t102; + --t012; + } + break; + + case 1: { /* 021C, 021U, 021D, 102 */ + ++t201; + --t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + ++t300; + --t201; + } + break; + + } + } + + } /* ****** move to next node3 ******** */ + } else { + t102 = t102 + (N_NODES - 2); + } + + for(j = 0; j < N_CHANGE_STATS; j++) { + triadtype = (Vertex)INPUT_PARAM[j]; + + switch(triadtype) { /* SEARCH_ON_THIS_TO_TRACK_DOWN_TRIADCENSUS_CHANGE + to undo triadcensus change, change - to plus in + next line: */ + case 1: t003 = -(t102+t201+t300); + CHANGE_STAT[j] += edgestate ? -(double)t003 : (double)t003; + break; + case 2: CHANGE_STAT[j] += edgestate ? -(double)t102 : (double)t102; + break; + case 3: CHANGE_STAT[j] += edgestate ? -(double)t201 : (double)t201; + break; + case 4: CHANGE_STAT[j] += edgestate ? -(double)t300 : (double)t300; + break; + } + } + } +} + + +/***************** + changestat: d_tripercent +*****************/ +C_CHANGESTAT_FN(c_tripercent) { + Edge e, e2; + Vertex node1, node2, node3; + int j; + Edge triwith, triwithout; + Edge degreewith, degreewithout, twostarwith, twostarwithout; + int ninputs = N_INPUT_PARAMS - N_NODES; + int MatchingOnAttribute = (ninputs>0); + double *attr=INPUT_PARAM, ratiowith, ratiowithout; + + if (MatchingOnAttribute) + attr = INPUT_PARAM + (ninputs-1); /* ptr to vertex attributes */ + + /* *** don't forget tail -> head */ + if (!edgestate) TOGGLE(tail, head); /* turn on the edge if it's missing */ + for (j=0; j < MAX(1, ninputs); j++) { + /* Count triangles with and without proposed edge */ + /* Simultaneously, find degree (use matching if necessary) with and without */ + triwith = triwithout = twostarwith = twostarwithout = 0; + for (node1 = 1; node1 <= N_NODES; node1++) { + degreewith = degreewithout = 0; + if (ninputs < 2 || EQUAL(attr[node1],INPUT_PARAM[j])) { + STEP_THROUGH_OUTEDGES(node1, e, node2) { + /* inside this loop, node1 < node2 always */ + if (!MatchingOnAttribute || EQUAL(attr[node1],attr[node2])) { + /* increment degree counter */ + ++degreewith; + if (node1!=tail || node2!=head) ++degreewithout; + STEP_THROUGH_OUTEDGES(node2, e2, node3) { + /* inside this loop, node1 < node2 < node3 always */ + if (!MatchingOnAttribute || EQUAL(attr[node2],attr[node3])) { + if (IS_OUTEDGE(node1, node3)) { + ++triwith; + if ((tail!=node1||head!=node2)&&(tail!=node2||head!=node3)&&(tail!=node1||head!=node3)) + ++triwithout; + } + } + } + } + } + STEP_THROUGH_INEDGES(node1, e, node2) { + /* inside this loop, node2 < node1 always. */ + /* We only do this to find correct degree for node1; */ + /* for triangles, node1 < node2 head */ + edgemult = edgestate ? -1.0 : 1.0; + change = 0; + if(N_INPUT_PARAMS > 0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]) { + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_INEDGE(node3, tail); + } + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_OUTEDGE(node3, tail) + IS_INEDGE(node3, tail); + } + if(N_CHANGE_STATS > 1) { /* diff = TRUE; matches must be tabled */ + for (j=0; j= d) - (b1deg >= d); + } +} + + +/***************** + changestat: c_b2mindegree +*****************/ +C_CHANGESTAT_FN(c_b2mindegree) { + int echange = edgestate ? -1 : 1; + Vertex b2deg = IN_DEG[head]; + for(unsigned int j = 0; j < N_CHANGE_STATS; j++) { + Vertex d = INPUT_PARAM[j]; + CHANGE_STAT[j] += (b2deg + echange >= d) - (b2deg >= d); + } +} From 6784e2ffc768af951935be6a58bb2be1f1bca73e Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:35:38 -0400 Subject: [PATCH 08/18] --Set up 'changestat_spheredist.c' to be appended to the 'changestat.c' placeholder. --- src/{changestat_spheredist.c => concat-changestat.c} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestat_spheredist.c => concat-changestat.c} (100%) diff --git a/src/changestat_spheredist.c b/src/concat-changestat.c similarity index 100% rename from src/changestat_spheredist.c rename to src/concat-changestat.c From e372b2998abcaf41ae034e8c7f52bf9eaffa3642 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:35:38 -0400 Subject: [PATCH 09/18] --Set up 'changestat.c' to be appended to the 'changestat.c' placeholder. --- src/{changestat.c => concat-changestat.c} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestat.c => concat-changestat.c} (100%) diff --git a/src/changestat.c b/src/concat-changestat.c similarity index 100% rename from src/changestat.c rename to src/concat-changestat.c From 3b263a7559b32458846690bb4859afc97d668af0 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:36:00 -0400 Subject: [PATCH 10/18] Concatenate files 'changestat.c' 'changestat_spheredist.c' into 'changestat.c'. --- src/{concat-changestat.c => changestat.c} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{concat-changestat.c => changestat.c} (100%) diff --git a/src/concat-changestat.c b/src/changestat.c similarity index 100% rename from src/concat-changestat.c rename to src/changestat.c From e450f7d5018f06da988e9083e33eaed450738d6e Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:37:17 -0400 Subject: [PATCH 11/18] --Duplicate changestats.c history into changestats_distance.c --- src/{changestats.c => changestats_distance.c} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestats.c => changestats_distance.c} (100%) diff --git a/src/changestats.c b/src/changestats_distance.c similarity index 100% rename from src/changestats.c rename to src/changestats_distance.c From fc3a868b8da5d742d5a66c68f8c52bbb63404d43 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:37:17 -0400 Subject: [PATCH 12/18] --Restore changestats.c --- src/changestats.c | 3491 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 3491 insertions(+) create mode 100644 src/changestats.c diff --git a/src/changestats.c b/src/changestats.c new file mode 100644 index 000000000..8a0297e0f --- /dev/null +++ b/src/changestats.c @@ -0,0 +1,3491 @@ +/* File src/changestats.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-2026 Statnet Commons + */ +#include "changestats.h" +#include "ergm_storage.h" +#include "ergm_dyad_hashmap.h" +#include "ergm_edgelist.h" + +/******************** changestats: A ***********/ + +/***************** + changestat: d_adegcor +*****************/ +D_CHANGESTAT_FN(d_adegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i),HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ +// CHANGE_STAT[0] = mtp->dstats[0] - current; +// Rprintf("c %f p %f",current,mtp->dstats[0]); + mtp->dstats[0] -= current; +// Rprintf(" p-c %f\n",mtp->dstats[0]); +FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_adegcor) { + Vertex tail, head, taildeg, headdeg; + Edge e; + double mu, mu2, sigma2, cross; + + mu = 0.0; + mu2 = 0.0; + cross = 0.0; + for(tail=1; tail <= N_NODES; tail++) { + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + mu += (double)(taildeg + headdeg); + mu2 += (double)(taildeg*taildeg + headdeg*headdeg); + cross += 2.0*taildeg*headdeg; + } + } + mu = mu / (2.0*N_EDGES); + sigma2 = mu2/(2.0*N_EDGES) - mu*mu; + CHANGE_STAT[0] = (cross / (2.0*N_EDGES) - mu*mu) / sigma2; +} + +/***************** + changestat: d_altkstar +*****************/ +C_CHANGESTAT_FN(c_altkstar) { + double lambda, oneexpl, change; + Vertex taild, headd=0; + + change = 0.0; + lambda = INPUT_PARAM[0]; + oneexpl = 1.0-1.0/lambda; + + /* *** don't forget tail -> head */ + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + if(taild!=0){ + change += (edgestate?-1:+1)*(1.0-pow(oneexpl,(double)taild)); + } + if(headd!=0){ + change += (edgestate?-1:+1)*(1.0-pow(oneexpl,(double)headd)); + } + CHANGE_STAT[0] = change*lambda; +} + +/***************** + changestat: d_asymmetric +*****************/ +C_CHANGESTAT_FN(c_asymmetric) { + double matchval, change; + int j, ninputs, noattr; + + ninputs = N_INPUT_PARAMS - N_NODES; + noattr = (N_INPUT_PARAMS == 0); + + /* *** don't forget tail -> head */ + change = (edgestate==IS_OUTEDGE(head, tail) ? 1.0 : -1.0) ; + if (noattr) { /* "plain vanilla" asymmetric, without node attributes */ + CHANGE_STAT[0] += change; + } else { /* Only consider asymmetrics where node attributes match */ + matchval = INPUT_PARAM[tail+ninputs-1]; + if (matchval == INPUT_PARAM[head+ninputs-1]) { /* We have a match! */ + if (ninputs==0) {/* diff=F in network statistic specification */ + CHANGE_STAT[0] += change; + } else { /* diff=T */ + for (j=0; j head */ + b1 = tail; + echange = IS_OUTEDGE(b1,head) ? -1 : 1; + b1deg = OUT_DEG[b1]; + CHANGE_STAT[0] += (b1deg + echange > 1) - (b1deg > 1); +} + +/***************** + changestat: d_b1concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_b1concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. */ + int j, echange, b1attr; + Vertex b1, b1deg; + + /* *** don't forget tail -> head */ + b1 = tail; + echange = IS_OUTEDGE(b1,head) ? -1 : 1; + b1deg = OUT_DEG[b1]; + b1attr = INPUT_PARAM[N_CHANGE_STATS + b1 - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (b1attr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (b1deg + echange > 1) - (b1deg > 1); + } + } +} + +/***************** + changestat: d_b1nodematch +*****************/ +C_CHANGESTAT_FN(c_b1nodematch) { + + Vertex node3, node4, ninputs; + int count, exponenttype, matchval, b2attrsize, attrval1, attrval2, diffstatus; + /* int j, numofstats; */ + Edge e, e2; + double beta, alpha, change=0.0, exponent; + const int BetaType=1, AlphaType=2; + + b2attrsize = INPUT_PARAM[0]; + + if(b2attrsize > 0){ + ninputs = N_INPUT_PARAMS - N_NODES - b2attrsize;/*have 2 sets of node attributes and b2attrvals */ + } + else{ + ninputs = N_INPUT_PARAMS - BIPARTITE; + } + + diffstatus = !(ninputs == 3); /* 1 if Diff = T and 0 if Diff = F */ + /* numofstats = diffstatus ? (b2attrsize == 0 ? (ninputs - 3): (ninputs - 3) * b2attrsize) : (b2attrsize == 0 ? 1 : b2attrsize); */ + + exponent = beta = INPUT_PARAM[1]; /* exponent on nodematch count */ + exponenttype = BetaType; + alpha = INPUT_PARAM[2]; + + if (beta >= 1.0 && alpha < 1.0) { + exponent = alpha; + exponenttype = AlphaType; + } + // Rprintf("N_INPUT_PARAMS = %d, N_NODES=%d\n", N_INPUT_PARAMS, N_NODES); + // Rprintf("ninputs = %d, beta=%f, alpha=%f, exponenttype=%d, exponent=%f\n", + // ninputs, beta, alpha, exponenttype, exponent); + + matchval = INPUT_PARAM[tail + ninputs - 1]; + + /* Now count the neighbors of head whose attribute value equals matchval */ + /* All neighbors of head are inedges because this is a bipartite network */ + count = 0; + change = 0.0; + + if(b2attrsize == 0){ + + STEP_THROUGH_INEDGES(head, e, node3) { + if (INPUT_PARAM[node3 + ninputs - 1] == matchval && tail != node3) { /* match! */ + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + count = 0; + + STEP_THROUGH_OUTEDGES(tail, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != head) { /* RPB */ + count += IS_OUTEDGE(node3, node4); /* add 1 if node4 connects node3 with tail */ + } + } + + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + change += (count==0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + + if (exponenttype == BetaType && count>0) { + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, exponent); + change -= 0.5*count*(exponent==0.0? (count==1? 0.0 : 1.0) : pow((count-1), exponent)); + } + + if (diffstatus) { /* diff=T */ + if(ninputs==4) /* keep != NULL*/ + CHANGE_STAT[0] += edgestate ? -change : change; + else + CHANGE_STAT[matchval-1] += edgestate ? -change : change; + + } else { /* diff=F */ + CHANGE_STAT[0] += edgestate ? -change : change; + } + + } else { + + attrval1 = INPUT_PARAM[head + ninputs + b2attrsize - 1]; + + STEP_THROUGH_INEDGES(head, e, node3) { + + if (INPUT_PARAM[node3 + ninputs - 1] == matchval && tail != node3) { /* match! */ + + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + + count = 0; + + STEP_THROUGH_OUTEDGES(tail, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != head) { + attrval2 = INPUT_PARAM[node4 + ninputs + b2attrsize - 1]; + if(attrval2 == attrval1) count += IS_OUTEDGE(node3, node4); + } + } + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + /* setting the change stat for each parameter */ + change += (count== 0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + if (exponenttype == BetaType && count > 0) { + + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + if(diffstatus){ + CHANGE_STAT[b2attrsize*(matchval-1) + attrval1 - 1] += edgestate ? -change : change; + } else{ + CHANGE_STAT[attrval1 - 1] += edgestate ? -change : change; + } + + } +} + +/***************** + changestat: d_b1starmix +*****************/ +C_CHANGESTAT_FN(c_b1starmix) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, taild; + int nstats; + double tailattr, headattr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + taild = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + if(headattr == INPUT_ATTRIB[node3-1]){++taild;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr && + INPUT_ATTRIB[nnodes+nstats+j] == headattr) { + change = CHOOSE(taild, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b1starmixhomophily +*****************/ +C_CHANGESTAT_FN(c_b1starmixhomophily) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, taild; + double tailattr, headattr; + + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + taild = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + if(headattr == INPUT_ATTRIB[node3-1]){++taild;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr) { + change = CHOOSE(taild, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b1twostar +*****************/ +C_CHANGESTAT_FN(c_b1twostar) { + double change; + int j; + Edge e; + Vertex node3, nnodes; + int nstats; + double tailattr, headattr, n3attr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + + /* *** don't forget tail -> head */ + change = IS_OUTEDGE(tail = tail, head)? -1.0 : 1.0 ; + tailattr = INPUT_PARAM[tail-1]; + headattr = INPUT_PARAM[head-1]; + + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of tail */ + n3attr = INPUT_PARAM[node3-1]; + for(j=0; j < N_CHANGE_STATS; j++) { + if (node3 != head && INPUT_PARAM[nnodes + j] == tailattr && + INPUT_PARAM[nnodes + nstats + j] == MIN(headattr, n3attr) && + INPUT_PARAM[nnodes + 2*nstats + j] == MAX(headattr, n3attr)) { + CHANGE_STAT[j] += change; + } + } + } +} + +/***************** + changestat: d_b2concurrent +*****************/ +C_CHANGESTAT_FN(c_b2concurrent) { + int echange; + Vertex b2, b2deg; + + /* *** don't forget tail -> head */ + b2 = head; + echange = IS_OUTEDGE(tail, b2) ? -1 : 1; + b2deg = IN_DEG[b2]; + CHANGE_STAT[0] += (b2deg + echange > 1) - (b2deg > 1); +} + +/***************** + changestat: d_b2concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_b2concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes.*/ + int j, echange, b2attr; + Vertex b2, b2deg; + + + /* *** don't forget tail -> head */ + b2 = head; + echange = IS_OUTEDGE(tail, b2) ? -1 : 1; + b2deg = IN_DEG[b2]; + b2attr = INPUT_PARAM[N_CHANGE_STATS + b2 - 1 - BIPARTITE]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (b2attr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (b2deg + echange > 1) - (b2deg > 1); + } + } +} + +/***************** + changestat: d_b2nodematch +*****************/ +C_CHANGESTAT_FN(c_b2nodematch) { + + Vertex node3, node4, ninputs; + int count, exponenttype, matchval, b1attrsize, attrval1, attrval2, diffstatus; + /* int j, ind, numofstats; */ + Edge e, e2; + double beta, alpha, change=0.0, exponent; + const int BetaType=1, AlphaType=2; + + b1attrsize = INPUT_PARAM[0]; + + if(b1attrsize > 0){ + ninputs = N_INPUT_PARAMS - N_NODES - b1attrsize;/*have 2 sets of node attributes and b2attrvals */ + } + else{ + ninputs = N_INPUT_PARAMS - N_NODES + BIPARTITE; + } + + diffstatus = !(ninputs == 3); /* 1 if Diff = T and o if Diff = F - RPB */ + /* numofstats = diffstatus ? (b1attrsize == 0 ? (ninputs - 3): (ninputs - 3) * b1attrsize) : (b1attrsize == 0 ? 1 : b1attrsize); */ + + exponent = beta = INPUT_PARAM[1]; /* exponent on nodematch count */ + exponenttype = BetaType; + alpha = INPUT_PARAM[2]; + if (beta >= 1.0 && alpha < 1.0) { + exponent = alpha; + exponenttype = AlphaType; + } + // Rprintf("N_INPUT_PARAMS = %d, N_NODES=%d\n", N_INPUT_PARAMS, N_NODES); + // Rprintf("ninputs = %d, beta=%f, alpha=%f, exponenttype=%d, exponent=%f\n", + // ninputs, beta, alpha, exponenttype, exponent); + + matchval = INPUT_PARAM[head + ninputs - BIPARTITE - 1]; + /* Now count the neighbors of tail whose attribute value equals matchval */ + /* All neighbors of tail are outedges because this is a bipartite network */ + count=0; + change = 0.0; + + /* RPB */ + /* double CHANGE[b1attrsize]; */ + + + if(b1attrsize == 0){ + + STEP_THROUGH_OUTEDGES(tail, e, node3) { + if (INPUT_PARAM[node3 + ninputs - BIPARTITE - 1] == matchval && head != node3) { /* match! */ + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting head and node3 */ + count = 0; + + STEP_THROUGH_INEDGES(head, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != tail) { + count += IS_OUTEDGE(node4, node3); /* add 1 if node4 connects node3 with head */ + } + } + /* if count==0, then the statistic is always 1 */ + // Rprintf("count is %d\n", count); + change += (count==0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + + if (exponenttype == BetaType && count>0) { + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + + if (diffstatus) { /* diff=T */ + if(ninputs==4) /* keep != NULL*/ + CHANGE_STAT[0] += edgestate ? -change : change; + else + CHANGE_STAT[matchval-1] += edgestate ? -change : change; + + } else { /* diff=F */ + CHANGE_STAT[0] += edgestate ? -change : change; + } + } else { + + attrval1 = INPUT_PARAM[tail + ninputs + N_NODES + b1attrsize - BIPARTITE - 1]; + + STEP_THROUGH_OUTEDGES(tail, e, node3) { + + if (INPUT_PARAM[node3 + ninputs - BIPARTITE - 1] == matchval && head != node3) { /* match! */ + + ++count; + + // Rprintf("Matching twostar found! %d and %d connect to %d\n==================\n", tail, node3, head); + if (exponenttype == AlphaType) { + /* calculate alpha change stat instead of beta change stat. */ + /* Look for number of two-paths connecting tail and node3, not via head */ + + count = 0; + + STEP_THROUGH_INEDGES(head, e2, node4) { + // Rprintf("node3=%d, node4=%d, alpha=%f\n", node3,node4,alpha); + if (node4 != tail) { + attrval2 = INPUT_PARAM[node4 + ninputs + N_NODES + b1attrsize - BIPARTITE - 1]; + if(attrval2 == attrval1) count += IS_OUTEDGE(node4, node3); + } + } + /* if count==0, then the statistic is always (plus or minus) 1 */ + // Rprintf("count is %d\n", count); + /* setting the change stat for each parameter */ + change += (count== 0 ? 1 : pow(count+1, exponent) - pow(count, exponent)); + } + } + } + /* If count==0 then the statistic cannot change; it is the same with or */ + /* without the proposed toggle */ + if (exponenttype == BetaType && count > 0) { + + /* Now raise count and count+1 to beta, find the difference */ + change = 0.5*(count+1)*pow(count, beta); + change -= 0.5*count*(beta==0.0? (count==1? 0.0 : 1.0) : pow((count-1), beta)); + } + + if(diffstatus){ + CHANGE_STAT[b1attrsize*(matchval-1) + attrval1 - 1] += edgestate ? -change : change; + } else{ + CHANGE_STAT[attrval1 - 1] += edgestate ? -change : change; + } + + } +} + +/***************** + changestat: d_b2starmix +*****************/ +C_CHANGESTAT_FN(c_b2starmix) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, headd; + int nstats; + double tailattr, headattr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + headd = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == tailattr && + INPUT_ATTRIB[nnodes+nstats+j] == headattr) { + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b2starmixhomophily +*****************/ +C_CHANGESTAT_FN(c_b2starmixhomophily) { + double change; + int j, kmo; + Edge e; + Vertex node3, nnodes, headd; + double tailattr, headattr; + + nnodes = N_NODES; + kmo = (int)INPUT_PARAM[0] - 1; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + headattr = INPUT_ATTRIB[head-1]; + headd = -(int)edgestate; /* if edge exists set to -1 because it will be recounted */ + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + if (INPUT_ATTRIB[nnodes+j] == headattr) { + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_b2twostar +*****************/ +C_CHANGESTAT_FN(c_b2twostar) { + double change; + int j; + Edge e; + Vertex node3, nnodes; + int nstats; + double tailattr, headattr, n3attr; + + nstats = (int)N_CHANGE_STATS; + nnodes = N_NODES; + + + /* *** don't forget tail -> head */ + change = edgestate? -1.0 : 1.0 ; + tailattr = INPUT_PARAM[tail-1]; + headattr = INPUT_PARAM[head-1]; + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + n3attr = INPUT_PARAM[node3-1]; + for(j=0; j < N_CHANGE_STATS; j++) { + if (node3 != tail && INPUT_PARAM[nnodes + j] == headattr && + INPUT_PARAM[nnodes + nstats + j] == MIN(tailattr, n3attr) && + INPUT_PARAM[nnodes + 2*nstats + j] == MAX(tailattr, n3attr)) { + CHANGE_STAT[j] += change; + } + } + } +} + +/***************** + changestat: d_balance +*****************/ +C_CHANGESTAT_FN(c_balance) { + int a, b, c, d, e, edgecount, t300, + t210, t120C, t120U, t120D, t201, t030C, t030T, t111U, + t111D, t021C, t021U, t021D, t102, t012; /* , t003; */ + Vertex node3; + + + /* *** don't forget tail -> head */ + if (DIRECTED) { /* directed version */ + t300 = 0; + t210 = 0; + t120C = 0; t120U = 0; t120D = 0; t201 = 0; + t030C = 0; t030T = 0; t111U = 0; t111D = 0; + t021C = 0; t021U = 0; t021D = 0; t102 = 0; + t012 = 0; + + if (MIN_OUTEDGE(head)!=0 || MIN_INEDGE(head)!=0 || + MIN_OUTEDGE(tail)!=0 || MIN_INEDGE(tail)!=0) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = IS_OUTEDGE(head, tail); + b = IS_OUTEDGE(head, node3); + c = IS_OUTEDGE(node3, head); + d = IS_OUTEDGE(node3, tail); + e = IS_OUTEDGE(tail, node3); + edgecount = (a + b + c + d + e); + + switch(edgecount) { + case 0: /* 012 */ + ++t012; + + case 1: { /* 021C, 021U, 021D, 102 */ + if ((b == 1) || (d == 1)) + ++t021C; + if (c == 1) + ++t021U; + if (e == 1) + ++t021D; + if (a == 1) + ++t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if ((b + d) == 2) + ++t030C; + if (((b + e) == 2) || ((c + d) == 2) || ((c + e) == 2)) + ++t030T; + if (((a + b) == 2) || ((a + e) == 2) || ((d + e) == 2)) + ++t111U; + if (((a + c) == 2) || ((a + d) == 2) || ((b + c) == 2)) + ++t111D; + } + break; + + case 3: { /* 120C, 120U, 120D, 201 */ + if (a == 1) { + if (((b + d) == 2) || ((c + e) == 2)) + ++t120C; + if ((b + e) == 2) + ++t120U; + if ((c + d) == 2) + ++t120D; + if (((b + c) == 2) || ((d + e) == 2)) + ++t201; + } else { + if (b == 1) { + if (((c + d) == 2) || ((d + e) == 2)) + ++t120C; + if ((c + e) == 2) + ++t120D; + } else { + ++t120U; + } + } + } + break; + + case 4: /* 210 */ + ++t210; + break; + + case 5: /* 300 */ + ++t300; + break; + } + + switch(edgecount) { + case 1: /* 102, 021D, 021U, 021C */ + --t012; + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if (((a + c) == 2) || ((a + e) == 2) || ((b + d) == 2) || + ((c + e) == 2)) + --t021C; + if (((a + d) == 2) || ((b + e) == 2)) + --t021U; + if (((a + b) == 2) || ((c + d) == 2)) + --t021D; + if (((b + c) == 2) || ((d + e) == 2)) + --t102; + } + break; + + case 3: { /* 201, 120D, 120U, 120C */ + if (a == 1) { + if ((c + e) == 2) + --t030C; + if (((c + d) == 2) || ((b + e) == 2) || ((b + d) == 2)) + --t030T; + if ((b + c) == 2) + --t111U; + if ((d + e) == 2) + --t111D; + } else { + if (b == 1) { + if ((c + d) == 2) + --t111U; + if (((c + e) == 2) || ((d + e) == 2)) + --t111D; + } + else + --t111U; + } + } + break; + + case 4: { /* 210 */ + if (a == 1) { + if (((b + c + e) == 3) || ((c + d + e) == 3)) + --t120C; + if ((b + c + d) == 3) + --t120U; + if ((b + d + e) == 3) + --t120D; + } else { + if ((b + c + d + e) == 4) + --t201; + } + } + break; + + case 5: /* 300 */ + --t210; + break; + } + } + } /* ****** move to next node3 ******** */ + } + else + t012 = t012 + (N_NODES - 2); + + /* t003 = (t300+t210+t120C+t120U+t120D+t201+t030C+t030T); + t003 = t003+(t111U+t111D+t021C+t021U+t021D+t102+t012); */ + b = t102 + t300; + CHANGE_STAT[0] += edgestate ? -(double)b : (double)b; + + /* *** don't forget tail -> head */ + }else{ /* undirected */ + t300 = 0; t201 = 0; t102 = 0; t012 = 0; + + if (MIN_OUTEDGE(head)!=0 || MIN_INEDGE(head)!=0 || + MIN_OUTEDGE(tail)!=0 || MIN_INEDGE(tail)!=0) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = IS_UNDIRECTED_EDGE(node3, head); + b = IS_UNDIRECTED_EDGE(node3, tail); + edgecount = (a + b); + + switch(edgecount){ + case 0: { /* 012 */ + ++t102; + --t012; + } + break; + + case 1: { /* 021C, 021U, 021D, 102 */ + ++t201; + --t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + ++t300; + --t201; + } + break; + } + } + + } /* ****** move to next node3 ******** */ + } else + t102 = t102 + (N_NODES - 2); + + /* t003 = (t102+t201+t300); */ + b = t102 + t300; + CHANGE_STAT[0] += edgestate ? -(double)b : (double)b; + } +} + +/***************** + changestat: d_boundeddegree +*****************/ +C_CHANGESTAT_FN(c_boundeddegree) { + int j, echange; + Vertex taild, headd=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = DEG(tail); + headd = DEG(head); + for(j = 0; j+1 < nstats; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + CHANGE_STAT[j] += (headd + echange == deg) - (headd == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); + CHANGE_STAT[nstats-1] += (headd + echange >= bound) - (headd >= bound); +} + +/***************** + changestat: d_boundedidegree +*****************/ +C_CHANGESTAT_FN(c_boundedidegree) { + int j, echange; + Vertex taild=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = IN_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); +} + +/***************** + changestat: d_boundedistar +*****************/ +C_CHANGESTAT_FN(c_boundedistar) { + double change, headod; + double newheadod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + headod = IN_DEG[head]; + newheadod = headod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = MIN(bound,CHOOSE(newheadod, k))-MIN(bound,CHOOSE(headod, k)); + CHANGE_STAT[j] += change; + } +} + +/***************** + changestat: d_boundedkstar +*****************/ +C_CHANGESTAT_FN(c_boundedkstar) { + double change, tailod, headod; + double newtailod, newheadod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + tailod = DEG(tail); + newtailod = tailod + (edgestate ? -1 : 1); + headod = DEG(head); + newheadod = headod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = (MIN(bound,CHOOSE(newtailod, k))-MIN(bound,CHOOSE(tailod, k))) + + (MIN(bound,CHOOSE(newheadod, k))-MIN(bound,CHOOSE(headod, k))); + + CHANGE_STAT[j] += change; /* (edgestate ? - change : change); */ + } +} + +/***************** + changestat: d_boundedodegree +*****************/ +C_CHANGESTAT_FN(c_boundedodegree) { + int j, echange; + Vertex taild=0, deg; + int nstats = (int)N_CHANGE_STATS; + Vertex bound = (Vertex)INPUT_PARAM[nstats-1]; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = OUT_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taild + echange == deg) - (taild == deg); + } + CHANGE_STAT[nstats-1] += (taild + echange >= bound) - (taild >= bound); +} + +/***************** + changestat: d_boundedostar +*****************/ +C_CHANGESTAT_FN(c_boundedostar) { + double change, tailod; + double newtailod; + int j, k, bound; + int p = N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + /* is there an edge for this toggle */ + tailod = OUT_DEG[tail]; + newtailod = tailod + (edgestate ? -1 : 1); + for(j=0; j < p; j++) { + k = ((int)INPUT_PARAM[j]); + bound = (int)INPUT_PARAM[j+p]; + change = MIN(bound,CHOOSE(newtailod, k))-MIN(bound,CHOOSE(tailod, k)); + CHANGE_STAT[j] += change; + } + } + +/***************** + changestat: d_boundedtriangle +*****************/ +Vertex CountTriangles (Vertex tail, Vertex head, int outcount, + int incount, Network *nwp); +C_CHANGESTAT_FN(c_boundedtriangle) { + Edge e; + Vertex node3; + double boundedchange, htcount; + Vertex tailtri, headtri; + int bound = (int)INPUT_PARAM[0]; + + /* *** don't forget tail -> head */ + tailtri=0; + headtri=0; + STEP_THROUGH_OUTEDGES(tail, e, node3) { + tailtri += CountTriangles(tail, node3, 1, 1, nwp); + } + STEP_THROUGH_INEDGES(tail, e, node3) { + tailtri += CountTriangles(tail, node3, 1, 1, nwp); + } + STEP_THROUGH_OUTEDGES(head, e, node3) { + headtri += CountTriangles(head, node3, 1, 1, nwp); + } + STEP_THROUGH_INEDGES(head, e, node3) { + headtri += CountTriangles(head, node3, 1, 1, nwp); + } + tailtri = tailtri/2; + headtri = headtri/2; + htcount = CountTriangles(tail, head, 1, 1, nwp); + boundedchange = (MIN(headtri+(edgestate ? -1:1)*htcount,bound)-MIN(headtri,bound)+ + MIN(tailtri+(edgestate ? -1:1)*htcount,bound)-MIN(tailtri,bound)); + CHANGE_STAT[0] += boundedchange; +} + +/***************** + CountTriangles: called by d_boundedtriangle +*****************/ +Vertex CountTriangles (Vertex tail, Vertex head, int outcount, int incount, + Network *nwp) { + Edge e; + Vertex change; + Vertex k; + + /* *** don't forget tail -> head */ + change=0; + if(outcount){ + STEP_THROUGH_OUTEDGES(head, e, k) /* step through outedges of head */ + { + if (IS_UNDIRECTED_EDGE(k,tail)) + ++change; + } + } + + if(incount){ + STEP_THROUGH_INEDGES(head, e, k) /* step through inedges of head */ + { + if (IS_UNDIRECTED_EDGE(k,tail)) + ++change; + } + } + return(change); +} + + + +/******************** changestats: C ***********/ +/***************** + changestat: d_concurrent +*****************/ +C_CHANGESTAT_FN(c_concurrent) { + int echange; + Vertex taildeg, headdeg; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + if(!DIRECTED){ + taildeg += IN_DEG[tail]; + headdeg += OUT_DEG[head]; + } + CHANGE_STAT[0] += (taildeg + echange > 1) - (taildeg > 1); + CHANGE_STAT[0] += (headdeg + echange > 1) - (headdeg > 1); +} + +/***************** + changestat: d_concurrent_by_attr +*****************/ +C_CHANGESTAT_FN(c_concurrent_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex taildeg, headdeg; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + if(!DIRECTED){ + taildeg += IN_DEG[tail]; + headdeg += OUT_DEG[head]; + } + tailattr = INPUT_PARAM[N_CHANGE_STATS + tail - 1]; + headattr = INPUT_PARAM[N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + if (tailattr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (taildeg + echange > 1) - (taildeg > 1); + } + if (headattr == INPUT_PARAM[j]) { /* we have attr match */ + CHANGE_STAT[j] += (headdeg + echange > 1) - (headdeg > 1); + } + } +} + +/***************** + changestat: d_ctriple +*****************/ +C_CHANGESTAT_FN(c_ctriple) { + Edge e; + Vertex change, node3; + int j; + double tailattr, edgemult; + + /* *** don't forget tail -> head */ + edgemult = edgestate ? -1.0 : 1.0; + change = 0; + if(N_INPUT_PARAMS > 0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]) { + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_OUTEDGE(node3, tail); + } + if(N_CHANGE_STATS > 1) { /* diff = TRUE; matches must be tabled */ + for (j=0; j=(from) && (x)<(to)) + +/***************** + changestat: d_degrange +*****************/ +C_CHANGESTAT_FN(c_degrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex taildeg = DEG(tail), headdeg = DEG(head); + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } +} + +/***************** + changestat: d_degrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_degrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 3*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex taildeg = DEG(tail), headdeg = DEG(head); + int tailattr = INPUT_PARAM[3*N_CHANGE_STATS + tail - 1], + headattr = INPUT_PARAM[3*N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } +} + +/***************** + changestat: d_degrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_degrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of degrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + Vertex taildeg, headdeg, v; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + taildeg=headdeg=-1; /* since tailattr==headattr, subtract the automatic match */ + taildeg=headdeg=0; + STEP_THROUGH_OUTEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_INEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_OUTEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + STEP_THROUGH_INEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(taildeg + echange, from, to) - FROM_TO(taildeg, from, to); + CHANGE_STAT[j] += FROM_TO(headdeg + echange, from, to) - FROM_TO(headdeg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_degree +*****************/ +C_CHANGESTAT_FN(c_degree) { + int j, echange; + Vertex taildeg, headdeg, deg; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + taildeg = DEG(tail); + headdeg = DEG(head); + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } +} + +/***************** + changestat: c_degdist +*****************/ +C_CHANGESTAT_FN(c_degdist) { + int echange = edgestate ? -1:+1; + + Vertex otd = DEG(tail), ohd = DEG(head), + ntd = otd + echange, nhd = ohd + echange; + + if(ntd > N_CHANGE_STATS || nhd > N_CHANGE_STATS) cutoff_error(mtp); + + if(otd) CHANGE_STAT[otd-1]--; + if(ohd) CHANGE_STAT[ohd-1]--; + if(ntd) CHANGE_STAT[ntd-1]++; + if(nhd) CHANGE_STAT[nhd-1]++; +} + +/***************** + changestat: d_degreepopularity +*****************/ +C_CHANGESTAT_FN(c_degreepopularity) { + double change; + + /* *** don't forget tail -> head */ + change = 0.0; + Vertex tdeg = DEG(tail); + Vertex hdeg = DEG(head); + if(edgestate){ + change -= sqrt(tdeg); + change += (tdeg-1.0)*(sqrt(tdeg-1.0)-sqrt(tdeg)); + change -= sqrt(hdeg); + change += (hdeg-1.0)*(sqrt(hdeg-1.0)-sqrt(hdeg)); + }else{ + change += sqrt(tdeg+1.0); + change += tdeg*(sqrt(tdeg+1.0)-sqrt(tdeg)); + change += sqrt(hdeg+1.0); + change += hdeg*(sqrt(hdeg+1.0)-sqrt(hdeg)); + } + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_degree_by_attr +*****************/ +C_CHANGESTAT_FN(c_degree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr, testattr; + Vertex taildeg, headdeg, d; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:1; + taildeg = DEG(tail); + headdeg = DEG(head); + tailattr = INPUT_PARAM[2*N_CHANGE_STATS + tail - 1]; + headattr = INPUT_PARAM[2*N_CHANGE_STATS + head - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += (taildeg + echange == d) - (taildeg == d); + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += (headdeg + echange == d) - (headdeg == d); + } +} + +/***************** + changestat: d_degree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_degree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex taildeg, headdeg, deg, v; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + tailattr = (int)nodeattr[tail]; + headattr = (int)nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + echange = edgestate ? -1:1; + taildeg=headdeg=-1; /* since tailattr==headattr, subtract the automatic match */ + taildeg=headdeg=0; + STEP_THROUGH_OUTEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_INEDGES(tail, e, v) { taildeg += (nodeattr[v]==tailattr); } + STEP_THROUGH_OUTEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + STEP_THROUGH_INEDGES(head, e, v) { headdeg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } + } +} + + +/***************** + changestat: c_distance +*****************/ +C_CHANGESTAT_FN(c_distance) { + Vertex t,h; + int j; + int nv,dim,logd,sphd; + double dis,dexp,*coord,sphr,logmind,logdoff,scale,mdis,pw; + + /*Set things up*/ + nv = (int)INPUT_PARAM[0]; /*Number of vertices*/ + dim = (int)INPUT_PARAM[1]; /*Dimension of coordinate space*/ + dexp = INPUT_PARAM[2]; /*Minkowski exponent*/ + logd = (int)INPUT_PARAM[3]; /*Should we use the log distance?*/ + sphd = (int)INPUT_PARAM[4]; /*Should we use lat/lon distances on the geosphere?*/ + sphr = INPUT_PARAM[5]; /*Radius to use for spherical coordinates*/ + logmind = INPUT_PARAM[6]; /*Distance minimum for the log case*/ + logdoff = INPUT_PARAM[7]; /*Distance offset for the log case*/ + scale = INPUT_PARAM[8]; /*Scaling adjustment for raw distance*/ + pw = INPUT_PARAM[9]; /*Power to which distance should be raised*/ + coord = INPUT_PARAM+10; /*Pointer to the coordinate matrix*/ + + /*Compute the changescore*/ + t=tail-1; /*0-indexed values, for C convenience*/ + h=head-1; + if(sphd){ /*Compute great circle distances on the sphere*/ + dis=spheredist(coord[t],coord[t+nv],coord[h],coord[h+nv],sphr); + }else{ + dis=0.0; + } + /*Compute Minkowski distances in free space (adding extra dimensions to + great circle distance if sphd and we were given dim>2)*/ + mdis=0.0; + for(j=0+2*sphd;j 0){ + nrow = (N_NODES)-(long int)(INPUT_PARAM[0]); + }else{ + nrow = (long int)(INPUT_PARAM[0]); + } + +/* Rprintf("nrow %d noffset %d\n",nrow, noffset); + Rprintf("attrib: "); + for(i=0;i<1000;i++) + Rprintf("%1.0f",INPUT_ATTRIB[i]); + + Rprintf("\n;"); */ + + if(DIRECTED){ + /* directed version */ + /*Get the initial state of the edge and its reflection*/ + refedgestate = (IS_OUTEDGE(head, tail)); + + /* *** don't forget tail -> head */ + + /*Get the dyadic covariate*/ + /* val = INPUT_ATTRIB[(head-1-nrow)+(tail-1)*ncols]; */ + index = (head-1-noffset)*nrow+(tail-1); + if(index >= 0 && index <= nrow*nrow){ + val = INPUT_ATTRIB[(head-1-noffset)*nrow+(tail-1)]; + /* Rprintf("tail %d head %d nrow %d ncols %d val %f\n",tail, head, nrow, ncols, val); */ + + /*Update the change statistics, as appropriate*/ + if(refedgestate){ /* Reflected edge is present */ + if(edgestate){ /* Toggled edge _was_ present */ + if(head>tail){ /* Mut to low->high */ + CHANGE_STAT[0] -= val; + CHANGE_STAT[1] += val; + }else{ /* Mut to high->low */ + CHANGE_STAT[0] -= val; + CHANGE_STAT[2] += val; + } + }else{ /* Toggled edge _was not_ present */ + if(head>tail){ /* Low->high to mut */ + CHANGE_STAT[1] -= val; + CHANGE_STAT[0] += val; + }else{ /* High->low to mut */ + CHANGE_STAT[2] -= val; + CHANGE_STAT[0] += val; + } + } + }else{ /* Reflected edge is absent */ + if(edgestate){ /* Toggled edge _was_ present */ + if(head>tail){ /* High->low to null */ + CHANGE_STAT[2] -= val; + }else{ /* Low->high to null */ + CHANGE_STAT[1] -= val; + } + }else{ /* Toggled edge _was not_ present */ + if(head>tail){ /* Null to high->low */ + CHANGE_STAT[2] += val; + }else{ /* Null to low->high */ + CHANGE_STAT[1] += val; + } + } + } + } +}else{ + /* undirected case (including bipartite) */ + + /* *** don't forget tail -> head */ + /*Get the initial edge state*/ + /*Get the covariate value*/ + /* val = INPUT_ATTRIB[(head-1-nrow)+(tail-1)*ncols]; */ + index = (head-1-noffset)*nrow+(tail-1); + if(index >= 0 && index <= nrow*((long int)(INPUT_PARAM[0]))){ + val = INPUT_ATTRIB[(head-1-noffset)*nrow+(tail-1)]; + /*Update the change statistic, based on the toggle type*/ + /* Rprintf("tail %d head %d nrow %d noffset %d val %f\n",tail, head, nrow, noffset, val); */ + /*Update the change statistic, based on the toggle type*/ + CHANGE_STAT[0] += edgestate ? -val : val; + } + } +} + + +/******************** changestats: G ***********/ + +/***************** + changestat: d_gwdegree +*****************/ + +#define GWD0(d0) ((decay) ? exp(loneexpd*(d0)) : (d0)==0) + +C_CHANGESTAT_FN(c_gwdegree) { + int echange=0; + double decay, loneexpd, change; + Vertex taild, headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + change = 0.0; + echange = edgestate ? -1:+1; + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + change += echange*(GWD0(taild) + GWD0(headd)); + + CHANGE_STAT[0] = change; + +} + +/***************** + changestat: d_gwdegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwdegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 200?) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int tailattr, headattr, echange=0; + double decay, loneexpd; + Vertex taild, headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:+1; + taild = DEG(tail) - edgestate; + tailattr = INPUT_PARAM[tail]; + CHANGE_STAT[tailattr-1] += echange*GWD0(taild); + + headd = DEG(head) - edgestate; + headattr = INPUT_PARAM[head]; + CHANGE_STAT[headattr-1] += echange*GWD0(headd); + +} + +/***************** + changestat: d_gwidegree +*****************/ +C_CHANGESTAT_FN(c_gwidegree) { + double decay, loneexpd, change; + Vertex headd=0; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + change = 0.0; + + /* *** don't forget tail -> head */ + headd = IN_DEG[head] - edgestate; + change += (edgestate? -1.0 : 1.0) * GWD0(headd); + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_gwidegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwidegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 2008) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int headattr, echange; + double decay, loneexpd; + Vertex headd; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + headd = IN_DEG[head] - edgestate; + headattr = INPUT_PARAM[head - BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + CHANGE_STAT[headattr-1] += echange*GWD0(headd); +} + +/***************** + changestat: d_gwodegree +*****************/ +C_CHANGESTAT_FN(c_gwodegree) { + double decay, loneexpd, change; + Vertex taild; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + change = 0.0; + + /* *** don't forget tail -> head */ + taild = OUT_DEG[tail] - edgestate; + change += (edgestate? -1 : 1) * GWD0(taild); + CHANGE_STAT[0] = change; +} + +/***************** + changestat: d_gwodegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_gwodegree_by_attr) { + /*The inputparams are assumed to be set up as follows: + The first value is the decay parameter (as in Hunter et al, JASA 2008) + The next sequence of values is the nodal attributes, coded as integers + from 1 through N_CHANGE_STATS + */ + int tailattr, echange; + double decay, loneexpd; + Vertex taild; + + decay = INPUT_PARAM[0]; + loneexpd = log1mexp(decay); + + /* *** don't forget tail -> head */ + echange = edgestate ? -1 : 1; + taild = OUT_DEG[tail] - edgestate; + tailattr = INPUT_PARAM[tail]; + CHANGE_STAT[tailattr-1] += echange*GWD0(taild); +} + +/******************** changestats: H ***********/ +/***************** + changestat: d_hamming +*****************/ +/* This function must be passed two networks in the forms of edgelists: + One is the network from which hamming distances are calculated and the + other is the network of weights on the dyads, which is an edgelist with a + third column of weights. Note that all non-edges in the second network + will have the value given by defaultval; thus, the "unweighted" hamming + distance is obtained when the default is 1.0 and the second network is + empty. */ +C_CHANGESTAT_FN(c_hamming) { + int discord; + + Edge wt_net_start= INPUT_PARAM[0]*2+2; + double defaultval = INPUT_PARAM[wt_net_start-1]; /* Hamming wt for non-edges in cov nwp */ + double *wt_net = INPUT_PARAM+wt_net_start; + + /* *** don't forget tail -> head */ + + discord = XOR(dEdgeListSearch(tail, head, INPUT_PARAM), edgestate); + + /* Second, search second network to see if the weight is different from + defaultval. In unweighted case, this network is empty. */ + + Edge wt_pos = dEdgeListSearch(tail, head, wt_net); + double val = wt_pos ? wt_net[wt_pos+2*(unsigned int)wt_net[0]] : defaultval; + + CHANGE_STAT[0] += (discord ? -val : val); + +} + +/******************** changestats: I ***********/ + +// A macro indicating whether x is in [from,to) +#define FROM_TO(x, from, to) ((x)>=(from) && (x)<(to)) + +/***************** + changestat: d_idegrange +*****************/ +C_CHANGESTAT_FN(c_idegrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex headideg = IN_DEG[head]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } +} + +/***************** + changestat: d_idegrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_idegrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex headideg = IN_DEG[head]; + int headattr = INPUT_PARAM[3*N_CHANGE_STATS + head - 1 - BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + for(j = 0; j < N_CHANGE_STATS; j++){ + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (headattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } +} + +/***************** + changestat: d_idegrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_idegrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of idegrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (headattr == tailattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + Vertex headideg=0, v; + STEP_THROUGH_INEDGES(head, e, v) { headideg += (nodeattr[v]==headattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(headideg + echange, from, to) - FROM_TO(headideg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_idegree +*****************/ +C_CHANGESTAT_FN(c_idegree) { + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1 : +1; + Vertex headd = IN_DEG[head]; + + for(j=0; j < N_CHANGE_STATS; j++){ + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] += (headd + echange == deg) - (headd == deg); + } +} + + +/***************** + changestat: d_idegdist +*****************/ +C_CHANGESTAT_FN(c_idegdist) { + int echange = edgestate ? -1 : +1; + Vertex ohd = IN_DEG[head], + nhd = ohd + echange; + + if(nhd > N_CHANGE_STATS) cutoff_error(mtp); + + if(ohd) CHANGE_STAT[ohd-1]--; + if(nhd) CHANGE_STAT[nhd-1]++; +} + + +/***************** + changestat: d_idegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_idegree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, headattr, testattr; + Vertex headdeg, d; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1 : +1; + headdeg = IN_DEG[head]; + headattr = INPUT_PARAM[2*N_CHANGE_STATS + head - 1- BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (headattr == testattr) /* we have head attr match */ + CHANGE_STAT[j] += (headdeg + echange == d) - (headdeg == d); + } +} + +/***************** + changestat: d_idegree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_idegree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j, echange, tailattr, headattr; + Vertex headdeg, deg, tmp; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + tailattr = (int)nodeattr[tail]; + headattr = (int)nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + echange=edgestate ? -1 : +1; + headdeg=0; + STEP_THROUGH_INEDGES(head, e, tmp){ + headdeg += (nodeattr[tmp]==headattr); + } + for(j = 0; j < N_CHANGE_STATS; j++) { + deg = (Vertex)INPUT_PARAM[j]; + CHANGE_STAT[j] += (headdeg + echange == deg) - (headdeg == deg); + } + } +} + +/***************** + changestat: d_idegreepopularity +*****************/ +C_CHANGESTAT_FN(c_idegreepopularity) { + double change; + Vertex deg=0; + + /* *** don't forget tail -> head */ + change = 0.0; + deg = (double)(IN_DEG[head]); + if(edgestate){ + change -= sqrt(deg); + change += (deg-1.0)*(sqrt(deg-1.0)-sqrt(deg)); + }else{ + change += sqrt(deg+1.0); + change += deg*(sqrt(deg+1.0)-sqrt(deg)); + } + CHANGE_STAT[0]=change; +} + +/***************** + changestat: d_intransitive +*****************/ +C_CHANGESTAT_FN(c_intransitive) { + Edge e; + Vertex node2; + double change; + + /* *** don't forget tail -> head */ + change = 0.0; + STEP_THROUGH_OUTEDGES(head, e, node2) { + if (node2 != tail){ + if (!IS_OUTEDGE(tail,node2)){ + change = change + 1.0; + } + } + } + STEP_THROUGH_INEDGES(head, e, node2) { + if (node2 != tail){ + if (IS_OUTEDGE(tail, node2)){ + change = change - 1.0; + } + } + } + STEP_THROUGH_INEDGES(tail, e, node2) { + if (node2 != head){ + if (!IS_OUTEDGE(node2,head)){ + change = change + 1.0; + } + } + } + CHANGE_STAT[0] += edgestate ? -change : change; +/* Rprintf("tail %d head %d edgestate %d change %f\n",tail,head, change); */ +} + +/***************** +changestat: d_isolatededges +*****************/ +D_CHANGESTAT_FN(d_isolatededges) { + int i, edgestate; + Vertex tail, head, neighbor, taild, headd; + Edge e; + + /* *** don't forget tail -> head */ + ZERO_ALL_CHANGESTATS(i); + FOR_EACH_TOGGLE(i) { + // is there an edge tail -> head? + edgestate = IS_OUTEDGE(tail=TAIL(i), head=HEAD(i)); + + taild = DEG(tail); + headd = DEG(head); + + if(edgestate) { // we are removing an edge + + // if head and tail both have degree one, then + // we are removing an isolated edge + if(taild == 1 && headd == 1) + CHANGE_STAT[0] -= 1; + + // if tail has degree 2 and has a degree one node other than head as a neighbor, + // then we are making a non-isolated edge into an isolated edge by removing + // the edge tail -> head + if(taild == 2) { + STEP_THROUGH_OUTEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != head) + CHANGE_STAT[0] += 1; + } + STEP_THROUGH_INEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != head) + CHANGE_STAT[0] += 1; + } + } + + // ditto head + if(headd == 2) { + STEP_THROUGH_OUTEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != tail) + CHANGE_STAT[0] += 1; + } + STEP_THROUGH_INEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1 && neighbor != tail) + CHANGE_STAT[0] += 1; + } + } + } else { // we are adding an edge + + // if head and tail both have degree zero, then + // we are adding an isolated edge + if(taild == 0 && headd == 0) + CHANGE_STAT[0] += 1; + + // if tail has degree 1 and so does its neighbor, then we + // are making an isolated edge into a non-isolated edge; + // note that for undirected graphs, this neighbor cannot + // be head, as the current toggle is to turn on the edge + // tail -> head + if(taild == 1) { + STEP_THROUGH_OUTEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + STEP_THROUGH_INEDGES(tail, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + } + + // ditto head + if(headd == 1) { + STEP_THROUGH_OUTEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + STEP_THROUGH_INEDGES(head, e, neighbor) { + if(DEG(neighbor) == 1) + CHANGE_STAT[0] -= 1; + } + } + } + + TOGGLE_IF_MORE_TO_COME(i); + } + + UNDO_PREVIOUS_TOGGLES(i); +} + +/***************** + changestat: d_isolates +*****************/ +C_CHANGESTAT_FN(c_isolates) { + int echange; + Vertex taild, headd=0; + + /* *** don't forget tail -> head */ + echange = edgestate ? -1:+1; + taild = DEG(tail); + headd = DEG(head); + CHANGE_STAT[0] += (taild + echange == 0) - (taild == 0); + CHANGE_STAT[0] += (headd + echange == 0) - (headd == 0); + +} + +S_CHANGESTAT_FN(s_isolates) { + /* *** don't forget tail -> head */ + CHANGE_STAT[0] = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++){ + if(DEG(tail) == 0) + CHANGE_STAT[0] ++; + } +} + +/***************** + changestat: d_istar +*****************/ +C_CHANGESTAT_FN(c_istar) { + double change, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double tailattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + headd = -(int)edgestate; + STEP_THROUGH_INEDGES(head, e, node3) {/* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headd = IN_DEG[head] - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/******************** changestats: K ***********/ +/***************** + changestat: d_kstar +*****************/ +C_CHANGESTAT_FN(c_kstar) { + double change, taild, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double tailattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + taild = -(int)edgestate; + STEP_THROUGH_OUTEDGES(tail, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++taild;} + } + STEP_THROUGH_INEDGES(tail, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++taild;} + } + headd = -(int)edgestate; + STEP_THROUGH_OUTEDGES(head, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + STEP_THROUGH_INEDGES(head, e, node3) { + if(tailattr == INPUT_ATTRIB[node3-1]){++headd;} + } + + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; +/* if (kmo==0) { + change=1; + } else { */ + change = CHOOSE(taild, kmo) + CHOOSE(headd, kmo); +/* } uncomment these few lines to define 1-stars as equivalent to + edges (currently, each edge is counted as two 1-stars) */ + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* *** don't forget tail -> head */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + taild = DEG(tail) - edgestate; + headd = DEG(head) - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) + { + kmo = ((int)INPUT_PARAM[j]) - 1; +/* if (kmo==0) { + change=1; + } else { */ + change = CHOOSE(taild, kmo) + CHOOSE(headd, kmo); +/* } uncomment these few lines to define 1-stars as equivalent to + edges (currently, each edge is counted as two 1-stars) */ + CHANGE_STAT[j] += (edgestate ? - change : change); + } + + } +} + + +/******************** changestats: L ***********/ +/***************** + changestat: d_localtriangle +*****************/ +C_CHANGESTAT_FN(c_localtriangle) { + Edge e; + Vertex node3, nmat; + double change; + + nmat = (Vertex)(INPUT_PARAM[0]); + + /* *** don't forget tail -> head */ + change = 0.0; + + if(INPUT_PARAM[1+(head-1)+(tail-1)*nmat] == 1.0){ + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(INPUT_PARAM[1+(node3-1)+(tail-1)*nmat] == 1.0 && + INPUT_PARAM[1+(node3-1)+(head-1)*nmat] == 1.0 ){ + if (DIRECTED){ + if (IS_INEDGE(node3,tail) ) ++change; + if (IS_OUTEDGE(node3,tail)) ++change; + }else{ + if (IS_UNDIRECTED_EDGE(node3,tail)) ++change; + } + } + } + + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(INPUT_PARAM[1+(node3-1)+(tail-1)*nmat] == 1.0 && + INPUT_PARAM[1+(node3-1)+(head-1)*nmat] == 1.0 ){ + if (DIRECTED) + { + if (IS_INEDGE(node3,tail) ) ++change; + if (IS_OUTEDGE(node3,tail)) ++change; + } + else + { + if (IS_UNDIRECTED_EDGE(node3,tail)) ++change; + } + } + } + + CHANGE_STAT[0] += edgestate ? - change : change; + + } +} + +/******************** changestats: M ***********/ +/***************** + changestat: d_m2star +*****************/ +C_CHANGESTAT_FN(c_m2star) { + int tailid, headod, change; + int backedgestate; + + + /* *** don't forget tail -> head */ + /* edgestate is 1 if the edge from tail to head */ + /* exists and will disappear */ + /* edgestate is 0 if the edge does not exist */ + backedgestate = (IS_OUTEDGE(head, tail)); + + tailid = IN_DEG[tail]; + headod = OUT_DEG[head]; + change = tailid + headod - 2*backedgestate; + CHANGE_STAT[0] += (edgestate ? -change : change); + +} + +/***************** + changestat: d_mutual + + (1,1) -> anything = -1 + anything -> (1,1) = +1 +*****************/ +C_CHANGESTAT_FN(c_mutual) { + double matchval, change; + int j, ninputs, noattr; + + ninputs = N_INPUT_PARAMS - N_NODES; + noattr = (N_INPUT_PARAMS == 0); + + /* *** don't forget tail -> head */ + if (IS_OUTEDGE(head,tail)) { /* otherwise, no change occurs */ + change = edgestate ? -1.0 : 1.0 ; + if (noattr) { /* "plain vanilla" mutual, without node attributes */ + CHANGE_STAT[0] += change; + } else { /* Only consider mutuals where node attributes match */ + matchval = INPUT_PARAM[tail+ninputs-1]; + if (matchval == INPUT_PARAM[head+ninputs-1]) { /* We have a match! */ + if (ninputs==0) {/* diff=F in network statistic specification */ + CHANGE_STAT[0] += change; + } else { /* diff=T */ + for (j=0; j head */ + if (IS_OUTEDGE(head,tail)) { /* otherwise, no change occurs */ + change = edgestate ? -1.0 : 1.0 ; + for (j=0; j head */ + edgestateth = (!IS_OUTEDGE(head,tail)); + + for(node3=1;node3<=N_NODES;node3++){ + if((node3!=tail)&&(node3!=head)){ + sc = edgestateth + (!IS_OUTEDGE(node3,tail)); + if(sc < 2){ + sc += (!IS_OUTEDGE(tail,node3)); + if(sc < 2){ + sc += (!IS_OUTEDGE(node3,head)); + if(sc < 2){ + sc += (!IS_OUTEDGE(head,node3)); + if(sc < 2){ + change=0.0; + if (sc == 0 && edgestate == 0 ){--change;} + if (sc == 0 && edgestate == 1 ){++change;} + if (sc == 1 && edgestate == 0 ){++change;} + if (sc == 1 && edgestate == 1 ){--change;} + CHANGE_STAT[0] += change; + } + } + } + } + } + } + +} + +/******************** changestats: O ***********/ + +// A macro indicating whether x is in [from,to) +#define FROM_TO(x, from, to) ((x)>=(from) && (x)<(to)) + +/***************** + changestat: d_odegrange +*****************/ +C_CHANGESTAT_FN(c_odegrange) { + int j, echange; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1:+1; + Vertex tailodeg = OUT_DEG[tail]; + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } +} + +/***************** + changestat: d_odegrange_by_attr +*****************/ +C_CHANGESTAT_FN(c_odegrange_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 3*nstats values are in triples: (from, to, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1:1; + Vertex tailodeg = OUT_DEG[tail]; + int tailattr = INPUT_PARAM[3*N_CHANGE_STATS + tail - 1]; + for(j = 0; j < N_CHANGE_STATS; j++){ + Vertex from = INPUT_PARAM[3*j], to = INPUT_PARAM[3*j + 1]; + int testattr = INPUT_PARAM[3*j + 2]; + if (tailattr == testattr) /* we have tail attr match */ + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } +} + +/***************** + changestat: d_odegrange_w_homophily +*****************/ +C_CHANGESTAT_FN(c_odegrange_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are the values of odegrange + The values following the first 2*nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS*2 - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange = edgestate ? -1:1; + Vertex tailodeg=0, v; + STEP_THROUGH_OUTEDGES(tail, e, v) { tailodeg += (nodeattr[v]==tailattr); } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex from = INPUT_PARAM[2*j], to = INPUT_PARAM[2*j+1]; + CHANGE_STAT[j] += FROM_TO(tailodeg + echange, from, to) - FROM_TO(tailodeg, from, to); + } + } +} + +#undef FROM_TO + +/***************** + changestat: d_odegree +*****************/ +C_CHANGESTAT_FN(c_odegree) { + int j; + + /* *** don't forget tail -> head */ + int echange = edgestate ? -1 : 1; + Vertex taild = OUT_DEG[tail]; + + for(j=0; j < N_CHANGE_STATS; j++) { + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] = (taild + echange == deg) - (taild == deg); + } +} + + +/***************** + changestat: d_odegdist +*****************/ +C_CHANGESTAT_FN(c_odegdist) { + int echange = edgestate ? -1 : +1; + Vertex otd = OUT_DEG[tail], + ntd = otd + echange; + + if(ntd > N_CHANGE_STATS) cutoff_error(mtp); + + if(otd) CHANGE_STAT[otd-1]--; + if(ntd) CHANGE_STAT[ntd-1]++; +} + + +/***************** + changestat: d_odegree_by_attr +*****************/ +C_CHANGESTAT_FN(c_odegree_by_attr) { + /* The inputparams are assumed to be set up as follows: + The first 2*nstats values are in pairs: (degree, attrvalue) + The values following the first 2*nstats values are the nodal attributes. + */ + int j, echange, tailattr, testattr; + Vertex taildeg, d; + + /* *** don't forget tail -> head */ + echange=edgestate ? -1 : +1; + taildeg = OUT_DEG[tail]; + tailattr = INPUT_PARAM[2*N_CHANGE_STATS + tail - 1]; + for(j = 0; j < N_CHANGE_STATS; j++) { + d = (Vertex)INPUT_PARAM[2*j]; + testattr = INPUT_PARAM[2*j + 1]; + if (tailattr == testattr) { /* we have tail attr match */ + CHANGE_STAT[j] += (taildeg + echange == d) - (taildeg == d); + } + } +} + +/***************** + changestat: d_odegree_w_homophily +*****************/ +C_CHANGESTAT_FN(c_odegree_w_homophily) { + /* The inputparams are assumed to be set up as follows: + The first nstats values are the values of degree + The values following the first nstats values are the nodal attributes. + */ + int j; + double *nodeattr; + Edge e; + + nodeattr = mtp->inputparams + N_CHANGE_STATS - 1; + + /* *** don't forget tail -> head */ + int tailattr = nodeattr[tail], headattr = nodeattr[head]; + if (tailattr == headattr) { /* They match; otherwise don't bother */ + int echange=edgestate ? -1 : +1; + Vertex taildeg=0, tmp; + STEP_THROUGH_OUTEDGES(tail, e, tmp){ + taildeg += (nodeattr[tmp]==tailattr); + } + for(j = 0; j < N_CHANGE_STATS; j++) { + Vertex deg = INPUT_PARAM[j]; + CHANGE_STAT[j] += (taildeg + echange == deg) - (taildeg == deg); + } + } +} + +/***************** + changestat: d_opentriad +*****************/ +C_CHANGESTAT_FN(c_opentriad) { + + /* *** don't forget tail -> head */ + Vertex node3; + Edge change = 0, e; + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + + // -3 * triangles + + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + change += IS_UNDIRECTED_EDGE(node3,tail); + } + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + change += IS_UNDIRECTED_EDGE(node3,tail); + } + CHANGE_STAT[0] += change * (edgestate ? 3.0 : -3.0); + + + // +1 * 2-stars + + Vertex taild = DEG(tail) - edgestate; + Vertex headd = DEG(head) - edgestate; + change = taild + headd; + CHANGE_STAT[0] += (edgestate ? -change : change); + +} + +/***************** + changestat: d_ostar +*****************/ +C_CHANGESTAT_FN(c_ostar) { + double change, headd=0.0; + int j, kmo; + Edge e; + Vertex node3; + int ninputs, nstats; + double headattr; + + ninputs = (int)N_INPUT_PARAMS; + nstats = (int)N_CHANGE_STATS; + + /* *** don't forget tail -> head */ + if(ninputs>nstats){ + /* match on attributes */ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headattr = INPUT_ATTRIB[head-1]; + if(headattr == INPUT_ATTRIB[tail-1]){ + headd = -(int)edgestate; + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* step through outedges of head */ + if(headattr == INPUT_ATTRIB[node3-1]){++headd;} + } + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } + }else{ + /* edgestate is 1 if edge exists and will disappear + edgestate is 0 if edge DNE and will appear */ + headd = OUT_DEG[tail] - edgestate; + for(j=0; j < N_CHANGE_STATS; j++) { + kmo = ((int)INPUT_PARAM[j]) - 1; + change = CHOOSE(headd, kmo); + CHANGE_STAT[j] += (edgestate ? - change : change); + } + } +} + +/***************** + changestat: d_odegreepopularity +*****************/ +C_CHANGESTAT_FN(c_odegreepopularity) { + double change; + Vertex deg=0; + + /* *** don't forget tail -> head */ + change = 0.0; + deg = (double)(OUT_DEG[tail]); + if(edgestate){ + change -= sqrt(deg); + change += (deg-1.0)*(sqrt(deg-1.0)-sqrt(deg)); + }else{ + change += sqrt(deg+1.0); + change += deg*(sqrt(deg+1.0)-sqrt(deg)); + } +CHANGE_STAT[0]=change; +} + +/******************** changestats: P ***********/ +/***************** + changestat: d_pdegcor +*****************/ +D_CHANGESTAT_FN(d_pdegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + mtp->dstats[0] -= current; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_pdegcor) { + Vertex taildeg, headdeg; + Edge e; + double mu, mu2, mutail, mutail2, sigma2, sigmatail2, cross; + + mu = 0.0; + mu2 = 0.0; + mutail = 0.0; + mutail2 = 0.0; + cross = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = OUT_DEG[tail]; + headdeg = IN_DEG[head]; + mu += (double)(headdeg); + mutail += (double)(taildeg); + mu2 += (double)(headdeg*headdeg); + mutail2 += (double)(taildeg*taildeg); + cross += taildeg*headdeg; + } + } + mu = mu / (N_EDGES); + mutail = mutail / (N_EDGES); + sigma2 = mu2/(N_EDGES) - mu*mu; + sigmatail2 = mutail2/(N_EDGES) - mutail*mutail; + CHANGE_STAT[0] = (cross / (N_EDGES) - mutail*mu) / sqrt(sigma2*sigmatail2); +} + +/******************** changestats: R ***********/ +/***************** + changestat: d_rdegcor +*****************/ +D_CHANGESTAT_FN(d_rdegcor) { + int i; + double current; + + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ + current = mtp->dstats[0]; + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } + (*(mtp->s_func))(mtp, nwp); /* Call s_??? function */ +// CHANGE_STAT[0] = mtp->dstats[0] - current; +// Rprintf("c %f p %f",current,mtp->dstats[0]); + mtp->dstats[0] -= current; +// Rprintf(" p-c %f\n",mtp->dstats[0]); + FOR_EACH_TOGGLE(i) { TOGGLE(TAIL(i), HEAD(i)); } +} +S_CHANGESTAT_FN(s_rdegcor) { + Vertex taildeg, headdeg; + Edge e; + double mu, mu2, sigma2, cross; + Vertex tailrank, headrank; + Vertex *ndeg=R_Calloc(N_NODES+1, Vertex); + + for(Vertex tail=0; tail <= N_NODES; tail++) { ndeg[tail]=0; } + for(Vertex tail=0; tail < N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + ndeg[taildeg+1]++; + ndeg[headdeg+1]++; + } + } +for(Vertex tail=1; tail <= N_NODES; tail++) { + ndeg[tail] += ndeg[tail-1]; +} +// Rprintf("tail %d taildeg[tail] %d \n",tail,ndeg[tail]);} + + mu = 0.0; + mu2 = 0.0; + cross = 0.0; + for(Vertex tail=1; tail <= N_NODES; tail++) { + Vertex head; + STEP_THROUGH_OUTEDGES(tail, e, head) { /* step through outedges of tail */ + taildeg = DEG(tail); + headdeg = DEG(head); + tailrank = (ndeg[taildeg+1]+ndeg[taildeg+2]+1)*0.5; + headrank = (ndeg[headdeg+1]+ndeg[headdeg+2]+1)*0.5; + mu += (double)(tailrank + headrank); + mu2 += (double)(tailrank*tailrank + headrank*headrank); + cross += 2.0*tailrank*headrank; + } + } + mu = mu / (2.0*N_EDGES); + sigma2 = mu2/(2.0*N_EDGES) - mu*mu; + CHANGE_STAT[0] = (cross / (2.0*N_EDGES) - mu*mu) / sigma2; + R_Free(ndeg); +} + +/***************** + changestat: d_simmelian +*****************/ +C_CHANGESTAT_FN(c_simmelian) { + Edge e; + Vertex change, node3; + + /* *** don't forget tail -> head */ + + if(IS_OUTEDGE(head, tail)){ + change = 0; + + STEP_THROUGH_OUTEDGES(head, e, node3) /* step through outedges of head */ + { + if (node3 != tail + && IS_OUTEDGE(node3, tail) + && IS_OUTEDGE(tail, node3) + && IS_OUTEDGE(node3, head) + ){++change;} + } + + CHANGE_STAT[0] += edgestate ? -(double)change : (double)change; + } + +} + +/***************** + changestat: d_simmelianties +*****************/ +C_CHANGESTAT_FN(c_simmelianties) { + Edge e, e2; + Vertex change, node3, node4, first, htflag; + + /* *** don't forget tail -> head */ + + if(IS_OUTEDGE(head, tail)){ + change = htflag = 0; + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if (node3 != tail + && IS_OUTEDGE(node3, tail) && IS_OUTEDGE(tail, node3) && IS_OUTEDGE(node3, head)){ + htflag=1; /* tail, head is itself in a simmelian triple (along with head, tail)*/ + first = 1; + /* Find out whether (tail, node3) is in any other simmelian triple */ + STEP_THROUGH_OUTEDGES (tail, e2, node4) { /* step through outedges of tail */ + if (node4 != head && node4 != node3 && IS_OUTEDGE(node4, tail) + && IS_OUTEDGE(node4, node3) && IS_OUTEDGE(node3, node4)){ + first = 0; + } + } + if(first){++change;} + first = 1; + /* Find out whether (head, node3) is in any other simmelian triple */ + STEP_THROUGH_OUTEDGES (head, e2, node4) { /* step through outedges of head */ + if (node4 != tail && node4 != node3 && IS_OUTEDGE(node4, head) + && IS_OUTEDGE(node4, node3) && IS_OUTEDGE(node3, node4)) { + first = 0; + } + } + if(first){++change;} + } + } + change += htflag; + change = 2*change; /* All changes must happen in pairs here; no tie can + be counted without its opposite */ + CHANGE_STAT[0] += edgestate ? -(double)change : (double)change; + } +} + +/******************** changestats: T ***********/ + +/***************** + changestat: d_threetrail +*****************/ +C_CHANGESTAT_FN(c_threetrail) { + int j, k, change, dchange[4]; + Edge e; + Vertex node3; + /* The four values of dchange represent the four different types of + directed threetrails oriented so that the middle step is always + "right" (R). In order: RRR, RRL, LRR, LRL + i.e., >>> >>< <>> <>< */ + + + /* *** don't forget tail -> head */ + /* Step A: Count threetrails in which tail->head is the middle edge */ + dchange[0] = IN_DEG[tail] * OUT_DEG[head]; /* R then R; may count head->tail->head->tail */ + dchange[1] = IN_DEG[tail] * (IN_DEG[head]-edgestate); /* R then L */ + dchange[2] = (OUT_DEG[tail]-edgestate) * OUT_DEG[head]; /* L then R */ + dchange[3] = (OUT_DEG[tail]-edgestate) * (IN_DEG[head]-edgestate); /* L then L */ + /* Step B: Count threetrails where tail is one endpoint */ + STEP_THROUGH_OUTEDGES(head, e, node3) { /* tail->head->node3-x which means -RL */ + dchange[1] += IN_DEG[node3]-1; /* RRL; subtract 1 for head itself */ + dchange[0] += OUT_DEG[node3]; /* RRR; possibly counted tail->head->tail->head */ + } + STEP_THROUGH_INEDGES(head, e, node3) { /* x-node3->head<-tail which means -RL*/ + if (node3 != tail) { + dchange[3] += OUT_DEG[node3]-1; /* LRL; subtract 1 for head itself */ + dchange[1] += IN_DEG[node3]; /* RRL */ + } + } + /* Step C: Count threetrails where head is one endpoint */ + STEP_THROUGH_INEDGES(tail, e, node3) { /* x-node3->tail->head which means -RR */ + dchange[2] += OUT_DEG[node3]-1; /* LRR; subtract 1 for tail itself */ + dchange[0] += IN_DEG[node3]; /* RRR; possibly counted tail->head->tail->head */ + } + STEP_THROUGH_OUTEDGES(tail, e, node3) { /* head<-tail->node3-x which means LR- */ + if (node3 != head) { + dchange[3] += IN_DEG[node3]-1; /* LRL; subtract 1 for tail itself */ + dchange[2] += OUT_DEG[node3]; /* LRR */ + } + } + /* Finally, correct for overcounted head->tail->head->tail and tail->head->tail->head */ + if (DIRECTED) { + dchange[0] -= IS_INEDGE(tail, head) * (1 + 2 * edgestate); + /* head->tail->head->tail is counted in A whenever IS_INEDGE(tail,head) but + TT->head->tail->head is only counted in B and C when also edgestate */ + for (j = 0; j < N_INPUT_PARAMS; j++) { + k = (int) INPUT_PARAM[j]; + CHANGE_STAT[j] += (edgestate ? -dchange[k-1] : dchange[k-1]); + } + } + else { /* Undirected case; don't need head->tail->head->tail correction */ + change = dchange[0] + dchange[1] + dchange[2] + dchange[3]; + CHANGE_STAT[0] += (edgestate ? -change : change); + } +} + +/***************** + changestat: d_transitive +*****************/ +C_CHANGESTAT_FN(c_transitive) { + Edge e; + Vertex node2; + double change; + + /* *** don't forget tail -> head */ + change = 0.0; /* change should become the number of transitive triples + a->b, b->c, a->c in which tail->head is found */ + + STEP_THROUGH_OUTEDGES(head, e, node2) { /* step through outedges of head */ + if (tail != node2 && IS_OUTEDGE(tail, node2)){ + change = change + 1.0; /* Here we have tail->head, head->node2, tail->node2 */ + } + } + STEP_THROUGH_INEDGES(head, e, node2) { /* step through inedges of head */ + if (tail != node2) { + change = change + IS_OUTEDGE(tail, node2) + IS_OUTEDGE(node2, tail); + /* Here we have tail->head and node2->head, with either node2->tail or tail->node2 */ + } + } +// STEP_THROUGH_INEDGES(tail, e, node2) { /* step through inedges of tail */ +// if (node2 != head){ +// if (!IS_OUTEDGE(node2, head)){ +// change = change - 1.0; +// } +// } +// } + CHANGE_STAT[0] += edgestate ? -change : change; +// Rprintf("tail %d head %d edgestate %d change %f C_S[0]=%f\n", tail, head, change,CHANGE_STAT[0]); +} + +C_CHANGESTAT_FN(c_transitiveties) { + int echange, ochange; + int L2th, L2tu, L2uh; + double cumchange; + double tailattr; + + + /* *** don't forget tail -> head */ + cumchange=0.0; + L2th=0; + ochange = GETWT(tail, head) ? -1 : 0; + echange = 2*ochange + 1; + if(N_INPUT_PARAMS>0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(tail, u) && (tailattr == INPUT_ATTRIB[u-1])){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v) && (tailattr == INPUT_ATTRIB[v-1])){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_INEDGES(head, e, u, { + if (GETWT(tail, u) && (tailattr == INPUT_ATTRIB[u-1])){ + L2th++; + } + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head) && (tailattr == INPUT_ATTRIB[v-1])){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + }else{ /* no attributes */ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(tail, u)){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v)){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_INEDGES(head, e, u, { + if (GETWT(tail, u)){ + L2th++; + } + if (GETWT(u, tail)){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head)){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + + cumchange += (L2th>0) ; +// Rprintf("L2th %d echange %d cumchange %f tail %d head %d\n", L2th, echange, cumchange,tail,head); + cumchange = echange*cumchange; + (CHANGE_STAT[0]) += cumchange; +} + +C_CHANGESTAT_FN(c_cyclicalties) { + int echange, ochange; + int L2th, L2tu, L2uh; + double cumchange; + double tailattr; + + + /* *** don't forget tail -> head */ + cumchange=0.0; + L2th=0; + ochange = GETWT(tail, head) ? -1 : 0; + echange = 2*ochange + 1; + if(N_INPUT_PARAMS>0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]){ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v) && (tailattr == INPUT_ATTRIB[v-1])){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through inedges of head */ + + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2th++; + } + if (GETWT(u, tail) && (tailattr == INPUT_ATTRIB[u-1])){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head) && (tailattr == INPUT_ATTRIB[v-1])){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + }else{ /* no attributes */ + /* step through outedges of head */ + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail)){ + L2tu=ochange; + /* step through inedges of u */ + EXEC_THROUGH_INEDGES(u, f, v, { + if(GETWT(tail, v)){ + L2tu++; + if(L2tu>0) {break;} + } + }); + cumchange += (L2tu==0); + } + }); + /* step through outedges of head */ + + EXEC_THROUGH_OUTEDGES(head, e, u, { + if (GETWT(u, tail)){ + L2th++; + } + if (GETWT(u, tail)){ + L2uh=ochange; + /* step through outedges of u */ + EXEC_THROUGH_OUTEDGES(u, f, v, { + if(GETWT(v, head)){ + L2uh++; + if(L2uh>0) {break;} + } + }); + cumchange += (L2uh==0) ; + } + }); + } + + cumchange += (L2th>0) ; +// Rprintf("L2th %d echange %d cumchange %f tail %d head %d\n", L2th, echange, cumchange,tail,head); + cumchange = echange*cumchange; + (CHANGE_STAT[0]) += cumchange; +} + +/***************** + changestat: d_triadcensus +*****************/ +C_CHANGESTAT_FN(c_triadcensus) { + int j, a, b, c, d, e, edgecount, t300, + t210, t120C, t120U, t120D, t201, t030C, t030T, t111U, + t111D, t021C, t021U, t021D, t102, t012, t003; + Vertex triadtype, node3; + + /* *** don't forget tail -> head */ + if (DIRECTED) { + /* directed version */ + t300 = 0; + t210 = 0; + t120C = 0; t120U = 0; t120D = 0; t201 = 0; + t030C = 0; t030T = 0; t111U = 0; t111D = 0; + t021C = 0; t021U = 0; t021D = 0; t102 = 0; + t012 = 0; + + if ( (MIN_OUTEDGE(head) != 0) || + (MIN_INEDGE(head) != 0) || + (MIN_OUTEDGE(tail) != 0) || + (MIN_INEDGE(tail) != 0) ) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = (IS_OUTEDGE(head, tail)); + b = (IS_OUTEDGE(head, node3)); + c = (IS_OUTEDGE(node3, head)); + d = (IS_OUTEDGE(node3, tail)); + e = (IS_OUTEDGE(tail, node3)); + edgecount = (a + b + c + d + e); + + switch(edgecount) { + case 0: /* 012 */ + ++t012; + + case 1: { /* 021C, 021U, 021D, 102 */ + if ((b == 1) || (d == 1)) + ++t021C; + if (c == 1) + ++t021U; + if (e == 1) + ++t021D; + if (a == 1) + ++t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if ((b + d) == 2) + ++t030C; + if (((b + e) == 2) || ((c + d) == 2) || ((c + e) == 2)) + ++t030T; + if (((a + b) == 2) || ((a + e) == 2) || ((d + e) == 2)) + ++t111U; + if (((a + c) == 2) || ((a + d) == 2) || ((b + c) == 2)) + ++t111D; + } + break; + + case 3: { /* 120C, 120U, 120D, 201 */ + if (a == 1) { + if (((b + d) == 2) || ((c + e) == 2)) + ++t120C; + if ((b + e) == 2) + ++t120U; + if ((c + d) == 2) + ++t120D; + if (((b + c) == 2) || ((d + e) == 2)) + ++t201; + } else { + if (b == 1) { + if (((c + d) == 2) || ((d + e) == 2)) + ++t120C; + if ((c + e) == 2) + ++t120D; + } else { + ++t120U; + } + } + } + break; + + case 4: /* 210 */ + ++t210; + break; + + case 5: /* 300 */ + ++t300; + break; + } + + switch(edgecount) { + case 1: /* 102, 021D, 021U, 021C */ + --t012; + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + if (((a + c) == 2) || ((a + e) == 2) || ((b + d) == 2) || + ((c + e) == 2)) + --t021C; + if (((a + d) == 2) || ((b + e) == 2)) + --t021U; + if (((a + b) == 2) || ((c + d) == 2)) + --t021D; + if (((b + c) == 2) || ((d + e) == 2)) + --t102; + } + break; + + case 3: { /* 201, 120D, 120U, 120C */ + if (a == 1) { + if ((c + e) == 2) + --t030C; + if (((c + d) == 2) || ((b + e) == 2) || ((b + d) == 2)) + --t030T; + if ((b + c) == 2) + --t111U; + if ((d + e) == 2) + --t111D; + } else { + if (b == 1) { + if ((c + d) == 2) + --t111U; + if (((c + e) == 2) || ((d + e) == 2)) + --t111D; + } else { + --t111U; + } + } + } + break; + + case 4: { /* 210 */ + if (a == 1) + { + if (((b + c + e) == 3) || ((c + d + e) == 3)) + --t120C; + if ((b + c + d) == 3) + --t120U; + if ((b + d + e) == 3) + --t120D; + } else { + if ((b + c + d + e) == 4) + --t201; + } + } + break; + + case 5: /* 300 */ + --t210; + break; + } + } + } /* ****** move to next node3 ******** */ + }else{ + t012 = t012 + (N_NODES - 2); + } + + for(j = 0; j < N_CHANGE_STATS; j++) { + triadtype = (Vertex)INPUT_PARAM[j]; + + switch(triadtype) { /* SEARCH_ON_THIS_TO_TRACK_DOWN_TRIADCENSUS_CHANGE + to undo triadcensus change, change - to plus in + next two lines: */ + case 1: t003 = -(t300+t210+t120C+t120U+t120D+t201+t030C+t030T); + t003 = t003-(t111U+t111D+t021C+t021U+t021D+t102+t012); + CHANGE_STAT[j] += edgestate ? -(double)t003 : (double)t003; + break; + case 2: CHANGE_STAT[j] += edgestate ? -(double)t012 : (double)t012; + break; + case 3: CHANGE_STAT[j] += edgestate ? -(double)t102 : (double)t102; + break; + case 4: CHANGE_STAT[j] += edgestate ? -(double)t021D : (double)t021D; + break; + case 5: CHANGE_STAT[j] += edgestate ? -(double)t021U : (double)t021U; + break; + case 6: CHANGE_STAT[j] += edgestate ? -(double)t021C : (double)t021C; + break; + case 7: CHANGE_STAT[j] += edgestate ? -(double)t111D : (double)t111D; + break; + case 8: CHANGE_STAT[j] += edgestate ? -(double)t111U : (double)t111U; + break; + case 9: CHANGE_STAT[j] += edgestate ? -(double)t030T : (double)t030T; + break; + case 10: CHANGE_STAT[j] += edgestate ? -(double)t030C : (double)t030C; + break; + case 11: CHANGE_STAT[j] += edgestate ? -(double)t201 : (double)t201; + break; + case 12: CHANGE_STAT[j] += edgestate ? -(double)t120D : (double)t120D; + break; + case 13: CHANGE_STAT[j] += edgestate ? -(double)t120U : (double)t120U; + break; + case 14: CHANGE_STAT[j] += edgestate ? -(double)t120C : (double)t120C; + break; + case 15: CHANGE_STAT[j] += edgestate ? -(double)t210 : (double)t210; + break; + case 16: CHANGE_STAT[j] += edgestate ? -(double)t300 : (double)t300; + break; + } + } + } else { + /* undirected */ + + /* *** don't forget tail -> head */ + t300 = 0; t201 = 0; t102 = 0; t012 = 0; + + if ( (MIN_OUTEDGE(head) != 0) || + (MIN_INEDGE(head) != 0) || + (MIN_OUTEDGE(tail) != 0) || + (MIN_INEDGE(tail) != 0) ) { + + /* ****** loop through node3 ****** */ + for (node3=1; node3 <= N_NODES; node3++) { + if (node3 != tail && node3 != head) { + a = (IS_UNDIRECTED_EDGE(node3, head)); + b = (IS_UNDIRECTED_EDGE(node3, tail)); + edgecount = (a + b); + + switch(edgecount) { + case 0: { /* 012 */ + ++t102; + --t012; + } + break; + + case 1: { /* 021C, 021U, 021D, 102 */ + ++t201; + --t102; + } + break; + + case 2: { /* 030C, 030T, 111U, 111D */ + ++t300; + --t201; + } + break; + + } + } + + } /* ****** move to next node3 ******** */ + } else { + t102 = t102 + (N_NODES - 2); + } + + for(j = 0; j < N_CHANGE_STATS; j++) { + triadtype = (Vertex)INPUT_PARAM[j]; + + switch(triadtype) { /* SEARCH_ON_THIS_TO_TRACK_DOWN_TRIADCENSUS_CHANGE + to undo triadcensus change, change - to plus in + next line: */ + case 1: t003 = -(t102+t201+t300); + CHANGE_STAT[j] += edgestate ? -(double)t003 : (double)t003; + break; + case 2: CHANGE_STAT[j] += edgestate ? -(double)t102 : (double)t102; + break; + case 3: CHANGE_STAT[j] += edgestate ? -(double)t201 : (double)t201; + break; + case 4: CHANGE_STAT[j] += edgestate ? -(double)t300 : (double)t300; + break; + } + } + } +} + + +/***************** + changestat: d_tripercent +*****************/ +C_CHANGESTAT_FN(c_tripercent) { + Edge e, e2; + Vertex node1, node2, node3; + int j; + Edge triwith, triwithout; + Edge degreewith, degreewithout, twostarwith, twostarwithout; + int ninputs = N_INPUT_PARAMS - N_NODES; + int MatchingOnAttribute = (ninputs>0); + double *attr=INPUT_PARAM, ratiowith, ratiowithout; + + if (MatchingOnAttribute) + attr = INPUT_PARAM + (ninputs-1); /* ptr to vertex attributes */ + + /* *** don't forget tail -> head */ + if (!edgestate) TOGGLE(tail, head); /* turn on the edge if it's missing */ + for (j=0; j < MAX(1, ninputs); j++) { + /* Count triangles with and without proposed edge */ + /* Simultaneously, find degree (use matching if necessary) with and without */ + triwith = triwithout = twostarwith = twostarwithout = 0; + for (node1 = 1; node1 <= N_NODES; node1++) { + degreewith = degreewithout = 0; + if (ninputs < 2 || EQUAL(attr[node1],INPUT_PARAM[j])) { + STEP_THROUGH_OUTEDGES(node1, e, node2) { + /* inside this loop, node1 < node2 always */ + if (!MatchingOnAttribute || EQUAL(attr[node1],attr[node2])) { + /* increment degree counter */ + ++degreewith; + if (node1!=tail || node2!=head) ++degreewithout; + STEP_THROUGH_OUTEDGES(node2, e2, node3) { + /* inside this loop, node1 < node2 < node3 always */ + if (!MatchingOnAttribute || EQUAL(attr[node2],attr[node3])) { + if (IS_OUTEDGE(node1, node3)) { + ++triwith; + if ((tail!=node1||head!=node2)&&(tail!=node2||head!=node3)&&(tail!=node1||head!=node3)) + ++triwithout; + } + } + } + } + } + STEP_THROUGH_INEDGES(node1, e, node2) { + /* inside this loop, node2 < node1 always. */ + /* We only do this to find correct degree for node1; */ + /* for triangles, node1 < node2 head */ + edgemult = edgestate ? -1.0 : 1.0; + change = 0; + if(N_INPUT_PARAMS > 0){ /* match on attributes */ + tailattr = INPUT_ATTRIB[tail-1]; + if(tailattr == INPUT_ATTRIB[head-1]) { + STEP_THROUGH_OUTEDGES(head, e, node3) { /* step through outedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_INEDGE(node3, tail); + } + STEP_THROUGH_INEDGES(head, e, node3) { /* step through inedges of head */ + if(tailattr == INPUT_ATTRIB[node3-1]) + change += IS_OUTEDGE(node3, tail) + IS_INEDGE(node3, tail); + } + if(N_CHANGE_STATS > 1) { /* diff = TRUE; matches must be tabled */ + for (j=0; j= d) - (b1deg >= d); + } +} + + +/***************** + changestat: c_b2mindegree +*****************/ +C_CHANGESTAT_FN(c_b2mindegree) { + int echange = edgestate ? -1 : 1; + Vertex b2deg = IN_DEG[head]; + for(unsigned int j = 0; j < N_CHANGE_STATS; j++) { + Vertex d = INPUT_PARAM[j]; + CHANGE_STAT[j] += (b2deg + echange >= d) - (b2deg >= d); + } +} From bcd881262d18d9b1c6f2c64333fd2fb267fb92da Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:45:45 -0400 Subject: [PATCH 13/18] --Set up 'changestats_distance.c' to be appended to the 'changestats_dyad_ind.c.template.do_not_include_directly.h' placeholder. --- ...cat-changestats_dyad_ind.c.template.do_not_include_directly.h} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestats_distance.c => concat-changestats_dyad_ind.c.template.do_not_include_directly.h} (100%) diff --git a/src/changestats_distance.c b/src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h similarity index 100% rename from src/changestats_distance.c rename to src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h From 9d12e26ccfe1f1c9dc09cd6e9af3c6738fd474a4 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 19:45:45 -0400 Subject: [PATCH 14/18] --Set up 'changestats_dyad_ind.c.template.do_not_include_directly.h' to be appended to the 'changestats_dyad_ind.c.template.do_not_include_directly.h' placeholder. --- ...cat-changestats_dyad_ind.c.template.do_not_include_directly.h} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{changestats_dyad_ind.c.template.do_not_include_directly.h => concat-changestats_dyad_ind.c.template.do_not_include_directly.h} (100%) diff --git a/src/changestats_dyad_ind.c.template.do_not_include_directly.h b/src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h similarity index 100% rename from src/changestats_dyad_ind.c.template.do_not_include_directly.h rename to src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h From 7a963c83a1c630eff0958349df7ccafd49bfa039 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 20:12:28 -0400 Subject: [PATCH 15/18] Concatenate files 'changestats_dyad_ind.c.template.do_not_include_directly.h' 'changestats_distance.c' into 'changestats_dyad_ind.c.template.do_not_include_directly.h'. --- ...estats_dyad_ind.c.template.do_not_include_directly.h} | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) rename src/{concat-changestats_dyad_ind.c.template.do_not_include_directly.h => changestats_dyad_ind.c.template.do_not_include_directly.h} (98%) diff --git a/src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h b/src/changestats_dyad_ind.c.template.do_not_include_directly.h similarity index 98% rename from src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h rename to src/changestats_dyad_ind.c.template.do_not_include_directly.h index a6df03509..078991b85 100644 --- a/src/concat-changestats_dyad_ind.c.template.do_not_include_directly.h +++ b/src/changestats_dyad_ind.c.template.do_not_include_directly.h @@ -439,7 +439,7 @@ ETYPE(C_CHANGESTAT_FN)(SVARIANT(c_sociality)) { /***************** changestat: c_distance *****************/ -C_CHANGESTAT_FN(c_distance) { +ETYPE(C_CHANGESTAT_FN)(SVARIANT(c_distance)) { Vertex t,h; int j; int nv,dim,logd,sphd; @@ -495,9 +495,6 @@ C_CHANGESTAT_FN(c_distance) { from below by logmind (after adding the offset).*/ if(logd) dis=log(MAX(logmind,dis+logdoff)); - /*Check for an edge, and update accordingly*/ - if(DIRECTED) - CHANGE_STAT[0] += IS_OUTEDGE(tail,head) ? -dis : dis; - else - CHANGE_STAT[0] += IS_UNDIRECTED_EDGE(tail,head) ? -dis : dis; + + CHANGE_STAT[0] = ECHANGE(dis); } From 83c73f62d8483f4a6f708c68e0a7350fa133baee Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 3 Aug 2026 20:56:35 -0400 Subject: [PATCH 16/18] Added valued distance() term initializer. --- R/InitWtErgmTerm.R | 15 +++++++++++++++ ...-9bc71012.Rd => distance-ergmTerm-d5141986.Rd} | 14 +++++++++++++- 2 files changed, 28 insertions(+), 1 deletion(-) rename man/{distance-ergmTerm-9bc71012.Rd => distance-ergmTerm-d5141986.Rd} (93%) diff --git a/R/InitWtErgmTerm.R b/R/InitWtErgmTerm.R index df2666d70..3878b7f05 100644 --- a/R/InitWtErgmTerm.R +++ b/R/InitWtErgmTerm.R @@ -301,6 +301,21 @@ InitWtErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm") binary_dind_wrap("diff", nw, a, ..., version=version) } +#' @templateVar name distance +#' @template ergmTerm-rdname +#' @usage +#' # valued: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +#' # log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1, form="sum") +#' @template ergmTerm-form +InitWtErgmTerm.distance <- function(nw, arglist, ...) { + a <- check.ErgmTerm(nw, arglist, + varnames = c("coord", "metric", "sphere", "radius", "log", "mindist", "distoff", "scale", "pow", "form"), + vartypes = c("numeric,matrix,data.frame,character", "numeric", "logical", "numeric", "logical", "numeric", "numeric", "numeric", "numeric", "character"), + required = c(TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE), + defaultvalues = list(NULL,2,FALSE,6371.0087714,TRUE,1e-5,0,1,1, "sum")) + + binary_dind_wrap("distance", nw, a, ...) +} #' @templateVar name edgecov #' @template ergmTerm-rdname #' @usage diff --git a/man/distance-ergmTerm-9bc71012.Rd b/man/distance-ergmTerm-d5141986.Rd similarity index 93% rename from man/distance-ergmTerm-9bc71012.Rd rename to man/distance-ergmTerm-d5141986.Rd index f48420796..3033d8b5f 100644 --- a/man/distance-ergmTerm-9bc71012.Rd +++ b/man/distance-ergmTerm-d5141986.Rd @@ -1,12 +1,16 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/InitErgmTerm.R +% Please edit documentation in R/InitErgmTerm.R, R/InitWtErgmTerm.R \name{distance-ergmTerm} \alias{distance-ergmTerm} \alias{InitErgmTerm.distance} +\alias{InitWtErgmTerm.distance} \title{Inter-point distances} \usage{ # binary: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, # log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1) + +# valued: distance(coord, metric=2, sphere=FALSE, radius=6371.0087714, +# log=TRUE, mindist=1e-5, distoff=0, scale=1, pow=1, form="sum") } \arguments{ \item{coord}{node by dimension coordinate matrix, name of a network @@ -40,6 +44,14 @@ operations (notably, offsetting, logging, and thresholding)} \item{pow}{power to which scaled distances should be raised prior to other other operations (notably, offsetting, logging, and thresholding)} + +\item{form}{how to aggregate tie values in a valued ERGM: \code{"sum"} +(the default) for a statistic of the form \eqn{\sum_{i,j} x_{i,j} +y_{i,j}}{sum[i,j] x[i,j]*y[i,j]}, where \eqn{y_{i,j}}{y[i,j]} is +the value of dyad \eqn{(i,j)} and \eqn{x_{i,j}}{x[i,j]} is the +term's covariate associated with it; and \code{"nonzero"} with the +edge considered to be present if its value is not 0. See +\code{\link{ergmTerm}} for more information.} } \description{ This term adds a single statistic to the model whose value is From 1903ad726eb8f5785a062a456068e73e9bdd51dc Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Fri, 21 Aug 2026 12:32:52 -0400 Subject: [PATCH 17/18] In the distance() term example, replaced \dontrun{} with \donttest{}. --- R/InitErgmTerm.R | 2 +- man/distance-ergmTerm-d5141986.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index a02c43af5..23784ebdd 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -2526,7 +2526,7 @@ InitErgmTerm.diff <- function(nw, arglist, ..., version=packageVersion("ergm")) #' This term can be used for directed or undirected networks. #' #' @examples -#' \dontrun{ #We turn these off by default, b/c they are a bit slow +#' \donttest{ #' #Create an example network, in a two-dimensional space #' # Effective SIF is 1/(1 + exp(-3 + 2 log(d))); apx inverse square #' n <- 300 diff --git a/man/distance-ergmTerm-d5141986.Rd b/man/distance-ergmTerm-d5141986.Rd index 3033d8b5f..2fc18b7d0 100644 --- a/man/distance-ergmTerm-d5141986.Rd +++ b/man/distance-ergmTerm-d5141986.Rd @@ -136,7 +136,7 @@ Note that both \code{distoff} and \code{mindist} are ignored when This term can be used for directed or undirected networks. } \examples{ -\dontrun{ #We turn these off by default, b/c they are a bit slow +\donttest{ #Create an example network, in a two-dimensional space # Effective SIF is 1/(1 + exp(-3 + 2 log(d))); apx inverse square n <- 300 From 15a386c7a646b0eaf1828599df9e80582307bbab Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Fri, 21 Aug 2026 17:33:05 -0700 Subject: [PATCH 18/18] Bumped the expected number of valued terms in the tests. --- tests/testthat/test-ergm-term-doc.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-ergm-term-doc.R b/tests/testthat/test-ergm-term-doc.R index b62c354a1..5677a5618 100644 --- a/tests/testthat/test-ergm-term-doc.R +++ b/tests/testthat/test-ergm-term-doc.R @@ -44,7 +44,7 @@ test_that("test search ergm term", { expect_equal(length(search.ergmTerms(keywords = 'bipartite', packages='ergm')), 36) ## expect_gt(length(search.ergmTerms(keywords = 'valued')), 44) - expect_equal(length(search.ergmTerms(keywords = 'valued', packages='ergm')), 50) + expect_equal(length(search.ergmTerms(keywords = 'valued', packages='ergm')), 51) ## expect_gt(length(search.ergmTerms(keywords = 'valued', packages=c('ergm', 'ergm.count'))), 44) })