Skip to content

Commit

Permalink
Improved the precision of GW decay calculations in *gw*() and in R by…
Browse files Browse the repository at this point in the history 
… doing more of them on log scale.
  • Loading branch information
@krivit
krivit committed Nov 28, 2024
1 parent baa6f5b commit b455e74
Show file tree
Hide file tree
Showing 7 changed files with 59 additions and 61 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -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="[email protected]"),
Expand All @@ -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),
Expand Down
8 changes: 4 additions & 4 deletions R/InitErgmTerm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
)
Expand Down
38 changes: 19 additions & 19 deletions src/changestats.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -1904,40 +1904,40 @@ 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);

}

/*****************
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;
}

Expand All @@ -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;
}

Expand All @@ -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 ***********/
Expand Down
22 changes: 10 additions & 12 deletions src/changestats_dgw_sp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; \
}
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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;

Expand Down
42 changes: 21 additions & 21 deletions src/changestats_experimental.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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);
}
Expand Down Expand Up @@ -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) {
Expand All @@ -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);
}
Expand All @@ -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);
}
Expand Down Expand Up @@ -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;
Expand All @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down
Loading

0 comments on commit b455e74

Please sign in to comment.