ENH: Make derivatives and integrals work on multidimensional array.

An axis keyword was added to the function signatures of xxxder and
xxxint, where xxx is any of poly, cheb, leg, lag, herm, herme. The
evaluation method for the Chebeshev series was also changed to avoid
using z_series and to more closely resemble the other implementations.

At some point the z_series will be removed from the chebyshev module
and only used for trigonometric series.

1 files changed, 87 insertions, 52 deletions

diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 39463b0cb..07f943bb1 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -621,26 +621,33 @@ def lagpow(cs, pow, maxpower=16) :
         return prd
 
 
-def lagder(cs, m=1, scl=1) :
+def lagder(c, m=1, scl=1, axis=0) :
     """
     Differentiate a Laguerre series.
 
-    Returns the series `cs` differentiated `m` times.  At each iteration the
-    result is multiplied by `scl` (the scaling factor is for use in a linear
-    change of variable).  The argument `cs` is the sequence of coefficients
-    from lowest order "term" to highest, e.g., [1,2,3] represents the series
-    ``P_0 + 2*P_1 + 3*P_2``.
+    Returns the Laguerre series coefficients `c` differentiated `m` times
+    along `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
+    while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
+    2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of Laguerre series coefficients ordered from low to high.
+    c: array_like
+        Array of Laguerre series coefficients. If `c` is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     m : int, optional
         Number of derivatives taken, must be non-negative. (Default: 1)
     scl : scalar, optional
         Each differentiation is multiplied by `scl`.  The end result is
         multiplication by ``scl**m``.  This is for use in a linear change of
         variable. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -667,49 +674,66 @@ def lagder(cs, m=1, scl=1) :
     array([ 1.,  2.,  3.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[:1]*0
     else :
         for i in range(cnt):
-            n = len(cs) - 1
-            cs *= scl
-            der = np.empty(n, dtype=cs.dtype)
-            for j in range(n, 0, -1):
-                der[j - 1] = -cs[j]
-                cs[j - 1] += cs[j]
-            cs = der
-        return cs
-
-
-def lagint(cs, m=1, k=[], lbnd=0, scl=1):
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
+            for j in range(n, 1, -1):
+                der[j - 1] = -c[j]
+                c[j - 1] += c[j]
+            der[0] = -c[1]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
+
+
+def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     Integrate a Laguerre series.
 
-    Returns a Laguerre series that is the Laguerre series `cs`, integrated
-    `m` times from `lbnd` to `x`.  At each iteration the resulting series
-    is **multiplied** by `scl` and an integration constant, `k`, is added.
+    Returns the Laguerre series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. At each iteration the resulting series is
+    **multiplied** by `scl` and an integration constant, `k`, is added.
     The scaling factor is for use in a linear change of variable.  ("Buyer
     beware": note that, depending on what one is doing, one may want `scl`
     to be the reciprocal of what one might expect; for more information,
-    see the Notes section below.)  The argument `cs` is a sequence of
-    coefficients, from lowest order Laguerre series "term" to highest,
-    e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
+    represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
+    2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
+
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of Laguerre series coefficients, ordered from low to high.
+    c : array_like
+        Array of Laguerre series coefficients. If `c` is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -723,6 +747,8 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -767,9 +793,12 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
     array([ 11.16666667,  -5.        ,  -3.        ,   2.        ])
 
     """
-    cnt = int(m)
-    if np.isscalar(k) :
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -777,28 +806,34 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     k = list(k) + [0]*(cnt - len(k))
     for i in range(cnt) :
-        n = len(cs)
-        cs *= scl
-        if n == 1 and cs[0] == 0:
-            cs[0] += k[i]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            tmp = np.empty(n + 1, dtype=cs.dtype)
-            tmp[0] = cs[0]
-            tmp[1] = -cs[0]
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]
+            tmp[1] = -c[0]
             for j in range(1, n):
-                tmp[j] += cs[j]
-                tmp[j + 1] = -cs[j]
+                tmp[j] += c[j]
+                tmp[j + 1] = -c[j]
             tmp[0] += k[i] - lagval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def lagval(x, c, tensor=True):
@@ -871,8 +906,8 @@ def lagval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor: