diff --git a/doc/index.rst b/doc/index.rst index 7dfe6cd25b14dbf97ee83edef8e4a41ca091115d..91e2e26b1bf724b779fe436a2ad23c2c1a7bed7f 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -45,6 +45,7 @@ Contents modes nodes quadrature + symbolic tools misc diff --git a/doc/symbolic.rst b/doc/symbolic.rst new file mode 100644 index 0000000000000000000000000000000000000000..2d5f58c390070dc8a33f8f82306e7720423961f1 --- /dev/null +++ b/doc/symbolic.rst @@ -0,0 +1,46 @@ +Symbolic +======== + +Mathematical formulation +------------------------ + +.. automodule:: modepy.symbolic + +Symbolic primitives +------------------- + +.. currentmodule:: modepy.symbolic + +.. autoclass:: TemplateQuadratureFamily + +.. autodata:: LEGENDRE_GAUSS_FAMILY + +.. autodata:: CHEBYSHEV_GAUSS_1_FAMILY + +.. autodata:: CHEBYSHEV_GAUSS_2_FAMILY + +.. autodata:: CLENSHAW_CURTIS_FAMILY + +.. autodata:: FEJER_1_FAMILY + +.. autodata:: FEJER_2_FAMILY + +.. autoclass:: QuadratureRingElement + + .. automethod:: __add__ + .. automethod:: __mul__ + .. automethod:: __rmul__ + .. automethod:: __pow__ + .. automethod:: get_exactness_diagram + .. automethod:: get_quadrature_rule + +.. autoclass:: QuadratureRingGenerator + +Tensor product rules +-------------------- + +.. currentmodule:: modepy.symbolic.tensor_product + +.. autoclass:: TensorProductQuadratureFactory + +.. vim: sw=4 diff --git a/modepy/__init__.py b/modepy/__init__.py index 5931c6f9bcf052b8334a1b27d978c755668e19c6..b29364c61a24044835dc685b6e0b8f60b846273f 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -27,8 +27,11 @@ THE SOFTWARE. from modepy.modes import ( jacobi, grad_jacobi, simplex_onb, grad_simplex_onb, simplex_monomial_basis, grad_simplex_monomial_basis, + tensor_product_basis, simplex_onb_with_mode_ids, + simplex_monomial_basis_with_mode_ids, simplex_best_available_basis, grad_simplex_best_available_basis) -from modepy.nodes import equidistant_nodes, warp_and_blend_nodes +from modepy.nodes import ( + equidistant_nodes, warp_and_blend_nodes, tensor_product_nodes) from modepy.matrices import (vandermonde, resampling_matrix, differentiation_matrices, inverse_mass_matrix, mass_matrix, @@ -54,9 +57,10 @@ __all__ = [ "simplex_onb", "grad_simplex_onb", "simplex_monomial_basis", "grad_simplex_monomial_basis", "simplex_best_available_basis", "grad_simplex_best_available_basis", - "tensor_product_basis", + "tensor_product_basis", "simplex_onb_with_mode_ids", + "simplex_monomial_basis_with_mode_ids", - "equidistant_nodes", "warp_and_blend_nodes", + "equidistant_nodes", "warp_and_blend_nodes", "tensor_product_nodes", "vandermonde", "resampling_matrix", "differentiation_matrices", "inverse_mass_matrix", "mass_matrix", "modal_face_mass_matrix", diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..34670f7d9404766514f624d8c6fb4ed10b942a88 --- /dev/null +++ b/modepy/symbolic/__init__.py @@ -0,0 +1,994 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__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 numpy as np +from modepy.quadrature import Quadrature + +import logging +logger = logging.getLogger(name=__name__) + +__doc__ = r""" +The symbolic manipulation of quadrature rules enables construction of tensor +product quadrature rules and more importantly sparse grid quadrature rules +in general dimensions. + +Terminologies +^^^^^^^^^^^^^ + +Dimension + Denote :math:`d \geq 0` as the dimension. + +Domains + We restrict our discussion to :math:`d`-dimensional template + domains :math:`\Omega_d = [-1, 1]^d`. + + If :math:`d = 0`, the domain is a single point, denoted as + :math:`\Omega_0 = \{0\}`. + +Function spaces + Let :math:`\omega: \Omega_d \rightarrow [0,+\infty)` be a measurable + function. :math:`L^2(\Omega_d, \omega)` denotes the + :math:`\omega`-weighted :math:`L^2` space. + + We consider a subspace that is Riemann integrable, + + .. math:: + + V_{d,\omega} := \{f \in L^2(\Omega_d, \omega) | + \sup_{x\in\Omega_d} f(x)\omega(x) < \infty, + \text{ and } \lim_{x\rightarrow x_0} f(x)\omega(x) = f(x_0)\omega(x_0) + \text{ a.e.}\}. + + If :math:`d = 0`, :math:`V_{0, \omega}` degenerates to the set of scalar + constants :math:`\mathbb{R}`. + +Integral + The Riemann integral is a linear functional defined as + :math:`I_\omega(f) := \int_{\Omega_d}f(x)\omega(x)dx, + \forall f\in V_{d, \omega}`. + + In the case of :math:`d = 0`, + :math:`I_\omega(f) := f(0)\omega(0)`, so that + :math:`I_1(1) = 2^d` holds in all dimensions. + +Evaluation + The evaluation functional is defined as :math:`ev_x(f) := f(x)`. + +Quadrature + A quadrature rule is a finite dimensional linear combination of + evaluation functionals + + .. math:: + + quad(f, \omega) := (\sum_{i=1}^{n} w_i ev_{x_i})(f\omega), + \forall f\in V_{d, \omega}. + + It is identified by the set of nodes and corresponding weights + :math:`\{(x_i, w_i), 1\leq i \leq n, w_i \neq 0\}`. + + A special rule is when the set of nodes and weights is empty. + It returns :math:`0` for any integrand and is thus "non-dimensional", + but we assign its dimension as :math:`-1` for convenience. + We denote it as :math:`0` as well, as will be shown later, it + is the additive identity. + + In addition, a :math:`0`-dimensional rule behaves like a + constant (the weight of its single node). We thus denote it as + its weight value. And as will be shown later, :math:`1` is the + multiplicative identity. + + In practice, the empty rule can be treated as any dimensional rule + + .. math:: + 0(f) = 0, \forall f \in V_{d,\omega}, \forall d, \forall \omega. + + Note that nodes are not allowed to have zero weights, but are + allowed to have negative weights. + +Quadrature chain + A sequence of quadrature rules :math:`\{quad_i\}_{i=1}^{N}` form a + quadrature chain if all quadrature rules in the sequence share the + same dimension and weight function. + +Consistency + An unbounded quadrature chain :math:`s_{d, \omega} = \{quad_i\}_{i=1}^{\infty}` + is said to be consistent, if + :math:`\lim_{i\rightarrow\infty} quad_i(f, \omega) = I_\omega(f), + \forall f\in V_{d, \omega}`. + +Embedding + Given :math:`1 \leq d_1 \leq d_2`, and a multi-index + :math:`\alpha = (\alpha_k)_{k=1}^{d_1}`, + with :math:`1\leq \alpha_1 < \alpha_2 < \ldots < \alpha_{d_1} \leq d_2`, + define the natural embedding of + :math:`X = [0, 1]^{d_1}` into :math:`Y = [0, 1]^{d_2}` as + + .. math:: + + em^\alpha: (x_i)_{i=1}^{d_1} \rightarrow (y_j)_{j=1}^{d_2}, + + where :math:`y_j = x_i \text{ if } j = \alpha_i`, and + :math:`y_j = 0` otherwise. + + Note that, :math:`Jacobian(em^\alpha) = 1` since :math:`\alpha_k` are + in increasing order. + + The embedding of evaluation operator is thus defined as + :math:`em^\alpha(ev_x) := ev(y_\alpha)`. And the embedding of quadrature + rules is defined as + :math:`em^\alpha(quad)(f) := + (\sum_{i=1}^{n} em^\alpha(w_i) em^\alpha(ev_{x_i}))(f\omega)`, + where :math:`em^\alpha(w_i) = w_i` if :math:`d_1 = d_2`, + otherwise :math:`em^\alpha(w_i) = 0` (thus removing the node). + Therefore, :math:`em^\alpha(quad) = 0` if :math:`d_1 < d_2`. + + If :math:`d_1=0`, there is only one natural embedding that maps the point + to the :math:`d_2`-dimensional origin. For simplicity, unless indicated + specifically, :math:`em` is used as a shorthand of + :math:`em^{(1, 2, \ldots, d_2)}`. + + In all cases, :math:`em^\alpha` acts like identity when :math:`d_1 = d_2`. + +Tensor product + Given :math:`d_1, d_2`, define the tensor product of + :math:`\Omega_{d_1}` and :math:`\Omega_{d_2}` as + + .. math:: + \Omega_{d_1} \otimes \Omega_{d_2} = \Omega_{d_1 + d_2}, + + associated with the bilinear map + + .. math:: + + \phi: \Omega_{d_1} \times \Omega_{d_2} \rightarrow + \Omega_{d_1} \otimes \Omega_{d_2}, + ((x_1, \ldots, x_{d_1}), (y_1, \ldots, y_{d_2})) + \mapsto (x_1, \ldots, x_{d_1}, y_1, \ldots, y_{d_2}). + + Specially, :math:`\Omega_d \otimes \Omega_0 = \Omega_0 \otimes + \Omega_d = \Omega_d`, + with :math:`\phi` being identity map. + + Naturally the tensor product of evaluation operators can be defined as + + .. math:: + + ev_x \otimes ev_y := ev_{\phi(x, y)}, + + as well as the tensor product of functions + + .. math:: + + (f \otimes g)(\phi(x, y)) := f(x) g(y), + +Algebra of quadrature rules +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Addition + Given two quadrature rules :math:`a, b`, both + in :math:`d` dimensions, the addition is defined as + + .. math:: + + (a+b)(f) := a(f) + b(f) = + (\sum_{i=1}^{n_a} w_{a,i} ev_{x_{a,i}})(f\omega_a) + + (\sum_{i=1}^{n_b} w_{b,j} ev_{x_{b,j}})(f\omega_b), + \forall f\in V_{d, \omega} \\ + = + (\sum_{i=1}^{n_a} + \frac{w_{a,i} \omega_a(x_{a, i})}{\omega_{a+b}(x_{a, i})} + ev_{x_{a,i}} + + \sum_{i=1}^{n_b} + \frac{w_{b,j} \omega_b(x_{b, j})}{\omega_{a+b}(x_{b, j})} + ev_{x_{b,j}}) + (f\omega_{a+b}), + \forall f\in V_{d, \omega} \\ + + where, :math:`\omega_{a+b}` is the weight function of the sum. + If :math:`\omega_a = \omega_b = \omega`, + then :math:`\omega_{a+b} = \omega`; otherwise the sum has + weight function constant one, :math:`\omega_{a+b} = 1`. + + More generally, if :math:`a, b` have different dimensions, define + :math:`d = max(dim(a), dim(b))`, and define + + .. math:: + + a + b := em(a) + em(b) = argmax_{quad\in\{a, b\}}dim(quad), + + Properties of addition: + + - :math:`dim(a + b) \in \{-1, max( dim(a), dim(b) )\}` + - (identity): :math:`a + 0 = a` + - (commutative): :math:`a + b = b + a` + - (associative): :math:`(a + b) + c = a + (b + c)` + +Multiplication + The product of two quadrature rules :math:`a, b` is defined by + their tensor product, + + .. math:: + + (ab)(f) := (a\otimes b)(f) = (\sum_{i=1}^{n_a}\sum_{j=1}^{n_b} + w_{a,i} w_{b,j} ev_{x_{a,i}} \otimes ev_{x_{b,j}})(f\omega), + \forall f\in V_{d_a + d_b, \omega}, + + where :math:`\omega = \omega_a \otimes \omega_b`. + + Properties of multiplication: + + - :math:`dim(a + b) \in \{-1, dim(a) + dim(b)\}` + - (identity): :math:`1 a = a 1 = a` + - (associative): :math:`(a b) c = a (b c)` + + Obviously, the + multiplication of two quadrature rules is not commutative in + general, unless one of the quadrature rules is of dimension + :math:`0` or less. + +Multiplication distributes over addition + + It can be shown by plugging in the definitions and checking that, + + .. math:: + + (a + b) c = a c + b c, + + and that + + .. math:: + + a (b + c) = a b + a c. + + +Ring construction +^^^^^^^^^^^^^^^^^ + +We want to construct quadrature rings from a set of generators, using +the operations defined above. The generators are classic families of +quadrature rules for 1D interval, such as Gauss-Legendre quadrature and +Gauss-Chebyshev quadrature. + +Denote the set of 0-th dimensional rules as :math:`Q^0`, +then :math:`\{0\} \cup Q^0` is a field isomorphic to :math:`\mathbb{R}`. + +For example, denote the n-th degree second kind Gauss-Chebyshev quadrature +rule as :math:`quad^{GC2}_n`, :math:`n\geq 0`. +Then :math:`Q^{GC2} := \{quad^{GC2}_n\}_{n=0}^{\infty}` is a +consistent quadrature chain. +And the quadrature ring generated by :math:`Q^{GC2}` is obtained as follows + +.. math:: + + R^{GC2} = span(\{0\} \cup Q^0 \cup Q^{GC2}) + := \left\{ \sum_{i=1}^{N} \prod_{j=1}^{p_i} q_{i,j} \right\}, + \text{ } q_{i,j} \in \{0\} \cup Q^0 \cup Q^{GC2}. + +We refer to each term in the summation as a **word**, and each +:math:`q_{i,j}` as a **letter**. Naturally, the set that letters +are drawn from (:math:`\{0\} \cup Q^0 \cup Q^{GC2}`) is called +the **alphabet**. + +:math:`R^{GC2}` is isomorphic to the non-commutative +polynomial ring on infinitely many indeterminates +(variables), where an isomorphism can be simply defined by mapping +:math:`Q^{GC2}` to the list of variables, while mapping :math:`Q^0` to nonzero +constants, and the empty rule :math:`0` to the constant :math:`0`. + +In other words, +:math:`R^{GC2}` is the monoid algebra over :math:`\{0\} \cup Q^0` +for the free (associative, unital) monoid on :math:`Q^{GC2}`. +""" + +# {{{ class TemplateQuadratureFamily + + +class TemplateQuadratureFamily(): + r"""Family of 1-D quadrature rules over [-1, 1], indexed by an integral + parameter *level*. + + .. attribute :: family_id + + An id number used to distinguish different families of rules when + added / multiplied together. + + .. attribute :: name + + (Optional) name of the family of quadrature rules. + + .. attribute :: is_nested + + A bool value indicating whether the quadrature rules can be nested. + + .. attribute :: weight_func + + A lambda / function object that gives the weight function omega. + + .. attribute :: effective_measure + + The sum of quadrature weights. For Gauss quadrature w.r.t weight + function :math:`\omega(x)`, this equals to + :math:`\int_{[-1,1]} \omega(x) dx`. + + .. attribute :: find_node_at_level + + A lambda / function object that when given an node and a level, + calculates the level-local id of the node inside the given level. + It should return -1 if the node does not belong to the specified level. + The input node's information is given as the node's first appearing + level and its index within that level. + + .. attribute :: level_to_nqpts + + A lambda / function object that maps level to number of quadrature + points. + + .. attribute :: level_to_exact_to + + A lambda / function object that maps level to the polynomial degree up + to which the quadrature is exact. + + .. attribute :: quad_class + + Quadrature class that is used to generate nodes and weights. + + .. attribute :: quad_class_extra_kwargs + + Extra keyword arguments for quad_class. + + """ + next_family_id = 0 + + def __init__(self, name=None, is_nested=False, weight_func=None, + find_node_at_level=None, effective_measure=None, + level_to_nqpts=None, level_to_exact_to=None, + quad_class=None, quad_class_extra_kwargs=None): + self.family_id = TemplateQuadratureFamily.next_family_id + TemplateQuadratureFamily.next_family_id += 1 + if name is None and quad_class is None: + self.name = "Quadrature" + elif quad_class is not None: + assert issubclass(quad_class, Quadrature) + self.name = quad_class.__name__ + else: + self.name = name + + self.is_nested = is_nested + + if level_to_nqpts is None: + # family of empty rules + self.level_to_nqpts = lambda l: 0 + else: + self.level_to_nqpts = level_to_nqpts + + if level_to_exact_to is None: + self.level_to_exact_to = lambda l: -1 + else: + self.level_to_exact_to = level_to_exact_to + + if weight_func is None: + self.weight_func = lambda x: 1 + assert effective_measure is None or np.allclose(effective_measure, 2) + else: + self.weight_func = weight_func + + if effective_measure is None: + self.effective_measure = 2 + else: + self.effective_measure = effective_measure + + if find_node_at_level is None: + if self.is_nested: + raise ValueError( + "Find_node_at_level must be supplied " + "for nested quadrature rules.") + + # default behavior: + # any node only appears once (strictly non-nested) + # + # this may, however, produce duplicate nodes + + def fnal(nlev, nlid, lev): + if nlev == lev: + return nlid + return -1 + self.find_node_at_level = fnal + else: + self.find_node_at_level = find_node_at_level + + if quad_class is None: + logger.warn("Abstract quad rule %s cannot be intantiated " + "with quad_class unspecified") + self.quad_class = None + self.quad_clsss_extra_kwargs = dict() + else: + assert issubclass(quad_class, Quadrature) + self.quad_class = quad_class + if quad_class_extra_kwargs is None: + self.quad_class_extra_kwargs = dict() + else: + self.quad_class_extra_kwargs = quad_class_extra_kwargs + + def __repr__(self): + return str(self.family_id) + + def __str__(self): + if self.is_nested: + prefix = "Nested" + else: + prefix = "" + return "%s (family_id = %d)" % (prefix + self.name, self.family_id) + + def get_node_info(self, node_index): + r"""Returns some basic information about a node, given its global index, + including the lowest levels containing the node, and its index within + that level. + """ + level = 0 + lev_node_id = node_index + while lev_node_id >= self.level_to_nqpts(level): + lev_node_id -= self.level_to_nqpts(level) + level += 1 + return level, lev_node_id + + def member_node_indices(self, n): + r"""Returns node indices of the n-th level quadrature rule within + this family, in the list :math:`X = [x_i]` of nodes for the whole + family. + + The node list :math:`X` is ordered in a way such that + :math:`X_n = X[s_n:s_n+level_to_nqpts(n)]`, + where :math:`X_n` is the set of nodes for the n-th level quadrature + rule, :math:`s_n = \sum_{i self.quad_desc.family[iaxis].level_to_nqpts(cur_level) + ): + cur_level_capacity = self.quad_desc.family[iaxis].level_to_nqpts( + cur_level) + next_level_nodes = (node_local_indices[iaxis] >= cur_level_capacity) + node_local_indices[iaxis][next_level_nodes] -= cur_level_capacity + node_levels[iaxis][next_level_nodes] += 1 + cur_level += 1 + return node_levels, node_local_indices + + def eliminate_duplicates(self, nodes, weights, tol=None): + """Performs numerical reduction on the generated quadrature rule. + """ + if tol is None: + tol = np.finfo(nodes.dtype).eps * 8 + + if all([fam.is_nested for fam in self.quad_desc.family]): + # nested nodes are ordered such that no numerical + # reduction is needed + return nodes, weights + + nnodes = len(weights) + + # distance matrix, with self-distance = 2 + distance_mat = np.max(np.abs( + nodes[:, :, None] - nodes[:, None, :] + ), axis=0) + 2 * np.diag(np.ones(nnodes)) + + dup_indices_pre = np.where(distance_mat < tol) + if not (dup_indices_pre[0].size > 0): + # no dups detected + return nodes, weights + + dup_indices = np.array( + [ + [dup_indices_pre[0][did], dup_indices_pre[1][did]] + for did in range(dup_indices_pre[0].shape[0]) + if dup_indices_pre[0][did] <= dup_indices_pre[1][did] + ]).transpose() + n_dup_nodes = dup_indices.shape[1] + print("Eliminated duplicate %d nodes" % n_dup_nodes) + + nred_nodes = nnodes - n_dup_nodes + red_indices = [i for i in range(nnodes) if i not in dup_indices[1]] + red_nodes = np.array(nodes[:, red_indices]) + red_weights = np.array(weights[red_indices]) + print("Number of reduced nodes: %d" % nred_nodes) + + return red_nodes, red_weights + + def __call__(self): + dim = self.quad_desc.dim + nnodes = len(self.quad_desc) + quad_classes, quad_classes_extra_kwargs = \ + self.get_quadrature_construction_params() + levels, level_indices = self.get_node_level_local_indices() + nodes = np.zeros([dim, nnodes]) + weights = np.ones(nnodes) + for iaxis, (quad_class, quad_extra_kwargs) in enumerate( + zip(quad_classes, quad_classes_extra_kwargs)): + axis_max_level = np.max(levels[iaxis]) + base_quadratures = [quad_class(l, **quad_extra_kwargs) + for l in range(axis_max_level + 1)] + w_ax = np.zeros(nnodes) + for nid in range(nnodes): + base_quad_rule = base_quadratures[levels[iaxis][nid]] + nodes[iaxis][nid] = base_quad_rule.nodes[level_indices[iaxis][nid]] + for l in range(axis_max_level + 1): + nlid = self.quad_desc.family[iaxis].find_node_at_level( + levels[iaxis][nid], level_indices[iaxis][nid], + l) + if nlid < 0: + # node not in this level + if not np.allclose( + self.quad_desc.weight_coeffs[iaxis, nid, l], 0): + raise RuntimeError("Node %d is not on the level %d " + "in the %d-th axis, but the " + "corresponding weight is nonzero." + % (nid, l, iaxis)) + continue + w_ax[nid] += (self.quad_desc.weight_coeffs[iaxis, nid, l] + * base_quadratures[l].weights[nlid]) + weights *= w_ax + + quad_rule = Quadrature( + *self.eliminate_duplicates(nodes, weights, tol=None)) + + if isinstance(self.quad_desc, QuadratureRingGenerator): + quad_rule.exact_to = self.quad_desc.family[0].level_to_exact_to( + self.quad_desc.max_level[0]) + else: + # TODO: get exactness diagram for non-generators + quad_rule.exact_to = -1 + + return quad_rule + +# }}} End class QuadratureFactory + +# {{{ class QuadratureRingElement + + +class QuadratureRingElement(): + r"""Element (member) of the ring of dim-D quadrature rules that are + generated by 1-D quadrature rules over [-1, 1]^dim. Only one family + of quadrature rules is allowed in each dimension. + + Note that, the ring generated by the 1D quadratures has the "valid" + quadrature rules (that actually integrate functions) as one of its + multiplicative subgroups, identified by the sum of quadrature + weights being certain constant that depends on the weight function + of the generating rules. + + .. attribute :: dim + + Dimension of the template quadrature rule. + + .. attribute :: nnodes + + Number of nodes in the quadrature rule. + + .. attribute :: family + + A tuple of length dim, containing TemplateQuadratureFamily objects + for each axis. + + .. attribute :: max_level + + A tuple of length dim, containing the maximum levels for each axis. + Note that this information is only used for memory management. It + does not necessarily reflect the actual tight bounds of levels + in each dimension. + + .. attribute :: effective_measure + + The sum of quadrature weights. + + .. attribute :: node_indices + + Numpy array of shape (dim, nnodes) + Indicies of the quadrature nodes that contains information + about family and level. + The nodes are sorted in lexicographical order. + + .. attribute :: weight_coeffs + + Numpy array of shape (dim, nnodes, max(max_level) + 1) + Coefficients :math:`a_{d, l}` for evaluating weights at each node from the + weights of the generator elements (1-D rules used for construction). + The weight :math:`w_k = \prod_d (\sum_l a_l w^l_{k, d})`, where `w^l_k` + is the weight of the k-th node inside the l-th level base quadrature rule. + + In effect, weight_coeffs records the full information about how the current + quadrature rule is generated using the base rules (generators). + """ + def __init__(self, dim, family=None, max_level=None, + node_indices=None, weight_coeffs=None, + effective_measure=None): + self.dim = dim + + if family is None: + self.family = (TemplateQuadratureFamily(), ) * dim + elif isinstance(family, (list, tuple)): + if len(family) == 1: + # uniform rules in each axis + self.family = tuple(family) * dim + elif len(family) == dim: + self.family = tuple(family) + else: + raise ValueError() + elif isinstance(family, TemplateQuadratureFamily): + # uniform rules in each axis + self.family = (family, ) * dim + else: + raise ValueError() + + if max_level is None: + self.max_level = (0, ) * dim + else: + if isinstance(max_level, (tuple, list)): + # uniform degrees in each axis + if len(max_level) == 1: + self.max_level = tuple(max_level) * dim + elif len(max_level) == dim: + self.max_level = tuple(max_level) + else: + raise ValueError() + elif isinstance(max_level, (int, np.uint)): + self.max_level = (max_level, ) * dim + else: + raise ValueError() + + if node_indices is None: + self.node_indices = np.array([[], ] * dim) + else: + if node_indices.shape[0] != dim: + raise ValueError("Invalid shape of indices array") + self.node_indices = node_indices + + nnodes = node_indices.shape[1] + degree = max(max_level) + 1 + if weight_coeffs is None: + self.weight_coeffs = np.array([[[], ] * nnodes, ] * dim) + else: + if weight_coeffs.shape != (dim, nnodes, degree): + raise ValueError("Invalid shape of weight coefficients array") + self.weight_coeffs = weight_coeffs + + if effective_measure is None: + self.effective_measure = 2**dim + else: + self.effective_measure = effective_measure + + @property + def nnodes(self): + return self.node_indices.shape[1] + + def __len__(self): + return len(self.node_indices[0]) + + def __add__(self, other): + """Add two rules. The set of nodes are combined, and the weights + for repeated nodes are added together. + """ + if not isinstance(other, QuadratureRingElement): + raise ValueError( + "Cannot add objects of QuadratureRingElement and %s" + % type(other)) + + if self.dim > other.dim: + return self + elif self.dim < other.dim: + return other + + if self.family != other.family: + raise ValueError("Adding different quadrature families (%s, %s) is not " + "currently supported." + % (self.family, other.family)) + + dim = self.dim + family = self.family + max_level = [max(mlthis, mlthat) + for mlthis, mlthat in zip(self.max_level, other.max_level)] + + combined_node_set = set([ + tuple(self.node_indices.transpose()[nid]) + for nid in range(len(self)) + ] + [ + tuple(other.node_indices.transpose()[nid]) + for nid in range(len(other)) + ]) + combined_nodes = sorted(combined_node_set) + nnodes = len(combined_nodes) + node_indices = np.zeros( + (dim, nnodes), + dtype=self.node_indices.dtype) + for nid in range(nnodes): + node_indices[:, nid] = combined_nodes[nid] + + weight_coefs_dict = { + node: np.zeros( + (dim, max(max_level) + 1), dtype=self.weight_coeffs.dtype) + for node in combined_nodes + } + for nid in range(len(self)): + node = tuple(self.node_indices.transpose()[nid]) + weight_coefs_dict[node][:, :max(self.max_level) + 1] \ + += self.weight_coeffs[:, nid, :] + for nid in range(len(other)): + node = tuple(other.node_indices.transpose()[nid]) + weight_coefs_dict[node][:, :max(other.max_level) + 1] \ + += other.weight_coeffs[:, nid, :] + + weight_coeffs = np.zeros( + (dim, nnodes, max(max_level) + 1), + dtype=self.weight_coeffs.dtype) + for nid in range(nnodes): + node = tuple(combined_nodes[nid]) + weight_coeffs[:, nid, :] = weight_coefs_dict[node] + + effective_measure = self.effective_measure + other.effective_measure + + return QuadratureRingElement(dim, family, max_level=max_level, + node_indices=node_indices, weight_coeffs=weight_coeffs, + effective_measure=effective_measure) + + def __mul__(self, other): + """When other is a number, returns a copy of self with weights + rescaled; when other is another TemplateQuadratureRule object, + returns the tensor product of the two rules (self x other). + """ + if np.isscalar(other): + return QuadratureRingElement(self.dim, self.family, self.max_level, + node_indices=self.node_indices.copy(), + weight_coeffs=self.weight_coeffs.copy() * other, + effective_measure=self.effective_measure * other,) + elif not isinstance(other, QuadratureRingElement): + raise ValueError( + "Cannot multiply objects of QuadratureRingElement and %s" + % type(other)) + + dim = self.dim + other.dim + family = self.family + other.family + max_level = self.max_level + other.max_level + + nnodes = len(self) * len(other) + node_indices = np.zeros((dim, nnodes), dtype=self.node_indices.dtype) + weight_coeffs = np.zeros((dim, nnodes, max(max_level) + 1), + dtype=self.weight_coeffs.dtype) + + node_indices[:self.dim, :] = np.repeat( + self.node_indices, repeats=len(other), axis=1) + node_indices[self.dim:, :] = np.tile( + other.node_indices, reps=(1, len(self))) + + weight_coeffs[:self.dim, :, :max(self.max_level) + 1] = np.repeat( + self.weight_coeffs, repeats=len(other), axis=1) + weight_coeffs[self.dim:, :, :max(other.max_level) + 1] = np.tile( + other.weight_coeffs, reps=(1, len(self), 1)) + + effective_measure = self.effective_measure * other.effective_measure + + return QuadratureRingElement(dim, family, max_level=max_level, + node_indices=node_indices, weight_coeffs=weight_coeffs, + effective_measure=effective_measure) + + def __rmul__(self, other): + """This is called only when multiplying a TemplateQuadratureRule + object from the right with a number, in which case the multiplication + is commutative. + """ + if not np.isscalar(other): + raise ValueError( + "Cannot multiply an object of QuadratureRingElement from " + "the right with an object of type %s" + % type(other)) + return self.__mul__(other) + + def __pow__(self, power): + """Shorthand for non-negative integer power of a rule. + """ + assert np.isscalar(power) + assert np.allclose(int(power), power) + assert power >= 0 + if power == 0: + return self.effective_measure + res = 1 + for i in range(int(power)): + res *= self + return res + + def get_exactness_diagram(self): + """Returns the exactness diagram as a list of tuples, where each + tuple contains degrees of polynomials in each axis direction that is + exactly integrated by the quadrature rule. + """ + raise NotImplementedError() + + def get_normalizing_constant(self): + """Returns the normalizing contant. When the weights are rescaled by + this constant, the quadrature rule is guranteed to be valid. + """ + target_measure = 1 + for rule in self.family: + target_measure *= rule.effective_measure + return target_measure / self.effective_measure + + def get_quadrature_rule(self, normalize=True): + """Returns a modepy.quadrature.Quadrature object containing actual + nodes and weights. + + If normalize if True, the weights are rescaled to ensure a valid + quadrature. + """ + quad_rule = QuadratureFactory(self)() + if normalize: + quad_rule.weights *= self.get_normalizing_constant() + return quad_rule + +# }}} End class QuadratureRingElement + +# {{{ class QuadratureRingGenerator + + +class QuadratureRingGenerator(QuadratureRingElement): + """Generator element of the quadrature ring. + By construction, such a quadrature rule belongs to a family of rules + (e.g. Chebyshev-Gauss quadrature, Clenshaw-Curtis quadrature). + + .. attribute :: level + + The index of the quadrature rule within its family, starting from 0. + """ + def __init__(self, family, level): + + assert level >= 0 + self.level = level + + weight_coeffs = np.zeros((1, family.level_to_nqpts(level), level + 1)) + weight_coeffs[0, :, level] = 1 + + QuadratureRingElement.__init__( + self, dim=1, family=family, max_level=(level, ), + node_indices=np.array([family.member_node_indices(level), ]), + weight_coeffs=weight_coeffs, + effective_measure=family.effective_measure) + +# }}} End class QuadratureRingGenerator diff --git a/modepy/symbolic/tensor_product.py b/modepy/symbolic/tensor_product.py new file mode 100644 index 0000000000000000000000000000000000000000..404a032565fbdf459f624ac48f87e7bd1a7c16ff --- /dev/null +++ b/modepy/symbolic/tensor_product.py @@ -0,0 +1,49 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__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 modepy.symbolic as sym + + +class TensorProductQuadratureFactory(): + """Generate tensor product rules. Nodes are ordered lexicographically. + + .. attribute :: generator + + The 1D quadrature rule, instance of + :class:`modepy.symbolic.QuadratureRingGenerator`. + + .. attribute :: dim + + Dimension of tensor product. + + """ + def __init__(self, generator, dim): + assert dim >= 1 + assert isinstance(generator, sym.QuadratureRingGenerator) + self.generator = generator + self.dim = dim + + def __call__(self): + quad = self.generator**self.dim + return quad.get_quadrature_rule() diff --git a/test/test_symbolic.py b/test/test_symbolic.py new file mode 100644 index 0000000000000000000000000000000000000000..2f5c8e792c7fa293724b522a7a89dce5108d4b4e --- /dev/null +++ b/test/test_symbolic.py @@ -0,0 +1,192 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2019 Xiaoyu Wei" + +__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 six.moves import range, zip # noqa + +import numpy as np +import pytest +import modepy.symbolic as sym +import modepy.symbolic.tensor_product as tensor_product + + +@pytest.mark.parametrize('fam', [ + sym.LEGENDRE_GAUSS_FAMILY, + sym.CHEBYSHEV_GAUSS_1_FAMILY, + sym.CHEBYSHEV_GAUSS_2_FAMILY]) +@pytest.mark.parametrize('lev', list(range(5))) +def test_ring_generators(fam, lev): + quad_sym = sym.QuadratureRingGenerator(family=fam, level=lev) + assert quad_sym.dim == 1 + + assert isinstance(quad_sym.family, tuple) + assert len(quad_sym.family) == quad_sym.dim + assert quad_sym.family[0] is fam + + assert isinstance(quad_sym.max_level, tuple) + assert len(quad_sym.max_level) == quad_sym.dim + assert quad_sym.max_level[0] == lev + + quad_rule = quad_sym.get_quadrature_rule(False) + assert len(quad_sym.node_indices) == quad_sym.dim + assert len(quad_sym.node_indices[0]) == len(quad_rule.weights) + print(quad_rule.nodes) + print(quad_rule.weights) + assert np.allclose(sum(quad_rule.weights), fam.effective_measure) + + assert quad_sym.weight_coeffs.shape == ( + quad_sym.dim, len(quad_rule.weights), + max(quad_sym.max_level) + 1) + for l in range(lev): + assert np.allclose(quad_sym.weight_coeffs[0, :, l], 0) + assert np.allclose(quad_sym.weight_coeffs[0, :, lev], 1) + + assert quad_rule.exact_to == quad_sym.family[0].level_to_exact_to(lev) + assert len(quad_rule.weights) == quad_sym.family[0].level_to_nqpts(lev) + assert quad_rule.nodes.shape == (quad_sym.dim, len(quad_rule.weights)) + + +@pytest.mark.parametrize(('fam', 'expect_nested'), [ + (sym.LEGENDRE_GAUSS_FAMILY, False), + (sym.CHEBYSHEV_GAUSS_1_FAMILY, False), + (sym.CHEBYSHEV_GAUSS_2_FAMILY, False), + ]) +def test_family_templates(fam, expect_nested): + assert fam.name[-10:] == 'Quadrature' + assert fam.is_nested == expect_nested + +# {{{ test additions + + +def test_ring_addition_1(): + a = sym.QuadratureRingGenerator(sym.LEGENDRE_GAUSS_FAMILY, 5) + b = sym.QuadratureRingGenerator(sym.LEGENDRE_GAUSS_FAMILY, 10) + c = a + b + + assert c.dim == 1 + assert c.family == a.family + assert c.max_level == (10, ) + + assert np.allclose(c.node_indices[0][:len(a.node_indices[0])], a.node_indices[0]) + assert np.allclose(c.node_indices[0][len(a.node_indices[0]):], b.node_indices[0]) + + for lev in range(11): + if lev == 5: + assert np.allclose(c.weight_coeffs[:, :len(a.node_indices[0]), lev], 1) + assert np.allclose(c.weight_coeffs[:, len(a.node_indices[0]):, lev], 0) + elif lev == 10: + assert np.allclose(c.weight_coeffs[:, :len(a.node_indices[0]), lev], 0) + assert np.allclose(c.weight_coeffs[:, len(a.node_indices[0]):, lev], 1) + else: + assert np.allclose(c.weight_coeffs[0, :, lev], 0) + + quad_rule_unormalized = c.get_quadrature_rule(normalize=False) + assert np.allclose(sum(quad_rule_unormalized.weights), c.effective_measure) + + quad_rule_normalized = c.get_quadrature_rule() + assert np.allclose(sum(quad_rule_normalized.weights), + sym.LEGENDRE_GAUSS_FAMILY.effective_measure) + + +@pytest.mark.xfail +def test_ring_addition_2(): + a = sym.QuadratureRingGenerator(sym.LEGENDRE_GAUSS_FAMILY, 5) + b = sym.QuadratureRingGenerator(sym.CHEBYSHEV_GAUSS_1_FAMILY, 10) + c = a + b # noqa + + +def test_ring_addition_3(): + a = sym.QuadratureRingGenerator(sym.LEGENDRE_GAUSS_FAMILY, 10) + c = a + a + + assert c.dim == 1 + assert c.family == a.family + assert c.max_level == (10, ) + + assert np.allclose(c.node_indices[0], a.node_indices[0]) + for lev in range(11): + if lev == 10: + assert np.allclose(c.weight_coeffs[:, :, lev], 2) + else: + assert np.allclose(c.weight_coeffs[:, :, lev], 0) + +# }}} End test additions + +# {{{ test multiplications + + +def test_ring_multiplication_1(): + a = sym.QuadratureRingGenerator(sym.LEGENDRE_GAUSS_FAMILY, 5) + b = sym.QuadratureRingGenerator(sym.CHEBYSHEV_GAUSS_1_FAMILY, 3) + c = a * b + + assert c.dim == 2 + assert c.family == (sym.LEGENDRE_GAUSS_FAMILY, + sym.CHEBYSHEV_GAUSS_1_FAMILY) + assert c.max_level == (5, 3) + + assert c.nnodes == a.nnodes * b.nnodes + assert np.allclose(c.node_indices[0, ::b.nnodes], a.node_indices) + assert np.allclose(c.node_indices[1, :b.nnodes], b.node_indices) + + assert c.weight_coeffs.shape == (2, c.nnodes, 6) + + +# }}} End test multiplications + +# {{{ test tensor product rules + +def drive_test_tensor_prod_quad(dim, level): + family = sym.LEGENDRE_GAUSS_FAMILY + generator = sym.QuadratureRingGenerator(family, level) + quad_factory = tensor_product.TensorProductQuadratureFactory(generator, dim) + + oned_rule = generator.get_quadrature_rule() + quad_rule = quad_factory() + + assert np.allclose(sum(quad_rule.weights), family.effective_measure**dim) + + assert np.allclose(quad_rule.nodes[-1][:family.level_to_nqpts(level)], + oned_rule.nodes) + + +def test_tensor_prod_quad(): + drive_test_tensor_prod_quad(2, 3) + drive_test_tensor_prod_quad(3, 3) + + +# }}} End test tensor product rules + + +# You can test individual routines by typing +# $ python test_nodes.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker