Merge pull request #34 from numericalEFT/pchou

add 3 examples

Author: kunyuan

examples/example_1.py

@@ -0,0 +1,70 @@
+# Example 1: Unit Circle and Half-Sphere Integrands Comparison
+
+import torch
+from MCintegration import MonteCarlo, MarkovChainMonteCarlo, Vegas, set_seed, get_device
+
+set_seed(42)
+device = get_device()
+
+
+def unit_circle_integrand(x, f):
+    f[:, 0] = (x[:, 0] ** 2 + x[:, 1] ** 2 < 1).double()
+    return f[:, 0]
+
+
+def half_sphere_integrand(x, f):
+    f[:, 0] = torch.clamp(1 - (x[:, 0] ** 2 + x[:, 1] ** 2), min=0) * 2
+    return f[:, 0]
+
+
+dim = 2
+bounds = [(-1, 1), (-1, 1)]
+n_eval = 6400000
+batch_size = 10000
+n_therm = 20
+
+vegas_map = Vegas(dim, device=device, ninc=10)
+
+# Monte Carlo and MCMC for Unit Circle
+mc_integrator = MonteCarlo(
+    f=unit_circle_integrand, bounds=bounds, batch_size=batch_size
+)
+mcmc_integrator = MarkovChainMonteCarlo(
+    f=unit_circle_integrand, bounds=bounds, batch_size=batch_size, nburnin=n_therm
+)
+
+print("Unit Circle Integration Results:")
+print("Plain MC:", mc_integrator(n_eval))
+print("MCMC:", mcmc_integrator(n_eval, mix_rate=0.5))
+
+# Train VEGAS map for Unit Circle
+vegas_map.adaptive_training(batch_size, unit_circle_integrand, alpha=0.5)
+vegas_integrator = MonteCarlo(
+    bounds, f=unit_circle_integrand, maps=vegas_map, batch_size=batch_size
+)
+vegasmcmc_integrator = MarkovChainMonteCarlo(
+    bounds,
+    f=unit_circle_integrand,
+    maps=vegas_map,
+    batch_size=batch_size,
+    nburnin=n_therm,
+)
+
+print("VEGAS:", vegas_integrator(n_eval))
+print("VEGAS-MCMC:", vegasmcmc_integrator(n_eval, mix_rate=0.5))
+
+# Monte Carlo and MCMC for Half-Sphere
+mc_integrator.f = half_sphere_integrand
+mcmc_integrator.f = half_sphere_integrand
+
+print("\nHalf-Sphere Integration Results:")
+print("Plain MC:", mc_integrator(n_eval))
+print("MCMC:", mcmc_integrator(n_eval, mix_rate=0.5))
+
+vegas_map.make_uniform()
+vegas_map.adaptive_training(batch_size, half_sphere_integrand, epoch=10, alpha=0.5)
+vegas_integrator.f = half_sphere_integrand
+vegasmcmc_integrator.f = half_sphere_integrand
+
+print("VEGAS:", vegas_integrator(n_eval))
+print("VEGAS-MCMC:", vegasmcmc_integrator(n_eval, mix_rate=0.5))

examples/example_2.py

@@ -0,0 +1,54 @@
+# Example 2: Sharp Peak Integrand in Higher Dimensions
+
+import torch
+from MCintegration import MonteCarlo, MarkovChainMonteCarlo, Vegas, set_seed, get_device
+
+set_seed(42)
+device = get_device()
+
+
+def sharp_integrands(x, f):
+    f[:, 0] = torch.sum((x - 0.5) ** 2, dim=-1)
+    f[:, 0] *= -200
+    f[:, 0].exp_()
+    f[:, 1] = f[:, 0] * x[:, 0]
+    f[:, 2] = f[:, 0] * x[:, 0] ** 2
+    return f.mean(dim=-1)
+
+
+dim = 4
+bounds = [(0, 1)] * dim
+n_eval = 500000
+batch_size = 10000
+n_therm = 20
+
+vegas_map = Vegas(dim, device=device, ninc=1000)
+
+# Plain MC and MCMC
+mc_integrator = MonteCarlo(
+    f=sharp_integrands, f_dim=3, bounds=bounds, batch_size=batch_size
+)
+mcmc_integrator = MarkovChainMonteCarlo(
+    f=sharp_integrands, f_dim=3, bounds=bounds, batch_size=batch_size, nburnin=n_therm
+)
+
+print("Sharp Peak Integration Results:")
+print("Plain MC:", mc_integrator(n_eval))
+print("MCMC:", mcmc_integrator(n_eval, mix_rate=0.5))
+
+# Train VEGAS map
+vegas_map.adaptive_training(batch_size, sharp_integrands, f_dim=3, epoch=10, alpha=2.0)
+vegas_integrator = MonteCarlo(
+    bounds, f=sharp_integrands, f_dim=3, maps=vegas_map, batch_size=batch_size
+)
+vegasmcmc_integrator = MarkovChainMonteCarlo(
+    bounds,
+    f=sharp_integrands,
+    f_dim=3,
+    maps=vegas_map,
+    batch_size=batch_size,
+    nburnin=n_therm,
+)
+
+print("VEGAS:", vegas_integrator(n_eval))
+print("VEGAS-MCMC:", vegasmcmc_integrator(n_eval, mix_rate=0.5))

examples/example_3.py

@@ -0,0 +1,54 @@
+# Example 3: Integration of log(x)/sqrt(x) using VEGAS
+
+import torch
+from MCintegration import MonteCarlo, MarkovChainMonteCarlo
+from MCintegration import Vegas, set_seed, get_device
+
+set_seed(42)
+device = get_device()
+
+
+def func(x, f):
+    f[:, 0] = torch.log(x[:, 0]) / torch.sqrt(x[:, 0])
+    return f[:, 0]
+
+
+dim = 1
+bounds = [[0, 1]] * dim
+n_eval = 500000
+batch_size = 10000
+alpha = 2.0
+ninc = 1000
+n_therm = 20
+
+vegas_map = Vegas(dim, device=device, ninc=ninc)
+
+# Train VEGAS map
+print("Training VEGAS map for log(x)/sqrt(x)... \n")
+vegas_map.adaptive_training(batch_size, func, epoch=10, alpha=alpha)
+
+print("Integration Results for log(x)/sqrt(x):")
+
+
+# Plain MC Integration
+mc_integrator = MonteCarlo(bounds, func, batch_size=batch_size)
+print("Plain MC Integral Result:", mc_integrator(n_eval))
+
+# MCMC Integration
+mcmc_integrator = MarkovChainMonteCarlo(
+    bounds, func, batch_size=batch_size, nburnin=n_therm
+)
+print("MCMC Integral Result:", mcmc_integrator(n_eval, mix_rate=0.5))
+
+# Perform VEGAS integration
+vegas_integrator = MonteCarlo(bounds, func, maps=vegas_map, batch_size=batch_size)
+res = vegas_integrator(n_eval)
+
+print("VEGAS Integral Result:", res)
+
+# VEGAS-MCMC Integration
+vegasmcmc_integrator = MarkovChainMonteCarlo(
+    bounds, func, maps=vegas_map, batch_size=batch_size, nburnin=n_therm
+)
+res_vegasmcmc = vegasmcmc_integrator(n_eval, mix_rate=0.5)
+print("VEGAS-MCMC Integral Result:", res_vegasmcmc)