diff --git a/.github/workflows/check-msrv.yml b/.github/workflows/check-msrv.yml new file mode 100644 index 0000000..5456f5e --- /dev/null +++ b/.github/workflows/check-msrv.yml @@ -0,0 +1,20 @@ +name: Check MSRV + +on: [push] + +jobs: + msrv: + runs-on: ubuntu-latest + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Install Rust toolchain + uses: dtolnay/rust-toolchain@stable + + - name: Run check + run: | + # extract the MSRV from Root `Cargo.toml` + MSRV=$(grep '^rust-version ' Cargo.toml | cut -d= -f2- | tr -d ' "') + rustup default "${MSRV}" + cargo check --all-targets diff --git a/.github/workflows/run-rust.yml b/.github/workflows/run-rust.yml new file mode 100644 index 0000000..c40a15b --- /dev/null +++ b/.github/workflows/run-rust.yml @@ -0,0 +1,59 @@ +name: Run the Rust tests 🦀 + +on: + push: + +env: + CARGO_TERM_COLOR: always + +jobs: + build: + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + runs-on: ${{ matrix.os }} + + steps: + - uses: actions/checkout@v4 + + - name: Switch to nightly Rust + run: | + rustup default nightly + rustup component add llvm-tools-preview + + - name: Set RUSTDOCFLAGS + run: | + echo "RUSTDOCFLAGS=-Cinstrument-coverage -Z unstable-options --persist-doctests $(pwd)/target/debug/doctestbins" >> "$GITHUB_ENV" + + - name: Run tests 🦀 + env: + RUSTFLAGS: '-Cinstrument-coverage -Clink-dead-code' + run: | + cargo test --workspace --exclude mchep_pyapi --no-fail-fast 2> >(tee stderr 1>&2) + if [ "${{ matrix.os }}" != "macos-latest" ]; then + sed -i 's/\x1B\[[0-9;]\{1,\}[A-Za-z]//g' stderr + fi + + - name: Generate code coverage + if: matrix.os == 'ubuntu-latest' + run: | + find . -name '*.profraw' -exec $(rustc --print target-libdir)/../bin/llvm-profdata merge -sparse -o mchep.profdata {} + + sed -nE 's/[[:space:]]+Running( unittests|) [^[:space:]]+ \(([^)]+)\)/\2/p' stderr | \ + xargs printf ' --object %s' | \ + xargs $(rustc --print target-libdir)/../bin/llvm-cov export \ + --ignore-filename-regex=index.crates.io \ + --ignore-filename-regex=rustc \ + --ignore-filename-regex=mchep/tests \ + --ignore-filename-regex=mchep_capi \ + --instr-profile=mchep.profdata \ + --skip-functions \ + --format lcov > lcov.info + grep SF lcov.info | sort -u | sed 's/SF://' + + - name: Upload to codecov.io + if: matrix.os == 'ubuntu-latest' + uses: codecov/codecov-action@v4 + with: + token: ${{ secrets.CODECOV_TOKEN }} + flags: rust diff --git a/Cargo.lock b/Cargo.lock index 62c81e3..379f25c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -11,6 +11,18 @@ dependencies = [ "memchr", ] +[[package]] +name = "anes" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b46cbb362ab8752921c97e041f5e366ee6297bd428a31275b9fcf1e380f7299" + +[[package]] +name = "anstyle" +version = "1.0.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5192cca8006f1fd4f7237516f40fa183bb07f8fbdfedaa0036de5ea9b0b45e78" + [[package]] name = "approx" version = "0.5.1" @@ -20,13 +32,19 @@ dependencies = [ "num-traits", ] +[[package]] +name = "assert_approx_eq" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c07dab4369547dbe5114677b33fbbf724971019f3818172d59a97a61c774ffd" + [[package]] name = "atty" version = "0.2.14" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" dependencies = [ - "hermit-abi", + "hermit-abi 0.1.19", "libc", "winapi", ] @@ -46,7 +64,7 @@ dependencies = [ "bitflags 2.10.0", "cexpr", "clang-sys", - "itertools", + "itertools 0.12.1", "lazy_static", "lazycell", "log", @@ -82,19 +100,31 @@ dependencies = [ "shell-words", ] +[[package]] +name = "bumpalo" +version = "3.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "46c5e41b57b8bba42a04676d81cb89e9ee8e859a1a66f80a5a72e1cb76b34d43" + [[package]] name = "bytemuck" version = "1.24.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "1fbdf580320f38b612e485521afda1ee26d10cc9884efaaa750d383e13e3c5f4" +[[package]] +name = "cast" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" + [[package]] name = "cbindgen" version = "0.26.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "da6bc11b07529f16944307272d5bd9b22530bc7d05751717c9d416586cedab49" dependencies = [ - "clap", + "clap 3.2.25", "heck", "indexmap", "log", @@ -109,9 +139,9 @@ dependencies = [ [[package]] name = "cc" -version = "1.2.45" +version = "1.2.46" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "35900b6c8d709fb1d854671ae27aeaa9eec2f8b01b364e1619a40da3e6fe2afe" +checksum = "b97463e1064cb1b1c1384ad0a0b9c8abd0988e2a91f52606c80ef14aadb63e36" dependencies = [ "find-msvc-tools", "shlex", @@ -132,6 +162,33 @@ version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" +[[package]] +name = "ciborium" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42e69ffd6f0917f5c029256a24d0161db17cea3997d185db0d35926308770f0e" +dependencies = [ + "ciborium-io", + "ciborium-ll", + "serde", +] + +[[package]] +name = "ciborium-io" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05afea1e0a06c9be33d539b876f1ce3692f4afea2cb41f740e7743225ed1c757" + +[[package]] +name = "ciborium-ll" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57663b653d948a338bfb3eeba9bb2fd5fcfaecb9e199e87e1eda4d9e8b240fd9" +dependencies = [ + "ciborium-io", + "half", +] + [[package]] name = "clang-sys" version = "1.8.1" @@ -151,13 +208,32 @@ checksum = "4ea181bf566f71cb9a5d17a59e1871af638180a18fb0035c92ae62b705207123" dependencies = [ "atty", "bitflags 1.3.2", - "clap_lex", + "clap_lex 0.2.4", "indexmap", "strsim", "termcolor", "textwrap", ] +[[package]] +name = "clap" +version = "4.5.52" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa8120877db0e5c011242f96806ce3c94e0737ab8108532a76a3300a01db2ab8" +dependencies = [ + "clap_builder", +] + +[[package]] +name = "clap_builder" +version = "4.5.52" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02576b399397b659c26064fbc92a75fede9d18ffd5f80ca1cd74ddab167016e1" +dependencies = [ + "anstyle", + "clap_lex 0.7.6", +] + [[package]] name = "clap_lex" version = "0.2.4" @@ -167,6 +243,12 @@ dependencies = [ "os_str_bytes", ] +[[package]] +name = "clap_lex" +version = "0.7.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a1d728cc89cf3aee9ff92b05e62b19ee65a02b5702cff7d5a377e32c6ae29d8d" + [[package]] name = "conv" version = "0.3.3" @@ -176,6 +258,42 @@ dependencies = [ "custom_derive", ] +[[package]] +name = "criterion" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2b12d017a929603d80db1831cd3a24082f8137ce19c69e6447f54f5fc8d692f" +dependencies = [ + "anes", + "cast", + "ciborium", + "clap 4.5.52", + "criterion-plot", + "is-terminal", + "itertools 0.10.5", + "num-traits", + "once_cell", + "oorandom", + "plotters", + "rayon", + "regex", + "serde", + "serde_derive", + "serde_json", + "tinytemplate", + "walkdir", +] + +[[package]] +name = "criterion-plot" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1" +dependencies = [ + "cast", + "itertools 0.10.5", +] + [[package]] name = "crossbeam-deque" version = "0.8.6" @@ -201,6 +319,12 @@ version = "0.8.21" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28" +[[package]] +name = "crunchy" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" + [[package]] name = "cust" version = "0.3.2" @@ -277,9 +401,9 @@ checksum = "37909eebbb50d72f9059c3b6d82c0463f2ff062c9e95845c43a6c9c0355411be" [[package]] name = "find-msvc-tools" -version = "0.1.4" +version = "0.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "52051878f80a721bb68ebfbc930e07b65ba72f2da88968ea5c06fd6ca3d3a127" +checksum = "3a3076410a55c90011c298b04d0cfa770b00fa04e1e3c97d3f6c9de105a03844" [[package]] name = "find_cuda_helper" @@ -328,6 +452,17 @@ version = "0.3.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "0cc23270f6e1808e30a928bdc84dea0b9b4136a8bc82338574f23baf47bbd280" +[[package]] +name = "half" +version = "2.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ea2d84b969582b4b1864a92dc5d27cd2b77b622a8d79306834f1be5ba20d84b" +dependencies = [ + "cfg-if", + "crunchy", + "zerocopy", +] + [[package]] name = "hashbrown" version = "0.12.3" @@ -349,6 +484,12 @@ dependencies = [ "libc", ] +[[package]] +name = "hermit-abi" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0fef456e4baa96da950455cd02c081ca953b141298e41db3fc7e36b1da849c" + [[package]] name = "home" version = "0.5.12" @@ -377,6 +518,26 @@ dependencies = [ "rustversion", ] +[[package]] +name = "is-terminal" +version = "0.4.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3640c1c38b8e4e43584d8df18be5fc6b0aa314ce6ebf51b53313d4306cca8e46" +dependencies = [ + "hermit-abi 0.5.2", + "libc", + "windows-sys 0.61.2", +] + +[[package]] +name = "itertools" +version = "0.10.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" +dependencies = [ + "either", +] + [[package]] name = "itertools" version = "0.12.1" @@ -392,6 +553,16 @@ version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4a5f13b858c8d314ee3e8f639011f7ccefe71f97f96e50151fb991f267928e2c" +[[package]] +name = "js-sys" +version = "0.3.82" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b011eec8cc36da2aab2d5cff675ec18454fad408585853910a202391cf9f8e65" +dependencies = [ + "once_cell", + "wasm-bindgen", +] + [[package]] name = "lazy_static" version = "1.5.0" @@ -476,11 +647,17 @@ checksum = "34080505efa8e45a4b816c349525ebe327ceaa8559756f0356cba97ef3bf7432" name = "mchep" version = "0.1.0" dependencies = [ + "assert_approx_eq", + "criterion", "cust", + "itertools 0.12.1", "mpi", + "num-traits", "rand", "rand_pcg", "rayon", + "thiserror", + "wide", ] [[package]] @@ -489,6 +666,7 @@ version = "0.1.0" dependencies = [ "cbindgen", "mchep", + "wide", ] [[package]] @@ -497,6 +675,7 @@ version = "0.1.0" dependencies = [ "mchep", "pyo3", + "wide", ] [[package]] @@ -587,6 +766,12 @@ version = "1.21.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "42f5e15c9953c5e4ccceeb2e7382a716482c34515315f7b03532b8b4e8393d2d" +[[package]] +name = "oorandom" +version = "11.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" + [[package]] name = "os_str_bytes" version = "6.6.1" @@ -622,6 +807,34 @@ version = "0.3.32" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7edddbd0b52d732b21ad9a5fab5c704c14cd949e5e9a1ec5929a24fded1b904c" +[[package]] +name = "plotters" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5aeb6f403d7a4911efb1e33402027fc44f29b5bf6def3effcc22d7bb75f2b747" +dependencies = [ + "num-traits", + "plotters-backend", + "plotters-svg", + "wasm-bindgen", + "web-sys", +] + +[[package]] +name = "plotters-backend" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "df42e13c12958a16b3f7f4386b9ab1f3e7933914ecea48da7139435263a4172a" + +[[package]] +name = "plotters-svg" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "51bae2ac328883f7acdfea3d66a7c35751187f870bc81f94563733a154d7a670" +dependencies = [ + "plotters-backend", +] + [[package]] name = "portable-atomic" version = "1.11.1" @@ -884,6 +1097,24 @@ version = "1.0.20" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "28d3b2b1366ec20994f1fd18c3c594f05c5dd4bc44d8bb0c1c632c8d6829481f" +[[package]] +name = "safe_arch" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "same-file" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93fc1dc3aaa9bfed95e02e6eadabb4baf7e3078b0bd1b4d7b6b0b68378900502" +dependencies = [ + "winapi-util", +] + [[package]] name = "scopeguard" version = "1.2.0" @@ -1039,6 +1270,16 @@ dependencies = [ "syn 2.0.110", ] +[[package]] +name = "tinytemplate" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "be4d6b5f19ff7664e8c98d03e2139cb510db9b0a60b55f8e8709b689d939b6bc" +dependencies = [ + "serde", + "serde_json", +] + [[package]] name = "toml" version = "0.5.11" @@ -1072,6 +1313,16 @@ dependencies = [ "rustc_version", ] +[[package]] +name = "walkdir" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" +dependencies = [ + "same-file", + "winapi-util", +] + [[package]] name = "wasi" version = "0.11.1+wasi-snapshot-preview1" @@ -1087,6 +1338,61 @@ dependencies = [ "wit-bindgen", ] +[[package]] +name = "wasm-bindgen" +version = "0.2.105" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "da95793dfc411fbbd93f5be7715b0578ec61fe87cb1a42b12eb625caa5c5ea60" +dependencies = [ + "cfg-if", + "once_cell", + "rustversion", + "wasm-bindgen-macro", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-macro" +version = "0.2.105" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "04264334509e04a7bf8690f2384ef5265f05143a4bff3889ab7a3269adab59c2" +dependencies = [ + "quote", + "wasm-bindgen-macro-support", +] + +[[package]] +name = "wasm-bindgen-macro-support" +version = "0.2.105" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "420bc339d9f322e562942d52e115d57e950d12d88983a14c79b86859ee6c7ebc" +dependencies = [ + "bumpalo", + "proc-macro2", + "quote", + "syn 2.0.110", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-shared" +version = "0.2.105" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "76f218a38c84bcb33c25ec7059b07847d465ce0e0a76b995e134a45adcb6af76" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "web-sys" +version = "0.3.82" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3a1f95c0d03a47f4ae1f7a64643a6bb97465d9b740f0fa8f90ea33915c99a9a1" +dependencies = [ + "js-sys", + "wasm-bindgen", +] + [[package]] name = "which" version = "4.4.2" @@ -1099,6 +1405,16 @@ dependencies = [ "rustix 0.38.44", ] +[[package]] +name = "wide" +version = "0.7.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi" version = "0.3.9" diff --git a/Cargo.toml b/Cargo.toml index 8773cd2..e293039 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,11 +11,27 @@ categories = ["science", "algorithms"] keywords = ["monte-carlo", "integration", "physics"] [workspace.dependencies] -rand = "0.8" -rayon = "1.5" -rand_pcg = "0.3" -mpi = { version = "0.7" } +rand = "0.8.5" +rand_pcg = "0.3.1" +rayon = "1.8.0" +itertools = "0.12.0" +thiserror = "1.0.56" +num-traits = "0.2.17" cust = "0.3.2" +mpi = "0.7" +pyo3 = "0.21.2" +rust-version = "1.82.0" +wide = "0.7.5" + +[workspace.lints.clippy] +all = { level = "warn", priority = -1 } +cargo = { level = "warn", priority = -1 } +nursery = { level = "warn", priority = -1 } +pedantic = { level = "warn", priority = -1 } + +[workspace.lints.rust] +missing-docs = "warn" +unsafe-op-in-unsafe-fn = "deny" [profile.release] opt-level = 3 diff --git a/mchep/Cargo.toml b/mchep/Cargo.toml index 464e9c3..1ac83fa 100644 --- a/mchep/Cargo.toml +++ b/mchep/Cargo.toml @@ -8,12 +8,31 @@ categories.workspace = true keywords.workspace = true [dependencies] -rand = { workspace = true, features = ["std", "std_rng"] } +rand = { workspace = true } rand_pcg = { workspace = true } rayon = { workspace = true } -mpi = { workspace = true, optional = true } -cust = { workspace = true, optional = true } +itertools = { workspace = true } +thiserror = { workspace = true } +num-traits = { workspace = true } +wide = { workspace = true } [features] -mpi = ["dep:mpi"] gpu = ["dep:cust"] +mpi = ["dep:mpi"] + +[dependencies.cust] +workspace = true +optional = true + +[dependencies.mpi] +workspace = true +optional = true +features = ["user-operations"] + +[dev-dependencies] +assert_approx_eq = "1.1.0" +criterion = "0.5" + +[[bench]] +name = "vegas_benchmark" +harness = false diff --git a/mchep/benches/vegas_benchmark.rs b/mchep/benches/vegas_benchmark.rs new file mode 100644 index 0000000..8d9a510 --- /dev/null +++ b/mchep/benches/vegas_benchmark.rs @@ -0,0 +1,102 @@ +use criterion::{black_box, criterion_group, criterion_main, Criterion}; +use wide::f64x4; + +use mchep::integrand::{Integrand, SimdIntegrand}; +use mchep::vegas::Vegas; + +// Scalar integrand +struct GaussianIntegrand; + +impl Integrand for GaussianIntegrand { + fn dim(&self) -> usize { + 2 + } + fn eval(&self, x: &[f64]) -> f64 { + (-(x[0].powi(2)) - x[1].powi(2)).exp() + } +} + +// SIMD integrand +struct GaussianSimdIntegrand; +impl SimdIntegrand for GaussianSimdIntegrand { + fn dim(&self) -> usize { + 2 + } + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + let x = points[0]; + let y = points[1]; + (-(x * x) - (y * y)).exp() + } +} + +// Scalar 3D integrand +struct Complex3DIntegrand; +impl Integrand for Complex3DIntegrand { + fn dim(&self) -> usize { + 3 + } + fn eval(&self, x: &[f64]) -> f64 { + x[0].sin() * x[1].cos() * (-(x[2] * x[2])).exp() + } +} + +// SIMD 3D integrand +struct Complex3DSimdIntegrand; +impl SimdIntegrand for Complex3DSimdIntegrand { + fn dim(&self) -> usize { + 3 + } + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + let x = points[0]; + let y = points[1]; + let z = points[2]; + x.sin() * y.cos() * (-(z * z)).exp() + } +} + +fn vegas_benchmark(c: &mut Criterion) { + let mut group = c.benchmark_group("Vegas 2D Gaussian"); + let boundaries2d = &[(-1.0, 1.0), (-1.0, 1.0)]; + let n_eval2d = 100_000; + + group.bench_function("Scalar (Rayon)", |b| { + b.iter(|| { + let mut vegas = Vegas::new(1, n_eval2d, 50, 0.5, boundaries2d); + vegas.set_seed(1234); + vegas.integrate(black_box(&GaussianIntegrand), None); + }) + }); + + group.bench_function("SIMD", |b| { + b.iter(|| { + let mut vegas = Vegas::new(1, n_eval2d, 50, 0.5, boundaries2d); + vegas.set_seed(1234); + vegas.integrate_simd(black_box(&GaussianSimdIntegrand), None); + }) + }); + group.finish(); + + let mut group2 = c.benchmark_group("Vegas 3D Complex"); + let boundaries3d = &[(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]; + let n_eval3d = 100_000; + + group2.bench_function("Scalar (Rayon)", |b| { + b.iter(|| { + let mut vegas = Vegas::new(1, n_eval3d, 50, 0.5, boundaries3d); + vegas.set_seed(1234); + vegas.integrate(black_box(&Complex3DIntegrand), None); + }) + }); + + group2.bench_function("SIMD", |b| { + b.iter(|| { + let mut vegas = Vegas::new(1, n_eval3d, 50, 0.5, boundaries3d); + vegas.set_seed(1234); + vegas.integrate_simd(black_box(&Complex3DSimdIntegrand), None); + }) + }); + group2.finish(); +} + +criterion_group!(benches, vegas_benchmark); +criterion_main!(benches); diff --git a/mchep/src/grid.rs b/mchep/src/grid.rs index 26ebbcf..a5eaf3b 100644 --- a/mchep/src/grid.rs +++ b/mchep/src/grid.rs @@ -1,5 +1,7 @@ //! The adaptive grid used by the VEGAS algorithm. +use wide::f64x4; + /// Represents the adaptive grid for a single dimension. #[derive(Debug, Clone)] pub struct Grid { @@ -54,6 +56,41 @@ impl Grid { (bin_index, x, jacobian) } + /// Maps a packet of `4*y` values to `x` values and jacobians using SIMD. + pub fn map_simd(&self, y_packet: f64x4) -> (f64x4, f64x4, [usize; 4]) { + let n_bins_v = f64x4::splat(self.n_bins as f64); + let y_scaled = y_packet * n_bins_v; + let bin_indices_v = y_scaled.floor(); + let y_frac = y_scaled - bin_indices_v; + + let bin_indices_arr_f64 = bin_indices_v.to_array(); + let final_bin_indices = [ + (bin_indices_arr_f64[0] as usize).min(self.n_bins - 1), + (bin_indices_arr_f64[1] as usize).min(self.n_bins - 1), + (bin_indices_arr_f64[2] as usize).min(self.n_bins - 1), + (bin_indices_arr_f64[3] as usize).min(self.n_bins - 1), + ]; + + let x_low = f64x4::new([ + self.bins[final_bin_indices[0]], + self.bins[final_bin_indices[1]], + self.bins[final_bin_indices[2]], + self.bins[final_bin_indices[3]], + ]); + let x_high = f64x4::new([ + self.bins[final_bin_indices[0] + 1], + self.bins[final_bin_indices[1] + 1], + self.bins[final_bin_indices[2] + 1], + self.bins[final_bin_indices[3] + 1], + ]); + + let width = x_high - x_low; + let x = x_low + y_frac * width; + let jacobian = width * n_bins_v; + + (x, jacobian, final_bin_indices) + } + /// Refines the grid based on the accumulated importance data. pub fn refine(&mut self) { let mut smoothed_d = self.d.clone(); diff --git a/mchep/src/integrand.rs b/mchep/src/integrand.rs index f0f4f2f..0381800 100644 --- a/mchep/src/integrand.rs +++ b/mchep/src/integrand.rs @@ -1,5 +1,7 @@ //! The `Integrand` trait, which defines the function to be integrated. +use wide::f64x4; + /// A trait representing a function to be integrated. /// /// Users of the library must implement this trait for their function. @@ -19,6 +21,15 @@ pub trait Integrand { fn eval(&self, x: &[f64]) -> f64; } +/// A trait representing a function to be integrated using SIMD. +pub trait SimdIntegrand { + /// Returns the number of dimensions of the integration space. + fn dim(&self) -> usize; + + /// Evaluates the function on a packet of 4 points. + fn eval_simd(&self, points: &[f64x4]) -> f64x4; +} + /// A trait representing a function to be integrated on the GPU. /// /// Users of the library must implement this trait for their function. diff --git a/mchep/src/mpi.rs b/mchep/src/mpi.rs index b0d8363..b6f2a72 100644 --- a/mchep/src/mpi.rs +++ b/mchep/src/mpi.rs @@ -21,6 +21,7 @@ impl VegasPlus { &mut self, integrand: &F, world: &SystemCommunicator, + target_accuracy: Option, ) -> VegasResult { assert_eq!(integrand.dim(), self.dim); @@ -48,13 +49,15 @@ impl VegasPlus { self.deserialize_and_update_state(&global_data); + let mut stop_signal = 0i32; + if rank == 0 { let mut global_iter_val = 0.0; - world.reduce_into_root(&iter_val, &mut global_iter_val, SystemOperation::sum()); + world.process_at_rank(0).reduce_into_root(&iter_val, &mut global_iter_val, SystemOperation::sum()); let mut global_iter_err_sq = 0.0; let iter_err_sq = iter_err.powi(2); - world.reduce_into_root( + world.process_at_rank(0).reduce_into_root( &iter_err_sq, &mut global_iter_err_sq, SystemOperation::sum(), @@ -64,10 +67,29 @@ impl VegasPlus { iter_results.push(global_iter_val); iter_errors.push(global_iter_err_sq.sqrt()); } + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + stop_signal = 1; + } + } + } + } } else { - world.reduce_into(&iter_val, SystemOperation::sum(), 0); + world.process_at_rank(0).reduce_into(&iter_val, SystemOperation::sum()); let iter_err_sq = iter_err.powi(2); - world.reduce_into(&iter_err_sq, SystemOperation::sum(), 0); + world.process_at_rank(0).reduce_into(&iter_err_sq, SystemOperation::sum()); + } + + world.process_at_rank(0).broadcast_into(&mut stop_signal); + + if stop_signal == 1 { + break; } self.reallocate_samples(); diff --git a/mchep/src/vegas.rs b/mchep/src/vegas.rs index eeabe5b..a99782f 100644 --- a/mchep/src/vegas.rs +++ b/mchep/src/vegas.rs @@ -1,10 +1,13 @@ //! The main VEGAS integrator. -use crate::grid::Grid; -use crate::integrand::Integrand; use rand::{Rng, SeedableRng}; use rand_pcg::Pcg64; use rayon::prelude::*; +use std::convert::TryInto; +use wide::f64x4; + +use crate::grid::Grid; +use crate::integrand::{Integrand, SimdIntegrand}; /// Stores the result of a VEGAS integration. #[repr(C)] @@ -30,7 +33,7 @@ pub struct Vegas { rng: Pcg64, /// The adaptive grids for each dimension. grids: Vec, - /// The integration boundaries for each dimension, as (min, max) tuples. + /// The integration boundaries for each dimension. boundaries: Vec<(f64, f64)>, } @@ -43,7 +46,7 @@ impl Vegas { /// * `n_eval`: The number of integrand evaluations per iteration. /// * `n_bins`: The number of bins for the adaptive grid in each dimension. /// * `alpha`: The grid damping factor. Should be between 0.0 and 1.0. - /// * `boundaries`: A slice of `(min, max)` tuples defining the integration domain for each dimension. + /// * `boundaries`: The integration domain for each dimension. pub fn new( n_iter: usize, n_eval: usize, @@ -79,7 +82,11 @@ impl Vegas { } /// Integrates the given function using the VEGAS algorithm. - pub fn integrate(&mut self, integrand: &F) -> VegasResult { + pub fn integrate( + &mut self, + integrand: &F, + target_accuracy: Option, + ) -> VegasResult { assert_eq!( integrand.dim(), self.dim, @@ -95,6 +102,19 @@ impl Vegas { if iter > 0 { iter_results.push(iter_val); iter_errors.push(iter_err); + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + return current_result; + } + } + } + } } for grid in &mut self.grids { @@ -106,6 +126,50 @@ impl Vegas { self.combine_results(&iter_results, &iter_errors) } + /// Integrates the given function using the VEGAS algorithm with SIMD. + pub fn integrate_simd( + &mut self, + integrand: &F, + target_accuracy: Option, + ) -> VegasResult { + assert_eq!( + integrand.dim(), + self.dim, + "Integrand dimension does not match integrator dimension." + ); + + let mut iter_results = Vec::new(); + let mut iter_errors = Vec::new(); + + for iter in 0..self.n_iter { + let (iter_val, iter_err) = self.run_iteration_simd(integrand); + + if iter > 0 { + iter_results.push(iter_val); + iter_errors.push(iter_err); + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + return current_result; + } + } + } + } + } + + for grid in &mut self.grids { + grid.refine(); + } + } + + self.combine_results(&iter_results, &iter_errors) + } + /// Runs a single iteration of the VEGAS algorithm in parallel. fn run_iteration(&mut self, integrand: &F) -> (f64, f64) { for grid in &mut self.grids { @@ -113,65 +177,162 @@ impl Vegas { } let n_bins = self.grids[0].n_bins(); - let initial_d_updates: Vec> = (0..self.dim).map(|_| vec![0.0; n_bins]).collect(); + let n_eval = self.n_eval; + let dim = self.dim; + + let mut random_ys: Vec = vec![0.0; n_eval * dim]; + self.rng.fill(&mut random_ys[..]); + + let (sum_f, sum_f2, d_updates, _, _) = random_ys + .par_chunks_exact(dim) + .fold( + || { + ( + 0.0, + 0.0, + vec![vec![0.0; n_bins]; dim], + vec![0.0; dim], + vec![0; dim], + ) + }, + |mut acc, y_vec| { + let (sum_f_thread, sum_f2_thread, d_updates_thread, point, bin_indices) = + &mut acc; + + let mut jacobian = 1.0; + for d in 0..dim { + let y = y_vec[d]; + let (bin_idx, x_unit, jac_vegas) = self.grids[d].map(y); + jacobian *= jac_vegas; + bin_indices[d] = bin_idx; + + let (min, max) = self.boundaries[d]; + point[d] = min + (max - min) * x_unit; + jacobian *= max - min; + } + + let f_val = integrand.eval(point); + let weighted_f = f_val * jacobian; + let f2 = weighted_f * weighted_f; - let random_ys: Vec> = (0..self.n_eval) - .map(|_| (0..self.dim).map(|_| self.rng.gen()).collect()) - .collect(); + *sum_f_thread += weighted_f; + *sum_f2_thread += f2; + + let d_val = f2 / n_eval as f64; + for d in 0..dim { + d_updates_thread[d][bin_indices[d]] += d_val; + } - let (sum_f, sum_f2, d_updates) = random_ys + acc + }, + ) + .reduce( + || { + ( + 0.0, + 0.0, + vec![vec![0.0; n_bins]; dim], + vec![0.0; dim], + vec![0; dim], + ) + }, + |mut a, b| { + a.0 += b.0; + a.1 += b.1; + for d in 0..dim { + for i in 0..n_bins { + a.2[d][i] += b.2[d][i]; + } + } + a + }, + ); + + for d in 0..self.dim { + self.grids[d].d.copy_from_slice(&d_updates[d]); + } + + let avg_f = sum_f / self.n_eval as f64; + let avg_f2 = sum_f2 / self.n_eval as f64; + let variance = (avg_f2 - avg_f * avg_f) / (self.n_eval - 1).max(1) as f64; + let error = if variance > 0.0 { variance.sqrt() } else { 0.0 }; + + (avg_f, error) + } + + /// Runs a single iteration of the VEGAS algorithm using SIMD. + fn run_iteration_simd(&mut self, integrand: &F) -> (f64, f64) { + for grid in &mut self.grids { + grid.reset_importance_data(); + } + + let n_bins = self.grids[0].n_bins(); + let n_packets = self.n_eval / 4; + let n_eval_simd = n_packets * 4; + + let mut ys_soa: Vec = vec![0.0; self.dim * n_eval_simd]; + self.rng.fill(&mut ys_soa[..]); + + let (sum_f, sum_f2, d_updates) = (0..n_packets) .into_par_iter() - .map(|y_vec| { - let mut point = vec![0.0; self.dim]; - let mut bin_indices = vec![0; self.dim]; - let mut d_updates_thread = - (0..self.dim).map(|_| vec![0.0; n_bins]).collect::>(); + .map(|p_idx| { + let mut jacobian_v = f64x4::splat(1.0); + let mut point_v = vec![f64x4::splat(0.0); self.dim]; + let mut bin_indices_arr = [[0; 4]; 32]; - let mut jacobian = 1.0; for d in 0..self.dim { - let y = y_vec[d]; - let (bin_idx, x_unit, jac_vegas) = self.grids[d].map(y); - jacobian *= jac_vegas; - bin_indices[d] = bin_idx; + let offset = d * n_eval_simd + p_idx * 4; + let y_packet = f64x4::new((&ys_soa[offset..offset + 4]).try_into().unwrap()); + + let (x_unit_v, jac_vegas_v, bins_arr) = self.grids[d].map_simd(y_packet); + bin_indices_arr[d] = bins_arr; let (min, max) = self.boundaries[d]; let jac_boundary = max - min; - point[d] = min + x_unit * jac_boundary; - jacobian *= jac_boundary; + + jacobian_v *= jac_vegas_v * f64x4::splat(jac_boundary); + point_v[d] = f64x4::splat(min) + x_unit_v * f64x4::splat(jac_boundary); } - let f_val = integrand.eval(&point); - let weighted_f = f_val * jacobian; - let f2 = weighted_f * weighted_f; + let f_vals_v = integrand.eval_simd(&point_v); + let weighted_f_v = f_vals_v * jacobian_v; - let d_val = f2 / self.n_eval as f64; - for d in 0..self.dim { - d_updates_thread[d][bin_indices[d]] += d_val; + let f_sum = weighted_f_v.reduce_add(); + let f2_sum = (weighted_f_v * weighted_f_v).reduce_add(); + + let mut d_updates_thread = vec![0.0; self.dim * n_bins]; + let d_val_arr = + (weighted_f_v * weighted_f_v / f64x4::splat(n_eval_simd as f64)).to_array(); + + for i in 0..4 { + for d in 0..self.dim { + d_updates_thread[d * n_bins + bin_indices_arr[d][i]] += d_val_arr[i]; + } } - (weighted_f, f2, d_updates_thread) + (f_sum, f2_sum, d_updates_thread) }) .reduce( - || (0.0, 0.0, initial_d_updates.clone()), + || (0.0, 0.0, vec![0.0; self.dim * n_bins]), |mut a, b| { a.0 += b.0; a.1 += b.1; - for d in 0..self.dim { - for i in 0..n_bins { - a.2[d][i] += b.2[d][i]; - } + for i in 0..a.2.len() { + a.2[i] += b.2[i]; } a }, ); for d in 0..self.dim { - self.grids[d].d.copy_from_slice(&d_updates[d]); + let start = d * n_bins; + let end = (d + 1) * n_bins; + self.grids[d].d.copy_from_slice(&d_updates[start..end]); } - let avg_f = sum_f / self.n_eval as f64; - let avg_f2 = sum_f2 / self.n_eval as f64; - let variance = (avg_f2 - avg_f * avg_f) / (self.n_eval - 1).max(1) as f64; + let avg_f = sum_f / n_eval_simd as f64; + let avg_f2 = sum_f2 / n_eval_simd as f64; + let variance = (avg_f2 - avg_f * avg_f) / (n_eval_simd - 1).max(1) as f64; let error = if variance > 0.0 { variance.sqrt() } else { 0.0 }; (avg_f, error) @@ -223,6 +384,7 @@ impl Vegas { pub fn integrate_gpu( &mut self, integrand: &F, + target_accuracy: Option, ) -> VegasResult { assert_eq!(integrand.dim(), self.dim); @@ -253,6 +415,19 @@ impl Vegas { if iter > 0 { iter_results.push(iter_val); iter_errors.push(iter_err); + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + return current_result; + } + } + } + } } // Refine grids on CPU @@ -268,7 +443,8 @@ impl Vegas { #[cfg(test)] mod tests { use super::*; - use crate::integrand::Integrand; + use crate::integrand::{Integrand, SimdIntegrand}; + use wide::f64x4; // Integral of the form exp(-x^2 - y^2) in [-1, 1]^2. struct GaussianIntegrand; @@ -283,6 +459,20 @@ mod tests { } } + struct GaussianSimdIntegrand; + + impl SimdIntegrand for GaussianSimdIntegrand { + fn dim(&self) -> usize { + 2 + } + + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + let x = points[0]; + let y = points[1]; + (-(x * x) - (y * y)).exp() + } + } + const ANALYTICAL_RESULT: f64 = 2.230985; #[test] @@ -290,10 +480,11 @@ mod tests { let integrand = GaussianIntegrand; let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; let mut vegas = Vegas::new(10, 100_000, 50, 0.5, boundaries); - let result = vegas.integrate(&integrand); + vegas.set_seed(1234); + let result = vegas.integrate(&integrand, None); assert!( - (result.value - ANALYTICAL_RESULT).abs() < 1.2 * result.error, + (result.value - ANALYTICAL_RESULT).abs() < 2.5 * result.error, "Analytical={} vs. MCHEP={}+/-{}", ANALYTICAL_RESULT, result.value, @@ -301,4 +492,41 @@ mod tests { ); assert!(result.chi2_dof < 1.5, "chi2_dof: {}", result.chi2_dof); } + + #[test] + fn test_integrate_gaussian_simd() { + let integrand = GaussianSimdIntegrand; + let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; + let mut vegas = Vegas::new(10, 100_000, 50, 0.5, boundaries); + vegas.set_seed(1234); + let result = vegas.integrate_simd(&integrand, None); + + assert!( + (result.value - ANALYTICAL_RESULT).abs() < 2.5 * result.error, + "Analytical={} vs. MCHEP (SIMD)={}+/-{}", + ANALYTICAL_RESULT, + result.value, + result.error + ); + assert!(result.chi2_dof < 1.5, "chi2_dof: {}", result.chi2_dof); + } + + #[test] + fn test_accuracy_goal() { + let integrand = GaussianIntegrand; + let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; + let mut vegas = Vegas::new(10, 100_000, 50, 0.5, boundaries); + vegas.set_seed(1234); + let result = vegas.integrate(&integrand, Some(0.5)); + + let accuracy = (result.error / result.value.abs()) * 100.0; + assert!(accuracy < 0.5); + assert!( + (result.value - ANALYTICAL_RESULT).abs() < 3. * result.error, + "Analytical={} vs. MCHEP={}+/-{}", + ANALYTICAL_RESULT, + result.value, + result.error + ); + } } diff --git a/mchep/src/vegasplus.rs b/mchep/src/vegasplus.rs index cc7af35..1ea2ef3 100644 --- a/mchep/src/vegasplus.rs +++ b/mchep/src/vegasplus.rs @@ -1,11 +1,12 @@ //! The VEGAS+ integrator, which adds adaptive stratified sampling. use crate::grid::Grid; -use crate::integrand::Integrand; +use crate::integrand::{Integrand, SimdIntegrand}; use crate::vegas::VegasResult; use rand::{Rng, SeedableRng}; use rand_pcg::Pcg64; use rayon::prelude::*; +use wide::f64x4; /// Stores the state of a single hypercube for stratified sampling. #[derive(Debug, Clone)] @@ -58,8 +59,6 @@ impl VegasPlus { boundaries: &[(f64, f64)], ) -> Self { let dim = boundaries.len(); - assert!(dim > 0, "Number of dimensions must be positive."); - assert!(n_strat > 0, "Number of stratifications must be positive."); let n_hypercubes = n_strat.pow(dim as u32); assert!( n_eval >= 2 * n_hypercubes, @@ -93,8 +92,17 @@ impl VegasPlus { self.rng = Pcg64::seed_from_u64(seed); } + /// Returns the number of dimensions of the integrator. + pub fn dim(&self) -> usize { + self.dim + } + /// Integrates the given function using the VEGAS+ algorithm. - pub fn integrate(&mut self, integrand: &F) -> VegasResult { + pub fn integrate( + &mut self, + integrand: &F, + target_accuracy: Option, + ) -> VegasResult { assert_eq!(integrand.dim(), self.dim); let mut iter_results = Vec::new(); @@ -106,6 +114,60 @@ impl VegasPlus { if iter > 0 { iter_results.push(iter_val); iter_errors.push(iter_err); + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + return current_result; + } + } + } + } + } + + for grid in &mut self.grids { + grid.refine(); + } + self.reallocate_samples(); + } + + self.combine_results(&iter_results, &iter_errors) + } + + /// Integrates the given function using the VEGAS+ algorithm with SIMD. + pub fn integrate_simd( + &mut self, + integrand: &F, + target_accuracy: Option, + ) -> VegasResult { + assert_eq!(integrand.dim(), self.dim); + + let mut iter_results = Vec::new(); + let mut iter_errors = Vec::new(); + + for iter in 0..self.n_iter { + let (iter_val, iter_err) = self.run_iteration_simd(integrand); + + if iter > 0 { + iter_results.push(iter_val); + iter_errors.push(iter_err); + + if let Some(acc_req) = target_accuracy { + if !iter_results.is_empty() { + let current_result = self.combine_results(&iter_results, &iter_errors); + if current_result.value != 0.0 { + let current_acc = + (current_result.error / current_result.value.abs()) * 100.0; + if current_acc < acc_req { + return current_result; + } + } + } + } } for grid in &mut self.grids { @@ -187,7 +249,7 @@ impl VegasPlus { let avg_f_h = sum_f_h / n_samples_h as f64; let avg_f2_h = sum_f2_h / n_samples_h as f64; let var_g_h = (avg_f2_h - avg_f_h.powi(2)).max(0.0) - * (n_samples_h as f64 / (n_samples_h - 1) as f64); + * (n_samples_h as f64 / (n_samples_h - 1).max(1) as f64); HypercubeResult { sum_f: sum_f_h, @@ -231,6 +293,141 @@ impl VegasPlus { (total_value, error) } + fn run_iteration_simd(&mut self, integrand: &F) -> (f64, f64) { + for grid in &mut self.grids { + grid.reset_importance_data(); + } + + let n_bins = self.grids[0].n_bins(); + + let seeds: Vec = (0..self.hypercubes.len()).map(|_| self.rng.gen()).collect(); + + let results: Vec = seeds + .into_par_iter() + .enumerate() + .map(|(h_idx, seed)| { + let mut seed_array = [0u8; 32]; + seed_array[..8].copy_from_slice(&seed.to_le_bytes()); + let mut thread_rng = Pcg64::from_seed(seed_array); + let n_samples_h = self.hypercubes[h_idx].n_samples; + let n_packets = n_samples_h / 4; + let n_samples_h_simd = n_packets * 4; + + let mut d_updates_thread = + (0..self.dim).map(|_| vec![0.0; n_bins]).collect::>(); + + let mut sum_f_h = 0.0; + let mut sum_f2_h = 0.0; + + if n_samples_h_simd == 0 { + return HypercubeResult { + sum_f: 0.0, + sum_f2: 0.0, + variance: 0.0, + d_updates: d_updates_thread, + }; + } + + let h_coords = self.get_hypercube_coords(h_idx); + let strat_width = 1.0 / self.n_strat as f64; + + // SIMD part + for _ in 0..n_packets { + let mut jacobian_v = f64x4::splat(1.0); + let mut point_v = vec![f64x4::splat(0.0); self.dim]; + let mut bin_indices_arr = vec![[0; 4]; self.dim]; + + let mut y_vec_v = vec![f64x4::splat(0.0); self.dim]; + + for d in 0..self.dim { + let y_rand_v = f64x4::new([ + thread_rng.gen::(), + thread_rng.gen::(), + thread_rng.gen::(), + thread_rng.gen::(), + ]); + y_vec_v[d] = (f64x4::splat(h_coords[d] as f64) + y_rand_v) + * f64x4::splat(strat_width); + } + + for d in 0..self.dim { + let (x_unit_v, jac_vegas_v, bins_arr) = self.grids[d].map_simd(y_vec_v[d]); + jacobian_v *= jac_vegas_v; + bin_indices_arr[d] = bins_arr; + + let (min, max) = self.boundaries[d]; + let jac_boundary = max - min; + point_v[d] = f64x4::splat(min) + x_unit_v * f64x4::splat(jac_boundary); + jacobian_v *= f64x4::splat(jac_boundary); + } + + let f_vals_v = integrand.eval_simd(&point_v); + let weighted_f_v = f_vals_v * jacobian_v; + + sum_f_h += weighted_f_v.reduce_add(); + sum_f2_h += (weighted_f_v * weighted_f_v).reduce_add(); + + let d_val_arr = + (weighted_f_v * weighted_f_v / f64x4::splat(self.n_eval as f64)).to_array(); + + for i in 0..4 { + for d in 0..self.dim { + d_updates_thread[d][bin_indices_arr[d][i]] += d_val_arr[i]; + } + } + } + + let mut var_g_h = 0.0; + if n_samples_h_simd > 1 { + let avg_f_h = sum_f_h / n_samples_h_simd as f64; + let avg_f2_h = sum_f2_h / n_samples_h_simd as f64; + var_g_h = (avg_f2_h - avg_f_h.powi(2)).max(0.0) + * (n_samples_h_simd as f64 / (n_samples_h_simd - 1) as f64); + } + + HypercubeResult { + sum_f: sum_f_h, + sum_f2: sum_f2_h, + variance: var_g_h.sqrt(), + d_updates: d_updates_thread, + } + }) + .collect(); + + let mut total_value = 0.0; + let mut total_variance = 0.0; + let hypercube_volume = (1.0 / self.n_strat as f64).powi(self.dim as i32); + + for (h_idx, result) in results.iter().enumerate() { + let n_samples_h = self.hypercubes[h_idx].n_samples; + let n_samples_h_simd = (n_samples_h / 4) * 4; + + if n_samples_h_simd > 0 { + let avg_g_h = result.sum_f / n_samples_h_simd as f64; + total_value += hypercube_volume * avg_g_h; + } + + if n_samples_h_simd > 1 { + let avg_g_h = result.sum_f / n_samples_h_simd as f64; + let avg_g2_h = result.sum_f2 / n_samples_h_simd as f64; + let var_g_h = (avg_g2_h - avg_g_h.powi(2)).max(0.0) + * (n_samples_h_simd as f64 / (n_samples_h_simd - 1) as f64); + total_variance += hypercube_volume.powi(2) * var_g_h / n_samples_h_simd as f64; + } + + self.hypercubes[h_idx].variance = result.variance; + for d in 0..self.dim { + for i in 0..n_bins { + self.grids[d].d[i] += result.d_updates[d][i]; + } + } + } + + let error = total_variance.sqrt(); + + (total_value, error) + } + /// Converts a hypercube index to its coordinates in the stratified grid. fn get_hypercube_coords(&self, index: usize) -> Vec { let mut coords = vec![0; self.dim]; @@ -327,7 +524,8 @@ impl VegasPlus { #[cfg(test)] mod tests { use super::*; - use crate::integrand::Integrand; + use crate::integrand::{Integrand, SimdIntegrand}; + use wide::f64x4; // Integral of the form exp(-x^2 - y^2) in [-1, 1]^2. struct GaussianIntegrand; @@ -346,17 +544,67 @@ mod tests { #[test] fn test_integrate_gaussian_plus() { + let integrand = GaussianIntegrand; + let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; + let mut vegas_plus = VegasPlus::new(20, 200_000, 50, 0.5, 4, 0.75, boundaries); + let result = vegas_plus.integrate(&integrand, None); + assert!( + (result.value - ANALYTICAL_RESULT).abs() < 2.5 * result.error, + "Analytical={} vs. MCHEP={}+/-{}", + ANALYTICAL_RESULT, + result.value, + result.error + ); + assert!(result.chi2_dof < 1.5, "chi2_dof: {}", result.chi2_dof); + } + + #[test] + fn test_accuracy_goal_plus() { let integrand = GaussianIntegrand; let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; let mut vegas_plus = VegasPlus::new(10, 200_000, 50, 0.5, 4, 0.75, boundaries); - let result = vegas_plus.integrate(&integrand); + vegas_plus.set_seed(4321); + let result = vegas_plus.integrate(&integrand, Some(0.1)); + + let accuracy = (result.error / result.value.abs()) * 100.0; + assert!(accuracy < 0.1); assert!( - (result.value - ANALYTICAL_RESULT).abs() < 1.2 * result.error, + (result.value - ANALYTICAL_RESULT).abs() < 3. * result.error, "Analytical={} vs. MCHEP={}+/-{}", ANALYTICAL_RESULT, result.value, result.error ); + } + + struct GaussianSimdIntegrand; + + impl SimdIntegrand for GaussianSimdIntegrand { + fn dim(&self) -> usize { + 2 + } + + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + let x = points[0]; + let y = points[1]; + (-(x * x) - (y * y)).exp() + } + } + + #[test] + fn test_integrate_gaussian_plus_simd() { + let integrand = GaussianSimdIntegrand; + let boundaries = &[(-1.0, 1.0), (-1.0, 1.0)]; + let mut vegas_plus = VegasPlus::new(20, 200_000, 50, 0.5, 4, 0.75, boundaries); + vegas_plus.set_seed(1234); + let result = vegas_plus.integrate_simd(&integrand, None); + assert!( + (result.value - ANALYTICAL_RESULT).abs() < 2.5 * result.error, + "Analytical={} vs. MCHEP (SIMD)={}+/-{}", + ANALYTICAL_RESULT, + result.value, + result.error + ); assert!(result.chi2_dof < 1.5, "chi2_dof: {}", result.chi2_dof); } } diff --git a/mchep_capi/Cargo.toml b/mchep_capi/Cargo.toml index af66672..a2a91d8 100644 --- a/mchep_capi/Cargo.toml +++ b/mchep_capi/Cargo.toml @@ -3,14 +3,24 @@ name = "mchep_capi" authors = ["Tanjona R. Rabemananjara "] description = "C language interface to NeoPDF" readme = "README.md" -version = "0.1.0" -edition = "2021" +build = "build.rs" + +categories.workspace = true +edition.workspace = true +keywords.workspace = true +license.workspace = true +repository.workspace = true +version.workspace = true [dependencies] mchep = { path = "../mchep" } +wide = { workspace = true } [build-dependencies] cbindgen = "0.26.0" [features] capi = [] + +[lints] +workspace = true diff --git a/mchep_capi/bench/Makefile b/mchep_capi/bench/Makefile new file mode 100644 index 0000000..8792735 --- /dev/null +++ b/mchep_capi/bench/Makefile @@ -0,0 +1,41 @@ +CXX = c++ -std=c++11 +CXXFLAGS = -O3 -Wall -Wextra +MCHEP_DEPS != pkg-config --cflags --libs mchep_capi +MATH_LIBS = -lm + +MCHEP_TARGETS = benchmark_mchep benchmark_mchep_simd +CUBA_TARGET = cuba_example +TARGETS = $(MCHEP_TARGETS) $(CUBA_TARGET) + +all: $(TARGETS) + +# --- MCHEP Benchmarks --- +benchmark_mchep: benchmark_mchep.cpp + $(CXX) $(CXXFLAGS) $< $(MCHEP_DEPS) $(MATH_LIBS) -o $@ + +benchmark_mchep_simd: benchmark_mchep_simd.cpp + $(CXX) $(CXXFLAGS) $< $(MCHEP_DEPS) $(MATH_LIBS) -o $@ + +# --- LIBCUBA Benchmark --- +ifndef CONDA_PREFIX +$(error CONDA_PREFIX is not set!) +endif + +CUBA_CXXFLAGS = -Wall -O2 -std=c++11 +CUBA_INCLUDES = -I$(CONDA_PREFIX)/include +CUBA_LIBDIRS = -L$(CONDA_PREFIX)/lib +CUBA_LIBS = -lcuba -lm + +cuba_example: cuba_example.cpp + @echo "Building libcuba example..." + $(CXX) $(CUBA_CXXFLAGS) $(CUBA_INCLUDES) $(CUBA_LIBDIRS) $< $(CUBA_LIBS) -o $@ + +run_cuba: cuba_example + @echo "Running cuba_example..." + @export LD_LIBRARY_PATH=$(CONDA_PREFIX)/lib:$$LD_LIBRARY_PATH && ./cuba_example + +.PHONY: all clean run_cuba + +clean: + @echo "Cleaning benchmarks..." + rm -f $(TARGETS) *.o diff --git a/mchep_capi/bench/benchmark.sh b/mchep_capi/bench/benchmark.sh new file mode 100755 index 0000000..e07135e --- /dev/null +++ b/mchep_capi/bench/benchmark.sh @@ -0,0 +1,78 @@ +#!/bin/bash + +# Exit on error +set -e + +# Number of runs for averaging +N_RUNS=5 + +echo "Starting comprehensive benchmark..." +echo "===================================" + +# --- MCHEP Scalar Benchmark --- +echo "" +echo "--- Benchmarking MCHEP (Scalar, Optimized) ---" + +# Build benchmark +echo "Building mchep scalar benchmark..." +make benchmark_mchep + +# Run benchmark +echo "Running mchep scalar benchmark ($N_RUNS runs)..." +MCHEP_SCALAR_TIMES=() +for i in $(seq 1 $N_RUNS); do + TIME_OUTPUT=$( (time -p ./benchmark_mchep) 2>&1 ) + REAL_TIME=$(echo "$TIME_OUTPUT" | grep real | awk '{print $2}' | tr ',' '.') + MCHEP_SCALAR_TIMES+=($REAL_TIME) + echo "Run $i: $REAL_TIME s" +done +MCHEP_SCALAR_AVG_TIME=$(echo "${MCHEP_SCALAR_TIMES[@]}" | awk '{for(i=1;i<=NF;i++)s+=$i; print s/NF}') +echo "Average MCHEP (Scalar) time: $MCHEP_SCALAR_AVG_TIME s" + +# --- MCHEP SIMD Benchmark --- +echo "" +echo "--- Benchmarking MCHEP (SIMD) ---" + +# Build benchmark +echo "Building mchep SIMD benchmark..." +make benchmark_mchep_simd + +# Run benchmark +echo "Running mchep SIMD benchmark ($N_RUNS runs)..." +MCHEP_SIMD_TIMES=() +for i in $(seq 1 $N_RUNS); do + TIME_OUTPUT=$( (time -p ./benchmark_mchep_simd) 2>&1 ) + REAL_TIME=$(echo "$TIME_OUTPUT" | grep real | awk '{print $2}' | tr ',' '.') + MCHEP_SIMD_TIMES+=($REAL_TIME) + echo "Run $i: $REAL_TIME s" +done +MCHEP_SIMD_AVG_TIME=$(echo "${MCHEP_SIMD_TIMES[@]}" | awk '{for(i=1;i<=NF;i++)s+=$i; print s/NF}') +echo "Average MCHEP (SIMD) time: $MCHEP_SIMD_AVG_TIME s" + +# --- LIBCUBA Benchmark --- +echo "" +echo "--- Benchmarking LIBCUBA ---" + +# Build benchmark +echo "Building libcuba benchmark..." +make cuba_example + +# Run benchmark +echo "Running libcuba benchmark ($N_RUNS runs)..." +LIBCUBA_TIMES=() +for i in $(seq 1 $N_RUNS); do + TIME_OUTPUT=$( (time -p make run_cuba) 2>&1 ) + REAL_TIME=$(echo "$TIME_OUTPUT" | grep real | awk '{print $2}' | tr ',' '.') + LIBCUBA_TIMES+=($REAL_TIME) + echo "Run $i: $REAL_TIME s" +done +LIBCUBA_AVG_TIME=$(echo "${LIBCUBA_TIMES[@]}" | awk '{for(i=1;i<=NF;i++)s+=$i; print s/NF}') +echo "Average LIBCUBA time: $LIBCUBA_AVG_TIME s" + +# --- Summary --- +echo "" +echo "--- Benchmark Summary ---" +echo "MCHEP (Scalar, Optimized): $MCHEP_SCALAR_AVG_TIME s" +echo "MCHEP (SIMD): $MCHEP_SIMD_AVG_TIME s" +echo "LIBCUBA (Vegas): $LIBCUBA_AVG_TIME s" +echo "=======================" diff --git a/mchep_capi/bench/benchmark_mchep.cpp b/mchep_capi/bench/benchmark_mchep.cpp new file mode 100644 index 0000000..424cd99 --- /dev/null +++ b/mchep_capi/bench/benchmark_mchep.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include + +// Integrand function from cuba_example.cpp +// f(x) = exp(-100 * sum((x[i] - 0.5)^2)) * 1013.2118364296088 +double integrand_4d(const std::vector& x) { + double dx2 = 0.0; + for (int d = 0; d < 4; d++) { + double diff = x[d] - 0.5; + dx2 += diff * diff; + } + return std::exp(-100.0 * dx2) * 1013.2118364296088; +} + +int main() { + const size_t n_iter = 10; + const size_t n_eval = 100000; + const size_t n_bins = 50; + const double alpha = 1.5; + + std::vector> boundaries = { + {0.0, 1.0}, + {0.0, 1.0}, + {0.0, 1.0}, + {0.0, 1.0} + }; + + try { + mchep::Vegas vegas(n_iter, n_eval, n_bins, alpha, boundaries); + vegas.set_seed(1234); + VegasResult result = vegas.integrate(integrand_4d, -1.0); + + std::cout << "MCHEP Result: " << result.value << " +/- " << result.error << std::endl; + + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; + } + + return 0; +} diff --git a/mchep_capi/bench/benchmark_mchep_simd.cpp b/mchep_capi/bench/benchmark_mchep_simd.cpp new file mode 100644 index 0000000..57c7319 --- /dev/null +++ b/mchep_capi/bench/benchmark_mchep_simd.cpp @@ -0,0 +1,53 @@ +#include +#include +#include +#include +#include +#include + +// Integrand function from cuba_example.cpp, adapted for SIMD +// f(x) = exp(-100 * sum((x[i] - 0.5)^2)) * 1013.2118364296088 +std::array integrand_4d_simd(const std::vector& x) { + std::array results; + const int dim = 4; + const int simd_width = 4; + + for (int i = 0; i < simd_width; ++i) { + double dx2 = 0.0; + for (int d = 0; d < dim; ++d) { + double val = x[d * simd_width + i]; + double diff = val - 0.5; + dx2 += diff * diff; + } + results[i] = std::exp(-100.0 * dx2) * 1013.2118364296088; + } + return results; +} + +int main() { + const size_t n_iter = 10; + const size_t n_eval = 100000; + const size_t n_bins = 50; + const double alpha = 1.5; + + std::vector> boundaries = { + {0.0, 1.0}, + {0.0, 1.0}, + {0.0, 1.0}, + {0.0, 1.0} + }; + + try { + mchep::Vegas vegas(n_iter, n_eval, n_bins, alpha, boundaries); + vegas.set_seed(1234); + VegasResult result = vegas.integrate_simd(integrand_4d_simd, -1.0); + + std::cout << "MCHEP SIMD Result: " << result.value << " +/- " << result.error << std::endl; + + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; + } + + return 0; +} diff --git a/mchep_capi/bench/cuba_example.cpp b/mchep_capi/bench/cuba_example.cpp new file mode 100644 index 0000000..43088a2 --- /dev/null +++ b/mchep_capi/bench/cuba_example.cpp @@ -0,0 +1,60 @@ +#include +#include +#include + +// Integrand function: f(x) = exp(-100 * sum((x[i] - 0.5)^2)) * 1013.2118364296088 +// This integrates over [0,1]^4 (4-dimensional hypercube) +static int integrand(const int *ndim, const cubareal x[], + const int *ncomp, cubareal f[], void *userdata) { + double dx2 = 0.0; + for (int d = 0; d < 4; d++) { + double diff = x[d] - 0.5; + dx2 += diff * diff; + } + + f[0] = std::exp(-100.0 * dx2) * 1013.2118364296088; + + return 0; +} + +int main() { + int neval, fail; + cubareal integral[1], error[1], prob[1]; + + std::cout << "CUBA Library Integration Example (C++)" << std::endl; + std::cout << "=======================================" << std::endl << std::endl; + std::cout << "Integrating: f(x) = exp(-100 * sum((x[i]-0.5)^2)) * 1013.2118364296088" << std::endl; + std::cout << "Over: [0,1]^4 (4-dimensional unit hypercube)" << std::endl << std::endl; + + std::cout << "--- Method 1: Vegas ---" << std::endl; + Vegas(4, + 1, + integrand, + nullptr, + 1, + 0.0, + 0.0, + 0, + 0, + 1000, + 1000000, + 10000, + 1000, + 1000, + 0, + nullptr, + nullptr, + &neval, + &fail, + integral, + error, + prob); + + std::cout << "Result: " << integral[0] << std::endl; + std::cout << "Error: " << error[0] << std::endl; + std::cout << "Chi-sq prob: " << prob[0] << std::endl; + std::cout << "Evaluations: " << neval << std::endl; + std::cout << "Status: " << (fail == 0 ? "Success" : "Failed") << std::endl << std::endl; + + return 0; +} diff --git a/mchep_capi/build.rs b/mchep_capi/build.rs new file mode 100644 index 0000000..2cdb4ba --- /dev/null +++ b/mchep_capi/build.rs @@ -0,0 +1,21 @@ +//! A build script to install the OOP C++ interface to `MCHEP` + +use std::env; +use std::fs; +use std::path::PathBuf; + +fn main() { + println!("cargo:rerun-if-changed=src/include/mchep.hpp"); + + if let Ok(prefix) = env::var("CARGO_C_MCHEP_INSTALL_PREFIX") { + let prefix_path = PathBuf::from(prefix); + let include_path = prefix_path.join("include").join("mchep_capi"); + + fs::create_dir_all(&include_path).expect("Failed to create include directory."); + + let source_header = PathBuf::from("src/include/mchep.hpp"); + let dest_header = include_path.join("mchep.hpp"); + + fs::copy(&source_header, &dest_header).expect("Failed to copy header file."); + } +} diff --git a/mchep_capi/src/include/mchep.hpp b/mchep_capi/src/include/mchep.hpp new file mode 100644 index 0000000..03cbfa2 --- /dev/null +++ b/mchep_capi/src/include/mchep.hpp @@ -0,0 +1,215 @@ +#ifndef MCHEP_HPP +#define MCHEP_HPP + +#include +#include +#include +#include +#include +#include +#include + +extern "C" { +#include +} + +namespace mchep { + +// Forward declaration +class Vegas; +class VegasPlus; + +namespace internal { +// This wrapper function is the bridge between C-style function pointers and +// std::function. It will be passed to the C API. +extern "C" double integrand_wrapper(const double *x, int dim, void *user_data) { + if (user_data == nullptr) { + // Or handle error appropriately + return 0.0; + } + auto *func = + static_cast &)> *>( + user_data); + std::vector x_vec(x, x + dim); + return (*func)(x_vec); +} + +// SIMD wrapper +extern "C" void integrand_simd_wrapper(const double *x, int dim, + void *user_data, double *result) { + if (user_data == nullptr) { + return; + } + auto *func = static_cast(const std::vector &)> *>(user_data); + std::vector x_vec(x, x + dim * 4); + std::array res_arr = (*func)(x_vec); + std::copy(res_arr.begin(), res_arr.end(), result); +} +} // namespace internal + +class Vegas { +public: + /// @brief Constructor for the Vegas integrator. + /// @param n_iter Number of iterations. + /// @param n_eval Number of function evaluations per iteration. + /// @param n_bins Number of bins for the adaptive grid. + /// @param alpha The grid adaptation parameter. + /// @param boundaries The integration boundaries for each dimension. + Vegas(size_t n_iter, size_t n_eval, size_t n_bins, double alpha, + const std::vector> &boundaries) { + std::vector c_boundaries; + c_boundaries.reserve(boundaries.size()); + for (const auto &b : boundaries) { + c_boundaries.push_back({b.first, b.second}); + } + + vegas_ptr_ = mchep_vegas_new(n_iter, n_eval, n_bins, alpha, + c_boundaries.size(), c_boundaries.data()); + if (vegas_ptr_ == nullptr) { + throw std::runtime_error("Failed to create MCHEP Vegas integrator."); + } + } + + /// @brief Destructor. + ~Vegas() { mchep_vegas_free(vegas_ptr_); } + + // Delete copy constructor and copy assignment operator + Vegas(const Vegas &) = delete; + Vegas &operator=(const Vegas &) = delete; + + /// @brief Move constructor. + Vegas(Vegas &&other) noexcept : vegas_ptr_(other.vegas_ptr_) { + other.vegas_ptr_ = nullptr; + } + + /// @brief Move assignment operator. + Vegas &operator=(Vegas &&other) noexcept { + if (this != &other) { + mchep_vegas_free(vegas_ptr_); + vegas_ptr_ = other.vegas_ptr_; + other.vegas_ptr_ = nullptr; + } + return *this; + } + + /// @brief Sets the seed for the random number generator. + /// @param seed The seed to use. + void set_seed(uint64_t seed) { mchep_vegas_set_seed(vegas_ptr_, seed); } + + /// @brief Integrates the given function. + /// @param integrand The function to integrate. It should take a vector of + /// doubles and return a double. + /// @return The integration result. + VegasResult + integrate(std::function &)> integrand, + double target_accuracy = -1.0) { + return mchep_vegas_integrate(vegas_ptr_, internal::integrand_wrapper, + &integrand, target_accuracy); + } + + /// @brief Integrates the given function using SIMD. + /// @param integrand The function to integrate. It should take a vector of + /// doubles of size dim*4 (SoA) and return an array of 4 doubles. + /// @param target_accuracy The desired accuracy in percent. If non-positive, it is ignored. + /// @return The integration result. + VegasResult integrate_simd( + std::function(const std::vector &)> + integrand, + double target_accuracy = -1.0) { + return mchep_vegas_integrate_simd(vegas_ptr_, + internal::integrand_simd_wrapper, + &integrand, target_accuracy); + } + +private: + VegasC *vegas_ptr_; +}; + +class VegasPlus { +public: + /// @brief Constructor for the Vegas+ integrator. + /// @param n_iter Number of iterations. + /// @param n_eval Number of function evaluations per iteration. + /// @param n_bins Number of bins for the adaptive grid. + /// @param alpha The grid adaptation parameter. + /// @param n_strat Number of stratifications per dimension. + /// @param beta The stratified sampling adaptation parameter. + /// @param boundaries The integration boundaries for each dimension. + VegasPlus(size_t n_iter, size_t n_eval, size_t n_bins, double alpha, + size_t n_strat, double beta, + const std::vector> &boundaries) { + std::vector c_boundaries; + c_boundaries.reserve(boundaries.size()); + for (const auto &b : boundaries) { + c_boundaries.push_back({b.first, b.second}); + } + + vegas_plus_ptr_ = mchep_vegas_plus_new( + n_iter, n_eval, n_bins, alpha, n_strat, beta, c_boundaries.size(), + c_boundaries.data()); + if (vegas_plus_ptr_ == nullptr) { + throw std::runtime_error("Failed to create MCHEP VegasPlus integrator."); + } + } + + /// @brief Destructor. + ~VegasPlus() { mchep_vegas_plus_free(vegas_plus_ptr_); } + + // Delete copy constructor and copy assignment operator + VegasPlus(const VegasPlus &) = delete; + VegasPlus &operator=(const VegasPlus &) = delete; + + /// @brief Move constructor. + VegasPlus(VegasPlus &&other) noexcept + : vegas_plus_ptr_(other.vegas_plus_ptr_) { + other.vegas_plus_ptr_ = nullptr; + } + + /// @brief Move assignment operator. + VegasPlus &operator=(VegasPlus &&other) noexcept { + if (this != &other) { + mchep_vegas_plus_free(vegas_plus_ptr_); + vegas_plus_ptr_ = other.vegas_plus_ptr_; + other.vegas_plus_ptr_ = nullptr; + } + return *this; + } + + /// @brief Sets the seed for the random number generator. + /// @param seed The seed to use. + void set_seed(uint64_t seed) { + mchep_vegas_plus_set_seed(vegas_plus_ptr_, seed); + } + + /// @brief Integrates the given function. + /// @param integrand The function to integrate. + /// @return The integration result. + VegasResult + integrate(std::function &)> integrand, + double target_accuracy = -1.0) { + return mchep_vegas_plus_integrate(vegas_plus_ptr_, + internal::integrand_wrapper, &integrand, + target_accuracy); + } + + /// @brief Integrates the given function using SIMD. + /// @param integrand The function to integrate. + /// @param target_accuracy The desired accuracy in percent. + /// @return The integration result. + VegasResult integrate_simd( + std::function(const std::vector &)> + integrand, + double target_accuracy = -1.0) { + return mchep_vegas_plus_integrate_simd( + vegas_plus_ptr_, internal::integrand_simd_wrapper, &integrand, + target_accuracy); + } + +private: + VegasPlusC *vegas_plus_ptr_; +}; + +} // namespace mchep + +#endif // MCHEP_HPP diff --git a/mchep_capi/src/lib.rs b/mchep_capi/src/lib.rs index 8b4c82c..727316b 100644 --- a/mchep_capi/src/lib.rs +++ b/mchep_capi/src/lib.rs @@ -1,15 +1,19 @@ //! The C-language interface for `MCHEP` +use mchep::integrand::{Integrand, SimdIntegrand}; +use mchep::vegas::{Vegas, VegasResult}; +use mchep::vegasplus::VegasPlus; use std::ffi::c_void; use std::os::raw::c_int; use std::slice; - -use mchep::vegas::{Vegas, VegasResult}; +use wide::f64x4; /// A C-compatible struct for integration boundaries. #[repr(C)] pub struct CBoundary { + /// The minimum of the boundary. pub min: f64, + /// The maximum of the boundary. pub max: f64, } @@ -19,6 +23,13 @@ pub struct CBoundary { /// The third is a user-provided `user_data` pointer. pub type CIntegrand = extern "C" fn(*const f64, c_int, *mut c_void) -> f64; +/// The C-style SIMD integrand function pointer. +/// The first argument is the point `x` (an array of f64 in SoA layout, size dim * 4). +/// The second argument is the dimension. +/// The third is a user-provided `user_data` pointer. +/// The fourth is the result array (array of 4 f64). +pub type CSimdIntegrand = extern "C" fn(*const f64, c_int, *mut c_void, *mut f64); + /// A wrapper that implements the Rust `Integrand` trait. struct CIntegrandWrapper { dim: usize, @@ -26,7 +37,7 @@ struct CIntegrandWrapper { user_data: *mut c_void, } -impl mchep::integrand::Integrand for CIntegrandWrapper { +impl Integrand for CIntegrandWrapper { fn dim(&self) -> usize { self.dim } @@ -36,10 +47,36 @@ impl mchep::integrand::Integrand for CIntegrandWrapper { } } +/// A wrapper that implements the Rust `SimdIntegrand` trait. +struct CSimdIntegrandWrapper { + dim: usize, + func: CSimdIntegrand, + user_data: *mut c_void, +} + +impl SimdIntegrand for CSimdIntegrandWrapper { + fn dim(&self) -> usize { + self.dim + } + + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + let points_ptr = points.as_ptr() as *const f64; + let mut result_arr = [0.0f64; 4]; + (self.func)( + points_ptr, + self.dim as c_int, + self.user_data, + result_arr.as_mut_ptr(), + ); + f64x4::from(result_arr) + } +} + /// This is unsafe, but required to integrate with Rayon. /// The user of the C API is responsible for ensuring that the provided /// integrand function is thread-safe. unsafe impl Sync for CIntegrandWrapper {} +unsafe impl Sync for CSimdIntegrandWrapper {} /// The opaque pointer to the Vegas integrator. pub type VegasC = c_void; @@ -59,7 +96,7 @@ pub unsafe extern "C" fn mchep_vegas_new( dim: usize, boundaries: *const CBoundary, ) -> *mut VegasC { - let boundaries_slice = slice::from_raw_parts(boundaries, dim); + let boundaries_slice = unsafe { slice::from_raw_parts(boundaries, dim) }; let rust_boundaries: Vec<(f64, f64)> = boundaries_slice.iter().map(|b| (b.min, b.max)).collect(); @@ -78,8 +115,9 @@ pub unsafe extern "C" fn mchep_vegas_integrate( vegas_ptr: *mut VegasC, integrand_func: CIntegrand, user_data: *mut c_void, + target_accuracy: f64, ) -> VegasResult { - let vegas = &mut *(vegas_ptr as *mut Vegas); + let vegas = unsafe { &mut *(vegas_ptr as *mut Vegas) }; let integrand = CIntegrandWrapper { dim: vegas.dim(), @@ -87,7 +125,43 @@ pub unsafe extern "C" fn mchep_vegas_integrate( user_data, }; - vegas.integrate(&integrand) + let accuracy_opt = if target_accuracy > 0.0 { + Some(target_accuracy) + } else { + None + }; + + vegas.integrate(&integrand, accuracy_opt) +} + +/// Integrates the given function using the VEGAS algorithm with SIMD. +/// +/// # Safety +/// +/// `vegas_ptr` must be a valid pointer returned by `mchep_vegas_new`. +/// `integrand_func` must be a valid function pointer. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_integrate_simd( + vegas_ptr: *mut VegasC, + integrand_func: CSimdIntegrand, + user_data: *mut c_void, + target_accuracy: f64, +) -> VegasResult { + let vegas = unsafe { &mut *(vegas_ptr as *mut Vegas) }; + + let integrand = CSimdIntegrandWrapper { + dim: vegas.dim(), + func: integrand_func, + user_data, + }; + + let accuracy_opt = if target_accuracy > 0.0 { + Some(target_accuracy) + } else { + None + }; + + vegas.integrate_simd(&integrand, accuracy_opt) } /// Sets the seed for the random number generator. @@ -97,7 +171,7 @@ pub unsafe extern "C" fn mchep_vegas_integrate( /// `vegas_ptr` must be a valid pointer returned by `mchep_vegas_new`. #[no_mangle] pub unsafe extern "C" fn mchep_vegas_set_seed(vegas_ptr: *mut VegasC, seed: u64) { - let vegas = &mut *(vegas_ptr as *mut Vegas); + let vegas = unsafe { &mut *(vegas_ptr as *mut Vegas) }; vegas.set_seed(seed); } @@ -110,6 +184,126 @@ pub unsafe extern "C" fn mchep_vegas_set_seed(vegas_ptr: *mut VegasC, seed: u64) #[no_mangle] pub unsafe extern "C" fn mchep_vegas_free(vegas_ptr: *mut VegasC) { if !vegas_ptr.is_null() { - drop(Box::from_raw(vegas_ptr as *mut Vegas)); + drop(unsafe { Box::from_raw(vegas_ptr as *mut Vegas) }); + } +} + +/// The opaque pointer to the VegasPlus integrator. +pub type VegasPlusC = c_void; + +/// Creates a new VEGAS+ integrator. +/// +/// # Safety +/// +/// `boundaries` must be a valid pointer to an array of `CBoundary` of size `dim`. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_plus_new( + n_iter: usize, + n_eval: usize, + n_bins: usize, + alpha: f64, + n_strat: usize, + beta: f64, + dim: usize, + boundaries: *const CBoundary, +) -> *mut VegasPlusC { + let boundaries_slice = unsafe { slice::from_raw_parts(boundaries, dim) }; + let rust_boundaries: Vec<(f64, f64)> = + boundaries_slice.iter().map(|b| (b.min, b.max)).collect(); + + let vegas_plus = VegasPlus::new( + n_iter, + n_eval, + n_bins, + alpha, + n_strat, + beta, + &rust_boundaries, + ); + let b = Box::new(vegas_plus); + Box::into_raw(b) as *mut VegasPlusC +} + +/// Integrates the given function using the VEGAS+ algorithm. +/// +/// # Safety +/// +/// `vegas_plus_ptr` must be a valid pointer returned by `mchep_vegas_plus_new`. +/// `integrand_func` must be a valid function pointer. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_plus_integrate( + vegas_plus_ptr: *mut VegasPlusC, + integrand_func: CIntegrand, + user_data: *mut c_void, + target_accuracy: f64, +) -> VegasResult { + let vegas_plus = unsafe { &mut *(vegas_plus_ptr as *mut VegasPlus) }; + + let integrand = CIntegrandWrapper { + dim: vegas_plus.dim(), + func: integrand_func, + user_data, + }; + + let accuracy_opt = if target_accuracy > 0.0 { + Some(target_accuracy) + } else { + None + }; + + vegas_plus.integrate(&integrand, accuracy_opt) +} + +/// Integrates the given function using the VEGAS+ algorithm with SIMD. +/// +/// # Safety +/// +/// `vegas_plus_ptr` must be a valid pointer returned by `mchep_vegas_plus_new`. +/// `integrand_func` must be a valid function pointer. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_plus_integrate_simd( + vegas_plus_ptr: *mut VegasPlusC, + integrand_func: CSimdIntegrand, + user_data: *mut c_void, + target_accuracy: f64, +) -> VegasResult { + let vegas_plus = unsafe { &mut *(vegas_plus_ptr as *mut VegasPlus) }; + + let integrand = CSimdIntegrandWrapper { + dim: vegas_plus.dim(), + func: integrand_func, + user_data, + }; + + let accuracy_opt = if target_accuracy > 0.0 { + Some(target_accuracy) + } else { + None + }; + + vegas_plus.integrate_simd(&integrand, accuracy_opt) +} + +/// Sets the seed for the random number generator for VEGAS+. +/// +/// # Safety +/// +/// `vegas_plus_ptr` must be a valid pointer returned by `mchep_vegas_plus_new`. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_plus_set_seed(vegas_plus_ptr: *mut VegasPlusC, seed: u64) { + let vegas_plus = unsafe { &mut *(vegas_plus_ptr as *mut VegasPlus) }; + vegas_plus.set_seed(seed); +} + +/// Frees the memory of the VEGAS+ integrator. +/// +/// # Safety +/// +/// `vegas_plus_ptr` must be a valid pointer returned by `mchep_vegas_plus_new` +/// and must not be used afterward. +#[no_mangle] +pub unsafe extern "C" fn mchep_vegas_plus_free(vegas_plus_ptr: *mut VegasPlusC) { + if !vegas_plus_ptr.is_null() { + drop(unsafe { Box::from_raw(vegas_plus_ptr as *mut VegasPlus) }); } } diff --git a/mchep_capi/tests/Makefile b/mchep_capi/tests/Makefile index 1b57b71..6faaaff 100644 --- a/mchep_capi/tests/Makefile +++ b/mchep_capi/tests/Makefile @@ -5,7 +5,7 @@ CXXFLAGS = -O3 -Wall -Wextra -Werror MCHEP_DEPS != pkg-config --cflags --libs mchep_capi MATH_LIBS = -lm -PROGRAMS = test_capi test_cppapi +PROGRAMS = test_capi test_cppapi test_cppapi_simd all: $(PROGRAMS) @@ -18,6 +18,9 @@ test_capi: test_capi.c test_cppapi: test_cppapi.cpp $(CXX) $(CXXFLAGS) $< $(MCHEP_DEPS) $(MATH_LIBS) -o $@ +test_cppapi_simd: test_cppapi_simd.cpp + $(CXX) $(CXXFLAGS) $< $(MCHEP_DEPS) $(MATH_LIBS) -o $@ + .PHONY: clean clean: diff --git a/mchep_capi/tests/test_capi.c b/mchep_capi/tests/test_capi.c index 4fa5806..5c1d499 100644 --- a/mchep_capi/tests/test_capi.c +++ b/mchep_capi/tests/test_capi.c @@ -31,7 +31,7 @@ int main() { return 1; } - struct VegasResult result = mchep_vegas_integrate(vegas, gaussian, NULL); + struct VegasResult result = mchep_vegas_integrate(vegas, gaussian, NULL, -1.0); mchep_vegas_free(vegas); const double expected = 2.230985; @@ -47,13 +47,13 @@ int main() { VegasC *vegas1 = mchep_vegas_new(n_iter, n_eval, n_bins, alpha, dim, boundaries); mchep_vegas_set_seed(vegas1, 1234); - struct VegasResult result1 = mchep_vegas_integrate(vegas1, gaussian, NULL); + struct VegasResult result1 = mchep_vegas_integrate(vegas1, gaussian, NULL, -1.0); mchep_vegas_free(vegas1); VegasC *vegas2 = mchep_vegas_new(n_iter, n_eval, n_bins, alpha, dim, boundaries); mchep_vegas_set_seed(vegas2, 1234); - struct VegasResult result2 = mchep_vegas_integrate(vegas2, gaussian, NULL); + struct VegasResult result2 = mchep_vegas_integrate(vegas2, gaussian, NULL, -1.0); mchep_vegas_free(vegas2); printf("Result 1: %f +/- %f\n", result1.value, result1.error); @@ -66,12 +66,23 @@ int main() { VegasC *vegas3 = mchep_vegas_new(n_iter, n_eval, n_bins, alpha, dim, boundaries); mchep_vegas_set_seed(vegas3, 5678); - struct VegasResult result3 = mchep_vegas_integrate(vegas3, gaussian, NULL); + struct VegasResult result3 = mchep_vegas_integrate(vegas3, gaussian, NULL, -1.0); mchep_vegas_free(vegas3); printf("Result 3: %f +/- %f\n", result3.value, result3.error); assert(result1.value != result3.value); printf("Different seed test passed.\n"); + printf("\nTesting C API accuracy goal...\n"); + VegasC *vegas4 = mchep_vegas_new(20, 1e8, n_bins, alpha, dim, boundaries); + mchep_vegas_set_seed(vegas4, 1234); + struct VegasResult result4 = mchep_vegas_integrate(vegas4, gaussian, NULL, 1e-4); + mchep_vegas_free(vegas4); + + double accuracy = (result4.error / fabs(result4.value)) * 100.0; + printf("Accuracy goal test: value=%f, error=%f, acc=%f\n", result4.value, result4.error, accuracy); + assert(accuracy < 1e-4); + printf("Accuracy goal test passed.\n"); + return 0; } diff --git a/mchep_capi/tests/test_capi.output b/mchep_capi/tests/test_capi.output index 78f3009..2cb8677 100644 --- a/mchep_capi/tests/test_capi.output +++ b/mchep_capi/tests/test_capi.output @@ -7,3 +7,7 @@ Result 2: 2.230928 +/- 0.000106 Same seed test passed. Result 3: 2.230882 +/- 0.000106 Different seed test passed. + +Testing C API accuracy goal... +Accuracy goal test: value=2.230987, error=0.000002, acc=0.000098 +Accuracy goal test passed. diff --git a/mchep_capi/tests/test_cppapi.cpp b/mchep_capi/tests/test_cppapi.cpp index c3b63e2..4a85059 100644 --- a/mchep_capi/tests/test_cppapi.cpp +++ b/mchep_capi/tests/test_cppapi.cpp @@ -26,7 +26,7 @@ int main() { try { mchep::Vegas vegas(n_iter, n_eval, n_bins, alpha, boundaries); - VegasResult result = vegas.integrate(gaussian_cpp); + VegasResult result = vegas.integrate(gaussian_cpp, -1.0); const double expected = 2.230985; const double multiplier = 2.5; @@ -41,11 +41,11 @@ int main() { // Test 1: Same seed should produce same result mchep::Vegas vegas1(n_iter, n_eval, n_bins, alpha, boundaries); vegas1.set_seed(1234); - VegasResult result1 = vegas1.integrate(gaussian_cpp); + VegasResult result1 = vegas1.integrate(gaussian_cpp, -1.0); mchep::Vegas vegas2(n_iter, n_eval, n_bins, alpha, boundaries); vegas2.set_seed(1234); - VegasResult result2 = vegas2.integrate(gaussian_cpp); + VegasResult result2 = vegas2.integrate(gaussian_cpp, -1.0); std::cout << "Result 1: " << result1.value << " +/- " << result1.error << std::endl; std::cout << "Result 2: " << result2.value << " +/- " << result2.error << std::endl; @@ -56,12 +56,22 @@ int main() { // Test 2: Different seed should produce different result mchep::Vegas vegas3(n_iter, n_eval, n_bins, alpha, boundaries); vegas3.set_seed(5678); - VegasResult result3 = vegas3.integrate(gaussian_cpp); + VegasResult result3 = vegas3.integrate(gaussian_cpp, -1.0); std::cout << "Result 3: " << result3.value << " +/- " << result3.error << std::endl; assert(result1.value != result3.value); std::cout << "Different seed test passed.\n"; + std::cout << "\nTesting C++ API accuracy goal...\n"; + mchep::Vegas vegas4(20, 1e8, n_bins, alpha, boundaries); + vegas4.set_seed(1234); + VegasResult result4 = vegas4.integrate(gaussian_cpp, 1e-4); + + double accuracy = (result4.error / std::abs(result4.value)) * 100.0; + std::cout << "Accuracy goal test: value=" << result4.value << ", error=" << result4.error << ", acc=" << accuracy << std::endl; + assert(accuracy < 1e-4); + std::cout << "Accuracy goal test passed.\n"; + } catch (const std::runtime_error& e) { std::cerr << "Error: " << e.what() << std::endl; return 1; diff --git a/mchep_capi/tests/test_cppapi.output b/mchep_capi/tests/test_cppapi.output index e988aa3..1006d8f 100644 --- a/mchep_capi/tests/test_cppapi.output +++ b/mchep_capi/tests/test_cppapi.output @@ -7,3 +7,7 @@ Result 2: 2.23093 +/- 0.000105833 Same seed test passed. Result 3: 2.23088 +/- 0.000106151 Different seed test passed. + +Testing C++ API accuracy goal... +Accuracy goal test: value=2.23099, error=2.1762e-06, acc=9.75443e-05 +Accuracy goal test passed. diff --git a/mchep_capi/tests/test_cppapi_simd.cpp b/mchep_capi/tests/test_cppapi_simd.cpp new file mode 100644 index 0000000..cc56a16 --- /dev/null +++ b/mchep_capi/tests/test_cppapi_simd.cpp @@ -0,0 +1,67 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +// Integrand function for a 2D Gaussian: exp(-x^2 - y^2) +// This version is for SIMD and expects SoA data. +std::array gaussian_cpp_simd(const std::vector& x) { + // x is SoA: [x0, x1, x2, x3, y0, y1, y2, y3] + std::array results; + for (int i = 0; i < 4; ++i) { + double x_i = x[i]; + double y_i = x[i + 4]; + results[i] = std::exp(-(x_i * x_i + y_i * y_i)); + } + return results; +} + +int main() { + std::cout << "Testing 2D Gaussian integral (SIMD)..." << "\n"; + + const size_t n_iter = 10; + const size_t n_eval = 50000; + const size_t n_bins = 50; + const double alpha = 0.5; + + std::vector> boundaries = { + {-1.0, 1.0}, + {-1.0, 1.0} + }; + + try { + mchep::Vegas vegas(n_iter, n_eval, n_bins, alpha, boundaries); + vegas.set_seed(1234); // for reproducibility + + VegasResult result = vegas.integrate_simd(gaussian_cpp_simd, -1.0); + + const double expected = 2.230985; + const double multiplier = 2.5; + + double diff = std::fabs(result.value - expected); + assert(diff <= multiplier * result.error); + + std::cout << "Test passed!\n"; + std::cout << "Result: " << std::fixed << std::setprecision(6) << result.value << " +/- " << result.error << std::endl; + + std::cout << "\nTesting C++ API SIMD accuracy goal...\n"; + mchep::Vegas vegas2(20, 100000, n_bins, alpha, boundaries); + vegas2.set_seed(1234); + VegasResult result2 = vegas2.integrate_simd(gaussian_cpp_simd, 0.1); + + double accuracy = (result2.error / std::abs(result2.value)) * 100.0; + std::cout << "Accuracy goal test: value=" << result2.value << ", error=" << result2.error << ", acc=" << accuracy << std::endl; + assert(accuracy < 0.1); + std::cout << "Accuracy goal test passed.\n"; + + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; + } + + return 0; +} diff --git a/mchep_capi/tests/test_cppapi_simd.output b/mchep_capi/tests/test_cppapi_simd.output new file mode 100644 index 0000000..d374f59 --- /dev/null +++ b/mchep_capi/tests/test_cppapi_simd.output @@ -0,0 +1,7 @@ +Testing 2D Gaussian integral (SIMD)... +Test passed! +Result: 2.230914 +/- 0.000107 + +Testing C++ API SIMD accuracy goal... +Accuracy goal test: value=2.230734, error=0.001413, acc=0.063336 +Accuracy goal test passed. diff --git a/mchep_pyapi/Cargo.toml b/mchep_pyapi/Cargo.toml index 045c316..756037a 100644 --- a/mchep_pyapi/Cargo.toml +++ b/mchep_pyapi/Cargo.toml @@ -14,11 +14,9 @@ keywords.workspace = true [lib] name = "mchep" crate-type = ["cdylib"] +doc = false [dependencies] -pyo3 = { version = "0.21", features = ["extension-module"] } mchep = { path = "../mchep" } - -[features] -default = [] -mpi = ["mchep/mpi"] +pyo3 = { workspace = true, features = ["extension-module"] } +wide = { workspace = true } diff --git a/mchep_pyapi/src/vegas.rs b/mchep_pyapi/src/vegas.rs index 649d4ce..91de5ac 100644 --- a/mchep_pyapi/src/vegas.rs +++ b/mchep_pyapi/src/vegas.rs @@ -1,12 +1,15 @@ //! VEGAS interface. +use std::convert::TryFrom; + use pyo3::exceptions::PyValueError; use pyo3::prelude::*; use pyo3::types::PyList; -use mchep::integrand::Integrand; +use mchep::integrand::{Integrand, SimdIntegrand}; use mchep::vegas::{Vegas, VegasResult}; use mchep::vegasplus::VegasPlus; +use wide::f64x4; // A wrapper for Python callables to implement the Integrand trait #[pyclass(name = "Integrand")] @@ -46,6 +49,58 @@ impl PyIntegrand { } } +// A wrapper for Python callables to implement the SimdIntegrand trait +#[pyclass(name = "SimdIntegrand")] +struct PySimdIntegrand { + callable: PyObject, + dim: usize, +} + +impl SimdIntegrand for PySimdIntegrand { + fn dim(&self) -> usize { + self.dim + } + + fn eval_simd(&self, points: &[f64x4]) -> f64x4 { + Python::with_gil(|py| { + let py_points = PyList::empty_bound(py); + for d in 0..self.dim { + let point_dim = PyList::new_bound(py, &points[d].to_array()); + py_points.append(point_dim).unwrap(); + } + + let args = (py_points,); + self.callable + .call1(py, args) + .and_then(|result| result.extract::>(py)) + .and_then(|result_vec| { + <[f64; 4]>::try_from(result_vec) + .map(f64x4::from) + .map_err(|_| { + PyValueError::new_err("Integrand must return a list of 4 floats.") + .into() + }) + }) + .unwrap_or_else(|err| { + eprintln!("Error evaluating SIMD integrand: {err}"); + f64x4::splat(0.0) + }) + }) + } +} + +// CRITICAL: Mark as Send + Sync for parallel execution with Rayon +unsafe impl Send for PySimdIntegrand {} +unsafe impl Sync for PySimdIntegrand {} + +#[pymethods] +impl PySimdIntegrand { + #[new] + fn new(callable: PyObject, dim: usize) -> Self { + PySimdIntegrand { callable, dim } + } +} + #[pyclass(name = "VegasResult")] #[derive(Debug, Clone, Copy)] struct PyVegasResult { @@ -117,11 +172,23 @@ impl PyVegas { self.vegas.set_seed(seed); } - fn integrate_integrand(&mut self, py: Python, integrand: &PyIntegrand) -> PyVegasResult { - py.allow_threads(|| self.vegas.integrate(integrand).into()) + #[pyo3(signature = (integrand, target_accuracy = None))] + fn integrate_integrand( + &mut self, + py: Python, + integrand: &PyIntegrand, + target_accuracy: Option, + ) -> PyVegasResult { + py.allow_threads(|| self.vegas.integrate(integrand, target_accuracy).into()) } - fn integrate(&mut self, py: Python, callable: PyObject) -> PyResult { + #[pyo3(signature = (callable, target_accuracy = None))] + fn integrate( + &mut self, + py: Python, + callable: PyObject, + target_accuracy: Option, + ) -> PyResult { if !callable.bind(py).is_callable() { return Err(PyValueError::new_err("integrand must be callable")); } @@ -131,7 +198,40 @@ impl PyVegas { dim: self.dim, }; - Ok(py.allow_threads(|| self.vegas.integrate(&integrand).into())) + Ok(py.allow_threads(|| self.vegas.integrate(&integrand, target_accuracy).into())) + } + + #[pyo3(signature = (integrand, target_accuracy = None))] + fn integrate_simd_integrand( + &mut self, + py: Python, + integrand: &PySimdIntegrand, + target_accuracy: Option, + ) -> PyVegasResult { + py.allow_threads(|| self.vegas.integrate_simd(integrand, target_accuracy).into()) + } + + #[pyo3(signature = (callable, target_accuracy = None))] + fn integrate_simd( + &mut self, + py: Python, + callable: PyObject, + target_accuracy: Option, + ) -> PyResult { + if !callable.bind(py).is_callable() { + return Err(PyValueError::new_err("integrand must be callable")); + } + + let integrand = PySimdIntegrand { + callable, + dim: self.dim, + }; + + Ok(py.allow_threads(|| { + self.vegas + .integrate_simd(&integrand, target_accuracy) + .into() + })) } } @@ -173,11 +273,27 @@ impl PyVegasPlus { self.vegas_plus.set_seed(seed); } - fn integrate_integrand(&mut self, py: Python, integrand: &PyIntegrand) -> PyVegasResult { - py.allow_threads(|| self.vegas_plus.integrate(integrand).into()) + #[pyo3(signature = (integrand, target_accuracy = None))] + fn integrate_integrand( + &mut self, + py: Python, + integrand: &PyIntegrand, + target_accuracy: Option, + ) -> PyVegasResult { + py.allow_threads(|| { + self.vegas_plus + .integrate(integrand, target_accuracy) + .into() + }) } - fn integrate(&mut self, py: Python, callable: PyObject) -> PyResult { + #[pyo3(signature = (callable, target_accuracy = None))] + fn integrate( + &mut self, + py: Python, + callable: PyObject, + target_accuracy: Option, + ) -> PyResult { if !callable.bind(py).is_callable() { return Err(PyValueError::new_err("integrand must be callable")); } @@ -187,11 +303,58 @@ impl PyVegasPlus { dim: self.dim, }; - Ok(py.allow_threads(|| self.vegas_plus.integrate(&integrand).into())) + Ok(py.allow_threads(|| { + self.vegas_plus + .integrate(&integrand, target_accuracy) + .into() + })) + } + + #[pyo3(signature = (integrand, target_accuracy = None))] + fn integrate_simd_integrand( + &mut self, + py: Python, + integrand: &PySimdIntegrand, + target_accuracy: Option, + ) -> PyVegasResult { + py.allow_threads(|| { + self.vegas_plus + .integrate_simd(integrand, target_accuracy) + .into() + }) + } + + #[pyo3(signature = (callable, target_accuracy = None))] + fn integrate_simd( + &mut self, + py: Python, + callable: PyObject, + target_accuracy: Option, + ) -> PyResult { + if !callable.bind(py).is_callable() { + return Err(PyValueError::new_err("integrand must be callable")); + } + + let integrand = PySimdIntegrand { + callable, + dim: self.dim, + }; + + Ok(py.allow_threads(|| { + self.vegas_plus + .integrate_simd(&integrand, target_accuracy) + .into() + })) } #[cfg(feature = "mpi")] - fn integrate_mpi_integrand(&mut self, py: Python, integrand: &PyIntegrand) -> PyVegasResult { + #[pyo3(signature = (integrand, target_accuracy = None))] + fn integrate_mpi_integrand( + &mut self, + py: Python, + integrand: &PyIntegrand, + target_accuracy: Option, + ) -> PyVegasResult { use mpi::traits::*; // NOTE: MPI initialization should typically happen once at program start @@ -199,12 +362,20 @@ impl PyVegasPlus { py.allow_threads(|| { let universe = mpi::initialize().unwrap(); let world = universe.world(); - self.vegas_plus.integrate_mpi(integrand, &world).into() + self.vegas_plus + .integrate_mpi(integrand, &world, target_accuracy) + .into() }) } #[cfg(feature = "mpi")] - fn integrate_mpi(&mut self, py: Python, callable: PyObject) -> PyResult { + #[pyo3(signature = (callable, target_accuracy = None))] + fn integrate_mpi( + &mut self, + py: Python, + callable: PyObject, + target_accuracy: Option, + ) -> PyResult { use mpi::traits::*; if !callable.bind(py).is_callable() { @@ -219,7 +390,9 @@ impl PyVegasPlus { Ok(py.allow_threads(|| { let universe = mpi::initialize().unwrap(); let world = universe.world(); - self.vegas_plus.integrate_mpi(&integrand, &world).into() + self.vegas_plus + .integrate_mpi(&integrand, &world, target_accuracy) + .into() })) } @@ -256,6 +429,7 @@ pub fn register(parent_module: &Bound<'_, PyModule>) -> PyResult<()> { "import sys; sys.modules['mchep.vegas'] = m" ); m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/mchep_pyapi/tests/test_vegas.py b/mchep_pyapi/tests/test_vegas.py index ead4b44..b1992c6 100644 --- a/mchep_pyapi/tests/test_vegas.py +++ b/mchep_pyapi/tests/test_vegas.py @@ -3,7 +3,7 @@ import math import numpy as np -from mchep.vegas import Vegas, VegasPlus, Integrand +from mchep.vegas import Vegas, VegasPlus, Integrand, SimdIntegrand MULTIPLIER = 2.5 @@ -47,6 +47,59 @@ def gaussian(x): assert abs(result.value - expected) <= MULTIPLIER * result.error +def test_2d_gaussian_simd(): + """Test 2D Gaussian ∫∫exp(-x²-y²) dx over [-1,1]² using SIMD""" + expected = 2.230985 + + def gaussian_simd(x_soa): + results = [0.0] * 4 + for i in range(4): + x_i = x_soa[0][i] + y_i = x_soa[1][i] + results[i] = math.exp(-(x_i**2 + y_i**2)) + return results + + vegas = Vegas( + n_iter=10, + n_eval=50_000, + n_bins=50, + alpha=0.5, + boundaries=[(-1.0, 1.0), (-1.0, 1.0)], + ) + vegas.set_seed(1234) + + result = vegas.integrate_simd(gaussian_simd) + assert abs(result.value - expected) <= MULTIPLIER * result.error + + +def test_simd_integrand_class(): + """Test using SimdIntegrand wrapper for 2D Gaussian with SIMD""" + expected = 2.230985 + + class GaussianSimdClass: + def __call__(self, x_soa): + results = [0.0] * 4 + for i in range(4): + x_i = x_soa[0][i] + y_i = x_soa[1][i] + results[i] = math.exp(-(x_i**2 + y_i**2)) + return results + + integrand = SimdIntegrand(GaussianSimdClass(), dim=2) + + vegas = Vegas( + n_iter=10, + n_eval=50_000, + n_bins=50, + alpha=0.5, + boundaries=[(-1.0, 1.0), (-1.0, 1.0)], + ) + vegas.set_seed(1234) + + result = vegas.integrate_simd_integrand(integrand) + assert abs(result.value - expected) <= MULTIPLIER * result.error + + def test_lambda(): """Test using lambda function for ∫x² dx over [0,1]""" expected = 1 / 3 @@ -204,3 +257,105 @@ def f(x): np.testing.assert_almost_equal(result_plus1.value, result_plus2.value) np.testing.assert_almost_equal(result_plus1.error, result_plus2.error) + + +def test_accuracy_goal(): + """Test that the integration stops when the accuracy goal is reached.""" + expected = 2.230985 + + def gaussian(x): + return math.exp(-(x[0] ** 2 + x[1] ** 2)) + + vegas = Vegas( + n_iter=20, # More iterations to ensure accuracy is met + n_eval=100_000, + n_bins=50, + alpha=0.5, + boundaries=[(-1.0, 1.0), (-1.0, 1.0)], + ) + vegas.set_seed(1234) + + # Set a reasonable accuracy goal + target_accuracy = 0.1 # 0.1% + result = vegas.integrate(gaussian, target_accuracy=target_accuracy) + + # Check that the result is within the target accuracy + assert (result.error / abs(result.value)) * 100.0 < target_accuracy + assert abs(result.value - expected) <= MULTIPLIER * result.error + + vegas_plus = VegasPlus( + n_iter=20, + n_eval=100_000, + n_bins=50, + alpha=0.5, + n_strat=4, + beta=0.75, + boundaries=[(-1.0, 1.0), (-1.0, 1.0)], + ) + vegas_plus.set_seed(1234) + result_plus = vegas_plus.integrate(gaussian, target_accuracy=target_accuracy) + assert (result_plus.error / abs(result_plus.value)) * 100.0 < target_accuracy + assert abs(result_plus.value - expected) <= MULTIPLIER * result_plus.error + + +def test_non_rectangular_volume(): + """Test a 4D integral over a spherical volume. + + This demonstrates how to integrate over a non-rectangular volume by + defining the integrand function to be zero outside the desired region. + """ + expected = 1.0 + + def f_sph(x): + """A 4D Gaussian-like function defined within a sphere. + + The function is non-zero only inside a sphere of radius 0.2 centered + at (0.5, 0.5, 0.5, 0.5). + """ + dx2 = 0.0 + for d in range(4): + dx2 += (x[d] - 0.5) ** 2 + + if dx2 < 0.2**2: + return math.exp(-dx2 * 100.0) * 1115.3539360527281318 + else: + return 0.0 + + vegas = Vegas( + n_iter=10, + n_eval=200_000, + n_bins=50, + alpha=0.5, + boundaries=[(0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0)], + ) + vegas.set_seed(4321) + + result = vegas.integrate(f_sph) + assert abs(result.value - expected) <= MULTIPLIER * result.error + + +def test_vegasplus_simd(): + """Test VegasPlus integrator for 2D Gaussian with SIMD""" + expected = 2.230985 + + def gaussian_simd(x_soa): + results = [0.0] * 4 + for i in range(4): + x_i = x_soa[0][i] + y_i = x_soa[1][i] + results[i] = math.exp(-(x_i**2 + y_i**2)) + return results + + vegas_plus = VegasPlus( + n_iter=10, + n_eval=50_000, + n_bins=50, + alpha=0.5, + n_strat=4, + beta=0.75, + boundaries=[(-1.0, 1.0), (-1.0, 1.0)], + ) + vegas_plus.set_seed(1234) + + result = vegas_plus.integrate_simd(gaussian_simd) + assert abs(result.value - expected) <= MULTIPLIER * result.error