Gaussian quadrature

README.md

@@ -3,7 +3,7 @@
 *** Based on https://github.com/othneildrew/Best-README-Template
 -->
 
-![Read the Docs (version)](https://img.shields.io/readthedocs/torchquad/main?style=flat-square) ![GitHub Workflow Status (branch)](https://img.shields.io/github/actions/workflow/status/esa/torchquad/.github/workflows/run_tests.yml/main?style=flat-square) ![GitHub last commit](https://img.shields.io/github/last-commit/esa/torchquad?style=flat-square)
+![Read the Docs (version)](https://img.shields.io/readthedocs/torchquad/main?style=flat-square) [![Tests](https://github.com/esa/torchquad/actions/workflows/run_tests.yml/badge.svg)](https://github.com/esa/torchquad/actions/workflows/run_tests.yml) ![GitHub last commit](https://img.shields.io/github/last-commit/esa/torchquad?style=flat-square)
 ![GitHub](https://img.shields.io/github/license/esa/torchquad?style=flat-square) ![Conda (channel only)](https://img.shields.io/conda/vn/conda-forge/torchquad?style=flat-square) ![PyPI](https://img.shields.io/pypi/v/torchquad?style=flat-square) ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/torchquad?style=flat-square)
 
 ![GitHub contributors](https://img.shields.io/github/contributors/esa/torchquad?style=flat-square)

docs/source/tutorial.rst

@@ -756,4 +756,40 @@ Now let's see how to do this a bit more simply, and in a way that provides signf
 
     torch.all(torch.isclose(result_vectorized, result)) # True!
 
+Custom Integrators
+------------------
+
+It is of course possible to extend our provided Integrators, perhaps for a special class of functions or for a new algorithm.
+
+.. code:: ipython3
+
+    class GaussHermite(Gaussian):
+    """Gauss Hermite quadrature rule in torch, for integrals of the form :math:`\\int_{-\\infty}^{+\\infty} e^{-x^{2}} f(x) dx`. It will correctly integrate
+    polynomials of degree :math:`2n - 1` or less over the interval
+    :math:`[-\\infty, \\infty]` with weight function :math:`f(x) = e^{-x^2}`. See https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature
+    """
+
+    def __init__(self):
+        super().__init__()
+        self.name = "Gauss-Hermite"
+        self.root_fn = scipy.special.roots_hermite
+        self.default_integration_domain = [[-1 * numpy.inf, numpy.inf]]
+        self.wrapper_func = None
+
+    @staticmethod
+    def _apply_composite_rule(cur_dim_areas, dim, hs, domain):
+        """Apply "composite" rule for gaussian integrals
+        cur_dim_areas will contain the areas per dimension
+        """
+        # We collapse dimension by dimension
+        for _ in range(dim):
+            cur_dim_areas = anp.sum(cur_dim_areas, axis=len(cur_dim_areas.shape) - 1)
+        return cur_dim_areas
+
+    gh=torchquad.GaussHermite()
+    integral=gh.integrate(lambda x: 1-x,dim=1,N=200) #integral from -inf to inf of np.exp(-(x**2))*(1-x)
+    # Computed integral was 1.7724538509055168.
+    # analytic result = sqrt(pi)
+
+
 

torchquad/__init__.py

@@ -9,6 +9,9 @@
 from .integration.simpson import Simpson
 from .integration.boole import Boole
 from .integration.vegas import VEGAS
+from .integration.gaussian import GaussLegendre
+from .integration.grid_integrator import GridIntegrator
+from .integration.base_integrator import BaseIntegrator
 
 from .integration.rng import RNG
 
@@ -22,16 +25,15 @@
 from .utils.deployment_test import _deployment_test
 
 __all__ = [
+    "GridIntegrator",
+    "BaseIntegrator",
     "IntegrationGrid",
     "MonteCarlo",
     "Trapezoid",
     "Simpson",
     "Boole",
     "VEGAS",
     "GaussLegendre",
-    "GaussJacobi",
-    "GaussLaguerre",
-    "GaussHermite",
     "RNG",
     "plot_convergence",
     "plot_runtime",

torchquad/integration/base_integrator.py

@@ -37,12 +37,14 @@ def _eval(self, points, weights=None, args=None):
             weights (backend tensor, optional): Integration weights. Defaults to None.
             args (list or tuple, optional): Any arguments required by the function. Defaults to None.
         """
-        result, num_points = self.evaluate_integrand(self._fn, points,weights=weights,args=args)
+        result, num_points = self.evaluate_integrand(
+            self._fn, points, weights=weights, args=args
+        )
         self._nr_of_fevals += num_points
         return result
 
     @staticmethod
-    def evaluate_integrand(fn, points, weights=None,args=None):
+    def evaluate_integrand(fn, points, weights=None, args=None):
         """Evaluate the integrand function at the passed points
 
         Args:
@@ -56,12 +58,12 @@ def evaluate_integrand(fn, points, weights=None,args=None):
             int: Number of evaluated points
         """
         num_points = points.shape[0]
-        
+
         if args is None:
             args = ()
-        
+
         result = fn(points, *args)
-        
+
         if infer_backend(result) != infer_backend(points):
             warnings.warn(
                 "The passed function's return value has a different numerical backend than the passed points. Will try to convert. Note that this may be slow as it results in memory transfers between CPU and GPU, if torchquad uses the GPU."
@@ -77,8 +79,15 @@ def evaluate_integrand(fn, points, weights=None,args=None):
             )
 
         if weights is not None:
+            if len(result.shape) > 1:
+                integrand_shape = anp.array(
+                    result.shape[1:], like=infer_backend(points)
+                )
+                weights = anp.repeat(
+                    anp.expand_dims(weights, axis=1), anp.prod(integrand_shape)
+                ).reshape((weights.shape[0], *(integrand_shape)))
             result *= weights
-            
+
         return result, num_points
 
     @staticmethod
@@ -105,3 +114,4 @@ def _check_inputs(dim=None, N=None, integration_domain=None):
             if dim is not None and dim != dim_domain:
                 raise ValueError(
                     "The dimension of the integration domain must match the passed function dimensionality dim."
+                )

torchquad/integration/boole.py

@@ -28,7 +28,7 @@ def integrate(self, fn, dim, N=None, integration_domain=None, backend=None):
         return super().integrate(fn, dim, N, integration_domain, backend)
 
     @staticmethod
-    def _apply_composite_rule(cur_dim_areas, dim, hs):
+    def _apply_composite_rule(cur_dim_areas, dim, hs, domain):
         """Apply composite Boole quadrature.
         cur_dim_areas will contain the areas per dimension
         """