|
| 1 | +import numpy |
| 2 | +from autoray import numpy as anp |
| 3 | +from .grid_integrator import GridIntegrator |
| 4 | + |
| 5 | + |
| 6 | +class Gaussian(GridIntegrator): |
| 7 | + """Gaussian quadrature methods inherit from this. Default behaviour is Gauss-Legendre quadrature on [-1,1].""" |
| 8 | + |
| 9 | + def __init__(self): |
| 10 | + super().__init__() |
| 11 | + self.name = "Gauss-Legendre" |
| 12 | + self.root_fn = numpy.polynomial.legendre.leggauss |
| 13 | + self.root_args = () |
| 14 | + self.default_integration_domain = [[-1, 1]] |
| 15 | + self.transform_interval = True |
| 16 | + self._cache = {} |
| 17 | + |
| 18 | + def integrate(self, fn, dim, N=8, integration_domain=None, backend=None): |
| 19 | + """Integrates the passed function on the passed domain using Simpson's rule. |
| 20 | +
|
| 21 | + Args: |
| 22 | + fn (func): The function to integrate over. |
| 23 | + dim (int): Dimensionality of the integration domain. |
| 24 | + N (int, optional): Total number of sample points to use for the integration. Should be odd. Defaults to 3 points per dimension if None is given. |
| 25 | + integration_domain (list or backend tensor, optional): Integration domain, e.g. [[-1,1],[0,1]]. Defaults to [-1,1]^dim. It also determines the numerical backend if possible. |
| 26 | + backend (string, optional): Numerical backend. This argument is ignored if the backend can be inferred from integration_domain. Defaults to the backend from the latest call to set_up_backend or "torch" for backwards compatibility. |
| 27 | +
|
| 28 | + Returns: |
| 29 | + backend-specific number: Integral value |
| 30 | + """ |
| 31 | + return super().integrate(fn, dim, N, integration_domain, backend) |
| 32 | + |
| 33 | + def _weights(self, N, dim, backend, requires_grad=False): |
| 34 | + """return the weights, broadcast across the dimensions, generated from the polynomial of choice |
| 35 | +
|
| 36 | + Args: |
| 37 | + N (int): number of nodes |
| 38 | + dim (int): number of dimensions |
| 39 | + backend (string): which backend array to return |
| 40 | +
|
| 41 | + Returns: |
| 42 | + backend tensor: the weights |
| 43 | + """ |
| 44 | + weights = anp.array(self._cached_points_and_weights(N)[1], like=backend) |
| 45 | + if backend == "torch": |
| 46 | + weights.requires_grad = requires_grad |
| 47 | + return anp.prod( |
| 48 | + anp.array( |
| 49 | + anp.stack( |
| 50 | + list(anp.meshgrid(*([weights] * dim))), like=backend, dim=0 |
| 51 | + ) |
| 52 | + ), |
| 53 | + axis=0, |
| 54 | + ).ravel() |
| 55 | + else: |
| 56 | + return anp.prod( |
| 57 | + anp.meshgrid(*([weights] * dim), like=backend), axis=0 |
| 58 | + ).ravel() |
| 59 | + |
| 60 | + def _roots(self, N, backend, requires_grad=False): |
| 61 | + """return the roots generated from the polynomial of choice |
| 62 | +
|
| 63 | + Args: |
| 64 | + N (int): number of nodes |
| 65 | + backend (string): which backend array to return |
| 66 | +
|
| 67 | + Returns: |
| 68 | + backend tensor: the roots |
| 69 | + """ |
| 70 | + roots = anp.array(self._cached_points_and_weights(N)[0], like=backend) |
| 71 | + if requires_grad: |
| 72 | + roots.requires_grad = True |
| 73 | + return roots |
| 74 | + |
| 75 | + @property |
| 76 | + def _grid_func(self): |
| 77 | + """ |
| 78 | + function for generating a grid to be integrated over i.e., the polynomial roots, resized to the domain. |
| 79 | + """ |
| 80 | + |
| 81 | + def f(a, b, N, requires_grad, backend=None): |
| 82 | + return self._resize_roots(a, b, self._roots(N, backend, requires_grad)) |
| 83 | + |
| 84 | + return f |
| 85 | + |
| 86 | + def _resize_roots(self, a, b, roots): # scale from [-1,1] to [a,b] |
| 87 | + """resize the roots based on domain of [a,b] |
| 88 | +
|
| 89 | + Args: |
| 90 | + a (backend tensor): lower bound |
| 91 | + b (backend tensor): upper bound |
| 92 | + roots (backend tensor): polynomial nodes |
| 93 | +
|
| 94 | + Returns: |
| 95 | + backend tensor: rescaled roots |
| 96 | + """ |
| 97 | + return roots |
| 98 | + |
| 99 | + # credit for the idea https://github.com/scipy/scipy/blob/dde50595862a4f9cede24b5d1c86935c30f1f88a/scipy/integrate/_quadrature.py#L72 |
| 100 | + def _cached_points_and_weights(self, N): |
| 101 | + """wrap the calls to get weights/roots in a cache |
| 102 | +
|
| 103 | + Args: |
| 104 | + N (int): number of nodes to return |
| 105 | + backend (string): which backend to use |
| 106 | +
|
| 107 | + Returns: |
| 108 | + tuple: nodes and weights |
| 109 | + """ |
| 110 | + root_args = (N, *self.root_args) |
| 111 | + if not isinstance(N, int): |
| 112 | + if hasattr(N, "item"): |
| 113 | + root_args = (N.item(), *self.root_args) |
| 114 | + else: |
| 115 | + raise NotImplementedError( |
| 116 | + f"N {N} is not an int and lacks an `item` method" |
| 117 | + ) |
| 118 | + if root_args in self._cache: |
| 119 | + return self._cache[root_args] |
| 120 | + self._cache[root_args] = self.root_fn(*root_args) |
| 121 | + return self._cache[root_args] |
| 122 | + |
| 123 | + @staticmethod |
| 124 | + def _apply_composite_rule(cur_dim_areas, dim, hs, domain): |
| 125 | + """Apply "composite" rule for gaussian integrals |
| 126 | +
|
| 127 | + cur_dim_areas will contain the areas per dimension |
| 128 | + """ |
| 129 | + # We collapse dimension by dimension |
| 130 | + for cur_dim in range(dim): |
| 131 | + cur_dim_areas = ( |
| 132 | + 0.5 |
| 133 | + * (domain[cur_dim][1] - domain[cur_dim][0]) |
| 134 | + * anp.sum(cur_dim_areas, axis=len(cur_dim_areas.shape) - 1) |
| 135 | + ) |
| 136 | + return cur_dim_areas |
| 137 | + |
| 138 | + |
| 139 | +class GaussLegendre(Gaussian): |
| 140 | + """Gauss Legendre quadrature rule in torch. See https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature. |
| 141 | +
|
| 142 | + Examples |
| 143 | + -------- |
| 144 | + >>> gl=torchquad.GaussLegendre() |
| 145 | + >>> integral = gl.integrate(lambda x:np.sin(x), dim=1, N=101, integration_domain=[[0,5]]) #integral from 0 to 5 of np.sin(x) |
| 146 | + |TQ-INFO| Computed integral was 0.7163378000259399 #analytic result = 1-np.cos(5)""" |
| 147 | + |
| 148 | + def __init__(self): |
| 149 | + super().__init__() |
| 150 | + |
| 151 | + def _resize_roots(self, a, b, roots): # scale from [-1,1] to [a,b] |
| 152 | + return ((b - a) / 2) * roots + ((a + b) / 2) |
0 commit comments