Codebase cleanup: CI, docs, lints, dedup, tests

.github/workflows/ci.yml

@@ -14,6 +14,9 @@ jobs:
   test:
     name: Test
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        features: ["", "--all-features", "--no-default-features"]
     steps:
       - uses: actions/checkout@v6
 
@@ -24,10 +27,10 @@ jobs:
         uses: Swatinem/rust-cache@v2
 
       - name: Build
-        run: cargo build
+        run: cargo build ${{ matrix.features }}
 
       - name: Run tests
-        run: cargo test
+        run: cargo test ${{ matrix.features }}
 
   lint:
     name: Lint
@@ -47,7 +50,7 @@ jobs:
         run: cargo fmt --all -- --check
 
       - name: Clippy
-        run: cargo clippy -- -D warnings
+        run: cargo clippy --all-features -- -D warnings
 
   docs:
     name: Documentation
@@ -83,53 +86,3 @@ jobs:
 
       - name: Test (default features)
         run: cargo test
-
-  coverage:
-    name: Coverage
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v6
-
-      - name: Install Rust
-        uses: dtolnay/rust-toolchain@stable
-
-      - name: Cache cargo
-        uses: Swatinem/rust-cache@v2
-
-      - name: Install cargo-tarpaulin
-        run: cargo install cargo-tarpaulin
-
-      - name: Generate coverage
-        run: cargo tarpaulin --out xml --timeout 300
-        continue-on-error: true
-
-      - name: Upload coverage
-        uses: codecov/codecov-action@v5
-        with:
-          files: cobertura.xml
-          token: ${{ secrets.CODECOV_TOKEN }}
-          fail_ci_if_error: false
-
-  security:
-    name: Security Audit
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v6
-
-      - name: Install Rust
-        uses: dtolnay/rust-toolchain@stable
-
-      - name: Cache cargo
-        uses: Swatinem/rust-cache@v2
-
-      - name: Install cargo-audit
-        run: cargo install cargo-audit --locked
-
-      - name: Audit dependencies
-        run: cargo audit
-
-      - name: Install cargo-deny
-        run: cargo install cargo-deny --locked
-
-      - name: Check dependency policies
-        run: cargo deny check

CHANGELOG.md

@@ -5,6 +5,40 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## [0.1.0] - Unreleased
 
