diff --git a/examples/example_1.py b/examples/example_1.py new file mode 100644 index 0000000..ae0f7f1 --- /dev/null +++ b/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)) diff --git a/examples/example_2.py b/examples/example_2.py new file mode 100644 index 0000000..3af5719 --- /dev/null +++ b/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)) diff --git a/examples/example_3.py b/examples/example_3.py new file mode 100644 index 0000000..b2ac147 --- /dev/null +++ b/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) diff --git a/examples/vegas_test.py b/examples/vegas_test.py deleted file mode 100644 index 57fbe07..0000000 --- a/examples/vegas_test.py +++ /dev/null @@ -1,131 +0,0 @@ -# Integration tests for VEGAS + MonteCarlo/MarkovChainMonteCarlo integral methods. -import torch -from MCintegration import MonteCarlo, MarkovChainMonteCarlo -from MCintegration import Vegas, set_seed, get_device - -set_seed(42) -device = get_device() -# device = torch.device("mps") -dtype = torch.float32 - - -def sharp_peak(x, f): - f[:, 0] = torch.sum((x - 0.5) ** 2, dim=-1) - f[:, 0] *= -200 - f[:, 0].exp_() - return f[:, 0] - - -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) - - -def func(x, f): - f[:, 0] = torch.log(x[:, 0]) / torch.sqrt(x[:, 0]) - return f[:, 0] - - -alpha = 2.0 -ninc = 1000 -# n_eval = 1000000 -n_eval = 500000 -batch_size = 10000 -n_therm = 10 - -print("\nCalculate the integral log(x)/x^0.5 in the bounds [0, 1]") -dim = 1 -bounds = [[0, 1]] * dim -print("train VEGAS map") -vegas_map = Vegas(dim, device=device, ninc=ninc, dtype=dtype) -vegas_map.adaptive_training(100000, func, epoch=10, alpha=alpha) - -vegas_integrator = MonteCarlo( - bounds, - func, - maps=vegas_map, - batch_size=batch_size, -) -res = vegas_integrator(n_eval) -print("VEGAS Integral results: ", res) - -vegasmcmc_integrator = MarkovChainMonteCarlo( - bounds, - func, - maps=vegas_map, - batch_size=batch_size, - nburnin=n_therm, -) -res = vegasmcmc_integrator(n_eval, mix_rate=0.5) -print("VEGAS-MarkovChainMonteCarlo Integral results: ", res) -print(type(res)) -print(res.sum_neval) -print(res.itn_results) -print(res.nitn) - -# Start Monte Carlo integration, including plain-MC, MarkovChainMonteCarlo, vegas, and vegas-MarkovChainMonteCarlo -print("\nCalculate the integral [h(X), x1 * h(X), x1^2 * h(X)] in the bounds [0, 1]^4") -print("h(X) = exp(-200 * (x1^2 + x2^2 + x3^2 + x4^2))") - -dim = 4 -bounds = [(0, 1)] * dim -vegas_map = Vegas(dim, device=device, ninc=ninc, dtype=dtype) -print("train VEGAS map for h(X)...") -vegas_map.adaptive_training(20000, sharp_peak, epoch=10, alpha=alpha) -# print(vegas_map.extract_grid()) - -print("VEGAS Integral results:") -vegas_integrator = MonteCarlo( - bounds, - sharp_integrands, - f_dim=3, - maps=vegas_map, - batch_size=batch_size, -) -res = vegas_integrator(neval=500000) -print( - " I[0] =", - res[0], - " I[1] =", - res[1], - " I[2] =", - res[2], - " I[1]/I[0] =", - res[1] / res[0], -) -print(type(res)) -print(type(res[0])) -print(res[0].sum_neval) -print(res[0].itn_results) -print(res[0].nitn) - - -print("VEGAS-MarkovChainMonteCarlo Integral results:") -vegasmcmc_integrator = MarkovChainMonteCarlo( - bounds, - sharp_integrands, - f_dim=3, - maps=vegas_map, - batch_size=batch_size, - nburnin=n_therm, -) -res = vegasmcmc_integrator(neval=500000, mix_rate=0.5) -print( - " I[0] =", - res[0], - " I[1] =", - res[1], - " I[2] =", - res[2], - " I[1]/I[0] =", - res[1] / res[0], -) -print(type(res)) -print(type(res[0])) -print(res[0].sum_neval) -print(res[0].itn_results) -print(res[0].nitn)