Update statistics() in MonteCarlo and MCMC

examples/benchmark_vegas.py

@@ -1,7 +1,6 @@
 import vegas
 import numpy as np
 import gvar
-import torch
 
 dim = 4
 nitn = 10
@@ -14,14 +13,11 @@ def f_batch(x):
     for d in range(dim):
         dx2 += (x[:, d] - 0.5) ** 2
     return np.exp(-200 * dx2)
-    # ans = np.empty((x.shape[0], 3), float)
-    # dx2 = 0.0
-    # for d in range(dim):
-    #     dx2 += (x[:, d] - 0.5) ** 2
-    # ans[:, 0] = np.exp(-200 * dx2)
-    # ans[:, 1] = x[:, 0] * ans[:, 0]
-    # ans[:, 2] = x[:, 0] ** 2 * ans[:, 0]
-    # return ans
+
+
+@vegas.lbatchintegrand
+def f1_batch(x):
+    return np.log(x[:, 0]) / np.sqrt(x[:, 0])
 
 
 def smc(f, map, neval, dim):
@@ -41,26 +37,53 @@ def mc(f, neval, dim):
     return (np.average(fy), np.std(fy) / neval**0.5)
 
 
-m = vegas.AdaptiveMap(dim * [[0, 1]], ninc=ninc)
-ny = 20000
-# torch.manual_seed(0)
-# y = torch.rand((ny, dim), dtype=torch.float64).numpy()
-y = np.random.uniform(0.0, 1.0, (ny, dim))  # 1000 random y's
+def vegas_map(f, dim, ny, alpha=0.5, niter=10):
+    m = vegas.AdaptiveMap(dim * [[0, 1]], ninc=ninc)
+    y = np.random.uniform(0.0, 1.0, (ny, dim))  # 1000 random y's
+
+    x = np.empty(y.shape, float)  # work space
+    jac = np.empty(y.shape[0], float)
+    f2 = np.empty(y.shape[0], float)
+
+    for itn in range(niter):  # 5 iterations to adapt
+        m.map(y, x, jac)  # compute x's and jac
+
+        f2 = (jac * f(x)) ** 2
+        m.add_training_data(y, f2)  # adapt
+        m.adapt(alpha=alpha)
+    return m
+
 
-x = np.empty(y.shape, float)  # work space
-jac = np.empty(y.shape[0], float)
-f2 = np.empty(y.shape[0], float)
+alpha = 2.0
+ny = 100000
+m = vegas_map(f1_batch, 1, ny, alpha=alpha)
+
+neval = 1000000
+# with map
+r = smc(f1_batch, m, neval, dim)
+print("   SMC + map:", f"{r[0]} +- {r[1]}")
+
+# without map
+r = mc(f1_batch, neval, dim)
+print("SMC (no map):", f"{r[0]} +- {r[1]}")
+
+# vegas with adaptive stratified sampling
+print("VEGAS using adaptive stratified sampling")
+integ = vegas.Integrator(dim * [[0, 1]])
+training = integ(f1_batch, nitn=10, neval=ny, alpha=alpha)  # adapt grid
+# print(integ.map)
+r = smc(f1_batch, integ.map, neval, dim)
+print("   SMC + map:", f"{r[0]} +- {r[1]}")
 
-for itn in range(10):  # 5 iterations to adapt
-    m.map(y, x, jac)  # compute x's and jac
+result = integ(f1_batch, nitn=10, neval=ny, adapt=False)
+# print(result)
+print("       vegas:", f"{result.mean} +- {result.sdev}")
 
-    f2 = (jac * f_batch(x)) ** 2
-    m.add_training_data(y, f2)  # adapt
-    # if itn == 0:
-    #     print(np.array(memoryview(m.sum_f)))
-    #     print(np.array(memoryview(m.n_f)))
-    m.adapt(alpha=0.5)
 
+### benchmark exp(-200 * (x1^2 + x2^2 + x3^2 + x4^2)) integral
+print("\nCalculate the integral h[X] in the bounds [0, 1]^4")
+print("h(X) = exp(-200 * (x1^2 + x2^2 + x3^2 + x4^2))")
+m = vegas_map(f_batch, dim, ny)
 
 # with map
 r = smc(f_batch, m, 50_000, dim)

src/base.py

@@ -1,6 +1,11 @@
 import torch
 from torch import nn
 import numpy as np
+import sys
+
+MINVAL = 10 ** (sys.float_info.min_10_exp + 50)
+MAXVAL = 10 ** (sys.float_info.max_10_exp - 50)
+EPSILON = 1e-16  # Small value to ensure numerical stability
 
 
 class BaseDistribution(nn.Module):