From b455e746be2a37a2baff0a47a671fa92a089e914 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Sun, 24 Nov 2024 22:56:01 +1100 Subject: [PATCH] Improved the precision of GW decay calculations in *gw*() and in R by doing more of them on log scale. --- DESCRIPTION | 6 ++--- R/InitErgmTerm.R | 8 +++---- src/changestats.c | 38 +++++++++++++++--------------- src/changestats_dgw_sp.c | 22 ++++++++---------- src/changestats_experimental.c | 42 +++++++++++++++++----------------- tests/testthat/test-miss.CD.R | 2 +- tests/testthat/test-miss.R | 2 +- 7 files changed, 59 insertions(+), 61 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index c3fbdaa81..a89fa7bfb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ergm -Version: 4.7.5-7458 -Date: 2024-11-04 +Version: 4.7.5-7485 +Date: 2024-11-24 Title: Fit, Simulate and Diagnose Exponential-Family Models for Networks Authors@R: c( person(c("Mark", "S."), "Handcock", role=c("aut"), email="handcock@stat.ucla.edu"), @@ -27,7 +27,7 @@ Imports: coda (>= 0.19-4.1), trust (>= 0.1-8), lpSolveAPI (>= 5.5.2.0-17.12), - statnet.common (>= 4.10.0), + statnet.common (>= 4.11.0), rle (>= 0.9.2), purrr (>= 1.0.2), rlang (>= 1.1.4), diff --git a/R/InitErgmTerm.R b/R/InitErgmTerm.R index 32ca87dea..3e1542b21 100644 --- a/R/InitErgmTerm.R +++ b/R/InitErgmTerm.R @@ -72,13 +72,13 @@ ergm_GWDECAY <- list( map = function(x,n,...) { i <- 1:n - x[1] * ifelse(i==1, 1, (exp(x[2])*(1-(1-exp(-x[2]))^i))) + x[1] * exp(x[2] + log1mexp(-log1mexp(x[2])*i)) }, gradient = function(x,n,...) { i <- 1:n - e2 <- exp(x[2]) - a <- 1-exp(-x[2]) - rbind((1-a^i)*e2, ifelse(i==1, 0, x[1] * ( (1-a^i)*e2 - i*a^(i-1) ) ) ) + a <- log1mexp(x[2]) + w <- exp(x[2] + log1mexp(-a*i)) + rbind(w, x[1] * (w - i*exp(a*(i-1)))) }, minpar = c(-Inf, 0) ) diff --git a/src/changestats.c b/src/changestats.c index 3af49d849..1dd1e9317 100644 --- a/src/changestats.c +++ b/src/changestats.c @@ -1875,20 +1875,20 @@ S_CHANGESTAT_FN(s_edges) { *****************/ C_CHANGESTAT_FN(c_gwdegree) { int echange=0; - double decay, oneexpd, change; + double decay, loneexpd, change; Vertex taild, headd=0, *id, *od; id=IN_DEG; od=OUT_DEG; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); /* *** don't forget tail -> head */ change = 0.0; echange = edgestate ? -1:+1; taild = od[tail] + id[tail] + (echange - 1)/2; headd = od[head] + id[head] + (echange - 1)/2; - change += echange*(pow(oneexpd,(double)taild)+pow(oneexpd,(double)headd)); + change += echange*(exp(loneexpd*taild)+exp(loneexpd*headd)); CHANGE_STAT[0] = change; @@ -1904,23 +1904,23 @@ C_CHANGESTAT_FN(c_gwdegree_by_attr) { from 1 through N_CHANGE_STATS */ int tailattr, headattr, echange=0; - double decay, oneexpd; + double decay, loneexpd; Vertex taild, headd=0, *id, *od; id=IN_DEG; od=OUT_DEG; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); /* *** don't forget tail -> head */ echange = edgestate ? -1:+1; taild = od[tail] + id[tail] + (echange - 1)/2; tailattr = INPUT_PARAM[tail]; - CHANGE_STAT[tailattr-1] += echange*(pow(oneexpd,(double)taild)); + CHANGE_STAT[tailattr-1] += echange*exp(loneexpd*taild); headd = od[head] + id[head] + (echange - 1)/2; headattr = INPUT_PARAM[head]; - CHANGE_STAT[headattr-1] += echange*(pow(oneexpd,(double)headd)); + CHANGE_STAT[headattr-1] += echange*exp(loneexpd*headd); } @@ -1928,16 +1928,16 @@ C_CHANGESTAT_FN(c_gwdegree_by_attr) { changestat: d_gwidegree *****************/ C_CHANGESTAT_FN(c_gwidegree) { - double decay, oneexpd, change; + double decay, loneexpd, change; Vertex headd=0; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); change = 0.0; /* *** don't forget tail -> head */ headd = IN_DEG[head] - edgestate; - change += (edgestate? -1.0 : 1.0) * pow(oneexpd,(double)headd); + change += (edgestate? -1.0 : 1.0) * exp(loneexpd*headd); CHANGE_STAT[0]=change; } @@ -1951,33 +1951,33 @@ C_CHANGESTAT_FN(c_gwidegree_by_attr) { from 1 through N_CHANGE_STATS */ int headattr, echange; - double decay, oneexpd; + double decay, loneexpd; Vertex headd; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); /* *** don't forget tail -> head */ echange = edgestate ? -1 : 1; headd = IN_DEG[head] + (echange - 1)/2; headattr = INPUT_PARAM[head - BIPARTITE]; /* BIPARTITE to make the b2 version a special case. */ - CHANGE_STAT[headattr-1] += echange*(pow(oneexpd,(double)headd)); + CHANGE_STAT[headattr-1] += echange*exp(loneexpd*headd); } /***************** changestat: d_gwodegree *****************/ C_CHANGESTAT_FN(c_gwodegree) { - double decay, oneexpd, change; + double decay, loneexpd, change; Vertex taild; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); change = 0.0; /* *** don't forget tail -> head */ taild = OUT_DEG[tail] - edgestate; - change += (edgestate? -1 : 1) * pow(oneexpd,(double)taild); + change += (edgestate? -1 : 1) * exp(loneexpd*taild); CHANGE_STAT[0] = change; } @@ -1991,17 +1991,17 @@ C_CHANGESTAT_FN(c_gwodegree_by_attr) { from 1 through N_CHANGE_STATS */ int tailattr, echange; - double decay, oneexpd; + double decay, loneexpd; Vertex taild; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); /* *** don't forget tail -> head */ echange = edgestate ? -1 : 1; taild = OUT_DEG[tail] + (echange - 1)/2; tailattr = INPUT_PARAM[tail]; - CHANGE_STAT[tailattr-1] += echange*(pow(oneexpd,(double)taild)); + CHANGE_STAT[tailattr-1] += echange*exp(loneexpd*taild); } /******************** changestats: H ***********/ diff --git a/src/changestats_dgw_sp.c b/src/changestats_dgw_sp.c index fd9862d8d..406f99886 100644 --- a/src/changestats_dgw_sp.c +++ b/src/changestats_dgw_sp.c @@ -88,29 +88,27 @@ } -#define gwsp_args tail,head,mtp,nwp,edgestate,spcache,alpha,oneexpa +#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 oneexpa) { \ + 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 += pow(oneexpa, L2-edgestate); \ + cumchange += alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0; \ },{ \ - if(alpha < 100.0) cumchange += exp(alpha)*(1-pow(oneexpa, L2)); \ - else cumchange += L2; \ + 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 oneexpa) { \ + 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 += pow(oneexpa, L2-edgestate)*2; \ + cumchange += (alpha ? exp(loneexpa*(L2-edgestate)) : L2-edgestate == 0) * 2; \ },{ \ - if(alpha < 100.0) cumchange += exp(alpha)*(1-pow(oneexpa, L2)); \ - else cumchange += L2; \ + cumchange += alpha ? exp(alpha + log1mexp(-loneexpa*L2)) : L2 != 0; \ }); \ return cumchange; \ } @@ -198,7 +196,7 @@ C_CHANGESTAT_FN(c_dgwdsp) { /*Set things up*/ StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double oneexpa = 1.0-exp(-alpha); /*Precompute (1-exp(-alpha))*/ + double loneexpa = log1mexp(alpha); /*Precompute log(1-exp(-alpha))*/ int type = IINPUT_PARAM[0]; /*Get the ESP type code to be used*/ double cumchange = 0; @@ -299,7 +297,7 @@ C_CHANGESTAT_FN(c_dgwesp) { /*Set things up*/ StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double oneexpa = 1.0-exp(-alpha); /*Precompute (1-exp(-alpha))*/ + double loneexpa = log1mexp(alpha); /*Precompute (1-exp(-alpha))*/ int type = IINPUT_PARAM[0]; /*Get the ESP type code to be used*/ double cumchange = 0; @@ -441,7 +439,7 @@ C_CHANGESTAT_FN(c_dgwnsp) { /*Set things up*/ StoreStrictDyadMapUInt *spcache = N_AUX ? AUX_STORAGE : NULL; double alpha = INPUT_PARAM[0]; /*Get alpha*/ - double oneexpa = 1.0-exp(-alpha); /*Precompute (1-exp(-alpha))*/ + double loneexpa = log1mexp(alpha); /*Precompute (1-exp(-alpha))*/ int type = IINPUT_PARAM[0]; /*Get the ESP type code to be used*/ double cumchange = 0; diff --git a/src/changestats_experimental.c b/src/changestats_experimental.c index 71d2a3b8a..c05e76a3d 100644 --- a/src/changestats_experimental.c +++ b/src/changestats_experimental.c @@ -624,7 +624,7 @@ D_CHANGESTAT_FN(d_geospartner) { *****************/ D_CHANGESTAT_FN(d_gwb1) { int i, echange=0; - double alpha, oneexpa, change; + double alpha, loneexpa, change; Vertex tail, head, taild=0, *od; /* int nb1 , nb2 @@ -635,14 +635,14 @@ D_CHANGESTAT_FN(d_gwb1) { od=OUT_DEG; change = 0.0; alpha = INPUT_PARAM[1]; - oneexpa = 1.0-exp(-alpha); + loneexpa = log1mexp(alpha); CHANGE_STAT[0] = 0.0; FOR_EACH_TOGGLE(i) { echange = IS_OUTEDGE(tail=TAIL(i), head=HEAD(i)) ? -1 : 1; taild = od[tail] + (echange - 1)/2; if(taild!=0){ - change += echange*(1.0-pow(oneexpa,(double)taild)); + change += echange*exp(log1mexp(-loneexpa*taild)); } TOGGLE_IF_MORE_TO_COME(i); } @@ -679,13 +679,13 @@ D_CHANGESTAT_FN(d_gwd) { D_CHANGESTAT_FN(d_gwdegree706) { /* Slight modification to the parameterization in d_gwdegree */ int i, echange=0; - double decay, oneexpd, change; + double decay, loneexpd, change; Vertex tail, head, taild, headd=0, *id, *od; id=IN_DEG; od=OUT_DEG; decay = INPUT_PARAM[0]; - oneexpd = 1.0-exp(-decay); + loneexpd = log1mexp(decay); change = 0.0; FOR_EACH_TOGGLE(i) { @@ -694,10 +694,10 @@ D_CHANGESTAT_FN(d_gwdegree706) { taild = od[tail] + id[tail] + (echange - 1)/2; headd = od[head] + id[head] + (echange - 1)/2; if(taild!=0){ - change -= echange*(1.0-pow(oneexpd,(double)taild)); + change -= echange*exp(log1mexp(-loneexpd*taild)); } if(headd!=0){ - change -= echange*(1.0-pow(oneexpd,(double)headd)); + change -= echange*exp(log1mexp(-loneexpd*headd)); } TOGGLE_IF_MORE_TO_COME(i); } @@ -710,24 +710,24 @@ D_CHANGESTAT_FN(d_gwdegree706) { *****************/ D_CHANGESTAT_FN(d_gwdegreealpha) { int i, echange=0; - double alpha, oneexpa, change; + double alpha, loneexpa, change; Vertex tail, head, taild, headd=0, *id, *od; id=IN_DEG; od=OUT_DEG; change = 0.0; alpha = INPUT_PARAM[0]; - oneexpa = 1.0-exp(-alpha); + loneexpa = log1mexp(alpha); FOR_EACH_TOGGLE(i) { echange = IS_OUTEDGE(tail=TAIL(i), head=HEAD(i)) ? -1 : 1; taild = od[tail] + id[tail] + (echange - 1)/2; headd = od[head] + id[head] + (echange - 1)/2; if(taild!=0){ - change += echange*(1.0-pow(oneexpa,(double)taild)); + change += echange*exp(log1mexp(-loneexpa*taild)); } if(headd!=0){ - change += echange*(1.0-pow(oneexpa,(double)headd)); + change += echange*exp(log1mexp(-loneexpa*headd)); } TOGGLE_IF_MORE_TO_COME(i); } @@ -772,7 +772,7 @@ D_CHANGESTAT_FN(d_gwdegreelambda) { *****************/ D_CHANGESTAT_FN(d_gwb2){ int i, echange=0; - double alpha, oneexpa, change; + double alpha, loneexpa, change; Vertex tail, head, headd=0, *id; /* int nb1, nb2; @@ -783,17 +783,17 @@ D_CHANGESTAT_FN(d_gwb2){ //od=OUT_DEG; change = 0.0; alpha = INPUT_PARAM[1]; - oneexpa = 1.0-exp(-alpha); + loneexpa = log1mexp(alpha); FOR_EACH_TOGGLE(i) { echange = IS_OUTEDGE(tail=TAIL(i), head=HEAD(i)) ? -1 : 1; // taild = od[tail] + id[head] + (echange - 1)/2; headd = id[head] + (echange - 1)/2; // if(taild!=0){ - // change += echange*(1.0-pow(oneexpa,(double)taild)); + // change += echange*exp(log1mexp(-loneexpa*taild)); // } if(headd!=0){ - change += echange*(1.0-pow(oneexpa,(double)headd)); + change += echange*exp(log1mexp(-loneexpa*headd)); } //Rprintf("tail %d head %d headd %d echange %d change %f\n", tail, head, headd, echange,change); // Rprintf(" tail %d head %d taild %d echange %d change %f\n", tail, head, taild, echange,change); @@ -1824,11 +1824,11 @@ D_CHANGESTAT_FN(d_gwb2share) { int i, echange, ochange; int L2uh; Vertex tail, head, u, v; - double alpha, oneexpa, cumchange; + double alpha, loneexpa, cumchange; CHANGE_STAT[0] = 0.0; alpha = INPUT_PARAM[0]; - oneexpa = 1.0-exp(-alpha); + loneexpa = log1mexp(alpha); FOR_EACH_TOGGLE(i) { cumchange=0.0; @@ -1840,7 +1840,7 @@ D_CHANGESTAT_FN(d_gwb2share) { STEP_THROUGH_INEDGES(u, v, f) { if(IS_OUTEDGE(MIN(v,head),MAX(v,head))) L2uh++; } - cumchange += pow(oneexpa,(double)L2uh); + cumchange += exp(loneexpa*L2uh); } } cumchange = echange*cumchange; @@ -1858,11 +1858,11 @@ D_CHANGESTAT_FN(d_gwb1share) { int i, echange, ochange; int L2tu; Vertex tail, head, u, v; - double alpha, oneexpa, cumchange; + double alpha, loneexpa, cumchange; CHANGE_STAT[0] = 0.0; alpha = INPUT_PARAM[0]; - oneexpa = 1.0-exp(-alpha); + loneexpa = log1mexp(alpha); FOR_EACH_TOGGLE(i) { cumchange=0.0; @@ -1874,7 +1874,7 @@ D_CHANGESTAT_FN(d_gwb1share) { STEP_THROUGH_OUTEDGES(u, v, f) { if(IS_OUTEDGE(MIN(v,tail),MAX(v,tail))) L2tu++; } - cumchange += pow(oneexpa,(double)L2tu); + cumchange += exp(loneexpa*L2tu); } } cumchange = echange*cumchange; diff --git a/tests/testthat/test-miss.CD.R b/tests/testthat/test-miss.CD.R index 28f1ccaa6..9d18e4336 100644 --- a/tests/testthat/test-miss.CD.R +++ b/tests/testthat/test-miss.CD.R @@ -59,7 +59,7 @@ test_that("bipartite undirected network", { # Add the curved+missing test here for now test_that("curved+missing", { set.seed(321) - n <- 50 + n <- 30 y <- network.initialize(n, directed=FALSE) # Create an empty network y <- simulate(y~edges, coef=logit(0.12), control=control.simulate(MCMC.burnin=2*n^2)) y.miss <- simulate(y~edges, coef=logit(0.01)) diff --git a/tests/testthat/test-miss.R b/tests/testthat/test-miss.R index 01b195cdd..0e20f707e 100644 --- a/tests/testthat/test-miss.R +++ b/tests/testthat/test-miss.R @@ -75,7 +75,7 @@ test_that("Bipartite Undirected Network", { # Add the curved+missing test here for now test_that("curved+missing", { set.seed(321) - n <- 100 + n <- 30 y <- network.initialize(n, directed=FALSE) # Create an empty network y <- simulate(y~edges, coef=logit(0.12), control=control.simulate(MCMC.burnin=2*n^2)) y.miss <- simulate(y~edges, coef=logit(0.01))