From 27d44fb6bddb35d680de4232893c93dc3b74ed8f Mon Sep 17 00:00:00 2001 From: Jonathan Feinberg Date: Tue, 16 Mar 2021 20:22:55 +0100 Subject: [PATCH 1/4] refactor quadrature function layout --- .../sampler/sequences/chebyshev.py | 4 +- .../distributions/sampler/sequences/grid.py | 4 +- chaospy/quadrature/__init__.py | 89 +++- chaospy/quadrature/chebyshev.py | 108 +++++ chaospy/quadrature/clenshaw_curtis.py | 205 ++------- chaospy/quadrature/discrete.py | 107 ++--- chaospy/quadrature/fejer.py | 159 ------- chaospy/quadrature/fejer_1.py | 91 ++++ chaospy/quadrature/fejer_2.py | 76 ++++ chaospy/quadrature/frontend.py | 115 ++--- chaospy/quadrature/gauss_legendre.py | 153 ------- chaospy/quadrature/gaussian.py | 17 +- chaospy/quadrature/gegenbauer.py | 59 +++ chaospy/quadrature/genz_keister/__init__.py | 11 - chaospy/quadrature/genz_keister/frontend.py | 56 --- chaospy/quadrature/genz_keister/gk16.py | 374 ---------------- chaospy/quadrature/genz_keister/gk18.py | 190 -------- chaospy/quadrature/genz_keister/gk22.py | 198 --------- chaospy/quadrature/genz_keister/gk24.py | 202 --------- chaospy/quadrature/grid.py | 44 +- chaospy/quadrature/hermite.py | 59 +++ chaospy/quadrature/hypercube.py | 411 ++++++++++++++++++ chaospy/quadrature/jacobi.py | 61 +++ .../{gauss_kronrod.py => kronrod.py} | 27 +- chaospy/quadrature/laguerre.py | 54 +++ chaospy/quadrature/legendre.py | 163 +++++++ chaospy/quadrature/leja.py | 12 +- .../{gauss_lobatto.py => lobatto.py} | 35 +- chaospy/quadrature/newton_cotes.py | 89 +--- .../{gauss_patterson.py => patterson.py} | 76 ++-- .../quadrature/{gauss_radau.py => radau.py} | 56 +-- chaospy/quadrature/sparse_grid.py | 9 +- chaospy/quadrature/{combine.py => utils.py} | 8 +- docs/reference/quadrature.rst | 64 ++- docs/user_guide/quadrature.rst | 383 ++++++++++------ pyproject.toml | 2 +- tests/recurrence/test_quadrature_creation.py | 4 +- 37 files changed, 1764 insertions(+), 2011 deletions(-) create mode 100644 chaospy/quadrature/chebyshev.py delete mode 100644 chaospy/quadrature/fejer.py create mode 100644 chaospy/quadrature/fejer_1.py create mode 100644 chaospy/quadrature/fejer_2.py delete mode 100644 chaospy/quadrature/gauss_legendre.py create mode 100644 chaospy/quadrature/gegenbauer.py delete mode 100644 chaospy/quadrature/genz_keister/__init__.py delete mode 100644 chaospy/quadrature/genz_keister/frontend.py delete mode 100644 chaospy/quadrature/genz_keister/gk16.py delete mode 100644 chaospy/quadrature/genz_keister/gk18.py delete mode 100644 chaospy/quadrature/genz_keister/gk22.py delete mode 100644 chaospy/quadrature/genz_keister/gk24.py create mode 100644 chaospy/quadrature/hermite.py create mode 100644 chaospy/quadrature/hypercube.py create mode 100644 chaospy/quadrature/jacobi.py rename chaospy/quadrature/{gauss_kronrod.py => kronrod.py} (92%) create mode 100644 chaospy/quadrature/laguerre.py create mode 100644 chaospy/quadrature/legendre.py rename chaospy/quadrature/{gauss_lobatto.py => lobatto.py} (84%) rename chaospy/quadrature/{gauss_patterson.py => patterson.py} (97%) rename chaospy/quadrature/{gauss_radau.py => radau.py} (78%) rename chaospy/quadrature/{combine.py => utils.py} (92%) diff --git a/chaospy/distributions/sampler/sequences/chebyshev.py b/chaospy/distributions/sampler/sequences/chebyshev.py index e373eddf..00c46a8c 100644 --- a/chaospy/distributions/sampler/sequences/chebyshev.py +++ b/chaospy/distributions/sampler/sequences/chebyshev.py @@ -42,7 +42,9 @@ """ import numpy + import chaospy +from chaospy.quadrature import utils def create_chebyshev_samples(order, dim=1): @@ -60,7 +62,7 @@ def create_chebyshev_samples(order, dim=1): ``[0, 1]^dim`` hyper-cube and ``shape == (dim, order)``. """ x_data = .5*numpy.cos(numpy.arange(order, 0, -1)*numpy.pi/(order+1)) + .5 - x_data = chaospy.quadrature.combine([x_data]*dim) + x_data = utils.combine([x_data]*dim) return x_data.T diff --git a/chaospy/distributions/sampler/sequences/grid.py b/chaospy/distributions/sampler/sequences/grid.py index 3db6e4c0..12f44f1f 100644 --- a/chaospy/distributions/sampler/sequences/grid.py +++ b/chaospy/distributions/sampler/sequences/grid.py @@ -42,7 +42,9 @@ """ import numpy + import chaospy +from chaospy.quadrature import utils def create_grid_samples(order, dim=1): @@ -59,7 +61,7 @@ def create_grid_samples(order, dim=1): Regular grid with ``shape == (dim, order)``. """ x_data = numpy.arange(1, order+1)/(order+1.) - x_data = chaospy.quadrature.combine([x_data]*dim) + x_data = utils.combine([x_data]*dim) return x_data.T diff --git a/chaospy/quadrature/__init__.py b/chaospy/quadrature/__init__.py index d166d691..afc93b32 100644 --- a/chaospy/quadrature/__init__.py +++ b/chaospy/quadrature/__init__.py @@ -1,18 +1,75 @@ r"""Collection of quadrature methods.""" +import logging +from functools import wraps + from .frontend import generate_quadrature -from .sparse_grid import construct_sparse_grid -from .combine import combine - -from .clenshaw_curtis import quad_clenshaw_curtis -from .discrete import quad_discrete -from .fejer import quad_fejer -from .gaussian import quad_gaussian -from .gauss_patterson import quad_gauss_patterson -from .gauss_legendre import quad_gauss_legendre -from .gauss_lobatto import quad_gauss_lobatto -from .gauss_kronrod import quad_gauss_kronrod, kronrod_jacobi -from .gauss_radau import quad_gauss_radau -from .genz_keister import quad_genz_keister -from .grid import quad_grid -from .leja import quad_leja -from .newton_cotes import quad_newton_cotes +from .sparse_grid import sparse_grid +from .utils import combine + +from .chebyshev import chebyshev_1, chebyshev_2 +from .clenshaw_curtis import clenshaw_curtis +from .discrete import discrete +from .fejer_1 import fejer_1 +from .fejer_2 import fejer_2 +from .gaussian import gaussian +from .gegenbauer import gegenbauer +from .grid import grid +from .hermite import hermite +from .jacobi import jacobi +from .kronrod import kronrod, kronrod_jacobi +from .laguerre import laguerre +from .legendre import legendre, legendre_proxy +from .leja import leja +from .lobatto import lobatto +from .newton_cotes import newton_cotes +from .patterson import patterson +from .radau import radau + + +INTEGRATION_COLLECTION = { + "clenshaw_curtis": clenshaw_curtis, + "discrete": discrete, + "fejer_1": fejer_1, + "fejer_2": fejer_2, + "gaussian": gaussian, + "grid": grid, + "kronrod": kronrod, + "legendre": legendre_proxy, + "leja": leja, + "lobatto": lobatto, + "newton_cotes": newton_cotes, + "patterson": patterson, + "radau": radau, +} + + +def quadrature_deprecation_warning(name, func=None): + """Announce deprecation warning for quad-func.""" + + if func is None: + func = globals()[name] + quad_name = "quad_%s" % name + + @wraps(func) + def wrapped(*args, **kwargs): + """Function wrapper adds warnings.""" + logger = logging.getLogger(__name__) + logger.warning("chaospy.%s name is to be deprecated; " + "Use chaospy.quadrature.%s instead", + quad_name, func.__name__) + return func(*args, **kwargs) + + globals()[quad_name] = wrapped + +quadrature_deprecation_warning("clenshaw_curtis", clenshaw_curtis) +quadrature_deprecation_warning("discrete", discrete) +quadrature_deprecation_warning("fejer", fejer_2) +quadrature_deprecation_warning("grid", grid) +quadrature_deprecation_warning("gaussian", gaussian) +quadrature_deprecation_warning("newton_cotes", newton_cotes) +quadrature_deprecation_warning("leja", leja) +quadrature_deprecation_warning("gauss_legendre", legendre_proxy) +quadrature_deprecation_warning("gauss_kronrod", kronrod) +quadrature_deprecation_warning("gauss_lobatto", lobatto) +quadrature_deprecation_warning("gauss_patterson", patterson) +quadrature_deprecation_warning("gauss_radau", radau) diff --git a/chaospy/quadrature/chebyshev.py b/chaospy/quadrature/chebyshev.py new file mode 100644 index 00000000..df5db499 --- /dev/null +++ b/chaospy/quadrature/chebyshev.py @@ -0,0 +1,108 @@ +"""Chebyshev-Gauss quadrature rule of the first kind.""" +import numpy +import chaospy + +from .hypercube import hypercube_quadrature + + +def chebyshev_1(order, lower=-1, upper=1, physicist=False): + r""" + Chebyshev-Gauss quadrature rule of the first kind. + + Compute the sample points and weights for Chebyshev-Gauss quadrature. The + sample points are the roots of the nth degree Chebyshev polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + first order Chebyshev-Gauss physicist means a weight function + :math:`1/\sqrt{1-x^2}` and weights that sum to :math`1/2`, and probabilist + means a weight function is :math:`1/\sqrt{x (1-x)}` and sum to 1. + + Args: + order (int): + The quadrature order. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.chebyshev_1(3) + >>> abscissas + array([[-0.92387953, -0.38268343, 0.38268343, 0.92387953]]) + >>> weights + array([0.25, 0.25, 0.25, 0.25]) + + See also: + :func:`chaospy.quadrature.chebyshev_2` + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Beta(0.5, 0.5, lower, upper)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights *= 0.5 if physicist else 1 + return abscissas[numpy.newaxis], weights + + + +def chebyshev_2(order, lower=-1, upper=1, physicist=False): + r""" + Chebyshev-Gauss quadrature rule of the second kind. + + Compute the sample points and weights for Chebyshev-Gauss quadrature. The + sample points are the roots of the nth degree Chebyshev polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + second order Chebyshev-Gauss physicist means a weight function + :math:`\sqrt{1-x^2}` and weights that sum to :math`2`, and probabilist + means a weight function is :math:`\sqrt{x (1-x)}` and sum to 1. + + Args: + order (int): + The quadrature order. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.chebyshev_2(3) + >>> abscissas + array([[-0.80901699, -0.30901699, 0.30901699, 0.80901699]]) + >>> weights + array([0.1381966, 0.3618034, 0.3618034, 0.1381966]) + + See also: + :func:`chaospy.quadrature.chebyshev_1` + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Beta(1.5, 1.5, lower, upper)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights *= 2 if physicist else 1 + return abscissas[numpy.newaxis], weights diff --git a/chaospy/quadrature/clenshaw_curtis.py b/chaospy/quadrature/clenshaw_curtis.py index 0ec804df..1d22bd6e 100644 --- a/chaospy/quadrature/clenshaw_curtis.py +++ b/chaospy/quadrature/clenshaw_curtis.py @@ -1,60 +1,16 @@ -""" -Generate the quadrature nodes and weights in Clenshaw-Curtis quadrature. - - -Example usage -------------- - -The first few orders with linear growth rule:: - - >>> distribution = chaospy.Uniform(0, 1) - >>> for order in [0, 1, 2, 3]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="clenshaw_curtis") - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0.5]] [1.] - 1 [[0. 1.]] [0.5 0.5] - 2 [[0. 0.5 1. ]] [0.167 0.667 0.167] - 3 [[0. 0.25 0.75 1. ]] [0.056 0.444 0.444 0.056] - -The first few orders with exponential growth rule:: - - >>> for order in [0, 1, 2]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="clenshaw_curtis", growth=True) - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0.5]] [1.] - 1 [[0. 0.5 1. ]] [0.167 0.667 0.167] - 2 [[0. 0.146 0.5 0.854 1. ]] [0.033 0.267 0.4 0.267 0.033] - -Applying the rule using Smolyak sparse grid:: - - >>> distribution = chaospy.Iid(chaospy.Uniform(0, 1), 2) - >>> abscissas, weights = chaospy.generate_quadrature( - ... 2, distribution, rule="clenshaw_curtis", - ... growth=True, sparse=True) - >>> abscissas.round(2) - array([[0. , 0. , 0. , 0.15, 0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.85, 1. , - 1. , 1. ], - [0. , 0.5 , 1. , 0.5 , 0. , 0.15, 0.5 , 0.85, 1. , 0.5 , 0. , - 0.5 , 1. ]]) - >>> weights.round(3) - array([ 0.028, -0.022, 0.028, 0.267, -0.022, 0.267, -0.089, 0.267, - -0.022, 0.267, 0.028, -0.022, 0.028]) -""" -from __future__ import division +"""Generate the quadrature nodes and weights in Clenshaw-Curtis quadrature.""" try: from functools import lru_cache -except ImportError: # pragma: no covere +except ImportError: # pragma: no coverage from functools32 import lru_cache import numpy import chaospy -from .combine import combine_quadrature +from .hypercube import hypercube_quadrature -def quad_clenshaw_curtis(order, domain, growth=False, segments=1): +def clenshaw_curtis(order, domain=(0., 1.), growth=False, segments=1): """ Generate the quadrature nodes and weights in Clenshaw-Curtis quadrature. @@ -71,133 +27,59 @@ def quad_clenshaw_curtis(order, domain, growth=False, segments=1): Args: order (int, numpy.ndarray): Quadrature order. - domain (chaospy.Distribution, numpy.ndarray): + domain (:class:`chaospy.Distribution`, numpy.ndarray): Either distribution or bounding of interval to integrate over. growth (bool): If True sets the growth rule for the quadrature rule to only include orders that enhances nested samples. segments (int): - Split intervals into N subintervals and create a patched + Split intervals into steps subintervals and create a patched quadrature based on the segmented quadrature. Can not be lower than `order`. If 0 is provided, default to square root of `order`. - Nested samples only exist when the number of segments are fixed. + Nested samples only appear when the number of segments are fixed. Returns: abscissas (numpy.ndarray): The quadrature points for where to evaluate the model function - with ``abscissas.shape == (len(dist), N)`` where ``N`` is the - number of samples. + with ``abscissas.shape == (len(dist), steps)`` where ``steps`` is + the number of samples. weights (numpy.ndarray): - The quadrature weights with ``weights.shape == (N,)``. + The quadrature weights with ``weights.shape == (steps,)``. Notes: Implemented as proposed by Waldvogel :cite:`waldvogel_fast_2006`. Example: - >>> abscissas, weights = quad_clenshaw_curtis(4, (0, 1)) + >>> abscissas, weights = chaospy.quadrature.clenshaw_curtis(4, (0, 1)) >>> abscissas.round(4) array([[0. , 0.1464, 0.5 , 0.8536, 1. ]]) >>> weights.round(4) array([0.0333, 0.2667, 0.4 , 0.2667, 0.0333]) - >>> abscissas, weights = quad_clenshaw_curtis(4, (0, 1), segments=0) - >>> abscissas.round(4) - array([[0. , 0.25, 0.5 , 0.75, 1. ]]) - >>> weights.round(4) - array([0.0833, 0.3333, 0.1667, 0.3333, 0.0833]) - - """ - if isinstance(domain, chaospy.Distribution): - abscissas, weights = quad_clenshaw_curtis( - order, (domain.lower, domain.upper), growth, segments) - - # Sometimes edge samples (inside the domain) falls out again from simple - # rounding errors. Edge samples needs to be adjusted. - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - order = numpy.asarray(order, dtype=int).flatten() - lower, upper = numpy.array(domain) - lower = numpy.asarray(lower).flatten() - upper = numpy.asarray(upper).flatten() - dim = max(lower.size, upper.size, order.size) + See also: + :func:`chaospy.quadrature.gaussian` + :func:`chaospy.quadrature.fejer_1` + :func:`chaospy.quadrature.fejer_2` - order = order*numpy.ones(dim, dtype=int) - lower = lower*numpy.ones(dim) - upper = upper*numpy.ones(dim) - segments = segments*numpy.ones(dim, dtype=int) - - if growth: - order = numpy.where(order > 0, 2**order, 0) - - abscissas, weights = zip(*[_clenshaw_curtis(order_, segment) - for order_, segment in zip(order, segments)]) - - return combine_quadrature(abscissas, weights, (lower, upper)) + """ + order = numpy.asarray(order) + order = numpy.where(growth, numpy.where(order > 0, 2**order, 0), order) + return hypercube_quadrature( + quad_func=clenshaw_curtis_simple, + order=order, + domain=domain, + segments=segments, + ) @lru_cache(None) -def _clenshaw_curtis(order, segments=1): - r""" - Backend method. - - Examples: - >>> abscissas, weights = _clenshaw_curtis(0, 0) - >>> abscissas - array([0.5]) - >>> weights - array([1.]) - >>> abscissas, weights = _clenshaw_curtis(2, 0) - >>> abscissas - array([0. , 0.5, 1. ]) - >>> weights - array([0.16666667, 0.66666667, 0.16666667]) - >>> abscissas, weights = _clenshaw_curtis(4, 0) - >>> abscissas - array([0. , 0.25, 0.5 , 0.75, 1. ]) - >>> weights - array([0.08333333, 0.33333333, 0.16666667, 0.33333333, 0.08333333]) - >>> abscissas, weights = _clenshaw_curtis(8, 0) - >>> abscissas.round(3) - array([0. , 0.073, 0.25 , 0.427, 0.5 , 0.573, 0.75 , 0.927, 1. ]) - >>> weights.round(3) - array([0.017, 0.133, 0.2 , 0.133, 0.033, 0.133, 0.2 , 0.133, 0.017]) - >>> abscissas, weights = _clenshaw_curtis(16, 0) - >>> abscissas.round(3) - array([0. , 0.037, 0.125, 0.213, 0.25 , 0.287, 0.375, 0.463, 0.5 , - 0.537, 0.625, 0.713, 0.75 , 0.787, 0.875, 0.963, 1. ]) - >>> weights.round(3) - array([0.008, 0.067, 0.1 , 0.067, 0.017, 0.067, 0.1 , 0.067, 0.017, - 0.067, 0.1 , 0.067, 0.017, 0.067, 0.1 , 0.067, 0.008]) +def clenshaw_curtis_simple(order): """ - if segments != 1 and order > 2: - if not segments: - segments = int(numpy.sqrt(order)) - assert segments < order, "few samples to distribute than intervals" - abscissas = [] - weights = [] - - nodes = numpy.linspace(0, 1, segments+1) - for idx, (lower, upper) in enumerate(zip(nodes[:-1], nodes[1:])): - - order_ = order//segments + (idx < (order%segments)) - abscissa, weight = _clenshaw_curtis(order_, segments=1) - abscissa = abscissa*(upper-lower) + lower - weight = weight*(upper-lower) - if abscissas: - weights[-1] += weight[0] - abscissa = abscissa[1:] - weight = weight[1:] - abscissas.extend(abscissa) - weights.extend(weight) - - assert len(abscissas) == order+1, (len(abscissas), order+1) - assert len(weights) == order+1 - return numpy.array(abscissas), numpy.array(weights) + Backend for Clenshaw-Curtis quadrature. + Use :func:`chaospy.quadrature.clenshaw_curtis` instead. + """ + order = int(order) if order == 0: return numpy.array([.5]), numpy.array([1.]) elif order == 1: @@ -206,21 +88,22 @@ def _clenshaw_curtis(order, segments=1): theta = (order-numpy.arange(order+1))*numpy.pi/order abscissas = 0.5*numpy.cos(theta)+0.5 - N = numpy.arange(1, order, 2) - length = len(N) - m = order-length + steps = numpy.arange(1, order, 2) + length = len(steps) + remains = order-length - v0 = numpy.concatenate([2./(N*(N-2)), [1./N[-1]], numpy.zeros(m)]) - v2 = -v0[:-1]-v0[:0:-1] - g0 = -numpy.ones(order) - g0[length] += order - g0[m] += order - g = g0/(order**2-1+(order%2)) + beta = numpy.hstack([2./(steps*(steps-2)), [1./steps[-1]], numpy.zeros(remains)]) + beta = -beta[:-1]-beta[:0:-1] - w = numpy.fft.ihfft(v2+g) - assert max(w.imag) < 1e-15 - w = w.real + gamma = -numpy.ones(order) + gamma[length] += order + gamma[remains] += order + gamma /= (order**2-1+(order%2)) - weights = numpy.concatenate([w, w[len(w)-2+(order%2)::-1]]) + weights = numpy.fft.ihfft(beta+gamma) + assert max(weights.imag) < 1e-15 + weights = weights.real + weights = numpy.hstack([weights, weights[len(weights)-2+(order%2)::-1]])/2 + assert numpy.isclose(numpy.sum(weights), 1) - return abscissas, weights/2 + return abscissas, weights diff --git a/chaospy/quadrature/discrete.py b/chaospy/quadrature/discrete.py index 65093fd3..3cc7f42e 100644 --- a/chaospy/quadrature/discrete.py +++ b/chaospy/quadrature/discrete.py @@ -1,62 +1,23 @@ -""" -Generate the quadrature abscissas and weights for simple grid. - -Available to ensure that discrete distributions works along side -continuous ones. - -Example usage -------------- - -The first few orders with linear growth rule:: - - >>> distribution = chaospy.DiscreteUniform(-2, 2) - >>> for order in [0, 1, 2, 3, 4, 5, 9]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="discrete") - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0]] [1.] - 1 [[-1 1]] [0.5 0.5] - 2 [[-2 0 2]] [0.333 0.333 0.333] - 3 [[-2 -1 1 2]] [0.25 0.25 0.25 0.25] - 4 [[-2 -1 0 1 2]] [0.2 0.2 0.2 0.2 0.2] - 5 [[-2 -1 0 1 2]] [0.2 0.2 0.2 0.2 0.2] - 9 [[-2 -1 0 1 2]] [0.2 0.2 0.2 0.2 0.2] - -As the accuracy of discrete distribution plateau when all contained values are -included, there is no reason to increase the number of nodes after this point. - -The first few orders with exponential growth rule where the nodes are nested:: - - >>> distribution = chaospy.DiscreteUniform(0, 10) - >>> for order in [0, 1, 2, 3, 4]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="discrete", growth=True) - ... print(order, abscissas) - 0 [[5]] - 1 [[1 5 9]] - 2 [[1 3 5 7 9]] - 3 [[ 0 1 3 4 5 6 7 9 10]] - 4 [[ 0 1 2 3 4 5 6 7 8 9 10]] -""" +"""Generate quadrature abscissas and weights for discrete distributions.""" import numpy import chaospy -from .combine import combine_quadrature -from .grid import quad_grid +from .hypercube import hypercube_quadrature -def quad_discrete(order, domain=(0, 1), growth=False, segments=1): +def discrete(order, domain=(0, 1), growth=False): """ Generate quadrature abscissas and weights for discrete distributions. - Same as regular grid, but `order` plateau at the `upper-lower-1`. - At this order, finite state discrete distributions are analytically - correct, and higher order will make the accuracy worsen. + A specialized quadrature designed for discrete distributions. It is defined + as an evenly spaced grid on the domain, rounded to nearest integer. Rule + will converge to where all integer values on the domain is covered. This + ensure that only necessary samples are evaluated. Args: order (int, numpy.ndarray): Quadrature order. - domain (chaospy.distributions.baseclass.Distribution, numpy.ndarray): + domain (:class:`chaospy.Distribution`, numpy.ndarray): Either distribution or bounding of interval to integrate over. growth (bool): if true sets the growth rule for the quadrature rule to only @@ -70,36 +31,36 @@ def quad_discrete(order, domain=(0, 1), growth=False, segments=1): Either distribution or bounding of interval to integrate over. Examples: - >>> distribution = chaospy.DiscreteUniform(-2, 2) - >>> abscissas, weights = chaospy.quad_discrete(4, distribution) - >>> abscissas.round(4) - array([[-2., -1., 0., 1., 2.]]) + >>> distribution = chaospy.Binomial(6, 0.4) + >>> abscissas, weights = chaospy.quadrature.discrete(4, distribution) + >>> abscissas + array([[0., 2., 3., 4., 6.]]) >>> weights.round(4) - array([0.2, 0.2, 0.2, 0.2, 0.2]) - >>> abscissas, weights = chaospy.quad_discrete(9, distribution) - >>> abscissas.round(4) - array([[-2., -1., 0., 1., 2.]]) + array([0.0601, 0.4006, 0.3561, 0.178 , 0.0053]) + >>> abscissas, weights = chaospy.quadrature.discrete(10, distribution) + >>> abscissas + array([[0., 1., 2., 3., 4., 5., 6.]]) >>> weights.round(4) - array([0.2, 0.2, 0.2, 0.2, 0.2]) + array([0.0467, 0.1866, 0.311 , 0.2765, 0.1382, 0.0369, 0.0041]) """ - if isinstance(domain, chaospy.Distribution): - abscissas, weights = quad_discrete(order, (domain.lower, domain.upper), growth=growth, segments=segments) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - del segments # Basically segments doesn't do much here - if growth: - order = numpy.where(order > 0, 2**(order), 0) + order = numpy.asarray(order) + order = numpy.where(growth, numpy.where(order > 0, 2**order, 0), order) + return hypercube_quadrature( + quad_func=discrete_simple, + order=order, + domain=domain, + auto_scale=False, + ) - order = numpy.atleast_1d(order) - order, lower, upper = numpy.broadcast_arrays(order, domain[0], domain[1]) - assert order.ndim == 1, "too many dimensions" - order_max = numpy.round(upper-lower).astype(int)-1 - order = numpy.where(order > order_max, order_max, order) +def discrete_simple(order, lower=-2, upper=2): + """ + Backend for discrete quadrature. - return quad_grid(order, (lower, upper)) + Use :func:`chaospy.quadrature.discrete` instead. + """ + order = min(order, round(upper-lower)-1) + abscissas = numpy.linspace(lower, upper, 2*order+3)[1::2].round() + weights = numpy.full(order+1, (upper-lower)/(order+1.)) + return abscissas, weights diff --git a/chaospy/quadrature/fejer.py b/chaospy/quadrature/fejer.py deleted file mode 100644 index ec76c426..00000000 --- a/chaospy/quadrature/fejer.py +++ /dev/null @@ -1,159 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Generate the quadrature abscissas and weights in Fejer quadrature. - -Example usage -------------- - -The first few orders with linear growth rule:: - - >>> distribution = chaospy.Uniform(0, 1) - >>> for order in [0, 1, 2, 3]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="fejer") - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0.5]] [1.] - 1 [[0.25 0.75]] [0.5 0.5] - 2 [[0.146 0.5 0.854]] [0.286 0.429 0.286] - 3 [[0.095 0.345 0.655 0.905]] [0.188 0.312 0.312 0.188] - -The first few orders with exponential growth rule:: - - >>> for order in [0, 1, 2]: # doctest: +NORMALIZE_WHITESPACE - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="fejer", growth=True) - ... print(order, abscissas.round(2), weights.round(2)) - 0 [[0.5]] [1.] - 1 [[0.15 0.5 0.85]] [0.29 0.43 0.29] - 2 [[0.04 0.15 0.31 0.5 0.69 0.85 0.96]] - [0.07 0.14 0.18 0.2 0.18 0.14 0.07] - -Applying the rule using Smolyak sparse grid:: - - >>> distribution = chaospy.Iid(chaospy.Uniform(0, 1), 2) - >>> abscissas, weights = chaospy.generate_quadrature( - ... 2, distribution, rule="fejer", growth=True, sparse=True) - >>> abscissas.round(3) - array([[0.038, 0.146, 0.146, 0.146, 0.309, 0.5 , 0.5 , 0.5 , 0.5 , - 0.5 , 0.5 , 0.5 , 0.691, 0.854, 0.854, 0.854, 0.962], - [0.5 , 0.146, 0.5 , 0.854, 0.5 , 0.038, 0.146, 0.309, 0.5 , - 0.691, 0.854, 0.962, 0.5 , 0.146, 0.5 , 0.854, 0.5 ]]) - >>> weights.round(3) - array([ 0.074, 0.082, -0.021, 0.082, 0.184, 0.074, -0.021, 0.184, - -0.273, 0.184, -0.021, 0.074, 0.184, 0.082, -0.021, 0.082, - 0.074]) -""" -from __future__ import division -try: - from functools import lru_cache -except ImportError: # pragma: no cover - from functools32 import lru_cache - -import numpy -import chaospy - -from .combine import combine_quadrature -from .clenshaw_curtis import _clenshaw_curtis - - -def quad_fejer(order, domain=(0, 1), growth=False, segments=1): - """ - Generate the quadrature abscissas and weights in Fejér quadrature. - - Fejér proposed two quadrature rules very similar to - :func:`quad_clenshaw_curtis`. The only difference is that the endpoints are - removed. That is, Fejér only used the interior extrema of the Chebyshev - polynomials, i.e. the true stationary points. This makes this a better - method for performing quadrature on infinite intervals, as the evaluation - does not contain illegal values. - - Args: - order (int, numpy.ndarray): - Quadrature order. - domain (chaospy.Distribution, numpy.ndarray): - Either distribution or bounding of interval to integrate over. - growth (bool): - If True sets the growth rule for the quadrature rule to only - include orders that enhances nested samples. - segments (int): - Split intervals into N subintervals and create a patched - quadrature based on the segmented quadrature. Can not be lower than - `order`. If 0 is provided, default to square root of `order`. - Nested samples only exist when the number of segments are fixed. - - Returns: - abscissas (numpy.ndarray): - The quadrature points for where to evaluate the model function - with ``abscissas.shape == (len(dist), N)`` where ``N`` is the - number of samples. - weights (numpy.ndarray): - The quadrature weights with ``weights.shape == (N,)``. - - Notes: - Implemented as proposed by Waldvogel :cite:`waldvogel_fast_2006`. - - Example: - >>> abscissas, weights = quad_fejer(3, (0, 1)) - >>> abscissas.round(4) - array([[0.0955, 0.3455, 0.6545, 0.9045]]) - >>> weights.round(4) - array([0.1804, 0.2996, 0.2996, 0.1804]) - >>> abscissas, weights = quad_fejer(3, (0, 1), segments=2) - >>> abscissas.round(4) - array([[0.125, 0.375, 0.625, 0.875]]) - >>> weights.round(4) - array([0.2222, 0.2222, 0.2222, 0.2222]) - """ - if isinstance(domain, chaospy.Distribution): - abscissas, weights = quad_fejer( - order, (domain.lower, domain.upper), growth) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - order = numpy.asarray(order, dtype=int).flatten() - lower, upper = numpy.array(domain) - lower = numpy.asarray(lower).flatten() - upper = numpy.asarray(upper).flatten() - - dim = max(lower.size, upper.size, order.size) - - order = order*numpy.ones(dim, dtype=int) - lower = lower*numpy.ones(dim) - upper = upper*numpy.ones(dim) - segments = segments*numpy.ones(dim, dtype=int) - - if growth: - order = numpy.where(order > 0, 2**(order+1)-2, 0) - - abscissas, weights = zip(*[_fejer(order_, segment) - for order_, segment in zip(order, segments)]) - - return combine_quadrature(abscissas, weights, (lower, upper)) - - -@lru_cache(None) -def _fejer(order, segments=1): - """Backend method.""" - if segments != 1 and order > 2: - if not segments: - segments = int(numpy.sqrt(order)) - assert segments < order, "few samples to distribute than intervals" - abscissas = [] - weights = [] - - nodes = numpy.linspace(0, 1, segments+1) - for idx, (lower, upper) in enumerate(zip(nodes[:-1], nodes[1:])): - order_ = order//segments + (idx+1 < (order%segments)) - abscissa, weight = _fejer(order_, segments=1) - abscissas.extend(abscissa*(upper-lower) + lower) - weights.extend(weight*(upper-lower)) - - assert len(abscissas) == order+1, (len(abscissas), order+1) - assert len(weights) == order+1 - return numpy.array(abscissas), numpy.array(weights) - - abscissas, weights = _clenshaw_curtis(order+2, segments) - return abscissas[1:-1], weights[1:-1] diff --git a/chaospy/quadrature/fejer_1.py b/chaospy/quadrature/fejer_1.py new file mode 100644 index 00000000..b4e4e8b7 --- /dev/null +++ b/chaospy/quadrature/fejer_1.py @@ -0,0 +1,91 @@ +# -*- coding: utf-8 -*- +"""Generate the quadrature abscissas and weights in Fejér type I quadrature.""" +try: + from functools import lru_cache +except ImportError: # pragma: no cover + from functools32 import lru_cache +import numpy + +from .hypercube import hypercube_quadrature + + +def fejer_1(order, domain=(0, 1), growth=False, segments=1): + """ + Generate the quadrature abscissas and weights in Fejér type I quadrature. + + Fejér proposed two quadrature rules very similar to + :func:`chaospy.quadrature.clenshaw_curtis`, except that it does not + include endpoints, making it more suitable to integrate unbound interval. + + Args: + order (int, numpy.ndarray): + Quadrature order. + domain (chaospy.Distribution, numpy.ndarray): + Either distribution or bounding of interval to integrate over. + growth (bool): + If True sets the growth rule for the quadrature rule to only + include orders that enhances nested samples. + segments (int): + Split intervals into N subintervals and create a patched + quadrature based on the segmented quadrature. Can not be lower than + `order`. If 0 is provided, default to square root of `order`. + Nested samples only exist when the number of segments are fixed. + + Returns: + abscissas (numpy.ndarray): + The quadrature points for where to evaluate the model function + with ``abscissas.shape == (len(dist), N)`` where ``N`` is the + number of samples. + weights (numpy.ndarray): + The quadrature weights with ``weights.shape == (N,)``. + + Notes: + Implemented as proposed by Waldvogel :cite:`waldvogel_fast_2006`. + + Example: + >>> abscissas, weights = chaospy.quadrature.fejer_1(3, (0, 1)) + >>> abscissas.round(4) + array([[0.0381, 0.3087, 0.6913, 0.9619]]) + >>> weights.round(4) + array([0.1321, 0.3679, 0.3679, 0.1321]) + + See also: + :func:`chaospy.quadrature.gaussian` + :func:`chaospy.quadrature.clenshaw_curtis` + :func:`chaospy.quadrature.fejer_2` + + """ + order = numpy.asarray(order) + order = numpy.where(growth, 2*3**order-1, order) + return hypercube_quadrature( + quad_func=fejer_1_simple, + order=order, + domain=domain, + segments=segments, + ) + + +@lru_cache(None) +def fejer_1_simple(order): + """Backend for Fejer type I quadrature.""" + order = int(order) + if order == 0: + return numpy.array([.5]), numpy.array([1.]) + order += 1 + + abscissas = -0.5*numpy.cos(numpy.pi*(numpy.arange(order)+0.5)/order)+0.5 + + steps = numpy.arange(1, order, 2) + length = len(steps) + remains = order-length + + kappa = numpy.arange(remains) + beta = numpy.hstack([2*numpy.exp(1j*numpy.pi*kappa/order)/(1-4*kappa**2), + numpy.zeros(length+1)]) + beta = beta[:-1]+numpy.conjugate(beta[:0:-1]) + + weights = numpy.fft.ifft(beta) + assert max(weights.imag) < 1e-15 + weights = weights.real/2. + + return abscissas, weights diff --git a/chaospy/quadrature/fejer_2.py b/chaospy/quadrature/fejer_2.py new file mode 100644 index 00000000..73e63606 --- /dev/null +++ b/chaospy/quadrature/fejer_2.py @@ -0,0 +1,76 @@ +# -*- coding: utf-8 -*- +"""Generate the quadrature abscissas and weights in Fejer quadrature.""" +import numpy +import chaospy + +from .hypercube import hypercube_quadrature +from .clenshaw_curtis import clenshaw_curtis_simple + + +def fejer_2(order, domain=(0, 1), growth=False, segments=1): + """ + Generate the quadrature abscissas and weights in Fejér type II quadrature. + + Fejér proposed two quadrature rules very similar to + :func:`chaospy.quadrature.clenshaw_curtis`. The only difference is that the + endpoints are removed. That is, Fejér only used the interior extrema of the + Chebyshev polynomials, i.e. the true stationary points. This makes this a + better method for performing quadrature on infinite intervals, as the + evaluation does not contain endpoint values. + + Args: + order (int, numpy.ndarray): + Quadrature order. + domain (chaospy.Distribution, numpy.ndarray): + Either distribution or bounding of interval to integrate over. + growth (bool): + If True sets the growth rule for the quadrature rule to only + include orders that enhances nested samples. + segments (int): + Split intervals into N subintervals and create a patched + quadrature based on the segmented quadrature. Can not be lower than + `order`. If 0 is provided, default to square root of `order`. + Nested samples only exist when the number of segments are fixed. + + Returns: + abscissas (numpy.ndarray): + The quadrature points for where to evaluate the model function + with ``abscissas.shape == (len(dist), N)`` where ``N`` is the + number of samples. + weights (numpy.ndarray): + The quadrature weights with ``weights.shape == (N,)``. + + Notes: + Implemented as proposed by Waldvogel :cite:`waldvogel_fast_2006`. + + Example: + >>> abscissas, weights = chaospy.quadrature.fejer_2(3, (0, 1)) + >>> abscissas.round(4) + array([[0.0955, 0.3455, 0.6545, 0.9045]]) + >>> weights.round(4) + array([0.1804, 0.2996, 0.2996, 0.1804]) + + See also: + :func:`chaospy.quadrature.gaussian` + :func:`chaospy.quadrature.clenshaw_curtis` + :func:`chaospy.quadrature.fejer_1` + + """ + order = numpy.asarray(order) + order = numpy.where(growth, numpy.where(order > 0, 2**(order+1)-2, 0), order) + return hypercube_quadrature( + quad_func=fejer_2_simple, + order=order, + domain=domain, + segments=segments, + ) + + +def fejer_2_simple(order): + """ + Backend for Fejer type II quadrature. + + Same as Clenshaw-Curtis, but with the end nodes removed. + """ + abscissas, weights = clenshaw_curtis_simple(order+2) + return abscissas[1:-1], weights[1:-1] diff --git a/chaospy/quadrature/frontend.py b/chaospy/quadrature/frontend.py index a6274f57..ef1cb32d 100644 --- a/chaospy/quadrature/frontend.py +++ b/chaospy/quadrature/frontend.py @@ -1,50 +1,34 @@ """Numerical quadrature node and weight generator.""" +import logging import numpy import chaospy -from .clenshaw_curtis import quad_clenshaw_curtis -from .discrete import quad_discrete -from .fejer import quad_fejer -from .gaussian import quad_gaussian -from .gauss_patterson import quad_gauss_patterson -from .gauss_legendre import quad_gauss_legendre -from .gauss_lobatto import quad_gauss_lobatto -from .gauss_kronrod import quad_gauss_kronrod -from .gauss_radau import quad_gauss_radau -from .genz_keister import quad_genz_keister -from .grid import quad_grid -from .leja import quad_leja -from .newton_cotes import quad_newton_cotes - -QUAD_NAMES = { +from .utils import combine +from .sparse_grid import sparse_grid + +SHORT_NAME_TABLE = { "c": "clenshaw_curtis", "clenshaw_curtis": "clenshaw_curtis", - "f": "fejer", "fejer": "fejer", + "f1": "fejer_1", "fejer_1": "fejer_1", + "f2": "fejer_2", "fejer_2": "fejer_2", "g": "gaussian", "gaussian": "gaussian", - "e": "gauss_legendre", "gauss_legendre": "gauss_legendre", - "l": "gauss_lobatto", "gauss_lobatto": "gauss_lobatto", - "k": "gauss_kronrod", "gauss_kronrod": "gauss_kronrod", - "p": "gauss_patterson", "gauss_patterson": "gauss_patterson", - "r": "gauss_radau", "gauss_radau": "gauss_radau", - "z": "genz_keister", "genz_keister": "genz_keister", + "e": "legendre", "legendre": "legendre", + "l": "lobatto", "lobatto": "lobatto", + "k": "kronrod", "kronrod": "kronrod", + "p": "patterson", "patterson": "patterson", + "r": "radau", "radau": "radau", "j": "leja", "leja": "leja", "n": "newton_cotes", "newton_cotes": "newton_cotes", "d": "discrete", "discrete": "discrete", "i": "grid", "grid": "grid", } -QUAD_FUNCTIONS = { - "clenshaw_curtis": quad_clenshaw_curtis, - "fejer": quad_fejer, - "gaussian": quad_gaussian, - "gauss_kronrod": quad_gauss_kronrod, - "gauss_legendre": quad_gauss_legendre, - "gauss_lobatto": quad_gauss_lobatto, - "gauss_patterson": quad_gauss_patterson, - "gauss_radau": quad_gauss_radau, - "genz_keister": quad_genz_keister, - "leja": quad_leja, - "newton_cotes": quad_newton_cotes, - "discrete": quad_discrete, - "grid": quad_grid, +DEPRECATED_SHORT_NAMES = { + "f": "f2", + "fejer": "fejer_2", + "gauss_kronrod": "kronrod", + "gauss_lobatto": "lobatto", + "gauss_patterson": "patterson", + "gauss_radau": "radau", + "gauss_legendre": "legendre", } @@ -108,12 +92,27 @@ def generate_quadrature( Examples: >>> distribution = chaospy.Iid(chaospy.Normal(0, 1), 2) >>> abscissas, weights = generate_quadrature( - ... 1, distribution, rule=("gaussian", "fejer")) + ... 1, distribution, rule=("gaussian", "fejer_2")) >>> abscissas.round(3) array([[-1. , -1. , 1. , 1. ], [-4.11, 4.11, -4.11, 4.11]]) >>> weights.round(3) - array([0.25, 0.25, 0.25, 0.25]) + array([0.222, 0.222, 0.222, 0.222]) + + See also: + :func:`chaospy.quadrature.clenshaw_curtis` + :func:`chaospy.quadrature.fejer_1` + :func:`chaospy.quadrature.fejer_2` + :func:`chaospy.quadrature.gaussian` + :func:`chaospy.quadrature.legendre_proxy` + :func:`chaospy.quadrature.lobatto` + :func:`chaospy.quadrature.kronrod` + :func:`chaospy.quadrature.patterson` + :func:`chaospy.quadrature.radau` + :func:`chaospy.quadrature.leja` + :func:`chaospy.quadrature.newton_cotes` + :func:`chaospy.quadrature.discrete` + :func:`chaospy.quadrature.grid` """ if not rule: @@ -123,7 +122,7 @@ def generate_quadrature( ] if sparse: - return chaospy.sparse_grid.construct_sparse_grid( + return sparse_grid( order=order, dist=dist, growth=growth, @@ -135,9 +134,10 @@ def generate_quadrature( ) if len(dist) == 1 or dist.stochastic_dependent: - if not isinstance(rule, str) and len(rule) == 1: + if not isinstance(rule, str) and len(set(rule)) == 1: rule = rule[0] - assert isinstance(rule, str), "dependencies require rule consistency" + assert isinstance(rule, str), ( + "dependencies require rule consistency; %s provided" % rule) abscissas, weights = _generate_quadrature( order=order, dist=dist, @@ -172,8 +172,8 @@ def generate_quadrature( ) for order_, dist_, rule_ in zip(order, dist, rule) ]) - abscissas = chaospy.combine([abscissa.T for abscissa in abscissas]).T - weights = numpy.prod(chaospy.combine([weight.T for weight in weights]), -1) + abscissas = combine([abscissa.T for abscissa in abscissas]).T + weights = numpy.prod(combine([weight.T for weight in weights]), -1) assert abscissas.shape == (len(dist), len(weights)) if dist.interpret_as_integer: @@ -183,6 +183,7 @@ def generate_quadrature( def _generate_quadrature(order, dist, rule, **kwargs): + logger = logging.getLogger(__name__) if isinstance(dist, chaospy.OperatorDistribution): args = ("left", "right") @@ -205,16 +206,30 @@ def _generate_quadrature(order, dist, rule, **kwargs): abscissas = (abscissas.T*const.T).T return abscissas, weights - rule = QUAD_NAMES[rule.lower()] + rule = rule.lower() + if rule in DEPRECATED_SHORT_NAMES: + logger.warning("quadrature rule '%s' is renamed to '%s'; " + "error will be raised in the future", + rule, DEPRECATED_SHORT_NAMES[rule]) + rule = DEPRECATED_SHORT_NAMES[rule] + rule = SHORT_NAME_TABLE[rule] + parameters = {} - if rule in ("clenshaw_curtis", "fejer", "newton_cotes", "discrete"): - parameters.update(growth=kwargs["growth"], segments=kwargs["segments"]) + if rule in ("clenshaw_curtis", "fejer_1", "fejer_2", "newton_cotes", "discrete", "grid"): + parameters["growth"] = kwargs["growth"] - if rule in ("gaussian", "gauss_kronrod", "gauss_radau", "gauss_lobatto"): - parameters.update(tolerance=kwargs["tolerance"], scaling=kwargs["scaling"], - n_max=kwargs["n_max"], recurrence_algorithm=kwargs["recurrence_algorithm"]) + if rule in ("clenshaw_curtis", "fejer_1", "fejer_2", "newton_cotes", "grid", "legendre"): + parameters["segments"] = kwargs["segments"] + + if rule in ("gaussian", "kronrod", "radau", "lobatto"): + parameters.update( + n_max=kwargs["n_max"], + tolerance=kwargs["tolerance"], + scaling=kwargs["scaling"], + recurrence_algorithm=kwargs["recurrence_algorithm"], + ) - quad_function = QUAD_FUNCTIONS[rule] + quad_function = chaospy.quadrature.INTEGRATION_COLLECTION[rule] abscissas, weights = quad_function(order, dist, **parameters) return abscissas, weights diff --git a/chaospy/quadrature/gauss_legendre.py b/chaospy/quadrature/gauss_legendre.py deleted file mode 100644 index 35b45e41..00000000 --- a/chaospy/quadrature/gauss_legendre.py +++ /dev/null @@ -1,153 +0,0 @@ -r""" -The Gauss-Legendre quadrature rule is properly supported by in :ref:`gaussian`. -However, as Gauss-Legendre is a special case where the weight function is -constant, it can in principle be used to integrate any weighting function. In -other words, this is the same Gauss-Legendre integration rule, but only in the -context of uniform distribution as weight function. Normalization of the -weights will be used to achieve the general integration form. - -It is also worth noting that this specific implementation of Gauss-Legendre is -faster to compute than the general version in :ref:`gaussian`. - -Example usage -------------- - -The first few orders:: - - >>> distribution = chaospy.Uniform(0, 1) - >>> for order in [0, 1, 2, 3]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_legendre") - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0.5]] [1.] - 1 [[0.211 0.789]] [0.5 0.5] - 2 [[0.113 0.5 0.887]] [0.278 0.444 0.278] - 3 [[0.069 0.33 0.67 0.931]] [0.174 0.326 0.326 0.174] - -Using an alternative distribution:: - - >>> distribution = chaospy.Beta(2, 4) - >>> for order in [0, 1, 2, 3]: - ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_legendre") - ... print(order, abscissas.round(3), weights.round(3)) - 0 [[0.5]] [1.] - 1 [[0.211 0.789]] [0.933 0.067] - 2 [[0.113 0.5 0.887]] [0.437 0.556 0.007] - 3 [[0.069 0.33 0.67 0.931]] [0.195 0.647 0.157 0.001] - -The abscissas stays the same, but the weights are re-adjusted for the new -weight function. -""" -import numpy -import chaospy - -from .combine import combine_quadrature - - -def quad_gauss_legendre( - order, - domain=(0, 1), - recurrence_algorithm="stieltjes", - rule="clenshaw_curtis", - tolerance=1e-10, - scaling=3, - n_max=5000, -): - r""" - Generate the quadrature nodes and weights in Gauss-Legendre quadrature. - - Note that this rule exists to allow for integrating functions with weight - functions without actually adding the quadrature. Like: - - .. math: - \int_a^b p(x) f(x) dx \approx \sum_i p(X_i) f(X_i) W_i - - instead of the more traditional: - - .. math: - \int_a^b p(x) f(x) dx \approx \sum_i f(X_i) W_i - - To get the behavior where the weight function is taken into consideration, - use :func:`chaospy.quad_gaussian`. - - Args: - order (int, numpy.ndarray): - Quadrature order. - domain (chaospy.distributions.baseclass.Distribution, numpy.ndarray): - Either distribution or bounding of interval to integrate over. - recurrence_algorithm (str): - Name of the algorithm used to generate abscissas and weights. - rule (str): - In the case of ``lanczos`` or ``stieltjes``, defines the - proxy-integration scheme. - tolerance (float): - The allowed relative error in norm between two quadrature orders - before method assumes convergence. - scaling (float): - A multiplier the adaptive order increases with for each step - quadrature order is not converged. Use 0 to indicate unit - increments. - n_max (int): - The allowed number of quadrature points to use in approximation. - - Returns: - (numpy.ndarray, numpy.ndarray): - abscissas: - The quadrature points for where to evaluate the model function - with ``abscissas.shape == (len(dist), N)`` where ``N`` is the - number of samples. - weights: - The quadrature weights with ``weights.shape == (N,)``. - - Example: - >>> abscissas, weights = quad_gauss_legendre(3) - >>> abscissas.round(4) - array([[0.0694, 0.33 , 0.67 , 0.9306]]) - >>> weights.round(4) - array([0.1739, 0.3261, 0.3261, 0.1739]) - """ - from ..distributions.baseclass import Distribution - from ..distributions.collection import Uniform - if isinstance(domain, Distribution): - abscissas, weights = quad_gauss_legendre( - order=order, - domain=(domain.lower, domain.upper), - recurrence_algorithm=recurrence_algorithm, - rule=rule, - tolerance=tolerance, - scaling=scaling, - n_max=n_max, - ) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - order = numpy.asarray(order, dtype=int).flatten() - lower, upper = numpy.array(domain) - lower = numpy.asarray(lower).flatten() - upper = numpy.asarray(upper).flatten() - - dim = max(lower.size, upper.size, order.size) - order = numpy.ones(dim, dtype=int)*order - lower = numpy.ones(dim)*lower - upper = numpy.ones(dim)*upper - - coefficients = chaospy.construct_recurrence_coefficients( - order=numpy.max(order), - dist=Uniform(0, 1), - recurrence_algorithm=recurrence_algorithm, - rule=rule, - tolerance=tolerance, - scaling=scaling, - n_max=n_max, - ) - - abscissas, weights = zip(*[chaospy.coefficients_to_quadrature( - coefficients[:order_+1]) for order_ in order]) - abscissas = list(numpy.asarray(abscissas).reshape(dim, -1)) - weights = list(numpy.asarray(weights).reshape(dim, -1)) - - return combine_quadrature(abscissas, weights, (lower, upper)) diff --git a/chaospy/quadrature/gaussian.py b/chaospy/quadrature/gaussian.py index 0a03cb63..5d4776d4 100644 --- a/chaospy/quadrature/gaussian.py +++ b/chaospy/quadrature/gaussian.py @@ -1,10 +1,10 @@ r"""Create Gaussian quadrature nodes and weights.""" import chaospy -from .combine import combine_quadrature +from .utils import combine_quadrature -def quad_gaussian( +def gaussian( order, dist, recurrence_algorithm="stieltjes", @@ -51,9 +51,9 @@ def quad_gaussian( Raises: NotImplementedError: - In the case of recurrence algorithm ``analytical``, error is raised - if the distribution does not implement the three terms recurrence - algorithm analytically. + In the case of ``analytical`` three terms recurrence algorithm, + error is raised if the distribution does not implement the feature. + coefficients. numpy.linalg.LinAlgError: For non-canonical random variables, the construction might fail because of illegal numerical operations. @@ -64,20 +64,19 @@ def quad_gaussian( Examples: >>> distribution = chaospy.Normal(0, 1) - >>> abscissas, weights = chaospy.quad_gaussian( - ... 5, distribution, recurrence_algorithm="stieltjes") + >>> abscissas, weights = chaospy.quadrature.gaussian(5, distribution) >>> abscissas.round(4) array([[-3.3243, -1.8892, -0.6167, 0.6167, 1.8892, 3.3243]]) >>> weights.round(4) array([0.0026, 0.0886, 0.4088, 0.4088, 0.0886, 0.0026]) >>> distribution = chaospy.J(chaospy.Uniform(), chaospy.Normal()) - >>> abscissas, weights = chaospy.quad_gaussian( - ... 2, distribution, recurrence_algorithm="chebyshev") + >>> abscissas, weights = chaospy.quadrature.gaussian(2, distribution) >>> abscissas.round(2) array([[ 0.11, 0.11, 0.11, 0.5 , 0.5 , 0.5 , 0.89, 0.89, 0.89], [-1.73, 0. , 1.73, -1.73, 0. , 1.73, -1.73, 0. , 1.73]]) >>> weights.round(3) array([0.046, 0.185, 0.046, 0.074, 0.296, 0.074, 0.046, 0.185, 0.046]) + """ coefficients = chaospy.construct_recurrence_coefficients( order=order, diff --git a/chaospy/quadrature/gegenbauer.py b/chaospy/quadrature/gegenbauer.py new file mode 100644 index 00000000..d26fc5fe --- /dev/null +++ b/chaospy/quadrature/gegenbauer.py @@ -0,0 +1,59 @@ +"""Gauss-Gegenbauer quadrature rule.""" +import numpy +import chaospy + +from .hypercube import hypercube_quadrature + + +def gegenbauer(order, alpha, lower=-1, upper=1, physicist=False): + """ + Gauss-Gegenbauer quadrature rule. + + Compute the sample points and weights for Gauss-Gegenbauer quadrature. The + sample points are the roots of the nth degree Gegenbauer polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + Gauss-Gegenbauer physicist means a weight function + :math:`(1-x^2)^{\alpha-0.5}` and weights that sum to :math`2^{2\alpha-1}`, + and probabilist means a weight function is + :math:`B(\alpha+0.5, \alpha+0.5) (x-x^2)^{\alpha+1/2}` (where :math:`B` is + the beta normalizing constant) which sum to 1. + + Args: + order (int): + The quadrature order. + alpha (float): + Gegenbauer shape parameter. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.gegenbauer(3, alpha=2) + >>> abscissas + array([[-0.72741239, -0.26621648, 0.26621648, 0.72741239]]) + >>> weights + array([0.10452141, 0.39547859, 0.39547859, 0.10452141]) + + See also: + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Beta(alpha+0.5, alpha+0.5, lower, upper)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights *= 2**(2*alpha-1) if physicist else 1 + return abscissas[numpy.newaxis], weights diff --git a/chaospy/quadrature/genz_keister/__init__.py b/chaospy/quadrature/genz_keister/__init__.py deleted file mode 100644 index e810e8e4..00000000 --- a/chaospy/quadrature/genz_keister/__init__.py +++ /dev/null @@ -1,11 +0,0 @@ -""" -Hermite Genz-Keister quadrature rules - -Adapted from John Burkardt's implementation in Matlab -""" -from .frontend import quad_genz_keister, GENS_KEISTER_FUNCTIONS - -from .gk16 import quad_genz_keister_16 -from .gk18 import quad_genz_keister_18 -from .gk22 import quad_genz_keister_22 -from .gk24 import quad_genz_keister_24 diff --git a/chaospy/quadrature/genz_keister/frontend.py b/chaospy/quadrature/genz_keister/frontend.py deleted file mode 100644 index 25269faf..00000000 --- a/chaospy/quadrature/genz_keister/frontend.py +++ /dev/null @@ -1,56 +0,0 @@ -""" -Frontend for the Hermite Genz-Keister quadrature rule. -""" -import numpy -import scipy.special - -from ..combine import combine -from .gk16 import quad_genz_keister_16 -from .gk18 import quad_genz_keister_18 -from .gk22 import quad_genz_keister_22 -from .gk24 import quad_genz_keister_24 - -GENS_KEISTER_FUNCTIONS = { - 16: quad_genz_keister_16, - 18: quad_genz_keister_18, - 22: quad_genz_keister_22, - 24: quad_genz_keister_24, -} - - -def quad_genz_keister(order, dist, rule=24): - """ - Genz-Keister quadrature rule. - - Examples: - >>> abscissas, weights = quad_genz_keister( - ... order=1, dist=chaospy.Iid(chaospy.Uniform(0, 1), 2)) - >>> abscissas.round(2) - array([[0.04, 0.04, 0.04, 0.5 , 0.5 , 0.5 , 0.96, 0.96, 0.96], - [0.04, 0.5 , 0.96, 0.04, 0.5 , 0.96, 0.04, 0.5 , 0.96]]) - >>> weights.round(2) - array([0.03, 0.11, 0.03, 0.11, 0.44, 0.11, 0.03, 0.11, 0.03]) - - """ - assert isinstance(rule, int) - - if len(dist) > 1: - - if isinstance(order, int): - values = [quad_genz_keister(order, d, rule) for d in dist] - else: - values = [quad_genz_keister(order[i], dist[i], rule) - for i in range(len(dist))] - - abscissas = [_[0][0] for _ in values] - abscissas = combine(abscissas).T - weights = [_[1] for _ in values] - weights = numpy.prod(combine(weights), -1) - - return abscissas, weights - - foo = GENS_KEISTER_FUNCTIONS[rule] - abscissas, weights = foo(order) - abscissas = dist.inv(scipy.special.ndtr(abscissas)) - abscissas = abscissas.reshape(1, abscissas.size) - return abscissas, weights diff --git a/chaospy/quadrature/genz_keister/gk16.py b/chaospy/quadrature/genz_keister/gk16.py deleted file mode 100644 index 2b615bcc..00000000 --- a/chaospy/quadrature/genz_keister/gk16.py +++ /dev/null @@ -1,374 +0,0 @@ -"""Hermite Genz-Keister 16 rule.""" - -import numpy - -def quad_genz_keister_16(order): - """ - Hermite Genz-Keister 16 rule. - - Args: - order (int): - The quadrature order. Must be in the interval (0, 8). - - Returns: - (:py:data:typing.Tuple[numpy.ndarray, numpy.ndarray]): - Abscissas and weights - - Examples: - >>> abscissas, weights = quad_genz_keister_16(1) - >>> abscissas.round(4) - array([-1.7321, 0. , 1.7321]) - >>> weights.round(4) - array([0.1667, 0.6667, 0.1667]) - """ - order = sorted(GENZ_KEISTER_16.keys())[order] - - abscissas, weights = GENZ_KEISTER_16[order] - abscissas = numpy.array(abscissas) - weights = numpy.array(weights) - - weights /= numpy.sum(weights) - abscissas *= numpy.sqrt(2) - - return abscissas, weights - - -GENZ_KEISTER_16 = { - 1 : (( - 0.0000000000000000E+00, - ), ( - 1.7724538509055159E+00, - )), - 3 : (( - -1.2247448713915889E+00, - 0.0000000000000000E+00, - 1.2247448713915889E+00, - ), ( - 2.9540897515091930E-01, - 1.1816359006036772E+00, - 2.9540897515091930E-01, - )), - 7 : (( - -2.9592107790638380E+00, - -1.2247448713915889E+00, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 1.2247448713915889E+00, - 2.9592107790638380E+00, - ), ( - 1.2330680655153448E-03, - 2.4557928535031393E-01, - 2.3286251787386100E-01, - 8.1310410832613500E-01, - 2.3286251787386100E-01, - 2.4557928535031393E-01, - 1.2330680655153448E-03, - )), - 9 : (( - -2.9592107790638380E+00, - -2.0232301911005157E+00, - -1.2247448713915889E+00, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 1.2247448713915889E+00, - 2.0232301911005157E+00, - 2.9592107790638380E+00, - ), ( - 1.6708826306882348E-04, - 1.4173117873979098E-02, - 1.6811892894767771E-01, - 4.7869428549114124E-01, - 4.5014700975378197E-01, - 4.7869428549114124E-01, - 1.6811892894767771E-01, - 1.4173117873979098E-02, - 1.6708826306882348E-04, - )), - 17 : (( - -4.4995993983103881E+00, - -3.6677742159463378E+00, - -2.9592107790638380E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.9592107790638380E+00, - 3.6677742159463378E+00, - 4.4995993983103881E+00, - ), ( - 3.7463469943051758E-08, - -1.4542843387069391E-06, - 1.8723818949278350E-04, - 1.2466519132805918E-02, - 3.4840719346803800E-03, - 1.5718298376652240E-01, - 2.5155825701712934E-02, - 4.5119803602358544E-01, - 4.7310733504965385E-01, - 4.5119803602358544E-01, - 2.5155825701712934E-02, - 1.5718298376652240E-01, - 3.4840719346803800E-03, - 1.2466519132805918E-02, - 1.8723818949278350E-04, - -1.4542843387069391E-06, - 3.7463469943051758E-08, - )), - 19 : (( - -4.4995993983103881E+00, - -3.6677742159463378E+00, - -2.9592107790638380E+00, - -2.2665132620567876E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.2665132620567876E+00, - 2.9592107790638380E+00, - 3.6677742159463378E+00, - 4.4995993983103881E+00, - ), ( - 1.5295717705322357E-09, - 1.0802767206624762E-06, - 1.0656589772852267E-04, - 5.1133174390883855E-03, - -1.1232438489069229E-02, - 3.2055243099445879E-02, - 1.1360729895748269E-01, - 1.0838861955003017E-01, - 3.6924643368920851E-01, - 5.3788160700510168E-01, - 3.6924643368920851E-01, - 1.0838861955003017E-01, - 1.1360729895748269E-01, - 3.2055243099445879E-02, - -1.1232438489069229E-02, - 5.1133174390883855E-03, - 1.0656589772852267E-04, - 1.0802767206624762E-06, - 1.5295717705322357E-09, - )), - 31 : (( - -6.3759392709822356E+00, - -5.6432578578857449E+00, - -5.0360899444730940E+00, - -4.4995993983103881E+00, - -3.6677742159463378E+00, - -2.9592107790638380E+00, - -2.5705583765842968E+00, - -2.2665132620567876E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.5794121348467671E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - -1.7606414208200893E-01, - 0.0000000000000000E+00, - 1.7606414208200893E-01, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.5794121348467671E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.2665132620567876E+00, - 2.5705583765842968E+00, - 2.9592107790638380E+00, - 3.6677742159463378E+00, - 4.4995993983103881E+00, - 5.0360899444730940E+00, - 5.6432578578857449E+00, - 6.3759392709822356E+00, - ), ( - 2.2365645607044459E-15, - -2.6304696458548942E-13, - 9.0675288231679823E-12, - 1.4055252024722478E-09, - 1.0889219692128120E-06, - 1.0541662394746661E-04, - 2.6665159778939428E-05, - 4.8385208205502612E-03, - -9.8566270434610019E-03, - 2.9409427580350787E-02, - 3.1210210352682834E-03, - 1.0939325071860877E-01, - 1.1594930984853116E-01, - 3.5393889029580544E-01, - 4.9855761893293160E-02, - 4.5888839636756751E-01, - 4.9855761893293160E-02, - 3.5393889029580544E-01, - 1.1594930984853116E-01, - 1.0939325071860877E-01, - 3.1210210352682834E-03, - 2.9409427580350787E-02, - -9.8566270434610019E-03, - 4.8385208205502612E-03, - 2.6665159778939428E-05, - 1.0541662394746661E-04, - 1.0889219692128120E-06, - 1.4055252024722478E-09, - 9.0675288231679823E-12, - -2.6304696458548942E-13, - 2.2365645607044459E-15, - )), - 33 : (( - -6.3759392709822356E+00, - -5.6432578578857449E+00, - -5.0360899444730940E+00, - -4.4995993983103881E+00, - -4.0292201405043713E+00, - -3.6677742159463378E+00, - -2.9592107790638380E+00, - -2.5705583765842968E+00, - -2.2665132620567876E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.5794121348467671E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - -1.7606414208200893E-01, - 0.0000000000000000E+00, - 1.7606414208200893E-01, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.5794121348467671E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.2665132620567876E+00, - 2.5705583765842968E+00, - 2.9592107790638380E+00, - 3.6677742159463378E+00, - 4.0292201405043713E+00, - 4.4995993983103881E+00, - 5.0360899444730940E+00, - 5.6432578578857449E+00, - 6.3759392709822356E+00, - ), ( - -1.7602932805372496E-15, - 4.7219278666417693E-13, - -3.4281570530349562E-11, - 2.7547825138935901E-09, - -2.3903343382803510E-08, - 1.2245220967158438E-06, - 9.8710009197409173E-05, - 1.4753204901862772E-04, - 3.7580026604304793E-03, - -4.9118576123877555E-03, - 2.0435058359107205E-02, - 1.3032872699027960E-02, - 9.6913444944583621E-02, - 1.3726521191567551E-01, - 3.1208656194697448E-01, - 1.8411696047725790E-01, - 2.4656644932829619E-01, - 1.8411696047725790E-01, - 3.1208656194697448E-01, - 1.3726521191567551E-01, - 9.6913444944583621E-02, - 1.3032872699027960E-02, - 2.0435058359107205E-02, - -4.9118576123877555E-03, - 3.7580026604304793E-03, - 1.4753204901862772E-04, - 9.8710009197409173E-05, - 1.2245220967158438E-06, - -2.3903343382803510E-08, - 2.7547825138935901E-09, - -3.4281570530349562E-11, - 4.7219278666417693E-13, - -1.7602932805372496E-15, - )), - 35 : (( - -6.3759392709822356E+00, - -5.6432578578857449E+00, - -5.0360899444730940E+00, - -4.4995993983103881E+00, - -4.0292201405043713E+00, - -3.6677742159463378E+00, - -3.3491639537131945E+00, - -2.9592107790638380E+00, - -2.5705583765842968E+00, - -2.2665132620567876E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.5794121348467671E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - -1.7606414208200893E-01, - 0.0000000000000000E+00, - 1.7606414208200893E-01, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.5794121348467671E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.2665132620567876E+00, - 2.5705583765842968E+00, - 2.9592107790638380E+00, - 3.3491639537131945E+00, - 3.6677742159463378E+00, - 4.0292201405043713E+00, - 4.4995993983103881E+00, - 5.0360899444730940E+00, - 5.6432578578857449E+00, - 6.3759392709822356E+00, - ), ( - 1.8684014894510604E-18, - 9.6599466278563243E-15, - 5.4896836948499462E-12, - 8.1553721816916897E-10, - 3.7920222392319532E-08, - 4.3737818040926989E-07, - 4.8462799737020461E-06, - 6.3328620805617891E-05, - 4.8785399304443770E-04, - 1.4515580425155904E-03, - 4.0967527720344047E-03, - 5.5928828911469180E-03, - 2.7780508908535097E-02, - 8.0245518147390893E-02, - 1.6371221555735804E-01, - 2.6244871488784277E-01, - 3.3988595585585218E-01, - 9.1262675363737921E-04, - 3.3988595585585218E-01, - 2.6244871488784277E-01, - 1.6371221555735804E-01, - 8.0245518147390893E-02, - 2.7780508908535097E-02, - 5.5928828911469180E-03, - 4.0967527720344047E-03, - 1.4515580425155904E-03, - 4.8785399304443770E-04, - 6.3328620805617891E-05, - 4.8462799737020461E-06, - 4.3737818040926989E-07, - 3.7920222392319532E-08, - 8.1553721816916897E-10, - 5.4896836948499462E-12, - 9.6599466278563243E-15, - 1.8684014894510604E-18, - )), -} diff --git a/chaospy/quadrature/genz_keister/gk18.py b/chaospy/quadrature/genz_keister/gk18.py deleted file mode 100644 index e6d4a2fb..00000000 --- a/chaospy/quadrature/genz_keister/gk18.py +++ /dev/null @@ -1,190 +0,0 @@ -"""Hermite Genz-Keister 18 rule.""" -import numpy - - -def quad_genz_keister_18(order): - """ - Hermite Genz-Keister 18 rule. - - Args: - order (int): - The quadrature order. Must be in the interval (0, 8). - - Returns: - (:py:data:typing.Tuple[numpy.ndarray, numpy.ndarray]): - Abscissas and weights - - Examples: - >>> abscissas, weights = quad_genz_keister_18(1) - >>> abscissas.round(4) - array([-1.7321, 0. , 1.7321]) - >>> weights.round(4) - array([0.1667, 0.6667, 0.1667]) - """ - order = sorted(GENZ_KEISTER_18.keys())[order] - - abscissas, weights = GENZ_KEISTER_18[order] - abscissas = numpy.array(abscissas) - weights = numpy.array(weights) - - weights /= numpy.sum(weights) - abscissas *= numpy.sqrt(2) - - return abscissas, weights - - -GENZ_KEISTER_18 = { - 1 : (( - 0.0000000000000000E+00, - ), ( - 1.7724538509055159E+00, - )), - 3 : (( - -1.2247448713915889E+00, - 0.0000000000000000E+00, - 1.2247448713915889E+00, - ), ( - 2.9540897515091930E-01, - 1.1816359006036772E+00, - 2.9540897515091930E-01, - )), - 9 : (( - -2.9592107790638380E+00, - -2.0232301911005157E+00, - -1.2247448713915889E+00, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 1.2247448713915889E+00, - 2.0232301911005157E+00, - 2.9592107790638380E+00, - ), ( - 1.6708826306882348E-04, - 1.4173117873979098E-02, - 1.6811892894767771E-01, - 4.7869428549114124E-01, - 4.5014700975378197E-01, - 4.7869428549114124E-01, - 1.6811892894767771E-01, - 1.4173117873979098E-02, - 1.6708826306882348E-04, - )), - 19 : (( - -4.4995993983103881E+00, - -3.6677742159463378E+00, - -2.9592107790638380E+00, - -2.2665132620567876E+00, - -2.0232301911005157E+00, - -1.8357079751751868E+00, - -1.2247448713915889E+00, - -8.7004089535290285E-01, - -5.2403354748695763E-01, - 0.0000000000000000E+00, - 5.2403354748695763E-01, - 8.7004089535290285E-01, - 1.2247448713915889E+00, - 1.8357079751751868E+00, - 2.0232301911005157E+00, - 2.2665132620567876E+00, - 2.9592107790638380E+00, - 3.6677742159463378E+00, - 4.4995993983103881E+00, - ), ( - 1.5295717705322357E-09, - 1.0802767206624762E-06, - 1.0656589772852267E-04, - 5.1133174390883855E-03, - -1.1232438489069229E-02, - 3.2055243099445879E-02, - 1.1360729895748269E-01, - 1.0838861955003017E-01, - 3.6924643368920851E-01, - 5.3788160700510168E-01, - 3.6924643368920851E-01, - 1.0838861955003017E-01, - 1.1360729895748269E-01, - 3.2055243099445879E-02, - -1.1232438489069229E-02, - 5.1133174390883855E-03, - 1.0656589772852267E-04, - 1.0802767206624762E-06, - 1.5295717705322357E-09, - )), - 37 : (( - -6.853200069757519, - -6.124527854622158, - -5.521865209868350, - -4.986551454150765, - -4.499599398310388, - -4.057956316089741, - -3.667774215946338, - -3.315584617593290, - -2.959210779063838, - -2.597288631188366, - -2.266513262056788, - -2.023230191100516, - -1.835707975175187, - -1.561553427651873, - -1.224744871391589, - -0.870040895352903, - -0.524033547486958, - -0.214618180588171, - 0.000000000000000, - 0.214618180588171, - 0.524033547486958, - 0.870040895352903, - 1.224744871391589, - 1.561553427651873, - 1.835707975175187, - 2.023230191100516, - 2.266513262056788, - 2.597288631188366, - 2.959210779063838, - 3.315584617593290, - 3.667774215946338, - 4.057956316089741, - 4.499599398310388, - 4.986551454150765, - 5.521865209868350, - 6.124527854622158, - 6.853200069757519, - ), ( - 0.19030350940130498E-20, - 0.187781893143728947E-16, - 0.182242751549129356E-13, - 0.45661763676186859E-11, - 0.422525843963111041E-09, - 0.16595448809389819E-07, - 0.295907520230744049E-06, - 0.330975870979203419E-05, - 0.32265185983739747E-04, - 0.234940366465975222E-03, - 0.985827582996483824E-03, - 0.176802225818295443E-02, - 0.43334988122723492E-02, - 0.15513109874859354E-01, - 0.442116442189845444E-01, - 0.937208280655245902E-01, - 0.143099302896833389E+00, - 0.147655710402686249E+00, - 0.968824552928425499E-01, - 0.147655710402686249E+00, - 0.143099302896833389E+00, - 0.937208280655245902E-01, - 0.442116442189845444E-01, - 0.15513109874859354E-01, - 0.43334988122723492E-02, - 0.176802225818295443E-02, - 0.985827582996483824E-03, - 0.234940366465975222E-03, - 0.32265185983739747E-04, - 0.330975870979203419E-05, - 0.295907520230744049E-06, - 0.16595448809389819E-07, - 0.422525843963111041E-09, - 0.45661763676186859E-11, - 0.182242751549129356E-13, - 0.187781893143728947E-16, - 0.19030350940130498E-20, - )), -} diff --git a/chaospy/quadrature/genz_keister/gk22.py b/chaospy/quadrature/genz_keister/gk22.py deleted file mode 100644 index e95c921a..00000000 --- a/chaospy/quadrature/genz_keister/gk22.py +++ /dev/null @@ -1,198 +0,0 @@ -"""Hermite Genz-Keister 22 rule.""" -import numpy - - -def quad_genz_keister_22(order): - """ - Hermite Genz-Keister 22 rule. - - Args: - order (int): - The quadrature order. Must be in the interval (0, 8). - - Returns: - (:py:data:typing.Tuple[numpy.ndarray, numpy.ndarray]): - Abscissas and weights - - Examples: - >>> abscissas, weights = quad_genz_keister_22(1) - >>> abscissas.round(4) - array([-1.7321, 0. , 1.7321]) - >>> weights.round(4) - array([0.1667, 0.6667, 0.1667]) - """ - order = sorted(GENZ_KEISTER_22.keys())[order] - - abscissas, weights = GENZ_KEISTER_22[order] - abscissas = numpy.array(abscissas) - weights = numpy.array(weights) - - weights /= numpy.sum(weights) - abscissas *= numpy.sqrt(2) - - return abscissas, weights - - -GENZ_KEISTER_22 = { - 1 : (( - 0.0000000000000000, - ), ( - 1.7724538509055159E+00, - )), - 3 : (( - -1.2247448713915889, - 0.0000000000000000, - 1.2247448713915889, - ), ( - 2.9540897515091930E-01, - 1.1816359006036772E+00, - 2.9540897515091930E-01, - )), - 9 : (( - -2.9592107790638380, - -2.0232301911005157, - -1.2247448713915889, - -0.52403354748695763, - 0.0000000000000000, - 0.52403354748695763, - 1.2247448713915889, - 2.0232301911005157, - 2.9592107790638380, - ), ( - 1.6708826306882348E-04, - 1.4173117873979098E-02, - 1.6811892894767771E-01, - 4.7869428549114124E-01, - 4.5014700975378197E-01, - 4.7869428549114124E-01, - 1.6811892894767771E-01, - 1.4173117873979098E-02, - 1.6708826306882348E-04, - )), - 19 : (( - -4.4995993983103881, - -3.6677742159463378, - -2.9592107790638380, - -2.2665132620567876, - -2.0232301911005157, - -1.8357079751751868, - -1.2247448713915889, - -0.87004089535290285, - -0.52403354748695763, - 0.0000000000000000, - 0.52403354748695763, - 0.87004089535290285, - 1.2247448713915889, - 1.8357079751751868, - 2.0232301911005157, - 2.2665132620567876, - 2.9592107790638380, - 3.6677742159463378, - 4.4995993983103881, - ), ( - 1.5295717705322357E-09, - 1.0802767206624762E-06, - 1.0656589772852267E-04, - 5.1133174390883855E-03, - -1.1232438489069229E-02, - 3.2055243099445879E-02, - 1.1360729895748269E-01, - 1.0838861955003017E-01, - 3.6924643368920851E-01, - 5.3788160700510168E-01, - 3.6924643368920851E-01, - 1.0838861955003017E-01, - 1.1360729895748269E-01, - 3.2055243099445879E-02, - -1.1232438489069229E-02, - 5.1133174390883855E-03, - 1.0656589772852267E-04, - 1.0802767206624762E-06, - 1.5295717705322357E-09, - )), - 41 : (( - -7.251792998192644, - -6.547083258397540, - -5.961461043404500, - -5.437443360177798, - -4.953574342912980, - -4.4995993983103881, - -4.070919267883068, - -3.6677742159463378, - -3.296114596212218, - -2.9592107790638380, - -2.630415236459871, - -2.2665132620567876, - -2.043834754429505, - -2.0232301911005157, - -1.8357079751751868, - -1.585873011819188, - -1.2247448713915889, - -0.87004089535290285, - -0.52403354748695763, - -0.195324784415805, - 0.0000000000000000, - 0.195324784415805, - 0.52403354748695763, - 0.87004089535290285, - 1.2247448713915889, - 1.585873011819188, - 1.8357079751751868, - 2.0232301911005157, - 2.043834754429505, - 2.2665132620567876, - 2.630415236459871, - 2.9592107790638380, - 3.296114596212218, - 3.6677742159463378, - 4.070919267883068, - 4.4995993983103881, - 4.953574342912980, - 5.437443360177798, - 5.961461043404500, - 6.547083258397540, - 7.251792998192644, - ), ( - 0.664195893812757801E-23, - 0.860427172512207236E-19, - 0.1140700785308509E-15, - 0.408820161202505983E-13, - 0.581803393170320419E-11, - 0.400784141604834759E-09, - 0.149158210417831408E-07, - 0.315372265852264871E-06, - 0.381182791749177506E-05, - 0.288976780274478689E-04, - 0.189010909805097887E-03, - 0.140697424065246825E-02, - - 0.144528422206988237E-01, - 0.178852543033699732E-01, - 0.705471110122962612E-03, - 0.165445526705860772E-01, - 0.45109010335859128E-01, - 0.928338228510111845E-01, - 0.145966293895926429E+00, - 0.165639740400529554E+00, - 0.562793426043218877E-01, - 0.165639740400529554E+00, - 0.145966293895926429E+00, - 0.928338228510111845E-01, - 0.45109010335859128E-01, - 0.165445526705860772E-01, - 0.705471110122962612E-03, - 0.178852543033699732E-01, - - 0.144528422206988237E-01, - 0.140697424065246825E-02, - 0.189010909805097887E-03, - 0.288976780274478689E-04, - 0.381182791749177506E-05, - 0.315372265852264871E-06, - 0.149158210417831408E-07, - 0.400784141604834759E-09, - 0.581803393170320419E-11, - 0.408820161202505983E-13, - 0.1140700785308509E-15, - 0.860427172512207236E-19, - 0.664195893812757801E-23, - )), -} diff --git a/chaospy/quadrature/genz_keister/gk24.py b/chaospy/quadrature/genz_keister/gk24.py deleted file mode 100644 index 1769442d..00000000 --- a/chaospy/quadrature/genz_keister/gk24.py +++ /dev/null @@ -1,202 +0,0 @@ -"""Hermite Genz-Keister 24 rule.""" -import numpy - - -def quad_genz_keister_24(order): - """ - Hermite Genz-Keister 24 rule. - - Args: - order (int): - The quadrature order. Must be in the interval (0, 8). - - Returns: - (:py:data:typing.Tuple[numpy.ndarray, numpy.ndarray]): - Abscissas and weights - - Examples: - >>> abscissas, weights = quad_genz_keister_24(1) - >>> abscissas.round(4) - array([-1.7321, 0. , 1.7321]) - >>> weights.round(4) - array([0.1667, 0.6667, 0.1667]) - """ - order = sorted(GENZ_KEISTER_24.keys())[order] - - abscissas, weights = GENZ_KEISTER_24[order] - abscissas = numpy.array(abscissas) - weights = numpy.array(weights) - - weights /= numpy.sum(weights) - abscissas *= numpy.sqrt(2) - - return abscissas, weights - - -GENZ_KEISTER_24 = { - 1 : (( - 0.0000000000000000, - ), ( - 1.7724538509055159E+00, - )), - 3 : (( - -1.2247448713915889, - 0.0000000000000000, - 1.2247448713915889, - ), ( - 2.9540897515091930E-01, - 1.1816359006036772E+00, - 2.9540897515091930E-01, - )), - 9 : (( - -2.9592107790638380, - -2.0232301911005157, - -1.2247448713915889, - -0.52403354748695763, - 0.0000000000000000, - 0.52403354748695763, - 1.2247448713915889, - 2.0232301911005157, - 2.9592107790638380, - ), ( - 1.6708826306882348E-04, - 1.4173117873979098E-02, - 1.6811892894767771E-01, - 4.7869428549114124E-01, - 4.5014700975378197E-01, - 4.7869428549114124E-01, - 1.6811892894767771E-01, - 1.4173117873979098E-02, - 1.6708826306882348E-04, - )), - 19 : (( - -4.4995993983103881, - -3.6677742159463378, - -2.9592107790638380, - -2.2665132620567876, - -2.0232301911005157, - -1.8357079751751868, - -1.2247448713915889, - -0.87004089535290285, - -0.52403354748695763, - 0.0000000000000000, - 0.52403354748695763, - 0.87004089535290285, - 1.2247448713915889, - 1.8357079751751868, - 2.0232301911005157, - 2.2665132620567876, - 2.9592107790638380, - 3.6677742159463378, - 4.4995993983103881, - ), ( - 1.5295717705322357E-09, - 1.0802767206624762E-06, - 1.0656589772852267E-04, - 5.1133174390883855E-03, - -1.1232438489069229E-02, - 3.2055243099445879E-02, - 1.1360729895748269E-01, - 1.0838861955003017E-01, - 3.6924643368920851E-01, - 5.3788160700510168E-01, - 3.6924643368920851E-01, - 1.0838861955003017E-01, - 1.1360729895748269E-01, - 3.2055243099445879E-02, - -1.1232438489069229E-02, - 5.1133174390883855E-03, - 1.0656589772852267E-04, - 1.0802767206624762E-06, - 1.5295717705322357E-09, - )), - 43 : (( - -10.167574994881873, - -7.231746029072501, - -6.535398426382995, - -5.954781975039809, - -5.434053000365068, - -4.952329763008589, - -4.4995993983103881, - -4.071335874253583, - -3.6677742159463378, - -3.295265921534226, - -2.9592107790638380, - -2.633356763661946, - -2.2665132620567876, - -2.089340389294661, - -2.0232301911005157, - -1.8357079751751868, - -1.583643465293944, - -1.2247448713915889, - -0.87004089535290285, - -0.52403354748695763, - -0.196029453662011, - 0.0000000000000000, - 0.196029453662011, - 0.52403354748695763, - 0.87004089535290285, - 1.2247448713915889, - 1.583643465293944, - 1.8357079751751868, - 2.0232301911005157, - 2.089340389294661, - 2.2665132620567876, - 2.633356763661946, - 2.9592107790638380, - 3.295265921534226, - 3.6677742159463378, - 4.071335874253583, - 4.4995993983103881, - 4.952329763008589, - 5.434053000365068, - 5.954781975039809, - 6.535398426382995, - 7.231746029072501, - 10.167574994881873, - ), ( - 0.546191947478318097E-37, - 0.87544909871323873E-23, - 0.992619971560149097E-19, - 0.122619614947864357E-15, - 0.421921851448196032E-13, - 0.586915885251734856E-11, - 0.400030575425776948E-09, - 0.148653643571796457E-07, - 0.316018363221289247E-06, - 0.383880761947398577E-05, - 0.286802318064777813E-04, - 0.184789465688357423E-03, - 0.150909333211638847E-02, - - 0.38799558623877157E-02, - 0.67354758901013295E-02, - 0.139966252291568061E-02, - 0.163616873493832402E-01, - 0.450612329041864976E-01, - 0.928711584442575456E-01, - 0.145863292632147353E+00, - 0.164880913687436689E+00, - 0.579595986101181095E-01, - 0.164880913687436689E+00, - 0.145863292632147353E+00, - 0.928711584442575456E-01, - 0.450612329041864976E-01, - 0.163616873493832402E-01, - 0.139966252291568061E-02, - 0.67354758901013295E-02, - - 0.38799558623877157E-02, - 0.150909333211638847E-02, - 0.184789465688357423E-03, - 0.286802318064777813E-04, - 0.383880761947398577E-05, - 0.316018363221289247E-06, - 0.148653643571796457E-07, - 0.400030575425776948E-09, - 0.586915885251734856E-11, - 0.421921851448196032E-13, - 0.122619614947864357E-15, - 0.992619971560149097E-19, - 0.87544909871323873E-23, - 0.546191947478318097E-37, - )) -} diff --git a/chaospy/quadrature/grid.py b/chaospy/quadrature/grid.py index 78078aec..c3d5bb5f 100644 --- a/chaospy/quadrature/grid.py +++ b/chaospy/quadrature/grid.py @@ -7,10 +7,10 @@ import numpy import chaospy -from .combine import combine_quadrature +from .hypercube import hypercube_quadrature -def quad_grid(order, domain=(0, 1)): +def grid(order, domain=(0, 1), growth=False, segments=1): """ Generate the quadrature abscissas and weights for simple grid. @@ -27,12 +27,12 @@ def quad_grid(order, domain=(0, 1)): The weights are all equal to `1/len(weights[0])`. Example: - >>> abscissas, weights = chaospy.quad_grid(4, chaospy.Uniform(-1, 1)) + >>> abscissas, weights = chaospy.quadrature.grid(4, chaospy.Uniform(-1, 1)) >>> abscissas.round(4) array([[-0.8, -0.4, 0. , 0.4, 0.8]]) >>> weights.round(4) array([0.2, 0.2, 0.2, 0.2, 0.2]) - >>> abscissas, weights = chaospy.quad_grid([1, 1]) + >>> abscissas, weights = chaospy.quadrature.grid([1, 1]) >>> abscissas.round(4) array([[0.25, 0.25, 0.75, 0.75], [0.25, 0.75, 0.25, 0.75]]) @@ -40,19 +40,23 @@ def quad_grid(order, domain=(0, 1)): array([0.25, 0.25, 0.25, 0.25]) """ - if isinstance(domain, chaospy.Distribution): - abscissas, weights = quad_grid(order, (domain.lower, domain.upper)) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - order = numpy.atleast_1d(order) - order, lower, upper = numpy.broadcast_arrays(order, domain[0], domain[1]) - assert order.ndim == 1, "too many dimensions" - abscissas = tuple(numpy.linspace(0, 1, 2*order_+3)[1::2] for order_ in order) - weights = tuple(numpy.repeat(1./(order_+1), order_+1) for order_ in order) - - abscissas_, weights_ = combine_quadrature(abscissas, weights, (lower, upper)) - return abscissas_, weights_ + order = numpy.asarray(order) + order = numpy.where(growth, numpy.where(order > 0, 3**order-1, 0), order) + return hypercube_quadrature( + quad_func=grid_simple, + order=order, + domain=domain, + segments=segments, + ) + + +def grid_simple(order): + """ + Backend for grid quadrature. + + Use :func:`chaospy.quadrature.grid` instead. + """ + order = int(order) + abscissas = numpy.linspace(0, 1, 2*order+3)[1::2] + weights = numpy.full(order+1, 1./(order+1)) + return abscissas, weights diff --git a/chaospy/quadrature/hermite.py b/chaospy/quadrature/hermite.py new file mode 100644 index 00000000..16d9870b --- /dev/null +++ b/chaospy/quadrature/hermite.py @@ -0,0 +1,59 @@ +"""Gauss-Hermite quadrature rule.""" +import numpy +import chaospy + +from .hypercube import hypercube_quadrature + + +def hermite(order, mu=0., sigma=1., physicist=False): + r""" + Gauss-Hermite quadrature rule. + + Compute the sample points and weights for Gauss-Hermite quadrature. The + sample points are the roots of the nth degree Hermite polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + Gauss-Hermite physicist means a weight function :math:`e^{-x^2}` and + weights that sum to :math`\sqrt(\pi)`, and probabilist means a weight + function is :math:`e^{-x^2/2}` and sum to 1. + + Args: + order (int): + The quadrature order. + mu (float): + Non-centrality parameter. + sigma (float): + Scale parameter. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.hermite(3) + >>> abscissas + array([[-2.33441422, -0.74196378, 0.74196378, 2.33441422]]) + >>> weights + array([0.04587585, 0.45412415, 0.45412415, 0.04587585]) + + See also: + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + sigma = float(sigma*2**-0.5 if physicist else sigma) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Normal(0, sigma)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights = weights*numpy.pi**0.5 if physicist else weights + if order%2 == 0: + abscissas[len(abscissas)//2] = 0 + abscissas += mu + return abscissas[numpy.newaxis], weights diff --git a/chaospy/quadrature/hypercube.py b/chaospy/quadrature/hypercube.py new file mode 100644 index 00000000..c7169941 --- /dev/null +++ b/chaospy/quadrature/hypercube.py @@ -0,0 +1,411 @@ +from functools import partial + +import numpy +import chaospy + +from .utils import combine_quadrature + + +def hypercube_quadrature( + quad_func, + order, + domain, + segments=None, + auto_scale=True, +): + """ + Enhance simple 1-dimensional unit quadrature with extra features. + + These features include handling of: + + * Distribution as domain by embedding density into the weights + * Scale to any intervals + * Multivariate support + * Repeat quadrature into segments + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights. If + ``auto_scale`` is true, the function should be on the form: + ``abscissas, weights = quad_func(order)`` and be defined on the unit + interval. Otherwise the call signature should be + ``abscissas, weights = quad_func(order, lower, upper)`` and it + should be defined on interval bound by ``lower`` and ``upper``. + order (int, Sequence[int]): + The quadrature order passed to the quadrature function. + domain (Tuple[float, float], :class:`chaospy.Distribution`): + Either interval on the format ``(lower, upper)`` or a distribution + to integrate over. If the latter, weights are adjusted to + incorporate the density at abscissas. + segments (Optional[int], Sequence[float]): + The number segments to split the interval on. If sequence is + provided, use as segment edges instead. + kwargs (Any): + Extra keyword arguments passed to `quad_func`. + + Returns: + Same as ``quad_func`` but adjusted to incorporate extra features. + + Examples: + >>> def my_quad(order): + ... return (numpy.linspace(0, 1, order+1), + ... 1./numpy.full(order+1, order+2)) + >>> my_quad(2) + (array([0. , 0.5, 1. ]), array([0.25, 0.25, 0.25])) + >>> hypercube_quadrature(my_quad, 2, domain=chaospy.Uniform(-1, 1)) + (array([[-1., 0., 1.]]), array([0.25, 0.25, 0.25])) + >>> abscissas, weights = hypercube_quadrature(my_quad, (1, 1), domain=(0, 1)) + >>> abscissas + array([[0., 0., 1., 1.], + [0., 1., 0., 1.]]) + >>> weights.round(5) + array([0.11111, 0.11111, 0.11111, 0.11111]) + + """ + if segments is None: + order, domain = align_arguments(order, domain) + kwargs = dict(order=order) + else: + order, domain, segments = align_arguments(order, domain, segments) + kwargs = dict(order=order, segments=segments) + + quad_func = partial(ensure_output, quad_func=quad_func) + if auto_scale: + if segments is not None: + quad_func = partial(split_into_segments, quad_func=quad_func) + quad_func = partial(scale_quadrature, quad_func=quad_func) + quad_func = partial(univariate_to_multivariate, quad_func=quad_func) + if isinstance(domain, chaospy.Distribution): + quad_func = partial( + distribution_to_domain, quad_func=quad_func, distribution=domain) + else: + quad_func = partial( + quad_func, lower=numpy.asarray(domain[0]), upper=numpy.asarray(domain[1])) + + return quad_func(**kwargs) + + +def align_arguments( + order, + domain, + segments=None, +): + """ + Extract dimensions from input arguments and broadcast relevant parts. + + Args: + order (int, Sequence[int]): + The quadrature order passed to the quadrature function. + domain (Tuple[float, float], :func:`chaospy.Distribution`): + Either interval on the format ``(lower, upper)`` or a distribution + to integrate over. + segments (Optional[int], Sequence[float]): + The number segments to split the interval on. If sequence is + provided, use as segment edges instead. + + Examples: + >>> order, domain, segments = align_arguments(1, chaospy.Uniform(0, 1), 1) + >>> order, domain, segments + (array([1]), Uniform(), array([1])) + >>> distribution = chaospy.Iid(chaospy.Uniform(0, 1), 2) + >>> order, domain = align_arguments(1, distribution) + >>> order, domain + (array([1, 1]), Iid(Uniform(), 2)) + + """ + args = [numpy.asarray(order)] + if isinstance(domain, chaospy.Distribution): + args += [numpy.zeros(len(domain))] + else: + args += list(domain) + if segments is not None: + segments = numpy.atleast_1d(segments) + assert segments.ndim <= 2 + if segments.ndim == 2: + args.append(segments[:, 0]) + else: + args.append(segments) + args = numpy.broadcast_arrays(*args) + + output = [args.pop(0)] + if not isinstance(domain, chaospy.Distribution): + output += [(args.pop(0), args.pop(0))] + else: + output += [domain] + if segments is not None: + if segments.ndim == 2: + segments = numpy.broadcast_arrays(segments, order)[0].T + else: + segments = args.pop(-1) + output += [segments] + return tuple(output) + + +def ensure_output(quad_func, **kwargs): + """ + Converts arrays to python native types and ensure quadrature output sizes. + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights. + kwargs (Any): + Extra keyword arguments passed to `quad_func`. + + Returns: + Same as ``quad_func(order, **kwargs)`` except numpy elements in + ``kwargs`` is replaced with Python native counterparts and + ``abscissas`` is ensured to be at least 2-dimensional. + + Examples: + >>> def my_quad(order): + ... return (numpy.linspace(0, 1, order+1), + ... 1./numpy.full(order+1, order+2)) + >>> my_quad(2) + (array([0. , 0.5, 1. ]), array([0.25, 0.25, 0.25])) + >>> ensure_output(my_quad, order=numpy.array([2])) + (array([[0. , 0.5, 1. ]]), array([0.25, 0.25, 0.25])) + + """ + kwargs = {key: (value.item() if isinstance(value, numpy.ndarray) else value) + for key, value in kwargs.items()} + abscissas, weights = quad_func(**kwargs) + abscissas = numpy.atleast_2d(abscissas) + assert abscissas.ndim == 2 + assert weights.ndim == 1 + assert abscissas.shape[-1] == len(weights) + return abscissas, weights + + +def univariate_to_multivariate(quad_func, **kwargs): + """ + Turn a univariate quadrature rule into a multivariate rule. + + The one-dimensional quadrature functions are combined into a multivariate + through tensor-product. The dimensionality is inferred from the keyword + arguments. Weights are adjusted to correspond to the multivariate scheme. + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights on the unit + interval. + kwargs (Any): + Keyword arguments passed to `quad_func`. If numerical value is + provided, it is used to infer the dimensions of the multivariate + output. Non-numerical values are passed as is. + + Returns: + Same as ``quad_func(order, **kwargs)`` except with multivariate + supported. + + Examples: + >>> def my_quad(order): + ... return (numpy.linspace(0, 1, order+1)[numpy.newaxis], + ... 1./numpy.full(order+1, order+2)) + >>> my_quad(1) + (array([[0., 1.]]), array([0.33333333, 0.33333333])) + >>> my_quad(2) + (array([[0. , 0.5, 1. ]]), array([0.25, 0.25, 0.25])) + >>> abscissas, weights = univariate_to_multivariate( + ... my_quad, order=numpy.array([1, 2, 1])) + >>> abscissas + array([[0. , 0. , 0. , 0. , 0. , 0. , 1. , 1. , 1. , 1. , 1. , 1. ], + [0. , 0. , 0.5, 0.5, 1. , 1. , 0. , 0. , 0.5, 0.5, 1. , 1. ], + [0. , 1. , 0. , 1. , 0. , 1. , 0. , 1. , 0. , 1. , 0. , 1. ]]) + >>> weights.round(6) + array([0.027778, 0.027778, 0.027778, 0.027778, 0.027778, 0.027778, + 0.027778, 0.027778, 0.027778, 0.027778, 0.027778, 0.027778]) + >>> univariate_to_multivariate(my_quad, order=numpy.array([1])) + (array([[0., 1.]]), array([0.33333333, 0.33333333])) + + """ + sizables = {key: value for key, value in kwargs.items() + if isinstance(value, (int, float, numpy.ndarray))} + sizables["_"] = numpy.zeros(len(kwargs.get("domain", [0]))) + nonsizables = {key: value for key, value in kwargs.items() + if not isinstance(value, (int, float, numpy.ndarray))} + keys = list(sizables) + args = numpy.broadcast_arrays(*[sizables[key] for key in keys]) + assert args[0].ndim == 1 + sizables = {key: value for key, value in zip(keys, args)} + del sizables["_"] + + results = [] + for idx in range(args[0].size): + sizable = kwargs.copy() + sizable.update({key: value[idx].item() + for key, value in sizables.items()}) + abscissas, weights = quad_func(**sizable) + results.append((abscissas.ravel(), weights)) + + abscissas, weights = zip(*results) + return combine_quadrature(abscissas, weights) + + +def distribution_to_domain(quad_func, distribution, **kwargs): + """ + Integrate over a distribution domain. + + Adjust weights to account for probability density. + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights. Must accept + the arguments ``lower`` and ``upper`` to define the interval it is + integrating over. + distribution (:class:`chaospy.Distribution`): + Distribution to adjust quadrature scheme to. + kwargs (Any): + Extra keyword arguments passed to `quad_func`. Can not include the + arguments ``lower`` and ``upper`` as they are taken from + ``distribution``. + + Returns: + Same as ``quad_func(order, **kwargs)`` except arguments ``lower`` and + ``upper`` are now replaced with a new ``distribution`` argument. + + Examples: + >>> def my_quad(lower=0, upper=1): + ... return (numpy.linspace(lower, upper, 5).reshape(1, -1)[:, 1:-1], + ... 1./numpy.full(3, 4)) + >>> my_quad() + (array([[0.25, 0.5 , 0.75]]), array([0.25, 0.25, 0.25])) + >>> distribution_to_domain(my_quad, chaospy.Uniform(-1, 1)) + (array([[-0.5, 0. , 0.5]]), array([0.125, 0.125, 0.125])) + >>> distribution_to_domain(my_quad, chaospy.Beta(2, 2)) + (array([[0.25, 0.5 , 0.75]]), array([0.225, 0.3 , 0.225])) + >>> distribution_to_domain(my_quad, chaospy.Exponential(1)) # doctest: +NORMALIZE_WHITESPACE + (array([[ 8.05924772, 16.11849545, 24.17774317]]), + array([2.32578431e-02, 7.35330570e-06, 2.32485465e-09])) + + """ + assert isinstance(distribution, chaospy.Distribution) + assert "lower" not in kwargs + assert "upper" not in kwargs + lower = distribution.lower + upper = distribution.upper + abscissas, weights = quad_func(lower=lower, upper=upper, **kwargs) + + # Sometimes edge samples (inside the distribution domain) falls out again from simple + # rounding errors. Edge samples needs to be adjusted. + eps = 1e-14*(distribution.upper-distribution.lower) + abscissas_ = numpy.clip(abscissas.T, distribution.lower+eps, distribution.upper-eps).T + weights_ = weights*distribution.pdf(abscissas_).ravel() + weights = weights_*numpy.sum(weights)/(numpy.sum(weights_)*numpy.prod(upper-lower)) + return abscissas, weights + + +def split_into_segments(quad_func, order, segments, **kwargs): + """ + Split a quadrature rule on ta unit interval to multiple segments. + + If both quadrature function includes the abscissas endpoints 0 and 1, then + the endpoints of each subsequent interval is collapsed and their weights + added together. This to avoid nodes from being repeated. + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights on the unit + interval. + order (int): + The quadrature order passed to the quadrature function. + segments (int, Sequence[float]): + The number segments to split the interval on. If sequence is + provided, use as segment edges instead. + kwargs (Any): + Extra keyword arguments passed to `quad_func`. + + Returns: + Same as ``quad_func(order, **kwargs)`` except segmented into + subintervals. + + Examples: + >>> def my_quad(order): + ... return (numpy.linspace(0, 1, order+1)[numpy.newaxis], + ... 1./numpy.full(order+1, order+2)) + >>> my_quad(2) + (array([[0. , 0.5, 1. ]]), array([0.25, 0.25, 0.25])) + >>> split_into_segments(my_quad, 4, segments=2) # doctest: +NORMALIZE_WHITESPACE + (array([[0. , 0.25, 0.5 , 0.75, 1. ]]), + array([0.125, 0.125, 0.25 , 0.125, 0.125])) + >>> split_into_segments(my_quad, 4, segments=3) # doctest: +NORMALIZE_WHITESPACE + (array([[0. , 0.16666667, 0.33333333, 0.66666667, 1. ]]), + array([0.08333333, 0.08333333, 0.19444444, 0.22222222, 0.11111111])) + >>> split_into_segments(my_quad, 4, segments=[0.2, 0.8]) # doctest: +NORMALIZE_WHITESPACE + (array([[0. , 0.1, 0.2, 0.5, 0.8, 0.9, 1. ]]), + array([0.05, 0.05, 0.2 , 0.15, 0.2 , 0.05, 0.05])) + + """ + segments = numpy.array(segments) + if segments.size == 1: + segments = int(segments) + if segments == 1 or order <= 2: + return quad_func(order=order, **kwargs) + if not segments: + segments = int(numpy.sqrt(order)) + assert segments < order, "few samples to distribute than intervals" + nodes = numpy.linspace(0, 1, segments+1) + else: + nodes, segments = numpy.hstack([[0], segments, [1]]), len(segments) + + abscissas = [] + weights = [] + for idx, (lower, upper) in enumerate(zip(nodes[:-1], nodes[1:])): + + assert lower < upper + order_ = order//segments + (idx < (order%segments)) + abscissa, weight = quad_func(order=order_, **kwargs) + weight = weight*(upper-lower) + abscissa = (abscissa.T*(upper-lower)+lower).T + if abscissas and numpy.allclose(abscissas[-1][:, -1], lower): + weights[-1][-1] += weight[0] + abscissa = abscissa[:, 1:] + weight = weight[1:] + abscissas.append(abscissa) + weights.append(weight) + + abscissas = numpy.hstack(abscissas) + weights = numpy.hstack(weights) + assert abscissas.shape == (1, len(weights)) + return abscissas, weights + + +def scale_quadrature(quad_func, order, lower, upper, **kwargs): + """ + Scale quadrature rule designed for unit interval to an arbitrary interval. + + Args: + quad_func (Callable): + Function that creates quadrature abscissas and weights on the unit + interval. + order (int): + The quadrature order passed to the quadrature function. + lower (float): + The new lower limit for the quadrature function. + upper (float): + The new upper limit for the quadrature function. + kwargs (Any): + Extra keyword arguments passed to `quad_func`. + + Returns: + Same as ``quad_func(order, **kwargs)`` except scaled to a new interval. + + + Examples: + >>> def my_quad(order): + ... return (numpy.linspace(0, 1, order+1)[numpy.newaxis], + ... 1./numpy.full(order+1, order+2)) + >>> my_quad(2) + (array([[0. , 0.5, 1. ]]), array([0.25, 0.25, 0.25])) + >>> scale_quadrature(my_quad, 2, lower=0, upper=2) + (array([[0., 1., 2.]]), array([0.5, 0.5, 0.5])) + >>> scale_quadrature(my_quad, 2, lower=-0.5, upper=0.5) + (array([[-0.5, 0. , 0.5]]), array([0.25, 0.25, 0.25])) + + """ + abscissas, weights = quad_func(order=order, **kwargs) + assert numpy.all(abscissas >= 0) and numpy.all(abscissas <= 1) + assert numpy.sum(weights) <= 1+1e-10 + assert numpy.sum(weights > 0) + weights = weights*(upper-lower) + abscissas = (abscissas.T*(upper-lower)+lower).T + return abscissas, weights diff --git a/chaospy/quadrature/jacobi.py b/chaospy/quadrature/jacobi.py new file mode 100644 index 00000000..e9defbca --- /dev/null +++ b/chaospy/quadrature/jacobi.py @@ -0,0 +1,61 @@ +"""Gauss-Jakobi quadrature rule.""" +import numpy +import chaospy + +from .hypercube import hypercube_quadrature + + +def jacobi(order, alpha, beta, lower=-1, upper=1, physicist=False): + """ + Gauss-Jacobi quadrature rule. + + Compute the sample points and weights for Gauss-Jacobi quadrature. The + sample points are the roots of the nth degree Jacobi polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + Gauss-Jacobi physicist means a weight function + :math:`(1-x)^\alpha (1+x)^\beta` and + weights that sum to :math`2^{\alpha+\beta}`, and probabilist means a weight + function is :math:`B(\alpha, \beta) x^{\alpha-1}(1-x)^{\beta-1}` (where + :math:`B` is the beta normalizing constant) which sum to 1. + + Args: + order (int): + The quadrature order. + alpha (float): + First Jakobi shape parameter. + beta (float): + Second Jakobi shape parameter. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.jacobi(3, alpha=2, beta=2) + >>> abscissas + array([[-0.69474659, -0.25056281, 0.25056281, 0.69474659]]) + >>> weights + array([0.09535261, 0.40464739, 0.40464739, 0.09535261]) + + See also: + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Beta(alpha+1, beta+1, lower, upper)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights *= 2**(alpha+beta) if physicist else 1 + return abscissas[numpy.newaxis], weights diff --git a/chaospy/quadrature/gauss_kronrod.py b/chaospy/quadrature/kronrod.py similarity index 92% rename from chaospy/quadrature/gauss_kronrod.py rename to chaospy/quadrature/kronrod.py index 04e47416..ca3b23a8 100644 --- a/chaospy/quadrature/gauss_kronrod.py +++ b/chaospy/quadrature/kronrod.py @@ -7,10 +7,11 @@ Generate Gauss-Kronrod quadrature rules:: >>> distribution = chaospy.Beta(2, 2, lower=-1, upper=1) - >>> for order in range(4): # doctest: +NORMALIZE_WHITESPACE + >>> for order in range(5): # doctest: +NORMALIZE_WHITESPACE ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_kronrod") + ... order, distribution, rule="kronrod") ... print(abscissas.round(2), weights.round(2)) + [[0.]] [1.] [[-0.65 -0. 0.65]] [0.23 0.53 0.23] [[-0.82 -0.45 0. 0.45 0.82]] [0.07 0.26 0.34 0.26 0.07] [[-0.89 -0.65 -0.34 -0. 0.34 0.65 0.89]] @@ -26,7 +27,7 @@ ... abscissas1, weights1 = chaospy.generate_quadrature( ... order, distribution, rule="gaussian") ... abscissas2, weights2 = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_kronrod") + ... order+1, distribution, rule="kronrod") ... print(abscissas1.round(2), abscissas2[:, 1::2].round(2)) [[0.]] [[-0.]] [[-0.58 0.58]] [[-0.58 0.58]] @@ -41,7 +42,7 @@ ... abscissas1, weights1 = chaospy.generate_quadrature( ... order, distribution, rule="gaussian") ... abscissas2, weights2 = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_kronrod") + ... order+1, distribution, rule="kronrod") ... print(abscissas1.round(2), abscissas2.round(2)) [[0.]] [[-1.73 0. 1.73]] [[-1. 1.]] [[-2.45 -1. 0. 1. 2.45]] @@ -50,7 +51,7 @@ exist:: >>> chaospy.generate_quadrature( # doctest: +IGNORE_EXCEPTION_DETAIL - ... 5, distribution, rule="gauss_kronrod") + ... 5, distribution, rule="kronrod") Traceback (most recent call last): ... numpy.linalg.LinAlgError: \ @@ -61,7 +62,7 @@ >>> distribution = chaospy.J( ... chaospy.Uniform(0, 1), chaospy.Beta(4, 5)) >>> abscissas, weights = chaospy.generate_quadrature( - ... 1, distribution, rule="gauss_kronrod") + ... 2, distribution, rule="kronrod") >>> abscissas.round(3) array([[0.037, 0.037, 0.037, 0.037, 0.037, 0.211, 0.211, 0.211, 0.211, 0.211, 0.5 , 0.5 , 0.5 , 0.5 , 0.5 , 0.789, 0.789, 0.789, @@ -80,10 +81,10 @@ import numpy import chaospy -from .combine import combine_quadrature +from .utils import combine_quadrature -def quad_gauss_kronrod( +def kronrod( order, dist, recurrence_algorithm="stieltjes", @@ -158,14 +159,13 @@ def quad_gauss_kronrod( Example: >>> distribution = chaospy.Uniform(-1, 1) - >>> abscissas, weights = quad_gauss_kronrod(3, distribution) + >>> abscissas, weights = chaospy.quadrature.kronrod(4, distribution) >>> abscissas.round(2) array([[-0.98, -0.86, -0.64, -0.34, 0. , 0.34, 0.64, 0.86, 0.98]]) >>> weights.round(3) array([0.031, 0.085, 0.133, 0.163, 0.173, 0.163, 0.133, 0.085, 0.031]) - """ - assert not rule.startswith("gauss"), "recursive Gaussian quadrature call" + """ length = int(numpy.ceil(3*(order+1) / 2.0)) coefficients = chaospy.construct_recurrence_coefficients( order=length, @@ -177,7 +177,7 @@ def quad_gauss_kronrod( n_max=n_max, ) - coefficients = [kronrod_jacobi(order+1, coeffs) for coeffs in coefficients] + coefficients = [kronrod_jacobi(order, coeffs) for coeffs in coefficients] abscissas, weights = chaospy.coefficients_to_quadrature(coefficients) return combine_quadrature(abscissas, weights) @@ -202,7 +202,8 @@ def kronrod_jacobi(order, coeffs): Three terms recurrence coefficients of the Gauss-Kronrod quadrature rule. """ - assert len(coeffs[0]) == int(math.ceil(3*order/2.0))+1 + if not order: + return kronrod_jacobi(1, coeffs)[:, :1] bound = int(math.floor(3*order/2.0))+1 coeffs_a = numpy.zeros(2*order+1) diff --git a/chaospy/quadrature/laguerre.py b/chaospy/quadrature/laguerre.py new file mode 100644 index 00000000..fab351b2 --- /dev/null +++ b/chaospy/quadrature/laguerre.py @@ -0,0 +1,54 @@ +"""Generalized Gauss-Laguerre quadrature rule.""" +import numpy +from scipy.special import gamma +import chaospy + +from .hypercube import hypercube_quadrature + + +def laguerre(order, alpha=0., physicist=False): + r""" + Generalized Gauss-Laguerre quadrature rule. + + Compute the sample points and weights for Gauss-Laguerre quadrature. The + sample points are the roots of the nth degree Laguerre polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less. + + Gaussian quadrature come in two variants: physicist and probabilist. For + Gauss-Laguerre physicist means a weight function :math:`x^\alpha e^{-x}` + and weights that sum to :math`\Gamma(\alpha+1)`, and probabilist means a + weight function is :math:`x^\alpha e^{-x}` and sum to 1. + + Args: + order (int): + The quadrature order. + alpha (float): + Shape parameter. Defaults to non-generalized Laguerre if 0. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.laguerre(2) + >>> abscissas + array([[0.41577456, 2.29428036, 6.28994508]]) + >>> weights + array([0.71109301, 0.27851773, 0.01038926]) + + See also: + :func:`chaospy.quadrature.gaussian` + + """ + order = int(order) + coefficients = chaospy.construct_recurrence_coefficients( + order=order, dist=chaospy.Gamma(alpha+1)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + weights *= gamma(alpha+1) if physicist else 1 + return abscissas[numpy.newaxis], weights diff --git a/chaospy/quadrature/legendre.py b/chaospy/quadrature/legendre.py new file mode 100644 index 00000000..33a9c294 --- /dev/null +++ b/chaospy/quadrature/legendre.py @@ -0,0 +1,163 @@ +"""Gauss-Legendre quadrature rule.""" +try: + from functools import lru_cache +except ImportError: # pragma: no cover + from functools32 import lru_cache + +import numpy +import chaospy + +from .hypercube import hypercube_quadrature + + +def legendre(order, lower=-1., upper=1., physicist=False): + """ + Gauss-Legendre quadrature rule. + + Compute the sample points and weights for Gauss-Legendre quadrature. The + sample points are the roots of the N-th degree Legendre polynomial. These + sample points and weights correctly integrate polynomials of degree + :math:`2N-1` or less over the interval ``[lower, upper]``. + + Gaussian quadrature come in two variants: physicist and probabilist. For + Gauss-Legendre physicist means a weight function constant 1 and weights sum + to ``upper-lower``, and probabilist means weight function constant + ``1/(upper-lower)`` while weights sum to 1. + + Args: + order (int): + The quadrature order. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + abscissas (numpy.ndarray): + The ``order+1`` quadrature points for where to evaluate the model + function with. + weights (numpy.ndarray): + The quadrature weights associated with each abscissas. + + Examples: + >>> abscissas, weights = chaospy.quadrature.legendre(2) + >>> abscissas + array([[-0.77459667, 0. , 0.77459667]]) + >>> weights + array([0.27777778, 0.44444444, 0.27777778]) + + See also: + :func:`chaospy.quadrature.gaussian` + :func:`chaospy.quadrature.legendre_proxy` + + """ + abscissas, weights = hypercube_quadrature( + legendre_simple, + order=int(order), + domain=(float(lower), float(upper)), + ) + weights = weights if physicist else weights/(upper-lower) + return abscissas, weights + + +def legendre_proxy( + order, + domain=(0, 1), + segments=1, +): + r""" + Generate proxy abscissas and weights from Legendre quadrature. + + Legendre provides optimal abscissas :math:`X_i` and weights :math:`W_i` to + solve the integration problem: + + .. math:: + + \int f(x) dx \approx \sum W_i f(X_i) + + over a function :math:`f`, where the probability density function :math:`p` + is uniform. + + Since the weight function is constant, it can in principle be used to + integrate any density function by considering it a part of the function. In + other words: + + .. math:: + + \int p(x) f(x) \approx \sum W_i p(X_i) f(X_i) = \sum W_i^' f(X_i) + + So when providing non-uniform distribution as `domain`, the weights will be + adjusted with: + + .. math:: + + W_i^' = W_i p(X_i) + + Bounds of the Legendre schemes is chosen to be the same as the distribution + provided. This makes it a bad choice for unbound distributions. + + To get optimal abscissas and weights directly from a density, use + :func:`chaospy.quadrature.gaussian` instead. + + Args: + order (int, numpy.ndarray): + Quadrature order. + domain (:class:`chaospy.Distribution`, numpy.ndarray): + Either distribution or bounding of interval to integrate over. + segments (int): + Split intervals into steps subintervals and create a patched + quadrature based on the segmented quadrature. Can not be lower than + `order`. If 0 is provided, default to square root of `order`. + Nested samples only appear when the number of segments are fixed. + + Returns: + abscissas (numpy.ndarray): + The quadrature points for where to evaluate the model function with + ``abscissas.shape == (len(dist), N)`` where ``N`` is the number of + samples. + weights (numpy.ndarray): + The quadrature weights with ``weights.shape == (N,)``. + + Example: + >>> abscissas, weights = chaospy.quadrature.legendre_proxy(3) + >>> abscissas.round(4) + array([[0.0694, 0.33 , 0.67 , 0.9306]]) + >>> weights.round(4) + array([0.1739, 0.3261, 0.3261, 0.1739]) + >>> abscissas, weights = chaospy.quadrature.legendre_proxy(3, chaospy.Uniform(0, 1)) + >>> abscissas.round(4) + array([[0.0694, 0.33 , 0.67 , 0.9306]]) + >>> weights.round(4) + array([0.1739, 0.3261, 0.3261, 0.1739]) + >>> abscissas, weights = chaospy.quadrature.legendre_proxy(3, chaospy.Beta(2, 2)) + >>> abscissas.round(4) + array([[0.0694, 0.33 , 0.67 , 0.9306]]) + >>> weights.round(4) + array([0.0674, 0.4326, 0.4326, 0.0674]) + + See also: + :func:`chaospy.quadrature.legendre`, + :func:`chaospy.quadrature.gaussian` + + """ + return hypercube_quadrature( + legendre_simple, + order=order, + domain=domain, + segments=segments, + ) + + +@lru_cache(None) +def legendre_simple(order): + """ + Simple Legendre quadrature on the [0, 1] interval. + + Use :func:`chaospy.quadrature.legendre` instead. + """ + coefficients = chaospy.construct_recurrence_coefficients( + order=int(order), dist=chaospy.Uniform(-1, 1)) + [abscissas], [weights] = chaospy.coefficients_to_quadrature(coefficients) + return abscissas*0.5+0.5, weights diff --git a/chaospy/quadrature/leja.py b/chaospy/quadrature/leja.py index 9fedbe46..256d71a1 100644 --- a/chaospy/quadrature/leja.py +++ b/chaospy/quadrature/leja.py @@ -25,13 +25,13 @@ from scipy.optimize import fminbound import chaospy -from .combine import combine +from .utils import combine -def quad_leja( +def leja( order, dist, - rule="fejer", + rule="fejer_2", ): """ Generate Leja quadrature node. @@ -59,8 +59,8 @@ def quad_leja( :cite:`narayan_adaptive_2014`. Example: - >>> abscissas, weights = quad_leja( - ... 2, chaospy.Iid(chaospy.Normal(0, 1), 2)) + >>> distribution = chaospy.Iid(chaospy.Normal(0, 1), 2) + >>> abscissas, weights = chaospy.quadrature.leja(2, distribution) >>> abscissas.round(2) array([[-1.41, -1.41, -1.41, 0. , 0. , 0. , 1.76, 1.76, 1.76], [-1.41, 0. , 1.76, -1.41, 0. , 1.76, -1.41, 0. , 1.76]]) @@ -73,7 +73,7 @@ def quad_leja( raise chaospy.StochasticallyDependentError( "Leja quadrature do not supper distribution with dependencies.") order = numpy.broadcast_to(order, len(dist)) - out = [quad_leja(order[_], dist[_]) for _ in range(len(dist))] + out = [leja(order[_], dist[_]) for _ in range(len(dist))] abscissas = [_[0][0] for _ in out] weights = [_[1] for _ in out] abscissas = combine(abscissas).T diff --git a/chaospy/quadrature/gauss_lobatto.py b/chaospy/quadrature/lobatto.py similarity index 84% rename from chaospy/quadrature/gauss_lobatto.py rename to chaospy/quadrature/lobatto.py index b5ab9969..fab80369 100644 --- a/chaospy/quadrature/gauss_lobatto.py +++ b/chaospy/quadrature/lobatto.py @@ -1,13 +1,6 @@ +# -*- coding: utf-8 -*- """ -Gauss-Radau formula for numerical estimation of integrals. It requires -:math:`m+1` points and fits all Polynomials to degree :math:`2m`, so it -effectively fits exactly all Polynomials of degree :math:`2m-3`. - -Gauss-Radau is defined by having two abscissas to be fixed to the endpoints, -while the others are built around these points. So if a distribution is defined -on the interval ``(a, b)``, then both ``a`` and ``b`` are abscissas in this -scheme. Note though that this does not always possible to achieve in practice, -and an error might be raised. +Generate the abscissas and weights in Gauss-Loboto quadrature. Example usage ------------- @@ -17,7 +10,7 @@ >>> distribution = chaospy.Beta(2, 2, lower=-1, upper=1) >>> for order in range(4): # doctest: +NORMALIZE_WHITESPACE ... X, W = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_lobatto") + ... order, distribution, rule="lobatto") ... print(X.round(2), W.round(2)) [[-1.]] [1.] [[-1. 1.]] [0.5 0.5] @@ -29,7 +22,7 @@ >>> distribution = chaospy.J(chaospy.Uniform(0, 1), chaospy.Beta(4, 5)) >>> X, W = chaospy.generate_quadrature( - ... 2, distribution, rule="gauss_lobatto") + ... 2, distribution, rule="lobatto") >>> X.round(3) array([[-0. , -0. , -0. , -0. , 0.276, 0.276, 0.276, 0.276, 0.724, 0.724, 0.724, 0.724, 1. , 1. , 1. , 1. ], @@ -43,14 +36,14 @@ from scipy.linalg import solve_banded, solve import chaospy -from .combine import combine_quadrature +from .utils import combine_quadrature -def quad_gauss_lobatto( +def lobatto( order, dist, recurrence_algorithm="stieltjes", - rule="fejer", + rule="fejer_2", tolerance=1e-10, scaling=3, n_max=5000, @@ -58,6 +51,15 @@ def quad_gauss_lobatto( """ Generate the abscissas and weights in Gauss-Loboto quadrature. + Also known as Lobatto quadrature, named after Dutch mathematician Rehuel + Lobatto. It is similar to Gaussian quadrature with the following + differences: + + * The integration points include the end points of the integration + interval. + * It is accurate for polynomials up to degree :math:`2n–3`, where :math:`n` + is the number of integration points. + Args: order (int): Quadrature order. @@ -89,12 +91,13 @@ def quad_gauss_lobatto( The quadrature weights with ``weights.shape == (N,)``. Example: - >>> abscissas, weights = quad_gauss_lobatto( - ... 4, chaospy.Uniform(-1, 1)) + >>> distribution = chaospy.Uniform(-1, 1) + >>> abscissas, weights = chaospy.quadrature.lobatto(4, distribution) >>> abscissas.round(3) array([[-1. , -0.872, -0.592, -0.209, 0.209, 0.592, 0.872, 1. ]]) >>> weights.round(3) array([0.018, 0.105, 0.171, 0.206, 0.206, 0.171, 0.105, 0.018]) + """ assert not rule.startswith("gauss"), "recursive Gaussian quadrature call" if order == 0: diff --git a/chaospy/quadrature/newton_cotes.py b/chaospy/quadrature/newton_cotes.py index 06f49234..e93b5c15 100644 --- a/chaospy/quadrature/newton_cotes.py +++ b/chaospy/quadrature/newton_cotes.py @@ -1,6 +1,4 @@ """ -Newton-Cotes quadrature, are a group of formulas for numerical integration -based on evaluating the integrand at equally spaced points. Example usage ------------- @@ -51,19 +49,22 @@ except ImportError: # pragma: no cover from functools32 import lru_cache import numpy -from scipy.integrate import newton_cotes +from scipy import integrate -from .combine import combine_quadrature +from .hypercube import hypercube_quadrature -def quad_newton_cotes(order, domain=(0, 1), growth=False, segments=1): +def newton_cotes(order, domain=(0, 1), growth=False, segments=1): """ Generate the abscissas and weights in Newton-Cotes quadrature. + Newton-Cotes quadrature, are a group of formulas for numerical integration + based on evaluating the integrand at equally spaced points. + Args: - order (int, numpy.ndarray): + order (int, numpy.ndarray:): Quadrature order. - domain (chaospy.distributions.baseclass.Distribution, numpy.ndarray): + domain (:func:`chaospy.Distribution`, ;class:`numpy.ndarray`): Either distribution or bounding of interval to integrate over. growth (bool): If True sets the growth rule for the quadrature rule to only @@ -84,75 +85,31 @@ def quad_newton_cotes(order, domain=(0, 1), growth=False, segments=1): The quadrature weights with ``weights.shape == (N,)``. Examples: - >>> abscissas, weights = quad_newton_cotes(4) + >>> abscissas, weights = chaospy.quadrature.newton_cotes(4) >>> abscissas.round(4) array([[0. , 0.25, 0.5 , 0.75, 1. ]]) >>> weights.round(4) array([0.0778, 0.3556, 0.1333, 0.3556, 0.0778]) - >>> abscissas, weights = quad_newton_cotes(4, segments=2) + >>> abscissas, weights = chaospy.quadrature.newton_cotes(4, segments=2) >>> abscissas.round(4) array([[0. , 0.25, 0.5 , 0.75, 1. ]]) >>> weights.round(4) - array([0.1667, 0.6667, 0.3333, 0.6667, 0.1667]) + array([0.0833, 0.3333, 0.1667, 0.3333, 0.0833]) + """ - from ..distributions.baseclass import Distribution - if isinstance(domain, Distribution): - abscissas, weights = quad_newton_cotes( - order, (domain.lower, domain.upper), growth, segments) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - order = numpy.asarray(order, dtype=int).flatten() - lower, upper = domain - lower = numpy.asarray(lower).flatten() - upper = numpy.asarray(upper).flatten() - dim = max(lower.size, upper.size, order.size) - - order = order*numpy.ones(dim, dtype=int) - lower = lower*numpy.ones(dim) - upper = upper*numpy.ones(dim) - segments = segments*numpy.ones(dim, dtype=int) - - results = [_newton_cotes(*args, growth=growth) - for args in zip(order, lower, upper, segments)] - abscissas = [args[0] for args in results] - weights = [args[1] for args in results] - return combine_quadrature(abscissas, weights) + order = numpy.asarray(order) + order = numpy.where(growth, numpy.where(order, 2**order, 0), order) + return hypercube_quadrature( + _newton_cotes, + order=order, + domain=domain, + segments=segments, + ) @lru_cache(None) -def _newton_cotes(order, lower, upper, segments=1, growth=False): +def _newton_cotes(order): """Backend for Newton-Cotes quadrature rule.""" if order == 0: - return numpy.array([0.5*(lower+upper)]), numpy.ones(1) - order = 2**order if growth else order - - if segments != 1 and order > 2: - segments = segments if segments else int(numpy.sqrt(order)) - assert segments < order, "fewer samples to distribute than intervals" - abscissas = [] - weights = [] - - nodes = numpy.linspace(0, 1, segments+1) - for lower, upper in zip(nodes[:-1], nodes[1:]): - - abscissa, weight = _newton_cotes(order//segments, lower, upper) - weight = weight*(upper-lower) - if abscissas: - weights[-1] += weight[0] - abscissa = abscissa[1:] - weight = weight[1:] - abscissas.extend(abscissa) - weights.extend(weight) - - assert len(abscissas) == order+1 - assert len(weights) == order+1 - return numpy.array(abscissas), numpy.array(weights) - - return ( - numpy.linspace(lower, upper, order+1), - newton_cotes(order)[0]/(upper-lower)/order, - ) + return numpy.full((1, 1), 0.5), numpy.ones(1) + return numpy.linspace(0, 1, order+1), integrate.newton_cotes(order)[0]/order diff --git a/chaospy/quadrature/gauss_patterson.py b/chaospy/quadrature/patterson.py similarity index 97% rename from chaospy/quadrature/gauss_patterson.py rename to chaospy/quadrature/patterson.py index 8b3fd38e..e0c17892 100644 --- a/chaospy/quadrature/gauss_patterson.py +++ b/chaospy/quadrature/patterson.py @@ -7,22 +7,26 @@ With increasing order:: >>> distribution = chaospy.Beta(2, 2, lower=-1, upper=1) - >>> for order in range(3): # doctest: +NORMALIZE_WHITESPACE + >>> for order in range(3): ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_patterson") + ... order, distribution, rule="patterson") ... print(abscissas.round(2), weights.round(2)) [[0.]] [1.] [[-0.77 0. 0.77]] [0.17 0.67 0.17] - [[-0.96 -0.77 -0.43 0. 0.43 0.77 0.96]] - [0.01 0.08 0.24 0.34 0.24 0.08 0.01] + [[-0.96 -0.77 -0.43 0. 0.43 0.77 0.96]] [0.01 0.08 0.24 0.34 0.24 0.08 0.01] """ +try: + from functools import lru_cache +except ImportError: # pragma: no cover + from functools32 import lru_cache + import numpy import chaospy -from .combine import combine +from .hypercube import hypercube_quadrature -def quad_gauss_patterson(order, domain): +def patterson(order, domain): """ Generate Gauss-Patterson quadrature abscissa and weights. @@ -32,7 +36,7 @@ def quad_gauss_patterson(order, domain): Gauss-Patterson rules do not have the super-high precision of the Gauss-Legendre rules. They trade this precision in exchange for the advantages of nestedness. This means that Gauss-Patterson rules are only - available for orders of 1, 3, 7, 15, 31, 63, 127, 255 or 511. + available for orders of 0, 2, 6, 14, 30, 62, 126, 254 or 510. Args: order (int): @@ -55,8 +59,8 @@ def quad_gauss_patterson(order, domain): Points to Quadrature Formulae":cite:`patterson_optimum_1968`. Example: - >>> abscissas, weights = chaospy.quad_gauss_patterson( - ... 1, chaospy.Iid(chaospy.Uniform(0, 1), 2)) + >>> distribution = chaospy.Iid(chaospy.Uniform(0, 1), 2) + >>> abscissas, weights = chaospy.quadrature.patterson(1, distribution) >>> abscissas.round(3) array([[0.113, 0.113, 0.113, 0.5 , 0.5 , 0.5 , 0.887, 0.887, 0.887], [0.113, 0.5 , 0.887, 0.113, 0.5 , 0.887, 0.113, 0.5 , 0.887]]) @@ -64,48 +68,30 @@ def quad_gauss_patterson(order, domain): array([0.077, 0.123, 0.077, 0.123, 0.198, 0.123, 0.077, 0.123, 0.077]) """ - if isinstance(domain, chaospy.Distribution): - abscissas, weights = quad_gauss_patterson( - order, (domain.lower, domain.upper)) - eps = 1e-14*(domain.upper-domain.lower) - abscissas_ = numpy.clip(abscissas.T, domain.lower+eps, domain.upper-eps).T - weights *= domain.pdf(abscissas_).flatten() - weights /= numpy.sum(weights) - return abscissas, weights - - domain = numpy.array(numpy.broadcast_arrays(*domain)) - lower, upper = numpy.atleast_2d(domain.T).T - lower, upper, order = numpy.broadcast_arrays(lower, upper, order) - - if lower.size > 1: - values = [quad_gauss_patterson(order_, (lower_, upper_)) - for order_, lower_, upper_ in zip(order, lower, upper)] + abscissas, weights = hypercube_quadrature( + patterson_simple, + order=int(order), + domain=domain, + ) + return abscissas, weights - abscissas = [value[0][0] for value in values] - weights = [value[1] for value in values] - abscissas = combine(abscissas).T - weights = numpy.prod(combine(weights), -1) - return abscissas, weights - order = sorted(PATTERSON_VALUES.keys())[int(order)] +@lru_cache(None) +def patterson_simple(order): + assert order < 9 abscissas, weights = PATTERSON_VALUES[int(order)] - - abscissas = .5*(abscissas*(upper-lower)+upper+lower) - abscissas = abscissas.reshape(1, abscissas.size) - weights /= numpy.sum(weights) - - return abscissas, weights + return numpy.array(abscissas)*0.5+0.5, numpy.array(weights)/2. PATTERSON_VALUES = { 0 : ((0e+00,), (2.0e+00,)), - 3 : (( + 1 : (( -0.77459666924148337704e+00, 0.0e+00, 0.77459666924148337704e+00, ), ( 0.555555555555555555556e+00, 0.888888888888888888889e+00, 0.555555555555555555556e+00, )), - 7 : (( + 2 : (( -0.96049126870802028342e+00, -0.77459666924148337704e+00, -0.43424374934680255800e+00, 0.0e+00, 0.43424374934680255800e+00, 0.77459666924148337704e+00, 0.96049126870802028342e+00, @@ -115,7 +101,7 @@ def quad_gauss_patterson(order, domain): 0.401397414775962222905e+00, 0.268488089868333440729e+00, 0.104656226026467265194e+00, )), - 15 : (( + 3 : (( -0.99383196321275502221e+00, -0.96049126870802028342e+00, -0.88845923287225699889e+00, -0.77459666924148337704e+00, -0.62110294673722640294e+00, -0.43424374934680255800e+00, @@ -133,7 +119,7 @@ def quad_gauss_patterson(order, domain): 0.0929271953151245376859e+00, 0.0516032829970797396969e+00, 0.0170017196299402603390e+00, )), - 31 : (( + 4 : (( -0.99909812496766759766e+00, -0.99383196321275502221e+00, -0.98153114955374010687e+00, -0.96049126870802028342e+00, -0.92965485742974005667e+00, -0.88845923287225699889e+00, @@ -167,7 +153,7 @@ def quad_gauss_patterson(order, domain): 0.0164460498543878109338e+00, 0.00843456573932110624631e+00, 0.00254478079156187441540e+00, )), - 63 : (( + 5 : (( -0.99987288812035761194e+00, -0.99909812496766759766e+00, -0.99720625937222195908e+00, -0.99383196321275502221e+00, -0.98868475754742947994e+00, -0.98153114955374010687e+00, @@ -233,7 +219,7 @@ def quad_gauss_patterson(order, domain): 0.00257904979468568827243e+00, 0.00126515655623006801137e+00, 0.000363221481845530659694e+00, )), - 127 : (( + 6 : (( -0.99998243035489159858e+00, -0.99987288812035761194e+00, -0.99959879967191068325e+00, -0.99909812496766759766e+00, -0.99831663531840739253e+00, -0.99720625937222195908e+00, @@ -363,7 +349,7 @@ def quad_gauss_patterson(order, domain): 0.000377746646326984660274e+00, 0.000180739564445388357820e+00, 0.0000505360952078625176247e+00, )), - 255 : (( + 7 : (( -0.99999759637974846462e+00, -0.99998243035489159858e+00, -0.99994399620705437576e+00, -0.99987288812035761194e+00, -0.99976049092443204733e+00, -0.99959879967191068325e+00, @@ -621,7 +607,7 @@ def quad_gauss_patterson(order, domain): 0.53275293669780613125e-04, 0.25157870384280661489e-04, 0.69379364324108267170e-05, )), - 511 : (( + 8 : (( -0.999999672956734384381e+00, -0.999997596379748464620e+00, -0.999992298136257588028e+00, -0.999982430354891598580e+00, -0.999966730098486276883e+00, -0.999943996207054375764e+00, diff --git a/chaospy/quadrature/gauss_radau.py b/chaospy/quadrature/radau.py similarity index 78% rename from chaospy/quadrature/gauss_radau.py rename to chaospy/quadrature/radau.py index c74e8d2a..c415b44f 100644 --- a/chaospy/quadrature/gauss_radau.py +++ b/chaospy/quadrature/radau.py @@ -1,14 +1,5 @@ """ -Gauss-Radau formula for numerical estimation of integrals. It requires -:math:`m+1` points and fits all Polynomials to degree :math:`2m`, so it -effectively fits exactly all Polynomials of degree :math:`2m-1`. - -It allows for a single abscissas to be user defined, while the others are built -around this point. - -Canonically, Radau is built around Legendre weight function with the fixed -point at the left end. Not all distributions/fixed point combinations allows -for the building of a quadrature scheme. +Generate the quadrature nodes and weights in Gauss-Radau quadrature. Example usage ------------- @@ -18,7 +9,7 @@ >>> distribution = chaospy.Beta(2, 2, lower=-1, upper=1) >>> for order in range(4): # doctest: +NORMALIZE_WHITESPACE ... abscissas, weights = chaospy.generate_quadrature( - ... order, distribution, rule="gauss_radau") + ... order, distribution, rule="radau") ... print(abscissas.round(2), weights.round(2)) [[-1.]] [1.] [[-1. 0.2]] [0.17 0.83] @@ -31,7 +22,7 @@ >>> distribution = chaospy.J( ... chaospy.Uniform(0, 1), chaospy.Beta(4, 5)) >>> abscissas, weights = chaospy.generate_quadrature( - ... 1, distribution, rule="gauss_radau") + ... 1, distribution, rule="radau") >>> abscissas.round(3) array([[0. , 0. , 0.667, 0.667], [0. , 0.5 , 0. , 0.5 ]]) @@ -42,7 +33,7 @@ >>> distribution = chaospy.Uniform(lower=-1, upper=1) >>> for fixed_point in numpy.linspace(-1, 1, 6): - ... abscissas, weights = chaospy.quad_gauss_radau( + ... abscissas, weights = chaospy.quadrature.radau( ... 2, distribution, fixed_point) ... print(abscissas.round(2), weights.round(2)) [[-1. -0.58 0.18 0.82]] [0.06 0.33 0.39 0.22] @@ -54,7 +45,7 @@ However, a fixed point at 0 is not allowed:: - >>> chaospy.quad_gauss_radau( # doctest: +IGNORE_EXCEPTION_DETAIL + >>> chaospy.quadrature.radau( # doctest: +IGNORE_EXCEPTION_DETAIL ... 3, distribution, fixed_point=0) Traceback (most recent call last): ... @@ -64,10 +55,10 @@ import scipy.linalg import chaospy -from .combine import combine_quadrature +from .utils import combine_quadrature -def quad_gauss_radau( +def radau( order, dist, fixed_point=None, @@ -80,15 +71,26 @@ def quad_gauss_radau( """ Generate the quadrature nodes and weights in Gauss-Radau quadrature. + Gauss-Radau formula for numerical estimation of integrals. It requires + :math:`m+1` points and fits all Polynomials to degree :math:`2m`, so it + effectively fits exactly all Polynomials of degree :math:`2m-1`. + + It allows for a single abscissas to be user defined, while the others are + built around this point. + + Canonically, Radau is built around Legendre weight function with the fixed + point at the left end. Not all distributions/fixed point combinations + allows for the building of a quadrature scheme. + Args: order (int): Quadrature order. - dist (chaospy.distributions.baseclass.Distribution): + dist (:class:`chaospy.Distribution`): The distribution weights to be used to create higher order nodes from. fixed_point (float): Fixed point abscissas assumed to be included in the quadrature. If - imitted, use distribution lower point ``dist.range()[0]``. + omitted, use distribution lower bound. rule (str): In the case of ``lanczos`` or ``stieltjes``, defines the proxy-integration scheme. @@ -101,22 +103,22 @@ def quad_gauss_radau( if that fails. Returns: - (numpy.ndarray, numpy.ndarray): - abscissas: - The quadrature points for where to evaluate the model function - with ``abscissas.shape == (len(dist), N)`` where ``N`` is the - number of samples. - weights: - The quadrature weights with ``weights.shape == (N,)``. + abscissas (numpy.ndarray): + The quadrature points for where to evaluate the model function + with ``abscissas.shape == (len(dist), N)`` where ``N`` is the + number of samples. + weights (numpy.ndarray): + The quadrature weights with ``weights.shape == (N,)``. Example: - >>> abscissas, weights = quad_gauss_radau(4, chaospy.Uniform(-1, 1)) + >>> distribution = chaospy.Uniform(-1, 1) + >>> abscissas, weights = chaospy.quadrature.radau(4, distribution) >>> abscissas.round(3) array([[-1. , -0.887, -0.64 , -0.295, 0.094, 0.468, 0.771, 0.955]]) >>> weights.round(3) array([0.016, 0.093, 0.152, 0.188, 0.196, 0.174, 0.125, 0.057]) + """ - assert not rule.startswith("gauss"), "recursive Gaussian quadrature call" if fixed_point is None: fixed_point = dist.lower else: diff --git a/chaospy/quadrature/sparse_grid.py b/chaospy/quadrature/sparse_grid.py index 45e2c405..89b619ee 100644 --- a/chaospy/quadrature/sparse_grid.py +++ b/chaospy/quadrature/sparse_grid.py @@ -9,7 +9,7 @@ import chaospy -def construct_sparse_grid( +def sparse_grid( order, dist, growth=None, @@ -55,15 +55,14 @@ def construct_sparse_grid( that ``abscissas.shape == (len(dist), len(weights))``. Example: - >>> distribution = chaospy.J( - ... chaospy.Normal(0, 1), chaospy.Uniform(-1, 1)) - >>> abscissas, weights = construct_sparse_grid(1, distribution) + >>> distribution = chaospy.J(chaospy.Normal(0, 1), chaospy.Uniform(-1, 1)) + >>> abscissas, weights = chaospy.quadrature.sparse_grid(1, distribution) >>> abscissas.round(4) array([[-1. , 0. , 0. , 0. , 1. ], [ 0. , -0.5774, 0. , 0.5774, 0. ]]) >>> weights.round(4) array([ 0.5, 0.5, -1. , 0.5, 0.5]) - >>> abscissas, weights = construct_sparse_grid([2, 1], distribution) + >>> abscissas, weights = chaospy.quadrature.sparse_grid([2, 1], distribution) >>> abscissas.round(2) array([[-1.73, -1. , -1. , -1. , 0. , 1. , 1. , 1. , 1.73], [ 0. , -0.58, 0. , 0.58, 0. , -0.58, 0. , 0.58, 0. ]]) diff --git a/chaospy/quadrature/combine.py b/chaospy/quadrature/utils.py similarity index 92% rename from chaospy/quadrature/combine.py rename to chaospy/quadrature/utils.py index 71acb8a9..8a563851 100644 --- a/chaospy/quadrature/combine.py +++ b/chaospy/quadrature/utils.py @@ -1,5 +1,7 @@ -"""Function to combine two dataset together with a tensor product.""" +from __future__ import division +from functools import partial import numpy +import chaospy def combine(args): @@ -20,7 +22,7 @@ def combine(args): Examples: >>> A, B = [1,2], [[4,4],[5,6]] - >>> chaospy.quadrature.combine([A, B]) + >>> combine([A, B]) array([[1, 4, 4], [1, 5, 6], [2, 4, 4], @@ -74,6 +76,6 @@ def combine_quadrature( weights = numpy.prod(weights, -1) assert len(weights.shape) == 1 - assert abscissas.shape == (dim,) + weights.shape + assert abscissas.shape == (dim,) + weights.shape, (abscissas, weights) return abscissas, weights diff --git a/docs/reference/quadrature.rst b/docs/reference/quadrature.rst index 52e57fd3..7e92008d 100644 --- a/docs/reference/quadrature.rst +++ b/docs/reference/quadrature.rst @@ -1,42 +1,64 @@ .. _quadrature_collection: -Quadrature Integration +Quadrature integration ====================== -.. currentmodule:: chaospy +.. currentmodule:: chaospy.quadrature -Gaussian Quadrature +Standard library +---------------- + +.. autosummary:: + :toctree: api + + clenshaw_curtis + fejer_1 + fejer_2 + gaussian + grid + legendre_proxy + leja + newton_cotes + +Discrete densities +------------------ + +.. autosummary:: + :toctree: api + + discrete + +Gaussian extensions ------------------- .. autosummary:: :toctree: api - quad_gaussian - quad_gauss_legendre - quad_gauss_kronrod - quad_gauss_lobatto - quad_gauss_patterson - quad_gauss_radau + legendre + kronrod + lobatto + patterson + radau -Other Quadrature Rules ----------------------- +Gaussian predefined +------------------- .. autosummary:: :toctree: api - quad_fejer - quad_clenshaw_curtis - quad_discrete - quad_genz_keister - quad_leja - quad_newton_cotes + chebyshev_1 + chebyshev_2 + gegenbauer + hermite + jacobi + legendre + laguerre -Utility Functions ------------------ +Helper functions +---------------- .. autosummary:: :toctree: api + sparse_grid kronrod_jacobi - combine - construct_sparse_grid diff --git a/docs/user_guide/quadrature.rst b/docs/user_guide/quadrature.rst index fa6a1042..edc02679 100644 --- a/docs/user_guide/quadrature.rst +++ b/docs/user_guide/quadrature.rst @@ -1,8 +1,32 @@ .. _quadrature: -Quadrature Integration +Quadrature integration ====================== +To see see the various quadrature functions available, see +:ref:`quadrature_collection` for an overview. Alternative to these individual +function, it is possible to use the convenience function +:func:`chaospy.generate_quadrature` which works as a frontend wrapper for all +quadrature functions. Just pass the name of the quadrature rule as a flag. +For example, to make Clenshaw-Curtis quadratures one could either do: + +.. code:: python + + >>> distribution = chaospy.Uniform(-1, 1) + >>> chaospy.quadrature.clenshaw_curtis(2, distribution) + (array([[-1., 0., 1.]]), array([0.16666667, 0.66666667, 0.16666667])) + +Or through the convenience function: + +.. code:: python + + >>> chaospy.generate_quadrature(2, distribution, rule="clenshaw_curtis") + (array([[-1., 0., 1.]]), array([0.16666667, 0.66666667, 0.16666667])) + +In other words the name of a quadrature function ``chaospy.quadrature.`` +can also be used as a keyword argument ``rule=""`` in +:func:`chaospy.generate_quadrature`. + Introduction ------------ @@ -12,121 +36,45 @@ requiring an analytical definition. In the scope of ``chaospy`` we limit this scope to focus on methods that can be reduced to the following approximation: .. math:: + \int p(x) g(x) dx \approx \sum_{n=1}^N W_n g(X_n) Here :math:`p(x)` is an weight function, which is assumed to be an probability -distribution, and :math:`W_n` and :math:`X_n` are respectively quadrature +density function, and :math:`W_n` and :math:`X_n` are respectively quadrature weights and abscissas used to define the approximation. -The simplest example applying such an approximation is Monte Carlo integration. -In such a method, you only need to select :math:`W_n=1/N` and :math:`X_n` to be -independent identical distributed samples drawn from the distribution of -:math:`p(x)`. In practice:: +The simplest application of such an approximation is Monte Carlo integration. +In Monte Carlo you only need to select :math:`W_n=1/N` for all :math:`n` and +:math:`X_n` to be independent identical distributed samples drawn from the +distribution of :math:`p(x)`. In practice: - >>> distribution = chaospy.Uniform(-1, 1) +.. code:: python + + >>> numpy.random.seed(1234) + >>> distribution = chaospy.Uniform(0, 1) >>> abscissas = distribution.sample(1000) >>> weights = 1./len(abscissas) -However, except for very high dimensional problems, Monte Carlo is quite an -inefficient way to perform numerical integration, and there exist quite a few -methods that performs better in most low-dimensional settings. If however, -Monto Carlo is your best choice, it might be worth taking a look at -:ref:`sampling`. - -.. note:: - Most quadrature rules optimized to a given weight function is referred to - as the :ref:`gaussian` rules. It does the embedding of the weight function - automatically as that is what it is designed for. For most other quadrature - rules, including a weight function is typically not canonical. This however - isn't very compatible with the Gaussian quadrature rules which take - probability density functions into account as part of their implementation. - It also does not match well with ``chaospy`` which assumes the density - weight function to be defined and incorporated implicit. - - To address this issue, the weight functions are incorporated into the - weight terms by substituting :math:`W^*_i \leftarraow W_i p(X_i)`, giving - us: - - .. math:: - \int p(x) g(x) dx \approx - \sum_i W_i p(X_i) g(X_i) = \sum_i W^{*}_i g(X_i) - - Which is the same format as the Gaussian quadrature rules. - - The consequence of this is that non-Gaussian quadrature rules only produces - the canonical weights for the probability distribution - ``chaospy.Uniform(0, 1)``, everything else is custom. To get around this - limitation, there are few workarounds: - - * Use a uniform distribution on an arbitrary interval ``Uniform(a, b)``, - and multiply the weight terms with the interval length: ``W *= (b-a)`` - * Use the quadrature rules directly from ``chaospy.quadrature.collection``. - * Adjust weights afterwards: ``W /= dist.pdf(X)`` - -To create quadrature abscissas and weights, use the -:func:`chaospy.generate_quadrature` function. Which type of quadrature to use -is defined by the flag ``rule``. This argument can either be the full name, or -a single letter representing the rule. These are as follows. - -.. _sparsegrid: - -Smolyak Sparse-Grid -------------------- - -As the number of dimensions increases linear, the number of samples increases -exponentially. This is known as the curse of dimensionality. Except for -switching to Monte Carlo integration, the is no way to completely guard against -this problem. However, there are some possibility to mitigate the problem -personally. One such strategy is to employ Smolyak sparse-grid quadrature. This -method uses a quadrature rule over a combination of different orders to tailor -a scheme that uses fewer abscissas points than a full tensor-product approach. - -To use Smolyak sparse-grid in ``chaospy``, just pass the flag ``sparse=True`` -to the :func:`chaospy.generate_quadrature` function. For example:: - - >>> distribution = chaospy.J( - ... chaospy.Uniform(0, 4), chaospy.Uniform(0, 4)) - >>> abscissas, weights = chaospy.generate_quadrature( - ... 3, distribution, sparse=True) - >>> abscissas.round(4) - array([[0., 0., 0., 1., 2., 2., 2., 2., 2., 3., 4., 4., 4.], - [0., 2., 4., 2., 0., 1., 2., 3., 4., 2., 0., 2., 4.]]) - >>> weights.round(4) - array([-0.0833, 0.2222, -0.0833, 0.4444, 0.2222, 0.4444, -1.3333, - 0.4444, 0.2222, 0.4444, -0.0833, 0.2222, -0.0833]) - -This compared to the full tensor-product grid:: +As an example, consider the problem of integrating the sinus function. +Applying the Monte Carlo integration technique is then as simple as: - >>> abscissas, weights = chaospy.generate_quadrature(3, distribution, sparse=False) - >>> abscissas.round(4) - array([[0., 0., 0., 0., 1., 1., 1., 1., 3., 3., 3., 3., 4., 4., 4., 4.], - [0., 1., 3., 4., 0., 1., 3., 4., 0., 1., 3., 4., 0., 1., 3., 4.]]) - >>> weights.round(4) - array([0.0031, 0.0247, 0.0247, 0.0031, 0.0247, 0.1975, 0.1975, 0.0247, - 0.0247, 0.1975, 0.1975, 0.0247, 0.0031, 0.0247, 0.0247, 0.0031]) +.. code:: python -The method works with all quadrature rules, but is known to be quite -inefficient when applied to rules that can not be nested. For example using -Gauss-Legendre samples:: + >>> approximation = numpy.sum(weights*numpy.sin(abscissas)) + >>> approximation + 0.4664434154275602 - >>> abscissas, weights = chaospy.generate_quadrature( - ... 6, distribution, rule="gauss_legendre", sparse=True) - >>> len(weights) - 140 - >>> abscissas, weights = chaospy.generate_quadrature( - ... 6, distribution, rule="gauss_legendre", sparse=False) - >>> len(weights) - 49 + >>> ref_solution = 1-numpy.cos(1) + >>> ref_solution + 0.45969769413186023 + >>> abs(ref_solution-approximation) + 0.006745721295699947 -.. note:: - Some quadrature rules are only partially nested at certain orders. These - include e.g. :func:`chaospy.quad_clenshaw_curtis`, - :func:`chaospy.quad_fejer` and :func:`chaospy.quad_newton_cotes`. To - exploit this nested-nes, the default behavior is to only include orders - that are properly nested. This implies that flipping the flag ``sparse`` - will result in a somewhat different scheme. To fix the scheme one way or - the other, explicitly include the flag ``growth=False`` or ``growth=True`` - respectively. +Though Monte Carlo is easy to implement, it also suffers from slow convergence +rate. It is quite inefficient at perform numerical integration, and there a few +methods that performs better in most low-dimensional settings. If however, +Monto Carlo is your best choice, it might be worth taking a look at +:ref:`sampling`. .. _gaussian: @@ -150,7 +98,9 @@ using the e.g. discretized Stieltjes algorithm. For example for the tailored quadrature rules defined above: -* Gauss-Hermit quadrature is tailored to the normal (Gaussian) distribution:: +* Gauss-Hermit quadrature is tailored to the normal (Gaussian) distribution: + + .. code:: python >>> distribution = chaospy.Normal(0, 1) >>> abscissas, weights = chaospy.generate_quadrature( @@ -160,7 +110,9 @@ For example for the tailored quadrature rules defined above: >>> weights.round(4) array([0.0026, 0.0886, 0.4088, 0.4088, 0.0886, 0.0026]) -* Gauss-Legendre quadrature is tailored to the Uniform distributions:: +* Gauss-Legendre quadrature is tailored to the Uniform distributions: + + .. code:: python >>> distribution = chaospy.Uniform(-1, 1) >>> abscissas, weights = chaospy.generate_quadrature( @@ -170,7 +122,9 @@ For example for the tailored quadrature rules defined above: >>> weights.round(4) array([0.0857, 0.1804, 0.234 , 0.234 , 0.1804, 0.0857]) -* Gauss-Jacobi quadrature is tailored to the Beta distribution:: +* Gauss-Jacobi quadrature is tailored to the Beta distribution: + + .. code:: python >>> distribution = chaospy.Beta(2, 4, lower=-1, upper=1) >>> abscissas, weights = chaospy.generate_quadrature( @@ -180,7 +134,9 @@ For example for the tailored quadrature rules defined above: >>> weights.round(4) array([0.0749, 0.272 , 0.355 , 0.2253, 0.0667, 0.0062]) -* Gauss-Laguerre quadrature is tailored to the Exponential distribution:: +* Gauss-Laguerre quadrature is tailored to the Exponential distribution: + + .. code:: python >>> distribution = chaospy.Exponential() >>> abscissas, weights = chaospy.generate_quadrature( @@ -190,7 +146,9 @@ For example for the tailored quadrature rules defined above: >>> weights.round(4) array([4.590e-01, 4.170e-01, 1.134e-01, 1.040e-02, 3.000e-04, 0.000e+00]) -* Generalized Gauss-Laguerre quadrature is tailored to the Gamma distribution:: +* Generalized Gauss-Laguerre quadrature is tailored to the Gamma distribution: + + .. code:: python >>> distribution = chaospy.Gamma(2, 4) >>> abscissas, weights = chaospy.generate_quadrature( @@ -205,7 +163,9 @@ as the distribution does not provide three terms recursion coefficients. In this scenario, the discretized counterpart is used instead as an approximation. For example, to mention a few: -* The Triangle distribution:: +* The Triangle distribution: + + .. code:: python >>> distribution = chaospy.Triangle(-1, 0, 1) >>> abscissas, weights = chaospy.generate_quadrature( @@ -215,7 +175,9 @@ For example, to mention a few: >>> weights.round(4) array([0.0295, 0.1475, 0.323 , 0.323 , 0.1475, 0.0295]) -* The Laplace distribution:: +* The Laplace distribution: + + .. code:: python >>> distribution = chaospy.Laplace(0, 1) >>> abscissas, weights = chaospy.generate_quadrature( @@ -225,7 +187,9 @@ For example, to mention a few: >>> weights.round(4) array([1.000e-04, 2.180e-02, 4.781e-01, 4.781e-01, 2.180e-02, 1.000e-04]) -* The Weibull distribution:: +* The Weibull distribution: + + .. code:: python >>> distribution = chaospy.Weibull() >>> abscissas, weights = chaospy.generate_quadrature( @@ -235,7 +199,9 @@ For example, to mention a few: >>> weights.round(4) array([4.589e-01, 4.170e-01, 1.134e-01, 1.040e-02, 3.000e-04, 0.000e+00]) -* The Rayleigh distribution:: +* The Rayleigh distribution: + + .. code:: python >>> distribution = chaospy.Rayleigh() >>> abscissas, weights = chaospy.generate_quadrature( @@ -246,7 +212,7 @@ For example, to mention a few: array([9.600e-02, 3.592e-01, 3.891e-01, 1.412e-01, 1.430e-02, 2.000e-04]) Statistician vs physicists --------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~ One of the more popular integration schemes when dealing with orthogonal polynomials are known as Gaussian quadrature. These are specially tailored @@ -269,7 +235,9 @@ example. \int_{-1}^1 0.5 g(x) dx \approx \sum_i W_i g(X_i) So to use ``chaospy`` to create a "true" Gaussian quadrature rule, one often has -to multiply the weights :math:`W_i` with some adjustment scalar. For example:: +to multiply the weights :math:`W_i` with some adjustment scalar. For example: + +.. code:: python >>> distribution = chaospy.Uniform(-1, 1) >>> N = 3 @@ -288,21 +256,172 @@ should be used. The various constants and distributions to achieve the various quadrature rules are as follows. -==================== ======================= ========================= =================== -Scheme Weight function Distribution Adjustment -==================== ======================= ========================= =================== -Hermite :math:`e^{-x^2}` ``Normal(0, 2**-0.5)`` :math:`\sqrt{\pi}` -Legendre :math:`1` ``Uniform(-1, 1)`` :math:`2` -Jakobi :math:`(1-x)^a(1+x)^b` ``Beta(a+1, b+1, -1, 1)`` :math:`2^{a+b}` -1. order Chebyshev :math:`1/\sqrt{1-x^2}` ``Beta(0.5, 0.5, -1, 1)`` :math:`1/2` -2. order Chebyshev :math:`\sqrt{1-x^2}` ``Beta(1.5, 1.5, -1, 1)`` :math:`2` -Laguerre :math:`e^{-x}` ``Exponential()`` :math:`1` -Generalized Laguerre :math:`x^a e^{-x}` ``Gamma(a+1)`` :math:`\Gamma(a+1)` -Gegenbaur :math:`(1-x^2)^{a-0.5}` ``Beta(a+.5,a+.5,-1,1)`` :math:`2^{2a-1}` -==================== ======================= ========================= =================== - -However, the list is not limited to these cases. Any and all valid weight -function are supported this way. However, not all weight functions does not -work very well. E.g. using the log-normal probability density function as -a weight function is known to scale badly. Which one works or not, depends on -context, so any non-standard use has to be done with some care. +======================================= ======================= ========================= =================== +Scheme Weight function Distribution Adjustment +======================================= ======================= ========================= =================== +:func:`~chaospy.quadrature.hermite` :math:`e^{-x^2}` ``Normal(0, 2**-0.5)`` :math:`\sqrt{\pi}` +:func:`~chaospy.quadrature.legendre` :math:`1` ``Uniform(-1, 1)`` :math:`2` +:func:`~chaospy.quadrature.jacobi` :math:`(1-x)^a(1+x)^b` ``Beta(a+1, b+1, -1, 1)`` :math:`2^{a+b}` +:func:`~chaospy.quadrature.chebyshev_1` :math:`1/\sqrt{1-x^2}` ``Beta(0.5, 0.5, -1, 1)`` :math:`1/2` +:func:`~chaospy.quadrature.chebyshev_2` :math:`\sqrt{1-x^2}` ``Beta(1.5, 1.5, -1, 1)`` :math:`2` +:func:`~chaospy.quadrature.laguerre` :math:`x^a e^{-x}` ``Gamma(a+1)`` :math:`\Gamma(a+1)` +:func:`~chaospy.quadrature.gegenbauer` :math:`(1-x^2)^{a-0.5}` ``Beta(a+.5,a+.5,-1,1)`` :math:`2^{2a-1}` +======================================= ======================= ========================= =================== + +The list is not limited to these cases. Any and all valid weight function are +supported this way. For the schemes listed are also functions where the +adjustments can be added by flag. Beyond the specific list here, any +distribution can be used to create the probabilist version. However, not all +weight functions work equally well. E.g. using the log-normal probability +density function as a weight function is known to scale badly. Which one works +or not, depends on context, so any non-standard use has to be done with some +care. + +Density as weight function +-------------------------- + +Most quadrature rules optimized to a given weight function is referred to as +the :ref:`gaussian` rules. It does the embedding of the weight function +automatically as that is what it is designed for. For most other quadrature +rules, including a weight function is typically not canonical. This however +isn't very compatible with the Gaussian quadrature rules which take probability +density functions directly into account as part of their implementation. It +also does not match well with ``chaospy`` which also assumes the density weight +function to be defined and incorporated implicit. + +To address this issue, the weight functions are incorporated into the weight +terms by substituting :math:`W^*_i \leftarraow W_i p(X_i)`, giving us: + +.. math:: + \int p(x) g(x) dx \approx + \sum_i W_i p(X_i) g(X_i) = \sum_i W^{*}_i g(X_i) + +Which is the same format as the Gaussian quadrature rules. + +The consequence of this is that non-Gaussian quadrature rules only produces the +canonical weights for the probability distribution ``chaospy.Uniform(0, 1)``, +everything else is custom. To get around this limitation, there are few +workarounds: + +* Use a uniform distribution on an arbitrary interval ``Uniform(a, b)``, and + multiply the weight terms with the interval length: ``W *= (b-a)`` +* Use the quadrature rules directly from ``chaospy.quadrature.collection``. +* Adjust weights afterwards: ``W /= dist.pdf(X)``. + +.. _sparsegrid: + +Smolyak Sparse-Grid +------------------- + +As the number of dimensions increases linear, the number of samples increases +exponentially. This is known as the curse of dimensionality. Except for +switching to Monte Carlo integration, the is no way to completely guard against +this problem. However, there are some possibility to mitigate the problem +personally. One such strategy is to employ Smolyak sparse-grid quadrature. This +method uses a quadrature rule over a combination of different orders to tailor +a scheme that uses fewer abscissas points than a full tensor-product approach. + +To use Smolyak sparse-grid in ``chaospy``, just pass the flag ``sparse=True`` +to the :func:`chaospy.generate_quadrature` function. For example: + +.. code:: python + + >>> distribution = chaospy.J( + ... chaospy.Uniform(0, 4), chaospy.Uniform(0, 4)) + >>> abscissas, weights = chaospy.generate_quadrature( + ... 3, distribution, sparse=True) + >>> abscissas.round(4) + array([[0., 0., 0., 1., 2., 2., 2., 2., 2., 3., 4., 4., 4.], + [0., 2., 4., 2., 0., 1., 2., 3., 4., 2., 0., 2., 4.]]) + >>> weights.round(4) + array([-0.0833, 0.2222, -0.0833, 0.4444, 0.2222, 0.4444, -1.3333, + 0.4444, 0.2222, 0.4444, -0.0833, 0.2222, -0.0833]) + +This compared to the full tensor-product grid: + +.. code:: python + + >>> abscissas, weights = chaospy.generate_quadrature(3, distribution, sparse=False) + >>> abscissas.round(4) + array([[0., 0., 0., 0., 1., 1., 1., 1., 3., 3., 3., 3., 4., 4., 4., 4.], + [0., 1., 3., 4., 0., 1., 3., 4., 0., 1., 3., 4., 0., 1., 3., 4.]]) + >>> weights.round(4) + array([0.0031, 0.0247, 0.0247, 0.0031, 0.0247, 0.1975, 0.1975, 0.0247, + 0.0247, 0.1975, 0.1975, 0.0247, 0.0031, 0.0247, 0.0247, 0.0031]) + +The method works with all quadrature rules, but is known to be quite +inefficient when applied to rules that can not be nested. For example using +Gauss-Legendre samples: + +.. code:: python + + >>> abscissas, weights = chaospy.generate_quadrature( + ... 6, distribution, rule="legendre", sparse=True) + >>> len(weights) + 139 + >>> abscissas, weights = chaospy.generate_quadrature( + ... 6, distribution, rule="legendre", sparse=False) + >>> len(weights) + 49 + +.. note:: + Some quadrature rules are only partially nested at certain orders. These + include e.g. :func:`chaospy.quadrature.clenshaw_curtis`, + :func:`chaospy.quadrature.fejer_1`, :func:`chaospy.quadrature.fejer_2` and + :func:`chaospy.quadrature.newton_cotes`. To exploit this nested-nes, the + default behavior is to only include orders that are properly nested. This + implies that flipping the flag ``sparse`` will result in a somewhat + different scheme. To fix the scheme one way or the other, explicitly + include the flag ``growth=False`` or ``growth=True`` respectively. + +Discrete probability distribution +--------------------------------- + +Quadrature rules with discrete distribution is not as often considered, as the +approximation is assumed to be of an integral. The analog to doing integration +in the discrete setting is to make a sum: + +.. math:: + + \sum_x p(x) g(x) + +As the domain of :math:`x` might be large, we propose the following +approximation strategy: Select abscissas evenly spaced over the discrete domain +and select the weights to be the point probability at each abscissas, but +adjusted so that the weights still sum to 1. When the number of abscissas reach +the size of the discrete domain, then it makes no more sense to increase the +number of nodes as the quadrature is no longer an approximation, but instead +the true formula. + +In practice we can use :func:`chaospy.quadrature.discrete`: + +.. code:: python + + >>> distribution = chaospy.DiscreteUniform(-4, 3) + >>> for order in range(4, 9): + ... abscissas, weights = chaospy.generate_quadrature( + ... order, distribution, rule="discrete") + ... print(order, abscissas.round(3), weights.round(3)) + 4 [[-4 -2 0 1 3]] [0.2 0.2 0.2 0.2 0.2] + 5 [[-4 -2 -1 0 1 3]] [0.167 0.167 0.167 0.167 0.167 0.167] + 6 [[-4 -3 -2 0 1 2 3]] [0.143 0.143 0.143 0.143 0.143 0.143 0.143] + 7 [[-4 -3 -2 -1 0 1 2 3]] [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125] + 8 [[-4 -3 -2 -1 0 1 2 3]] [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125] + +As the accuracy of discrete distribution plateau when all contained values are +included, there is no reason to increase the number of nodes after this point. + +The first few orders with exponential growth rule where the nodes are nested: + +.. code:: python + + >>> distribution = chaospy.DiscreteUniform(0, 10) + >>> for order in range(5): + ... abscissas, weights = chaospy.generate_quadrature( + ... order, distribution, rule="discrete", growth=True) + ... print(order, abscissas) + 0 [[5]] + 1 [[1 5 9]] + 2 [[1 3 5 7 9]] + 3 [[ 0 1 3 4 5 6 7 9 10]] + 4 [[ 0 1 2 3 4 5 6 7 8 9 10]] diff --git a/pyproject.toml b/pyproject.toml index e1b1fd90..2978915b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.masonry.api" [tool.poetry] name = "chaospy" -version = "4.2.4" +version = "4.3.0" description = "Numerical tool for perfroming uncertainty quantification" license = "MIT" authors = ["Jonathan Feinberg"] diff --git a/tests/recurrence/test_quadrature_creation.py b/tests/recurrence/test_quadrature_creation.py index 8d593b70..3eaf994b 100644 --- a/tests/recurrence/test_quadrature_creation.py +++ b/tests/recurrence/test_quadrature_creation.py @@ -14,7 +14,7 @@ def test_1d_quadrature_creation( analytical_distribution, recurrence_algorithm): """Check 1-D quadrature rule.""" - abscissas, weights = chaospy.quad_gaussian( + abscissas, weights = chaospy.gaussian( order=8, dist=analytical_distribution, recurrence_algorithm=recurrence_algorithm, @@ -35,7 +35,7 @@ def test_3d_quadrature_creation( analytical_distribution, recurrence_algorithm): """Check 3-D quadrature rule.""" distribution = chaospy.Iid(analytical_distribution, 3) - abscissas, weights = chaospy.quad_gaussian( + abscissas, weights = chaospy.gaussian( order=3, dist=distribution, recurrence_algorithm=recurrence_algorithm, From 90500805138cc751aae6d14da551ea1d2f390795 Mon Sep 17 00:00:00 2001 From: Jonathan Feinberg Date: Thu, 18 Mar 2021 10:18:08 +0100 Subject: [PATCH 2/4] refactor orthogonal->expansion --- CHANGELOG.rst | 38 ++++++ chaospy/__init__.py | 4 +- chaospy/distributions/approximation.py | 6 +- chaospy/distributions/operators/joint.py | 15 --- .../sampler/sequences/chebyshev.py | 12 +- .../distributions/sampler/sequences/halton.py | 55 ++++---- .../sampler/sequences/hammersley.py | 39 +++--- .../distributions/sampler/sequences/sobol.py | 8 +- .../sampler/sequences/van_der_corput.py | 50 +++---- chaospy/expansion/__init__.py | 38 ++++++ chaospy/expansion/chebyshev.py | 99 ++++++++++++++ chaospy/{orthogonal => expansion}/cholesky.py | 4 +- chaospy/{orthogonal => expansion}/frontend.py | 15 +-- chaospy/expansion/gegenbauer.py | 51 +++++++ .../{orthogonal => expansion}/gram_schmidt.py | 6 +- chaospy/expansion/hermite.py | 55 ++++++++ chaospy/expansion/jacobi.py | 40 ++++++ chaospy/{orthogonal => expansion}/lagrange.py | 14 +- chaospy/expansion/laguerre.py | 33 +++++ chaospy/expansion/legendre.py | 36 +++++ .../stieltjes.py} | 10 +- chaospy/orthogonal/__init__.py | 6 - chaospy/quadrature/__init__.py | 8 +- chaospy/quadrature/hermite.py | 2 +- chaospy/recurrence/stieltjes.py | 8 +- docs/reference/index.rst | 1 - docs/reference/orthogonal.rst | 23 ---- docs/reference/polynomial.rst | 126 +++++++++++------- docs/reference/quadrature.rst | 2 +- docs/user_guide/index.rst | 13 -- docs/user_guide/orthogonality.rst | 95 ------------- docs/user_guide/polynomial.rst | 96 +++++++++++++ pyproject.toml | 2 +- tests/recurrence/test_quadrature_creation.py | 4 +- tests/test_orth.py | 14 +- tests/test_stress.py | 17 +-- 36 files changed, 691 insertions(+), 354 deletions(-) create mode 100644 chaospy/expansion/__init__.py create mode 100644 chaospy/expansion/chebyshev.py rename chaospy/{orthogonal => expansion}/cholesky.py (98%) rename chaospy/{orthogonal => expansion}/frontend.py (92%) create mode 100644 chaospy/expansion/gegenbauer.py rename chaospy/{orthogonal => expansion}/gram_schmidt.py (92%) create mode 100644 chaospy/expansion/hermite.py create mode 100644 chaospy/expansion/jacobi.py rename chaospy/{orthogonal => expansion}/lagrange.py (87%) create mode 100644 chaospy/expansion/laguerre.py create mode 100644 chaospy/expansion/legendre.py rename chaospy/{orthogonal/three_terms_recurrence.py => expansion/stieltjes.py} (93%) delete mode 100644 chaospy/orthogonal/__init__.py delete mode 100644 docs/reference/orthogonal.rst delete mode 100644 docs/user_guide/orthogonality.rst diff --git a/CHANGELOG.rst b/CHANGELOG.rst index ab304241..dadd7957 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,44 @@ Master Branch ============= +Version 4.3.1 (2021-03-18) +========================== + +Refactoring `orthogonal -> expansion` module. + +ADDED: + * Dedicated classical orthogonal expansion schemes: + `chaospy.expansion.{chebyshev_1,chebyshev_2,gegenbauer,hermite,jacobi,laguerre,legendre}` +CHANGED: + * Function renames: + `chaospy.{orth_ttr,orth_chol,orth_gs,lagrange_polynomial} -> + chaospy.expansion.{stieltjes,cholesky,gram_schmidt,lagrange}` + * Docs update. + +Version 4.3.0 (2021-01-20) +========================== + +Refactoring `quadrature` module. + +ADDED: + * `chaospy.quadrature.fejer_1` is added. + * Dedicated classical quadrature schemes: + `chaospy.quadrature.{chebyshev,gegenbauer,hermite,jacobi,laguerre,legendre}` +CHANGED: + * Bound checks for the triangle distribution. (Thanks to @yoelcortes.) + * Refactored hypercube quadrature to common backend. This gives lots of flags + like `seqments` and + * Function renames: + `chaospy.quad_{clenshaw_curtis,discrete,fejer,gaussian,grid,gauss_lengendre,gauss_kronrod,gauss_lobatto,gauss_patterson,gauss_radau} -> + chaospy.quadrature.{clenshaw_curtis,fejer_2,gaussian,grid,legendre_proxy,kronrod,lobatto,patterson,radau}` + * Patterson growth rule changed from `0, 3, 7, ...` to `0, 1, 2, ...` but + maps backwards. Defaults to not have growth parameter, as setting it false + makes no sense. + * Renamed: `chaospy.generate_sparse_grid -> chaospy.quadrature.sparse_grid` +REMOVED: + * Genz-Keister quadrature `quad_genz_keister` is deprecated as it does not + fit the `chaospy` scheme very well. + Version 4.2.4 (2021-02-23) ========================== diff --git a/chaospy/__init__.py b/chaospy/__init__.py index 76fff84b..71fdd0b4 100644 --- a/chaospy/__init__.py +++ b/chaospy/__init__.py @@ -12,7 +12,7 @@ import chaospy.descriptives import chaospy.distributions -import chaospy.orthogonal +import chaospy.expansion import chaospy.spectral import chaospy.quadrature import chaospy.saltelli @@ -20,7 +20,7 @@ import chaospy.recurrence from chaospy.distributions import * -from chaospy.orthogonal import * +from chaospy.expansion import * from chaospy.spectral import * from chaospy.quadrature import * from chaospy.saltelli import * diff --git a/chaospy/distributions/approximation.py b/chaospy/distributions/approximation.py index e84b67a8..3f3ae897 100644 --- a/chaospy/distributions/approximation.py +++ b/chaospy/distributions/approximation.py @@ -171,12 +171,12 @@ def approximate_moment( k_loc (Sequence[int, ...]): The exponents of the moments of interest with ``shape == (dim,)``. order (int): - The quadrature order used in approximation. If omitted, calculated - to be ``1000/log2(len(distribution)+1)``. + The quadrature order used in approximation. If omitted, defaults to + ``1e5**(1./len(distribution))``. rule (str): Quadrature rule for integrating moments. kwargs: - Extra args passed to `chaospy.generate_quadrature`. + Extra args passed to :func:`chaospy.generate_quadrature`. Examples: >>> distribution = chaospy.Uniform(1, 4) diff --git a/chaospy/distributions/operators/joint.py b/chaospy/distributions/operators/joint.py index 544ccbd1..1f255b14 100644 --- a/chaospy/distributions/operators/joint.py +++ b/chaospy/distributions/operators/joint.py @@ -221,18 +221,3 @@ def _cache(self, idx, cache, get): if isinstance(out, chaospy.Distribution): return self return out - - -# def J(*args, **kwargs): -# """ -# Joint random variable. - -# Too be deprecated, use `chaospy.J` instead. - -# Args: -# args (chaospy.Distribution): -# Distribution to join together. -# """ -# logger = logging.getLogger(__name__) -# logger.warning("DepricationWarning: J to be replaced with Joint.") -# return Joint(*args, **kwargs) diff --git a/chaospy/distributions/sampler/sequences/chebyshev.py b/chaospy/distributions/sampler/sequences/chebyshev.py index 00c46a8c..46d7b0fc 100644 --- a/chaospy/distributions/sampler/sequences/chebyshev.py +++ b/chaospy/distributions/sampler/sequences/chebyshev.py @@ -49,7 +49,7 @@ def create_chebyshev_samples(order, dim=1): """ - Chebyshev sampling function. + Generate Chebyshev pseudo-random samples. Args: order (int): @@ -60,6 +60,16 @@ def create_chebyshev_samples(order, dim=1): Returns: samples following Chebyshev sampling scheme mapped to the ``[0, 1]^dim`` hyper-cube and ``shape == (dim, order)``. + + Examples: + >>> samples = chaospy.create_chebyshev_samples(6, 1) + >>> samples.round(4) + array([[0.0495, 0.1883, 0.3887, 0.6113, 0.8117, 0.9505]]) + >>> samples = chaospy.create_chebyshev_samples(3, 2) + >>> samples.round(3) + array([[0.146, 0.146, 0.146, 0.5 , 0.5 , 0.5 , 0.854, 0.854, 0.854], + [0.146, 0.5 , 0.854, 0.146, 0.5 , 0.854, 0.146, 0.5 , 0.854]]) + """ x_data = .5*numpy.cos(numpy.arange(order, 0, -1)*numpy.pi/(order+1)) + .5 x_data = utils.combine([x_data]*dim) diff --git a/chaospy/distributions/sampler/sequences/halton.py b/chaospy/distributions/sampler/sequences/halton.py index 02da8d98..ae54e6cf 100644 --- a/chaospy/distributions/sampler/sequences/halton.py +++ b/chaospy/distributions/sampler/sequences/halton.py @@ -1,30 +1,4 @@ -""" -Create samples from the `Halton sequence`_. - -In statistics, Halton sequences are sequences used to generate points in space -for numerical methods such as Monte Carlo simulations. Although these sequences -are deterministic, they are of low discrepancy, that is, appear to be random -for many purposes. They were first introduced in 1960 and are an example of -a quasi-random number sequence. They generalise the one-dimensional van der -Corput sequences. - -Example usage -------------- - -Standard usage:: - - >>> distribution = chaospy.J(chaospy.Uniform(0, 1), chaospy.Uniform(0, 1)) - >>> samples = distribution.sample(3, rule="halton") - >>> samples.round(4) - array([[0.125 , 0.625 , 0.375 ], - [0.4444, 0.7778, 0.2222]]) - >>> samples = distribution.sample(4, rule="halton") - >>> samples.round(4) - array([[0.125 , 0.625 , 0.375 , 0.875 ], - [0.4444, 0.7778, 0.2222, 0.5556]]) - -.. _Halton sequence: https://en.wikipedia.org/wiki/Halton_sequence -""" +"""Create samples from the Halton sequence.""" import numpy from .van_der_corput import create_van_der_corput_samples @@ -33,9 +7,15 @@ def create_halton_samples(order, dim=1, burnin=-1, primes=()): """ - Create Halton sequence. + Create samples from the Halton sequence. - For ``dim == 1`` the sequence falls back to Van Der Corput sequence. + In statistics, Halton sequences are sequences used to generate points in + space for numerical methods such as Monte Carlo simulations. Although these + sequences are deterministic, they are of low discrepancy, that is, appear + to be random for many purposes. They were first introduced in 1960 and are + an example of a quasi-random number sequence. They generalise the + one-dimensional van der Corput sequences. For ``dim == 1`` the sequence + falls back to Van Der Corput sequence. Args: order (int): @@ -49,8 +29,21 @@ def create_halton_samples(order, dim=1, burnin=-1, primes=()): The (non-)prime base to calculate values along each axis. If empty, growing prime values starting from 2 will be used. - Returns (numpy.ndarray): - Halton sequence with ``shape == (dim, order)``. + Returns: + (numpy.ndarray): + Halton sequence with ``shape == (dim, order)``. + + Examples: + >>> distribution = chaospy.J(chaospy.Uniform(0, 1), chaospy.Uniform(0, 1)) + >>> samples = distribution.sample(3, rule="halton") + >>> samples.round(4) + array([[0.125 , 0.625 , 0.375 ], + [0.4444, 0.7778, 0.2222]]) + >>> samples = distribution.sample(4, rule="halton") + >>> samples.round(4) + array([[0.125 , 0.625 , 0.375 , 0.875 ], + [0.4444, 0.7778, 0.2222, 0.5556]]) + """ primes = list(primes) if not primes: diff --git a/chaospy/distributions/sampler/sequences/hammersley.py b/chaospy/distributions/sampler/sequences/hammersley.py index a82b6775..15eb9211 100644 --- a/chaospy/distributions/sampler/sequences/hammersley.py +++ b/chaospy/distributions/sampler/sequences/hammersley.py @@ -1,26 +1,4 @@ -""" -Create samples from the `Hammersley set`_. - -The Hammersley set is equivalent to the Halton sequence, except for one -dimension is replaced with a regular grid. - -Example usage -------------- - -Standard usage:: - - >>> distribution = chaospy.J(chaospy.Uniform(0, 1), chaospy.Uniform(0, 1)) - >>> samples = distribution.sample(3, rule="hammersley") - >>> samples.round(4) - array([[0.75 , 0.125, 0.625], - [0.25 , 0.5 , 0.75 ]]) - >>> samples = distribution.sample(4, rule="hammersley") - >>> samples.round(4) - array([[0.75 , 0.125, 0.625, 0.375], - [0.2 , 0.4 , 0.6 , 0.8 ]]) - -.. _Hammersley set: https://en.wikipedia.org/wiki/Low-discrepancy_sequence#Hammersley_set -""" +"""Create samples from the Hammersley set.""" import numpy from .halton import create_halton_samples @@ -30,7 +8,8 @@ def create_hammersley_samples(order, dim=1, burnin=-1, primes=()): """ Create samples from the Hammersley set. - For ``dim == 1`` the sequence falls back to Van Der Corput sequence. + The Hammersley set is equivalent to the Halton sequence, except for one + dimension is replaced with a regular grid. Args: order (int): @@ -47,6 +26,18 @@ def create_hammersley_samples(order, dim=1, burnin=-1, primes=()): Returns: (numpy.ndarray): Hammersley set with ``shape == (dim, order)``. + + Examples: + >>> distribution = chaospy.J(chaospy.Uniform(0, 1), chaospy.Uniform(0, 1)) + >>> samples = distribution.sample(3, rule="hammersley") + >>> samples.round(4) + array([[0.75 , 0.125, 0.625], + [0.25 , 0.5 , 0.75 ]]) + >>> samples = distribution.sample(4, rule="hammersley") + >>> samples.round(4) + array([[0.75 , 0.125, 0.625, 0.375], + [0.2 , 0.4 , 0.6 , 0.8 ]]) + """ if dim == 1: return create_halton_samples( diff --git a/chaospy/distributions/sampler/sequences/sobol.py b/chaospy/distributions/sampler/sequences/sobol.py index 8294b552..35735a27 100644 --- a/chaospy/distributions/sampler/sequences/sobol.py +++ b/chaospy/distributions/sampler/sequences/sobol.py @@ -1,10 +1,5 @@ """ - -Example usage -------------- - -Standard usage:: - +Generates samples from the Sobol sequence. Papers:: @@ -35,7 +30,6 @@ Preprint IPM Akad. Nauk SSSR, Number 40, Moscow 1976. -.. _Sobel sequence: https://en.wikipedia.org/wiki/Sobol_sequence """ import math diff --git a/chaospy/distributions/sampler/sequences/van_der_corput.py b/chaospy/distributions/sampler/sequences/van_der_corput.py index eb1a0909..fcab00e8 100644 --- a/chaospy/distributions/sampler/sequences/van_der_corput.py +++ b/chaospy/distributions/sampler/sequences/van_der_corput.py @@ -1,36 +1,20 @@ -""" -Create `Van Der Corput` low discrepancy sequence samples. - -A van der Corput sequence is an example of the simplest one-dimensional -low-discrepancy sequence over the unit interval; it was first described in 1935 -by the Dutch mathematician J. G. van der Corput. It is constructed by reversing -the base-n representation of the sequence of natural numbers (1, 2, 3, ...). - -In practice, use Halton sequence instead of Van Der Corput, as it is the -same, but generalized to work in multiple dimensions. - -Example usage -------------- - -Using base 10:: - - >>> create_van_der_corput_samples(range(11), number_base=10) - array([0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 0.01, 0.11]) - -Using base 2:: - - >>> create_van_der_corput_samples(range(8), number_base=2) - array([0.5 , 0.25 , 0.75 , 0.125 , 0.625 , 0.375 , 0.875 , 0.0625]) - -.. Van Der Corput: https://en.wikipedia.org/wiki/Van_der_Corput_sequence -""" +"""Create Van Der Corput low discrepancy sequence samples.""" from __future__ import division import numpy def create_van_der_corput_samples(idx, number_base=2): """ - Van der Corput samples. + Create Van Der Corput low discrepancy sequence samples. + + A van der Corput sequence is an example of the simplest one-dimensional + low-discrepancy sequence over the unit interval; it was first described in + 1935 by the Dutch mathematician J. G. van der Corput. It is constructed by + reversing the base-n representation of the sequence of natural numbers + :math:`(1, 2, 3, ...)`. + + In practice, use Halton sequence instead of Van Der Corput, as it is the + same, but generalized to work in multiple dimensions. Args: idx (int, numpy.ndarray): @@ -39,8 +23,16 @@ def create_van_der_corput_samples(idx, number_base=2): number_base (int): The numerical base from where to create the samples from. - Returns (float, numpy.ndarray): - Van der Corput samples. + Returns: + (numpy.ndarray): + Van der Corput samples. + + Examples: + #>>> chaospy.create_van_der_corput_samples(range(11), number_base=10) + #array([0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 0.01, 0.11]) + #>>> chaospy.create_van_der_corput_samples(range(8), number_base=2) + #array([0.5 , 0.25 , 0.75 , 0.125 , 0.625 , 0.375 , 0.875 , 0.0625]) + """ assert number_base > 1 diff --git a/chaospy/expansion/__init__.py b/chaospy/expansion/__init__.py new file mode 100644 index 00000000..cd5e9927 --- /dev/null +++ b/chaospy/expansion/__init__.py @@ -0,0 +1,38 @@ +r"""Collection of polynomial expansion constructors.""" +import logging +from functools import wraps + +from .chebyshev import chebyshev_1, chebyshev_2 +from .cholesky import cholesky +from .frontend import generate_expansion +from .gegenbauer import gegenbauer +from .gram_schmidt import gram_schmidt +from .hermite import hermite +from .jacobi import jacobi +from .stieltjes import stieltjes +from .lagrange import lagrange +from .laguerre import laguerre +from .legendre import legendre + +__all__ = ["generate_expansion"] + + +def expansion_deprecation_warning(name, func): + + @wraps(func) + def wrapped(*args, **kwargs): + """Function wrapper adds warnings.""" + logger = logging.getLogger(__name__) + logger.warning("chaospy.%s name is to be deprecated; " + "Use chaospy.expansion.%s instead", + name, func.__name__) + return func(*args, **kwargs) + + globals()[name] = wrapped + __all__.append(name) + + +expansion_deprecation_warning("orth_ttr", stieltjes) +expansion_deprecation_warning("orth_chol", cholesky) +expansion_deprecation_warning("orth_gs", gram_schmidt) +expansion_deprecation_warning("lagrange_polynomial", lagrange) diff --git a/chaospy/expansion/chebyshev.py b/chaospy/expansion/chebyshev.py new file mode 100644 index 00000000..00a239a8 --- /dev/null +++ b/chaospy/expansion/chebyshev.py @@ -0,0 +1,99 @@ +"""Chebyshev polynomials of the first kind.""" +import numpy +import chaospy + + +def chebyshev_1( + order, + lower=-1, + upper=1, + physicist=False, + normed=False, + retall=False, +): + """ + Chebyshev polynomials of the first kind. + + Args: + order (int): + The polynomial order. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + (numpoly.ndpoly, numpy.ndarray): + Chebyshev polynomial expansion. Norms of the orthogonal + expansion on the form ``E(orth**2, dist)``. + + Examples: + >>> polynomials, norms = chaospy.expansion.chebyshev_1(4, retall=True) + >>> polynomials + polynomial([1.0, q0, q0**2-0.5, q0**3-0.75*q0, q0**4-q0**2+0.125]) + >>> norms + array([1. , 0.5 , 0.125 , 0.03125 , 0.0078125]) + >>> chaospy.expansion.chebyshev_1(3, physicist=True) + polynomial([1.0, q0, 2.0*q0**2-1.0, 4.0*q0**3-2.5*q0]) + >>> chaospy.expansion.chebyshev_1(3, lower=0.5, upper=1.5, normed=True).round(3) + polynomial([1.0, 2.828*q0-2.828, 11.314*q0**2-22.627*q0+9.899, + 45.255*q0**3-135.765*q0**2+127.279*q0-36.77]) + + """ + multiplier = 1+numpy.arange(order).astype(bool) if physicist else 1 + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Beta(0.5, 0.5, lower, upper), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials + + +def chebyshev_2( + order, + lower=-1, + upper=1, + physicist=False, + normed=False, + retall=False, +): + """ + Chebyshev polynomials of the second kind. + + Args: + order (int): + The quadrature order. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Returns: + (numpoly.ndpoly, numpy.ndarray): + Chebyshev polynomial expansion. Norms of the orthogonal + expansion on the form ``E(orth**2, dist)``. + + Examples: + >>> polynomials, norms = chaospy.expansion.chebyshev_2(4, retall=True) + >>> polynomials + polynomial([1.0, q0, q0**2-0.25, q0**3-0.5*q0, q0**4-0.75*q0**2+0.0625]) + >>> norms + array([1. , 0.25 , 0.0625 , 0.015625 , 0.00390625]) + >>> chaospy.expansion.chebyshev_2(3, physicist=True) + polynomial([1.0, 2.0*q0, 4.0*q0**2-0.5, 8.0*q0**3-2.0*q0]) + >>> chaospy.expansion.chebyshev_2(3, lower=0.5, upper=1.5, normed=True).round(3) + polynomial([1.0, 4.0*q0-4.0, 16.0*q0**2-32.0*q0+15.0, + 64.0*q0**3-192.0*q0**2+184.0*q0-56.0]) + + """ + multiplier = 2 if physicist else 1 + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Beta(1.5, 1.5, lower, upper), multiplier=multiplier) + if normed: + polynomials= chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/orthogonal/cholesky.py b/chaospy/expansion/cholesky.py similarity index 98% rename from chaospy/orthogonal/cholesky.py rename to chaospy/expansion/cholesky.py index baf53ae2..db993e3d 100644 --- a/chaospy/orthogonal/cholesky.py +++ b/chaospy/expansion/cholesky.py @@ -29,7 +29,7 @@ -def orth_chol( +def cholesky( order, dist, normed=False, @@ -66,7 +66,7 @@ def orth_chol( Examples: >>> distribution = chaospy.Normal() - >>> expansion, norms = chaospy.orth_chol(3, distribution, retall=True) + >>> expansion, norms = chaospy.expansion.cholesky(3, distribution, retall=True) >>> expansion.round(4) polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0]) >>> norms diff --git a/chaospy/orthogonal/frontend.py b/chaospy/expansion/frontend.py similarity index 92% rename from chaospy/orthogonal/frontend.py rename to chaospy/expansion/frontend.py index 10588c76..fff45876 100644 --- a/chaospy/orthogonal/frontend.py +++ b/chaospy/expansion/frontend.py @@ -1,18 +1,17 @@ """Frontend function for generating polynomial expansions.""" -from .three_terms_recurrence import orth_ttr -from .cholesky import orth_chol -from .gram_schmidt import orth_gs -from .lagrange import lagrange_polynomial +from .stieltjes import stieltjes +from .cholesky import cholesky +from .gram_schmidt import gram_schmidt EXPANSION_NAMES = { - "ttr": "three_terms_recurrence", "three_terms_recurrence": "three_terms_recurrence", + "ttr": "stieltjes", "three_terms_recurrence": "stieltjes", "stieltjes": "stieltjes", "chol": "cholesky", "cholesky": "cholesky", "gs": "gram_schmidt", "gram_schmidt": "gram_schmidt", } EXPANSION_FUNCTIONS = { - "three_terms_recurrence": orth_ttr, - "cholesky": orth_chol, - "gram_schmidt": orth_gs, + "stieltjes": stieltjes, + "cholesky": cholesky, + "gram_schmidt": gram_schmidt, } diff --git a/chaospy/expansion/gegenbauer.py b/chaospy/expansion/gegenbauer.py new file mode 100644 index 00000000..583f4487 --- /dev/null +++ b/chaospy/expansion/gegenbauer.py @@ -0,0 +1,51 @@ +import numpy +import chaospy + + +def gegenbauer( + order, + alpha, + lower=-1, + upper=1, + physicist=False, + normed=False, + retall=False, +): + """ + Gegenbauer polynomials. + + Args: + order (int): + The polynomial order. + alpha (float): + Gegenbauer shape parameter. + lower (float): + Lower bound for the integration interval. + upper (float): + Upper bound for the integration interval. + physicist (bool): + Use physicist weights instead of probabilist. + + Examples: + >>> polynomials, norms = chaospy.expansion.gegenbauer(4, 1, retall=True) + >>> polynomials + polynomial([1.0, q0, q0**2-0.25, q0**3-0.5*q0, q0**4-0.75*q0**2+0.0625]) + >>> norms + array([1. , 0.25 , 0.0625 , 0.015625 , 0.00390625]) + >>> chaospy.expansion.gegenbauer(3, 1, physicist=True) + polynomial([1.0, 2.0*q0, 4.0*q0**2-0.5, 8.0*q0**3-2.0*q0]) + >>> chaospy.expansion.gegenbauer(3, 1, lower=0.5, upper=1.5, normed=True).round(3) + polynomial([1.0, 4.0*q0-4.0, 16.0*q0**2-32.0*q0+15.0, + 64.0*q0**3-192.0*q0**2+184.0*q0-56.0]) + + """ + multiplier = 1 + if physicist: + multiplier = numpy.arange(1, order+1) + multiplier = 2*(multiplier+alpha-1)/multiplier + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Beta(alpha+0.5, alpha+0.5, lower, upper), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/orthogonal/gram_schmidt.py b/chaospy/expansion/gram_schmidt.py similarity index 92% rename from chaospy/orthogonal/gram_schmidt.py rename to chaospy/expansion/gram_schmidt.py index f5b0af66..b0a419d4 100644 --- a/chaospy/orthogonal/gram_schmidt.py +++ b/chaospy/expansion/gram_schmidt.py @@ -6,7 +6,7 @@ import chaospy -def orth_gs(order, dist, normed=False, graded=True, reverse=True, +def gram_schmidt(order, dist, normed=False, graded=True, reverse=True, retall=False, cross_truncation=1., **kws): """ Gram-Schmidt process for generating orthogonal polynomials. @@ -41,12 +41,12 @@ def orth_gs(order, dist, normed=False, graded=True, reverse=True, Examples: >>> distribution = chaospy.J(chaospy.Normal(), chaospy.Normal()) - >>> polynomials, norms = chaospy.orth_gs(2, distribution, retall=True) + >>> polynomials, norms = chaospy.expansion.gram_schmidt(2, distribution, retall=True) >>> polynomials.round(10) polynomial([1.0, q1, q0, q1**2-1.0, q0*q1, q0**2-1.0]) >>> norms.round(10) array([1., 1., 1., 2., 1., 2.]) - >>> polynomials = chaospy.orth_gs(2, distribution, normed=True) + >>> polynomials = chaospy.expansion.gram_schmidt(2, distribution, normed=True) >>> polynomials.round(3) polynomial([1.0, q1, q0, 0.707*q1**2-0.707, q0*q1, 0.707*q0**2-0.707]) diff --git a/chaospy/expansion/hermite.py b/chaospy/expansion/hermite.py new file mode 100644 index 00000000..94cfcb43 --- /dev/null +++ b/chaospy/expansion/hermite.py @@ -0,0 +1,55 @@ +"""Hermite orthogonal polynomial expansion.""" +import numpy +import chaospy + + +def hermite( + order, + mu=0., + sigma=1., + physicist=False, + normed=False, + retall=False, +): + """ + Hermite orthogonal polynomial expansion. + + Args: + order (int): + The quadrature order. + mu (float): + Non-centrality parameter. + sigma (float): + Scale parameter. + physicist (bool): + Use physicist weights instead of probabilist variant. + normed (bool): + If True orthonormal polynomials will be used. + retall (bool): + If true return numerical stabilized norms as well. Roughly the same + as ``cp.E(orth**2, dist)``. + + Returns: + (numpoly.ndpoly, numpy.ndarray): + Hermite polynomial expansion. Norms of the orthogonal + expansion on the form ``E(orth**2, dist)``. + + Examples: + >>> polynomials, norms = chaospy.expansion.hermite(4, retall=True) + >>> polynomials + polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0, q0**4-6.0*q0**2+3.0]) + >>> norms + array([ 1., 1., 2., 6., 24.]) + >>> chaospy.expansion.hermite(3, physicist=True) + polynomial([1.0, 2.0*q0, 4.0*q0**2-2.0, 8.0*q0**3-12.0*q0]) + >>> chaospy.expansion.hermite(3, sigma=2.5, normed=True).round(3) + polynomial([1.0, 0.4*q0, 0.113*q0**2-0.707, 0.026*q0**3-0.49*q0]) + + """ + multiplier = 2 if physicist else 1 + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Normal(mu, sigma), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/expansion/jacobi.py b/chaospy/expansion/jacobi.py new file mode 100644 index 00000000..d7e46be9 --- /dev/null +++ b/chaospy/expansion/jacobi.py @@ -0,0 +1,40 @@ +import numpy +import chaospy + + +def jacobi( + order, + alpha, + beta, + lower=-1, + upper=1, + physicist=False, + normed=False, + retall=False, +): + """ + Jacobi polynomial expansion. + + Examples: + >>> polynomials, norms = chaospy.expansion.jacobi(4, 0.5, 0.5, retall=True) + >>> polynomials + polynomial([1.0, q0, q0**2-0.5, q0**3-0.75*q0, q0**4-q0**2+0.125]) + >>> norms + array([1. , 0.5 , 0.125 , 0.03125 , 0.0078125]) + >>> chaospy.expansion.jacobi(3, 0.5, 0.5, physicist=True).round(4) + polynomial([1.0, 1.5*q0, 2.5*q0**2-0.8333, 4.375*q0**3-2.1146*q0]) + >>> chaospy.expansion.jacobi(3, 1.5, 0.5, normed=True) + polynomial([1.0, 2.0*q0, 4.0*q0**2-1.0, 8.0*q0**3-4.0*q0]) + + """ + multiplier = 1 + if physicist: + multiplier = numpy.arange(1, order+1) + multiplier = ((2*multiplier+alpha+beta-1)*(2*multiplier+alpha+beta)/ + (2*multiplier*(multiplier+alpha+beta))) + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Beta(alpha, beta, lower=lower, upper=upper), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/orthogonal/lagrange.py b/chaospy/expansion/lagrange.py similarity index 87% rename from chaospy/orthogonal/lagrange.py rename to chaospy/expansion/lagrange.py index a7499a0b..fc3361bb 100644 --- a/chaospy/orthogonal/lagrange.py +++ b/chaospy/expansion/lagrange.py @@ -4,9 +4,9 @@ import numpoly -def lagrange_polynomial(abscissas, graded=True, reverse=True, sort=None): +def lagrange(abscissas, graded=True, reverse=True, sort=None): """ - Create Lagrange polynomials. + Create Lagrange polynomial expansion. Args: abscissas (numpy.ndarray): @@ -22,13 +22,13 @@ def lagrange_polynomial(abscissas, graded=True, reverse=True, sort=None): considered bigger than ``q0**3*q1``, instead of the opposite. Example: - >>> chaospy.lagrange_polynomial([4]).round(4) + >>> chaospy.expansion.lagrange([4]).round(4) polynomial([4.0]) - >>> chaospy.lagrange_polynomial([-10, 10]).round(4) + >>> chaospy.expansion.lagrange([-10, 10]).round(4) polynomial([-0.05*q0+0.5, 0.05*q0+0.5]) - >>> chaospy.lagrange_polynomial([-1, 0, 1]).round(4) + >>> chaospy.expansion.lagrange([-1, 0, 1]).round(4) polynomial([0.5*q0**2-0.5*q0, -q0**2+1.0, 0.5*q0**2+0.5*q0]) - >>> poly = chaospy.lagrange_polynomial([[1, 0, 1], [0, 1, 2]]) + >>> poly = chaospy.expansion.lagrange([[1, 0, 1], [0, 1, 2]]) >>> poly.round(4) polynomial([-0.5*q1+0.5*q0+0.5, -q0+1.0, 0.5*q1+0.5*q0-0.5]) >>> poly([1, 0, 1], [0, 1, 2]).round(14) @@ -37,7 +37,7 @@ def lagrange_polynomial(abscissas, graded=True, reverse=True, sort=None): [0., 0., 1.]]) >>> nodes = numpy.array([[ 0.17, 0.15, 0.17, 0.19], ... [14.94, 16.69, 16.69, 16.69]]) - >>> poly = chaospy.lagrange_polynomial(nodes) # doctest: +IGNORE_EXCEPTION_DETAIL + >>> poly = chaospy.expansion.lagrange(nodes) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): ... LinAlgError: Lagrange abscissas resulted in invertible matrix diff --git a/chaospy/expansion/laguerre.py b/chaospy/expansion/laguerre.py new file mode 100644 index 00000000..2622187a --- /dev/null +++ b/chaospy/expansion/laguerre.py @@ -0,0 +1,33 @@ +import numpy +import chaospy + + +def laguerre( + order, + alpha=0., + physicist=False, + normed=False, + retall=False, +): + """ + Examples: + >>> polynomials, norms = chaospy.expansion.laguerre(3, retall=True) + >>> polynomials + polynomial([1.0, q0-1.0, q0**2-4.0*q0+2.0, q0**3-9.0*q0**2+18.0*q0-6.0]) + >>> norms + array([ 1., 1., 4., 36.]) + >>> chaospy.expansion.laguerre(3, physicist=True).round(5) + polynomial([1.0, -q0+1.0, 0.5*q0**2-2.0*q0+2.0, + -0.16667*q0**3+1.5*q0**2-5.33333*q0+4.66667]) + >>> chaospy.expansion.laguerre(3, alpha=2, normed=True).round(3) + polynomial([1.0, 0.577*q0-1.732, 0.204*q0**2-1.633*q0+2.449, + 0.053*q0**3-0.791*q0**2+3.162*q0-3.162]) + + """ + multiplier = -1./numpy.arange(1, order+1) if physicist else 1. + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Gamma(alpha+1), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/expansion/legendre.py b/chaospy/expansion/legendre.py new file mode 100644 index 00000000..d274c3b3 --- /dev/null +++ b/chaospy/expansion/legendre.py @@ -0,0 +1,36 @@ +import numpy +import chaospy + + +def legendre( + order, + lower=-1, + upper=1, + physicist=False, + normed=False, + retall=False, +): + """ + Examples: + >>> polynomials, norms = chaospy.expansion.legendre(3, retall=True) + >>> polynomials.round(5) + polynomial([1.0, q0, q0**2-0.33333, q0**3-0.6*q0]) + >>> norms + array([1. , 0.33333333, 0.08888889, 0.02285714]) + >>> chaospy.expansion.legendre(3, physicist=True).round(3) + polynomial([1.0, 1.5*q0, 2.5*q0**2-0.556, 4.375*q0**3-1.672*q0]) + >>> chaospy.expansion.legendre(3, lower=0, upper=1, normed=True).round(3) + polynomial([1.0, 3.464*q0-1.732, 13.416*q0**2-13.416*q0+2.236, + 52.915*q0**3-79.373*q0**2+31.749*q0-2.646]) + + """ + multiplier = 1. + if physicist: + multiplier = numpy.arange(1, order+1) + multiplier = (2*multiplier+1)/(multiplier+1) + _, [polynomials], [norms] = chaospy.recurrence.analytical_stieltjes( + order, chaospy.Uniform(lower, upper), multiplier=multiplier) + if normed: + polynomials = chaospy.true_divide(polynomials, numpy.sqrt(norms)) + norms[:] = 1. + return (polynomials, norms) if retall else polynomials diff --git a/chaospy/orthogonal/three_terms_recurrence.py b/chaospy/expansion/stieltjes.py similarity index 93% rename from chaospy/orthogonal/three_terms_recurrence.py rename to chaospy/expansion/stieltjes.py index 23820068..2a0808b2 100644 --- a/chaospy/orthogonal/three_terms_recurrence.py +++ b/chaospy/expansion/stieltjes.py @@ -38,10 +38,10 @@ discretized Stieltjes method (described in the `paper by Golub and Welsch`_). In ``chaospy`` constructing orthogonal polynomial using the three term -recurrence scheme can be done through ``orth_ttr``. For example:: +recurrence scheme can be done through ``stieltjes``. For example:: >>> dist = chaospy.Iid(chaospy.Gamma(1), 2) - >>> orths = chaospy.orth_ttr(2, dist) + >>> orths = chaospy.expansion.stieltjes(2, dist) >>> orths.round(4) polynomial([1.0, q1-1.0, q0-1.0, q1**2-4.0*q1+2.0, q0*q1-q1-q0+1.0, q0**2-4.0*q0+2.0]) @@ -59,7 +59,7 @@ import chaospy -def orth_ttr(order, dist, normed=False, graded=True, reverse=True, +def stieltjes(order, dist, normed=False, graded=True, reverse=True, retall=False, cross_truncation=1.): """ Create orthogonal polynomial expansion from three terms recurrence formula. @@ -98,12 +98,12 @@ def orth_ttr(order, dist, normed=False, graded=True, reverse=True, Examples: >>> distribution = chaospy.J(chaospy.Normal(), chaospy.Normal()) - >>> polynomials, norms = chaospy.orth_ttr(2, distribution, retall=True) + >>> polynomials, norms = chaospy.expansion.stieltjes(2, distribution, retall=True) >>> polynomials.round(10) polynomial([1.0, q1, q0, q1**2-1.0, q0*q1, q0**2-1.0]) >>> norms.round(10) array([1., 1., 1., 2., 1., 2.]) - >>> polynomials = chaospy.orth_ttr(2, distribution, normed=True) + >>> polynomials = chaospy.expansion.stieltjes(2, distribution, normed=True) >>> polynomials.round(3) polynomial([1.0, q1, q0, 0.707*q1**2-0.707, q0*q1, 0.707*q0**2-0.707]) diff --git a/chaospy/orthogonal/__init__.py b/chaospy/orthogonal/__init__.py deleted file mode 100644 index 4d624d9f..00000000 --- a/chaospy/orthogonal/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -r"""Collection of polynomial expansion constructors.""" -from .frontend import generate_expansion -from .three_terms_recurrence import orth_ttr -from .lagrange import lagrange_polynomial -from .gram_schmidt import orth_gs -from .cholesky import orth_chol diff --git a/chaospy/quadrature/__init__.py b/chaospy/quadrature/__init__.py index afc93b32..82d5373d 100644 --- a/chaospy/quadrature/__init__.py +++ b/chaospy/quadrature/__init__.py @@ -25,6 +25,8 @@ from .patterson import patterson from .radau import radau +__all__ = ["generate_quadrature", "sparse_grid", "combine"] + INTEGRATION_COLLECTION = { "clenshaw_curtis": clenshaw_curtis, @@ -43,11 +45,8 @@ } -def quadrature_deprecation_warning(name, func=None): +def quadrature_deprecation_warning(name, func): """Announce deprecation warning for quad-func.""" - - if func is None: - func = globals()[name] quad_name = "quad_%s" % name @wraps(func) @@ -60,6 +59,7 @@ def wrapped(*args, **kwargs): return func(*args, **kwargs) globals()[quad_name] = wrapped + __all__.append(quad_name) quadrature_deprecation_warning("clenshaw_curtis", clenshaw_curtis) quadrature_deprecation_warning("discrete", discrete) diff --git a/chaospy/quadrature/hermite.py b/chaospy/quadrature/hermite.py index 16d9870b..fa4f276a 100644 --- a/chaospy/quadrature/hermite.py +++ b/chaospy/quadrature/hermite.py @@ -27,7 +27,7 @@ def hermite(order, mu=0., sigma=1., physicist=False): sigma (float): Scale parameter. physicist (bool): - Use physicist weights instead of probabilist. + Use physicist weights instead of probabilist variant. Returns: abscissas (numpy.ndarray): diff --git a/chaospy/recurrence/stieltjes.py b/chaospy/recurrence/stieltjes.py index dfa21d55..2ef86290 100644 --- a/chaospy/recurrence/stieltjes.py +++ b/chaospy/recurrence/stieltjes.py @@ -154,19 +154,21 @@ def discretized_stieltjes( return coeffs, orths, norms -def analytical_stieltjes(order, dist): +def analytical_stieltjes(order, dist, multiplier=1): """Analytical Stieltjes' method""" dimensions = len(dist) mom_order = numpy.arange(order+1).repeat(dimensions) mom_order = mom_order.reshape(order+1, dimensions).T coeffs = dist.ttr(mom_order) coeffs[1, :, 0] = 1. + orders = numpy.arange(order, dtype=int) + multiplier, orders = numpy.broadcast_arrays(multiplier, orders) var = numpoly.variable(dimensions) orth = [numpy.zeros(dimensions), numpy.ones(dimensions)] - for order_ in range(order): + for order_, multiplier_ in zip(orders, multiplier): orth.append( - (var-coeffs[0, :, order_])*orth[-1]-coeffs[1, :, order_]*orth[-2]) + multiplier_*((var-coeffs[0, :, order_])*orth[-1]-coeffs[1, :, order_]*orth[-2])) orth = numpoly.polynomial(orth[1:]).T norms = numpy.cumprod(coeffs[1], 1) diff --git a/docs/reference/index.rst b/docs/reference/index.rst index e46273bf..2b867488 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -14,7 +14,6 @@ public. descriptives distributions polynomial - orthogonal sampling quadrature recurrence diff --git a/docs/reference/orthogonal.rst b/docs/reference/orthogonal.rst deleted file mode 100644 index 91d57c87..00000000 --- a/docs/reference/orthogonal.rst +++ /dev/null @@ -1,23 +0,0 @@ -Polynomial Expansions -===================== - -.. currentmodule:: chaospy - -Orthogonal Expansions ---------------------- - -.. autosummary:: - :toctree: api - - orth_ttr - orth_chol - orth_gs - - -Non-orthogonal Expansions -------------------------- - -.. autosummary:: - :toctree: api - - lagrange_polynomial diff --git a/docs/reference/polynomial.rst b/docs/reference/polynomial.rst index 61d5ebc6..f4985a69 100644 --- a/docs/reference/polynomial.rst +++ b/docs/reference/polynomial.rst @@ -5,52 +5,79 @@ Polynomials .. currentmodule:: chaospy -Baseclass ---------- +Variable constructors +--------------------- .. autosummary:: - :template: ndpoly.rst - :toctree: api + :toctree: api - ndpoly + variable + polynomial + aspolynomial + symbols + polynomial_from_attributes + ndpoly.from_attributes + polynomial_from_roots -.. autosummary:: - :toctree: api - - ndpoly.coefficients - ndpoly.dtype - ndpoly.exponents - ndpoly.indeterminants - ndpoly.keys - ndpoly.names - ndpoly.values - ndpoly.KEY_OFFSET - -Exceptions +Expansions ---------- .. autosummary:: - :toctree: api + :toctree: api + + monomial + expansion.lagrange + +Orthogonal constructors +~~~~~~~~~~~~~~~~~~~~~~~ + +.. autosummary:: + :toctree: api - FeatureNotSupported + expansion.stieltjes + expansion.cholesky + expansion.gram_schmidt -Constructors ------------- +Pre-defined orthogonal +~~~~~~~~~~~~~~~~~~~~~~ .. autosummary:: - :toctree: api + :toctree: api - variable - polynomial - aspolynomial - monomial - symbols - polynomial_from_attributes - ndpoly.from_attributes - polynomial_from_roots + expansion.chebyshev_1 + expansion.chebyshev_2 + expansion.gegenbauer + expansion.hermite + expansion.jacobi + expansion.laguerre + expansion.legendre + +Baseclass +--------- + +.. autosummary:: + :template: ndpoly.rst + :toctree: api + + ndpoly + +.. autosummary:: + :toctree: api + + ndpoly.coefficients + ndpoly.dtype + ndpoly.exponents + ndpoly.indeterminants + ndpoly.keys + ndpoly.names + ndpoly.values + ndpoly.KEY_OFFSET + +Helper functions +---------------- Leading coefficient -------------------- +~~~~~~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -60,7 +87,7 @@ Leading coefficient sortable_proxy Polynomial specific -------------------- +~~~~~~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -74,7 +101,7 @@ Polynomial specific ndpoly.tonumpy Array creation --------------- +~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -89,7 +116,7 @@ Array creation zeros_like Arithmetics ------------ +~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -106,7 +133,7 @@ Arithmetics square Division --------- +~~~~~~~~ .. autosummary:: :toctree: api @@ -122,7 +149,7 @@ Division true_divide Logic ------ +~~~~~ .. autosummary:: :toctree: api @@ -142,7 +169,7 @@ Logic not_equal Rounding --------- +~~~~~~~~ .. autosummary:: :toctree: api @@ -155,7 +182,7 @@ Rounding round_ Sums/Products -------------- +~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -166,7 +193,7 @@ Sums/Products sum Differentiation ---------------- +~~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -178,7 +205,7 @@ Differentiation hessian Min/Max -------- +~~~~~~~ .. autosummary:: :toctree: api @@ -193,7 +220,7 @@ Min/Max minimum Conditionals ------------- +~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -204,7 +231,7 @@ Conditionals where Save/Load ---------- +~~~~~~~~~ .. autosummary:: :toctree: api @@ -217,7 +244,7 @@ Save/Load savez_compressed Stacking/Splitting ------------------- +~~~~~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -234,7 +261,7 @@ Stacking/Splitting vstack Shape manipulation ------------------- +~~~~~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -251,7 +278,7 @@ Shape manipulation transpose Miscellaneous -------------- +~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -267,7 +294,7 @@ Miscellaneous result_type Global options --------------- +~~~~~~~~~~~~~~ .. autosummary:: :toctree: api @@ -277,11 +304,12 @@ Global options set_options Utilities ---------- +~~~~~~~~~ .. autosummary:: :toctree: api + cross_truncate + FeatureNotSupported glexindex glexsort - cross_truncate diff --git a/docs/reference/quadrature.rst b/docs/reference/quadrature.rst index 7e92008d..fa34e502 100644 --- a/docs/reference/quadrature.rst +++ b/docs/reference/quadrature.rst @@ -40,7 +40,7 @@ Gaussian extensions patterson radau -Gaussian predefined +Predefined Gaussian ------------------- .. autosummary:: diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst index 4a3f67a8..b7b140e7 100644 --- a/docs/user_guide/index.rst +++ b/docs/user_guide/index.rst @@ -11,7 +11,6 @@ User guide sampling.rst quadrature.rst polynomial.rst - orthogonality.rst descriptive.rst zbibliography.rst @@ -72,18 +71,6 @@ The user guide is split into the following topics: :target: ./polynomial.html :align: middle -:ref:`orthogonality` --------------------- - -+-----------------+-----------------------------------------------------------+ -| |orthogonality| | Creation, manipulation and analysis of orthogonal | -| | polynomials for use in model approximations. | -+-----------------+-----------------------------------------------------------+ - -.. |orthogonality| image:: figures/orthogonality.png - :target: ./orthogonality.html - :align: middle - :ref:`quadrature` ----------------- diff --git a/docs/user_guide/orthogonality.rst b/docs/user_guide/orthogonality.rst deleted file mode 100644 index 6bfdb272..00000000 --- a/docs/user_guide/orthogonality.rst +++ /dev/null @@ -1,95 +0,0 @@ -.. _orthogonality: - -Orthogonal Polynomials -====================== - -The core idea of polynomial chaos expansions is that the polynomials used as an -expansion are all mutually orthogonal. The relation is typically written -mathematically as: - -.. math:: - \left\langle \Phi_n, \Phi_m \right\rangle = 0 \qquad n \neq m - -In practice this relation is instead expressed by the equivalent notation using -expected values: - -.. math:: - \mbox E\left(\Phi_n \Phi_m\right) = 0 \qquad n \neq m - -In ``chaospy`` this property can be tested by taking the outer product of two -expansions, and checking if the expected value of the resulting matrix is -diagonal. For example, for a basic monomial:: - - >>> expansion = chaospy.monomial(4) - >>> expansion - polynomial([1, q0, q0**2, q0**3]) - >>> outer_product = chaospy.outer(expansion, expansion) - >>> outer_product - polynomial([[1, q0, q0**2, q0**3], - [q0, q0**2, q0**3, q0**4], - [q0**2, q0**3, q0**4, q0**5], - [q0**3, q0**4, q0**5, q0**6]]) - >>> distribution = chaospy.Normal() - >>> chaospy.E(outer_product, distribution) - array([[ 1., 0., 1., 0.], - [ 0., 1., 0., 3.], - [ 1., 0., 3., 0.], - [ 0., 3., 0., 15.]]) - -In other words, the basic monomial (beyond polynomial order 1) are not -orthogonal. - -But if we replace the basic monomial with an explicit orthogonal polynomial -constructor, we get:: - - >>> expansion = chaospy.generate_expansion(3, distribution) - >>> expansion - polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0]) - >>> outer_product = chaospy.outer(expansion, expansion) - >>> chaospy.E(outer_product, distribution).round(15) - array([[1., 0., 0., 0.], - [0., 1., 0., 0.], - [0., 0., 2., 0.], - [0., 0., 0., 6.]]) - -A fully diagonal matrix, which implies all the polynomials in the expansion are -mutually orthogonal. - -Algorithms ----------- - -There are three algorithms available: - -+------------------------+--------------------------------------------------+ -| Algorithm | Description | -+------------------------+--------------------------------------------------+ -| three_terms_recurrence | Three terms recurrence coefficients generated | -| | using Stieltjes :cite:`stieltjes_quelques_1884` | -| | and Golub-Welsch method | -| | :cite:`golub_calculation_1967`. The most stable | -| | of the methods, but do not work on | -| | dependent distributions. | -+------------------------+--------------------------------------------------+ -| gram_schmidt | Gram-Schmidt orthogonalization method applied on | -| | polynomial expansions. Know for being | -| | numerically unstable. | -+------------------------+--------------------------------------------------+ -| cholesky | Orthogonalization through decorrelation of the | -| | covariance matrix. Uses Gill-King's Cholesky | -| | decomposition method for higher numerical | -| | stability. Still not scalable to high number of | -| | dimensions. | -+------------------------+--------------------------------------------------+ - -.. _lagrange: - -Lagrange Polynomials --------------------- - -Lagrange polynomials are not a method for creating orthogonal polynomials. -Instead it is an interpolation method for creating an polynomial expansion that -has the property that each polynomial interpolates exactly one point in space -with the value 1 and has the value 0 for all other interpolation values. -For more details, see this `article on Lagrange polynomials`_. - -.. _article on Lagrange polynomials: https://en.wikipedia.org/wiki/Lagrange_polynomial diff --git a/docs/user_guide/polynomial.rst b/docs/user_guide/polynomial.rst index 0e094ddf..3cdafbcf 100644 --- a/docs/user_guide/polynomial.rst +++ b/docs/user_guide/polynomial.rst @@ -324,6 +324,102 @@ E.g.: array([[2, 3, 0, 0], [0, 0, 2, 3]]) +.. _orthogonality: + +Orthogonal expansions +--------------------- + +The core idea of polynomial chaos expansions is that the polynomials used as an +expansion are all mutually orthogonal. The relation is typically written +mathematically as: + +.. math:: + \left\langle \Phi_n, \Phi_m \right\rangle = 0 \qquad n \neq m + +In practice this relation is instead expressed by the equivalent notation using +expected values: + +.. math:: + \mbox E\left(\Phi_n \Phi_m\right) = 0 \qquad n \neq m + +In ``chaospy`` this property can be tested by taking the outer product of two +expansions, and checking if the expected value of the resulting matrix is +diagonal. For example, for a basic monomial:: + + >>> expansion = chaospy.monomial(4) + >>> expansion + polynomial([1, q0, q0**2, q0**3]) + >>> outer_product = chaospy.outer(expansion, expansion) + >>> outer_product + polynomial([[1, q0, q0**2, q0**3], + [q0, q0**2, q0**3, q0**4], + [q0**2, q0**3, q0**4, q0**5], + [q0**3, q0**4, q0**5, q0**6]]) + >>> distribution = chaospy.Normal() + >>> chaospy.E(outer_product, distribution) + array([[ 1., 0., 1., 0.], + [ 0., 1., 0., 3.], + [ 1., 0., 3., 0.], + [ 0., 3., 0., 15.]]) + +In other words, the basic monomial (beyond polynomial order 1) are not +orthogonal. + +But if we replace the basic monomial with an explicit orthogonal polynomial +constructor, we get:: + + >>> expansion = chaospy.generate_expansion(3, distribution) + >>> expansion + polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0]) + >>> outer_product = chaospy.outer(expansion, expansion) + >>> chaospy.E(outer_product, distribution).round(15) + array([[1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., 2., 0.], + [0., 0., 0., 6.]]) + +A fully diagonal matrix, which implies all the polynomials in the expansion are +mutually orthogonal. + +Algorithms +~~~~~~~~~~ + +There are three algorithms available: + ++------------------------+--------------------------------------------------+ +| Algorithm | Description | ++------------------------+--------------------------------------------------+ +| three_terms_recurrence | Three terms recurrence coefficients generated | +| | using Stieltjes :cite:`stieltjes_quelques_1884` | +| | and Golub-Welsch method | +| | :cite:`golub_calculation_1967`. The most stable | +| | of the methods, but do not work on | +| | dependent distributions. | ++------------------------+--------------------------------------------------+ +| gram_schmidt | Gram-Schmidt orthogonalization method applied on | +| | polynomial expansions. Know for being | +| | numerically unstable. | ++------------------------+--------------------------------------------------+ +| cholesky | Orthogonalization through decorrelation of the | +| | covariance matrix. Uses Gill-King's Cholesky | +| | decomposition method for higher numerical | +| | stability. Still not scalable to high number of | +| | dimensions. | ++------------------------+--------------------------------------------------+ + +.. _lagrange: + +Lagrange polynomials +-------------------- + +Lagrange polynomials are not a method for creating orthogonal polynomials. +Instead it is an interpolation method for creating an polynomial expansion that +has the property that each polynomial interpolates exactly one point in space +with the value 1 and has the value 0 for all other interpolation values. +For more details, see this `article on Lagrange polynomials`_. + +.. _article on Lagrange polynomials: https://en.wikipedia.org/wiki/Lagrange_polynomial + .. _numpy_functions: Numpy functions diff --git a/pyproject.toml b/pyproject.toml index 2978915b..83d51dde 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "poetry.masonry.api" [tool.poetry] name = "chaospy" -version = "4.3.0" +version = "4.3.1" description = "Numerical tool for perfroming uncertainty quantification" license = "MIT" authors = ["Jonathan Feinberg"] diff --git a/tests/recurrence/test_quadrature_creation.py b/tests/recurrence/test_quadrature_creation.py index 3eaf994b..a782b94e 100644 --- a/tests/recurrence/test_quadrature_creation.py +++ b/tests/recurrence/test_quadrature_creation.py @@ -14,7 +14,7 @@ def test_1d_quadrature_creation( analytical_distribution, recurrence_algorithm): """Check 1-D quadrature rule.""" - abscissas, weights = chaospy.gaussian( + abscissas, weights = chaospy.quadrature.gaussian( order=8, dist=analytical_distribution, recurrence_algorithm=recurrence_algorithm, @@ -35,7 +35,7 @@ def test_3d_quadrature_creation( analytical_distribution, recurrence_algorithm): """Check 3-D quadrature rule.""" distribution = chaospy.Iid(analytical_distribution, 3) - abscissas, weights = chaospy.gaussian( + abscissas, weights = chaospy.quadrature.gaussian( order=3, dist=distribution, recurrence_algorithm=recurrence_algorithm, diff --git a/tests/test_orth.py b/tests/test_orth.py index 4a21d587..82ebd3e0 100644 --- a/tests/test_orth.py +++ b/tests/test_orth.py @@ -18,9 +18,9 @@ def test_operator_E(): assert np.allclose(cp.E(poly, dist), res) -def test_orth_ttr(): +def test_expansion_stieltjes(): dist = cp.Normal(0, 1) - orth = cp.orth_ttr(5, dist) + orth = cp.expansion.stieltjes(5, dist) outer = cp.outer(orth, orth) Cov1 = cp.E(outer, dist) Diatoric = Cov1 - np.diag(np.diag(Cov1)) @@ -30,16 +30,16 @@ def test_orth_ttr(): assert np.allclose(Cov1[1:,1:], Cov2) -def test_orth_chol(): +def test_expansion_cholesky(): dist = cp.Normal(0, 1) - orth1 = cp.orth_ttr(5, dist, normed=True) - orth2 = cp.orth_chol(5, dist, normed=True) + orth1 = cp.expansion.cholesky(5, dist, normed=True) + orth2 = cp.expansion.cholesky(5, dist, normed=True) eps = cp.sum((orth1-orth2)**2) assert np.allclose(eps(np.linspace(-100, 100, 5)), 0) -def test_orth_norms(): +def test_expansion_stieltjes_norms(): dist = cp.Normal(0, 1) - orth = cp.orth_ttr(5, dist, normed=True) + orth = cp.expansion.stieltjes(5, dist, normed=True) norms = cp.E(orth**2, dist) assert np.allclose(norms, 1) diff --git a/tests/test_stress.py b/tests/test_stress.py index e432b380..93b0a7af 100644 --- a/tests/test_stress.py +++ b/tests/test_stress.py @@ -50,16 +50,11 @@ def test_quasimc(): def test_orthogonals(): dist = cp.Iid(cp.Normal(), dim) - cp.orth_gs(order, dist) - cp.orth_ttr(order, dist) - cp.orth_chol(order, dist) + cp.expansion.gram_schmidt(order, dist) + cp.expansion.stieltjes(order, dist) + cp.expansion.cholesky(order, dist) -# def test_approx_orthogonals(): -# dist = cp.Iid(normal(), dim) -# cp.orth_ttr(order, dist) -# - def test_quadrature(): dist = cp.Iid(cp.Normal(), dim) gq = cp.generate_quadrature @@ -78,7 +73,7 @@ def test_approx_quadrature(): def test_integration(): dist = cp.Iid(cp.Normal(), dim) - orth, norms = cp.orth_ttr(order, dist, retall=1) + orth, norms = cp.expansion.stieltjes(order, dist, retall=1) gq = cp.generate_quadrature nodes, weights = gq(order, dist, rule="C") vals = np.zeros((len(weights), size)) @@ -87,7 +82,7 @@ def test_integration(): def test_regression(): dist = cp.Iid(cp.Normal(), dim) - orth, norms = cp.orth_ttr(order, dist, retall=1) + orth, norms = cp.expansion.stieltjes(order, dist, retall=1) data = dist.sample(samples) vals = np.zeros((samples, size)) cp.fit_regression(orth, data, vals) @@ -95,7 +90,7 @@ def test_regression(): def test_descriptives(): dist = cp.Iid(cp.Normal(), dim) - orth = cp.orth_ttr(order, dist) + orth = cp.expansion.stieltjes(order, dist) cp.E(orth, dist) cp.Var(orth, dist) cp.Cov(orth, dist) From 256cafe4a094069ea42ebaca3af1ca625e541ab9 Mon Sep 17 00:00:00 2001 From: Jonathan Feinberg Date: Tue, 6 Apr 2021 15:00:55 +0200 Subject: [PATCH 3/4] including genz-keister back again --- chaospy/distributions/approximation.py | 5 +- chaospy/quadrature/__init__.py | 7 + chaospy/quadrature/frontend.py | 10 + chaospy/quadrature/gaussian.py | 4 +- chaospy/quadrature/genz_keister.py | 300 +++++++++++++++++++++++++ docs/reference/quadrature.rst | 4 + 6 files changed, 325 insertions(+), 5 deletions(-) create mode 100644 chaospy/quadrature/genz_keister.py diff --git a/chaospy/distributions/approximation.py b/chaospy/distributions/approximation.py index 3f3ae897..1724365a 100644 --- a/chaospy/distributions/approximation.py +++ b/chaospy/distributions/approximation.py @@ -156,7 +156,7 @@ def approximate_inverse( def approximate_moment( distribution, k_loc, - order=None, + order=100000, rule="clenshaw_curtis", **kwargs ): @@ -171,8 +171,7 @@ def approximate_moment( k_loc (Sequence[int, ...]): The exponents of the moments of interest with ``shape == (dim,)``. order (int): - The quadrature order used in approximation. If omitted, defaults to - ``1e5**(1./len(distribution))``. + The quadrature order used in approximation. rule (str): Quadrature rule for integrating moments. kwargs: diff --git a/chaospy/quadrature/__init__.py b/chaospy/quadrature/__init__.py index 82d5373d..4dc5c730 100644 --- a/chaospy/quadrature/__init__.py +++ b/chaospy/quadrature/__init__.py @@ -12,6 +12,8 @@ from .fejer_1 import fejer_1 from .fejer_2 import fejer_2 from .gaussian import gaussian +from .genz_keister import ( + genz_keister_16, genz_keister_18, genz_keister_22, genz_keister_24) from .gegenbauer import gegenbauer from .grid import grid from .hermite import hermite @@ -34,6 +36,10 @@ "fejer_1": fejer_1, "fejer_2": fejer_2, "gaussian": gaussian, + "genz_keister_16": genz_keister_16, + "genz_keister_18": genz_keister_18, + "genz_keister_22": genz_keister_22, + "genz_keister_24": genz_keister_24, "grid": grid, "kronrod": kronrod, "legendre": legendre_proxy, @@ -73,3 +79,4 @@ def wrapped(*args, **kwargs): quadrature_deprecation_warning("gauss_lobatto", lobatto) quadrature_deprecation_warning("gauss_patterson", patterson) quadrature_deprecation_warning("gauss_radau", radau) +quadrature_deprecation_warning("genz_keister", genz_keister_24) diff --git a/chaospy/quadrature/frontend.py b/chaospy/quadrature/frontend.py index ef1cb32d..926b9991 100644 --- a/chaospy/quadrature/frontend.py +++ b/chaospy/quadrature/frontend.py @@ -20,6 +20,10 @@ "n": "newton_cotes", "newton_cotes": "newton_cotes", "d": "discrete", "discrete": "discrete", "i": "grid", "grid": "grid", + "z16": "genz_keister_16", "genz_keister_16": "genz_keister_16", + "z18": "genz_keister_18", "genz_keister_18": "genz_keister_18", + "z22": "genz_keister_22", "genz_keister_22": "genz_keister_22", + "z24": "genz_keister_24", "genz_keister_24": "genz_keister_24", } DEPRECATED_SHORT_NAMES = { "f": "f2", @@ -29,6 +33,8 @@ "gauss_patterson": "patterson", "gauss_radau": "radau", "gauss_legendre": "legendre", + "z": "genz_keister_24", + "genz_keister": "genz_keister_24", } @@ -113,6 +119,10 @@ def generate_quadrature( :func:`chaospy.quadrature.newton_cotes` :func:`chaospy.quadrature.discrete` :func:`chaospy.quadrature.grid` + :func:`chaospy.quadrature.genz_keister_16` + :func:`chaospy.quadrature.genz_keister_18` + :func:`chaospy.quadrature.genz_keister_22` + :func:`chaospy.quadrature.genz_keister_24` """ if not rule: diff --git a/chaospy/quadrature/gaussian.py b/chaospy/quadrature/gaussian.py index 5d4776d4..9bfcc615 100644 --- a/chaospy/quadrature/gaussian.py +++ b/chaospy/quadrature/gaussian.py @@ -25,10 +25,10 @@ def gaussian( algorithms exists. Args: - dist (chaospy.Distribution): - The distribution which density will be used as weight function. order (int): The order of the quadrature. + dist (chaospy.Distribution): + The distribution which density will be used as weight function. recurrence_algorithm (str): Name of the algorithm used to generate abscissas and weights. rule (str): diff --git a/chaospy/quadrature/genz_keister.py b/chaospy/quadrature/genz_keister.py new file mode 100644 index 00000000..223f13b2 --- /dev/null +++ b/chaospy/quadrature/genz_keister.py @@ -0,0 +1,300 @@ +""" +Hermite Genz-Keister quadrature rules + +Adapted from John Burkardt's implementation in Matlab +""" +import numpy +import scipy + +from .utils import combine_quadrature + +GENZ_KEISTER_STORE = { + 1: ((0.0000000000000000e+00,), (1.7724538509055159,)), + 3: ((0.0000000000000000e+00, 1.2247448713915889), + (1.1816359006036772, 0.29540897515091930)), + 7: ((0.0000000000000000, 0.52403354748695763, + 1.2247448713915889, 2.9592107790638380), + (0.81310410832613500, 0.23286251787386100, + 0.24557928535031393, 0.0012330680655153448)), + 9: ((0.0000000000000000, 0.52403354748695763, 1.2247448713915889, + 2.0232301911005157, 2.9592107790638380), + (0.45014700975378197, 0.47869428549114124, 0.16811892894767771, + 0.014173117873979098, 1.6708826306882348e-4)), + 17: ((0.0000000000000000, 0.52403354748695763, 0.87004089535290285, + 1.2247448713915889, 1.8357079751751868, 2.0232301911005157, + 2.9592107790638380, 3.6677742159463378, 4.4995993983103881), + (0.47310733504965385, 0.45119803602358544, 0.025155825701712934, + 0.15718298376652240, 0.0034840719346803800, 0.012466519132805918, + 1.8723818949278350e-04, -1.4542843387069391e-06, 3.7463469943051758e-08)), + 19: ((0.0000000000000000, 0.52403354748695763, 0.87004089535290285, + 1.2247448713915889, 1.8357079751751868, 2.0232301911005157, + 2.2665132620567876, 2.9592107790638380, + 3.6677742159463378, 4.4995993983103881), + (0.53788160700510168, 0.36924643368920851, 0.10838861955003017, + 0.11360729895748269, 0.032055243099445879, -0.011232438489069229, + 5.1133174390883855e-03, 1.0656589772852267e-04, + 1.0802767206624762e-06, 1.5295717705322357e-09)), + 31: ((0.0000000000000000, 0.17606414208200893, 0.52403354748695763, + 0.87004089535290285, 1.2247448713915889, 1.5794121348467671, + 1.8357079751751868, 2.0232301911005157, 2.2665132620567876, + 2.5705583765842968, 2.9592107790638380, 3.6677742159463378, + 4.4995993983103881, 5.0360899444730940, + 5.6432578578857449, 6.3759392709822356), + (0.45888839636756751, 0.049855761893293160, 0.35393889029580544, + 0.11594930984853116, 0.10939325071860877, 0.0031210210352682834, + 0.029409427580350787, -0.0098566270434610019, 0.0048385208205502612, + 2.6665159778939428e-05, 1.0541662394746661e-04, 1.0889219692128120e-06, + 1.4055252024722478e-09, 9.0675288231679823e-12, + -2.6304696458548942e-13, 2.2365645607044459e-15)), + 33: ((0.0000000000000000, 0.17606414208200893, 0.52403354748695763, + 0.87004089535290285, 1.2247448713915889, 1.5794121348467671, + 1.8357079751751868, 2.0232301911005157, 2.2665132620567876, + 2.5705583765842968, 2.9592107790638380, 3.6677742159463378, + 4.0292201405043713, 4.4995993983103881, 5.0360899444730940, + 5.6432578578857449, 6.3759392709822356), + (2.4656644932829619e-01, 1.8411696047725790e-01, 3.1208656194697448e-01, + 1.3726521191567551e-01, 9.6913444944583621e-02, 1.3032872699027960e-02, + 2.0435058359107205e-02, -4.9118576123877555e-03, 3.7580026604304793e-03, + 1.4753204901862772e-04, 9.8710009197409173e-05, 1.2245220967158438e-06, + -2.3903343382803510e-08, 2.7547825138935901e-09, -3.4281570530349562e-11, + 4.7219278666417693e-13, -1.7602932805372496e-15)), + 35: ((0.0000000000000000e+00, 1.7606414208200893e-01, 5.2403354748695763e-01, + 8.7004089535290285e-01, 1.2247448713915889e+00, 1.5794121348467671e+00, + 1.8357079751751868e+00, 2.0232301911005157e+00, 2.2665132620567876e+00, + 2.5705583765842968e+00, 2.9592107790638380e+00, 3.3491639537131945e+00, + 3.6677742159463378e+00, 4.0292201405043713e+00, 4.4995993983103881e+00, + 5.0360899444730940e+00, 5.6432578578857449e+00, 6.3759392709822356e+00), + (9.1262675363737921e-04, 3.3988595585585218e-01, 2.6244871488784277e-01, + 1.6371221555735804e-01, 8.0245518147390893e-02, 2.7780508908535097e-02, + 5.5928828911469180e-03, 4.0967527720344047e-03, 1.4515580425155904e-03, + 4.8785399304443770e-04, 6.3328620805617891e-05, 4.8462799737020461e-06, + 4.3737818040926989e-07, 3.7920222392319532e-08, 8.1553721816916897e-10, + 5.4896836948499462e-12, 9.6599466278563243e-15, 1.8684014894510604e-18)), + 37: ((0.000000000000000, 0.214618180588171, 0.524033547486958, + 0.870040895352903, 1.224744871391589, 1.561553427651873, + 1.835707975175187, 2.023230191100516, 2.266513262056788, + 2.597288631188366, 2.959210779063838, 3.315584617593290, + 3.667774215946338, 4.057956316089741, 4.499599398310388, + 4.986551454150765, 5.521865209868350, 6.124527854622158, + 6.853200069757519), + (0.968824552928425499e-01, 0.147655710402686249e+00, + 0.143099302896833389e+00, 0.937208280655245902e-01, + 0.442116442189845444e-01, 0.15513109874859354e-01, + 0.43334988122723492e-02, 0.176802225818295443e-02, + 0.985827582996483824e-03, 0.234940366465975222e-03, + 0.32265185983739747e-04, 0.330975870979203419e-05, + 0.295907520230744049e-06, 0.16595448809389819e-07, + 0.422525843963111041e-09, 0.45661763676186859e-11, + 0.182242751549129356e-13, 0.187781893143728947e-16, + 0.19030350940130498e-20)), + 41: ((0.0000000000000000, 0.195324784415805, 0.52403354748695763, + 0.87004089535290285, 1.2247448713915889, 1.585873011819188, + 1.8357079751751868, 2.0232301911005157, 2.043834754429505, + 2.2665132620567876, 2.630415236459871, 2.9592107790638380, + 3.296114596212218, 3.6677742159463378, 4.070919267883068, + 4.4995993983103881, 4.953574342912980, 5.437443360177798, + 5.961461043404500, 6.547083258397540, 7.251792998192644), + (0.562793426043218877e-01, 0.165639740400529554e+00, + 0.145966293895926429e+00, 0.928338228510111845e-01, + 0.45109010335859128e-01, 0.165445526705860772e-01, + 0.705471110122962612e-03, 0.178852543033699732e-01, + - 0.144528422206988237e-01, 0.140697424065246825e-02, + 0.189010909805097887e-03, 0.288976780274478689e-04, + 0.381182791749177506e-05, 0.315372265852264871e-06, + 0.149158210417831408e-07, 0.400784141604834759e-09, + 0.581803393170320419e-11, 0.408820161202505983e-13, + 0.1140700785308509e-15, 0.860427172512207236e-19, + 0.664195893812757801e-23)), + 43: ((0.0000000000000000, 0.196029453662011, 0.52403354748695763, + 0.87004089535290285, 1.2247448713915889, 1.583643465293944, + 1.8357079751751868, 2.0232301911005157, 2.089340389294661, + 2.2665132620567876, 2.633356763661946, 2.9592107790638380, + 3.295265921534226, 3.6677742159463378, 4.071335874253583, + 4.4995993983103881, 4.952329763008589, 5.434053000365068, + 5.954781975039809, 6.535398426382995, 7.231746029072501, + 10.167574994881873), + (0.579595986101181095e-01, 0.164880913687436689e+00, + 0.145863292632147353e+00, 0.928711584442575456e-01, + 0.450612329041864976e-01, 0.163616873493832402e-01, + 0.139966252291568061e-02, 0.67354758901013295e-02, + -0.38799558623877157e-02, 0.150909333211638847e-02, + 0.184789465688357423e-03, 0.286802318064777813e-04, + 0.383880761947398577e-05, 0.316018363221289247e-06, + 0.148653643571796457e-07, 0.400030575425776948e-09, + 0.586915885251734856e-11, 0.421921851448196032e-13, + 0.122619614947864357e-15, 0.992619971560149097e-19, + 0.87544909871323873e-23, 0.546191947478318097e-37)), +} + +RULES = { + 16: [1, 3, 7, 9, 17, 19, 31], + 18: [1, 3, 9, 19, 37], + 22: [1, 3, 9, 19, 41], + 24: [1, 3, 9, 19, 43], +} + +def genz_keister_16(order, dist=None): + """ + Create Genz-Keister variant 16 quadrature nodes and weights. + + Args: + order (int, Sequence[int]): + The order of the quadrature. + dist (Optional[chaospy.Distribution]): + The distribution which density will be used as weight function. + If omitted, standard Gaussian is assumed. + + Returns: + (numpy.ndarray, numpy.ndarray): + Genz-Keister quadrature abscissas and weights. + + Examples: + >>> nodes, weights = genz_keister_16(4) + >>> nodes.round(2) + array([[-6.36, -5.19, -4.18, -2.86, -2.6 , -1.73, -1.23, -0.74, 0. , + 0.74, 1.23, 1.73, 2.6 , 2.86, 4.18, 5.19, 6.36]]) + >>> weights.round(8) + array([ 2.0000000e-08, -8.2000000e-07, 1.0564000e-04, 7.0334800e-03, + 1.9656800e-03, 8.8681000e-02, 1.4192650e-02, 2.5456123e-01, + 2.6692223e-01, 2.5456123e-01, 1.4192650e-02, 8.8681000e-02, + 1.9656800e-03, 7.0334800e-03, 1.0564000e-04, -8.2000000e-07, + 2.0000000e-08]) + + """ + return genz_keister(order, dist, rule=16) + + +def genz_keister_18(order, dist=None): + """ + Create Genz-Keister variant 18 quadrature nodes and weights. + + Args: + order (int, Sequence[int]): + The order of the quadrature. + dist (Optional[chaospy.Distribution]): + The distribution which density will be used as weight function. + If omitted, standard Gaussian is assumed. + + Returns: + (numpy.ndarray, numpy.ndarray): + Genz-Keister quadrature abscissas and weights. + + Examples: + >>> nodes, weights = genz_keister_18(2) + >>> nodes.round(2) + array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]]) + >>> weights.round(8) + array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01, + 2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03, + 9.4270000e-05]) + + """ + return genz_keister(order, dist, rule=18) + + +def genz_keister_22(order, dist=None): + """ + Create Genz-Keister variant 22 quadrature nodes and weights. + + Args: + order (int, Sequence[int]): + The order of the quadrature. + dist (Optional[chaospy.Distribution]): + The distribution which density will be used as weight function. + If omitted, standard Gaussian is assumed. + + Returns: + (numpy.ndarray, numpy.ndarray): + Genz-Keister quadrature abscissas and weights. + + Examples: + >>> nodes, weights = genz_keister_22(2) + >>> nodes.round(2) + array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]]) + >>> weights.round(8) + array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01, + 2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03, + 9.4270000e-05]) + + """ + return genz_keister(order, dist, rule=22) + + +def genz_keister_24(order, dist=None): + """ + Create Genz-Keister variant 24 quadrature nodes and weights. + + Args: + order (int, Sequence[int]): + The order of the quadrature. + dist (Optional[chaospy.Distribution]): + The distribution which density will be used as weight function. + If omitted, standard Gaussian is assumed. + + Returns: + (numpy.ndarray, numpy.ndarray): + Genz-Keister quadrature abscissas and weights. + + Examples: + >>> nodes, weights = genz_keister_24(2) + >>> nodes.round(2) + array([[-4.18, -2.86, -1.73, -0.74, 0. , 0.74, 1.73, 2.86, 4.18]]) + >>> weights.round(8) + array([9.4270000e-05, 7.9963300e-03, 9.4850950e-02, 2.7007433e-01, + 2.5396825e-01, 2.7007433e-01, 9.4850950e-02, 7.9963300e-03, + 9.4270000e-05]) + + """ + return genz_keister(order, dist, rule=24) + + +def genz_keister(order, dist=None, rule=24): + """ + Create Genz-Keister quadrature nodes and weights. + + Args: + order (int, Sequence[int]): + The order of the quadrature. + dist (Optional[chaospy.Distribution]): + The distribution which density will be used as weight function. + If omitted, standard Gaussian is assumed. + rule (int, Sequence[int]): + The Genz-Keister rule name. Supported rules are 16, 18, 22 and 24. + + Returns: + (numpy.ndarray, numpy.ndarray): + Genz-Keister quadrature abscissas and weights. + + Examples: + >>> genz_keister(0) + (array([[0.]]), array([1.])) + >>> genz_keister(1) # doctest: +NORMALIZE_WHITESPACE + (array([[-1.73205081, 0. , 1.73205081]]), + array([0.16666667, 0.66666667, 0.16666667])) + + """ + shape = (1,) if dist is None else (len(dist),) + order = numpy.broadcast_to(order, shape) + rule = numpy.broadcast_to(rule, shape) + nodes, weights = zip(*[_genz_keister(order_, rule_) + for order_, rule_ in zip(order, rule)]) + nodes, weights = combine_quadrature(nodes, weights) + if dist is not None: + nodes = dist.inv(scipy.special.ndtr(nodes)) + return nodes, weights + + +def _genz_keister(order, rule): + assert rule in RULES, "rule %d not in known rules: %s" % (rule, list(RULES)) + assert order <= len(RULES[rule]), ( + "rule genz_keister_%d limited at order %d" % (rule, order)) + order = RULES[rule][order] + nodes, weights = GENZ_KEISTER_STORE[order] + length = len(nodes) + nodes = numpy.array(nodes[::-1]+nodes[1:]) + nodes[:length-1] *= -1 + nodes *= numpy.sqrt(2) + weights = numpy.array(weights[::-1]+weights[1:]) + weights /= numpy.sum(weights) + + return nodes, weights diff --git a/docs/reference/quadrature.rst b/docs/reference/quadrature.rst index fa34e502..b7184ca4 100644 --- a/docs/reference/quadrature.rst +++ b/docs/reference/quadrature.rst @@ -34,6 +34,10 @@ Gaussian extensions .. autosummary:: :toctree: api + genz_keister_16 + genz_keister_18 + genz_keister_22 + genz_keister_24 legendre kronrod lobatto From 340cd8c57765480145f76f605548d21495913ed9 Mon Sep 17 00:00:00 2001 From: Jonathan Feinberg Date: Thu, 22 Apr 2021 08:58:52 +0200 Subject: [PATCH 4/4] some new stuff --- .../distributions/baseclass/distribution.py | 14 +++++-- chaospy/quadrature/frontend.py | 2 +- docs/reference/high_level_interface.rst | 37 ------------------- docs/reference/index.rst | 3 +- docs/reference/polynomial.rst | 10 +++++ docs/reference/quadrature.rst | 5 +++ docs/reference/recurrence.rst | 1 + docs/reference/sampling.rst | 4 +- docs/user_guide/index.rst | 34 ++++++++--------- docs/user_guide/sampling.rst | 6 +-- pyproject.toml | 3 +- 11 files changed, 51 insertions(+), 68 deletions(-) delete mode 100644 docs/reference/high_level_interface.rst diff --git a/chaospy/distributions/baseclass/distribution.py b/chaospy/distributions/baseclass/distribution.py index 5d3ddccf..72d5fcc0 100644 --- a/chaospy/distributions/baseclass/distribution.py +++ b/chaospy/distributions/baseclass/distribution.py @@ -1,4 +1,5 @@ """Abstract baseclass for all distributions.""" +import typing import logging import numpy @@ -15,8 +16,9 @@ class Distribution(object): interpret_as_integer = False """ - Flag indicating that return value from the methods sample, and inv - should be interpreted as integers instead of floating point. + Flag indicating that return value from the methods + :func:`Distribution.sample`, and :func:`Distribution.inv` should be + interpreted as integers instead of floating point. """ @property @@ -204,6 +206,10 @@ def fwd(self, x_data): q_data = q_data.reshape(shape) return q_data + def _cdf(self, x_data, **parameters): + raise NotImplementedError( + "%s most define _cdf method." % self.__class__.__name__) + def _get_fwd(self, x_data, idx, cache): """In-process function for getting cdf-values.""" logger = logging.getLogger(__name__) @@ -691,9 +697,9 @@ def _cache(self, idx, cache, get): """Backend function of retrieving cache values.""" return self - def __getitem__(self, index): + def __getitem__(self, index: typing.Union[int, numpy.ndarray, slice]): if isinstance(index, numpy.number): - assert index.dtype == int + assert numpy.asarray(index).dtype == int index = int(index) if isinstance(index, int): if not -len(self) < index < len(self): diff --git a/chaospy/quadrature/frontend.py b/chaospy/quadrature/frontend.py index 926b9991..fdddb139 100644 --- a/chaospy/quadrature/frontend.py +++ b/chaospy/quadrature/frontend.py @@ -26,7 +26,7 @@ "z24": "genz_keister_24", "genz_keister_24": "genz_keister_24", } DEPRECATED_SHORT_NAMES = { - "f": "f2", + "f": "fejer_2", "fejer": "fejer_2", "gauss_kronrod": "kronrod", "gauss_lobatto": "lobatto", diff --git a/docs/reference/high_level_interface.rst b/docs/reference/high_level_interface.rst deleted file mode 100644 index 6aee9e5a..00000000 --- a/docs/reference/high_level_interface.rst +++ /dev/null @@ -1,37 +0,0 @@ -High-level interfaces -===================== - -.. currentmodule:: chaospy - -Orthogonal expansion --------------------- - -.. autosummary:: - :toctree: api - - generate_expansion - -Quadrature ----------- - -.. autosummary:: - :toctree: api - - generate_quadrature - -Fit chaos expansion -------------------- - -.. autosummary:: - :toctree: api - - fit_regression - fit_quadrature - -Three terms recurrence coefficients ------------------------------------ - -.. autosummary:: - :toctree: api - - construct_recurrence_coefficients diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 2b867488..62d75e61 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -10,10 +10,9 @@ public. .. toctree:: :maxdepth: 1 - high_level_interface - descriptives distributions polynomial sampling quadrature recurrence + descriptives diff --git a/docs/reference/polynomial.rst b/docs/reference/polynomial.rst index f4985a69..c97a86c5 100644 --- a/docs/reference/polynomial.rst +++ b/docs/reference/polynomial.rst @@ -34,10 +34,20 @@ Orthogonal constructors .. autosummary:: :toctree: api + generate_expansion expansion.stieltjes expansion.cholesky expansion.gram_schmidt +Model builder +~~~~~~~~~~~~~ + +.. autosummary:: + :toctree: api + + fit_quadrature + fit_regression + Pre-defined orthogonal ~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/reference/quadrature.rst b/docs/reference/quadrature.rst index b7184ca4..002ad942 100644 --- a/docs/reference/quadrature.rst +++ b/docs/reference/quadrature.rst @@ -3,6 +3,11 @@ Quadrature integration ====================== +.. autosummary:: + :toctree: api + + chaospy.generate_quadrature + .. currentmodule:: chaospy.quadrature Standard library diff --git a/docs/reference/recurrence.rst b/docs/reference/recurrence.rst index 62f63871..0a750056 100644 --- a/docs/reference/recurrence.rst +++ b/docs/reference/recurrence.rst @@ -9,6 +9,7 @@ Recurrence algorithms .. autosummary:: :toctree: api + construct_recurrence_coefficients modified_chebyshev lanczos stieltjes diff --git a/docs/reference/sampling.rst b/docs/reference/sampling.rst index 17f2cd68..0f00da69 100644 --- a/docs/reference/sampling.rst +++ b/docs/reference/sampling.rst @@ -1,7 +1,7 @@ .. _sampling_collection: -Random and low-discrepency samples -================================== +Sampling and sequences +====================== .. currentmodule:: chaospy diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst index b7b140e7..a9a906fc 100644 --- a/docs/user_guide/index.rst +++ b/docs/user_guide/index.rst @@ -22,18 +22,6 @@ check out the :ref:`tutorial` section as well. The user guide is split into the following topics: -:ref:`chaos_expansion` ----------------------- - -+-----------------+-----------------------------------------------------------+ -| |chaos| | Overview over the fundamentals of using `chaospy` to make | -| | polynomial chaos expansions. | -+-----------------+-----------------------------------------------------------+ - -.. |chaos| image:: figures/recurrence.png - :target: ./chaos_expansion.html - :align: middle - :ref:`distributions` -------------------- @@ -59,6 +47,17 @@ The user guide is split into the following topics: :target: ./sampling.html :align: middle +:ref:`quadrature` +----------------- + ++-----------------+-----------------------------------------------------------+ +| |quadrature| | Quadrature rules for numerical integration. | ++-----------------+-----------------------------------------------------------+ + +.. |quadrature| image:: figures/quadrature.png + :target: ./quadrature.html + :align: middle + :ref:`polynomial` ----------------- @@ -71,15 +70,16 @@ The user guide is split into the following topics: :target: ./polynomial.html :align: middle -:ref:`quadrature` ------------------ +:ref:`chaos_expansion` +---------------------- +-----------------+-----------------------------------------------------------+ -| |quadrature| | Quadrature rules for numerical integration. | +| |chaos| | Overview over the fundamentals of using `chaospy` to make | +| | polynomial chaos expansions. | +-----------------+-----------------------------------------------------------+ -.. |quadrature| image:: figures/quadrature.png - :target: ./quadrature.html +.. |chaos| image:: figures/recurrence.png + :target: ./chaos_expansion.html :align: middle :ref:`descriptives` diff --git a/docs/user_guide/sampling.rst b/docs/user_guide/sampling.rst index b0ad7fdd..010c6248 100644 --- a/docs/user_guide/sampling.rst +++ b/docs/user_guide/sampling.rst @@ -1,7 +1,7 @@ .. _sampling: -Sampling and low-discrepency sequences -====================================== +Sampling and sequences +====================== Introduction ------------ @@ -135,7 +135,7 @@ advantage. .. _antithetic: -Antithetic Variates +Antithetic variates ------------------- Create `antithetic variates`_ from variables on the unit hyper-cube. diff --git a/pyproject.toml b/pyproject.toml index 83d51dde..f7deebdd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,11 +20,10 @@ classifiers = [ ] [tool.poetry.dependencies] -python = "^2.7||>=3.6" +python = ">=3.6" numpoly = "*" numpy = "*" scipy = "*" -functools32 = { version = "*", python = "^2.7" } [tool.poetry.dev-dependencies] pytest = "*"