feat: add opt-in VEGAS error reporting (return_error)

[gomezzz]

What

VEGAS already computes a standard deviation, chi-squared and goodness-of-fit
internally, then threw them away and returned only the bare integral. This adds
an opt-in return_error=True that surfaces them:

from torchquad import VEGAS
result = VEGAS().integrate(fn, dim=2, N=40000, integration_domain=dom, return_error=True)
print(result)        # VEGASResult(4.00144 +/- 0.0014, chi2/dof = 1, Q = 0.41)
result.integral, result.sdev, result.chi2, result.dof, result.Q

• New VEGASResult dataclass (integral, sdev, chi2, dof, Q,
  nr_of_fevals) with a compact value +/- sdev (chi2/dof, Q) repr; Q is the
  chi-squared goodness-of-fit p-value via SciPy's regularized upper incomplete
  gamma. Exposed as torchquad.VEGASResult.
• The tensor fields keep the integration backend. Only integral is on the
  gradient path; sdev/chi2 derive from VEGAS's detached variances and are
  not differentiable.
• return_error defaults to False → the existing bare-value return is
  unchanged (no breaking change).

Test plan

• vegas_result_test.py (numpy + torch): result fields well-formed and
  self-consistent, error estimate brackets the true error, backward-compatible
  default return; numpy bit-for-bit reproducibility (torch bit-reproducibility
  isn't asserted — GPU reductions aren't deterministic)
• existing vegas_test.py still green
• ruff / pydoclint / vulture clean; sphinx-build -W builds
  (VEGASResult auto-documented via __all__, per-field attribute docstrings)

Roadmap F2.

[github-actions]

Coverage Report

File	Stmts	Miss	Cover	Missing
torchquad
__init__.py	27	2	93%	67–68
torchquad/integration
base_integrator.py	49	7	86%	28, 71–74, 78, 113, 117, 122
boole.py	30	3	90%	54, 73–76
gaussian.py	51	3	94%	84, 108, 128
grid_integrator.py	97	3	97%	50, 164, 268
integration_grid.py	49	4	92%	64–65, 132, 135
monte_carlo.py	93	3	97%	79, 101, 230
qmc.py	20	1	95%	67
rng.py	60	7	88%	71–74, 93–94, 102
simpson.py	29	3	90%	52, 71–74
utils.py	131	11	92%	35–36, 102, 173, 177, 184, 196, 198, 204, 299, 306
vegas.py	167	4	98%	106–107, 212, 359
vegas_map.py	103	8	92%	247–257
vegas_result.py	25	1	96%	22
torchquad/utils
deployment_test.py	192	63	67%	53–55, 65–66, 69–70, 78–79, 91–93, 164, 175–176, 180–181, 183–184, 196–199, 211–213, 221, 232–234, 237, 247, 250, 253, 256–259, 267–275, 280–288, 297–306, 316–323
enable_cuda.py	15	6	60%	17–22
set_log_level.py	14	2	86%	45–46
set_precision.py	34	9	74%	35–38, 50–51, 70–72
TOTAL	1282	140	89%

Tests	Skipped	Failures	Errors	Time
125	0 💤	0 ❌	0 🔥	2m 34s ⏱️

[gomezzz]

Claude: Thanks — all three addressed in 467d25f:

1. Differentiability (Important): narrowed the claim in the VEGASResult docstring and the PR description — only integral is on the gradient path; sdev/chi2 come from VEGAS's detached per-iteration variances and are not differentiable.
2. _get_chisq zero-variance: added the and sig2 != 0.0 guard (matching _get_error/_get_result) so the now-user-facing chi2/Q can't become inf/nan on a zero-variance iteration.
3. Q oracle: the fields test now asserts result.Q == gammaincc(dof/2, chi2/2) from the returned chi2/dof, so a swapped-argument or wrong-formula regression fails instead of passing a range check.

Re-ran vegas_test.py + vegas_result_test.py: 8 passed.

[gomezzz]

Claude: ### Accuracy comparison: Monte Carlo vs Sobol QMC vs VEGAS

Relative error against the analytic value at an equal sample budget of N = 2¹⁴ = 16384, numpy backend, float64, fixed seed=0. Every integrand is separable so the exact integral is known in closed form.

Integrand	dim	Monte Carlo	Sobol QMC	VEGAS
smooth: ∏ cos(πx/2)	2	3.84e-03	4.53e-09	4.33e-04
smooth: ∏ e^(-x)	5	4.84e-04	8.55e-06	8.15e-05
sharply peaked: ∏ e^(-((x-0.5)/0.1)²)	3	2.48e-02	4.17e-03	2.46e-03
smooth: ∏ e^(-x)	8	1.34e-03	1.11e-04	3.46e-04

Takeaways

• On smooth integrands Sobol QMC is the clear winner — up to ~6 orders of magnitude better than plain MC (and consistently ahead of VEGAS), matching the expected O(1/N) vs O(1/√N) convergence.
• On the sharply peaked integrand, VEGAS's adaptive stratification pays off and it edges out Sobol; Sobol still beats plain MC ~6×, but its low-discrepancy advantage shrinks when the mass is concentrated in a small region.
• Plain Monte Carlo is the weakest in every case.

Caveat: VEGAS's N is an approximate max function-evaluation budget spread across warmup + refinement iterations, so its effective sample count isn't identical to MC/Sobol — this is a same-budget comparison, not an identical-sample-count one.

[gomezzz]

Claude: Added an error-calibration test against the analytic collection (c3b18c2) — thanks, the earlier tests only exercised one integrand.

vegas_result_test.py::_vegas_error_calibration_test now runs every real, scalar-integrand function from integration_test_functions (1-D and 3-D; VEGAS doesn't support complex or multi-integrand functions, so those are skipped) with return_error=True, and for each asserts:

• VEGAS converged to the closed-form value (worst observed rel error ~0.8%);
• its self-reported sdev brackets the true error: |integral − exact| < 8·sdev.

Measured true_error/sdev across the collection: max ~3.4, mean ~1 on both numpy and torch — i.e. the estimate is a genuine 1σ-scale error bar, not just a plausible positive number. The 8× bound leaves headroom for VEGAS's stochastic (and GPU-nondeterministic) spread.

Runs on numpy and torch. Related tracked guidance is in the REVIEW.md update on #244.

[gomezzz]

Claude: Merged develop in to resolve conflicts after #246 landed — VEGAS integrate() now carries both return_error (this PR) and args (from #246), verified together (vegas + result + args tests, 22 passed). CI green. Includes the error-calibration test against the analytic collection.