From df4cb1a0e9919f067b86becfe3ef02840aa08f00 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 28 Jan 2025 09:24:35 -0600 Subject: [PATCH] Add Vioreanu-Rokhlin quadrature finding Co-authored-by: Alex Fikl --- doc/conf.py | 7 + doc/index.rst | 1 + doc/quad_construction.rst | 4 + modepy/quadrature/construction.py | 451 ++++++++++++++++++++++++++ modepy/test/test_quad_construction.py | 284 ++++++++++++++++ 5 files changed, 747 insertions(+) create mode 100644 doc/quad_construction.rst create mode 100644 modepy/quadrature/construction.py create mode 100644 modepy/test/test_quad_construction.py diff --git a/doc/conf.py b/doc/conf.py index d5bf56ed..306fd786 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -23,3 +23,10 @@ } autodoc_member_order = "bysource" + +# FIXME: This would be nice to have, but I was not able to resolve the +# resulting mess of Sphinx errors. +# +# autodoc_type_aliases = { +# "Integrand": "modepy.quadrature.construction.Integrand", +# } diff --git a/doc/index.rst b/doc/index.rst index 573ea321..57636f9e 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -49,6 +49,7 @@ Contents modes nodes quadrature + quad_construction tools misc 🚀 Github diff --git a/doc/quad_construction.rst b/doc/quad_construction.rst new file mode 100644 index 00000000..2f218ff4 --- /dev/null +++ b/doc/quad_construction.rst @@ -0,0 +1,4 @@ +Support for construction of new quadrature rules +================================================ + +.. automodule:: modepy.quadrature.construction diff --git a/modepy/quadrature/construction.py b/modepy/quadrature/construction.py new file mode 100644 index 00000000..c9974ea7 --- /dev/null +++ b/modepy/quadrature/construction.py @@ -0,0 +1,451 @@ +""" +.. versionadded:: 2025.2 + +.. autodata:: Integrand +.. autoclass:: LinearCombinationIntegrand +.. autofunction:: linearly_combine +.. autofunction:: orthogonalize_basis +.. autofunction:: adapt_2d_integrands_to_complex_arg +.. autofunction:: guess_nodes_vioreanu_rokhlin +.. autofunction:: find_weights_undetermined_coefficients +.. autoclass:: QuadratureResidualJacobian +.. autofunction:: quad_residual_and_jacobian +.. autofunction:: quad_gauss_newton_increment +""" + +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 +from warnings import warn + +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: + """ + .. note:: + + New linear combinations should be created by :func:`linearly_combine`, + which flattens nested combinations by taking advantage of associativity, + substantially decreasing evaluation cost. + + .. autoattribute:: coefficients + .. autoattribute:: functions + .. automethod:: __call__ + """ + 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: + """Evaluate the linear combination of :attr:`functions` with + :attr:`coefficients` at *points*. + """ + 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] + ) -> LinearCombinationIntegrand: + r""" + Returns a linear combination of *functions* with coefficients. + If *functions* are themselves :class:`LinearCombinationIntegrand`\ s, + and all *functions* use a consistent basis, then the resulting + :class:`LinearCombinationIntegrand` will be flattened via associativity, + leading to faster execution time. + """ + + lcfunctions = [ + f for f in functions + if isinstance(f, LinearCombinationIntegrand) + ] + if len(lcfunctions) != len(functions): + if lcfunctions: + warn("Only some functions passed to linearly_combine were themselves " + "linear combination. This case is unhandled." + "Returning an un-flattened linear combination.", + stacklevel=2) + 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)) + if basis[:ncommon] != f.functions[:ncommon]: + warn("Not all functions passed to linearly_combine used the " + "same basis. Returning an un-flattened linear combination instead.", + stacklevel=2) + + return LinearCombinationIntegrand(coefficients, functions) + + basis.extend(f.functions[ncommon:]) + + ncoeff = len(f.coefficients) + matrix[i, :ncoeff] = f.coefficients + + return LinearCombinationIntegrand(coefficients @ matrix, basis) + + +def _mass_matrix( + integrate: Callable[[Integrand], np.inexact], + basis: Sequence[Integrand], + ) -> np.ndarray: + n = len(basis) + integrals = [ + [ + integrate(_ProductIntegrand((basis[i], _ConjugateIntegrand(basis[j])))) + for j in range(i+1) + ] + for i in range(n) + ] + + # Let np.array deal with figuring out the result dtype + return np.array([ + # Fill out the full matrix from the triangle computed above + [integrals[i][j] if i >= j else integrals[j][i] for j in range(n)] + for i in range(n) + ]) + + +def orthogonalize_basis( + integrate: Callable[[Integrand], np.inexact], + basis: Sequence[Integrand], + ) -> Sequence[LinearCombinationIntegrand]: + r""" + Let :math:`\Omega\subset\mathbb C^n` be a domain. (Domains + over the reals are allowable as well.) Returns linear combinations + of functions in *basis* that is orthogonal under the (complex-valued) + :math:`L^2` inner product induced by *integrate*. + + :arg integrate: Computes an integral of the passed integrand over + :math:`\Omega`. Must use integration nodes compatible with + *basis*. + :arg basis: A sequence of functions that accept an array of nodes + of shape either ``(n, nnodes)`` or ``(nnodes,)`` and return + an array of shape ``(nnodes,)`` with the value of the + basis function at the node. + """ + n = len(basis) + mass_mat = _mass_matrix(integrate, basis) + + l_factor = la.cholesky(mass_mat) + + try: + from scipy.linalg import solve_triangular + except ImportError: + # *shrug* I guess? The triangular version is also O(n**3). + l_inv = la.inv(l_factor) + else: + 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) + ] + + +@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]: + r"""Converts a list of :data:`Integrand`\ s taking nodes in real-valued + arrays of shape ``(2, nnodes)`` to ones accepting a single + complex-valued array of shape ``(nnodes,)``. See + :func:`guess_nodes_vioreanu_rokhlin` + for the main intended use case. + """ + return [_ComplexToNDAdapter(f) for f in functions] + + +def guess_nodes_vioreanu_rokhlin( + integrate: Callable[[Integrand], np.inexact], + onb: Sequence[Integrand], + ) -> np.ndarray: + r""" + Let :math:`\Omega\subset\mathbb C` be a convex domain. + Finds interpolation nodes for :math:`\Omega` based on the + multiplication-operator technique in + [Vioreanu2011]_. May be useful as an initial guess for a Gauss-Newton + process to find new quadrature rules, as driven by + :func:`quad_residual_and_jacobian`. + + :arg integrate: Must accurately integrate a product of two functions from *onb* and + a degree-1 monomial over :math:`\Omega`. + :arg onb: An orthonormal basis of functions. Each function takes + in an array of (complex) nodes and returns an array + of the same shape containing the function values. + Functions are expected to be orthogonal with respect to + the complex-valued inner product. + + Integrands accepting real-valued node coordinates in two dimensions + (shaped ``(2, nnodes)`` may be converted to functions acceptable by + this function via :func:`adapt_2d_integrands_to_complex_arg`.) + :returns: An array of shape ``(len(onb),)`` containing nodes, complex-valued. + + .. note:: + + This function is based on complex arithmetic and thus only + usable for domains $\Omega$ with one and two (real-valued) + dimensions. + + .. note:: + + (Empirically) it seems acceptable for the functions + in *onb* to be exclusively real-valued, i.e., particularly, + no assumptions on complex differentiability are needed. + + .. note:: + + As noted in [Vioreanu2011]_, the returned nodes should obey + the symmetry of the domain. + + """ + n = len(onb) + mat = np.array([ + [ + integrate( + _ProductIntegrand(( + _identity_integrand, + onb[i], + _ConjugateIntegrand(onb[j])))) + for j in range(n) + ] + for i in range(n) + ]) + + return la.eigvals(mat) + + +def find_weights_undetermined_coefficients( + integrands: Sequence[Integrand], + nodes: np.ndarray, + reference_integrals: np.ndarray, + ) -> np.ndarray: + """ + Uses the method of undetermined coefficients to find weights + for a quadrature rule using *nodes*, for the provided + *integrands*. + + :arg integrands: A sequence of functions that accept an array of nodes + of shape either ``(ndim, nnodes)`` or ``(nnodes,)`` and return + an array of shape ``(nnodes,)`` with the value of the + basis function at the node. + :arg nodes: An array with shape ``(ndim, nnodes)`` or ``(nnodes,)``, real-valued. + Must be compatible with *integrands*. + :arg reference_integrals: An array with shape ``(len(integrands),)`` + + .. note:: + + This function also supports overdetermined systems. In that case, it + will provide the least squares solution. + + + .. note:: + + Unlike :func:`guess_nodes_vioreanu_rokhlin`, domains of any dimensionality + are allowed. + """ + + 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: + """ + Contains residual value and Jacobian for a quadrature residual, see + :func:`quad_residual_and_jacobian`. May be reused for residuals + that include additional terms, such as those that penalize asymmetry + of the noddes. + + .. autoattribute:: residual + .. autoattribute:: dresid_dweights + .. autoattribute:: dresid_dnodes + """ + + residual: np.ndarray + """Shaped ``(nintegrands,)``.""" + + dresid_dweights: np.ndarray + """Shaped ``(nintegrands, nnodes)``.""" + + dresid_dnodes: np.ndarray + """Shaped ``(nintegrands, ndim*nnodes)``.""" + + +def quad_residual_and_jacobian( + nodes: np.ndarray, + weights: np.ndarray, + integrands: Sequence[Integrand], + integrand_derivatives: Sequence[Sequence[Integrand]], + reference_integrals: np.ndarray, + ) -> QuadratureResidualJacobian: + r""" + Computes the residual and Jacobian of the objective function + + .. math:: + + \begin{bmatrix} + \sum_{j=0}^{\text{nnodes-1}} \psi_0(\boldsymbol x_j) w_j - I_0\\ + \vdots\\ + \sum_{j=0}^{\text{nnodes-1}} \psi_{\text{nintegrands-1}} + (\boldsymbol x_j) w_j - I_{\text{nintegrands}-1} + \end{bmatrix}. + + Typically used with :meth:`quad_gauss_newton_increment` + to drive an iterative process for finding nodes and weights of a quadrature + rule. + + An initial guess for the nodes may be obtained via + :func:`guess_nodes_vioreanu_rokhlin`, + and an initial guess for the weights (given nodes) may be found + via :func:`find_weights_undetermined_coefficients`. + + .. note:: + + Unlike :func:`guess_nodes_vioreanu_rokhlin`, domains of any dimensionality + are allowed. + + :arg nodes: An array with shape ``(ndim, nnodes)`` or ``(nnodes,)``, real-valued. + Must be compatible with *integrands*. + :arg weights: An array with shape ``(nnodes,)``, real-valued + :arg integrands: A sequence of functions that accept an array of nodes + of shape either ``(ndim, nnodes)`` or ``(nnodes,)`` and return + an array of shape ``(nnodes,)`` with the value of the + basis function at the node. + :arg integrand_derivatives: Derivatives of *integrands* along the + coordinate axes, one list of + functions per axis, with each sub-list matching *integrands* + in structure. + :arg reference_integrals: Integrals of *integrands*, shaped ``(len(integrands),)``. + """ + nintegrands = len(integrands) + ndim, _nnodes = nodes.shape + + if __debug__: + if len(integrand_derivatives) != ndim: + raise ValueError("sequence of integrand derivatives must have ndim entries") + + for iaxis, ax_derivatives in enumerate(integrand_derivatives): + if len(ax_derivatives) != nintegrands: + raise ValueError( + "number of integrands must match number of integrand" + f"along axis {iaxis} (0-based)") + + from modepy import vandermonde + + vdm_t = vandermonde(integrands, nodes).T + residual = vdm_t @ weights - reference_integrals + + dresid_dnodes = np.hstack([ + vandermonde(ax_derivatives, nodes).T * weights + for ax_derivatives in integrand_derivatives + ]) + + return QuadratureResidualJacobian( + residual=residual, + dresid_dweights=vdm_t, + dresid_dnodes=dresid_dnodes, + ) + + +def quad_gauss_newton_increment( + qrj: QuadratureResidualJacobian + ) -> tuple[np.ndarray, np.ndarray]: + """Return the Gauss-Newton increment based on the residual and Jacobian + (see :func:`quad_residual_and_jacobian`), + separated into the weight increment and the nodes increment, + in the two entries of the returned tuple. The nodes increment + has shape ``(ndim, nnodes)``. + + .. note:: + + [Vioreanu2011]_ suggests that step size control should be used + in the Gauss-Newton process. No notion of step size control + is included in the returned increments. + """ + + full_jacobian = np.hstack([qrj.dresid_dweights, qrj.dresid_dnodes]) + full_increment = -la.lstsq(full_jacobian, qrj.residual)[0] + + _nintegrands, nnodes = qrj.dresid_dweights.shape + + return full_increment[:nnodes], full_increment[nnodes:].reshape(-1, nnodes) diff --git a/modepy/test/test_quad_construction.py b/modepy/test/test_quad_construction.py new file mode 100644 index 00000000..19644f2e --- /dev/null +++ b/modepy/test/test_quad_construction.py @@ -0,0 +1,284 @@ +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. +""" + +from collections.abc import Callable + +import numpy as np +import numpy.linalg as la + +import modepy as mp +import modepy.quadrature.construction as constr + + +def test_quad_finder_lincomb() -> None: + basis = constr.adapt_2d_integrands_to_complex_arg( + mp.monomial_basis_for_space( + mp.PN(2, order=5), + mp.Simplex(2)).functions, + ) + + quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + base_quad = mp.Quadrature( + nodes=quad.nodes[0] + quad.nodes[1]*1j, + weights=quad.weights, + ) + + n = len(basis) + rng = np.random.default_rng(15) + + m1 = rng.uniform(size=(n, n)) + b1 = [constr.linearly_combine(m1[i], basis) for i in range(n)] + b1_ref = [constr.LinearCombinationIntegrand(m1[i], basis) for i in range(n)] + + for test_f, ref_f in zip(b1, b1_ref, strict=True): + assert la.norm(test_f(base_quad.nodes) - ref_f(base_quad.nodes), 2) < 1e-15 + + m2 = rng.uniform(size=(n, n)) + b2 = [constr.linearly_combine(m2[i], b1) for i in range(n)] + b2_ref = [constr.LinearCombinationIntegrand(m2[i], b1_ref) for i in range(n)] + + for test_f in b2: + assert isinstance(test_f, constr.LinearCombinationIntegrand) + assert all(not isinstance(subf, constr.LinearCombinationIntegrand) + for subf in test_f.functions) + + for test_f, ref_f in zip(b2, b2_ref, strict=True): + assert la.norm(test_f(base_quad.nodes) - ref_f(base_quad.nodes), 2) < 1e-13 + + +def test_orthogonalize_basis() -> None: + orig_basis = mp.monomial_basis_for_space(mp.PN(2, order=7), mp.Simplex(2)) + basis = constr.adapt_2d_integrands_to_complex_arg(orig_basis.functions) + + quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + base_quad = mp.Quadrature( + nodes=quad.nodes[0] + quad.nodes[1]*1j, + weights=quad.weights, + ) + + def integrate(integrand: Callable[[np.ndarray], np.ndarray]) -> np.inexact: + res = base_quad(integrand) + assert isinstance(res, np.inexact) + return res + + onb = constr.orthogonalize_basis(integrate, basis) + onb = constr.orthogonalize_basis(integrate, onb) + + eye = constr._mass_matrix(integrate, onb) + assert la.norm(eye - np.eye(len(onb))) < 2e-12 + + +def test_guess_nodes_vr() -> None: + order = 7 + shape = mp.Simplex(2) + orig_basis = mp.orthonormal_basis_for_space( + mp.PN(2, order=order), + shape) + basis = constr.adapt_2d_integrands_to_complex_arg(orig_basis.functions) + + quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + base_quad = mp.Quadrature( + nodes=quad.nodes[0] + quad.nodes[1]*1j, + weights=quad.weights, + ) + + def integrate(integrand: Callable[[np.ndarray], np.ndarray]) -> np.inexact: + res = base_quad(integrand) + assert isinstance(res, np.inexact) + return res + + onb = constr.orthogonalize_basis(integrate, basis) + + eye = constr._mass_matrix(integrate, onb) + assert la.norm(eye - np.eye(len(onb))) < 1e-14 + + nodes_complex = constr.guess_nodes_vioreanu_rokhlin(integrate, onb) + nodes = np.array([nodes_complex.real, nodes_complex.imag]) + + from modepy.tools import estimate_lebesgue_constant + assert estimate_lebesgue_constant(order, nodes, shape) < 32 + + if 0: + import matplotlib.pyplot as plt + plt.plot(nodes[0], nodes[1], "x") + plt.gca().set_aspect("equal") + plt.show() + + +def test_undetermined_coeffs() -> None: + order = 7 + space = mp.PN(2, order=order) + shape = mp.Simplex(2) + basis = mp.orthonormal_basis_for_space(space, shape) + nodes = mp.edge_clustered_nodes_for_space(space, shape) + + ref_quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + + ref_integrals = np.array([ref_quad(f) for f in basis.functions]) + weights = constr.find_weights_undetermined_coefficients( + basis.functions, nodes, ref_integrals) + quad = mp.Quadrature(nodes, weights) + + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + ) + + from modepy.tools import Monomial + for comb in gnitstam(order, 2): + f = Monomial(comb) + i_f = quad(f) + ref = f.simplex_integral() + err = abs(i_f - ref) + print(order, repr(f), err) + assert err < 2e-15, (err, comb, i_f, ref) + + +def test_quad_residual_derivatives() -> None: + order = 7 + shape = mp.Simplex(2) + onb = mp.orthonormal_basis_for_space( + mp.PN(2, order=order), + shape) + complex_onb = constr.adapt_2d_integrands_to_complex_arg(onb.functions) + + base_quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + base_quad_complex = mp.Quadrature( + nodes=base_quad.nodes[0] + base_quad.nodes[1]*1j, + weights=base_quad.weights, + ) + + def integrate(integrand: Callable[[np.ndarray], np.ndarray]) -> np.inexact: + res = base_quad_complex(integrand) + assert isinstance(res, np.inexact) + return res + + eye = constr._mass_matrix(integrate, complex_onb) + assert la.norm(eye - np.eye(len(complex_onb))) < 7e-14 + + nodes_complex = constr.guess_nodes_vioreanu_rokhlin(integrate, complex_onb) + nodes = np.array([nodes_complex.real, nodes_complex.imag]) + + ref_integrals = np.array([base_quad(f) for f in onb.functions]) + weights = constr.find_weights_undetermined_coefficients( + onb.functions, nodes, ref_integrals) + + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + ) + + from modepy.tools import Monomial + + quad = mp.Quadrature(nodes, weights) + + for comb in gnitstam(order, 2): + f = Monomial(comb) + assert abs(quad(f) - f.simplex_integral()) < 2e-15 + + all_derivatives = [onb.derivatives(iaxis) for iaxis in range(shape.dim)] + + qrj = constr.quad_residual_and_jacobian( + nodes, weights, onb.functions, all_derivatives, ref_integrals) + + rng = np.random.default_rng(17) + + # check weight derivatives + w_direction = rng.uniform(size=weights.shape) + for h in [0.1]: + qrj_offset = constr.quad_residual_and_jacobian( + nodes, weights + h * w_direction, + onb.functions, all_derivatives, + ref_integrals) + + err = la.norm( + (qrj_offset.residual - qrj.residual) / h + - qrj.dresid_dweights @ w_direction, 2) + # Residual is linear in the weights + assert err < 1e-12 + + # check node derivatives + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + n_direction = rng.uniform(size=nodes.shape) + for h in [10**-i for i in range(1, 4)]: + qrj_offset = constr.quad_residual_and_jacobian( + nodes + h * n_direction, weights, + onb.functions, all_derivatives, + ref_integrals) + + err = la.norm( + (qrj_offset.residual - qrj.residual) / h + - qrj.dresid_dnodes @ n_direction.reshape(-1), 2) + eoc_rec.add_data_point(h, float(err)) + + print(eoc_rec) + assert eoc_rec.order_estimate() >= 0.95 + + +def test_quad_gauss_newton() -> None: + order = 7 + shape = mp.Simplex(2) + onb = mp.orthonormal_basis_for_space( + mp.PN(2, order=order), + shape) + complex_onb = constr.adapt_2d_integrands_to_complex_arg(onb.functions) + + base_quad = mp.XiaoGimbutasSimplexQuadrature(20, 2) + base_quad_complex = mp.Quadrature( + nodes=base_quad.nodes[0] + base_quad.nodes[1]*1j, + weights=base_quad.weights, + ) + + def integrate(integrand: Callable[[np.ndarray], np.ndarray]) -> np.inexact: + res = base_quad_complex(integrand) + assert isinstance(res, np.inexact) + return res + + eye = constr._mass_matrix(integrate, complex_onb) + assert la.norm(eye - np.eye(len(complex_onb))) < 7e-14 + + nodes_complex = constr.guess_nodes_vioreanu_rokhlin(integrate, complex_onb) + nodes = np.array([nodes_complex.real, nodes_complex.imag]) + + integrands = mp.orthonormal_basis_for_space( + mp.PN(2, order=11), + shape) + + ref_integrals = np.array([base_quad(f) for f in integrands.functions]) + weights = constr.find_weights_undetermined_coefficients( + integrands.functions, nodes, ref_integrals) + + all_derivatives = [integrands.derivatives(iaxis) for iaxis in range(shape.dim)] + + for _ in range(8): + qrj = constr.quad_residual_and_jacobian( + nodes, weights, integrands.functions, all_derivatives, ref_integrals) + weights_inc, nodes_inc = constr.quad_gauss_newton_increment(qrj) + weights += weights_inc + nodes += nodes_inc + + resid_norm = la.norm(qrj.residual) + + assert resid_norm < 1e-14