From f6bded4b1c3a369bd7b1f609d798ccac9461bb07 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 1 Jul 2026 08:02:38 +0000 Subject: [PATCH 1/7] Initial plan From da5f5acf39d179d4ab56e66674d51e938c7bef52 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 1 Jul 2026 08:10:26 +0000 Subject: [PATCH 2/7] Port shared-partner changestats to C++ --- src/MHproposals_triadic.cpp | 18 +- src/changestats_dgw_sp.c | 500 ------------------------------------ src/changestats_dgw_sp.cpp | 249 ++++++++++++++++++ src/changestats_dgw_sp.h | 45 ++++ 4 files changed, 296 insertions(+), 516 deletions(-) delete mode 100644 src/changestats_dgw_sp.c create mode 100644 src/changestats_dgw_sp.cpp diff --git a/src/MHproposals_triadic.cpp b/src/MHproposals_triadic.cpp index 91c02e40f..c650d783a 100644 --- a/src/MHproposals_triadic.cpp +++ b/src/MHproposals_triadic.cpp @@ -134,24 +134,10 @@ extern "C" MH_P_FN(Mp_SPDyad){ Dyad oldtd = kh_size(spcache), newtd = oldtd; - // The following is setting up to use macros developed for the *sp - // terms. Rboolean edgeflag = nw(p.tail[0], p.head[0]); - int echange = edgeflag ? -1 : +1; Vertex tail = p.tail[0], head = p.head[0]; - -#define sp_nonzero newtd += (L2 + echange != 0) - (L2 != 0); - - switch(type){ - case L2UTP: dspUTP_change(sp_nonzero, ); break; - case L2OTP: dspOTP_change(sp_nonzero, ); break; - case L2ITP: dspITP_change(sp_nonzero, ); break; - case L2OSP: dspOSP_change(sp_nonzero, ); break; - case L2ISP: dspISP_change(sp_nonzero, ); break; - default: error("In ergm:Mp_SPDyad(), an unsupported type of triad: %d.", type); - } - -#undef sp_nonzero + if(type == L2RTP) error("In ergm:Mp_SPDyad(), an unsupported type of triad: %d.", type); + newtd += ergm::sp::dsp_nonzero_change(type, tail, head, nwp, edgeflag, spcache); // q(y | y*) / q(y* | y) = 1/TD(y*) / (1/TD(y)) = TD(y) / TD(y*) p.logratio += log(oldtd) - log(newtd); diff --git a/src/changestats_dgw_sp.c b/src/changestats_dgw_sp.c deleted file mode 100644 index ab7f18b71..000000000 --- a/src/changestats_dgw_sp.c +++ /dev/null @@ -1,500 +0,0 @@ -/* File src/changestats_dgw_sp.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_dgw_sp.h" -#include "ergm_storage.h" -#include "changestats.h" - -#define all_calcs(term) \ - dvec_calc(term) \ - dist_calc(term) \ - gw_calc(term) - -#define all_calcs2(term) \ - dvec_calc2(term) \ - dist_calc2(term) \ - gw_calc2(term) - - -#define sp_args tail,head,mtp,nwp,edgestate,spcache,N_CHANGE_STATS,dvec,CHANGE_STAT - -#define dvec_calc(term) \ - static inline void term ## _calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, int nd, Vertex *dvec, double *cs) { \ - int echange = edgestate ? -1 : 1; \ - term ## _change({ \ - for(unsigned int j = 0; j < nd; j++){ \ - Vertex deg = dvec[j]; \ - cs[j] += ((L2+echange == deg) - (L2 == deg)); \ - } \ - },{ \ - for(unsigned int j = 0; j < nd; j++){ \ - Vertex deg = dvec[j]; \ - cs[j] += (echange)*(L2 == deg); \ - } \ - }); \ - } - -#define dvec_calc2(term) \ - static inline void term ## _calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, int nd, Vertex *dvec, double *cs) { \ - int echange = edgestate ? -1 : 1; \ - term ## _change({ \ - for(unsigned int j = 0; j < nd; j++){ \ - Vertex deg = (Vertex)dvec[j]; \ - cs[j] += ((L2+echange == deg) - (L2 == deg))*2; \ - } \ - },{ \ - for(unsigned int j = 0; j < nd; j++){ \ - Vertex deg = (Vertex)dvec[j]; \ - cs[j] += (echange)*(L2 == deg); \ - } \ - }); \ - } - - -#define spd_args tail,head,mtp,nwp,edgestate,spcache,N_CHANGE_STATS,CHANGE_STAT - -#define dist_calc(term) \ - static inline void term ## _dist_calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, int nd, double *cs) { \ - int echange = edgestate ? -1 : 1; \ - term ## _change({ \ - int nL2 = L2 + echange; \ - if(nL2 > nd) cutoff_error(mtp); \ - if(L2) cs[L2-1]--; \ - if(nL2) cs[nL2-1]++; \ - },{ \ - if(L2 > nd) cutoff_error(mtp); \ - if(L2) cs[L2-1] += echange; \ - }); \ - } - -#define dist_calc2(term) \ - static inline void term ## _dist_calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, int nd, double *cs) { \ - int echange = edgestate ? -1 : 1; \ - term ## _change({ \ - int nL2 = L2 + echange; \ - if(nL2 > nd) cutoff_error(mtp); \ - if(L2) cs[L2-1]-=2; \ - if(nL2) cs[nL2-1]+=2; \ - },{ \ - if(L2 > nd) cutoff_error(mtp); \ - if(L2) cs[L2-1] += echange; \ - }); \ - } - - -#define gwsp_args tail,head,mtp,nwp,edgestate,spcache,alpha,loneexpa - -#define gw_calc(term) \ - static inline double term ## _gw_calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa) { \ - double cumchange = 0; \ - term ## _change({ \ - cumchange += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; \ - },{ \ - cumchange += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; \ - }); \ - return cumchange; \ - } - - -#define gw_calc2(term) \ - static inline double term ## _gw_calc(Vertex tail, Vertex head, ModelTerm *mtp, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa) { \ - double cumchange = 0; \ - term ## _change({ \ - cumchange += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * 2; \ - },{ \ - cumchange += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; \ - }); \ - return cumchange; \ - } - - -all_calcs(dspUTP) -all_calcs(dspOTP) -all_calcs(dspITP) -all_calcs2(dspOSP) -all_calcs2(dspISP) -all_calcs2(dspRTP) - -/***************** - changestat: d_dsp -*****************/ -/* - Note that d_esp is a meta-function, dispatching actual changescore - calculation to one of the esp*_calc routines, based on the selected shared - partner type code. - - Type codes are as follows (where (i,j) is the focal edge): - - UTP - Undirected two-path (undirected graphs only) - OTP - Outgoing two-path (i->k->j) - ITP - Incoming two-path (i<-k<-j) - RTP - Reciprocated two-path (i<->k<->j) - OSP - Outgoing shared partner (i->k<-j) - ISP - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. UTP should always be used for undirected graphs; OTP is the traditional directed default. -*/ -C_CHANGESTAT_FN(c_ddsp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - Vertex *dvec = (Vertex*) IINPUT_PARAM+1; /*Get the pointer to the L2 stats list*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: dspUTP_calc(sp_args); break; - case L2OTP: dspOTP_calc(sp_args); break; - case L2ITP: dspITP_calc(sp_args); break; - case L2RTP: dspRTP_calc(sp_args); break; - case L2OSP: dspOSP_calc(sp_args); break; - case L2ISP: dspISP_calc(sp_args); break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -C_CHANGESTAT_FN(c_ddspdist) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: dspUTP_dist_calc(spd_args); break; - case L2OTP: dspOTP_dist_calc(spd_args); break; - case L2ITP: dspITP_dist_calc(spd_args); break; - case L2RTP: dspRTP_dist_calc(spd_args); break; - case L2OSP: dspOSP_dist_calc(spd_args); break; - case L2ISP: dspISP_dist_calc(spd_args); break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -/***************** - changestat: d_gwdsp -*****************/ - -/* - Note that d_gwesp is a meta-function for all geometrically weighted DSP stats; the specific type of DSP to be employed is determined by the type argument (INPUT_PARAM[1]). Type codes are as follows (where (i,j) is the focal edge): - - OTP (0) - Outgoing two-path (i->k->j) - ITP (1) - Incoming two-path (i<-k<-j) - RTP (2) - Reciprocated two-path (i<->k<->j) - OSP (3) - Outgoing shared partner (i->k<-j) - ISP (4) - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. The default, OTP, retains the original behavior of esp/gwesp. In the case of undirected graphs, OTP should be used (the others assume a directed network memory structure, and are not safe in the undirected case). -*/ -C_CHANGESTAT_FN(c_dgwdsp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double loneexpa = log1mexp(alpha); /*Precompute log(1-exp(-alpha))*/ - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - double cumchange = 0; - - /*Obtain the DSP changescores (by type)*/ - switch(type){ - case L2UTP: cumchange = dspUTP_gw_calc(gwsp_args); break; - case L2OTP: cumchange = dspOTP_gw_calc(gwsp_args); break; - case L2ITP: cumchange = dspITP_gw_calc(gwsp_args); break; - case L2RTP: cumchange = dspRTP_gw_calc(gwsp_args); break; - case L2OSP: cumchange = dspOSP_gw_calc(gwsp_args); break; - case L2ISP: cumchange = dspISP_gw_calc(gwsp_args); break; - } - - CHANGE_STAT[0] = edgestate ? -cumchange : cumchange; -} - - -all_calcs(espUTP) - all_calcs(espOTP) - all_calcs(espITP) - all_calcs(espOSP) - all_calcs(espISP) - all_calcs(espRTP) - - -/***************** - changestat: d_esp -*****************/ -/* - Note that d_esp is a meta-function, dispatching actual changescore - calculation to one of the esp*_calc routines, based on the selected shared - partner type code. - - Type codes are as follows (where (i,j) is the focal edge): - - UTP - Undirected two-path (undirected graphs only) - OTP - Outgoing two-path (i->k->j) - ITP - Incoming two-path (i<-k<-j) - RTP - Reciprocated two-path (i<->k<->j) - OSP - Outgoing shared partner (i->k<-j) - ISP - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. UTP should always be used for undirected graphs; OTP is the traditional directed default. -*/ - C_CHANGESTAT_FN(c_desp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - Vertex *dvec = (Vertex*) IINPUT_PARAM+1; /*Get the pointer to the ESP stats list*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: espUTP_calc(sp_args); break; - case L2OTP: espOTP_calc(sp_args); break; - case L2ITP: espITP_calc(sp_args); break; - case L2RTP: espRTP_calc(sp_args); break; - case L2OSP: espOSP_calc(sp_args); break; - case L2ISP: espISP_calc(sp_args); break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -C_CHANGESTAT_FN(c_despdist) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: espUTP_dist_calc(spd_args); break; - case L2OTP: espOTP_dist_calc(spd_args); break; - case L2ITP: espITP_dist_calc(spd_args); break; - case L2RTP: espRTP_dist_calc(spd_args); break; - case L2OSP: espOSP_dist_calc(spd_args); break; - case L2ISP: espISP_dist_calc(spd_args); break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -/***************** - changestat: d_gwesp -*****************/ - -/* - Note that d_gwesp is a meta-function for all geometrically weighted ESP stats; the specific type of ESP to be employed is determined by the type argument (INPUT_PARAM[1]). Type codes are as follows (where (i,j) is the focal edge): - - OTP (0) - Outgoing two-path (i->k->j) - ITP (1) - Incoming two-path (i<-k<-j) - RTP (2) - Reciprocated two-path (i<->k<->j) - OSP (3) - Outgoing shared partner (i->k<-j) - ISP (4) - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. The default, OTP, retains the original behavior of esp/gwesp. In the case of undirected graphs, OTP should be used (the others assume a directed network memory structure, and are not safe in the undirected case). -*/ -C_CHANGESTAT_FN(c_dgwesp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double loneexpa = log1mexp(alpha); /*Precompute (1-exp(-alpha))*/ - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - double cumchange = 0; - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: cumchange = espUTP_gw_calc(gwsp_args); break; - case L2OTP: cumchange = espOTP_gw_calc(gwsp_args); break; - case L2ITP: cumchange = espITP_gw_calc(gwsp_args); break; - case L2RTP: cumchange = espRTP_gw_calc(gwsp_args); break; - case L2OSP: cumchange = espOSP_gw_calc(gwsp_args); break; - case L2ISP: cumchange = espISP_gw_calc(gwsp_args); break; - } - - CHANGE_STAT[0] = edgestate ? -cumchange : cumchange; -} - - -/***************** - changestat: d_nsp -*****************/ -/* - Note that d_esp is a meta-function, dispatching actual changescore - calculation to one of the esp*_calc routines, based on the selected shared - partner type code. - - Type codes are as follows (where (i,j) is the focal edge): - - UTP - Undirected two-path (undirected graphs only) - OTP - Outgoing two-path (i->k->j) - ITP - Incoming two-path (i<-k<-j) - RTP - Reciprocated two-path (i<->k<->j) - OSP - Outgoing shared partner (i->k<-j) - ISP - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. UTP should always be used for undirected graphs; OTP is the traditional directed default. -*/ -#define NEGATE_CHANGE_STATS for(unsigned int i = 0; i < N_CHANGE_STATS; i++) CHANGE_STAT[i] *= -1; -C_CHANGESTAT_FN(c_dnsp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - Vertex *dvec = (Vertex*) IINPUT_PARAM+1; /*Get the pointer to the NSP stats list*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: - espUTP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspUTP_calc(sp_args); - break; - case L2OTP: - espOTP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspOTP_calc(sp_args); - break; - case L2ITP: - espITP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspITP_calc(sp_args); - break; - case L2RTP: - espRTP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspRTP_calc(sp_args); - break; - case L2OSP: - espOSP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspOSP_calc(sp_args); - break; - case L2ISP: - espISP_calc(sp_args); - NEGATE_CHANGE_STATS; - dspISP_calc(sp_args); - break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -C_CHANGESTAT_FN(c_dnspdist) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: - espUTP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspUTP_dist_calc(spd_args); - break; - case L2OTP: - espOTP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspOTP_dist_calc(spd_args); - break; - case L2ITP: - espITP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspITP_dist_calc(spd_args); - break; - case L2RTP: - espRTP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspRTP_dist_calc(spd_args); - break; - case L2OSP: - espOSP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspOSP_dist_calc(spd_args); - break; - case L2ISP: - espISP_dist_calc(spd_args); - NEGATE_CHANGE_STATS; - dspISP_dist_calc(spd_args); - break; - } - /*We're done! (Changestats were written in by the calc routine.)*/ -} - - -/***************** - changestat: d_gwnsp -*****************/ - -/* - Note that d_gwesp is a meta-function for all geometrically weighted NSP stats; the specific type of NSP to be employed is determined by the type argument (INPUT_PARAM[1]). Type codes are as follows (where (i,j) is the focal edge): - - OTP (0) - Outgoing two-path (i->k->j) - ITP (1) - Incoming two-path (i<-k<-j) - RTP (2) - Reciprocated two-path (i<->k<->j) - OSP (3) - Outgoing shared partner (i->k<-j) - ISP (4) - Incoming shared partner (i<-k->j) - - Only one type may be specified per esp term. The default, OTP, retains the original behavior of esp/gwesp. In the case of undirected graphs, OTP should be used (the others assume a directed network memory structure, and are not safe in the undirected case). -*/ -C_CHANGESTAT_FN(c_dgwnsp) { - /*Set things up*/ - StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; - double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double loneexpa = log1mexp(alpha); /*Precompute (1-exp(-alpha))*/ - L2Type type = (L2Type) IINPUT_PARAM[0]; /*Get the L2 type code to be used*/ - double cumchange = 0; - - /*Obtain the changescores (by type)*/ - switch(type){ - case L2UTP: - cumchange = dspUTP_gw_calc(gwsp_args) - espUTP_gw_calc(gwsp_args); - break; - case L2OTP: - cumchange = dspOTP_gw_calc(gwsp_args) - espOTP_gw_calc(gwsp_args); - break; - case L2ITP: - cumchange = dspITP_gw_calc(gwsp_args) - espITP_gw_calc(gwsp_args); - break; - case L2RTP: - cumchange = dspRTP_gw_calc(gwsp_args) - espRTP_gw_calc(gwsp_args); - break; - case L2OSP: - cumchange = dspOSP_gw_calc(gwsp_args) - espOSP_gw_calc(gwsp_args); - break; - case L2ISP: - cumchange = dspISP_gw_calc(gwsp_args) - espISP_gw_calc(gwsp_args); - break; - } - - CHANGE_STAT[0] = edgestate ? -cumchange : cumchange; -} - - -/***************** - changestat: c_ddspbwrap -*****************/ -C_CHANGESTAT_FN(c_ddspbwrap) { - c_ddsp(tail, head, mtp, nwp, edgestate); - - // correct for double counting of directed vs. undirected dyads - for(int ind = 0; ind < N_CHANGE_STATS; ind++) CHANGE_STAT[ind] /= 2.0; -} - - -C_CHANGESTAT_FN(c_ddspdistbwrap) { - c_ddspdist(tail, head, mtp, nwp, edgestate); - - // correct for double counting of directed vs. undirected dyads - for(int ind = 0; ind < N_CHANGE_STATS; ind++) CHANGE_STAT[ind] /= 2.0; -} - - -/***************** - changestat: c_dgwdspbwrap -*****************/ -C_CHANGESTAT_FN(c_dgwdspbwrap) { - c_dgwdsp(tail, head, mtp, nwp, edgestate); - - // correct for double counting of directed vs. undirected dyads - CHANGE_STAT[0] /= 2.0; -} - diff --git a/src/changestats_dgw_sp.cpp b/src/changestats_dgw_sp.cpp new file mode 100644 index 000000000..026567438 --- /dev/null +++ b/src/changestats_dgw_sp.cpp @@ -0,0 +1,249 @@ +/* File src/changestats_dgw_sp.cpp 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 + +#include "cpp/ergm_changestat.h" + +#include "changestats_dgw_sp.h" +#include "changestats.h" + +using ergm::ErgmCppModelTerm; + +namespace { + +inline StoreStrictDyadMapUInt* get_spcache(const ErgmCppModelTerm<>& mt){ + return mt.aux_storage.size() ? static_cast(mt.aux_storage[0]) : nullptr; +} + +inline int dsp_path_multiplier(L2Type type){ + switch(type){ + case L2OSP: + case L2ISP: + case L2RTP: + return 2; + default: + return 1; + } +} + +inline void negate_change_stats(ErgmCppModelTerm<>& mt){ + for(double& stat : mt.stat) stat *= -1.0; +} + +template +inline void dsp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ + int echange = edgestate ? -1 : 1; + ergm::sp::dsp_change(type, tail, head, nwp, spcache, + [&](int L2){ + for(std::size_t j = 0; j < mt.stat.size(); ++j){ + int deg = mt.iinput[j+1]; + mt.stat[j] += (((L2 + echange) == deg) - (L2 == deg)) * PathMultiplier; + } + }, + [&](int){}); +} + +inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ + int echange = edgestate ? -1 : 1; + ergm::sp::esp_change(type, tail, head, nwp, spcache, + [&](int L2){ + for(std::size_t j = 0; j < mt.stat.size(); ++j){ + int deg = mt.iinput[j+1]; + mt.stat[j] += ((L2 + echange) == deg) - (L2 == deg); + } + }, + [&](int L2){ + for(std::size_t j = 0; j < mt.stat.size(); ++j){ + int deg = mt.iinput[j+1]; + mt.stat[j] += echange * (L2 == deg); + } + }); +} + +template +inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ + int echange = edgestate ? -1 : 1; + int nd = static_cast(mt.stat.size()); + ergm::sp::dsp_change(type, tail, head, nwp, spcache, + [&](int L2){ + int nL2 = L2 + echange; + if(nL2 > nd) cutoff_error(mtp); + if(L2) mt.stat[L2-1] -= PathMultiplier; + if(nL2) mt.stat[nL2-1] += PathMultiplier; + }, + [&](int){}); +} + +inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ + int echange = edgestate ? -1 : 1; + int nd = static_cast(mt.stat.size()); + ergm::sp::esp_change(type, tail, head, nwp, spcache, + [&](int L2){ + int nL2 = L2 + echange; + if(nL2 > nd) cutoff_error(mtp); + if(L2) mt.stat[L2-1]--; + if(nL2) mt.stat[nL2-1]++; + }, + [&](int L2){ + if(L2 > nd) cutoff_error(mtp); + if(L2) mt.stat[L2-1] += echange; + }); +} + +template +inline double dsp_gw_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ + double cumchange = 0; + ergm::sp::dsp_change(type, tail, head, nwp, spcache, + [&](int L2){ + cumchange += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * PathMultiplier; + }, + [&](int){}); + return cumchange; +} + +inline double esp_gw_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ + double cumchange = 0; + ergm::sp::esp_change(type, tail, head, nwp, spcache, + [&](int L2){ + cumchange += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; + }, + [&](int L2){ + cumchange += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; + }); + return cumchange; +} + +} // namespace + +/***************** + changestat: d_dsp +*****************/ +C_CHANGESTAT_CPP(ddsp, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nwp, edgestate, spcache); + else dsp_vector_change<2>(type, tail, head, mt, nwp, edgestate, spcache); + }) + +C_CHANGESTAT_CPP(ddspdist, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + }) + +/***************** + changestat: d_gwdsp +*****************/ +C_CHANGESTAT_CPP(dgwdsp, { + auto *spcache = get_spcache(mt); + double alpha = mt.dinput[0]; + double loneexpa = log1mexp(alpha); + L2Type type = static_cast(mt.iinput[0]); + + double cumchange = dsp_path_multiplier(type) == 1 + ? dsp_gw_change<1>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa) + : dsp_gw_change<2>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + + mt.stat[0] = edgestate ? -cumchange : cumchange; + }) + +/***************** + changestat: d_esp +*****************/ +C_CHANGESTAT_CPP(desp, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + esp_vector_change(type, tail, head, mt, nwp, edgestate, spcache); + }) + +C_CHANGESTAT_CPP(despdist, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + esp_dist_change(type, tail, head, mtp, mt, nwp, edgestate, spcache); + }) + +/***************** + changestat: d_gwesp +*****************/ +C_CHANGESTAT_CPP(dgwesp, { + auto *spcache = get_spcache(mt); + double alpha = mt.dinput[0]; + double loneexpa = log1mexp(alpha); + L2Type type = static_cast(mt.iinput[0]); + + double cumchange = esp_gw_change(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + mt.stat[0] = edgestate ? -cumchange : cumchange; + }) + +/***************** + changestat: d_nsp +*****************/ +C_CHANGESTAT_CPP(dnsp, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + esp_vector_change(type, tail, head, mt, nwp, edgestate, spcache); + negate_change_stats(mt); + if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nwp, edgestate, spcache); + else dsp_vector_change<2>(type, tail, head, mt, nwp, edgestate, spcache); + }) + +C_CHANGESTAT_CPP(dnspdist, { + auto *spcache = get_spcache(mt); + L2Type type = static_cast(mt.iinput[0]); + + esp_dist_change(type, tail, head, mtp, mt, nwp, edgestate, spcache); + negate_change_stats(mt); + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + }) + +/***************** + changestat: d_gwnsp +*****************/ +C_CHANGESTAT_CPP(dgwnsp, { + auto *spcache = get_spcache(mt); + double alpha = mt.dinput[0]; + double loneexpa = log1mexp(alpha); + L2Type type = static_cast(mt.iinput[0]); + + double dspchange = dsp_path_multiplier(type) == 1 + ? dsp_gw_change<1>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa) + : dsp_gw_change<2>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + double cumchange = dspchange - esp_gw_change(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + + mt.stat[0] = edgestate ? -cumchange : cumchange; + }) + +/***************** + changestat: c_ddspbwrap +*****************/ +C_CHANGESTAT_CPP(ddspbwrap, { + c_ddsp(tail, head, mtp, nwp, edgestate); + for(double &stat : mt.stat) stat /= 2.0; + }) + +C_CHANGESTAT_CPP(ddspdistbwrap, { + c_ddspdist(tail, head, mtp, nwp, edgestate); + for(double &stat : mt.stat) stat /= 2.0; + }) + +/***************** + changestat: c_dgwdspbwrap +*****************/ +C_CHANGESTAT_CPP(dgwdspbwrap, { + c_dgwdsp(tail, head, mtp, nwp, edgestate); + mt.stat[0] /= 2.0; + }) diff --git a/src/changestats_dgw_sp.h b/src/changestats_dgw_sp.h index c7b4f00c4..06c4b14e9 100644 --- a/src/changestats_dgw_sp.h +++ b/src/changestats_dgw_sp.h @@ -612,4 +612,49 @@ typedef enum {L2UTP, L2OTP, L2ITP, L2RTP, L2OSP, L2ISP} L2Type; }); \ call_subroutine_focus(th, subroutine_focus); +#ifdef __cplusplus +namespace ergm { +inline namespace v1 { +namespace sp { + +template +inline void dsp_change(L2Type type, Vertex tail, Vertex head, Network *nwp, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ + switch(type){ + case L2UTP: { dspUTP_change(update_path(L2);, update_focus(L2);); break; } + case L2OTP: { dspOTP_change(update_path(L2);, update_focus(L2);); break; } + case L2ITP: { dspITP_change(update_path(L2);, update_focus(L2);); break; } + case L2RTP: { dspRTP_change(update_path(L2);, update_focus(L2);); break; } + case L2OSP: { dspOSP_change(update_path(L2);, update_focus(L2);); break; } + case L2ISP: { dspISP_change(update_path(L2);, update_focus(L2);); break; } + default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); + } +} + +template +inline void esp_change(L2Type type, Vertex tail, Vertex head, Network *nwp, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ + switch(type){ + case L2UTP: { espUTP_change(update_path(L2);, update_focus(L2);); break; } + case L2OTP: { espOTP_change(update_path(L2);, update_focus(L2);); break; } + case L2ITP: { espITP_change(update_path(L2);, update_focus(L2);); break; } + case L2RTP: { espRTP_change(update_path(L2);, update_focus(L2);); break; } + case L2OSP: { espOSP_change(update_path(L2);, update_focus(L2);); break; } + case L2ISP: { espISP_change(update_path(L2);, update_focus(L2);); break; } + default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); + } +} + +inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ + int echange = edgestate ? -1 : 1; + int delta = 0; + dsp_change(type, tail, head, nwp, spcache, + [&](int L2){ delta += (L2 + echange != 0) - (L2 != 0); }, + [&](int){}); + return delta; +} + +} // namespace sp +} // namespace v1 +} // namespace ergm +#endif + #endif // _CHANGESTATS_DGW_SP_H_ From ac5ae73a0a29f0f33c43da05e4dcf41e91b6985e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 1 Jul 2026 08:29:19 +0000 Subject: [PATCH 3/7] Refactor SP helper templates --- src/MHproposals_triadic.cpp | 2 +- src/changestats_dgw_sp.cpp | 61 +++++----- src/changestats_dgw_sp.h | 224 +++++++++++++++++++++++++++++++++--- 3 files changed, 238 insertions(+), 49 deletions(-) diff --git a/src/MHproposals_triadic.cpp b/src/MHproposals_triadic.cpp index c650d783a..869d5d8e0 100644 --- a/src/MHproposals_triadic.cpp +++ b/src/MHproposals_triadic.cpp @@ -137,7 +137,7 @@ extern "C" MH_P_FN(Mp_SPDyad){ Rboolean edgeflag = nw(p.tail[0], p.head[0]); Vertex tail = p.tail[0], head = p.head[0]; if(type == L2RTP) error("In ergm:Mp_SPDyad(), an unsupported type of triad: %d.", type); - newtd += ergm::sp::dsp_nonzero_change(type, tail, head, nwp, edgeflag, spcache); + newtd += ergm::sp::dsp_nonzero_change(type, tail, head, nw, edgeflag, spcache); // q(y | y*) / q(y* | y) = 1/TD(y*) / (1/TD(y)) = TD(y) / TD(y*) p.logratio += log(oldtd) - log(newtd); diff --git a/src/changestats_dgw_sp.cpp b/src/changestats_dgw_sp.cpp index 026567438..faffde547 100644 --- a/src/changestats_dgw_sp.cpp +++ b/src/changestats_dgw_sp.cpp @@ -15,6 +15,7 @@ #include "changestats.h" using ergm::ErgmCppModelTerm; +using ergm::ErgmCppNetwork; namespace { @@ -38,9 +39,9 @@ inline void negate_change_stats(ErgmCppModelTerm<>& mt){ } template -inline void dsp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void dsp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; - ergm::sp::dsp_change(type, tail, head, nwp, spcache, + ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ for(std::size_t j = 0; j < mt.stat.size(); ++j){ int deg = mt.iinput[j+1]; @@ -50,9 +51,9 @@ inline void dsp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppMode [&](int){}); } -inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; - ergm::sp::esp_change(type, tail, head, nwp, spcache, + ergm::sp::esp_change(type, tail, head, nw, spcache, [&](int L2){ for(std::size_t j = 0; j < mt.stat.size(); ++j){ int deg = mt.iinput[j+1]; @@ -68,10 +69,10 @@ inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppMode } template -inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int nd = static_cast(mt.stat.size()); - ergm::sp::dsp_change(type, tail, head, nwp, spcache, + ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ int nL2 = L2 + echange; if(nL2 > nd) cutoff_error(mtp); @@ -81,10 +82,10 @@ inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mt [&](int){}); } -inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int nd = static_cast(mt.stat.size()); - ergm::sp::esp_change(type, tail, head, nwp, spcache, + ergm::sp::esp_change(type, tail, head, nw, spcache, [&](int L2){ int nL2 = L2 + echange; if(nL2 > nd) cutoff_error(mtp); @@ -98,9 +99,9 @@ inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mt } template -inline double dsp_gw_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ +inline double dsp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ double cumchange = 0; - ergm::sp::dsp_change(type, tail, head, nwp, spcache, + ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ cumchange += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * PathMultiplier; }, @@ -108,9 +109,9 @@ inline double dsp_gw_change(L2Type type, Vertex tail, Vertex head, Network *nwp, return cumchange; } -inline double esp_gw_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ +inline double esp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ double cumchange = 0; - ergm::sp::esp_change(type, tail, head, nwp, spcache, + ergm::sp::esp_change(type, tail, head, nw, spcache, [&](int L2){ cumchange += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; }, @@ -129,16 +130,16 @@ C_CHANGESTAT_CPP(ddsp, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nwp, edgestate, spcache); - else dsp_vector_change<2>(type, tail, head, mt, nwp, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nw, edgestate, spcache); + else dsp_vector_change<2>(type, tail, head, mt, nw, edgestate, spcache); }) C_CHANGESTAT_CPP(ddspdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nwp, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nw, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mtp, mt, nw, edgestate, spcache); }) /***************** @@ -151,8 +152,8 @@ C_CHANGESTAT_CPP(dgwdsp, { L2Type type = static_cast(mt.iinput[0]); double cumchange = dsp_path_multiplier(type) == 1 - ? dsp_gw_change<1>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa) - : dsp_gw_change<2>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + ? dsp_gw_change<1>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa) + : dsp_gw_change<2>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); mt.stat[0] = edgestate ? -cumchange : cumchange; }) @@ -164,14 +165,14 @@ C_CHANGESTAT_CPP(desp, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_vector_change(type, tail, head, mt, nwp, edgestate, spcache); + esp_vector_change(type, tail, head, mt, nw, edgestate, spcache); }) C_CHANGESTAT_CPP(despdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_dist_change(type, tail, head, mtp, mt, nwp, edgestate, spcache); + esp_dist_change(type, tail, head, mtp, mt, nw, edgestate, spcache); }) /***************** @@ -183,7 +184,7 @@ C_CHANGESTAT_CPP(dgwesp, { double loneexpa = log1mexp(alpha); L2Type type = static_cast(mt.iinput[0]); - double cumchange = esp_gw_change(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + double cumchange = esp_gw_change(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); mt.stat[0] = edgestate ? -cumchange : cumchange; }) @@ -194,20 +195,20 @@ C_CHANGESTAT_CPP(dnsp, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_vector_change(type, tail, head, mt, nwp, edgestate, spcache); + esp_vector_change(type, tail, head, mt, nw, edgestate, spcache); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nwp, edgestate, spcache); - else dsp_vector_change<2>(type, tail, head, mt, nwp, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nw, edgestate, spcache); + else dsp_vector_change<2>(type, tail, head, mt, nw, edgestate, spcache); }) C_CHANGESTAT_CPP(dnspdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_dist_change(type, tail, head, mtp, mt, nwp, edgestate, spcache); + esp_dist_change(type, tail, head, mtp, mt, nw, edgestate, spcache); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nwp, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mtp, mt, nwp, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nw, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mtp, mt, nw, edgestate, spcache); }) /***************** @@ -220,9 +221,9 @@ C_CHANGESTAT_CPP(dgwnsp, { L2Type type = static_cast(mt.iinput[0]); double dspchange = dsp_path_multiplier(type) == 1 - ? dsp_gw_change<1>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa) - : dsp_gw_change<2>(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); - double cumchange = dspchange - esp_gw_change(type, tail, head, nwp, edgestate, spcache, alpha, loneexpa); + ? dsp_gw_change<1>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa) + : dsp_gw_change<2>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); + double cumchange = dspchange - esp_gw_change(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); mt.stat[0] = edgestate ? -cumchange : cumchange; }) diff --git a/src/changestats_dgw_sp.h b/src/changestats_dgw_sp.h index 06c4b14e9..6bce80eb3 100644 --- a/src/changestats_dgw_sp.h +++ b/src/changestats_dgw_sp.h @@ -613,40 +613,228 @@ typedef enum {L2UTP, L2OTP, L2ITP, L2RTP, L2OSP, L2ISP} L2Type; call_subroutine_focus(th, subroutine_focus); #ifdef __cplusplus +#include "cpp/ergm_network.h" + namespace ergm { inline namespace v1 { namespace sp { -template -inline void dsp_change(L2Type type, Vertex tail, Vertex head, Network *nwp, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ +template +inline int count_otp(NetworkView& nw, Vertex tail, Vertex head){ + int count = 0; + for(auto u: nw.in_neighbors(head)) count += nw(tail, u); + return count; +} + +template +inline int count_utp(NetworkView& nw, Vertex tail, Vertex head){ + int count = 0; + for(auto u: nw.neighbors(head)) count += nw(u, tail); + return count; +} + +template +inline int count_osp(NetworkView& nw, Vertex tail, Vertex head){ + int count = 0; + for(auto u: nw.out_neighbors(head)) + if(u != tail) count += nw(tail, u); + return count; +} + +template +inline int count_isp(NetworkView& nw, Vertex tail, Vertex head){ + int count = 0; + for(auto u: nw.in_neighbors(head)) + if(u != tail) count += nw(u, tail); + return count; +} + +template +inline int count_rtp(NetworkView& nw, Vertex tail, Vertex head, Vertex exclude1, Vertex exclude2){ + int count = 0; + for(auto u: nw.out_neighbors(tail)) + if(u != exclude1 && u != exclude2 && nw(u, tail)) + count += nw(u, head) && nw(head, u); + return count; +} + +template +inline void dsp_change(Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus){ + if constexpr(type == L2UTP){ + for(auto u: nw.neighbors(head)) + if(u != tail) + update_path(spcache ? GETUDMUI(tail, u, spcache) : count_utp(nw, tail, u)); + + for(auto u: nw.neighbors(tail)) + if(u != head) + update_path(spcache ? GETUDMUI(u, head, spcache) : count_utp(nw, u, head)); + }else if constexpr(type == L2OTP){ + for(auto k: nw.out_neighbors(head)) + if(k != tail) + update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); + + for(auto k: nw.in_neighbors(tail)) + if(k != head) + update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + }else if constexpr(type == L2ITP){ + for(auto k: nw.out_neighbors(head)) + if(k != tail) + update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); + + for(auto k: nw.in_neighbors(tail)) + if(k != head) + update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + }else if constexpr(type == L2RTP){ + if(nw(head, tail)){ + for(auto k: nw.out_neighbors(tail)) + if(k != head && nw(k, tail)) + update_path(spcache ? GETUDMUI(k, head, spcache) : count_rtp(nw, k, head, tail, head)); + + for(auto k: nw.out_neighbors(head)) + if(k != tail && nw(k, head)) + update_path(spcache ? GETUDMUI(k, tail, spcache) : count_rtp(nw, k, tail, head, tail)); + } + }else if constexpr(type == L2OSP){ + for(auto k: nw.in_neighbors(head)) + if(k != tail) + update_path(spcache ? GETUDMUI(tail, k, spcache) : count_osp(nw, tail, k)); + }else if constexpr(type == L2ISP){ + for(auto k: nw.out_neighbors(tail)) + if(k != head) + update_path(spcache ? GETUDMUI(k, head, spcache) : count_isp(nw, head, k)); + } +} + +template +inline void dsp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ switch(type){ - case L2UTP: { dspUTP_change(update_path(L2);, update_focus(L2);); break; } - case L2OTP: { dspOTP_change(update_path(L2);, update_focus(L2);); break; } - case L2ITP: { dspITP_change(update_path(L2);, update_focus(L2);); break; } - case L2RTP: { dspRTP_change(update_path(L2);, update_focus(L2);); break; } - case L2OSP: { dspOSP_change(update_path(L2);, update_focus(L2);); break; } - case L2ISP: { dspISP_change(update_path(L2);, update_focus(L2);); break; } + case L2UTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2OTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2ITP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2RTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2OSP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2ISP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } } -template -inline void esp_change(L2Type type, Vertex tail, Vertex head, Network *nwp, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ +template +inline void esp_change(Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ + if constexpr(type == L2UTP){ + int L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; + + for(auto u: nw.neighbors(head)) + if(nw(u, tail)){ + if(!spcache) L2th++; + update_path(spcache ? GETUDMUI(tail, u, spcache) : count_utp(nw, tail, u)); + update_path(spcache ? GETUDMUI(u, head, spcache) : count_utp(nw, u, head)); + } + + update_focus(L2th); + }else if constexpr(type == L2OTP){ + int L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; + + for(auto k: nw.out_neighbors(tail)){ + if(!spcache && k != head && nw(k, head)) L2th++; + if(k != head && nw(head, k)) + update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); + } + + for(auto k: nw.in_neighbors(head)) + if(k != tail && nw(k, tail)) + update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + + update_focus(L2th); + }else if constexpr(type == L2ITP){ + int L2th = spcache ? GETDDMUI(head, tail, spcache) : 0; + + for(auto k: nw.out_neighbors(head)) + if(k != tail && nw(k, tail)){ + if(!spcache) L2th++; + update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + } + + for(auto k: nw.in_neighbors(tail)) + if(k != head && nw(head, k)) + update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); + + update_focus(L2th); + }else if constexpr(type == L2RTP){ + int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + bool htedge = nw(head, tail); + + for(auto k: nw.in_neighbors(tail)){ + if(k != head){ + if(!spcache) L2th += nw(tail, k) && nw(head, k) && nw(k, head); + if(htedge && nw(head, k) && nw(k, head)) + update_path(spcache ? GETUDMUI(k, tail, spcache) : count_rtp(nw, k, tail, tail, 0)); + } + } + + for(auto k: nw.out_neighbors(tail)) + if(k != head && htedge && nw(head, k) && nw(k, head)) + update_path(spcache ? GETUDMUI(tail, k, spcache) : count_rtp(nw, k, tail, tail, 0)); + + for(auto k: nw.in_neighbors(head)) + if(k != tail && htedge && nw(tail, k) && nw(k, tail)) + update_path(spcache ? GETUDMUI(k, head, spcache) : count_rtp(nw, k, head, head, 0)); + + for(auto k: nw.out_neighbors(head)) + if(k != tail && htedge && nw(tail, k) && nw(k, tail)) + update_path(spcache ? GETUDMUI(head, k, spcache) : count_rtp(nw, k, head, head, 0)); + + update_focus(L2th); + }else if constexpr(type == L2OSP){ + int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + + for(auto k: nw.out_neighbors(tail)) + if(k != head){ + if(!spcache) L2th += nw(head, k); + if(nw(k, head)) + update_path(spcache ? GETUDMUI(tail, k, spcache) : count_osp(nw, tail, k)); + } + + for(auto k: nw.in_neighbors(tail)) + if(k != head && nw(k, head)) + update_path(spcache ? GETUDMUI(k, tail, spcache) : count_osp(nw, tail, k)); + + update_focus(L2th); + }else if constexpr(type == L2ISP){ + int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + + for(auto k: nw.in_neighbors(head)) + if(k != tail){ + if(!spcache) L2th += nw(k, tail); + if(nw(tail, k)) + update_path(spcache ? GETUDMUI(k, head, spcache) : count_isp(nw, head, k)); + } + + for(auto k: nw.out_neighbors(head)) + if(k != tail && nw(tail, k)) + update_path(spcache ? GETUDMUI(head, k, spcache) : count_isp(nw, head, k)); + + update_focus(L2th); + } +} + +template +inline void esp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ switch(type){ - case L2UTP: { espUTP_change(update_path(L2);, update_focus(L2);); break; } - case L2OTP: { espOTP_change(update_path(L2);, update_focus(L2);); break; } - case L2ITP: { espITP_change(update_path(L2);, update_focus(L2);); break; } - case L2RTP: { espRTP_change(update_path(L2);, update_focus(L2);); break; } - case L2OSP: { espOSP_change(update_path(L2);, update_focus(L2);); break; } - case L2ISP: { espISP_change(update_path(L2);, update_focus(L2);); break; } + case L2UTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2OTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2ITP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2RTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2OSP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; + case L2ISP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } } -inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, Network *nwp, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +template +inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int delta = 0; - dsp_change(type, tail, head, nwp, spcache, + dsp_change(type, tail, head, nw, spcache, [&](int L2){ delta += (L2 + echange != 0) - (L2 != 0); }, [&](int){}); return delta; From f13cb756e446ee482d51830139b1b3474fd9d612 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 1 Jul 2026 08:38:25 +0000 Subject: [PATCH 4/7] Remove obsolete SP header macros --- src/changestats_dgw_sp.h | 599 --------------------------------------- 1 file changed, 599 deletions(-) diff --git a/src/changestats_dgw_sp.h b/src/changestats_dgw_sp.h index 6bce80eb3..a4b1652ff 100644 --- a/src/changestats_dgw_sp.h +++ b/src/changestats_dgw_sp.h @@ -16,603 +16,6 @@ typedef enum {L2UTP, L2OTP, L2ITP, L2RTP, L2OSP, L2ISP} L2Type; -#define call_subroutine_path(count, subroutine_path) \ - {int L2 = (L2 ## count); \ - {subroutine_path}} - -#define call_subroutine_focus(count, subroutine_focus) \ - {int L2 = (L2 ## count); \ - {subroutine_focus}} - - -/************************** - dsp Calculation functions -**************************/ - -/* - Changescore for ESPs based on two-paths in undirected graphs i.e. configurations for edge i<->j such that i<->k<->j (where <-> here denotes an undirected edge). - - UTP: - L2th - count t<->k<->h - L2tk - for each t<->k neq h: k<->h, count u such that k<->u<->h - L2hk - for each h<->k neq t: k<->t, count u such that k<->u<->t - - This function will only work properly with undirected graphs, and should only be called in that case. -*/ - -#define dspUTP_change(subroutine_path, subroutine_focus) \ - /* step through edges of head */ \ - EXEC_THROUGH_EDGES(head,e,u, { \ - if (u!=tail){ \ - int L2tu; \ - if(spcache) L2tu = GETUDMUI(tail,u,spcache); \ - else{ \ - L2tu=0; \ - /* step through edges of u */ \ - EXEC_THROUGH_EDGES(u,f,v, { \ - if(IS_UNDIRECTED_EDGE(v,tail)!= 0) L2tu++; \ - }); \ - } \ - call_subroutine_path(tu, subroutine_path); \ - } \ - }); \ - EXEC_THROUGH_EDGES(tail,e,u, { \ - if (u!=head){ \ - int L2uh; \ - if(spcache) L2uh = GETUDMUI(u,head,spcache); \ - else{ \ - L2uh=0; \ - /* step through edges of u */ \ - EXEC_THROUGH_EDGES(u,f,v, { \ - if(IS_UNDIRECTED_EDGE(v,head)!= 0) L2uh++; \ - }); \ - } \ - call_subroutine_path(uh, subroutine_path); \ - } \ - }); - - -/* - Changescore for dsps based on outgoing two-paths, i.e. configurations for non-edge i->j such that i->k->j. - - This function should only be used in the directed case -*/ - -#define dspOTP_change(subroutine_path, subroutine_focus) \ - /* step through outedges of head (i.e., k: t->k)*/ \ - EXEC_THROUGH_OUTEDGES(head, e, k, { \ - if(k!=tail){ /*Only use contingent cases*/ \ - int L2tk; \ - if(spcache) L2tk = GETDDMUI(tail,k,spcache); \ - else{ \ - L2tk=0; \ - /* step through inedges of k, incl. (head,k) itself */ \ - EXEC_THROUGH_INEDGES(k, f, u, { \ - L2tk+=IS_OUTEDGE(tail,u); /*Increment if there is a trans edge*/ \ - }); \ - } \ - call_subroutine_path(tk, subroutine_path); \ - } \ - }); \ - /* step through inedges of tail (i.e., k: k->t)*/ \ - EXEC_THROUGH_INEDGES(tail, e, k, { \ - if (k!=head){ /*Only use contingent cases*/ \ - int L2kh; \ - if(spcache) L2kh = GETDDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /* step through outedges of k , incl. (k,tail) itself */ \ - EXEC_THROUGH_OUTEDGES(k, f, u, { \ - L2kh+=IS_OUTEDGE(u,head); /*Increment if there is a trans edge*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - }); - - -/* - Changescore for DSPs based on incoming two-paths, i.e. configurations for edge i->j such that j->k->i. - IE cyclical shared partners - - ITP: - L2th - count j->k->i - L2hk - for each j->k neq i: k->i, count u such that k->u->j - L2kt - for each k->i neq j: j->k, count u such that i->u->k - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define dspITP_change(subroutine_path, subroutine_focus) \ - /* step through outedges of head (i.e., k: h->k)*/ \ - EXEC_THROUGH_OUTEDGES(head, e, k, { \ - if((k!=tail)){ /*Only use contingent cases*/ \ - int L2kt; \ - /*We have a h->k->t two-path, so add it to our count.*/ \ - if(spcache) L2kt = GETDDMUI(tail,k,spcache); /* spcache is an OTP cache. */ \ - else{ \ - L2kt=0; \ - /*Now, count # u such that k->u->h (so that we know k's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k, f, u, { \ - L2kt+=IS_OUTEDGE(tail,u); /*Increment if there is a cyclic edge*/ \ - }); \ - } \ - call_subroutine_path(kt, subroutine_path); \ - } \ - }); \ - /* step through inedges of tail (i.e., k: k->t)*/ \ - EXEC_THROUGH_INEDGES(tail, e, k, { \ - if((k!=head)){ /*Only use contingent cases*/ \ - int L2hk; \ - if(spcache) L2hk = GETDDMUI(k,head,spcache); \ - else{ \ - L2hk=0; \ - /*Now, count # u such that t->u->k (so that we know k's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k, f, u, { \ - L2hk+=IS_OUTEDGE(u,head); /*Increment if there is a cyclic edge*/ \ - }); \ - } \ - call_subroutine_path(hk, subroutine_path); \ - } \ - }); - - -/* - Changescore for DSPs based on outgoing shared partners, i.e. configurations for edge i->j such that i->k and j->k (with k!=j). - - OSP: - L2th - count t->k, h->k - L2tk - for each t->k neq h: k->h, count u such that t->u, k->u - L2kt - for each k->t neq h: k->h, count u such that t->u, k->u - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define dspOSP_change(subroutine_path, subroutine_focus) \ - /* step through outedges of tail (i.e., k: t->k, k->h, k!=h)*/ \ - EXEC_THROUGH_INEDGES(head, e, k, { \ - if(k!=tail){ \ - int L2tk; \ - /*Do we have a t->k,h->k SP? If so, add it to our count.*/ \ - if(spcache) L2tk = GETUDMUI(tail,k,spcache); \ - else{ \ - L2tk=0; \ - /*Now, count # u such that t->u,k->u (to get t->k's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k, f, u, { \ - if (u != tail) \ - /*Increment if there is an OSP */ \ - L2tk+=(IS_OUTEDGE(tail,u)); \ - }); \ - } \ - call_subroutine_path(tk, subroutine_path); \ - } \ - }); - - -/* - Changescore for ESPs based on incoming shared partners, i.e. configurations for edge i->j such that i->k and j->k (with k!=j). - - ISP: - L2th - count k->t, k->h - L2hk - for each h->k neq t: t->k, count u such that u->h, u->k - L2kh - for each k->h neq t: t->k, count u such that u->h, u->k - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define dspISP_change(subroutine_path, subroutine_focus) \ - /* step through inedges of head (i.e., k: k->h, t->k, k!=t)*/ \ - EXEC_THROUGH_OUTEDGES(tail, e, k, { \ - int L2kh; \ - if(k!=head){ \ - if(spcache) L2kh = GETUDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /*Now, count # u such that u->h,u->k (to get h>k's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k, f, u, { \ - if(u!=head) \ - L2kh+=IS_OUTEDGE(u,head); /*Increment if there is an ISP*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - }); - - -/* - Changescore for DSPs based on reciprocated two-paths, i.e. configurations for edge i->j such that i<->k and j<->k (with k!=j). - - RTP: - L2kh - for each k->h neq t: h->t,k<->t, count u such that k<->u<->h - L2kt - for each k->t neq h: h->t,k<->h, count u such that k<->u<->t - - Thanks to the symmetries involved, this covers all cases. - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define dspRTP_change(subroutine_path, subroutine_focus) \ - int htedge=IS_OUTEDGE(head,tail); /*Is there an h->t (reciprocating) edge?*/ \ - if(htedge){ /* Otherwise, t->h doesn't make a difference. */ \ - /* step through reciprocated outedges of tail (t->k: k!=h,k<-t)*/ \ - EXEC_THROUGH_OUTEDGES(tail,e,k,{ \ - if(k!=head&&IS_OUTEDGE(k,tail)){ \ - int L2kh; \ - if(spcache) L2kh = GETUDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /*Now, count # u such that k<->u<->h (to get k->h's SP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u,{ \ - if(u!=tail&&u!=head&&(IS_OUTEDGE(u,k))) \ - L2kh+=(IS_OUTEDGE(u,head)&&IS_OUTEDGE(head,u)); /*k<->u<->h?*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - }); \ - /* step through reciprocated outedges of tail (t->k: k!=h,k<-t)*/ \ - EXEC_THROUGH_OUTEDGES(head,e,k,{ \ - if(k!=tail&&IS_OUTEDGE(k,head)){ \ - int L2kt; \ - if(spcache) L2kt = GETUDMUI(k,tail,spcache); \ - else{ \ - L2kt=0; \ - /*Now, count # u such that k<->u<->t (to get k->t's SP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u,{ \ - if(u!=head&&u!=tail&&(IS_OUTEDGE(u,k))) \ - L2kt+=(IS_OUTEDGE(u,tail)&&IS_OUTEDGE(tail,u)); /*k<->u<->t?*/ \ - }); \ - } \ - call_subroutine_path(kt, subroutine_path); \ - } \ - }); \ - } - - - -/************************** - ESP Calculation functions -**************************/ - -/* - Changescore for ESPs based on two-paths in undirected graphs i.e. configurations for edge i<->j such that i<->k<->j (where <-> here denotes an undirected edge). - - UTP: - L2th - count t<->k<->h - L2tk - for each t<->k neq h: k<->h, count u such that k<->u<->h - L2hk - for each h<->k neq t: k<->t, count u such that k<->u<->t - - This function will only work properly with undirected graphs, and should only be called in that case. -*/ -#define espUTP_change(subroutine_path, subroutine_focus) \ - int L2th; \ - if(spcache) L2th = GETDDMUI(tail,head,spcache); else L2th=0; \ - /* step through outedges of head */ \ - EXEC_THROUGH_EDGES(head,e,u, { \ - if (IS_UNDIRECTED_EDGE(u,tail) != 0){ \ - int L2tu; \ - int L2uh; \ - if(spcache){ \ - L2tu = GETUDMUI(tail,u,spcache); \ - L2uh = GETUDMUI(u,head,spcache); \ - }else{ \ - L2th++; \ - L2tu=0; \ - L2uh=0; \ - /* step through edges of u */ \ - EXEC_THROUGH_EDGES(u,f,v, { \ - if(IS_UNDIRECTED_EDGE(v,head)!= 0) L2uh++; \ - if(IS_UNDIRECTED_EDGE(v,tail)!= 0) L2tu++; \ - }); \ - } \ - call_subroutine_path(tu, subroutine_path); \ - call_subroutine_path(uh, subroutine_path); \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - - -/* - Changescore for ESPs based on outgoing two-paths, i.e. configurations for edge i->j such that i->k->j. - - OTP: - L2th - count i->k->j - L2tk - for each i->k neq j: j->k, count u such that i->u->k - L2kh - for each k->j neq i: k->i, count u such that k->u->j - - This function should only be used in the directed case, with espUTP being used in the undirected case. -*/ -#define espOTP_change(subroutine_path, subroutine_focus) \ - int L2th; \ - if(spcache) L2th = GETDDMUI(tail,head,spcache); else L2th=0; \ - /* step through outedges of tail (i.e., k: t->k)*/ \ - EXEC_THROUGH_OUTEDGES(tail,e,k, { \ - if(!spcache&&(k!=head)&&(IS_OUTEDGE(k,head))){ \ - /*We have a t->k->h two-path, so add it to our count.*/ \ - L2th++; \ - } \ - int L2tk; \ - if((k!=head)&&(IS_OUTEDGE(head,k))){ /*Only use contingent cases*/ \ - if(spcache) L2tk = GETDDMUI(tail,k,spcache); \ - else{ \ - L2tk=0; \ - /*Now, count # u such that t->u->k (to find t->k's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k,f,u, { \ - if(u!=tail) \ - L2tk+=IS_OUTEDGE(tail,u); /*Increment if there is a trans edge*/ \ - }); \ - } \ - call_subroutine_path(tk, subroutine_path); \ - } \ - }); \ - /* step through inedges of head (i.e., k: k->h)*/ \ - EXEC_THROUGH_INEDGES(head,e,k, { \ - int L2kh; \ - if((k!=tail)&&(IS_OUTEDGE(k,tail))){ /*Only use contingent cases*/ \ - if(spcache) L2kh = GETDDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /*Now, count # u such that k->u->j (to find k->h's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if(u!=head) \ - L2kh+=IS_OUTEDGE(u,head); /*Increment if there is a trans edge*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - - -/* - Changescore for ESPs based on incoming two-paths, i.e. configurations for edge i->j such that j->k->i. - - ITP: - L2th - count j->k->i - L2hk - for each j->k neq i: k->i, count u such that k->u->j - L2kt - for each k->i neq j: j->k, count u such that i->u->k - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define espITP_change(subroutine_path, subroutine_focus) \ - int L2th; \ - if(spcache) L2th = GETDDMUI(head,tail,spcache); else L2th=0; \ - /* step through outedges of head (i.e., k: h->k)*/ \ - EXEC_THROUGH_OUTEDGES(head,e,k, { \ - int L2hk; \ - if((k!=tail)&&(IS_OUTEDGE(k,tail))){ /*Only use contingent cases*/ \ - if(spcache) L2hk = GETDDMUI(k,head,spcache); \ - else{ \ - /*We have a h->k->t two-path, so add it to our count.*/ \ - L2th++; \ - L2hk=0; \ - /*Now, count # u such that k->u->h (so that we know k's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if(u!=head) \ - L2hk+=IS_OUTEDGE(u,head); /*Increment if there is a cyclic edge*/ \ - }); \ - } \ - call_subroutine_path(hk, subroutine_path); \ - } \ - }); \ - /* step through inedges of tail (i.e., k: k->t)*/ \ - EXEC_THROUGH_INEDGES(tail,e,k, { \ - int L2kt; \ - if((k!=head)&&(IS_OUTEDGE(head,k))){ /*Only use contingent cases*/ \ - if(spcache) L2kt = GETDDMUI(tail,k,spcache); \ - else{ \ - L2kt=0; \ - /*Now, count # u such that t->u->k (so that we know k's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k,f,u, { \ - if(u!=tail) \ - L2kt+=IS_OUTEDGE(tail,u); /*Increment if there is a cyclic edge*/ \ - }); \ - } \ - call_subroutine_path(kt, subroutine_path); \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - - -/* - Changescore for ESPs based on outgoing shared partners, i.e. configurations for edge i->j such that i->k and j->k (with k!=j). - - OSP: - L2th - count t->k, h->k - L2tk - for each t->k neq h: k->h, count u such that t->u, k->u - L2kt - for each k->t neq h: k->h, count u such that t->u, k->u - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define espOSP_change(subroutine_path, subroutine_focus) \ - int L2th; \ - if(spcache) L2th = GETUDMUI(tail,head,spcache); else L2th=0; \ - /* step through outedges of tail (i.e., k: t->k, k->h, k!=h)*/ \ - EXEC_THROUGH_OUTEDGES(tail,e,k, { \ - if(k!=head){ \ - if(!spcache) \ - /*Do we have a t->k,h->k SP? If so, add it to our count.*/ \ - L2th+=IS_OUTEDGE(head,k); \ - \ - if(IS_OUTEDGE(k,head)){ /*Only consider stats that could change*/ \ - int L2tk; \ - if(spcache) L2tk = GETUDMUI(tail,k,spcache); \ - else{ \ - L2tk=0; \ - /*Now, count # u such that t->u,k->u (to get t->k's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if(u!=tail) \ - L2tk+=IS_OUTEDGE(tail,u); /*Increment if there is an OSP*/ \ - }); \ - } \ - call_subroutine_path(tk, subroutine_path); \ - } \ - } \ - }); \ - /* step through inedges of tail (i.e., k: k->t, k->h, k!=h)*/ \ - EXEC_THROUGH_INEDGES(tail,e,k, { \ - if((k!=head)&&(IS_OUTEDGE(k,head))){ /*Only stats that could change*/ \ - int L2kt; \ - if(spcache) L2kt = GETUDMUI(k,tail,spcache); \ - else{ \ - L2kt=0; \ - /*Now, count # u such that t->u,k->u (to get k->t's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if(u!=tail) \ - L2kt+=IS_OUTEDGE(tail,u); /*Increment if there is an OSP*/ \ - }); \ - } \ - call_subroutine_path(kt, subroutine_path); \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - - -/* - Changescore for ESPs based on incoming shared partners, i.e. configurations for edge i->j such that i->k and j->k (with k!=j). - - ISP: - L2th - count k->t, k->h - L2hk - for each h->k neq t: t->k, count u such that u->h, u->k - L2kh - for each k->h neq t: t->k, count u such that u->h, u->k - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define espISP_change(subroutine_path, subroutine_focus) \ - int L2th; \ - if(spcache) L2th = GETUDMUI(tail,head,spcache); else L2th=0; \ - /* step through inedges of head (i.e., k: k->h, t->k, k!=t)*/ \ - EXEC_THROUGH_INEDGES(head,e,k, { \ - if(k!=tail){ \ - if(!spcache) \ - /*Do we have a k->t,k->h SP? If so, add it to our count.*/ \ - L2th+=IS_OUTEDGE(k,tail); \ - \ - if(IS_OUTEDGE(tail,k)){ /*Only consider stats that could change*/ \ - int L2kh; \ - if(spcache) L2kh = GETUDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /*Now, count # u such that u->h,u->k (to get h>k's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k,f,u, { \ - if(u!=head) \ - L2kh+=IS_OUTEDGE(u,head); /*Increment if there is an ISP*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - } \ - }); \ - /* step through outedges of head (i.e., k: h->k, t->k, k!=t)*/ \ - EXEC_THROUGH_OUTEDGES(head,e,k, { \ - if((k!=tail)&&(IS_OUTEDGE(tail,k))){ /*Only stats that could change*/ \ - int L2hk; \ - if(spcache) L2hk = GETUDMUI(head,k,spcache); \ - else{ \ - L2hk=0; \ - /*Now, count # u such that u->h,u->k (to get k->h's ESP value)*/ \ - EXEC_THROUGH_INEDGES(k,f,u, { \ - if(u!=head) \ - L2hk+=IS_OUTEDGE(u,head); /*Increment if there is an ISP*/ \ - }); \ - } \ - call_subroutine_path(hk, subroutine_path); \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - - -/* - Changescore for ESPs based on reciprocated two-paths, i.e. configurations for edge i->j such that i<->k and j<->k (with k!=j). - - RTP: - L2th - count t<->k<->h - L2kt - for each k->t neq h: h->t,k<->h, count u such that k<->u<->t - L2tk - for each t->k neq h: h->t,k<->h, count u such that k<->u<->t - L2kh - for each k->h neq t: h->t,k<->t, count u such that k<->u<->h - L2hk - for each h->k neq t: h->t,k<->t, count u such that k<->u<->h - - We assume that this is only called for directed graphs - otherwise, use the baseline espUTP function. -*/ -#define espRTP_change(subroutine_path, subroutine_focus) \ - int L2th; /*Two-path counts for various edges*/ \ - /* NB: RTP is for directed networks, but the focus dyad is undirected. */ \ - if(spcache) L2th = GETUDMUI(tail,head,spcache); else L2th=0; \ - int htedge=IS_OUTEDGE(head,tail); /*Is there an h->t (reciprocating) edge?*/ \ - /* step through inedges of tail (k->t: k!=h,h->t,k<->h)*/ \ - EXEC_THROUGH_INEDGES(tail,e,k, { \ - if(k!=head){ \ - if(!spcache) \ - /*Do we have a t<->k<->h TP? If so, add it to our count.*/ \ - L2th+=(IS_OUTEDGE(tail,k)&&IS_OUTEDGE(head,k)&&IS_OUTEDGE(k,head)); \ - if(htedge&&IS_OUTEDGE(head,k)&&IS_OUTEDGE(k,head)){ /*Only consider stats that could change*/ \ - int L2kt; \ - if(spcache) L2kt = GETUDMUI(k,tail,spcache); \ - else{ \ - L2kt=0; \ - /*Now, count # u such that k<->u<->t (to get (k,t)'s ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if((u!=tail)&&(IS_OUTEDGE(u,k))) \ - L2kt+=(IS_OUTEDGE(u,tail)&&IS_OUTEDGE(tail,u)); /*k<->u<->t?*/ \ - }); \ - } \ - call_subroutine_path(kt, subroutine_path); \ - } \ - } \ - }); \ - /* step through outedges of tail (t->k: k!=h,h->t,k<->h)*/ \ - EXEC_THROUGH_OUTEDGES(tail,e,k, { \ - if(k!=head){ \ - if(htedge&&IS_OUTEDGE(head,k)&&IS_OUTEDGE(k,head)){ /*Only consider stats that could change*/ \ - int L2tk; \ - if(spcache) L2tk = GETUDMUI(tail,k,spcache); \ - else{ \ - L2tk=0; \ - /*Now, count # u such that k<->u<->t (to get (tk)'s ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if((u!=tail)&&(IS_OUTEDGE(u,k))) \ - L2tk+=(IS_OUTEDGE(u,tail)&&IS_OUTEDGE(tail,u)); /*k<->u<->t?*/ \ - }); \ - } \ - call_subroutine_path(tk, subroutine_path); \ - } \ - } \ - }); \ - /* step through inedges of head (k->h: k!=t,h->t,k<->t)*/ \ - EXEC_THROUGH_INEDGES(head,e,k, { \ - if(k!=tail){ \ - if(htedge&&IS_OUTEDGE(tail,k)&&IS_OUTEDGE(k,tail)){ /*Only consider stats that could change*/ \ - int L2kh; \ - if(spcache) L2kh = GETUDMUI(k,head,spcache); \ - else{ \ - L2kh=0; \ - /*Now, count # u such that k<->u<->h (to get k->h's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if((u!=head)&&(IS_OUTEDGE(u,k))) \ - L2kh+=(IS_OUTEDGE(u,head)&&IS_OUTEDGE(head,u)); /*k<->u<->h?*/ \ - }); \ - } \ - call_subroutine_path(kh, subroutine_path); \ - } \ - } \ - }); \ - /* step through outedges of head (h->k: k!=t,h->t,k<->t)*/ \ - EXEC_THROUGH_OUTEDGES(head,e,k, { \ - if(k!=tail){ \ - if(htedge&&IS_OUTEDGE(tail,k)&&IS_OUTEDGE(k,tail)){ /*Only consider stats that could change*/ \ - int L2hk; \ - if(spcache) L2hk = GETUDMUI(head,k,spcache); \ - else{ \ - L2hk=0; \ - /*Now, count # u such that k<->u<->h (to get h->k's ESP value)*/ \ - EXEC_THROUGH_OUTEDGES(k,f,u, { \ - if((u!=head)&&(IS_OUTEDGE(u,k))) \ - L2hk+=(IS_OUTEDGE(u,head)&&IS_OUTEDGE(head,u)); /*k<->u<->h?*/ \ - }); \ - } \ - call_subroutine_path(hk, subroutine_path); \ - } \ - } \ - }); \ - call_subroutine_focus(th, subroutine_focus); - -#ifdef __cplusplus #include "cpp/ergm_network.h" namespace ergm { @@ -844,5 +247,3 @@ inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, NetworkView } // namespace v1 } // namespace ergm #endif - -#endif // _CHANGESTATS_DGW_SP_H_ From fb3d74734df7763adbe0f1494c3b9720f6326712 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Fri, 17 Jul 2026 10:41:28 +1000 Subject: [PATCH 5/7] Removed the no-longer-necessary mtp arguments. --- src/changestats_dgw_sp.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/changestats_dgw_sp.cpp b/src/changestats_dgw_sp.cpp index faffde547..7dd60ab1b 100644 --- a/src/changestats_dgw_sp.cpp +++ b/src/changestats_dgw_sp.cpp @@ -69,31 +69,31 @@ inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppMode } template -inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int nd = static_cast(mt.stat.size()); ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ int nL2 = L2 + echange; - if(nL2 > nd) cutoff_error(mtp); + if(nL2 > nd) cutoff_error(mt.ptr); if(L2) mt.stat[L2-1] -= PathMultiplier; if(nL2) mt.stat[nL2-1] += PathMultiplier; }, [&](int){}); } -inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ModelTerm *mtp, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int nd = static_cast(mt.stat.size()); ergm::sp::esp_change(type, tail, head, nw, spcache, [&](int L2){ int nL2 = L2 + echange; - if(nL2 > nd) cutoff_error(mtp); + if(nL2 > nd) cutoff_error(mt.ptr); if(L2) mt.stat[L2-1]--; if(nL2) mt.stat[nL2-1]++; }, [&](int L2){ - if(L2 > nd) cutoff_error(mtp); + if(L2 > nd) cutoff_error(mt.ptr); if(L2) mt.stat[L2-1] += echange; }); } @@ -138,8 +138,8 @@ C_CHANGESTAT_CPP(ddspdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nw, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mtp, mt, nw, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mt, nw, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mt, nw, edgestate, spcache); }) /***************** @@ -172,7 +172,7 @@ C_CHANGESTAT_CPP(despdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_dist_change(type, tail, head, mtp, mt, nw, edgestate, spcache); + esp_dist_change(type, tail, head, mt, nw, edgestate, spcache); }) /***************** @@ -205,10 +205,10 @@ C_CHANGESTAT_CPP(dnspdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - esp_dist_change(type, tail, head, mtp, mt, nw, edgestate, spcache); + esp_dist_change(type, tail, head, mt, nw, edgestate, spcache); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mtp, mt, nw, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mtp, mt, nw, edgestate, spcache); + if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mt, nw, edgestate, spcache); + else dsp_dist_change<2>(type, tail, head, mt, nw, edgestate, spcache); }) /***************** From 8f4e8f0343602056318f11de9d8ad3a2cdc05529 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Fri, 17 Jul 2026 14:07:20 +1000 Subject: [PATCH 6/7] Made the dgwesp C++ code more consistent for different calculation types and added comments explaining the algorithms. --- src/changestats_dgw_sp.cpp | 36 ++--- src/changestats_dgw_sp.h | 321 +++++++++++++++++++++---------------- 2 files changed, 201 insertions(+), 156 deletions(-) diff --git a/src/changestats_dgw_sp.cpp b/src/changestats_dgw_sp.cpp index 7dd60ab1b..59915b5b9 100644 --- a/src/changestats_dgw_sp.cpp +++ b/src/changestats_dgw_sp.cpp @@ -99,26 +99,22 @@ inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelT } template -inline double dsp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ - double cumchange = 0; +inline void dsp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ - cumchange += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * PathMultiplier; + mt.stat[0] += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * PathMultiplier; }, [&](int){}); - return cumchange; } -inline double esp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ - double cumchange = 0; +inline void esp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ ergm::sp::esp_change(type, tail, head, nw, spcache, [&](int L2){ - cumchange += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; + mt.stat[0] += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; }, [&](int L2){ - cumchange += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; + mt.stat[0] += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; }); - return cumchange; } } // namespace @@ -151,11 +147,10 @@ C_CHANGESTAT_CPP(dgwdsp, { double loneexpa = log1mexp(alpha); L2Type type = static_cast(mt.iinput[0]); - double cumchange = dsp_path_multiplier(type) == 1 - ? dsp_gw_change<1>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa) - : dsp_gw_change<2>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); + if(dsp_path_multiplier(type) == 1) dsp_gw_change<1>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + else dsp_gw_change<2>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); - mt.stat[0] = edgestate ? -cumchange : cumchange; + if(edgestate) negate_change_stats(mt); }) /***************** @@ -184,8 +179,9 @@ C_CHANGESTAT_CPP(dgwesp, { double loneexpa = log1mexp(alpha); L2Type type = static_cast(mt.iinput[0]); - double cumchange = esp_gw_change(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); - mt.stat[0] = edgestate ? -cumchange : cumchange; + esp_gw_change(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + + if(edgestate) negate_change_stats(mt); }) /***************** @@ -220,12 +216,12 @@ C_CHANGESTAT_CPP(dgwnsp, { double loneexpa = log1mexp(alpha); L2Type type = static_cast(mt.iinput[0]); - double dspchange = dsp_path_multiplier(type) == 1 - ? dsp_gw_change<1>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa) - : dsp_gw_change<2>(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); - double cumchange = dspchange - esp_gw_change(type, tail, head, nw, edgestate, spcache, alpha, loneexpa); + esp_gw_change(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + negate_change_stats(mt); + if(dsp_path_multiplier(type) == 1) dsp_gw_change<1>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + else dsp_gw_change<2>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); - mt.stat[0] = edgestate ? -cumchange : cumchange; + if(edgestate) negate_change_stats(mt); }) /***************** diff --git a/src/changestats_dgw_sp.h b/src/changestats_dgw_sp.h index a4b1652ff..f7cd9b104 100644 --- a/src/changestats_dgw_sp.h +++ b/src/changestats_dgw_sp.h @@ -22,94 +22,117 @@ namespace ergm { inline namespace v1 { namespace sp { -template -inline int count_otp(NetworkView& nw, Vertex tail, Vertex head){ +/* The following functions calculate or obtain the number of two-paths + of specified type from tail to head. */ + +inline int count_utp(ErgmCppNetwork& nw, Vertex tail, Vertex head, StoreStrictDyadMapUInt *spcache){ + if(spcache) return GETUDMUI(tail, head, spcache); int count = 0; - for(auto u: nw.in_neighbors(head)) count += nw(tail, u); + // h - v - t + for(auto v: nw.neighbors(head)) count += nw(v, tail); return count; } -template -inline int count_utp(NetworkView& nw, Vertex tail, Vertex head){ +inline int count_otp(ErgmCppNetwork& nw, Vertex tail, Vertex head, StoreStrictDyadMapUInt *spcache){ + if(spcache) return GETDDMUI(tail, head, spcache); int count = 0; - for(auto u: nw.neighbors(head)) count += nw(u, tail); + // t -> v -> h + for(auto v: nw.out_neighbors(tail)) count += nw(v, head); return count; } -template -inline int count_osp(NetworkView& nw, Vertex tail, Vertex head){ +/* itp == otp with tail and head swapped. */ + +inline int count_rtp(ErgmCppNetwork& nw, Vertex tail, Vertex head, Vertex exclude1, StoreStrictDyadMapUInt *spcache){ + if(spcache) return GETUDMUI(tail, head, spcache); int count = 0; - for(auto u: nw.out_neighbors(head)) - if(u != tail) count += nw(tail, u); + // t <-> v <-> h + for(auto v: nw.out_neighbors(tail)) + if(v != exclude1 && v != head && nw(v, tail)) + count += nw(v, head) && nw(head, v); return count; } -template -inline int count_isp(NetworkView& nw, Vertex tail, Vertex head){ +inline int count_osp(ErgmCppNetwork& nw, Vertex tail, Vertex head, StoreStrictDyadMapUInt *spcache){ + if(spcache) return GETUDMUI(tail, head, spcache); int count = 0; - for(auto u: nw.in_neighbors(head)) - if(u != tail) count += nw(u, tail); + // h -> v <- t + for(auto v: nw.out_neighbors(head)) + if(v != tail) count += nw(tail, v); return count; } -template -inline int count_rtp(NetworkView& nw, Vertex tail, Vertex head, Vertex exclude1, Vertex exclude2){ +inline int count_isp(ErgmCppNetwork& nw, Vertex tail, Vertex head, StoreStrictDyadMapUInt *spcache){ + if(spcache) return GETUDMUI(tail, head, spcache); int count = 0; - for(auto u: nw.out_neighbors(tail)) - if(u != exclude1 && u != exclude2 && nw(u, tail)) - count += nw(u, head) && nw(head, u); + // h <- v -> t + for(auto v: nw.in_neighbors(head)) + if(v != tail) count += nw(v, tail); return count; } -template -inline void dsp_change(Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus){ - if constexpr(type == L2UTP){ +template +inline void dsp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus){ + switch(type){ + case L2UTP: + // h - u - #v - t: focus dyad: t - u for(auto u: nw.neighbors(head)) if(u != tail) - update_path(spcache ? GETUDMUI(tail, u, spcache) : count_utp(nw, tail, u)); + update_path(count_utp(nw, tail, u, spcache)); + // t - u - #v - h: focus dyad: h - u for(auto u: nw.neighbors(tail)) if(u != head) - update_path(spcache ? GETUDMUI(u, head, spcache) : count_utp(nw, u, head)); - }else if constexpr(type == L2OTP){ - for(auto k: nw.out_neighbors(head)) - if(k != tail) - update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); - - for(auto k: nw.in_neighbors(tail)) - if(k != head) - update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); - }else if constexpr(type == L2ITP){ - for(auto k: nw.out_neighbors(head)) - if(k != tail) - update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); - - for(auto k: nw.in_neighbors(tail)) - if(k != head) - update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); - }else if constexpr(type == L2RTP){ - if(nw(head, tail)){ - for(auto k: nw.out_neighbors(tail)) - if(k != head && nw(k, tail)) - update_path(spcache ? GETUDMUI(k, head, spcache) : count_rtp(nw, k, head, tail, head)); + update_path(count_utp(nw, u, head, spcache)); + + break; + case L2OTP: + case L2ITP: // Identical without a focus dyad to give direction. + // h -> u <- #v <- t: focus dyad t - u + for(auto u: nw.out_neighbors(head)) + if(u != tail) + update_path(count_otp(nw, tail, u, spcache)); + + // t <- u -> #v -> h: focus dyad h - u + for(auto u: nw.in_neighbors(tail)) + if(u != head) + update_path(count_otp(nw, u, head, spcache)); - for(auto k: nw.out_neighbors(head)) - if(k != tail && nw(k, head)) - update_path(spcache ? GETUDMUI(k, tail, spcache) : count_rtp(nw, k, tail, head, tail)); + break; + case L2RTP: + if(nw(head, tail)){ + // t <-> u <-> #v <-> h: focus dyad h - u + for(auto u: nw.out_neighbors(tail)) + if(u != head && nw(u, tail)) + update_path(count_rtp(nw, u, head, tail, spcache)); + + // h <-> u <-> #v <-> t: focus dyad t - u + for(auto u: nw.out_neighbors(head)) + if(u != tail && nw(u, head)) + update_path(count_rtp(nw, u, tail, head, spcache)); } - }else if constexpr(type == L2OSP){ - for(auto k: nw.in_neighbors(head)) - if(k != tail) - update_path(spcache ? GETUDMUI(tail, k, spcache) : count_osp(nw, tail, k)); - }else if constexpr(type == L2ISP){ - for(auto k: nw.out_neighbors(tail)) - if(k != head) - update_path(spcache ? GETUDMUI(k, head, spcache) : count_isp(nw, head, k)); + + break; + case L2OSP: + // h <- u -> #v <- t: focus dyad t - u + for(auto u: nw.in_neighbors(head)) + if(u != tail) + update_path(count_osp(nw, tail, u, spcache)); + + break; + case L2ISP: + // t -> u <- #v -> h: focus dyad h - u + for(auto u: nw.out_neighbors(tail)) + if(u != head) + update_path(count_isp(nw, head, u, spcache)); + + break; + default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } } -template -inline void dsp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ +template +inline void dsp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ switch(type){ case L2UTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; case L2OTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; @@ -121,107 +144,134 @@ inline void dsp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, S } } -template -inline void esp_change(Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ - if constexpr(type == L2UTP){ - int L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; +template +inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ + int L2th; + bool htedge; + + switch(type){ + case L2UTP: + L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; + // h - u - t: focus dyad t - h + // h - u - t - #v - u: focus dyad t - u + // u - #v - h - u - t: focus dyad h - u for(auto u: nw.neighbors(head)) if(nw(u, tail)){ if(!spcache) L2th++; - update_path(spcache ? GETUDMUI(tail, u, spcache) : count_utp(nw, tail, u)); - update_path(spcache ? GETUDMUI(u, head, spcache) : count_utp(nw, u, head)); + update_path(count_utp(nw, tail, u, spcache)); + update_path(count_utp(nw, u, head, spcache)); } - update_focus(L2th); - }else if constexpr(type == L2OTP){ - int L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; + break; + case L2OTP: + L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; - for(auto k: nw.out_neighbors(tail)){ - if(!spcache && k != head && nw(k, head)) L2th++; - if(k != head && nw(head, k)) - update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); + // t -> u -> h: focus dyad t -> h + // u <- #v <- t -> u <- h: focus dyad t -> u + for(auto u: nw.out_neighbors(tail)){ + if(!spcache && u != head && nw(u, head)) L2th++; + if(u != head && nw(head, u)) + update_path(count_otp(nw, tail, u, spcache)); } - for(auto k: nw.in_neighbors(head)) - if(k != tail && nw(k, tail)) - update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + // u -> #v -> h <- u -> h: focus dyad u -> h + for(auto u: nw.in_neighbors(head)) + if(u != tail && nw(u, tail)) + update_path(count_otp(nw, u, head, spcache)); - update_focus(L2th); - }else if constexpr(type == L2ITP){ - int L2th = spcache ? GETDDMUI(head, tail, spcache) : 0; + break; + case L2ITP: + L2th = spcache ? GETDDMUI(head, tail, spcache) : 0; - for(auto k: nw.out_neighbors(head)) - if(k != tail && nw(k, tail)){ + // h -> u -> t: focus dyad t -> h + // u -> #v -> h -> u -> t: focus dyad h -> u + for(auto u: nw.out_neighbors(head)) + if(u != tail && nw(u, tail)){ if(!spcache) L2th++; - update_path(spcache ? GETDDMUI(k, head, spcache) : count_otp(nw, k, head)); + update_path(count_otp(nw, u, head, spcache)); } - for(auto k: nw.in_neighbors(tail)) - if(k != head && nw(head, k)) - update_path(spcache ? GETDDMUI(tail, k, spcache) : count_otp(nw, tail, k)); - - update_focus(L2th); - }else if constexpr(type == L2RTP){ - int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; - bool htedge = nw(head, tail); - - for(auto k: nw.in_neighbors(tail)){ - if(k != head){ - if(!spcache) L2th += nw(tail, k) && nw(head, k) && nw(k, head); - if(htedge && nw(head, k) && nw(k, head)) - update_path(spcache ? GETUDMUI(k, tail, spcache) : count_rtp(nw, k, tail, tail, 0)); + // u <- #v <- t <- u <- h: focus dyad u -> t + for(auto u: nw.in_neighbors(tail)) + if(u != head && nw(head, u)) + update_path(count_otp(nw, tail, u, spcache)); + + break; + case L2RTP: + L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + htedge = nw(head, tail); + + // t <-> u <-> h: focus dyad t -> h + // u <-> h -> t <- u <-> #v <-> t: focus dyad t <- u + for(auto u: nw.in_neighbors(tail)){ + if(u != head){ + if(!spcache) L2th += nw(tail, u) && nw(head, u) && nw(u, head); + if(htedge && nw(head, u) && nw(u, head)) + update_path(count_rtp(nw, u, tail, tail, spcache)); } } - for(auto k: nw.out_neighbors(tail)) - if(k != head && htedge && nw(head, k) && nw(k, head)) - update_path(spcache ? GETUDMUI(tail, k, spcache) : count_rtp(nw, k, tail, tail, 0)); - - for(auto k: nw.in_neighbors(head)) - if(k != tail && htedge && nw(tail, k) && nw(k, tail)) - update_path(spcache ? GETUDMUI(k, head, spcache) : count_rtp(nw, k, head, head, 0)); - - for(auto k: nw.out_neighbors(head)) - if(k != tail && htedge && nw(tail, k) && nw(k, tail)) - update_path(spcache ? GETUDMUI(head, k, spcache) : count_rtp(nw, k, head, head, 0)); - - update_focus(L2th); - }else if constexpr(type == L2OSP){ - int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; - - for(auto k: nw.out_neighbors(tail)) - if(k != head){ - if(!spcache) L2th += nw(head, k); - if(nw(k, head)) - update_path(spcache ? GETUDMUI(tail, k, spcache) : count_osp(nw, tail, k)); + // u <-> h -> t -> u <-> #v <-> t: focus dyad t -> u + for(auto u: nw.out_neighbors(tail)) + if(u != head && htedge && nw(head, u) && nw(u, head)) + update_path(count_rtp(nw, u, tail, tail, spcache)); + + // u <-> t <- h <- u <-> #v <-> h: focus dyad u -> h + for(auto u: nw.in_neighbors(head)) + if(u != tail && htedge && nw(tail, u) && nw(u, tail)) + update_path(count_rtp(nw, u, head, head, spcache)); + + // u <-> t -> h <- u <-> #v <-> h: focus dyad u <- h + for(auto u: nw.out_neighbors(head)) + if(u != tail && htedge && nw(tail, u) && nw(u, tail)) + update_path(count_rtp(nw, u, head, head, spcache)); + + break; + case L2OSP: + L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + + // t -> u <- h: focus dyad t -> h + // u -> #v <- t -> u -> h: focus dyad t -> u + for(auto u: nw.out_neighbors(tail)) + if(u != head){ + if(!spcache) L2th += nw(head, u); + if(nw(u, head)) + update_path(count_osp(nw, tail, u, spcache)); } - for(auto k: nw.in_neighbors(tail)) - if(k != head && nw(k, head)) - update_path(spcache ? GETUDMUI(k, tail, spcache) : count_osp(nw, tail, k)); - - update_focus(L2th); - }else if constexpr(type == L2ISP){ - int L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; - - for(auto k: nw.in_neighbors(head)) - if(k != tail){ - if(!spcache) L2th += nw(k, tail); - if(nw(tail, k)) - update_path(spcache ? GETUDMUI(k, head, spcache) : count_isp(nw, head, k)); + // u -> #v <- t <- u -> h: focus dyad u -> t + for(auto u: nw.in_neighbors(tail)) + if(u != head && nw(u, head)) + update_path(count_osp(nw, tail, u, spcache)); + + break; + case L2ISP: + L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; + + // h <- u -> t: focus dyad t -> h + // u <- #v -> h <- u <- t: focus dyad u -> h + for(auto u: nw.in_neighbors(head)) + if(u != tail){ + if(!spcache) L2th += nw(u, tail); + if(nw(tail, u)) + update_path(count_isp(nw, head, u, spcache)); } - for(auto k: nw.out_neighbors(head)) - if(k != tail && nw(tail, k)) - update_path(spcache ? GETUDMUI(head, k, spcache) : count_isp(nw, head, k)); + // u <- #v -> h -> u <- t: focus dyad h -> u + for(auto u: nw.out_neighbors(head)) + if(u != tail && nw(tail, u)) + update_path(count_isp(nw, head, u, spcache)); - update_focus(L2th); + break; + default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } + + update_focus(L2th); } -template -inline void esp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ +template +inline void esp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ switch(type){ case L2UTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; case L2OTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; @@ -233,8 +283,7 @@ inline void esp_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, S } } -template -inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, NetworkView& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ +inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int delta = 0; dsp_change(type, tail, head, nw, spcache, From 8a13eed03cd95b53ab6294f3b8a4027806c2f896 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Fri, 17 Jul 2026 14:58:21 +1000 Subject: [PATCH 7/7] Removed some unnecessary templating. --- src/changestats_dgw_sp.cpp | 31 ++++++++------------- src/changestats_dgw_sp.h | 56 +++++++++++--------------------------- 2 files changed, 27 insertions(+), 60 deletions(-) diff --git a/src/changestats_dgw_sp.cpp b/src/changestats_dgw_sp.cpp index 59915b5b9..a605b74c1 100644 --- a/src/changestats_dgw_sp.cpp +++ b/src/changestats_dgw_sp.cpp @@ -23,7 +23,7 @@ inline StoreStrictDyadMapUInt* get_spcache(const ErgmCppModelTerm<>& mt){ return mt.aux_storage.size() ? static_cast(mt.aux_storage[0]) : nullptr; } -inline int dsp_path_multiplier(L2Type type){ +inline constexpr int dsp_path_multiplier(L2Type type){ switch(type){ case L2OSP: case L2ISP: @@ -38,14 +38,13 @@ inline void negate_change_stats(ErgmCppModelTerm<>& mt){ for(double& stat : mt.stat) stat *= -1.0; } -template inline void dsp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ for(std::size_t j = 0; j < mt.stat.size(); ++j){ int deg = mt.iinput[j+1]; - mt.stat[j] += (((L2 + echange) == deg) - (L2 == deg)) * PathMultiplier; + mt.stat[j] += (((L2 + echange) == deg) - (L2 == deg)) * dsp_path_multiplier(type); } }, [&](int){}); @@ -68,7 +67,6 @@ inline void esp_vector_change(L2Type type, Vertex tail, Vertex head, ErgmCppMode }); } -template inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int nd = static_cast(mt.stat.size()); @@ -76,8 +74,8 @@ inline void dsp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelT [&](int L2){ int nL2 = L2 + echange; if(nL2 > nd) cutoff_error(mt.ptr); - if(L2) mt.stat[L2-1] -= PathMultiplier; - if(nL2) mt.stat[nL2-1] += PathMultiplier; + if(L2) mt.stat[L2-1] -= dsp_path_multiplier(type); + if(nL2) mt.stat[nL2-1] += dsp_path_multiplier(type); }, [&](int){}); } @@ -98,11 +96,10 @@ inline void esp_dist_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelT }); } -template inline void dsp_gw_change(L2Type type, Vertex tail, Vertex head, ErgmCppModelTerm<>& mt, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache, double alpha, double loneexpa){ ergm::sp::dsp_change(type, tail, head, nw, spcache, [&](int L2){ - mt.stat[0] += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * PathMultiplier; + mt.stat[0] += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * dsp_path_multiplier(type); }, [&](int){}); } @@ -126,16 +123,14 @@ C_CHANGESTAT_CPP(ddsp, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nw, edgestate, spcache); - else dsp_vector_change<2>(type, tail, head, mt, nw, edgestate, spcache); + dsp_vector_change(type, tail, head, mt, nw, edgestate, spcache); }) C_CHANGESTAT_CPP(ddspdist, { auto *spcache = get_spcache(mt); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mt, nw, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mt, nw, edgestate, spcache); + dsp_dist_change(type, tail, head, mt, nw, edgestate, spcache); }) /***************** @@ -147,8 +142,7 @@ C_CHANGESTAT_CPP(dgwdsp, { double loneexpa = log1mexp(alpha); L2Type type = static_cast(mt.iinput[0]); - if(dsp_path_multiplier(type) == 1) dsp_gw_change<1>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); - else dsp_gw_change<2>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + dsp_gw_change(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); if(edgestate) negate_change_stats(mt); }) @@ -193,8 +187,7 @@ C_CHANGESTAT_CPP(dnsp, { esp_vector_change(type, tail, head, mt, nw, edgestate, spcache); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_vector_change<1>(type, tail, head, mt, nw, edgestate, spcache); - else dsp_vector_change<2>(type, tail, head, mt, nw, edgestate, spcache); + dsp_vector_change(type, tail, head, mt, nw, edgestate, spcache); }) C_CHANGESTAT_CPP(dnspdist, { @@ -203,8 +196,7 @@ C_CHANGESTAT_CPP(dnspdist, { esp_dist_change(type, tail, head, mt, nw, edgestate, spcache); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_dist_change<1>(type, tail, head, mt, nw, edgestate, spcache); - else dsp_dist_change<2>(type, tail, head, mt, nw, edgestate, spcache); + dsp_dist_change(type, tail, head, mt, nw, edgestate, spcache); }) /***************** @@ -218,8 +210,7 @@ C_CHANGESTAT_CPP(dgwnsp, { esp_gw_change(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); negate_change_stats(mt); - if(dsp_path_multiplier(type) == 1) dsp_gw_change<1>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); - else dsp_gw_change<2>(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); + dsp_gw_change(type, tail, head, mt, nw, edgestate, spcache, alpha, loneexpa); if(edgestate) negate_change_stats(mt); }) diff --git a/src/changestats_dgw_sp.h b/src/changestats_dgw_sp.h index f7cd9b104..477145815 100644 --- a/src/changestats_dgw_sp.h +++ b/src/changestats_dgw_sp.h @@ -71,8 +71,9 @@ inline int count_isp(ErgmCppNetwork& nw, Vertex tail, Vertex head, StoreStrictDy return count; } -template -inline void dsp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus){ + +template +inline void dsp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus){ switch(type){ case L2UTP: // h - u - #v - t: focus dyad: t - u @@ -84,8 +85,8 @@ inline void dsp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.neighbors(tail)) if(u != head) update_path(count_utp(nw, u, head, spcache)); - break; + case L2OTP: case L2ITP: // Identical without a focus dyad to give direction. // h -> u <- #v <- t: focus dyad t - u @@ -97,8 +98,8 @@ inline void dsp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.in_neighbors(tail)) if(u != head) update_path(count_otp(nw, u, head, spcache)); - break; + case L2RTP: if(nw(head, tail)){ // t <-> u <-> #v <-> h: focus dyad h - u @@ -111,41 +112,29 @@ inline void dsp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict if(u != tail && nw(u, head)) update_path(count_rtp(nw, u, tail, head, spcache)); } - break; + case L2OSP: // h <- u -> #v <- t: focus dyad t - u for(auto u: nw.in_neighbors(head)) if(u != tail) update_path(count_osp(nw, tail, u, spcache)); - break; + case L2ISP: // t -> u <- #v -> h: focus dyad h - u for(auto u: nw.out_neighbors(tail)) if(u != head) update_path(count_isp(nw, head, u, spcache)); - break; - default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); - } -} -template -inline void dsp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ - switch(type){ - case L2UTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2OTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2ITP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2RTP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2OSP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2ISP: dsp_change(tail, head, nw, spcache, update_path, update_focus); break; default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } } -template -inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ + +template +inline void esp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ int L2th; bool htedge; @@ -162,8 +151,8 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict update_path(count_utp(nw, tail, u, spcache)); update_path(count_utp(nw, u, head, spcache)); } - break; + case L2OTP: L2th = spcache ? GETDDMUI(tail, head, spcache) : 0; @@ -179,8 +168,8 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.in_neighbors(head)) if(u != tail && nw(u, tail)) update_path(count_otp(nw, u, head, spcache)); - break; + case L2ITP: L2th = spcache ? GETDDMUI(head, tail, spcache) : 0; @@ -196,8 +185,8 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.in_neighbors(tail)) if(u != head && nw(head, u)) update_path(count_otp(nw, tail, u, spcache)); - break; + case L2RTP: L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; htedge = nw(head, tail); @@ -226,8 +215,8 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.out_neighbors(head)) if(u != tail && htedge && nw(tail, u) && nw(u, tail)) update_path(count_rtp(nw, u, head, head, spcache)); - break; + case L2OSP: L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; @@ -244,8 +233,8 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.in_neighbors(tail)) if(u != head && nw(u, head)) update_path(count_osp(nw, tail, u, spcache)); - break; + case L2ISP: L2th = spcache ? GETUDMUI(tail, head, spcache) : 0; @@ -262,27 +251,14 @@ inline void esp_change(Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrict for(auto u: nw.out_neighbors(head)) if(u != tail && nw(tail, u)) update_path(count_isp(nw, head, u, spcache)); - break; + default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); } update_focus(L2th); } -template -inline void esp_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, StoreStrictDyadMapUInt *spcache, UpdatePath update_path, UpdateFocus update_focus){ - switch(type){ - case L2UTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2OTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2ITP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2RTP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2OSP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - case L2ISP: esp_change(tail, head, nw, spcache, update_path, update_focus); break; - default: error("In ergm shared partner helper, an unsupported type of triad: %d.", type); - } -} - inline int dsp_nonzero_change(L2Type type, Vertex tail, Vertex head, ErgmCppNetwork& nw, Rboolean edgestate, StoreStrictDyadMapUInt *spcache){ int echange = edgestate ? -1 : 1; int delta = 0;