Add Vioreanu-Rokhlin quadrature finding

doc/quadrature.rst

@@ -62,4 +62,9 @@ Quadratures on the hypercube
 .. autoclass:: LegendreGaussTensorProductQuadrature
     :show-inheritance:
 
+Support for development of new quadrature rules
+-----------------------------------------------
+
+.. automodule:: modepy.quadrature.finding
+
 .. vim: sw=4

modepy/quadrature/finding.py

@@ -0,0 +1,261 @@
+"""
+.. autofunction:: adapt_2d_integrands_to_complex_arg
+.. autofunction:: orthogonalize_basis
+.. autofunction:: guess_nodes_vr
+.. autofunction:: find_weights_undetermined_coefficients
+.. autoclass:: QuadratureResidualJacobian
+.. autofunction:: quad_residual_and_jacobian
+"""
+
+from __future__ import annotations
+
+
+__copyright__ = """
+Copyright (C)  2024 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import operator
+from collections.abc import Callable, Sequence
+from dataclasses import dataclass
+from functools import reduce
+from typing import TypeAlias
+
+import numpy as np
+import numpy.linalg as la
+
+
+# FIXME: Better name?
+Integrand: TypeAlias = Callable[[np.ndarray], np.ndarray]
+
+
+@dataclass(frozen=True)
+class _ProductIntegrand:
+    functions: Sequence[Integrand]
+
+    def __call__(self, points: np.ndarray) -> np.ndarray:
+        return reduce(operator.mul, (f(points) for f in self.functions))
+
+
+@dataclass(frozen=True)
+class _ConjugateIntegrand:
+    function: Integrand
+
+    def __call__(self, points: np.ndarray) -> np.ndarray:
+        return self.function(points).conj()
+
+
+def _identity_integrand(points: np.ndarray) -> np.ndarray:
+    return points
+
+
+@dataclass(frozen=True)
+class _LinearCombinationIntegrand:
+    coefficients: np.ndarray
+    functions: Sequence[Integrand]
+
+    def __post_init__(self):
+        assert len(self.coefficients) == len(self.functions)
+
+    def __call__(self, points: np.ndarray) -> np.ndarray:
+        return sum(
+            (coeff * func(points)
+                for coeff, func in zip(self.coefficients, self.functions, strict=True)),
+            np.zeros(()))
+
+
+def linearly_combine(
+            coefficients: np.ndarray,
+            functions: Sequence[Integrand]
+        ) -> Integrand:
+    """
+    Takes advantage of associativity when linearly combining linear combinations.
+    """
+
+    lcfunctions = [
+        f for f in functions
+        if isinstance(f, _LinearCombinationIntegrand)
+    ]
+    if len(lcfunctions) != len(functions):
+        return _LinearCombinationIntegrand(coefficients, functions)
+
+    basis: list[Integrand] = []
+    n = len(lcfunctions)
+    matrix = np.zeros((n, n), dtype=np.complex128)
+    for i, f in enumerate(lcfunctions):
+        ncommon = min(len(basis), len(f.functions))
+        assert basis[:ncommon]  == f.functions[:ncommon]
+        basis.extend(f.functions[ncommon:])
+
+        ncoeff = len(f.coefficients)
+        matrix[i, :ncoeff] = f.coefficients
+
+    return _LinearCombinationIntegrand(coefficients @ matrix, basis)
+
+
+@dataclass(frozen=True)
+class _ComplexToNDAdapter:
+    function: Integrand
+
+    def __call__(self, points: np.ndarray) -> np.ndarray:
+        rpoints  = np.array([points.real, points.imag])
+        return self.function(rpoints)
+
+
+def adapt_2d_integrands_to_complex_arg(
+            functions: Sequence[Integrand]
+        ) -> Sequence[Integrand]:
+    return [_ComplexToNDAdapter(f) for f in functions]
+
+
+def _mass_matrix(
+            integrate: Callable[[Integrand], np.inexact],
+            basis: Sequence[Integrand],
+        ) -> np.ndarray:
+    n = len(basis)
+    mass_mat = np.zeros((n, n), dtype=np.complex128)
+    for i in range(n):
+        for j in range(i+1):
+            mass_mat[i, j] = integrate(
+                    _ProductIntegrand((basis[i], _ConjugateIntegrand(basis[j]))))
+    return mass_mat
+
+
+def orthogonalize_basis(
+            integrate: Callable[[Integrand], np.inexact],
+            basis: Sequence[Integrand],
+        ) -> Sequence[Integrand]:
+    r"""
+    Let :math:`\Omega\subset\mathbb C` be a convex domain.
+
+    :arg integrate: Computes an integral of the passed integrand over
+        :math:`\Omega`
+    """
+    n = len(basis)
+    mass_mat = _mass_matrix(integrate, basis)
+
+    l_factor = la.cholesky(mass_mat)
+
+    from scipy.linalg import solve_triangular
+    l_inv = solve_triangular(l_factor, np.eye(n), lower=True)
+
+    assert la.norm(np.triu(l_inv, 1), "fro") < 1e-14
+
+    return [
+        linearly_combine(l_inv[i, :i+1], basis[:i+1])
+        for i in range(n)
+    ]
+
+
+def guess_nodes_vr(
+            integrate: Callable[[Integrand], np.inexact],
+            onb: Sequence[Integrand],
+        ) -> np.ndarray:
+    """
+    Finds interpolation nodes based on the multiplication-operator technique in
+    [Vioreanu2011]_.
+
+    :arg integrate: must accurately integrate a product of two functions from *onb* and
+        a degree-1 monomial.
+    :arg onb: An orthonormal basis of functions
+    :returns: an array of shape ``(2, len(onb))`` containing nodes
+    """
+    n = len(onb)
+    mat = np.empty((n, n), dtype=np.complex128)
+    for i in range(n):
+        for j in range(n):
+            mat[i, j] = integrate(
+                _ProductIntegrand((
+                    _identity_integrand,
+                    onb[i],
+                    _ConjugateIntegrand(onb[j]))))
+
+    nodes_complex = la.eigvals(mat)
+    return np.array([nodes_complex.real, nodes_complex.imag])
+
+
+def find_weights_undetermined_coefficients(
+            integrands: Sequence[Integrand],
+            nodes: np.ndarray,
+            reference_integrals: np.ndarray,
+        ) -> np.ndarray:
+    """
+    :arg nodes: shaped ``(ndim, nnodes)``, real-valued
+    :arg reference_integrals: shaped ``(len(integrands),)``
+
+    .. note::
+
+        Tolerates overdetermined systems, will provide least squares solution.
+    """
+
+    if len(reference_integrals) != len(integrands):
+        raise ValueError(
+                "number of integrands must match number of reference integrals")
+    if len(reference_integrals) < len(nodes):
+        from warnings import warn
+        warn("Underdetermined quadrature system", stacklevel=2)
+
+    from modepy import vandermonde
+    return la.lstsq(vandermonde(integrands, nodes).T, reference_integrals)[0]
+
+
+@dataclass(frozen=True)
+class QuadratureResidualJacobian:
+    """
+    .. autoattribute:: residual
+    .. autoattribute:: dresid_dweights
+    """
+
+    residual: np.ndarray
+    dresid_dweights: np.ndarray
+
+
+def quad_residual_and_jacobian(
+            nodes: np.ndarray,
+            weights: np.ndarray,
+            integrands: Sequence[Integrand],
+            integrand_derivatives: Sequence[Integrand],
+            reference_integrals: np.ndarray,
+        ) -> QuadratureResidualJacobian:
+    """
+    :arg nodes: shaped ``(ndim, nnodes)``, real-valued
+    :arg weights: shaped ``(nnodes,)``, real-valued
+    :arg reference_integrals: shaped ``(len(integrands),)``
+    """
+    nintegrands = len(integrands)
+    _ndim, nnodes = nodes.shape
+
+    from modepy import vandermonde
+
+    vdm_t = vandermonde(integrands, nodes).T
+    residual = vdm_t @ weights - reference_integrals
+
+    dresid_dweights = np.empty((nintegrands, nnodes))
+    for inode in range(nnodes):
+        w_with_one = weights.copy()
+        w_with_one[inode] = 1
+        dresid_dweights = vdm_t @ w_with_one
+
+    return QuadratureResidualJacobian(
+        residual=residual,
+        dresid_dweights=dresid_dweights,
+    )