refactor orthogonal->expansion

Author: jonathf

chaospy/__init__.py

@@ -12,15 +12,15 @@
 
 import chaospy.descriptives
 import chaospy.distributions
-import chaospy.orthogonal
+import chaospy.expansion
 import chaospy.spectral
 import chaospy.quadrature
 import chaospy.saltelli
 import chaospy.regression
 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 *

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)

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

chaospy/orthogonal/cholesky.py 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

chaospy/orthogonal/frontend.py 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,
 }
 
 

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

chaospy/orthogonal/gram_schmidt.py 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])
 

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

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