From 5eabd8b9e13c0d7849938a855f72d9d9cbac4977 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 29 Aug 2019 21:48:13 -0500 Subject: [PATCH 01/25] Structuring sparse grid --- modepy/sparse_gird/__init__.py | 125 +++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 modepy/sparse_gird/__init__.py diff --git a/modepy/sparse_gird/__init__.py b/modepy/sparse_gird/__init__.py new file mode 100644 index 0000000..5e0ce5d --- /dev/null +++ b/modepy/sparse_gird/__init__.py @@ -0,0 +1,125 @@ +__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 + + +class TemplateQuadratureFamily(): + """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 :: level_to_nqpts + + A lambda / function object that maps level to number of quadrature + points. + """ + next_family_id = 0 + + def __init__(self, name=None, is_nested=False): + self.family_id = TemplateQuadratureFamily.next_family_id + TemplateQuadratureFamily.next_family_id += 1 + if name is None: + self.name = "Quadrature" + else: + self.name = name + self.is_nested = is_nested + + def __repr__(self): + return self.family_id + + def __str__(self): + if self.is_nested: + prefix = "Nested" + else: + prefix = "" + return "%s(%d)" % (prefix + self.name, self.family_id) + + +class QuadratureRingElement(): + """Element (member) of the ring of dim-D quadrature rules that are + generated by 1-D quadrature rules over [-1, 1]^dim. + + .. attribute :: dim + + Dimension of the template quadrature rule. + + """ + def __init__(self, dim, q, is_nested, level_to_nqpts): + self.dim = dim + self.is_nested = is_nested + + def __add__(self, other): + """Add two rules. The set of nodes are combined, and the weights + for repeated nodes are added together. + """ + pass + + def __radd__(self, other): + return self.add(other) + + 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. + """ + pass + + 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. + """ + assert np.isscalar(other) + return self.__mul__(other) + + +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 :: family + + A TemplateQuadratureFamily object. + + .. attribute :: level + + The index of the quadrature rule within its family. + + .. attribute :: node_indicies + + Unique indicies of the quadrature nodes that contains information + about family and level. + """ -- GitLab From 53af0927ab59d345f87dd3099118fd2f4b6e6a50 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 30 Aug 2019 18:44:21 -0500 Subject: [PATCH 02/25] Implement quadrature nesting --- modepy/__init__.py | 1 + modepy/sparse_gird/__init__.py | 125 ---------------- modepy/sparse_grid/__init__.py | 263 +++++++++++++++++++++++++++++++++ modepy/sparse_grid/nesting.py | 77 ++++++++++ 4 files changed, 341 insertions(+), 125 deletions(-) delete mode 100644 modepy/sparse_gird/__init__.py create mode 100644 modepy/sparse_grid/__init__.py create mode 100644 modepy/sparse_grid/nesting.py diff --git a/modepy/__init__.py b/modepy/__init__.py index 9fa1d31..797c7f2 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -27,6 +27,7 @@ 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_best_available_basis, grad_simplex_best_available_basis) from modepy.nodes import equidistant_nodes, warp_and_blend_nodes from modepy.matrices import (vandermonde, diff --git a/modepy/sparse_gird/__init__.py b/modepy/sparse_gird/__init__.py deleted file mode 100644 index 5e0ce5d..0000000 --- a/modepy/sparse_gird/__init__.py +++ /dev/null @@ -1,125 +0,0 @@ -__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 - - -class TemplateQuadratureFamily(): - """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 :: level_to_nqpts - - A lambda / function object that maps level to number of quadrature - points. - """ - next_family_id = 0 - - def __init__(self, name=None, is_nested=False): - self.family_id = TemplateQuadratureFamily.next_family_id - TemplateQuadratureFamily.next_family_id += 1 - if name is None: - self.name = "Quadrature" - else: - self.name = name - self.is_nested = is_nested - - def __repr__(self): - return self.family_id - - def __str__(self): - if self.is_nested: - prefix = "Nested" - else: - prefix = "" - return "%s(%d)" % (prefix + self.name, self.family_id) - - -class QuadratureRingElement(): - """Element (member) of the ring of dim-D quadrature rules that are - generated by 1-D quadrature rules over [-1, 1]^dim. - - .. attribute :: dim - - Dimension of the template quadrature rule. - - """ - def __init__(self, dim, q, is_nested, level_to_nqpts): - self.dim = dim - self.is_nested = is_nested - - def __add__(self, other): - """Add two rules. The set of nodes are combined, and the weights - for repeated nodes are added together. - """ - pass - - def __radd__(self, other): - return self.add(other) - - 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. - """ - pass - - 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. - """ - assert np.isscalar(other) - return self.__mul__(other) - - -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 :: family - - A TemplateQuadratureFamily object. - - .. attribute :: level - - The index of the quadrature rule within its family. - - .. attribute :: node_indicies - - Unique indicies of the quadrature nodes that contains information - about family and level. - """ diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py new file mode 100644 index 0000000..b26ec8e --- /dev/null +++ b/modepy/sparse_grid/__init__.py @@ -0,0 +1,263 @@ +__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 + +# {{{ class TemplateQuadratureFamily + + +class TemplateQuadratureFamily(): + """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 :: 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. + """ + next_family_id = 0 + + def __init__(self, name=None, is_nested=False, + level_to_nqpts=None, level_to_exact_to=None): + self.family_id = TemplateQuadratureFamily.next_family_id + TemplateQuadratureFamily.next_family_id += 1 + if name is None: + self.name = "Quadrature" + 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: + # any sane quadrature rule should integrate constants + # exactly even with no known algebraic exactness + self.level_to_exact_to = lambda l: 0 + else: + self.level_to_exact_to = level_to_exact_to + + 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 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= 0 + self.level = level + + QuadratureRingElement.__init__( + self, dim=1, family=family, + node_indices=np.array([family.member_node_indices(level), ]) + ) + +# }}} End class QuadratureRingGenerator diff --git a/modepy/sparse_grid/nesting.py b/modepy/sparse_grid/nesting.py new file mode 100644 index 0000000..d7d631a --- /dev/null +++ b/modepy/sparse_grid/nesting.py @@ -0,0 +1,77 @@ +__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 modepy.sparse_grid import TemplateQuadratureFamily + + +class NestedQuadratureFactory(): + r"""Builds nested quadrature rules based on a "nestable" quadrature, e.g., + Chebyshev-Gauss quadrature of the second kind, where the nodes are all of + the form + + :math:`x^(n)_k = cos(\frac{k}{n} \pi), 1 \leq k \leq n-1`. + + It is clear that when :math:`n` doubles, + :math:`x^(n)_k` and :math:`x^(2n)_{2k}` are the same nodes. + + The nested quadrature rule is then determined by the "initial conditions" and + the nesting structure of the base family. + + NOTE: this class performs the calculation without verifying if the base + family is indeed nestable under the parameters given. + """ + def __init__( + self, base_family, name=None, + initial_index=0, initial_spacing=2, nesting_multiple=2): + """Denote U_n the base family's n-th quadrature rule, the nested rule will + be [U_{n0}, U_{n0 + sp}, U_{n0 + sp * r}, U_{n0 + sp * r^2}, ...]. + + :arg initial_index: the index of the first nested rule (n0). + :arg initial_spacing: the distance between the first two nested rules (sp). + :arg nesting_multiple: the multiple by which the spacing grows (r). + """ + assert isinstance(base_family, TemplateQuadratureFamily) + self.base_family = base_family + if name is None: + self.name = "Nested%s(%d, %d)" % ( + base_family.name, initial_index, nesting_multiple) + else: + self.name = name + self.initial_index = initial_index + self.initial_spacing = initial_spacing + self.nesting_multiple = nesting_multiple + + def get_level_to_nqpts_map(self): + base_map = self.base_family.level_to_nqpts + + def level_to_nqpts(lev): + base_lev = self.initial_index + self.initial_spacing * ( + 2**lev - 1) + return base_map(base_lev) + + return level_to_nqpts + + def __call__(self): + newrule = TemplateQuadratureFamily( + name=self.name, is_nested=True, + level_to_nqpts=self.get_level_to_nqpts_map() + ) + return newrule -- GitLab From 09410a6b03ef818d9f2c28f6efd974032b7ff948 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 31 Aug 2019 11:11:58 -0500 Subject: [PATCH 03/25] Add QuadratureFactory --- modepy/sparse_grid/__init__.py | 90 ++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index b26ec8e..c20865a 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -135,6 +135,96 @@ CHEBYSHEV_GAUSS_2_FAMILY = TemplateQuadratureFamily( # }}} End 1d quadrature families +# {{{ class QuadratureFactory + + +class QuadratureFactory(): + """Maps algebraic quadrature descriptions to actual quadrature rules. + + .. attribute :: quad_desc + + QuadratureRingElement object that describes the quadrature rule + algebraically. + """ + supported_families = [ + LEGENDRE_GAUSS_FAMILY, + CHEBYSHEV_GAUSS_1_FAMILY, CHEBYSHEV_GAUSS_2_FAMILY, + ] + + def __init__(self, quad_desc): + assert isinstance(quad_desc, QuadratureRingElement) + self.quad_desc = quad_desc + + def get_quadrature_construction_params(self): + dim = self.quad_desc.dim + quad_classes = [None, ] * dim + quad_extra_kwargs = [dict(), ] * dim + + for iaxis, qfam in enumerate(self.quad_desc.family): + + if qfam not in self.__class__.supported_families: + raise ValueError("Unsupported quadrature family %s" % qfam) + + if qfam == LEGENDRE_GAUSS_FAMILY: + from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature + quad_classes[iaxis] = LegendreGaussQuadrature + + elif qfam == CHEBYSHEV_GAUSS_1_FAMILY: + from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature + quad_classes[iaxis] = ChebyshevGaussQuadrature + quad_extra_kwargs[iaxis]['kind'] = 1 + + elif qfam == CHEBYSHEV_GAUSS_2_FAMILY: + from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature + quad_classes[iaxis] = ChebyshevGaussQuadrature + quad_extra_kwargs[iaxis]['kind'] = 2 + + else: + raise RuntimeError("This line should not be reached") + + return quad_classes, quad_extra_kwargs + + def get_node_level_local_indices(self): + """Get levels and level-local indices in one axis dimension. + For nested nodes, the search gives the lowest levels containing them. + """ + dim = self.quad_desc.dim + node_local_indices = np.copy(self.quad_desc.node_indices) + node_levels = np.zeros_like(node_local_indices) + for iaxis in range(dim): + cur_level = 0 + while (np.max(node_local_indices[iaxis]) + > 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 + # pylint will give a false positive here + # https://github.com/PyCQA/pylint/issues/2747 + node_levels[iaxis][next_level_nodes] += 1 + cur_level += 1 + return node_levels, node_local_indices + + def __call__(self): + dim = self.quad_desc.dim + nnodes = len(self.quad_desc) + quad_cps = self.get_quadrature_construction_params() + levels, level_indices = self.get_node_level_local_indices() + nodes = np.zeros([dim, nnodes]) + weights = np.zeros(nnodes) + for iaxis, (quad_class, quad_extra_kwargs) in enumerate(quad_cps): + axis_max_level = np.max(levels[iaxis]) + base_quadratures = [quad_class(l, **quad_extra_kwargs) + for l in range(axis_max_level + 1)] + for nid in range(nnodes): + base_quad_rule = base_quadratures[levels[iaxis][nid]] + nodes[iaxis][nid] = base_quad_rule.nodes[level_indices[iaxis][nid]] + # TODO: compute weights from expressions + return Quadrature(nodes, weights) + +# }}} End class QuadratureFactory + # {{{ class QuadratureRingElement -- GitLab From 0584f8d2214f3300307badebeea08cd75c8af98e Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 31 Aug 2019 12:36:10 -0500 Subject: [PATCH 04/25] Add quad elements --- modepy/sparse_grid/__init__.py | 119 ++++++++++++++++++++++++++++----- modepy/sparse_grid/nesting.py | 2 +- 2 files changed, 102 insertions(+), 19 deletions(-) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index c20865a..039e400 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -242,18 +242,33 @@ class QuadratureRingElement(): 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 :: 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_expressions + .. attribute :: weight_coeffs - List of expressions for evaluating weights at each node from the + Numpy array of shape (dim, nnodes, max(max_level)) + 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, node_indices=None): + def __init__(self, dim, family=None, max_level=None, + node_indices=None, weight_coeffs=None): self.dim = dim if family is None: @@ -263,7 +278,7 @@ class QuadratureRingElement(): # uniform rules in each axis self.family = tuple(family) * dim elif len(family) == dim: - self.family = family + self.family = tuple(family) else: raise ValueError() elif isinstance(family, TemplateQuadratureFamily): @@ -272,6 +287,22 @@ class QuadratureRingElement(): 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: @@ -279,6 +310,15 @@ class QuadratureRingElement(): raise ValueError("Invalid shape of indices array") self.node_indices = node_indices + nnodes = node_indices.shape[1] + degree = max(max_level) + 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 + def __len__(self): return len(self.node_indices[0]) @@ -286,16 +326,61 @@ class QuadratureRingElement(): """Add two rules. The set of nodes are combined, and the weights for repeated nodes are added together. """ - pass - - def __radd__(self, other): - return self.add(other) + if self.dim != other.dim: + raise ValueError("Cannot add with different dims (%d, %d)" + % (self.dim, other.dim)) + if self.family != other.family: + raise ValueError("Cannot add with different quadrature families (%s, %s)" + % (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)), 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)] \ + += 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)] \ + += other.weight_coeffs[:, nid, :] + + weight_coeffs = np.zeros( + (dim, nnodes, max(max_level)), + dtype=self.weight_coeffs.dtype) + for nid in range(nnodes): + node = tuple(combined_nodes[nid]) + weight_coeffs[:, nid, :] = weight_coefs_dict[node] + + return self.__class__(self.dim, self.family, max_level=max_level, + node_indices=node_indices, weight_coeffs=None) 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. """ + # TODO pass def __rmul__(self, other): @@ -311,20 +396,14 @@ class QuadratureRingElement(): tuple contains degrees of polynomials in each axis direction that is exactly integrated by the quadrature rule. """ + # TODO return [(0, ) * self.dim] def get_quadrature_rule(self): """Returns a modepy.quadrature.Quadrature object containing actual nodes and weights. """ - nodes = None - weights = None - quad = Quadrature(nodes, weights) - - exact_diag = self.get_exactness_diagram() - quad.exact_to = max([sum(pows) for pows in exact_diag]) - - return quad + return QuadratureFactory(self)() # }}} End class QuadratureRingElement @@ -345,9 +424,13 @@ class QuadratureRingGenerator(QuadratureRingElement): 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, - node_indices=np.array([family.member_node_indices(level), ]) + self, dim=1, family=family, max_level=(level, ), + node_indices=np.array([family.member_node_indices(level), ]), + weight_coeffs=weight_coeffs): ) # }}} End class QuadratureRingGenerator diff --git a/modepy/sparse_grid/nesting.py b/modepy/sparse_grid/nesting.py index d7d631a..fee35b1 100644 --- a/modepy/sparse_grid/nesting.py +++ b/modepy/sparse_grid/nesting.py @@ -22,7 +22,7 @@ THE SOFTWARE. from modepy.sparse_grid import TemplateQuadratureFamily -class NestedQuadratureFactory(): +class NestedQuadratureFamilyFactory(): r"""Builds nested quadrature rules based on a "nestable" quadrature, e.g., Chebyshev-Gauss quadrature of the second kind, where the nodes are all of the form -- GitLab From b4cd4f76856d010e51d5df2e5bbf8f21713446f9 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 31 Aug 2019 13:04:33 -0500 Subject: [PATCH 05/25] Multiply quads --- modepy/sparse_grid/__init__.py | 69 +++++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 17 deletions(-) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index 039e400..3fe39de 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -229,7 +229,7 @@ class QuadratureFactory(): class QuadratureRingElement(): - """Element (member) of the ring of dim-D quadrature rules that are + 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. @@ -326,6 +326,10 @@ class QuadratureRingElement(): """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: raise ValueError("Cannot add with different dims (%d, %d)" % (self.dim, other.dim)) @@ -337,13 +341,13 @@ class QuadratureRingElement(): 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_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( @@ -372,23 +376,55 @@ class QuadratureRingElement(): node = tuple(combined_nodes[nid]) weight_coeffs[:, nid, :] = weight_coefs_dict[node] - return self.__class__(self.dim, self.family, max_level=max_level, - node_indices=node_indices, weight_coeffs=None) + return self.__class__(dim, family, max_level=max_level, + node_indices=node_indices, weight_coeffs=weight_coeffs) 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. + returns the tensor product of the two rules (self x other). """ - # TODO - pass + if np.isscalar(other): + return self.__class__(self.dim, self.family, self.max_level, + node_indices=self.node_indices.copy(), + weight_coeffs=self.weight_coeffs.copy() * 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)), + 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)] = np.repeat( + self.weight_coeffs, repeats=len(other), axis=1) + weight_coeffs[self.dim:, :, :max(other.max_level)] = np.tile( + other.weight_coeffs, reps=(1, len(self), 1)) + + return self.__class__(dim, family, max_level=max_level, + node_indices=node_indices, weight_coeffs=weight_coeffs) 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. """ - assert np.isscalar(other) + 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 get_exactness_diagram(self): @@ -424,13 +460,12 @@ class QuadratureRingGenerator(QuadratureRingElement): assert level >= 0 self.level = level - weight_coeffs=np.zeros((1, family.level_to_nqpts(level), level + 1)) + 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): - ) + weight_coeffs=weight_coeffs) # }}} End class QuadratureRingGenerator -- GitLab From 5ad75974a857916bad2dfcbf0f297cf138028569 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 31 Aug 2019 13:50:13 -0500 Subject: [PATCH 06/25] Redo weight eval --- modepy/sparse_grid/__init__.py | 57 +++++++++++++++++++++++++++++++--- 1 file changed, 53 insertions(+), 4 deletions(-) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index 3fe39de..f8423c4 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -43,6 +43,12 @@ class TemplateQuadratureFamily(): A bool value indicating whether the quadrature rules can be nested. + .. 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. + .. attribute :: level_to_nqpts A lambda / function object that maps level to number of quadrature @@ -56,6 +62,7 @@ class TemplateQuadratureFamily(): next_family_id = 0 def __init__(self, name=None, is_nested=False, + find_node_at_level=None, level_to_nqpts=None, level_to_exact_to=None): self.family_id = TemplateQuadratureFamily.next_family_id TemplateQuadratureFamily.next_family_id += 1 @@ -78,6 +85,24 @@ class TemplateQuadratureFamily(): else: self.level_to_exact_to = level_to_exact_to + 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) + + def fnal(nid, lev): + level_first, level_id = self.get_node_info(nid) + if level_first == lev: + return level_id + return -1 + self.find_node_at_level = fnal + else: + self.find_node_at_level = find_node_at_level + def __repr__(self): return str(self.family_id) @@ -88,6 +113,18 @@ class TemplateQuadratureFamily(): 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 @@ -187,6 +224,9 @@ class QuadratureFactory(): def get_node_level_local_indices(self): """Get levels and level-local indices in one axis dimension. For nested nodes, the search gives the lowest levels containing them. + + The results should be in line with TemplateQuadratureFamily.get_node_info(), + but in this function the calculation is done in batches. """ dim = self.quad_desc.dim node_local_indices = np.copy(self.quad_desc.node_indices) @@ -200,8 +240,6 @@ class QuadratureFactory(): cur_level) next_level_nodes = (node_local_indices[iaxis] >= cur_level_capacity) node_local_indices[iaxis][next_level_nodes] -= cur_level_capacity - # pylint will give a false positive here - # https://github.com/PyCQA/pylint/issues/2747 node_levels[iaxis][next_level_nodes] += 1 cur_level += 1 return node_levels, node_local_indices @@ -212,15 +250,26 @@ class QuadratureFactory(): quad_cps = self.get_quadrature_construction_params() levels, level_indices = self.get_node_level_local_indices() nodes = np.zeros([dim, nnodes]) - weights = np.zeros(nnodes) + weights = np.ones(nnodes) for iaxis, (quad_class, quad_extra_kwargs) in enumerate(quad_cps): 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]] - # TODO: compute weights from expressions + for l in range(axis_max_level + 1): + nlid = self.quad_desc.find_node_at_level(nid, l) + if nlid < 0: + # node not in this level + assert np.allclose( + self.quad_desc.weight_coeffs[iaxis, nid, l], 0) + continue + w_ax[nid] += (self.quad_desc.weight_coeffs[iaxis, nid, l] + * base_quadratures[l].weights[nlid]) + weights *= w_ax + return Quadrature(nodes, weights) # }}} End class QuadratureFactory -- GitLab From 0a26e32a44369f1d2a0aa23bbf8a092d616334ae Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 31 Aug 2019 15:38:12 -0500 Subject: [PATCH 07/25] Bug fixes --- modepy/sparse_grid/__init__.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index f8423c4..01633c5 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -247,11 +247,13 @@ class QuadratureFactory(): def __call__(self): dim = self.quad_desc.dim nnodes = len(self.quad_desc) - quad_cps = self.get_quadrature_construction_params() + 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(quad_cps): + 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)] @@ -260,11 +262,16 @@ class QuadratureFactory(): 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.find_node_at_level(nid, l) + nlid = self.quad_desc.family[iaxis].find_node_at_level( + self.quad_desc.node_indices[iaxis][nid], l) if nlid < 0: # node not in this level - assert np.allclose( - self.quad_desc.weight_coeffs[iaxis, nid, l], 0) + 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]) @@ -339,7 +346,7 @@ class QuadratureRingElement(): if max_level is None: self.max_level = (0, ) * dim else: - if isinstance(max_level, [tuple, list]): + if isinstance(max_level, (tuple, list)): # uniform degrees in each axis if len(max_level) == 1: self.max_level = tuple(max_level) * dim @@ -347,7 +354,7 @@ class QuadratureRingElement(): self.max_level = tuple(max_level) else: raise ValueError() - elif isinstance(max_level, [int, np.uint]): + elif isinstance(max_level, (int, np.uint)): self.max_level = (max_level, ) * dim else: raise ValueError() @@ -360,7 +367,7 @@ class QuadratureRingElement(): self.node_indices = node_indices nnodes = node_indices.shape[1] - degree = max(max_level) + degree = max(max_level) + 1 if weight_coeffs is None: self.weight_coeffs = np.array([[[], ] * nnodes, ] * dim) else: -- GitLab From 8187c0a2b922780416501039c83fb256d45c0aa6 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 1 Sep 2019 11:51:37 -0500 Subject: [PATCH 08/25] Fix QuadratureFactory --- modepy/sparse_grid/__init__.py | 86 ++++++++++++++++++++++++++++------ modepy/sparse_grid/nesting.py | 46 ++++++++++++++++-- 2 files changed, 112 insertions(+), 20 deletions(-) diff --git a/modepy/sparse_grid/__init__.py b/modepy/sparse_grid/__init__.py index 01633c5..0040ccb 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/sparse_grid/__init__.py @@ -48,6 +48,8 @@ class TemplateQuadratureFamily(): 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 @@ -62,8 +64,8 @@ class TemplateQuadratureFamily(): next_family_id = 0 def __init__(self, name=None, is_nested=False, - find_node_at_level=None, - level_to_nqpts=None, level_to_exact_to=None): + find_node_at_level=None, + level_to_nqpts=None, level_to_exact_to=None): self.family_id = TemplateQuadratureFamily.next_family_id TemplateQuadratureFamily.next_family_id += 1 if name is None: @@ -88,16 +90,17 @@ class TemplateQuadratureFamily(): if find_node_at_level is None: if self.is_nested: raise ValueError( - "Find_node_at_level must be supplied " - "for nested quadrature rules.") + "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(nid, lev): - level_first, level_id = self.get_node_info(nid) - if level_first == lev: - return level_id + def fnal(nlev, nlid, lev): + if nlev == lev: + return nlid return -1 self.find_node_at_level = fnal else: @@ -123,7 +126,7 @@ class TemplateQuadratureFamily(): while lev_node_id >= self.level_to_nqpts(level): lev_node_id -= self.level_to_nqpts(level) level += 1 - return level, lev_node_id, + return level, lev_node_id def member_node_indices(self, n): r"""Returns node indices of the n-th level quadrature rule within @@ -176,7 +179,17 @@ CHEBYSHEV_GAUSS_2_FAMILY = TemplateQuadratureFamily( class QuadratureFactory(): - """Maps algebraic quadrature descriptions to actual quadrature rules. + r"""Maps algebraic quadrature descriptions to actual quadrature rules. + + When the base family is nested, the algebraic description is guaranteed to + have no duplicate nodes. + + When the base family is not nested, the algebraic description may contain + duplicate nodes that are not known to the symbolic engine, e.g., + when :math:`x_k^n = cos((k+1)/n \pi)`, :math:`x_{k}^n = x_{2k+1}^{2n}`. + Then a numerical elimination is performed on the generated quadrature points, + where nodes close to each other as identified by a given tolerance + are combined. .. attribute :: quad_desc @@ -244,6 +257,46 @@ class QuadratureFactory(): 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]: + # 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) @@ -263,7 +316,8 @@ class QuadratureFactory(): 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( - self.quad_desc.node_indices[iaxis][nid], l) + levels[iaxis][nid], level_indices[iaxis][nid], + l) if nlid < 0: # node not in this level if not np.allclose( @@ -277,7 +331,9 @@ class QuadratureFactory(): * base_quadratures[l].weights[nlid]) weights *= w_ax - return Quadrature(nodes, weights) + return Quadrature( + *self.eliminate_duplicates(nodes, weights, tol=None) + ) # }}} End class QuadratureFactory @@ -324,7 +380,7 @@ class QuadratureRingElement(): quadrature rule is generated using the base rules (generators). """ def __init__(self, dim, family=None, max_level=None, - node_indices=None, weight_coeffs=None): + node_indices=None, weight_coeffs=None): self.dim = dim if family is None: @@ -384,8 +440,8 @@ class QuadratureRingElement(): """ if not isinstance(other, QuadratureRingElement): raise ValueError( - "Cannot add objects of QuadratureRingElement and %s" - % type(other)) + "Cannot add objects of QuadratureRingElement and %s" + % type(other)) if self.dim != other.dim: raise ValueError("Cannot add with different dims (%d, %d)" % (self.dim, other.dim)) diff --git a/modepy/sparse_grid/nesting.py b/modepy/sparse_grid/nesting.py index fee35b1..301dc7e 100644 --- a/modepy/sparse_grid/nesting.py +++ b/modepy/sparse_grid/nesting.py @@ -40,13 +40,19 @@ class NestedQuadratureFamilyFactory(): """ def __init__( self, base_family, name=None, - initial_index=0, initial_spacing=2, nesting_multiple=2): - """Denote U_n the base family's n-th quadrature rule, the nested rule will + initial_index=0, initial_spacing=2, nesting_multiple=2, + padding=1): + r"""Denote U_n the base family's n-th quadrature rule, the nested rule will be [U_{n0}, U_{n0 + sp}, U_{n0 + sp * r}, U_{n0 + sp * r^2}, ...]. :arg initial_index: the index of the first nested rule (n0). :arg initial_spacing: the distance between the first two nested rules (sp). :arg nesting_multiple: the multiple by which the spacing grows (r). + :arg padding: the padding parameter associated the following invariant: + :math:`x_i^n = x_{(i + padding) * r - padding}^{n+1}`, + where :math:`0 \leq i < nqpts(n)`. + In other words, padding is used for computing k given i, e.g., + with Chebyshev-Gauss nodes `padding = 1`, since k = i + 1. """ assert isinstance(base_family, TemplateQuadratureFamily) self.base_family = base_family @@ -58,20 +64,50 @@ class NestedQuadratureFamilyFactory(): self.initial_index = initial_index self.initial_spacing = initial_spacing self.nesting_multiple = nesting_multiple + self.padding = padding + + def get_base_level(self, lev): + """Computes the level in the base family given the level + in the nested sub-family. + """ + return self.initial_index + self.initial_spacing * ( + 2**lev - 1) def get_level_to_nqpts_map(self): base_map = self.base_family.level_to_nqpts def level_to_nqpts(lev): - base_lev = self.initial_index + self.initial_spacing * ( - 2**lev - 1) + base_lev = self.get_base_level(lev) return base_map(base_lev) return level_to_nqpts + def get_node_at_level_finder(self): + """Returns find_node_at_level needed for the TemplateQuadratureFamily + construction. + + Currently it makes the following assumptions on the base_family: + + 1. base_family's level_to_nqpts map is linear with unit slope + (:math:`f(n) = n + c`). + + 2. base_family is uniformly nested (as is the case with + (:math:`x_k = cos(k / n)`), i.e. new nodes and nodes from + the previous levels are staggered. + """ + def find_node_at_level(nlev, nlid, lev): + if nlev < lev: + return -1 + k = self.padding + nlid + k_lev = k * self.nesting_multiple**(lev - nlev) + return k_lev - self.padding + + return find_node_at_level + def __call__(self): newrule = TemplateQuadratureFamily( name=self.name, is_nested=True, - level_to_nqpts=self.get_level_to_nqpts_map() + level_to_nqpts=self.get_level_to_nqpts_map(), + find_node_at_level=self.get_node_at_level_finder(), ) return newrule -- GitLab From eaacfbf375a4cff1d557faed50db898474776efb Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 15 Sep 2019 13:08:02 -0500 Subject: [PATCH 09/25] Minor fixes --- modepy/{sparse_grid => symbolic}/__init__.py | 8 ++++---- modepy/{sparse_grid => symbolic}/nesting.py | 0 2 files changed, 4 insertions(+), 4 deletions(-) rename modepy/{sparse_grid => symbolic}/__init__.py (98%) rename modepy/{sparse_grid => symbolic}/nesting.py (100%) diff --git a/modepy/sparse_grid/__init__.py b/modepy/symbolic/__init__.py similarity index 98% rename from modepy/sparse_grid/__init__.py rename to modepy/symbolic/__init__.py index 0040ccb..4874729 100644 --- a/modepy/sparse_grid/__init__.py +++ b/modepy/symbolic/__init__.py @@ -511,7 +511,7 @@ class QuadratureRingElement(): 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)), + weight_coeffs = np.zeros((dim, nnodes, max(max_level) + 1), dtype=self.weight_coeffs.dtype) node_indices[:self.dim, :] = np.repeat( @@ -519,12 +519,12 @@ class QuadratureRingElement(): node_indices[self.dim:, :] = np.tile( other.node_indices, reps=(1, len(self))) - weight_coeffs[:self.dim, :, :max(self.max_level)] = np.repeat( + 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)] = np.tile( + weight_coeffs[self.dim:, :, :max(other.max_level) + 1] = np.tile( other.weight_coeffs, reps=(1, len(self), 1)) - return self.__class__(dim, family, max_level=max_level, + return QuadratureRingElement(dim, family, max_level=max_level, node_indices=node_indices, weight_coeffs=weight_coeffs) def __rmul__(self, other): diff --git a/modepy/sparse_grid/nesting.py b/modepy/symbolic/nesting.py similarity index 100% rename from modepy/sparse_grid/nesting.py rename to modepy/symbolic/nesting.py -- GitLab From 07e50964390ee206851eaa9bebe128679c581af4 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 19 Sep 2019 13:29:29 -0500 Subject: [PATCH 10/25] Add tests for generators --- modepy/symbolic/__init__.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index 4874729..0bd1d76 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -147,10 +147,13 @@ class TemplateQuadratureFamily(): return np.arange(sn, sn + self.level_to_nqpts(n)) + # }}} End class TemplateQuadratureFamily # {{{ 1d quadrature families +# Naming convention (for the name attribute): +# always end with the class name for the concrete quadrature class. # LegendreGaussQuadrature LEGENDRE_GAUSS_FAMILY = TemplateQuadratureFamily( @@ -276,7 +279,7 @@ class QuadratureFactory(): ), axis=0) + 2 * np.diag(np.ones(nnodes)) dup_indices_pre = np.where(distance_mat < tol) - if not dup_indices_pre[0]: + if not (dup_indices_pre[0].size > 0): # no dups detected return nodes, weights @@ -370,7 +373,7 @@ class QuadratureRingElement(): .. attribute :: weight_coeffs - Numpy array of shape (dim, nnodes, max(max_level)) + 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` @@ -544,8 +547,7 @@ class QuadratureRingElement(): tuple contains degrees of polynomials in each axis direction that is exactly integrated by the quadrature rule. """ - # TODO - return [(0, ) * self.dim] + raise NotImplementedError() def get_quadrature_rule(self): """Returns a modepy.quadrature.Quadrature object containing actual @@ -565,7 +567,7 @@ class QuadratureRingGenerator(QuadratureRingElement): .. attribute :: level - The index of the quadrature rule within its family. + The index of the quadrature rule within its family, starting from 0. """ def __init__(self, family, level): -- GitLab From b00ffd59193ffe414fc70593f42e494643e1a5d2 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 19 Sep 2019 14:15:34 -0500 Subject: [PATCH 11/25] Test additions --- modepy/symbolic/__init__.py | 37 +++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index 0bd1d76..08c1ce2 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -81,9 +81,7 @@ class TemplateQuadratureFamily(): self.level_to_nqpts = level_to_nqpts if level_to_exact_to is None: - # any sane quadrature rule should integrate constants - # exactly even with no known algebraic exactness - self.level_to_exact_to = lambda l: 0 + self.level_to_exact_to = lambda l: -1 else: self.level_to_exact_to = level_to_exact_to @@ -159,21 +157,24 @@ class TemplateQuadratureFamily(): LEGENDRE_GAUSS_FAMILY = TemplateQuadratureFamily( name='LegendreGaussQuadrature', is_nested=False, - level_to_nqpts=lambda n: n + 1 + level_to_nqpts=lambda n: n + 1, + level_to_exact_to=lambda n: 2 * n + 1, ) # ChebyshevGaussQuadrature(kind=1) CHEBYSHEV_GAUSS_1_FAMILY = TemplateQuadratureFamily( name='FirstKindChebyshevGaussQuadrature', is_nested=False, - level_to_nqpts=lambda n: n + 1 + level_to_nqpts=lambda n: n + 1, + level_to_exact_to=lambda n: 2 * n + 1, ) # ChebyshevGaussQuadrature(kind=2) CHEBYSHEV_GAUSS_2_FAMILY = TemplateQuadratureFamily( name='SecondKindChebyshevGaussQuadrature', is_nested=False, - level_to_nqpts=lambda n: n + 1 + level_to_nqpts=lambda n: n + 1, + level_to_exact_to=lambda n: 2 * n + 1, ) # }}} End 1d quadrature families @@ -334,9 +335,17 @@ class QuadratureFactory(): * base_quadratures[l].weights[nlid]) weights *= w_ax - return Quadrature( - *self.eliminate_duplicates(nodes, weights, tol=None) - ) + 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 @@ -472,26 +481,26 @@ class QuadratureRingElement(): node_indices[:, nid] = combined_nodes[nid] weight_coefs_dict = { - node: np.zeros((dim, max(max_level)), dtype=self.weight_coeffs.dtype) + 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)] \ + 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)] \ + weight_coefs_dict[node][:, :max(other.max_level) + 1] \ += other.weight_coeffs[:, nid, :] weight_coeffs = np.zeros( - (dim, nnodes, max(max_level)), + (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] - return self.__class__(dim, family, max_level=max_level, + return QuadratureRingElement(dim, family, max_level=max_level, node_indices=node_indices, weight_coeffs=weight_coeffs) def __mul__(self, other): -- GitLab From 7f5723b80dcc5b19240b825344a5d5a22af91f58 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 19 Sep 2019 14:32:19 -0500 Subject: [PATCH 12/25] Fixes --- modepy/symbolic/__init__.py | 3 +- test/test_symbolic.py | 137 ++++++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+), 1 deletion(-) create mode 100644 test/test_symbolic.py diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index 08c1ce2..c651a77 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -481,7 +481,8 @@ class QuadratureRingElement(): node_indices[:, nid] = combined_nodes[nid] weight_coefs_dict = { - node: np.zeros((dim, max(max_level) + 1), dtype=self.weight_coeffs.dtype) + node: np.zeros( + (dim, max(max_level) + 1), dtype=self.weight_coeffs.dtype) for node in combined_nodes } for nid in range(len(self)): diff --git a/test/test_symbolic.py b/test/test_symbolic.py new file mode 100644 index 0000000..1eabfbc --- /dev/null +++ b/test/test_symbolic.py @@ -0,0 +1,137 @@ +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 + + +@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() + assert len(quad_sym.node_indices) == quad_sym.dim + assert len(quad_sym.node_indices[0]) == len(quad_rule.weights) + + 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 in sym.QuadratureFactory.supported_families + 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) + + +@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 + + +# 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 -- GitLab From ae343db5647591bbe25c17f5b5421dbd767a3afc Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 19 Sep 2019 16:55:16 -0500 Subject: [PATCH 13/25] Use effective_measure to track normalizers --- modepy/symbolic/__init__.py | 77 +++++++++++++++++++++++++++++++++---- test/test_symbolic.py | 34 +++++++++++++++- 2 files changed, 102 insertions(+), 9 deletions(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index c651a77..de73ea9 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -43,6 +43,10 @@ class TemplateQuadratureFamily(): A bool value indicating whether the quadrature rules can be nested. + .. attribute :: effective_measure + + The sum of quadrature weights (depends on the weight function). + .. attribute :: find_node_at_level A lambda / function object that when given an node and a level, @@ -64,7 +68,7 @@ class TemplateQuadratureFamily(): next_family_id = 0 def __init__(self, name=None, is_nested=False, - find_node_at_level=None, + find_node_at_level=None, effective_measure=None, level_to_nqpts=None, level_to_exact_to=None): self.family_id = TemplateQuadratureFamily.next_family_id TemplateQuadratureFamily.next_family_id += 1 @@ -85,6 +89,11 @@ class TemplateQuadratureFamily(): else: self.level_to_exact_to = level_to_exact_to + 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( @@ -157,6 +166,7 @@ class TemplateQuadratureFamily(): LEGENDRE_GAUSS_FAMILY = TemplateQuadratureFamily( name='LegendreGaussQuadrature', is_nested=False, + effective_measure=2, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, ) @@ -165,6 +175,7 @@ LEGENDRE_GAUSS_FAMILY = TemplateQuadratureFamily( CHEBYSHEV_GAUSS_1_FAMILY = TemplateQuadratureFamily( name='FirstKindChebyshevGaussQuadrature', is_nested=False, + effective_measure=np.pi, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, ) @@ -173,6 +184,7 @@ CHEBYSHEV_GAUSS_1_FAMILY = TemplateQuadratureFamily( CHEBYSHEV_GAUSS_2_FAMILY = TemplateQuadratureFamily( name='SecondKindChebyshevGaussQuadrature', is_nested=False, + effective_measure=np.pi / 2, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, ) @@ -227,11 +239,13 @@ class QuadratureFactory(): from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature quad_classes[iaxis] = ChebyshevGaussQuadrature quad_extra_kwargs[iaxis]['kind'] = 1 + quad_extra_kwargs[iaxis]['backend'] = 'scipy' elif qfam == CHEBYSHEV_GAUSS_2_FAMILY: from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature quad_classes[iaxis] = ChebyshevGaussQuadrature quad_extra_kwargs[iaxis]['kind'] = 2 + quad_extra_kwargs[iaxis]['backend'] = 'builtin' else: raise RuntimeError("This line should not be reached") @@ -357,10 +371,20 @@ class QuadratureRingElement(): 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 @@ -373,6 +397,10 @@ class QuadratureRingElement(): 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) @@ -392,7 +420,8 @@ class QuadratureRingElement(): quadrature rule is generated using the base rules (generators). """ def __init__(self, dim, family=None, max_level=None, - node_indices=None, weight_coeffs=None): + node_indices=None, weight_coeffs=None, + effective_measure=None): self.dim = dim if family is None: @@ -443,6 +472,15 @@ class QuadratureRingElement(): 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]) @@ -501,8 +539,11 @@ class QuadratureRingElement(): 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) + 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 @@ -512,7 +553,8 @@ class QuadratureRingElement(): if np.isscalar(other): return self.__class__(self.dim, self.family, self.max_level, node_indices=self.node_indices.copy(), - weight_coeffs=self.weight_coeffs.copy() * other) + 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" @@ -537,8 +579,11 @@ class QuadratureRingElement(): 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) + node_indices=node_indices, weight_coeffs=weight_coeffs, + effective_measure=effective_measure) def __rmul__(self, other): """This is called only when multiplying a TemplateQuadratureRule @@ -559,11 +604,26 @@ class QuadratureRingElement(): """ raise NotImplementedError() - def get_quadrature_rule(self): + 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. """ - return QuadratureFactory(self)() + quad_rule = QuadratureFactory(self)() + if normalize: + quad_rule.weights *= self.get_normalizing_constant() + return quad_rule # }}} End class QuadratureRingElement @@ -590,6 +650,7 @@ class QuadratureRingGenerator(QuadratureRingElement): QuadratureRingElement.__init__( self, dim=1, family=family, max_level=(level, ), node_indices=np.array([family.member_node_indices(level), ]), - weight_coeffs=weight_coeffs) + weight_coeffs=weight_coeffs, + effective_measure=family.effective_measure) # }}} End class QuadratureRingGenerator diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 1eabfbc..7c03aba 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -46,9 +46,12 @@ def test_ring_generators(fam, lev): assert len(quad_sym.max_level) == quad_sym.dim assert quad_sym.max_level[0] == lev - quad_rule = quad_sym.get_quadrature_rule() + 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), @@ -97,6 +100,13 @@ def test_ring_addition_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(): @@ -122,6 +132,28 @@ def test_ring_addition_3(): # }}} 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 + # You can test individual routines by typing # $ python test_nodes.py 'test_routine()' -- GitLab From b03a56b4f5424b189e7f73db4b95e68094559398 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 19 Sep 2019 17:08:30 -0500 Subject: [PATCH 14/25] Py2 fix --- modepy/symbolic/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index de73ea9..fdaafb4 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -1,3 +1,5 @@ +from __future__ import division + __copyright__ = "Copyright (C) 2019 Xiaoyu Wei" __license__ = """ -- GitLab From 1feda0cebefcd6dc99ac8f3d9692e5ec70c76f03 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 20 Sep 2019 09:34:23 -0500 Subject: [PATCH 15/25] Fix docstring --- modepy/symbolic/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index fdaafb4..f04577d 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -47,7 +47,9 @@ class TemplateQuadratureFamily(): .. attribute :: effective_measure - The sum of quadrature weights (depends on the weight function). + 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 -- GitLab From ee67d40a15f66b1d2cd498f834f9d22ab13d263c Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 20 Sep 2019 09:37:02 -0500 Subject: [PATCH 16/25] Save nesting for later --- modepy/symbolic/nesting.py | 113 ------------------------------------- 1 file changed, 113 deletions(-) delete mode 100644 modepy/symbolic/nesting.py diff --git a/modepy/symbolic/nesting.py b/modepy/symbolic/nesting.py deleted file mode 100644 index 301dc7e..0000000 --- a/modepy/symbolic/nesting.py +++ /dev/null @@ -1,113 +0,0 @@ -__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 modepy.sparse_grid import TemplateQuadratureFamily - - -class NestedQuadratureFamilyFactory(): - r"""Builds nested quadrature rules based on a "nestable" quadrature, e.g., - Chebyshev-Gauss quadrature of the second kind, where the nodes are all of - the form - - :math:`x^(n)_k = cos(\frac{k}{n} \pi), 1 \leq k \leq n-1`. - - It is clear that when :math:`n` doubles, - :math:`x^(n)_k` and :math:`x^(2n)_{2k}` are the same nodes. - - The nested quadrature rule is then determined by the "initial conditions" and - the nesting structure of the base family. - - NOTE: this class performs the calculation without verifying if the base - family is indeed nestable under the parameters given. - """ - def __init__( - self, base_family, name=None, - initial_index=0, initial_spacing=2, nesting_multiple=2, - padding=1): - r"""Denote U_n the base family's n-th quadrature rule, the nested rule will - be [U_{n0}, U_{n0 + sp}, U_{n0 + sp * r}, U_{n0 + sp * r^2}, ...]. - - :arg initial_index: the index of the first nested rule (n0). - :arg initial_spacing: the distance between the first two nested rules (sp). - :arg nesting_multiple: the multiple by which the spacing grows (r). - :arg padding: the padding parameter associated the following invariant: - :math:`x_i^n = x_{(i + padding) * r - padding}^{n+1}`, - where :math:`0 \leq i < nqpts(n)`. - In other words, padding is used for computing k given i, e.g., - with Chebyshev-Gauss nodes `padding = 1`, since k = i + 1. - """ - assert isinstance(base_family, TemplateQuadratureFamily) - self.base_family = base_family - if name is None: - self.name = "Nested%s(%d, %d)" % ( - base_family.name, initial_index, nesting_multiple) - else: - self.name = name - self.initial_index = initial_index - self.initial_spacing = initial_spacing - self.nesting_multiple = nesting_multiple - self.padding = padding - - def get_base_level(self, lev): - """Computes the level in the base family given the level - in the nested sub-family. - """ - return self.initial_index + self.initial_spacing * ( - 2**lev - 1) - - def get_level_to_nqpts_map(self): - base_map = self.base_family.level_to_nqpts - - def level_to_nqpts(lev): - base_lev = self.get_base_level(lev) - return base_map(base_lev) - - return level_to_nqpts - - def get_node_at_level_finder(self): - """Returns find_node_at_level needed for the TemplateQuadratureFamily - construction. - - Currently it makes the following assumptions on the base_family: - - 1. base_family's level_to_nqpts map is linear with unit slope - (:math:`f(n) = n + c`). - - 2. base_family is uniformly nested (as is the case with - (:math:`x_k = cos(k / n)`), i.e. new nodes and nodes from - the previous levels are staggered. - """ - def find_node_at_level(nlev, nlid, lev): - if nlev < lev: - return -1 - k = self.padding + nlid - k_lev = k * self.nesting_multiple**(lev - nlev) - return k_lev - self.padding - - return find_node_at_level - - def __call__(self): - newrule = TemplateQuadratureFamily( - name=self.name, is_nested=True, - level_to_nqpts=self.get_level_to_nqpts_map(), - find_node_at_level=self.get_node_at_level_finder(), - ) - return newrule -- GitLab From b1a41c275f8b84c1447c1a3d8f3d0de499fe6058 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 20 Sep 2019 09:37:35 -0500 Subject: [PATCH 17/25] Fix docstring --- modepy/symbolic/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index f04577d..676294f 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -29,7 +29,7 @@ from modepy.quadrature import Quadrature class TemplateQuadratureFamily(): - """Family of 1-D quadrature rules over [-1, 1], indexed by an integral + r"""Family of 1-D quadrature rules over [-1, 1], indexed by an integral parameter *level*. .. attribute :: family_id -- GitLab From 2928931ce24ffffbe02ef13f595949070444028f Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 20 Sep 2019 17:51:25 -0500 Subject: [PATCH 18/25] Add documentation --- doc/index.rst | 1 + modepy/symbolic/__init__.py | 116 ++++++++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+) diff --git a/doc/index.rst b/doc/index.rst index 7dfe6cd..91e2e26 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -45,6 +45,7 @@ Contents modes nodes quadrature + symbolic tools misc diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index 676294f..0cb63ad 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -25,6 +25,122 @@ THE SOFTWARE. import numpy as np from modepy.quadrature import Quadrature +__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` as the dimension. + +Domains + We restrict our discussion to :math:`d`-dimensional template + domains :math:`\Omega_d = [-1, 1]^d`. + +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.}\}. + +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}`. + +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\}`. + + If the set of nodes and + weights is empty, we call it an empty rule, and denote it by + :math:`0`, as it will be the additive identity for any group of + quadrature rules under the addition operation defined below. + 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}`. + + +Algebra of quadrature rules +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Addition + Given a quadrature chain with two elements :math:`a, b`, define + + .. math:: + + (a+b)(f) := (\sum_{i=1}^{n_a} w_{a,i} ev_{x_{a,i}} + + \sum_{i=1}^{n_b} w_{b,i} ev_{x_{b,j}})(f\omega), + \forall f\in V_{d, \omega}. + + Alternatively, if :math:`a, b` have different dimensions, define + + .. math:: + + a + b := argmax_{quad\in\{a, b\}}dim(quad). + + The rationale of this definition is that, the rule with lower + dimension is tensored with empty rules to "match" the higher + dimension before added together. + + Obviously, in either case the addition is commutative + :math:`a + b = b + a`. + + Note that addition of two rules with same dimension but different + weight functions is not allowed. + +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,i}})(f\omega), + \forall f\in V_{d_a + d_b, \omega}, + + where :math:`\omega = \omega_a \otimes \omega_b`. Obviously, the + multiplication of two quadrature rules is not commutative. + +""" + # {{{ class TemplateQuadratureFamily -- GitLab From 7b0d8f3475e8cb9cbae898afcb613084ca7caa56 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 20 Sep 2019 17:55:25 -0500 Subject: [PATCH 19/25] Add rst doc source --- doc/symbolic.rst | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 doc/symbolic.rst diff --git a/doc/symbolic.rst b/doc/symbolic.rst new file mode 100644 index 0000000..becf18a --- /dev/null +++ b/doc/symbolic.rst @@ -0,0 +1,6 @@ +Symbolic +======== + +.. automodule:: modepy.symbolic + +.. vim: sw=4 -- GitLab From f633f6ff446f066a0853404faa3481e601947d9e Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Sep 2019 13:50:17 -0500 Subject: [PATCH 20/25] Make addition and multiplication definitions self-consistent --- modepy/symbolic/__init__.py | 166 +++++++++++++++++++++++++++++++----- 1 file changed, 146 insertions(+), 20 deletions(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index 0cb63ad..b6f3b81 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -34,12 +34,15 @@ Terminologies ^^^^^^^^^^^^^ Dimension - Denote :math:`d` as the 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 @@ -54,11 +57,18 @@ Function spaces \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)`. @@ -74,11 +84,18 @@ Quadrature 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\}`. - If the set of nodes and - weights is empty, we call it an empty rule, and denote it by - :math:`0`, as it will be the additive identity for any group of - quadrature rules under the addition operation defined below. - The empty rule can be treated as any dimensional rule + 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. @@ -97,35 +114,102 @@ Consistency :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) := + (\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 a quadrature chain with two elements :math:`a, b`, define + Given a quadrature chain in :math:`d` dimensions with two elements + :math:`a, b`, the addition is defined as .. math:: - (a+b)(f) := (\sum_{i=1}^{n_a} w_{a,i} ev_{x_{a,i}} + (a+b)(f) := + (\sum_{i=1}^{n_a} w_{a,i} ev_{x_{a,i}} + \sum_{i=1}^{n_b} w_{b,i} ev_{x_{b,j}})(f\omega), \forall f\in V_{d, \omega}. - Alternatively, if :math:`a, b` have different dimensions, define + More generally, if :math:`a, b` have different dimensions, define + :math:`d = max(dim(a), dim(b))`, and define .. math:: - a + b := argmax_{quad\in\{a, b\}}dim(quad). - - The rationale of this definition is that, the rule with lower - dimension is tensored with empty rules to "match" the higher - dimension before added together. - - Obviously, in either case the addition is commutative - :math:`a + b = b + a`. + a + b := em(a) + em(b) = argmax_{quad\in\{a, b\}}dim(quad), Note that addition of two rules with same dimension but different weight functions is not allowed. + 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, @@ -133,11 +217,53 @@ Multiplication .. 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,i}})(f\omega), + 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`. Obviously, the - multiplication of two quadrature rules is not commutative. + 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. + +For example, dentoe the (n-1)-th degree second kind Gauss-Chebyshev quadrature +rule as :math:`quad^{GC2}_n`, :math:`n\geq 1`, and let :math:`quad^{GC2}_0 = 0`. +Then :math:`Q := \{quad^{GC2}_n\}_{n=0}^{\infty}` is a consistent quadrature chain. +Let + +.. math:: + + R^{GC2} := \left\{ \sum_i \prod_j quad_{i,j} \right\}, \text{ } q_{i,j} \in Q. """ -- GitLab From be0fbdf2478625ed9970b01133e854b26126a7a7 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Sep 2019 14:10:19 -0500 Subject: [PATCH 21/25] Add tensor product quadrature factory, update docs --- doc/symbolic.rst | 40 +++++++ modepy/symbolic/__init__.py | 226 +++++++++++++++++++++++++----------- test/test_symbolic.py | 25 +++- 3 files changed, 223 insertions(+), 68 deletions(-) diff --git a/doc/symbolic.rst b/doc/symbolic.rst index becf18a..2d5f58c 100644 --- a/doc/symbolic.rst +++ b/doc/symbolic.rst @@ -1,6 +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/symbolic/__init__.py b/modepy/symbolic/__init__.py index b6f3b81..d86ccba 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -25,6 +25,9 @@ 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 @@ -134,7 +137,7 @@ Embedding 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) := + :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). @@ -183,15 +186,29 @@ Algebra of quadrature rules ^^^^^^^^^^^^^^^^^^^^^^^^^^^ Addition - Given a quadrature chain in :math:`d` dimensions with two elements - :math:`a, b`, the addition is defined as + Given two quadrature rules :math:`a, b`, both + in :math:`d` dimensions, the addition is defined as .. math:: - (a+b)(f) := - (\sum_{i=1}^{n_a} w_{a,i} ev_{x_{a,i}} - + \sum_{i=1}^{n_b} w_{b,i} ev_{x_{b,j}})(f\omega), - \forall f\in V_{d, \omega}. + (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 @@ -200,9 +217,6 @@ Addition a + b := em(a) + em(b) = argmax_{quad\in\{a, b\}}dim(quad), - Note that addition of two rules with same dimension but different - weight functions is not allowed. - Properties of addition: - :math:`dim(a + b) \in \{-1, max( dim(a), dim(b) )\}` @@ -256,15 +270,35 @@ the operations defined above. The generators are classic families of quadrature rules for 1D interval, such as Gauss-Legendre quadrature and Gauss-Chebyshev quadrature. -For example, dentoe the (n-1)-th degree second kind Gauss-Chebyshev quadrature -rule as :math:`quad^{GC2}_n`, :math:`n\geq 1`, and let :math:`quad^{GC2}_0 = 0`. -Then :math:`Q := \{quad^{GC2}_n\}_{n=0}^{\infty}` is a consistent quadrature chain. -Let +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} := \left\{ \sum_i \prod_j quad_{i,j} \right\}, \text{ } q_{i,j} \in Q. + 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 @@ -287,6 +321,10 @@ class TemplateQuadratureFamily(): 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 @@ -310,18 +348,32 @@ class TemplateQuadratureFamily(): 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, + 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): + 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: + 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: @@ -335,6 +387,12 @@ class TemplateQuadratureFamily(): 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: @@ -359,6 +417,19 @@ class TemplateQuadratureFamily(): 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) @@ -400,39 +471,69 @@ class TemplateQuadratureFamily(): return np.arange(sn, sn + self.level_to_nqpts(n)) + def get_quadrature_construction_params(self): + return self.quad_class, self.quad_class_extra_kwargs # }}} End class TemplateQuadratureFamily # {{{ 1d quadrature families -# Naming convention (for the name attribute): -# always end with the class name for the concrete quadrature class. -# LegendreGaussQuadrature +from modepy.quadrature.jacobi_gauss import ( + LegendreGaussQuadrature, ChebyshevGaussQuadrature) +from modepy.quadrature.clenshaw_curtis import ( + ClenshawCurtisQuadrature, FejerQuadrature) + +#: Gauss-Legendre quadrature family LEGENDRE_GAUSS_FAMILY = TemplateQuadratureFamily( - name='LegendreGaussQuadrature', - is_nested=False, - effective_measure=2, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, + quad_class=LegendreGaussQuadrature, ) -# ChebyshevGaussQuadrature(kind=1) +#: First kind Chebyshev-Gauss quadrature family CHEBYSHEV_GAUSS_1_FAMILY = TemplateQuadratureFamily( name='FirstKindChebyshevGaussQuadrature', - is_nested=False, + weight_func=lambda x: (1 - x**2)**(-1 / 2), effective_measure=np.pi, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, + quad_class=ChebyshevGaussQuadrature, + quad_class_extra_kwargs={'kind': 1, 'backend': 'scipy'}, ) -# ChebyshevGaussQuadrature(kind=2) +#: Second kind Chebyshev-Gauss quadrature family CHEBYSHEV_GAUSS_2_FAMILY = TemplateQuadratureFamily( name='SecondKindChebyshevGaussQuadrature', - is_nested=False, + weight_func=lambda x: (1 - x**2)**(1 / 2), effective_measure=np.pi / 2, level_to_nqpts=lambda n: n + 1, level_to_exact_to=lambda n: 2 * n + 1, + quad_class=ChebyshevGaussQuadrature, + quad_class_extra_kwargs={'kind': 2, 'backend': 'builtin'}, + ) + +#: Clenshaw-Curtis quadrature family +CLENSHAW_CURTIS_FAMILY = TemplateQuadratureFamily( + level_to_nqpts=lambda n: n + 1, + level_to_exact_to=lambda n: n, + quad_class=ClenshawCurtisQuadrature, + ) + +#: First kind Fejer quadrature family +FEJER_1_FAMILY = TemplateQuadratureFamily( + name='FirstKindFejerQuadrature', + level_to_nqpts=lambda n: n, + quad_class=FejerQuadrature, + quad_class_extra_kwargs={'kind': 1}, + ) + +#: Second kind Fejer quadrature family +FEJER_2_FAMILY = TemplateQuadratureFamily( + name='SecondKindFejerQuadrature', + level_to_nqpts=lambda n: n - 1, + quad_class=FejerQuadrature, + quad_class_extra_kwargs={'kind': 2}, ) # }}} End 1d quadrature families @@ -458,45 +559,18 @@ class QuadratureFactory(): QuadratureRingElement object that describes the quadrature rule algebraically. """ - supported_families = [ - LEGENDRE_GAUSS_FAMILY, - CHEBYSHEV_GAUSS_1_FAMILY, CHEBYSHEV_GAUSS_2_FAMILY, - ] - def __init__(self, quad_desc): assert isinstance(quad_desc, QuadratureRingElement) self.quad_desc = quad_desc def get_quadrature_construction_params(self): - dim = self.quad_desc.dim - quad_classes = [None, ] * dim - quad_extra_kwargs = [dict(), ] * dim - - for iaxis, qfam in enumerate(self.quad_desc.family): - - if qfam not in self.__class__.supported_families: - raise ValueError("Unsupported quadrature family %s" % qfam) - - if qfam == LEGENDRE_GAUSS_FAMILY: - from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature - quad_classes[iaxis] = LegendreGaussQuadrature - - elif qfam == CHEBYSHEV_GAUSS_1_FAMILY: - from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature - quad_classes[iaxis] = ChebyshevGaussQuadrature - quad_extra_kwargs[iaxis]['kind'] = 1 - quad_extra_kwargs[iaxis]['backend'] = 'scipy' - - elif qfam == CHEBYSHEV_GAUSS_2_FAMILY: - from modepy.quadrature.jacobi_gauss import ChebyshevGaussQuadrature - quad_classes[iaxis] = ChebyshevGaussQuadrature - quad_extra_kwargs[iaxis]['kind'] = 2 - quad_extra_kwargs[iaxis]['backend'] = 'builtin' - - else: - raise RuntimeError("This line should not be reached") - - return quad_classes, quad_extra_kwargs + classes = [] + params = [] + for quadrule in self.quad_desc.family: + clss, prms = quadrule.get_quadrature_construction_params() + classes.append(clss) + params.append(prms) + return classes, params def get_node_level_local_indices(self): """Get levels and level-local indices in one axis dimension. @@ -738,12 +812,17 @@ class QuadratureRingElement(): raise ValueError( "Cannot add objects of QuadratureRingElement and %s" % type(other)) - if self.dim != other.dim: - raise ValueError("Cannot add with different dims (%d, %d)" - % (self.dim, other.dim)) + + if self.dim > other.dim: + return self + elif self.dim < other.dim: + return other + if self.family != other.family: - raise ValueError("Cannot add with different quadrature families (%s, %s)" + 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) @@ -797,7 +876,7 @@ class QuadratureRingElement(): returns the tensor product of the two rules (self x other). """ if np.isscalar(other): - return self.__class__(self.dim, self.family, self.max_level, + 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,) @@ -843,6 +922,19 @@ class QuadratureRingElement(): % 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 diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 7c03aba..bbe3a03 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -27,6 +27,7 @@ from six.moves import range, zip # noqa import numpy as np import pytest import modepy.symbolic as sym +from modepy.symbolic.tensor_product import TensorProductQuadratureFactory @pytest.mark.parametrize('fam', [ @@ -71,7 +72,6 @@ def test_ring_generators(fam, lev): (sym.CHEBYSHEV_GAUSS_2_FAMILY, False), ]) def test_family_templates(fam, expect_nested): - assert fam in sym.QuadratureFactory.supported_families assert fam.name[-10:] == 'Quadrature' assert fam.is_nested == expect_nested @@ -154,6 +154,29 @@ def test_ring_multiplication_1(): # }}} 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 = 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()' -- GitLab From 3e7caf6b7eabd4a897663aee9189f961132acda8 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Sep 2019 14:16:28 -0500 Subject: [PATCH 22/25] Fix imports --- modepy/symbolic/__init__.py | 2 +- test/test_symbolic.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/modepy/symbolic/__init__.py b/modepy/symbolic/__init__.py index d86ccba..34670f7 100644 --- a/modepy/symbolic/__init__.py +++ b/modepy/symbolic/__init__.py @@ -1,4 +1,4 @@ -from __future__ import division +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2019 Xiaoyu Wei" diff --git a/test/test_symbolic.py b/test/test_symbolic.py index bbe3a03..2f5c8e7 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -27,7 +27,7 @@ from six.moves import range, zip # noqa import numpy as np import pytest import modepy.symbolic as sym -from modepy.symbolic.tensor_product import TensorProductQuadratureFactory +import modepy.symbolic.tensor_product as tensor_product @pytest.mark.parametrize('fam', [ @@ -159,7 +159,7 @@ def test_ring_multiplication_1(): def drive_test_tensor_prod_quad(dim, level): family = sym.LEGENDRE_GAUSS_FAMILY generator = sym.QuadratureRingGenerator(family, level) - quad_factory = TensorProductQuadratureFactory(generator, dim) + quad_factory = tensor_product.TensorProductQuadratureFactory(generator, dim) oned_rule = generator.get_quadrature_rule() quad_rule = quad_factory() -- GitLab From b0a4d7b72054a842bd6af5266426e30c16d9ff06 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Sep 2019 14:17:36 -0500 Subject: [PATCH 23/25] Add new py file --- modepy/symbolic/tensor_product.py | 49 +++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 modepy/symbolic/tensor_product.py diff --git a/modepy/symbolic/tensor_product.py b/modepy/symbolic/tensor_product.py new file mode 100644 index 0000000..404a032 --- /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() -- GitLab From a54366ca381d5c7918709e141ec883fa58e8392c Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 26 Sep 2019 10:33:10 -0500 Subject: [PATCH 24/25] Fix doc generation --- modepy/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index e276d94..02137c2 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -27,9 +27,10 @@ THE SOFTWARE. from modepy.modes import ( jacobi, grad_jacobi, simplex_onb, grad_simplex_onb, simplex_monomial_basis, grad_simplex_monomial_basis, - tensor_product_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, -- GitLab From fed4319ffdecb975ee5421e8a6163e96bd959d5c Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 26 Sep 2019 11:09:11 -0500 Subject: [PATCH 25/25] flake8 fixes --- modepy/__init__.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/modepy/__init__.py b/modepy/__init__.py index 02137c2..b29364c 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -30,7 +30,8 @@ from modepy.modes import ( 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, tensor_product_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, @@ -56,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", -- GitLab