ENH: Modify the various polynomial series so that multidimensional coefficient arrays can be used. Add functions for evaluation of 2D and 3D polynomial series evaluated either on a specified set of points or on a cartesian product of 1D points.

The new functions have names polyval2d, polygrid2d, polyval3d, and
polygrid3d, where 'poly' can be replaced by any of 'leg', 'cheb', 'lag',
'herm', or 'herme'. These additional functions should cover the common
multidimensional cases and provide examples for anyone who wants to go to
higher dimensions.

1 files changed, 222 insertions, 36 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index b4f50d90e..6212f2bc5 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -25,6 +25,10 @@ Arithmetic
 - `chebdiv` -- divide one Chebyshev series by another.
 - `chebpow` -- raise a Chebyshev series to an positive integer power
 - `chebval` -- evaluate a Chebyshev series at given points.
+- `chebval2d` -- evaluate a 2D Chebyshev series at given points.
+- `chebval3d` -- evaluate a 3D Chebyshev series at given points.
+- `chebgrid2d` -- evaluate a 2D Chebyshev series on a Cartesian product.
+- `chebgrid3d` -- evaluate a 3D Chebyshev series on a Cartesian product.
 
 Calculus
 --------
@@ -78,18 +82,19 @@ References
 """
 from __future__ import division
 
-__all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
-        'chebadd', 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow',
-        'chebval', 'chebder', 'chebint', 'cheb2poly', 'poly2cheb',
-        'chebfromroots', 'chebvander', 'chebfit', 'chebtrim', 'chebroots',
-        'chebpts1', 'chebpts2', 'Chebyshev']
-
 import numpy as np
 import numpy.linalg as la
 import polyutils as pu
 import warnings
 from polytemplate import polytemplate
 
+__all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
+    'chebadd', 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow',
+    'chebval', 'chebval2d', 'chebval3d', 'chebgrid2d', 'chebgrid3d',
+    'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots',
+    'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1',
+    'chebpts2', 'Chebyshev']
+
 chebtrim = pu.trimcoef
 
 #
@@ -1025,36 +1030,59 @@ def chebint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def chebval(x, cs):
-    """Evaluate a Chebyshev series.
+def chebval(x, c, tensor=True):
+    """
+    Evaluate a Chebyshev series.
+
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    If `cs` is of length `n`, this function returns :
+    ``p(x) = c[0]*T_0(x) + c[1]*T_1(x) + ... + c[n]*T_n(x)``
 
-    ``p(x) = cs[0]*T_0(x) + cs[1]*T_1(x) + ... + cs[n-1]*T_{n-1}(x)``
+    If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
 
-    If x is a sequence or array then p(x) will have the same shape as x.
-    If r is a ring_like object that supports multiplication and addition
-    by the values in `cs`, then an object of the same type is returned.
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     Parameters
     ----------
-    x : array_like, ring_like
-        Array of numbers or objects that support multiplication and
-        addition with themselves and with the elements of `cs`.
-    cs : array_like
-        1-d array of Chebyshev coefficients ordered from low to high.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Chebyshev series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    chebfit
-
-    Examples
-    --------
+    chebval2d, chebgrid2d, chebval3d, chebgrid3d
 
     Notes
     -----
@@ -1064,28 +1092,186 @@ def chebval(x, cs):
     --------
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
         x2 = 2*x
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
-            c0 = cs[-i] - c1
+            c0 = c[-i] - c1
             c1 = tmp + c1*x2
     return c0 + c1*x
 
 
+def chebval2d(x, y, c):
+    """
+    Evaluate 2D Chebyshev series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * T_i(x) * T_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Chebyshev seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Chebyshev series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    chebval, chebgrid2d, chebval3d, chebgrid3d
+
+    """
+    return chebval(y, chebval(x, c), False)
+
+
+def chebgrid2d(x, y, c):
+    """
+    Evaluate 2D Chebyshev series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * T_i(a) * T_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Chebyshev series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Chebyshev series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebval3d, chebgrid3d
+
+    """
+    return chebval(y, chebval(x, c))
+
+
+def chebval3d(x, y, z, c):
+    """
+    Evaluate 3D Chebyshev series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * T_i(x) * T_j(y) * T_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Chebyshev seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `chebval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Chebyshev series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebgrid2d, chebgrid3d
+
+    """
+    return chebval(z, chebval2d(x, y, c), False)
+
+
+def chebgrid3d(x, y, z, c):
+    """
+    Evaluate 3D Chebyshev series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * T_i(a) * T_j(b) * T_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Chebyshev seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Chebyshev series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebgrid2d, chebval3d
+
+    """
+    return chebval(z, chebgrid2d(x, y, c))
+
+
 def chebvander(x, deg) :
     """Vandermonde matrix of given degree.