-[Unreleased]: https://github.com/Entrolution/bilby/compare/v0.1.0...HEAD
+Initial release of bilby.
+
+### Added
+
+- **Gauss-Legendre** quadrature with Bogaert (2014) asymptotic expansion for large n and Newton iteration for small n
+- **Gauss-Kronrod** embedded pairs (G7K15, G10K21, G15K31, G20K41, G25K51, G30K61) for error estimation
+- **Classical families** via Golub-Welsch eigenvalue solver:
+  - Gauss-Jacobi (arbitrary alpha, beta)
+  - Gauss-Hermite (physicists' convention)
+  - Gauss-Laguerre (generalised, arbitrary alpha)
+  - Gauss-Chebyshev Types I and II (closed-form)
+  - Gauss-Lobatto (both endpoints)
+  - Gauss-Radau (left or right endpoint)
+  - Clenshaw-Curtis (nested, Chebyshev extrema)
+- **Adaptive integration** -- QUADPACK-style global subdivision with configurable GK pair, break points, absolute/relative tolerance
+- **Domain transforms** -- `integrate_infinite`, `integrate_semi_infinite_lower`, `integrate_semi_infinite_upper` for unbounded domains
+- **Tanh-sinh** (double-exponential) quadrature for endpoint singularities
+- **Cauchy principal value** integration via the subtraction technique
+- **Oscillatory integration** via Filon-Clenshaw-Curtis with adaptive fallback for small omega
+- **Weighted integration** -- unified API for product integration with Jacobi, Laguerre, Hermite, Chebyshev, and log weight functions
+- **Multi-dimensional cubature**:
+  - Tensor product rules (isotropic and anisotropic)
+  - Smolyak sparse grids with Clenshaw-Curtis basis
+  - Genz-Malik adaptive cubature (degree-7/degree-5 embedded rule)
+  - Monte Carlo (plain pseudo-random with Welford variance estimation)
+  - Quasi-Monte Carlo via Sobol and Halton low-discrepancy sequences
+- **Generic over `F: Float`** via num-traits for core rule types
+- **`no_std` compatible** -- works with `alloc` only (disable default `std` feature)
+- **`parallel` feature** -- rayon-based parallel variants:
+  - `GaussLegendre::new_par` for parallel node generation
+  - `integrate_composite_par` for parallel composite rules
+  - `integrate_par` / `integrate_box_par` for cubature rules
+  - `MonteCarloIntegrator::integrate_par` for parallel MC/QMC
+- **Precomputed rule cache** (`cache` module, requires `std`) with `GL5`, `GL10`, `GL15`, `GL20`, `GL50`, `GL100` lazy singletons
+- **Benchmarks** via criterion for node generation, 1D integration, and cubature

CONTRIBUTING.md

@@ -56,7 +56,7 @@ Before submitting a PR, ensure:
 cargo fmt
 
 # Run clippy
-cargo clippy -- -D warnings
+cargo clippy --all-features -- -D warnings
 
 # Check documentation
 cargo doc --no-deps

src/adaptive.rs

@@ -97,24 +97,28 @@ impl Default for AdaptiveIntegrator {
 
 impl AdaptiveIntegrator {
     /// Set the Gauss-Kronrod pair used for panel evaluation.
+    #[must_use]
     pub fn with_pair(mut self, pair: GKPair) -> Self {
         self.pair = pair;
         self
     }
 
     /// Set the absolute error tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set the relative error tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Set the maximum number of function evaluations.
+    #[must_use]
     pub fn with_max_evals(mut self, n: usize) -> Self {
         self.max_evals = n;
         self
@@ -125,6 +129,11 @@ impl AdaptiveIntegrator {
     /// Returns a [`QuadratureResult`] with the best estimate and error bound.
     /// If the tolerance could not be achieved within the evaluation budget,
     /// `converged` will be `false` but the best estimate is still returned.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if both tolerances are non-positive
+    /// or if `max_evals` is less than one Gauss-Kronrod panel evaluation.
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -142,6 +151,11 @@ impl AdaptiveIntegrator {
     /// Break points split the interval at known discontinuities or singularities.
     /// The adaptive algorithm then refines each sub-interval independently in a
     /// global priority queue.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if any break point is outside
+    /// `(a, b)`, is NaN, or if both tolerances are non-positive.
     pub fn integrate_with_breaks<G>(
         &self,
         a: f64,
@@ -189,6 +203,9 @@ impl AdaptiveIntegrator {
 
         // Seed the priority queue with initial intervals
         for &(ia, ib) in intervals {
+            // Exact comparison is intentional: degenerate zero-width
+            // intervals are an edge case, not a floating-point tolerance issue.
+            #[allow(clippy::float_cmp)]
             if ia == ib {
                 continue;
             }
@@ -210,10 +227,7 @@ impl AdaptiveIntegrator {
         while !self.tolerance_met(total_estimate, total_error)
             && num_evals + 2 * evals_per_call <= self.max_evals
         {
-            let worst = match heap.pop() {
-                Some(s) => s,
-                None => break,
-            };
+            let Some(worst) = heap.pop() else { break };
 
             // Bisect the worst interval
             let mid = 0.5 * (worst.a + worst.b);
@@ -340,6 +354,16 @@ where
 /// Break points indicate locations of discontinuities or singularities.
 /// The interval is split at these points before adaptive refinement begins.
 ///
+/// # Example
+///
+/// ```
+/// use bilby::adaptive_integrate_with_breaks;
+///
+/// // Integrate |x| over [-1, 1] with a break at the cusp
+/// let result = adaptive_integrate_with_breaks(|x: f64| x.abs(), -1.0, 1.0, &[0.0], 1e-12).unwrap();
+/// assert!((result.value - 1.0).abs() < 1e-12);
+/// ```
+///
 /// # Errors
 ///
 /// Returns [`QuadratureError::InvalidInput`] if any break point is outside `(a, b)` or NaN.
@@ -507,4 +531,62 @@ mod tests {
         assert!(r.num_evals > 0);
         assert!(r.num_evals <= 200, "num_evals={}", r.num_evals);
     }
+
+    /// Reversed bounds for exp(x): ∫₁⁰ exp(x) dx = -(e - 1).
+    #[test]
+    fn reversed_bounds_exp() {
+        let forward = adaptive_integrate(f64::exp, 0.0, 1.0, 1e-12).unwrap();
+        let reverse = adaptive_integrate(f64::exp, 1.0, 0.0, 1e-12).unwrap();
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
+
+    /// Reversed bounds for a polynomial: ∫₁⁰ x² dx = -1/3.
+    #[test]
+    fn reversed_bounds_polynomial() {
+        let r = adaptive_integrate(|x| x * x, 1.0, 0.0, 1e-12).unwrap();
+        assert!(r.converged, "err={}", r.error_estimate);
+        assert!(
+            (r.value + 1.0 / 3.0).abs() < 1e-12,
+            "value={}, expected -1/3",
+            r.value
+        );
+    }
+
+    /// Reversed bounds with break points: ∫₁⁻¹ |x| dx = -1.
+    #[test]
+    fn reversed_bounds_with_breaks() {
+        let forward =
+            adaptive_integrate_with_breaks(|x: f64| x.abs(), -1.0, 1.0, &[0.0], 1e-12).unwrap();
+        let reverse =
+            adaptive_integrate_with_breaks(|x: f64| x.abs(), 1.0, -1.0, &[0.0], 1e-12).unwrap();
+        assert!(forward.converged, "forward err={}", forward.error_estimate);
+        assert!(reverse.converged, "reverse err={}", reverse.error_estimate);
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
+
+    /// Reversed bounds with break points for step function.
+    #[test]
+    fn reversed_bounds_step_with_breaks() {
+        let step = |x: f64| if x < 2.0 { 1.0 } else { 3.0 };
+        let forward = adaptive_integrate_with_breaks(step, 0.0, 4.0, &[2.0], 1e-12).unwrap();
+        let reverse = adaptive_integrate_with_breaks(step, 4.0, 0.0, &[2.0], 1e-12).unwrap();
+        assert!(forward.converged, "forward err={}", forward.error_estimate);
+        assert!(reverse.converged, "reverse err={}", reverse.error_estimate);
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
 }

src/cauchy_pv.rs

@@ -48,18 +48,21 @@ impl Default for CauchyPV {
 
 impl CauchyPV {
     /// Set absolute tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set relative tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Set maximum number of function evaluations.
+    #[must_use]
     pub fn with_max_evals(mut self, n: usize) -> Self {
         self.max_evals = n;
         self
@@ -68,6 +71,13 @@ impl CauchyPV {
     /// Compute PV ∫ₐᵇ f(x)/(x-c) dx.
     ///
     /// The singularity c must be strictly inside (a, b).
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN.
+    /// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside
+    /// `(a, b)` or if `f(c)` is not finite.
+    #[allow(clippy::many_single_char_names)] // a, b, c, f, g are conventional in quadrature
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -128,6 +138,23 @@ impl CauchyPV {
 }
 
 /// Convenience: Cauchy principal value integration with default settings.
+///
+/// # Example
+///
+/// ```
+/// use bilby::pv_integrate;
+///
+/// // PV integral of x^2/(x - 0.3) over [0, 1]
+/// let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln();
+/// let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap();
+/// assert!((result.value - exact).abs() < 1e-7);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN.
+/// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside
+/// `(a, b)` or if `f(c)` is not finite.
 pub fn pv_integrate<G>(
     f: G,
     a: f64,