From 0359bdeef1e5f7d769daa2ce4c8c0fe79d3b4d5a Mon Sep 17 00:00:00 2001
From: "Pavel N. Krivitsky" <p.krivitsky@unsw.edu.au>
Date: Sun, 24 Nov 2024 22:56:01 +1100
Subject: [PATCH] Improved the precision of GW decay calculations in dgw*sp()
 and in R by doing more of them on log scale.

---
 DESCRIPTION              |  4 ++--
 R/InitErgmTerm.R         |  8 ++++----
 R/ergm.utility.R         |  3 +++
 src/changestats_dgw_sp.c | 22 ++++++++++------------
 4 files changed, 19 insertions(+), 18 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index c3fbdaa81..05b6f8e23 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"),
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/R/ergm.utility.R b/R/ergm.utility.R
index cbf16b115..c7a3e4d06 100644
--- a/R/ergm.utility.R
+++ b/R/ergm.utility.R
@@ -275,3 +275,6 @@ is.const.sample <- function(x){
   if(is.mcmc.list(x)) x <- as.matrix(x)
   isTRUE(all.equal(apply(x,2,stats::sd), rep(0,ncol(x)), check.names=FALSE))
 }
+
+## Implementation of the C macro in R.
+log1mexp <- function(x) ifelse(-x > -log(2), log(-expm1(-x)), log1p(-exp(-x)))
diff --git a/src/changestats_dgw_sp.c b/src/changestats_dgw_sp.c
index fd9862d8d..e7cb78d07 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 += exp(loneexpa*(L2-edgestate));                      \
       },{                                                               \
-        if(alpha < 100.0) cumchange += exp(alpha)*(1-pow(oneexpa, L2)); \
-        else cumchange += L2;                                           \
+        cumchange += exp(alpha + log1mexp(-loneexpa*L2));                \
       });                                                               \
     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 += exp(loneexpa*(L2-edgestate))*2;                    \
       },{                                                               \
-        if(alpha < 100.0) cumchange += exp(alpha)*(1-pow(oneexpa, L2)); \
-        else cumchange += L2;                                           \
+        cumchange += exp(alpha + log1mexp(-loneexpa*L2));                \
       });                                                               \
     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;