From 5d7065465e13ac0b1bf404ba12cd0572f9e76c23 Mon Sep 17 00:00:00 2001 From: Jonathan Feinberg Date: Thu, 18 Mar 2021 10:18:08 +0100 Subject: [PATCH] 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 eee68ce2..77b8a4f1 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 4564db1b..39d86804 100644 --- a/chaospy/distributions/operators/joint.py +++ b/chaospy/distributions/operators/joint.py @@ -225,18 +225,3 @@ def _cache(self, idx, cache): 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..2e9ed721 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)