int_part.hh: Improve numerical stability

src/graph/inference/blockmodel/graph_blockmodel.cc

@@ -136,6 +136,8 @@ simple_degs_t copy_simple_degs(simple_degs_t& degs)
 
 void export_sbm_state();
 
+double spence(double);
+
 void export_blockmodel_state()
 {
     using namespace boost::python;
@@ -192,6 +194,7 @@ void export_blockmodel_state()
     def("log_q_approx", log_q_approx);
     def("log_q_approx_big", log_q_approx_big);
     def("log_q_approx_small", log_q_approx_small);
+    def("spence", spence);
 
     def("positive_w_log_P", positive_w_log_P<size_t>);
     def("signed_w_log_P", signed_w_log_P<size_t>);

src/graph/inference/support/int_part.cc

@@ -89,9 +89,9 @@ double q_rec_memo(int n, int k)
     return res;
 }
 
-double log_q_approx_big(int n, int k)
+double log_q_approx_big(size_t n, size_t k)
 {
-    double C = M_PI * sqrt(2/3.);
+    constexpr double C = M_PI * sqrt(2/3.);
     double S = C * sqrt(n) - log(4 * sqrt(3) * n);
     if (k < n)
     {
@@ -101,7 +101,7 @@ double log_q_approx_big(int n, int k)
     return S;
 }
 
-double log_q_approx_small(int n, int k)
+double log_q_approx_small(size_t n, size_t k)
 {
     return lbinom_fast(n - 1, k - 1) - lgamma_fast(k + 1);
 }
@@ -112,22 +112,23 @@ double get_v(double u, double epsilon=1e-8)
     double delta = 1;
     while (delta > epsilon)
     {
-        double n_v = u * sqrt(-v*v/2 - spence(exp(v)));
+        // spence(exp(v)) = -spence(exp(-v)) - (v*v)/2
+        double n_v = u * sqrt(spence(exp(-v)));
         delta = abs(n_v - v);
         v = n_v;
     }
     return v;
 }
 
-double log_q_approx(int n, int k)
+double log_q_approx(size_t n, size_t k)
 {
     if (k < pow(n, 1/4.))
         return log_q_approx_small(n, k);
     double u = k / sqrt(n);
     double v = get_v(u);
-    double lf = log(v) - log(1 - exp(-v) * (1 + u * u/2)) / 2 - log(2) * 3 / 2
+    double lf = log(v) - log1p(- exp(-v) * (1 + u * u/2)) / 2 - log(2) * 3 / 2
         - log(u) - log(M_PI);
-    double g = 2 * v / u - u * log(1-exp(-v));
+    double g = 2 * v / u - u * log1p(-exp(-v));
     return lf - log(n) + sqrt(n) * g;
 }
 

src/graph/inference/support/int_part.hh

@@ -23,9 +23,9 @@ using namespace boost;
 void init_q_cache(size_t n_max);
 double q_rec(int n, int k);
 double q_rec_memo(int n, int k);
-double log_q_approx(int n, int k);
-double log_q_approx_big(int n, int k);
-double log_q_approx_small(int n, int k);
+double log_q_approx(size_t n, size_t k);
+double log_q_approx_big(size_t n, size_t k);
+double log_q_approx_small(size_t n, size_t k);
 
 extern boost::multi_array<double, 2> __q_cache;
 

src/graph/inference/support/spence.cc

@@ -143,7 +143,7 @@ double spence(double x)
     y = -w * polevl(w, A, 7) / polevl(w, B, 7);
 
     if (flag & 1)
-        y = (M_PI * M_PI) / 6.0 - std::log(x) * std::log(1.0 - x) - y;
+        y = (M_PI * M_PI) / 6.0 - std::log(x) * std::log1p(-x) - y;
 
     if (flag & 2) {
         z = std::log(x);