diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 3827e33..2c1c3eb 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -39,6 +39,8 @@ jobs: python -m pip install --upgrade pip pip install coverage pip install pytest pytest-cov + pip install -e .[test] + pip install pytest-xdist pip install -r requirements.txt - name: Run tests @@ -49,7 +51,10 @@ jobs: export MASTER_ADDR=localhost export MASTER_PORT=12345 export PYTHONPATH=MCintegration + export COVERAGE_PROCESS_START=.coveragerc pytest --cov --cov-report=xml --ignore=examples + coverage combine + coverage xml - name: Upload coverage to Codecov uses: codecov/codecov-action@v4 diff --git a/MCintegration/maps_test.py b/MCintegration/maps_test.py index 7863828..7024cfe 100644 --- a/MCintegration/maps_test.py +++ b/MCintegration/maps_test.py @@ -1,11 +1,45 @@ import unittest import torch - -# import numpy as np +import numpy as np from maps import Map, CompositeMap, Vegas, Configuration from base import LinearMap +class TestConfiguration(unittest.TestCase): + def setUp(self): + self.batch_size = 5 + self.dim = 3 + self.f_dim = 2 + self.device = "cpu" + self.dtype = torch.float64 + + def test_configuration_initialization(self): + config = Configuration( + batch_size=self.batch_size, + dim=self.dim, + f_dim=self.f_dim, + device=self.device, + dtype=self.dtype, + ) + + self.assertEqual(config.batch_size, self.batch_size) + self.assertEqual(config.dim, self.dim) + self.assertEqual(config.f_dim, self.f_dim) + self.assertEqual(config.device, self.device) + + self.assertEqual(config.u.shape, (self.batch_size, self.dim)) + self.assertEqual(config.x.shape, (self.batch_size, self.dim)) + self.assertEqual(config.fx.shape, (self.batch_size, self.f_dim)) + self.assertEqual(config.weight.shape, (self.batch_size,)) + self.assertEqual(config.detJ.shape, (self.batch_size,)) + + self.assertEqual(config.u.dtype, self.dtype) + self.assertEqual(config.x.dtype, self.dtype) + self.assertEqual(config.fx.dtype, self.dtype) + self.assertEqual(config.weight.dtype, self.dtype) + self.assertEqual(config.detJ.dtype, self.dtype) + + class TestMap(unittest.TestCase): def setUp(self): self.device = "cpu" @@ -24,6 +58,35 @@ def test_inverse_not_implemented(self): with self.assertRaises(NotImplementedError): self.map.inverse(torch.tensor([0.5, 0.5], dtype=self.dtype)) + def test_forward_with_detJ(self): + # Create a simple linear map for testing: x = u * A + b + # With A=[1, 1] and b=[0, 0], we have x = u + linear_map = LinearMap([1, 1], [0, 0], device=self.device) + + # Test forward_with_detJ method + u = torch.tensor([[0.5, 0.5]], dtype=torch.float64, device=self.device) + x, detJ = linear_map.forward_with_detJ(u) + + # Since it's a linear map from [0,0] to [1,1], x should equal u + self.assertTrue(torch.allclose(x, u)) + + # Determinant of Jacobian should be 1 for linear map with slope 1 + # forward_with_detJ returns actual determinant, not log + self.assertAlmostEqual(detJ.item(), 1.0) + + # Test with a different linear map: x = u * [2, 3] + [1, 1] + # So u = [0.5, 0.5] should give x = [0.5*2+1, 0.5*3+1] = [2, 2.5] + linear_map2 = LinearMap([2, 3], [1, 1], device=self.device) + u2 = torch.tensor([[0.5, 0.5]], dtype=torch.float64, device=self.device) + x2, detJ2 = linear_map2.forward_with_detJ(u2) + expected_x2 = torch.tensor( + [[2.0, 2.5]], dtype=torch.float64, device=self.device + ) + self.assertTrue(torch.allclose(x2, expected_x2)) + + # Determinant should be 2 * 3 = 6 + self.assertAlmostEqual(detJ2.item(), 6.0) + class TestCompositeMap(unittest.TestCase): def setUp(self): @@ -99,6 +162,32 @@ def test_initialization(self): self.assertTrue(torch.equal(self.vegas.grid, self.init_grid)) self.assertEqual(self.vegas.inc.shape, (2, self.ninc)) + def test_ninc_initialization_types(self): + # Test ninc initialization with int + vegas_int = Vegas(self.dim, ninc=5) + self.assertEqual(vegas_int.ninc.tolist(), [5, 5]) + + # Test ninc initialization with list + vegas_list = Vegas(self.dim, ninc=[5, 10]) + self.assertEqual(vegas_list.ninc.tolist(), [5, 10]) + + # Test ninc initialization with numpy array + vegas_np = Vegas(self.dim, ninc=np.array([3, 7])) + self.assertEqual(vegas_np.ninc.tolist(), [3, 7]) + + # Test ninc initialization with torch tensor + vegas_tensor = Vegas(self.dim, ninc=torch.tensor([4, 6])) + self.assertEqual(vegas_tensor.ninc.tolist(), [4, 6]) + + # Test ninc initialization with invalid type + with self.assertRaises(ValueError): + Vegas(self.dim, ninc="invalid") + + def test_ninc_shape_validation(self): + # Test ninc shape validation + with self.assertRaises(ValueError): + Vegas(self.dim, ninc=[1, 2, 3]) # Wrong length + def test_add_training_data(self): # Test adding training data self.vegas.add_training_data(self.sample) @@ -137,6 +226,16 @@ def test_forward(self): self.assertEqual(x.shape, u.shape) self.assertEqual(log_jac.shape, (u.shape[0],)) + def test_forward_with_detJ(self): + # Test forward_with_detJ transformation + u = torch.tensor([[0.1, 0.2], [0.3, 0.4]], dtype=torch.float64) + x, det_jac = self.vegas.forward_with_detJ(u) + self.assertEqual(x.shape, u.shape) + self.assertEqual(det_jac.shape, (u.shape[0],)) + + # Determinant should be positive + self.assertTrue(torch.all(det_jac > 0)) + def test_forward_out_of_bounds(self): # Test forward transformation with out-of-bounds u values u = torch.tensor( diff --git a/MCintegration/mc_multicpu_test.py b/MCintegration/mc_multicpu_test.py index 1be7b9f..242c1a2 100644 --- a/MCintegration/mc_multicpu_test.py +++ b/MCintegration/mc_multicpu_test.py @@ -79,7 +79,10 @@ def two_integrands(x, f): if dist.is_initialized(): dist.destroy_process_group() - +def test_mcmc_singlethread(): + world_size = 1 + init_process(rank=0, world_size=world_size, fn=run_mcmc, backend=backend) + def test_mcmc(world_size=2): # Use fewer processes than CPU cores to avoid resource contention world_size = min(world_size, mp.cpu_count()) diff --git a/MCintegration/utils_test.py b/MCintegration/utils_test.py index d186774..636d8a3 100644 --- a/MCintegration/utils_test.py +++ b/MCintegration/utils_test.py @@ -164,17 +164,6 @@ def test_converged_criteria(self): self.assertTrue(ravg.converged(0.1, 0.1)) self.assertFalse(ravg.converged(0.001, 0.001)) - def test_multiplication_with_another_ravg(self): - ravg1 = RAvg(weighted=True) - ravg1.update(2.0, 0.1) - ravg2 = RAvg(weighted=True) - ravg2.update(3.0, 0.1) - - result = ravg1 * ravg2 - self.assertAlmostEqual(result.mean, 6.0) - sdev = (0.1 / 2**2 + 0.1 / 3**2) ** 0.5 * 6.0 - self.assertAlmostEqual(result.sdev, sdev) - def test_multiplication(self): ravg1 = RAvg(weighted=True) # Test multiplication by another RAvg object @@ -216,16 +205,19 @@ def test_multiplication(self): np.allclose([r.sdev for r in result], [2.0 * ravg1.sdev, 3.0 * ravg1.sdev]) ) - def test_division_with_another_ravg(self): - ravg1 = RAvg(weighted=True) - ravg1.update(6.0, 0.1) - ravg2 = RAvg(weighted=True) - ravg2.update(3.0, 0.1) + # Test multiplication with unweighted RAvg + ravg_unweighted = RAvg(weighted=False) + ravg_unweighted.update(2.0, 0.1) + result = ravg_unweighted * 3.0 + self.assertAlmostEqual(result.mean, 6.0) + self.assertAlmostEqual(result.sdev, ravg_unweighted.sdev * 3) - result = ravg1 / ravg2 - self.assertAlmostEqual(result.mean, 2.0) - sdev = (0.1 / 6.0**2 + 0.1 / 3.0**2) ** 0.5 * 2.0 - self.assertAlmostEqual(result.sdev, sdev) + # Test multiplication with zero variance + ravg_zero_var = RAvg(weighted=True) + ravg_zero_var.update(2.0, 0.0) + result = ravg_zero_var * 3.0 + self.assertAlmostEqual(result.mean, 6.0) + self.assertAlmostEqual(result.sdev, 0.0) def test_division(self): ravg1 = RAvg(weighted=True) @@ -271,6 +263,53 @@ def test_division(self): np.allclose([r.sdev for r in result], [ravg1.sdev / 2.0, ravg1.sdev / 3.0]) ) + # Test division with unweighted RAvg + ravg_unweighted = RAvg(weighted=False) + ravg_unweighted.update(6.0, 0.1) + result = ravg_unweighted / 3.0 + self.assertAlmostEqual(result.mean, 2.0) + self.assertAlmostEqual(result.sdev, ravg_unweighted.sdev / 3) + + # Test division with zero variance + ravg_zero_var = RAvg(weighted=True) + ravg_zero_var.update(6.0, 0.0) + result = ravg_zero_var / 3.0 + self.assertAlmostEqual(result.mean, 2.0) + self.assertAlmostEqual(result.sdev, 0.0) + + # Test division of zero by RAvg + zero_ravg = RAvg(weighted=True) + zero_ravg.update(0.0, 0.1) + divisor_ravg = RAvg(weighted=True) + divisor_ravg.update(3.0, 0.1) + result = zero_ravg / divisor_ravg + self.assertAlmostEqual(result.mean, 0.0) + # sdev = (0.1 / 0.0**2 + 0.1 / 3.0**2) ** 0.5 * 0.0 # This would be NaN + # For 0/x, the error propagation gives 0 * sqrt((0.1/0.0^2) + (0.1/3.0^2)) + # But since we're dividing by zero, we need to be careful + # In practice, gvar handles this appropriately + + def test_vector_operations_not_implemented(self): + # Test that NotImplemented is returned for vector operations + ravg = RAvg(weighted=True) + ravg.update(2.0, 0.1) + + # Test multiplication with list (should return NotImplemented) + result = ravg.__mul__([1, 2, 3]) + self.assertEqual(result, NotImplemented) + + # Test division with list (should return NotImplemented) + result = ravg.__truediv__([1, 2, 3]) + self.assertEqual(result, NotImplemented) + + # Test multiplication with numpy array (should return NotImplemented) + result = ravg.__mul__(np.array([1, 2, 3])) + self.assertEqual(result, NotImplemented) + + # Test division with numpy array (should return NotImplemented) + result = ravg.__truediv__(np.array([1, 2, 3])) + self.assertEqual(result, NotImplemented) + class TestUtils(unittest.TestCase): def setUp(self): diff --git a/README.md b/README.md index bccdd05..2e7618c 100644 --- a/README.md +++ b/README.md @@ -2,3 +2,144 @@ [![alpha](https://img.shields.io/badge/docs-alpha-blue.svg)](https://numericaleft.github.io/MCintegration.py/) [![Build Status](https://github.com/numericalEFT/MCIntegration.py/workflows/CI/badge.svg)](https://github.com/numericalEFT/MCIntegration.py/actions) [![codecov](https://codecov.io/gh/numericalEFT/MCintegration.py/graph/badge.svg?token=851N2CNOTN)](https://codecov.io/gh/numericalEFT/MCintegration.py) +A Python library for Monte Carlo integration with support for multi-CPU and GPU computations. + +## Overview + +MCintegration is a specialized library designed for numerical integration using Monte Carlo methods. It provides efficient implementations of various integration algorithms with focus on applications in computational physics and effective field theories (EFT). + +The library offers: +- Multiple Monte Carlo integration algorithms +- Support for multi-CPU parallelization +- GPU acceleration capabilities +- Integration with PyTorch for tensor-based computations + +## Installation + +```bash +pip install mcintegration +``` + +Or install from source: + +```bash +python setup.py install +``` + +## Usage + +### Example 1: Unit Circle Integration + +This example demonstrates different Monte Carlo methods for integrating functions over [-1,1]×[-1,1]: + +```python +from MCintegration import MonteCarlo, MarkovChainMonteCarlo, Vegas +import torch + +# Define integrand function +def unit_circle(x, f): + r2 = x[:, 0]**2 + x[:, 1]**2 + f[:, 0] = (r2 <= 1).float() + return f.mean(dim=-1) + +# Set up integration parameters +dim = 2 +bounds = [(-1, 1)] * dim +n_eval = 6400000 +batch_size = 10000 +n_therm = 100 + +# Create integrator instances +mc = MonteCarlo(f=unit_circle, bounds=bounds, batch_size=batch_size) +mcmc = MarkovChainMonteCarlo(f=unit_circle, bounds=bounds, batch_size=batch_size, nburnin=n_therm) + +# Perform integration +result_mc = mc(n_eval) +result_mcmc = mcmc(n_eval) +``` + +### Example 2: Singular Function Integration + +This example shows integration of a function with a singularity at x=0: + +```python +# Integrate log(x)/sqrt(x) which has a singularity at x=0 +def singular_func(x, f): + f[:, 0] = torch.log(x[:, 0]) / torch.sqrt(x[:, 0]) + return f[:, 0] + +# Set up integration parameters +dim = 1 +bounds = [(0, 1)] +n_eval = 6400000 +batch_size = 10000 +n_therm = 100 + +# Use VEGAS algorithm which adapts to the singular structure +vegas_map = Vegas(dim, ninc=1000) +vegas_map.adaptive_training(batch_size, singular_func) + +# Create integrator instances using the adapted vegas map +vegas_mc = MonteCarlo(f=singular_func, bounds=bounds, batch_size=batch_size, maps=vegas_map) +vegas_mcmc = MarkovChainMonteCarlo(f=singular_func, bounds=bounds, batch_size=batch_size, nburnin=n_therm, maps=vegas_map) + +# Perform integration +result_vegas = vegas_mc(n_eval) +result_vegas_mcmc = vegas_mcmc(n_eval) +``` + +### Example 3: Multiple Sharp Peak Integrands in Higher Dimensions + +This example demonstrates integration of a sharp Gaussian peak and its moments in 4D space: + +```python +# Define a sharp peak and its moments integrands +# This represents a Gaussian peak centered at (0.5, 0.5, 0.5, 0.5) +def sharp_integrands(x, f): + f[:, 0] = torch.sum((x - 0.5) ** 2, dim=-1) # Distance from center + f[:, 0] *= -200 # Scale by width parameter + f[:, 0].exp_() # Exponentiate to create Gaussian + f[:, 1] = f[:, 0] * x[:, 0] # First moment + f[:, 2] = f[:, 0] * x[:, 0] ** 2 # Second moment + return f.mean(dim=-1) + +# Set up 4D integration with sharp peak +dim = 4 +bounds = [(0, 1)] * dim +n_eval = 6400000 +batch_size = 10000 +n_therm = 100 + +# Use VEGAS algorithm which adapts to the peak structure +vegas_map = Vegas(dim, ninc=1000) +vegas_map.adaptive_training(batch_size, sharp_integrands, f_dim=3) + +# Create integrator instances using the adapted vegas map +vegas_mc = MonteCarlo(f=sharp_integrands, f_dim=3, bounds=bounds, batch_size=batch_size, maps=vegas_map) +vegas_mcmc = MarkovChainMonteCarlo(f=sharp_integrands, f_dim=3, bounds=bounds, batch_size=batch_size, nburnin=n_therm, maps=vegas_map) + +# Perform integration +result_vegas = vegas_mc(n_eval) +result_vegas_mcmc = vegas_mcmc(n_eval) +``` + +## Features + +- **Base integration methods**: Core Monte Carlo algorithms in `MCintegration/base.py` +- **Integrator implementations**: Various MC integration strategies in `MCintegration/integrators.py` +- **Variable transformations**: Coordinate mapping utilities in `MCintegration/maps.py` +- **Utility functions**: Helper functions for numerical computations in `MCintegration/utils.py` +- **Multi-CPU support**: Parallel processing capabilities demonstrated in `MCintegration/mc_multicpu_test.py` +- **GPU acceleration**: CUDA-enabled functions through PyTorch in the examples directory + + +## Requirements + +- Python 3.7+ +- NumPy +- PyTorch +- gvar + +## Acknowledgements and Related Packages +The development of `MCIntegration.py` has been greatly inspired and influenced by `vegas` package. We would like to express our appreciation to the following: +- [vegas](https://github.com/gplepage/vegas) A Python package offering Monte Carlo estimations of multidimensional integrals, with notable improvements on the original Vegas algorithm. It's been a valuable reference for us. Learn more from the vegas [documentation](https://vegas.readthedocs.io/). **Reference: G. P. Lepage, J. Comput. Phys. 27, 192 (1978) and G. P. Lepage, J. Comput. Phys. 439, 110386 (2021) [arXiv:2009.05112](https://arxiv.org/abs/2009.05112)**. \ No newline at end of file diff --git a/license.md b/license.md new file mode 100644 index 0000000..b964bdf --- /dev/null +++ b/license.md @@ -0,0 +1,7 @@ +Copyright (c) 2025: Pengcheng Hou, Tao Wang, Caiyu Fan, and Kun Chen. + +Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. \ No newline at end of file